I am trying to get the heating by ultrasound example working but in a 2D heterogeneous medium., but it does not work. I think the problem is that defining time arrey (kgrid.setTime) depends on Nt and dt and they also depend on ppp which is a matrix and not a single integer anymore. Is there a way to fix this? I appreciate any hint. Here is the code I use but does not work:
clearvars;
% =========================================================================
% ACOUSTIC SIMULATION
% =========================================================================
% define the PML size
pml_size = 20; % [grid points]
% define the grid parameters
Nx = 256 - 2 * pml_size; % [grid points]
Ny = 256 - 2 * pml_size; % [grid points]
dx = 0.25e-3; % [m]
dy = 0.25e-3; % [m]
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dy);
% define the properties of the propagation medium
medium.alpha_coeff = 0.75; % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;
medium.sound_speed = 1500 * ones(Nx, Ny); % [m/s]
medium.sound_speed(1:Nx/2, :) = 1800; % [m/s]
medium.density = 1000 * ones(Nx, Ny); % [kg/m^3]
medium.density(:, Ny/4:Ny) = 1200; % [kg/m^3]
% define the source parameters
diameter = 45e-3; % [m]
radius = 35e-3; % [m]
freq = 1e6; % [Hz]
amp = 0.5e6; % [Pa]
% define a focused ultrasound transducer
source.p_mask = makeArc([Nx, Ny], [1, Ny/2], round(radius / dx), round(diameter / dx) + 1, [Nx/2, Ny/2]);
% calculate the time step using an integer number of points per period
ppw = 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 ) / 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);
sensor.record = {'p', 'p_max_all'};
% 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 = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});