% Comparison between the on-axis pressure from a circular piston source
% simulated using k-Wave with the analytical expression from [1]. The
% errors in the k-Wave simulation arise due to the stair-cased
% representation of the circular source. These decrease as the scale
% parameter is increased.
%
% [1] A. D. Pierce, Acoustics: An Introduction to its Physical Principles
% and Applications. New York: Acoustical Society of America, 1989.
%
% author: Bradley Treeby
% date: 21 March 2014
% last update: 20th November 2014
%       
% This function is part of the k-Wave Toolbox (https://bug.medphys.ucl.ac.uk/kwave)
% Copyright (C) 2009-2014 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 LITERALS
% =========================================================================
    
% select which code to run
%   1: MATLAB CPU code
%   2: MATLAB GPU code
%   3: C++ code
MODEL           = 3;

% scale parameter (changes the resolution of the simulation)
SC              = 1;

% grid parameters
PML_SIZE        = 10;                    % size of the perfectly matched layer [grid points]
Nx              = 256*SC - 2*PML_SIZE;   % number of grid points in the x direction
Ny              = 128*SC - 2*PML_SIZE;   % number of grid points in the y direction
Nz              = 128*SC - 2*PML_SIZE;   % number of grid points in the z direction

% sampling parameters
PPW             = 12*SC;    % points per wavelength (this defines the grid size)
PPP             = 50*SC;    % points per period (this defines the CFL)
T_END           = 40e-6;    % total compute time [s] (this must be long enough to reach steady state)
RECORD_PERIODS  = 2;        % number of periods to record

% medium parameters
C0              = 1500;     % sound speed [m/s]
RHO0            = 1000;     % density [kg/m^3]

% source parameters
F0              = 1e6;      % source frequency [Hz]
SOURCE_RADIUS   = 40*SC;    % piston radius [grid points]
SOURCE_MAG      = 0.5e6;    % source pressure [Pa]

% =========================================================================
% CREATE GRID
% =========================================================================

% calculate the grid spacing based on the PPW and F0
dx = C0 / (PPW * F0);   % [m]
dt = 1 / (PPP*F0);      % [s]

% calculate the CFL
disp(['CFL = ' num2str(C0*dt/dx)]);

% create the computational grid
kgrid = makeGrid(Nx, dx, Ny, dx, Nz, dx);

% create the time array
kgrid.t_array = 0:dt:T_END;

% assign medium properties
medium.sound_speed = C0;
medium.density = RHO0;

% =========================================================================
% CREATE SOURCE
% =========================================================================

% create piston source mask
piston = makeDisc(Ny, Nz, Ny/2, Nz/2, SOURCE_RADIUS);
source.p_mask = zeros(Nx, Ny, Nz);
source.p_mask(1, :, :) = piston;

% create time varying source
source.p = SOURCE_MAG * sin(2*pi*F0*kgrid.t_array);

% filter source with up-ramp
ramp_length = round((2*pi/F0)/dt);
ramp = (-cos( (0:(ramp_length-1))*pi/ramp_length ) + 1)/2;
source.p(1:ramp_length) = ramp .* source.p(1:ramp_length);

% =========================================================================
% CREATE SENSOR MASK
% =========================================================================

% set sensor mask to record central axis, not including the source point
sensor.mask = zeros(Nx, Ny, Nz);
sensor.mask(2:end, Ny/2, Nz/2) = 1;

% record the maximum pressure
sensor.record = {'p_max'};

% average only the final few periods when the field is in steady state
sensor.record_start_index = kgrid.Nt - RECORD_PERIODS*PPP + 1;

% =========================================================================
% RUN k-WAVE SIMULATION
% =========================================================================

% run code
switch MODEL
    case 1
        % MATLAB CPU
        sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...
            'DataCast', 'single', 'PMLInside', false, ...
            'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG);
    case 2
        % MATLAB GPU
        sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...
            'DataCast', 'gpuArray-single', 'PMLInside', false, ...
            'DisplayMask', source.p_mask, 'PlotScale', [-1, 1]*SOURCE_MAG);
    case 3
        % C++
        sensor_data = kspaceFirstOrder3DC(kgrid, medium, source, sensor, 'PMLInside', false);
end

% =========================================================================
% ANALYTICAL SOLUTION
% =========================================================================

% calculate the wavenumber
k = 2*pi*F0./C0;

% define radius and axis
a = SOURCE_RADIUS*dx;
x = (kgrid.x_vec(2:end, :) - min(kgrid.x_vec(:)));

% calculate the analytical solution for a piston in an infinite baffle
% for comparison (Eq 5-7.3 in Pierce)
r_a = sqrt(x.^2 + a^2);
p_ref = SOURCE_MAG*abs(-2i*exp(1i*k*(x+r_a)/2).*sin((k*r_a - k*x)/2));

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

% plot the pressure along the focal axis of the piston
figure;
plot(x/a, p_ref, 'k-');
hold on;
plot(x/a, sensor_data.p_max, 'bx')
xlabel('Normalised Distance');
legend('Exact', 'k-Wave');