Dear Xiong,

Your question has made me think of a few points. I hope they’re helpful.

1) k-Wave will apply the time integration that you mention automatically. It will multiply each sample of the time-varying source.p by the timestep dt. I think you therefore just need to write:

source.p = Gruneisen coefficient * H1;

where this H1 is the sampled version of your known time varying H1(t). See comment (3) below about checking the amplitude is correct.

2) Just to clear up a possible confusion. I said in my last post that "if your EM pulse is Gaussian, then the pressure source will be the time derivative of a Gaussian". That is not true for k-Wave. The source term for the underlying wave equation is the time derivative, but that is taken care of in the code. The input source.p should be the actual shape of the EM pulse you are simulating. For instance, take the example "Point Source In A Homogeneous Propagation Medium Example" as a starting point, comment out the definition of the sinusoidal source and replace it with a Gaussian. Note the time offset.

% define a time varying photoacoustic Gaussian pulse

% with the Gaussian normalised so its integral sum(gaussian)*dt = 1

source_mag = 1;

gruneisen = 0.11;

pulsewidth = 1e-6;

t_offset = 5e-6;

gaussian = exp(-((kgrid.t_array-t_offset)/pulsewidth).^2)/(pulsewidth*sqrt(pi));

source.p = source_mag*gruneisen*gaussian;

Also adjust the following line to make sure the plot scaling is sensible:

[sensor_data p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor,... 'PlotScale',[-1e4 1e4]);

Just to be tidy, you might also want to delete the instructions that "plot the final wave-field".

3) It might be a good idea to check that you are getting the amplitude you expect by putting in a plane source (so there’s no geometrical spreading). Adjusting the above example a little:

% define a line of source elements to generate a plane wave

source.p_mask = zeros(Nz, Nx);

source.p_mask(Nz/2, :) = 1;

Also set PMLAlpha to turn off the PML at the sides so the plane wave remains truly plane:

[sensor_data p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor,... 'PlotScale',[-7e4 7e4],'PMLAlpha',[0 4]);

4) Final thing. There’s a bug, which we hope to fix for the next release, which means that if you want quantitatively accurate results when you are using a time-varying source you should set dx = dz to be equal.

Kind regards,

Ben