% Simulations In Three Dimensions Example
%
% This example provides a demonstration of using k-Wave for the simulation
% and detection of photoacoustic waves within a three-dimensional
% heterogeneous propagation medium. It builds on the Homogeneous
% Propagation Medium, Heterogeneous Propagation Medium, and Controlling The
% Absorbing Boundary Layer examples.
%
% author: Bradley Treeby
% date: 1st July 2009
% last update: 16th July 2009
%  
% This example is part of the k-Wave Toolbox (http://www.k-wave.org)

clear all;

% =========================================================================
% SIMULATION
% =========================================================================

% create the computational grid
Nx = 64;            % number of pixels in the x direction
Ny = 64;            % number of pixels in the y direction
Nz = 64;            % number of pixels in the z direction
dx = 0.1e-3;        % pixel width [m]
dy = 0.1e-3;        % pixel width [m]
dz = 0.1e-3;        % pixel height [m]
kgrid = makeGrid(Nx, dx, Ny, dy, Nz, dz);

% define the properties of the propagation medium
c0 = 1500;                      % [m/s]
c = c0*ones(Nz, Nx, Ny);        % [m/s]
c(1:Nz/2, :, :) = c0*1.2;

rho0 = 1000;                    % [kg/m^3]
rho = rho0*ones(Nz, Nx, Ny);    % [kg/m^3]
rho(:, Nx/4:end, :) = rho0*1.2;

% create initial pressure distribution using makeBall
ball_magnitude = 10;
ball_x_pos = 36;    % pixels
ball_y_pos = 36;    % pixels
ball_z_pos = 36;    % pixels
ball_radius = 5;    % pixels
ball_1 = ball_magnitude*makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);

ball_magnitude = 10;
ball_x_pos = 20;    % pixels
ball_y_pos = 20;    % pixels
ball_z_pos = 20;    % pixels
ball_radius = 3;    % pixels
ball_2 = ball_magnitude*makeBall(Nx, Ny, Nz, ball_x_pos, ball_y_pos, ball_z_pos, ball_radius);

p0 = ball_1 + ball_2;

% define a series of Cartesian points to collect the data
z = (-22:2:22)*dz;          % [m]
x = (-22:2:22)*dx;          % [m]
y = 22*dy*ones(size(z));    % [m]
sensor_mask = [x; y; z];

% input arguments
input_args = {'PlotLayout', true, 'PMLSize', 10};

% run the simulation
sensor_data = kspaceFirstOrder3D(p0, kgrid, c, rho, 'auto', sensor_mask, input_args{:});

% =========================================================================
% VISUALISATION
% =========================================================================

% plot the simulated sensor data
figure;
imagesc(sensor_data, [-1, 1]);
colormap(getColorMap);
ylabel('Sensor Position');
xlabel('Time Step');
colorbar;