Dear,
I really appreciate you've distributed this good and useful simulation tool for everyone.
It seems that K-wave tool misses the scaling factor that must be considered between discrete and continuous Fast Fourier transform in Matlab.
Please run the following simple K-wave code with n=128 and n=256. (or you could try n=64 and n=128. However, for this setup, the difference of the ball surface structure between n=64 and 128 is intensified, which causes additional difference in results.)
Since all grid values in this simulation code are set by given physical quantities, different n does not affect simulation result much.
You can verify that the magnitude of the simulated photoacoustic(PA) signals of n=128 and 256 is quite different each other. The PA magnitude of n=256 is almost 4 times larger than that of n=128. As mentioned above, the PA magnitude must be almost the same. (Although there might be slight difference due to grid size difference, for example, maximum frequency values and so on)
Since Dt(=kgrid.t_array(2)-kgrid.t_array(1)) of n=256 case is twice times smaller than Dt of n=128, the PA signal magnitudes are the same if you consider PA*Dt^2.
Related to this PA signal magnitude discrepancy issue, I wonder how much the absolute values of PA signals simulated by K-wave are reliable. We cannot simply consider PA*Dt^2 that is far different from the real PA signal value. I think some factor containing Dt must be missed on converting your theoretical equations to K-wave code.
-------------------------code------------------------------------
WinSize = 30e-3; % [m] physical length of the 3D cube where PA imaging happens
% Ball
B_radius = 2e-3; %[m]
% create the computational grid
n = 256; % grid number in WinSize
Nx = n; % number of grid points in the x direction
Ny = n; % number of grid points in the y direction
Nz = n; % number of grid points in the z direction
dx = WinSize/Nx; % grid point spacing in the x direction [m]
dy = dx; % grid point spacing in the y direction [m]
dz = dx; % grid point spacing in the z direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
medium.sound_speed = 1500; % ultrasound speed in soft tissue[m/s]
kgrid.makeTime(medium.sound_speed);
% create a concave sensor
sphere_offset = 10; % [grid points]
f0 = 0.02; % focal length of the focused ultrasound transducer
NA = 0.375;
RadiusofC = ceil(f0/dx); % [grid points; raduis of curvarture = focal length]
diameter = 2*ceil(NA*RadiusofC)+1; % [grid points; must be an odd number]
bowl_pos = [1 + sphere_offset, Ny/2, Nz/2]; % [grid points]
focus_pos = [Nx/2, Ny/2, Nz/2]; % [grid points; focal axis direction: from bowl_pos to focus_pos]
sensor.mask = makeBowl([Nx, Ny, Nz], bowl_pos, RadiusofC, diameter, focus_pos);
% create a ball PA object with exponentially attenuated irradiance
ball_radius = ceil(B_radius/dx);
ball_pos_x = 1 + sphere_offset + RadiusofC + ball_radius;
ball_pos_y = Ny/2;
ball_pos_z = Nz/2;
source.p0 = makeBall(Nx, Ny, Nz, ball_pos_x, ball_pos_y, ball_pos_z, ball_radius);
sensor_data = kspaceFirstOrder3DG(kgrid, medium, source, sensor);
sensor_data1 = sum(sensor_data, 1);
% figure
plot(kgrid.t_array, sensor_data1, '-');
----------------------------------------------------------