% Using A Binary Sensor Mask Example
%
% This example demonstrates how to use a binary sensor mask for the
% detection of photoacoustic waves within a two-dimensional homogeneous
% propagation medium. It builds on the Homogeneous Propagation Medium
% Example.
%
% author: Bradley Treeby
% date: 29th June 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 = 128;           % number of pixels in the x (column) direction
Nz = 128;           % number of pixels in the z (row) direction
dx = 0.1e-3;        % pixel width [m]
dz = 0.1e-3;        % pixel height [m]
kgrid = makeGrid(Nx, dx, Nz, dz);

% define the properties of the propagation medium
c = 1500;           % [m/s]
rho = 1000;         % [kg/m^3]

% create initial pressure distribution using makeDisc
disc_magnitude = 5;
disc_x_pos = 50;    % pixels
disc_z_pos = 50;    % pixels
disc_radius = 8;    % pixels
disc_1 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius);

disc_magnitude = 3;
disc_x_pos = 60;    % pixels
disc_z_pos = 80;    % pixels
disc_radius = 5;    % pixels
disc_2 = disc_magnitude*makeDisc(Nx, Nz, disc_x_pos, disc_z_pos, disc_radius);

p0 = disc_1 + disc_2;

% define a binary sensor mask
sensor_x_pos = Nx/2;        % pixels
sensor_z_pos = Nz/2;        % pixels
sensor_radius = Nx/2 - 20;  % pixels
sensor_arc_angle = 3*pi/2;  % radians
sensor_mask = makeCircle(Nx, Nz, Nx/2, Nz/2, sensor_radius, sensor_arc_angle);

% run the simulation
sensor_data = kspaceFirstOrder2D(p0, kgrid, c, rho, 'auto', sensor_mask);

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

% plot the initial pressure and sensor distribution
figure;
imagesc(kgrid.x(1,:)*1e3, kgrid.z(:,1)*1e3, p0 + disc_magnitude*sensor_mask, [-1 1]);
colormap(getColorMap);
ylabel('z-position [mm]');
xlabel('x-position [mm]');
axis image;

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