% Modelling a Focussed Detector Example
%
% This example demonstrates how to model a focussed detector, when the
% focussing arises due to spatial averaging over the detector's surface.
% It uses photoacoustic sources, but the focussing principle is general and
% will hold for all source types.
%
% author: Ben Cox
% date: 20th January 2010
% last update: 22nd January 2010
%  
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009, 2010 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
% 
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details. 
% 
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>. 

clear all;

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

% create the computational grid
Nx = 180;           % number of pixels in the x (column) direction
Nz = 180;           % 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
medium.sound_speed = 1500;    % [m/s]
medium.density = 1000;        % [kg/m^3]

% define a sensor as part of a circle centred on the grid
sensor_radius = 65; % [pixels]
arc_angle = pi;     % [rad]
sensor.mask = makeCircle(Nx, Nz, Nx/2, Nz/2, sensor_radius, arc_angle);

% define the array of temporal points [s]
[kgrid.t_array, dt] = makeTime(kgrid, medium.sound_speed);
kgrid.t_array = (0:round(length(kgrid.t_array)*0.75))*dt;

% place a disc-shaped source near the focus of the detector
source.p0 = makeDisc(Nx, Nz, Nx/2, Nz/2, 4)*2;

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

% place a disc-shaped source away from the focus of the detector
source.p0 = makeDisc(Nx, Nz, Nx/2+20, Nz/2, 4)*2;

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

% average the outputs of the simulations over the points making up the sensor
sensor_output1 = sum(sensor_data1,1)/sum(sensor.mask(:));
sensor_output2 = sum(sensor_data2,1)/sum(sensor.mask(:));

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

% plot the recorded time signals
figure
plot(kgrid.t_array, sensor_output1, 'k', kgrid.t_array, sensor_output2, 'r')
xlabel('Time (s)')
ylabel('Pressure measured over semi-circular detector') 
legend('Source at focus, sum(output^2) = 0.57', 'Source off focus, sum(output^2) = 0.04')
