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{:});