% Example script for simulation of time series data using k-Wave
% This script is compiled from example snippets from the official
% documentation of the k-Wave toolbox.
% Please look at the k-Wave most recent documentation to make sure that 
% you are using k-Wave correctly.

% create the computational grid
Nx = 256;           % number of grid points in the x (row) direction
Ny = 256;           % number of grid points in the y (column) direction
dx = 0.5e-3;        % grid point spacing in the x direction [m]
dy = 0.5e-3;        % grid point spacing in the y direction [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]
medium.alpha_coeff = 0.75;  % [dB/(MHz^y cm)]
medium.alpha_power = 1.5;

% create initial pressure distribution using makeDisc
disc_magnitude = 5; % [Pa]
disc_x_pos = 100;    % [grid points]
disc_y_pos = 100;    % [grid points]
disc_radius = 8;    % [grid points]
disc_1 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_magnitude = 3; % [Pa]
disc_x_pos = 160;    % [grid points]
disc_y_pos = 120;    % [grid points]
disc_radius = 4;    % [grid points]
disc_2 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

disc_magnitude = 18; % [Pa]
disc_x_pos = 130;    % [grid points]
disc_y_pos = 150;    % [grid points]
disc_radius = 2;    % [grid points]
disc_3 = disc_magnitude * makeDisc(Nx, Ny, disc_x_pos, disc_y_pos, disc_radius);

source.p0 = disc_1 + disc_2 + disc_3;

% define a centered circular sensor
num_sensor_points = 50;
mask = zeros(2, num_sensor_points);
for i = 1:num_sensor_points
   mask(1, i) = -0.04+i*(0.08/num_sensor_points);
   mask(2, i) = -0.04 + 0.005;
end
sensor.mask = mask;

% define the frequency response of the sensor elements
center_freq = 3e6;      % [Hz]
bandwidth = 200;         % [%]
sensor.frequency_response = [center_freq, bandwidth];

% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);

% Saving the simulation in the IPASC format using the PACFISH MATLAB tool
% Please make sure to take a look at the most up-to-date examples on the
% PACFISH Github repository (https://github.com/IPASC/PACFISH/) to make
% sure that the code to export a k-Wave simulation has not changed.
pf_adapter = pacfish.kwave_adapter(sensor, sensor_data, medium, kgrid, [-0.03, 0.03, 0, 0, -0.03, 0.03]);
pa_data = pf_adapter.get_pa_data();
pacfish.write_data("kwave_2Dsim_linear_array.hdf5", pa_data, 1);

% plot the simulated sensor data (as a sanity tst that the simulations were
% successful)
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;