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);