Hi Kevin,
To get the same results as k-Wave, you can set w0 to zero, and c0 to the local sound speed. See below for an example based on the Modelling Power Law Absorption example.
Two things to keep in mind. First, when using kspaceFirstOrderND there will be some error introduced in the calculation of the absorption term as discussed here. Second, if using high values of absorption, second order effects mean the absorption and dispersion deviate slightly from a power law - see Eq. (38) and (40) in this paper.
Brad.
--
clear all;
% create the computational grid
Nx = 1024; % number of grid points in the x (row) direction
x = 12.8e-3; % grid size in the x direction [m]
dx = x/Nx; % grid point spacing in the x direction [m]
kgrid = makeGrid(Nx, dx);
% define the properties of the propagation medium
medium.sound_speed = 1500; % [m/s]
medium.alpha_coeff = 0.2;
medium.alpha_power = 1.5;
% define time array
t_end = 4e-6;
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed, [], t_end);
% create a spatial delta pulse
source_pos = Nx/4; % [grid points]
source.p0 = zeros(Nx, 1);
source.p0(source_pos) = 1;
% define the sensor positions
source_sensor_dist = 0.5e-3; % [m]
sensor_sensor_dist = 1e-3; % [m]
sensor_pos_1 = source_pos + round(source_sensor_dist/dx);
sensor_pos_2 = source_pos + round((source_sensor_dist + sensor_sensor_dist)/dx);
% calculate discrete distance between the sensor positions
d = (sensor_pos_2 - sensor_pos_1)*dx; % [m]
d_cm = d*100;
% create a Binary sensor mask
sensor.mask = zeros(Nx, 1);
sensor.mask(sensor_pos_1) = 1;
sensor.mask(sensor_pos_2) = 1;
% run the simulation without visualisation
sensor_data = kspaceSecondOrder(kgrid, medium, source, sensor, 'PlotSim', false);
% calculate the amplitude and phase spectrum at the two sensor
% positions
[f, as1, ps1] = spect(sensor_data(1, :), 1/dt);
[~, as2, ps2] = spect(sensor_data(2, :), 1/dt);
% calculate the attenuation from the amplitude spectrums
attenuation = -20*log10(as2./as1)./d_cm;
% calculate the corresponding theoretical attenuation in dB/cm
attenuation_th = medium.alpha_coeff.*(f./1e6).^medium.alpha_power;
% calculate the dispersion (dependence of the sound speed on frequency)
% from the phase spectrums
cp = 2*pi.*f.*d./(unwrap(ps1) - unwrap(ps2));
% calculate the corresponding theoretical dispersion using the
% Kramers-Kronig relation for power law absorption
cp_kk = powerLawKramersKronig(2*pi*f, 0, medium.sound_speed, db2neper(medium.alpha_coeff, medium.alpha_power), medium.alpha_power);
% plot the attenuation
figure;
ds = 8;
f_max = 50;
plot(f(1:ds:end)./1e6, attenuation(:, 1:ds:end), 'ko', f./1e6, attenuation_th, 'k-');
set(gca, 'XLim', [0 f_max]);
box on;
xlabel('Frequency [MHz]');
ylabel('\alpha [dB/cm]');
% plot the dispersion
figure
plot(f(1:ds:end)./1e6, cp(:, 1:ds:end), 'ko', f./1e6, cp_kk, 'k-');
set(gca, 'XLim', [0 f_max]);
box on;
xlabel('Frequency [MHz]');
ylabel('C_p [m/s]');