Hi Bradley, thank you for your interest in this problem. I have tested the following matlab script and it should run needing only kwave on the path. It inlines all the necessary functions, and replicates the condition I stated. You can run it by calling
G = kwave_inversion_demo;
then inspecting the object G which is the shear modulus. The correct answer is 4000 and it is mostly very close, however you will see some smooth artifact comes through that is far off, and needless to say, the material is supposed to be homogeneous.
Okay here is the code:
-----
function G = kwave_inversion_demo
K = 5; G = 2; dx = .002; f = 50; ncycles = 4';
prefs = kwave_prefs('K', K, 'G', G, 'dx', dx, 'f', f, 'ncycles', ncycles);
u = kwave_eb(prefs);
u = cat(5, u.ux, u.uy, u.uz);
u_t = size(u,4);
final_cycle_start = round(u_t*(ncycles-1)/ncycles);
u_cycle = u(:,:,:,final_cycle_start:end,:);
% temporal fourier transform, retain driving frequency bin
u_ft = fft(u_cycle, [], 4);
u_ft = squeeze(u_ft(:,:,:,2,:));
% apply Helmholtz equation to complex wavefield
u_lap = lap(u_ft) ./ (prefs.dx)^2;
G = prefs.rho*(2*pi*prefs.f)^2*sum(abs(u_ft), 4) ./ sum(abs(u_lap), 4);
end
function u = kwave_eb(prefs)
SIZE = 96;
PML_SIZE = 16;
SOURCE_POS = SIZE / 2;
Nx = SIZE; % number of grid points in the x direction
Ny = SIZE; % number of grid points in the y direction
Nz = SIZE; % number of grid points in the z direction
dx = prefs.dx; % grid point spacing in the x direction [m]
dy = prefs.dx; % grid point spacing in the y direction [m]
dz = prefs.dx; % grid point spacing in the z direction [m]
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);
medium.sound_speed_compression = prefs.K*ones(Nx, Ny, Nz); % [m/s]
medium.sound_speed_shear = prefs.G*ones(Nx, Ny, Nz); % [m/s]
medium.density = prefs.rho*ones(Nx, Ny, Nz); % [kg/m^3]
% define the absorption properties
medium.alpha_coeff_compression = 0.1; % [dB/(MHz^2 cm)]
medium.alpha_coeff_shear = 0.1; % [dB/(MHz^2 cm)]
% create the time array
t_end = prefs.t;
kgrid.t_array = makeTime(kgrid, max(medium.sound_speed_compression(:)), [], t_end);
source.u_mask = zeros(Nx, Ny, Nz);
source.u_mask(SOURCE_POS, SOURCE_POS, SOURCE_POS) = 1;
% define source to be a velocity source
source_freq = prefs.f; % [Hz]
source_mag = 1; % [Pa]
source.ux = source_mag*sin(2*pi*source_freq*kgrid.t_array);
source.uy = source_mag*sin(2*pi*source_freq*kgrid.t_array);
source.uz = source_mag*sin(2*pi*source_freq*kgrid.t_array);
%source.ux = filterTimeSeries(kgrid, medium, source.ux);
%source.uy = filterTimeSeries(kgrid, medium, source.uy);
%source.uz = filterTimeSeries(kgrid, medium, source.uz);
% define sensor mask in x-y plane using cuboid corners, where a rectangular
% mask is defined using the xyz coordinates of two opposing corners in the
% form [x1, y1, z1, x2, y2, z2].'
sensor.mask = [1 + PML_SIZE, 1 + PML_SIZE, 1 + PML_SIZE, Nx - PML_SIZE, Ny - PML_SIZE, Nz - PML_SIZE].';
% record the maximum pressure in the plane
sensor.record = {'u'};
% define input arguments
input_args = {'PMLSize', PML_SIZE};
% run the simulation with PML inside
u = pstdElastic3D(kgrid, medium, source, sensor, input_args{:});
end
function prefs = kwave_prefs(varargin)
prefs = default_prefs;
fields = get_fields;
for n = 1:2:nargin
name = varargin{n};
value = varargin{n+1};
for f = 1:numel(fields)
if strcmpi(fields{f}, name)
evalc(['prefs.',name,'=',num2str(value),';']);
end
end
end
end
function prefs = default_prefs
prefs.K = 1540;
prefs.G = 1.54;
prefs.dx = 0.002;
prefs.rho = 1000;
prefs.f = 50;
prefs.ncycles = 4;
prefs.t = (1/prefs.f)*prefs.ncycles;
end
function fields = get_fields
fields = {'K', 'G', 'dx', 'rho', 'f', 'ncycles', 't'};
end
function y = lap(x)
laplacian = zeros(3,3,3);
laplacian(2,2,1) = 1;
laplacian(2,2,3) = 1;
laplacian(:,:,2) = [0 1 0; 1 -6 1; 0 1 0];
y = convn(x, laplacian, 'same');
end