Hi,
First of all thank you for building this amazing toolbox.
I have simulated the bowel transducer in 3D and then compared the focus pressure amplitude with the analytical model of a concave spherical transducer (from Theory of Focusing Radiators by H.T. O'Neil). I also checked this analytical model results from the Sonic Concepts transducer (H-276) and it was correct. I couldn't find the mistake in my code I think I have done everything correctly but the results shouldn't deviate that much (like the peak pressure at the focus must be 58 MHz but the simulated pressure is 46.8 MHz ). I have posted the code below:
(`% define the grid parameters 3D
Nx = 150; % [grid points]
Ny = 150; % [grid points]
Nz = 92;
dx = 0.1e-3; % [m]
dy = 0.1e-3; % [m]
dz = 0.1e-3; % [m]
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
% define the properties of the propagation medium water
medium.sound_speed = 1500; % [m/s]
medium.density = 1000; % [kg/m^3]
medium.alpha_coeff = 0.025; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
% define parameters
grid_size = [150,150,92];
bowl_pos = [Nx/2, Ny/2, 2];
diameter = 15e-3; % [m]
radius = 8.5e-3; % [m]
focus_pos = [Nx/2, Ny/2, 87];
freq = 1.5e6; % [Hz]
amp = 2.06e6; % [Pa]
% create bowl
source.p_mask = makeBowl(grid_size, bowl_pos, round(radius / dx), round(diameter / dx)+1, focus_pos, 'Plot', true);
% calculate the time step using an integer number of points per period
ppw = max(medium.sound_speed(:)) / (freq * dx); % points per wavelength
cfl = 0.3; % cfl number
ppp = ceil(ppw / cfl); % points per period
T = 1 / freq; % period [s]
dt = T / ppp; % time step [s]
% calculate the number of time steps to reach steady state
t_end = sqrt( kgrid.x_size.^2 + kgrid.y_size.^2 + kgrid.z_size.^2 ) / max(medium.sound_speed(:));
Nt = round(t_end / dt);
% create the time array
kgrid.setTime(Nt, dt);
% define the input signal
source.p = createCWSignals(kgrid.t_array, freq, amp, 0);
% set the sensor mask to cover the entire grid
sensor.mask = ones(Nx, Ny, Nz);
sensor.record = {'p','p_rms','u','u_rms', 'I','I_avg'};
% record the last 3 cycles in steady state
num_periods = 3;
T_points = round(num_periods * T / kgrid.dt);
sensor.record_start_index = Nt - T_points + 1;
% set the input arguements
input_args = {'PMLInside', false, 'PlotPML', false, 'DisplayMask', ...
'off', 'PlotScale', [-1, 1] * amp};
% run the acoustic simulation
sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, input_args{:});
% convert the absorption coefficient to nepers/m
alpha_np = db2neper(medium.alpha_coeff, medium.alpha_power) * ...
(2 * pi * freq).^medium.alpha_power;
% extract the pressure amplitude at each position
p = extractAmpPhase(sensor_data.p, 1/kgrid.dt, freq);
% reshape the data, and calculate the volume rate of heat deposition
p = reshape(p, Nx, Ny, Nz);