% Sensor Element Directivity in 2D Example
%
% This example demonstrates how a single-element sensor, or the individual
% elements of a sensor array, can be made directional, ie. the sensitivity
% depends on the angle of incidence of the incoming wave.
%
% author: Ben Cox
% date: 20th January 2010
% last update: 1st February 2011
%  
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009, 2010, 2011 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

% =========================================================================
% DEFINE THE GRID AND MEDIUM PROPERTIES
% =========================================================================

% create the computational grid
Nx = 64;        % number of pixels in the x (column) direction
Nz = 64;       % number of pixels in the z (row) direction
dx = 1e-3/Nx;   % pixel width [m]
dz = dx/2;     	% pixel height [m]
kgrid = makeGrid(Nx, dx, Nz, dz);

% define the properties of the propagation medium
medium.sound_speed = 1500;  % [m/s]

% create the time array, then halve it to make the example end sooner
[kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed);
kgrid.t_array = (0:0.5*length(kgrid.t_array))*dt;

% =========================================================================
% DEFINE THE DIRECTIONAL SENSOR ARRAY
% =========================================================================

% define a line of single-element sensors
sensor.mask = zeros(Nz,Nx);
sensor.mask(24,2:2:63) = 1;

% define the angle of max directivity for each single-element sensor 
% 0 = max sensitivity in z direction (up/down)
% pi/2 or -pi/2 = max sensitivity in x direction (left/right)
sensor.directivity_angle = zeros(Nz,Nx);
sensor.directivity_angle(24,2:2:63) = (-1:1/15:1)*pi/2;

% define the directivity size
sensor.directivity_size = 16*kgrid.dx;

% =========================================================================
% SIMULATION AND VISUALISATION FOR AN INITIAL VALUE PROBLEM
% =========================================================================

% define the initial pressure distribution
source.p0 = zeros(Nz,Nx);
source.p0(39:41 ,:) = 2;

% run the simulation
[sensor_data p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor, 'PMLAlpha', [0 4]);

% plot the largest value of the output for each sensor
figure, plot(((-1:1/15:1)*pi/2),max(sensor_data,[],2),'o')
xlabel('sensor directivity angle (radians)')
ylabel('maxima of single-element sensors'' outputs')

% =========================================================================
% SIMULATION AND VISUALISATION FOR TIME-VARYING SOURCE PROBLEM
% =========================================================================

% define a time varying sinusoidal source instead of an initial pressure
source.p_mask = source.p0;
source.p0 = [];
source_freq = 12e6;
source_mag = 1;
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);
    
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, 'PMLAlpha', [0 4]);

% plot the largest value of the output for each sensor
figure, plot(((-1:1/15:1)*pi/2), max(sensor_data(:,100:length(kgrid.t_array)),[],2),'o')
xlabel('sensor directivity angle (radians)')
ylabel('maxima of single-element sensors'' outputs')
axis([-pi/2 pi/2 0 1.02])
figure, polar(((-1:1/15:1)'*pi/2), max(sensor_data(:,100:length(kgrid.t_array)),[],2))

