Hi,
Is it possible in k-wave to calculate the sound pressure distribution after N time steps by running simulations in a loop after 1 time step? The loop runs kspaceFirstOrder where at each grid node the wave source is the value from the previous time step in the loop. It seems to me that this should work, but maybe I'm missing something and it's not working. I have attached my code below.
Kind regards
pml_size = 20; % [grid points]
% define the grid parameters
Nx = 512 - 2 * pml_size; % [grid points]
Ny = 512 - 2 * pml_size; % [grid points]
dx = 10e-6; % [m]
dy = 10e-6; % [m]
% create the computational grid
kgrid = kWaveGrid(Nx, dx, Ny, dy);
%%
% define the properties of the propagation medium
medium.sound_speed = 1480; % [m/s]
medium.density = 1000; % [kg/m^3]
medium.alpha_coeff = 0.02; % [dB/(MHz^y cm)]
medium.alpha_power = 1.1;
% define the source parameters
diameter = 4e-3; % [m]
radius = 3e-3; % [m]
freq = 1.05e6; % [Hz]
amp = 0.1e6; % [Pa]
mask= logical(zeros(Nx,Ny));
mask(Nx/2,Ny/2) = 1;
source.p_mask = zeros(Nx,Ny);
source.p_mask(mask) = 1;
% 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);
kgrid.setTime(Nt, dt);
% define the input signal
input_signal = createCWSignals(kgrid.t_array, freq, amp, 0);
% create the time array
kgrid.setTime(1, dt);
source.p = zeros(1,1);
source.p = input_signal(:,1);
% set the sensor mask to cover the entire grid
sensor.mask = ones(Nx, Ny);
sensor.record = {'p'};
% record one point
sensor.record_start_index = 1;
% set the input arguements
input_args = {'PMLInside', false, 'PlotPML', false, 'DisplayMask', ...
'off', 'PlotScale', [-1, 1] * amp, ...
'PlotSim', false, 'DataCast', 'gpuArray-single'};
% run the acoustic simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% % reshape the data
p = reshape(sensor_data.p, Nx, Ny,1);
%%
source.p_mask = ones(Nx,Ny);
f = waitbar(0);
for i=2:length(input_signal)
source.p = p(:,:,i-1);
source.p = source.p(:);
source.p(mask(:)) = input_signal(:,i);
% run the acoustic simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:});
% % reshape the data
p (:,:,i) = reshape(sensor_data.p, Nx, Ny,1);
waitbar(i/(length(input_signal)),f)
end
close(f)