function [sensor_data, field_data] = kspaceFirstOrder3D(kgrid, medium, source, sensor, varargin)
%KSPACEFIRSTORDER3D     3D time-domain simulation of wave propagation.
%
% DESCRIPTION:
%       kspaceFirstOrder3D simulates the time-domain propagation of linear
%       compressional waves through a three-dimensional homogeneous or
%       heterogeneous acoustic medium given four input structures: kgrid,
%       medium, source, and sensor. The computation is based on a
%       first-order k-space model which allows power law absorption and a
%       heterogeneous sound speed and density. At each time-step (defined
%       by kgrid.t_array), the pressure at the positions defined by
%       sensor.mask are recorded and stored. If kgrid.t_array is set to
%       'auto', this array is automatically generated using makeTime. An
%       absorbing boundary condition called a perfectly matched layer (PML)
%       is implemented to prevent waves that leave one side of the domain
%       being reintroduced from the opposite side (a consequence of using
%       the FFT to compute the spatial derivatives in the wave equation).
%       This allows infinite domain simulations to be computed using small
%       computational grids.
%
%       For a homogeneous medium the formulation is exact and the
%       time-steps are only limited by the effectiveness of the perfectly
%       matched layer. For a heterogeneous medium, the solution represents
%       a leap-frog pseudospectral method with a Laplacian correction that
%       improves the accuracy of computing the temporal derivatives. This
%       allows larger time-steps to be taken without instability compared
%       to conventional pseudospectral time-domain methods. The
%       computational grids are staggered both spatially and temporally.
%
%       An initial pressure distribution can be specified by assigning a
%       matrix (the same size as the computational grid) of arbitrary
%       numeric values to source.p0. A time varying pressure source can
%       similarly be specified by assigning a binary matrix (i.e., a matrix
%       of 1's and 0's with the same dimensions as the computational grid)
%       to source.p_mask where the 1's represent the pixels that form part
%       of the source. The time varying input signals are then assigned to
%       source.p. This must be the same length as kgrid.t_array and can be
%       a single time series (in which case it is applied to all source
%       elements), or a matrix of time series following the source elements
%       using MATLAB's standard column-wise linear matrix index ordering. A
%       time varying velocity source can be specified in an analogous
%       fashion, where the source location is specified by source.u_mask,
%       and the time varying input velocity is assigned to source.ux,
%       source.uy, and source.uz.
%
%       The pressure is returned as an array of time series at the sensor
%       locations defined by sensor.mask. This can be given either as a
%       binary grid (i.e., a matrix of 1's and 0's with the same dimensions
%       as the computational grid) representing the pixels within the
%       computational grid that will collect the data, or as a series of
%       arbitrary Cartesian coordinates within the grid at which the
%       pressure values are calculated at each time step via interpolation.
%       The Cartesian points must be given as a 3 by N matrix corresponding
%       to the x, y and z positions, respectively. The final pressure field
%       over the complete computational grid can also be obtained using the
%       output field_data. If no output is required, the sensor input can
%       be replaced with an empty array [].
%
%       If sensor.mask is given as a set of Cartesian coordinates, the
%       computed sensor_data is returned in the same order. If sensor.mask
%       is given as a binary grid, sensor_data is returned using MATLAB's
%       standard column-wise linear matrix index ordering. In both cases,
%       the recorded data is indexed as sensor_data(sensor_position, time).
%       For a binary sensor mask, the pressure values at a particular time
%       can be restored to the sensor positions within the computation grid
%       using unmaskSensorData.
%
%       By default, the recorded pressure field is passed directly to the
%       output arguments sensor_data and field_data. However, the particle
%       velocity can also be recorded by setting the optional input
%       parameter 'ReturnVelocity' to true. In this case, the output
%       arguments sensor_data and field_data are returned as structures
%       with the pressure and particle velocity appended as separate
%       fields. In two dimensions, these fields are given by sensor_data.p,
%       sensor_data.ux, sensor_data.uy, and sensor_data.uz, and
%       field_final.p, field_final.ux, field_final.uy, and field_final.uz.
%
%       kspaceFirstOrder3D may also be used for time reversal image
%       reconstruction by assigning the time varying pressure recorded over
%       an arbitrary sensor surface to the input field
%       sensor.time_reversal_boundary_data. This data is then enforced in
%       time reversed order as a time varying Dirichlet boundary condition
%       over the sensor surface given by sensor.mask. The boundary data
%       must be indexed as
%       sensor.time_reversal_boundary_data(sensor_position, time). If
%       sensor.mask is given as a set of Cartesian coordinates, the
%       boundary data must be given in the same order. An equivalent binary
%       sensor mask (computed using nearest neighbour interpolation) is
%       then used to place the pressure values into the computational grid
%       at each time step. If sensor.mask is given as a binary grid of
%       sensor points, the boundary data must be ordered using MATLAB's
%       standard column-wise linear matrix indexing. If no additional
%       inputs are required, the source input can be replaced with an empty
%       array [].
%
%       Acoustic attenuation compensation can also be included during time
%       reversal image reconstruction by assigning the absorption
%       parameters medium.alpha_coeff and medium.alpha_power and reversing
%       the sign of the absorption term by setting medium.alpha_sign = [-1,
%       1]. This forces the propagating waves to grow according to the
%       absorption parameters instead of decay. The reconstruction should
%       then be regularised by assigning a filter to medium.alpha_filter
%       (this can be created using getAlphaFilter).
%
%       Note: To run a simple reconstruction example using time reversal
%       (that commits the 'inverse crime' of using the same numerical
%       parameters and model for data simulation and image reconstruction),
%       the sensor_data returned from a k-Wave simulation can be passed
%       directly to sensor.time_reversal_boundary_data  with the input
%       fields source.p0 and source.p removed or set to zero.
%
% USAGE:
%       sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor)
%       sensor_data = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...) 
%
%       [sensor_data, field_data] = kspaceFirstOrder3D(kgrid, medium, source, sensor)
%       [sensor_data, field_data] = kspaceFirstOrder3D(kgrid, medium, source, sensor, ...) 
%
%       kspaceFirstOrder3D(kgrid, medium, source, [])
%       kspaceFirstOrder3D(kgrid, medium, source, [], ...) 
%
% INPUTS:
% The minimum fields that must be assigned to run an initial value problem
% (for example, a photoacoustic forward simulation) are marked with a *. 
%
%       kgrid*              - k-space grid structure returned by makeGrid
%                             containing Cartesian and k-space grid fields  
%       kgrid.t_array*      - evenly spaced array of time values [s] (set
%                             to 'auto' by makeGrid) 
%
%
%       medium.sound_speed* - sound speed distribution within the acoustic
%                             medium [m/s] 
%       medium.density*     - density distribution within the acoustic
%                             medium [kg/m^3] 
%       medium.alpha_power  - power law absorption exponent
%       medium.alpha_coeff  - power law absorption coefficient 
%                             [dB/(MHz^y cm)] 
%       medium.alpha_mode   - optional input to force either the absorption
%                             or dispersion terms in the equation of state
%                             to be excluded; valid inputs are
%                             'no_absorption' or 'no_dispersion' 
%       medium.alpha_filter - frequency domain filter applied to the
%                             absorption and dispersion terms in the
%                             equation of state 
%       medium.alpha_sign   - two element array used to control the sign of
%                             absorption and dispersion terms in the
%                             equation of state  
%
%
%       source.p0*          - initial pressure within the acoustic medium
%       source.p            - time varying pressure at each of the source
%                             positions given by source.p_mask 
%       source.p_mask       - binary grid specifying the positions of the
%                             time varying pressure source distribution
%       source.ux           - time varying particle velocity in the
%                             x-direction at each of the source positions
%                             given by source.u_mask 
%       source.uy           - time varying particle velocity in the
%                             y-direction at each of the source positions
%                             given by source.u_mask 
%       source.uz           - time varying particle velocity in the
%                             z-direction at each of the source positions
%                             given by source.u_mask  
%       source.u_mask       - binary grid specifying the positions of the
%                             time varying particle velocity distribution 
%
%
%       sensor.mask*        - binary grid or a set of Cartesian points
%                             where the pressure is recorded at each
%                             time-step  
%       sensor.time_reversal_boundary_data - time varying pressure
%                             enforced as a Dirichlet boundary condition
%                             over sensor.mask  
%       sensor.frequency response - two element array specifying the center
%                             frequency and percentage bandwidth of a
%                             frequency domain Gaussian filter applied to
%                             the sensor_data    
%
% Note: For heterogeneous medium parameters, medium.sound_speed and
% medium.density must be given in matrix form with the same dimensions as
% kgrid. For homogeneous medium parameters, these can be given as single
% numeric values. If the medium is homogeneous and velocity inputs or
% outputs are not required, it is not necessary to specify medium.density.
%
% OPTIONAL INPUTS:
%       Optional 'string', value pairs that may be used to modify the
%       default computational settings.
%
%       'CartInterp'- Interpolation mode used to extract the pressure when
%                     a Cartesian sensor mask is given. If set to 'nearest'
%                     and more than one Cartesian point maps to the same
%                     pixel, duplicated data points are discarded and
%                     sensor_data will be returned with less points than
%                     that specified by sensor.mask (default = 'nearest').
%       'CreateLog' - Boolean controlling whether the command line output
%                     is saved using the diary function with a date and
%                     time stamped filename (default = false). 
%       'DataCast'  - String input of the data type that variables are cast
%                     to before computation. For example, setting to
%                     'single' will speed up the computation time (due to
%                     the improved efficiency of fftn and ifftn for this
%                     data type) at the expense of a loss in precision.
%                     This variable is also useful for utilising GPU
%                     parallelisation through libraries such as GPUmat or
%                     AccelerEyesJacket by setting 'DataCast' to
%                     'GPUsingle' or 'gsingle' (default = 'off').
%       'DisplayMask' - Binary matrix overlayed onto the animated
%                     simulation display. Elements set to 1 within the
%                     display mask are set to black within the display
%                     (default = sensor.mask).
%       'LogScale'  - Boolean controlling whether the pressure field is log
%                     compressed before display (default = false). The data
%                     is compressed by scaling both the positive and
%                     negative values between 0 and 1 (truncating the data
%                     to the given plot scale), adding a scalar value
%                     (compression factor) and then using the corresponding
%                     portion of a log10 plot for the compression (the
%                     negative parts are remapped to be negative thus the
%                     default color scale will appear unchanged). The
%                     amount of compression can be controlled by adjusting
%                     the compression factor which can be given in place of
%                     the Boolean input. The closer the compression factor
%                     is to zero, the steeper the corresponding part of the
%                     log10 plot used, and the greater the compression (the
%                     default compression factor is 0.02).
%       'MovieArgs' - Settings for movie2avi. Parameters must be given as
%                     {param, value, ...} pairs within a cell array
%                     (default = {}).
%       'MovieName' - Name of the movie produced when 'RecordMovie' is set
%                     to true (default = 'date-time-kspaceFirstOrder2D').
%       'PlotFreq'  - The number of iterations which must pass before the
%                     simulation plot is updated (default = 10).
%       'PlotLayout'- Boolean controlling whether three plots are produced
%                     of the initial simulation layout (initial pressure,
%                     sound speed, density) (default = false).
%       'PlotPML'   - Boolean controlling whether the perfectly matched
%                     layer is shown in the simulation plots. If set to
%                     false, the PML is not displayed (default = true).
%       'PlotScale' - [min, max] values used to control the scaling for
%                     imagesc (visualisation). If set to 'auto', a
%                     symmetric plot scale is chosen automatically for each
%                     plot frame.
%       'PlotSim'   - Boolean controlling whether the simulation iterations
%                     are progressively plotted (default = true).
%       'PMLAlpha'  - Absorption within the perfectly matched layer in
%                     Nepers per metre (default = 4).
%       'PMLInside' - Boolean controlling whether the perfectly matched
%                     layer is inside or outside the grid. If set to false,
%                     the input grids are enlarged by PMLSize before
%                     running the simulation (default = true). 
%       'PMLSize'   - size of the perfectly matched layer in voxels. By
%                     default, the PML is added evenly to all sides of the
%                     grid, however, both PMLSize and PMLAlpha can be given
%                     as three element arrays to specify the x, y, and z
%                     properties, respectively. To remove the PML, set the
%                     appropriate PMLAlpha to zero rather than forcing the
%                     PML to be of zero size (default = 10).
%       'RecordMovie' - Boolean controlling whether the displayed image
%                     frames are captured and stored as a movie using
%                     movie2avi (default = false). 
%       'ReturnVelocity' - Boolean controlling whether the acoustic
%                     particle velocity at the positions defined by
%                     sensor.mask are also returned. If set to true, the
%                     output argument sensor_data is returned as a
%                     structure with the pressure and particle velocity
%                     appended as separate fields. In three dimensions,
%                     these fields are given by sensor_data.p,
%                     sensor_data.ux, sensor_data.uy, and sensor_data.uz.
%                     The field data is similarly returned as field_data.p,
%                     field_data.ux, field_data.uy, and field_data.uz.  
%       'Smooth'    - Boolean controlling whether source.p0,
%                     medium.sound_speed, and medium.density are smoothed
%                     using smooth before computation. 'Smooth' can either
%                     be given as a single Boolean value or as a 3 element
%                     array to control the smoothing of source.p0,
%                     medium.sound_speed, and medium.density,
%                     independently.  
%
% OUTPUTS:
% If 'ReturnVelocity' is set to false:
%       sensor_data - time varying pressure recorded at the sensor
%                     positions given by sensor.mask
%       field_data  - final pressure field
%
% If 'ReturnVelocity' is set to true:
%       sensor_data.p   - time varying pressure recorded at the sensor
%                         positions given by sensor.mask 
%       sensor_data.ux  - time varying particle velocity in the x-direction
%                         recorded at the sensor positions given by
%                         sensor.mask  
%       sensor_data.uy  - time varying particle velocity in the y-direction
%                         recorded at the sensor positions given by
%                         sensor.mask  
%       sensor_data.uz  - time varying particle velocity in the z-direction
%                         recorded at the sensor positions given by
%                         sensor.mask  
%       field_data.p    - final pressure field
%       field_data.ux   - final field of the particle velocity in the
%                         x-direction 
%       field_data.uy   - final field of the particle velocity in the
%                         y-direction 
%       field_data.uz   - final field of the particle velocity in the
%                         z-direction 
%
% ABOUT:
%       author      - Bradley Treeby and Ben Cox
%       date        - 7th April 2009
%       last update - 10th 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
%
% See also fftn, ifftn, imagesc, kspaceFirstOrder1D, kspaceFirstOrder2D,
% makeGrid, makeTime, smooth, unmaskSensorData 

% 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/>.

% suppress mlint warnings that arise from using sub-scripts
%#ok<*NASGU>
%#ok<*COLND>
%#ok<*NODEF>
%#ok<*INUSL>

% start the timer and store the start time
start_time = clock;
tic;

% =========================================================================
% DEFINE LITERALS
% =========================================================================

% minimum number of input variables
NUM_REQ_INPUT_VARIABLES = 4;

% optional input defaults
CARTESIAN_INTERP_DEF = 'nearest';
CREATE_LOG_DEF = false;
DATA_CAST_DEF = 'off';
DISPLAY_MASK_DEF = 'default';
LOG_SCALE_DEF = false;
LOG_SCALE_COMPRESSION_FACTOR_DEF = 0.02;
MOVIE_ARGS_DEF = {};
MOVIE_NAME_DEF = [getDateString '-kspaceFirstOrder3D'];
PLOT_FREQ_DEF = 10;
PLOT_LAYOUT_DEF = false;
PLOT_SCALE_DEF = [-1 1];
PLOT_SCALE_LOG_DEF = false;
PLOT_SIM_DEF = true;
PLOT_PML_DEF = true;
PML_ALPHA_DEF = 4;
PML_INSIDE_DEF = true;
PML_SIZE_DEF = 10;
RECORD_MOVIE_DEF = false;
RETURN_VELOCITY_DEF = false;
SMOOTH_P0_DEF = true;
SMOOTH_C0_DEF = false;
SMOOTH_RHO0_DEF = false;
USE_KSPACE_DEF = true;
USE_SG_DEF = true;

% set default movie compression
MOVIE_COMP_WIN = 'Cinepak';
MOVIE_COMP_LNX = 'None';
MOVIE_COMP_64B = 'None';

% set additional literals
MFILE = mfilename;
COLOR_MAP = getColorMap;
DT_WARNING_CFL = 0.45;  
ESTIMATE_SIM_TIME_STEPS = 10;
LOG_NAME = ['k-Wave-Log-' getDateString];
PLOT_SCALE_WARNING = 20;
LOG_SCALE_COMPRESSION_FACTOR = 0.01;

% =========================================================================
% CHECK INPUTS STRUCTURES AND OPTIONAL INPUTS
% =========================================================================

kspaceFirstOrder_inputChecking;

% forced exit for release B.0.1 inputs
if early_exit
    return
end

% =========================================================================
% UPDATE COMMAND LINE STATUS
% =========================================================================

disp(['  start time: ' datestr(start_time)]);
disp(['  dt: ' scaleSI(dt) 's, t_end: ' scaleSI(t_array(end)) 's, time steps: ' num2str(length(t_array))]);
[x_sc scale prefix] = scaleSI(min([kgrid.x_size, kgrid.y_size, kgrid.z_size])); %#ok<ASGLU>
disp(['  input grid size: ' num2str(kgrid.Nx) ' by ' num2str(kgrid.Ny) ' by ' num2str(kgrid.Nz) ' pixels (' num2str(kgrid.x_size*scale) ' by ' num2str(kgrid.y_size*scale) ' by ' num2str(kgrid.z_size*scale) prefix 'm)']); 
if (kgrid.kx_max == kgrid.kz_max) && (kgrid.kx_max == kgrid.ky_max)
    disp(['  maximum supported frequency: ' scaleSI( kgrid.k_max * min(c(:)) / (2*pi) ) 'Hz']);
else
    disp(['  maximum supported frequency: ' scaleSI( kgrid.kx_max * min(c(:)) / (2*pi) ) 'Hz by ' scaleSI( kgrid.ky_max * min(c(:)) / (2*pi) ) 'Hz by ' scaleSI( kgrid.kz_max * min(c(:)) / (2*pi) ) 'Hz']);
end

% =========================================================================
% SMOOTH AND ENLARGE INPUT GRIDS
% =========================================================================

% smooth p0 distribution if required and then restore maximum magnitude
if isfield(source, 'p0') && smooth_p0
    disp('  smoothing p0 distribution...');      
    source.p0 = smooth(kgrid, source.p0, true);
end

% expand grid if the PML is set to be outside the input grid
if ~PML_inside

    % expand the computational grid
    disp('  expanding computational grid...');
    kgrid = makeGrid(kgrid.Nx + 2*PML_x_size, kgrid.dx, kgrid.Ny + 2*PML_y_size, kgrid.dy, kgrid.Nz + 2*PML_z_size, kgrid.dz);
               
    % expand the grid matrices
    expand_size = [PML_x_size, PML_y_size, PML_z_size]; %#ok<NASGU>
    kspaceFirstOrder_expandGridMatrices;
    clear expand_size;
    
    % update command line status
    disp(['  computational grid size: ' num2str(kgrid.Nx) ' by ' num2str(kgrid.Ny) ' by ' num2str(kgrid.Nz) ' pixels']);
end

if ~plot_PML
    % create indexes to allow inputs to be placed into the larger
    % simulation grid
    x1 = (PML_x_size + 1);
    x2 = kgrid.Nx - PML_x_size;
    y1 = (PML_y_size + 1);
    y2 = kgrid.Ny - PML_y_size;    
    z1 = (PML_z_size + 1);
    z2 = kgrid.Nz - PML_z_size;
else
    % create indexes to place the source input exactly into the simulation
    % grid
    x1 = 1;
    x2 = kgrid.Nx;
    y1 = 1;
    y2 = kgrid.Ny;
    z1 = 1;
    z2 = kgrid.Nz;
end    

% select reference sound speed based on heterogeneity maps
c_max = max(c(:));

% smooth c distribution if required
if smooth_c && numDim(c) == 3
    disp('  smoothing sound speed distribution...');      
    c = smooth(kgrid, c);
end
    
% smooth rho0 distribution if required
if smooth_rho0 && numDim(rho0) == 3
    disp('  smoothing density distribution...');      
    rho0 = smooth(kgrid, rho0);
end

% =========================================================================
% PREPARE STAGGERED COMPUTATIONAL GRIDS AND OPERATORS
% =========================================================================

% create the staggered grids
x_sg = kgrid.x + kgrid.dx/2;
y_sg = kgrid.y + kgrid.dy/2;
z_sg = kgrid.z + kgrid.dz/2;

% interpolate the values of the density at the staggered grid locations
% where r1 = (x + dx/2, y, z), r2 = (x, y + dy/2, z), r3 = (x, y, z + dz/2)
% values outside of the interpolation range are replaced with their
% original values 
if numDim(rho0) == 3   && use_sg
    
    % rho0 is heterogeneous and staggered grids are used
    rho0_r1 = interpn(kgrid.z, kgrid.x, kgrid.y, rho0, kgrid.z, x_sg, kgrid.y, '*linear');
    rho0_r2 = interpn(kgrid.z, kgrid.x, kgrid.y, rho0, kgrid.z, kgrid.x, y_sg, '*linear');
    rho0_r3 = interpn(kgrid.z, kgrid.x, kgrid.y, rho0, z_sg, kgrid.x, kgrid.y, '*linear');
    
    % set values outside of the interpolation range to original values
    rho0_r1(isnan(rho0_r1)) = rho0(isnan(rho0_r1));
    rho0_r2(isnan(rho0_r2)) = rho0(isnan(rho0_r2));    
    rho0_r3(isnan(rho0_r3)) = rho0(isnan(rho0_r3));
    
else
    % rho0 is homogeneous or staggered grids are not used
    rho0_r1 = rho0;
    rho0_r2 = rho0;
    rho0_r3 = rho0;
end

% define the location of the perfectly matched layer within the grid
x0_min = kgrid.x(1) + PML_x_size*kgrid.dx;
x0_max = kgrid.x(end) - PML_x_size*kgrid.dx;
y0_min = kgrid.y(1) + PML_y_size*kgrid.dy;
y0_max = kgrid.y(end) - PML_y_size*kgrid.dy;
z0_min = kgrid.z(1) + PML_z_size*kgrid.dz;
z0_max = kgrid.z(end) - PML_z_size*kgrid.dz;

% set the PML attenuation over the pressure (regular) grid
ax = PML_x_alpha*(c_max/kgrid.dx)*((kgrid.x - x0_max)./(kgrid.x(end) - x0_max)).^4.*(kgrid.x >= x0_max)...
    + PML_x_alpha*(c_max/kgrid.dx)*((kgrid.x - x0_min)./(kgrid.x(1) - x0_min)).^4.*(kgrid.x <= x0_min);

ay = PML_y_alpha*(c_max/kgrid.dy)*((kgrid.y - y0_max)./(kgrid.y(end) - y0_max)).^4.*(kgrid.y >= y0_max)...
    + PML_y_alpha*(c_max/kgrid.dy)*((kgrid.y - y0_min)./(kgrid.y(1) - y0_min)).^4.*(kgrid.y <= y0_min);

az = PML_z_alpha*(c_max/kgrid.dz)*((kgrid.z - z0_max)./(kgrid.z(end) - z0_max)).^4.*(kgrid.z >= z0_max)...
    + PML_z_alpha*(c_max/kgrid.dz)*((kgrid.z - z0_min)./(kgrid.z(1) - z0_min)).^4.*(kgrid.z <= z0_min);

% set the PML attenuation over the velocity (staggered) grid
ax_sg = PML_x_alpha*(c_max/kgrid.dx)*((x_sg - x0_max)./(kgrid.x(end) - x0_max)).^4.*(x_sg >= x0_max)...
    + PML_x_alpha*(c_max/kgrid.dx)*((x_sg - x0_min)./(kgrid.x(1) - x0_min)).^4.*(x_sg <= x0_min);

ay_sg = PML_y_alpha*(c_max/kgrid.dy)*((y_sg - y0_max)./(kgrid.y(end) - y0_max)).^4.*(y_sg >= y0_max)...
    + PML_y_alpha*(c_max/kgrid.dy)*((y_sg - y0_min)./(kgrid.y(1) - y0_min)).^4.*(y_sg <= y0_min);

az_sg = PML_z_alpha*(c_max/kgrid.dz)*((z_sg - z0_max)./(kgrid.z(end) - z0_max)).^4.*(z_sg >= z0_max)...
    + PML_z_alpha*(c_max/kgrid.dz)*((z_sg - z0_min)./(kgrid.z(1) - z0_min)).^4.*(z_sg <= z0_min);

% precompute absorbing boundary condition operators
abc_x = exp(-ax_sg*dt/2);
abc_y = exp(-ay_sg*dt/2);
abc_z = exp(-az_sg*dt/2);
abc_x_alt = exp(-ax*dt);
abc_y_alt = exp(-ay*dt);
abc_z_alt = exp(-az*dt);

% define the modified first order k-space derivative operators and
% multiply by the staggered grid shift operators (options use_kspace and
% use_sg exist for debugging)
if use_kspace
    % k-space
    if use_sg
        ddx_k_shift_pos = 1i*kgrid.kx.*sinc(c_max*dt*kgrid.k/2) .* exp(1i*kgrid.kx*kgrid.dx/2);
        ddx_k_shift_neg = 1i*kgrid.kx.*sinc(c_max*dt*kgrid.k/2) .* exp(-1i*kgrid.kx*kgrid.dx/2);
        ddy_k_shift_pos = 1i*kgrid.ky.*sinc(c_max*dt*kgrid.k/2) .* exp(1i*kgrid.ky*kgrid.dy/2);
        ddy_k_shift_neg = 1i*kgrid.ky.*sinc(c_max*dt*kgrid.k/2) .* exp(-1i*kgrid.ky*kgrid.dy/2);
        ddz_k_shift_pos = 1i*kgrid.kz.*sinc(c_max*dt*kgrid.k/2) .* exp(1i*kgrid.kz*kgrid.dz/2);
        ddz_k_shift_neg = 1i*kgrid.kz.*sinc(c_max*dt*kgrid.k/2) .* exp(-1i*kgrid.kz*kgrid.dz/2);
    else
        ddx_k_shift_pos = 1i*kgrid.kx.*sinc(c_max*dt*kgrid.k/2);
        ddx_k_shift_neg = 1i*kgrid.kx.*sinc(c_max*dt*kgrid.k/2);
        ddy_k_shift_pos = 1i*kgrid.ky.*sinc(c_max*dt*kgrid.k/2);
        ddy_k_shift_neg = 1i*kgrid.ky.*sinc(c_max*dt*kgrid.k/2);
        ddz_k_shift_pos = 1i*kgrid.kz.*sinc(c_max*dt*kgrid.k/2);
        ddz_k_shift_neg = 1i*kgrid.kz.*sinc(c_max*dt*kgrid.k/2);        
    end
else
    % pseudospectral
    if use_sg
        ddx_k_shift_pos = 1i*kgrid.kx .* exp(1i*kgrid.kx*kgrid.dx/2);
        ddx_k_shift_neg = 1i*kgrid.kx .* exp(-1i*kgrid.kx*kgrid.dx/2);
        ddy_k_shift_pos = 1i*kgrid.ky .* exp(1i*kgrid.ky*kgrid.dy/2);
        ddy_k_shift_neg = 1i*kgrid.ky .* exp(-1i*kgrid.ky*kgrid.dy/2);
        ddz_k_shift_pos = 1i*kgrid.kz .* exp(1i*kgrid.kz*kgrid.dz/2);
        ddz_k_shift_neg = 1i*kgrid.kz .* exp(-1i*kgrid.kz*kgrid.dz/2);
    else
        ddx_k_shift_pos = 1i*kgrid.kx;
        ddx_k_shift_neg = 1i*kgrid.kx;
        ddy_k_shift_pos = 1i*kgrid.ky;
        ddy_k_shift_neg = 1i*kgrid.ky;
        ddz_k_shift_pos = 1i*kgrid.kz;
        ddz_k_shift_neg = 1i*kgrid.kz;        
    end
end

% pre-shift variables used as frequency domain multipliers
ddx_k_shift_pos = ifftshift(ddx_k_shift_pos);
ddx_k_shift_neg = ifftshift(ddx_k_shift_neg);
ddy_k_shift_pos = ifftshift(ddy_k_shift_pos);
ddy_k_shift_neg = ifftshift(ddy_k_shift_neg);
ddz_k_shift_pos = ifftshift(ddz_k_shift_pos);
ddz_k_shift_neg = ifftshift(ddz_k_shift_neg);

% cleanup unused variables
clear ax* ay* az* x0_min x0_max x_sg PML_x_alpha y0_min y0_max y_sg PML_y_alpha z0_min z0_max z_sg PML_z_alpha;

% =========================================================================
% PREPARE DATA MASKS AND STORAGE VARIABLES
% =========================================================================

kspaceFirstOrder_createAbsorptionAndStorageVariables;

% =========================================================================
% SET INITIAL CONDITIONS
% =========================================================================

% add the initial pressure to rho as a mass source
if isfield(source, 'p0')
    rhox = source.p0./(3*c.^2);
    rhoy = source.p0./(3*c.^2);
    rhoz = source.p0./(3*c.^2);
else
    rhox = zeros(kgrid.Nz, kgrid.Nx, kgrid.Ny);
    rhoy = zeros(kgrid.Nz, kgrid.Nx, kgrid.Ny);
    rhoz = zeros(kgrid.Nz, kgrid.Nx, kgrid.Ny);   
end

% add p_source(t = t1) to rho as a mass source and create source variables
if p_source 
        
    % scale the pressure source by 2*dt*c/dx divided by 3*c^2 (because the
    % source is added to the density and divided into three components) 
    if numel(c) == 1
        source.p = source.p .* (2*dt./(3*c*kgrid.dx));
    else
        % compute the scale parameter seperately for each source position
        % based on the sound speed at that position
        for p_index = 1:length(source.p(:,1))
            source.p(p_index, :) = source.p(p_index, :) .* (2.*dt./(3*c(p_index).*kgrid.dx));
        end
    end
    
    % extract the initial source pressure and add to rho as a mass source
    rhox(ps_index) = rhox(ps_index) + source.p(p_element_index, 1);
    rhoy(ps_index) = rhoy(ps_index) + source.p(p_element_index, 1);
    rhoz(ps_index) = rhoz(ps_index) + source.p(p_element_index, 1);    
    
end

% calculate initial pressure using an adiabatic equation of state
p = c.^2.*(rhox + rhoy + rhoz);

% compute u(t = t1 + dt/2) based on the assumption u(dt/2) = -u(-dt/2)
% which forces u(t = t1) = 0
ux_r1 = dt./rho0_r1 .* real(ifftn(ddx_k_shift_pos .* fftn(p)));
uy_r2 = dt./rho0_r2 .* real(ifftn(ddy_k_shift_pos .* fftn(p)));
uz_r3 = dt./rho0_r3 .* real(ifftn(ddz_k_shift_pos .* fftn(p)));
    
% add u_source(t = t1) to u
if u_source
    
    % scale the velocity source by 2*dt*c/dx
    if numel(c) == 1
        source.ux = source.ux .* (2*c*dt./kgrid.dx);
        source.uy = source.uy .* (2*c*dt./kgrid.dy);
        source.uz = source.uz .* (2*c*dt./kgrid.dz);
    else
        % compute the scale parameter seperately for each source position
        % based on the sound speed at that position
        for u_index = 1:length(source.ux(:,1))        
            source.ux(u_index, :) = source.ux(u_index, :) .* (2.*c(u_index).*dt./kgrid.dx);
            source.uy(u_index, :) = source.uy(u_index, :) .* (2.*c(u_index).*dt./kgrid.dy);
            source.uz(u_index, :) = source.uz(u_index, :) .* (2.*c(u_index).*dt./kgrid.dz);
        end
    end
    
    % extract the initial source velocity and add to u
    ux_r1(us_index) = ux_r1(us_index) + source.ux(u_element_index, 1);
    uy_r2(us_index) = uy_r2(us_index) + source.uy(u_element_index, 1);
    uz_r3(us_index) = uz_r3(us_index) + source.uz(u_element_index, 1);
    
end

% setup the time index variable
if ~time_rev
    index_start = 2;
    index_step = 1;
    index_end = length(t_array); 
else
    % reverse the order of the input data
    sensor.time_reversal_boundary_data = fliplr(sensor.time_reversal_boundary_data);
    index_start = 1;
    index_step = 1;
    
    % stop one time point before the end so the last points are not
    % propagated
    index_end = length(t_array) - 1;
end

% store p(t = t1) in forward mode
if use_sensor && ~time_rev
    
    % precomputate the triangulation points
    if ~binary_sensor_mask
        % try to use TriScatteredInterp (only supported post R2009a)
        try
            disp('  calculating Delaunay triangulation (TriScatteredInterp)...');
            F_interp = TriScatteredInterp(reshape(kgrid.z, [], 1), reshape(kgrid.x, [], 1), reshape(kgrid.y, [], 1), reshape(p, [], 1));
        catch ME
            % if TriScatteredInterp doesn't exist, return an error
            if strcmp(ME.identifier, 'MATLAB:UndefinedFunction')
                error('The inbuilt function TriScatteredInterp is not supported in this version of MATLAB. Try running the simulation with ''CartInterp'' set to ''nearest''');
            else
                rethrow(ME);
            end
        end
    end    
   
    % store p(t = t1)
    switch extract_data_case
        case 1
            % return velocity = false
            % binary sensor mask = false

            F_interp.V = reshape(p, [], 1);
            sensor_data(:, 1) = F_interp(sensor_z, sensor_y, sensor_x);
            
        case 2
            % return velocity = false
            % binary sensor mask = true  
            
            sensor_data(:, 1) = p(sensor_mask_ind);
            
        case 3
            % return velocity = true
            % binary sensor mask = false            

            F_interp.V = reshape(p, [], 1);
            sensor_data.p(:, 1) = F_interp(sensor_z, sensor_y, sensor_x);
            F_interp.V = reshape(ux_r1, [], 1);
            sensor_data.ux(:, 1) = F_interp(sensor_z, sensor_y, sensor_x);
            F_interp.V = reshape(uy_r2, [], 1);
            sensor_data.uy(:, 1) = F_interp(sensor_z, sensor_y, sensor_x);            
            F_interp.V = reshape(uz_r3, [], 1);
            sensor_data.uz(:, 1) = F_interp(sensor_z, sensor_y, sensor_x);            

        case 4
            % return velocity = true
            % binary sensor mask = true
            
            sensor_data.p(:, 1) = p(sensor_mask_ind);
            sensor_data.ux(:, 1) = ux_r1(sensor_mask_ind);
            sensor_data.uy(:, 1) = uy_r2(sensor_mask_ind);            
            sensor_data.uz(:, 1) = uz_r3(sensor_mask_ind);                        
            
    end        
end


% =========================================================================
% PREPARE VISUALISATIONS
% =========================================================================

% pre-compute suitable axes scaling factor
if plot_layout || plot_sim
    [x_sc, scale, prefix] = scaleSI(max([kgrid.x(1,:,1), squeeze(kgrid.y(1,1,:))', kgrid.z(:,1,1)']));  %#ok<ASGLU>
end

% plot the simulation layout
if plot_layout
    
    % initial pressure
    if isfield(source, 'p0')
        
        % create plot variable
        p_plot = p(z1:z2, x1:x2, y1:y2);  
        
        % update plot scale if set to automatic or log
        if plot_scale_auto || plot_scale_log
            kspaceFirstOrder_adjustPlotScale;
        end        

        % plot the initial pressure distribution
        figure;
        planeplot(scale*kgrid.x_vec(x1:x2), scale*kgrid.y_vec(y1:y2), scale*kgrid.z_vec(z1:z2), p_plot, 'Initial Pressure: ', plot_scale, prefix, COLOR_MAP);

    end
    
    % plot c if heterogeneous
    if numDim(c) == 3
        c_plot_scale = [0.5*min(c(:)), max(c(:))/0.5];
        figure;
        planeplot(scale*kgrid.x_vec(x1:x2), scale*kgrid.y_vec(y1:y2), scale*kgrid.z_vec(z1:z2), c(z1:z2, x1:x2, y1:y2), 'c: ', c_plot_scale, prefix, COLOR_MAP);
    end

    % plot rho0 if heterogeneous
    if numDim(rho0) == 3    
        rho0_plot_scale = [0.5*min(rho0(:)), max(rho0(:))/0.5];
        figure;
        planeplot(scale*kgrid.x_vec(x1:x2), scale*kgrid.y_vec(y1:y2), scale*kgrid.z_vec(z1:z2),  rho0(z1:z2, x1:x2, y1:y2), 'rho0: ', rho0_plot_scale, prefix, COLOR_MAP);
    end

end

% initialise the figures used for animation
if plot_sim
    img = figure;
    if ~time_rev
        pbar = waitbar(0, 'Computing Pressure Field');
    else
        pbar = waitbar(0, 'Computing Time Reversed Field');
    end
end 

% set movie index variable
if record_movie
    
    % force getframe compatability with dual monitors
    movegui(img);
    
    % define frame index
    frame_index = 1;
    
    % save the first frame if in forward mode
    if ~time_rev
        
        % create plot variable
        p_plot = p(z1:z2, x1:x2, y1:y2);       

        % update plot scale if set to automatic or log
        if plot_scale_auto || plot_scale_log
            kspaceFirstOrder_adjustPlotScale;
        end         
        
        % add display mask onto plot
        if strcmp(display_mask, 'default')
            p_plot(sensor.mask(z1:z2, x1:x2, y1:y2) ~= 0) = plot_scale(2);
        elseif ~strcmp(display_mask, 'off')
            p_plot(display_mask(z1:z2, x1:x2, y1:y2) ~= 0) = plot_scale(2);
        end

        % update plot
        planeplot(scale*kgrid.x_vec(x1:x2), scale*kgrid.y_vec(y1:y2), scale*kgrid.z_vec(z1:z2), p_plot, '', plot_scale, prefix, COLOR_MAP);

        % force plot update
        drawnow;        
        
        % set background color to white
        set(gcf, 'Color', [1 1 1]);

        % save the movie frame
        movie_frames(frame_index) = getframe(gcf);

        % update frame index
        frame_index  = frame_index  + 1;
    end
end

% =========================================================================
% DATA CASTING
% =========================================================================

kspaceFirstOrder_dataCast;

% =========================================================================
% LOOP THROUGH TIME STEPS
% =========================================================================

% update command line status
disp(['  precomputation completed in ' scaleTime(toc)]);
disp('  starting time loop...');

% restart timing variable
tic;

% precompute fft of p
p_k = fftn(p);

% start time loop
for t_index = index_start:index_step:index_end

    % enforce time reversal bounday condition
    if time_rev

        % load pressure value and enforce as a Dirichlet boundary condition
        p(sensor_mask_ind) = sensor.time_reversal_boundary_data(:, t_index);

        % update p_k
        p_k = fftn(p);

        % compute rhox and rhoz using an adiabatic equation of state
        rhox_mod = p./(3*c.^2);
        rhoy_mod = p./(3*c.^2);
        rhoz_mod = p./(3*c.^2);
        rhox(sensor_mask_ind) = rhox_mod(sensor_mask_ind);
        rhoy(sensor_mask_ind) = rhoy_mod(sensor_mask_ind);
        rhoz(sensor_mask_ind) = rhoz_mod(sensor_mask_ind);
           
    end    

    % calculate ux, uy and uz at t + dt/2 using dp/dx, dp/dy and dp/dz at t
    ux_r1 = abc_x .* (  abc_x.*ux_r1 - dt./rho0_r1 .* real(ifftn(ddx_k_shift_pos .* p_k))  );
    uy_r2 = abc_y .* (  abc_y.*uy_r2 - dt./rho0_r2 .* real(ifftn(ddy_k_shift_pos .* p_k))  );
    uz_r3 = abc_z .* (  abc_z.*uz_r3 - dt./rho0_r3 .* real(ifftn(ddz_k_shift_pos .* p_k))  );

    % add in the velocity source term
    if u_source
        ux_r1(us_index) = ux_r1(us_index) + source.ux(u_element_index, t_index);
        uy_r2(us_index) = uy_r2(us_index) + source.uy(u_element_index, t_index);
        uz_r3(us_index) = uz_r3(us_index) + source.uz(u_element_index, t_index);
    end    
    
    % calculate dux/dx, duydy and duz/dz at t + dt/2
    duxdx_k = ddx_k_shift_neg .* fftn(ux_r1);
    duydy_k = ddy_k_shift_neg .* fftn(uy_r2);
    duzdz_k = ddz_k_shift_neg .* fftn(uz_r3);   

    % calculate rhox, rhoy and rhoz at t + dt
    rhox = abc_x_alt .* (  rhox - dt.*rho0 .* real(ifftn(duxdx_k))  );
    rhoy = abc_y_alt .* (  rhoy - dt.*rho0 .* real(ifftn(duydy_k))  );        
    rhoz = abc_z_alt .* (  rhoz - dt.*rho0 .* real(ifftn(duzdz_k))  );      
    
    % add in the pressure source term as a mass source   
    if p_source
        rhox(ps_index) = rhox(ps_index) + source.p(p_element_index, t_index);
        rhoy(ps_index) = rhoy(ps_index) + source.p(p_element_index, t_index);  
        rhoz(ps_index) = rhoz(ps_index) + source.p(p_element_index, t_index);          
    end
    
    switch equation_of_state
        case 'lossless';
            % calculate p using an adiabatic equation of state
            p = c.^2.*(rhox + rhoy + rhoz);
        case 'absorbing';
            % calculate p using an absorbing equation of state 
            p = c.^2.*(  (rhox + rhoy + rhoz) + real(ifftn( absorb_param.*(duxdx_k + duydy_k + duzdz_k) + dispers_param.*fftn(rhox + rhoy + rhoz) ))  );            
        case 'stokes'
            % calculate p using stokes' equation of state
            p = c.^2.*(  (rhox + rhoy + rhoz) + absorb_param .* real(ifftn(duxdx_k + duydy_k + duzdz_k))  );
    end        

    % precompute fft of p here so it can be modified for visualisation
    p_k = fftn(p);        

    % extract required data
    switch extract_data_case
        case 1
            % return velocity = false
            % binary sensor mask = false
            
            F_interp.V = reshape(p, [], 1);
            sensor_data(:, t_index) = F_interp(sensor_z, sensor_y, sensor_x);
            
        case 2
            % return velocity = false
            % binary sensor mask = true  
            
            sensor_data(:, t_index) = p(sensor_mask_ind);
                        
        case 3
            % return velocity = true
            % binary sensor mask = false            

            F_interp.V = reshape(p, [], 1);
            sensor_data.p(:, t_index) = F_interp(sensor_z, sensor_y, sensor_x);
            F_interp.V = reshape(ux_r1, [], 1);
            sensor_data.ux(:, t_index) = F_interp(sensor_z, sensor_y, sensor_x);
            F_interp.V = reshape(uy_r2, [], 1);
            sensor_data.uy(:, t_index) = F_interp(sensor_z, sensor_y, sensor_x);            
            F_interp.V = reshape(uz_r3, [], 1);
            sensor_data.uz(:, t_index) = F_interp(sensor_z, sensor_y, sensor_x);            

        case 4
            % return velocity = true
            % binary sensor mask = true
            
            sensor_data.p(:, t_index) = p(sensor_mask_ind);
            sensor_data.ux(:, t_index) = ux_r1(sensor_mask_ind);
            sensor_data.uy(:, t_index) = uy_r2(sensor_mask_ind);            
            sensor_data.uz(:, t_index) = uz_r3(sensor_mask_ind);                    
    end        

    % estimate the time to run the simulation
    if t_index == ESTIMATE_SIM_TIME_STEPS
        disp(['  estimated simulation time ' scaleTime(toc*index_end/t_index) '...']);
    end
    
    % plot data if required
    if plot_sim && rem(t_index, plot_freq) == 0           

        % update progress bar
        waitbar(t_index/length(t_array), pbar);
        drawnow;

        % ensure p is cast as a CPU variable
        p_plot = double(p(z1:z2, x1:x2, y1:y2));       

        % update plot scale if set to automatic or log
        if plot_scale_auto || plot_scale_log
            kspaceFirstOrder_adjustPlotScale;
        end         
        
        % add display mask onto plot
        if strcmp(display_mask, 'default')
            p_plot(sensor.mask(z1:z2, x1:x2, y1:y2) ~= 0) = plot_scale(2);
        elseif ~strcmp(display_mask, 'off')
            p_plot(display_mask(z1:z2, x1:x2, y1:y2) ~= 0) = plot_scale(2);
        end     
        
        % update plot
        planeplot(scale*kgrid.x_vec(x1:x2), scale*kgrid.y_vec(y1:y2), scale*kgrid.z_vec(z1:z2), p_plot, '',  plot_scale, prefix, COLOR_MAP);

        % save movie frames if required
        if record_movie
            
            % set background color to white
            set(gcf, 'Color', [1 1 1]);

            % save the movie frame
            movie_frames(frame_index) = getframe(gcf);

            % update frame index
            frame_index  = frame_index  + 1;
            
        end
    end
end

% assign the final time reversal values
if time_rev
    p(sensor_mask_ind) = sensor.time_reversal_boundary_data(:, index_end + 1);
end

% update command line status
disp(['  simulation completed in ' scaleTime(toc)]);

% =========================================================================
% CLEAN UP
% =========================================================================

% store the final movie frame and then save the movie
if record_movie
    
    % update command line status
    disp('  saving movie file...');
    
    % ensure p is cast as a CPU variable
    p_plot = double(p(z1:z2, x1:x2, y1:y2));       

    % update plot scale if set to automatic or log
    if plot_scale_auto || plot_scale_log
        kspaceFirstOrder_adjustPlotScale;
    end    
    
    % add display mask onto plot
    if strcmp(display_mask, 'default')
        p_plot(sensor.mask(z1:z2, x1:x2, y1:y2) ~= 0) = plot_scale(2);
    elseif ~strcmp(display_mask, 'off')
        p_plot(display_mask(z1:z2, x1:x2, y1:y2) ~= 0) = plot_scale(2);
    end

    % update plot
    planeplot(scale*kgrid.x_vec(x1:x2), scale*kgrid.y_vec(y1:y2), scale*kgrid.z_vec(z1:z2), p_plot, '',  plot_scale, prefix, COLOR_MAP);

    % force plot update
    drawnow;
        
    % set background color to white
    set(gcf, 'Color', [1 1 1]);

    % save the movie frame
    movie_frames(frame_index) = getframe(gcf);
    
    % save movie file
    kspaceFirstOrder_saveMovieFile;

end

% clean up used figures
if plot_sim
    close(img);
    close(pbar);
end

% reset the indexing variables to allow original grid size to be maintained
if (~plot_PML && PML_inside)
    x1 = x1 - PML_x_size;
    x2 = x2 + PML_x_size;
    y1 = y1 - PML_y_size;
    y2 = y2 + PML_y_size;    
    z1 = z1 - PML_z_size;
    z2 = z2 + PML_z_size;
elseif (plot_PML && ~PML_inside)
    x1 = x1 + PML_x_size;
    x2 = x2 - PML_x_size;
    y1 = y1 + PML_y_size;
    y2 = y2 - PML_y_size;
    z1 = z1 + PML_z_size;
    z2 = z2 - PML_z_size;    
end

% save the final pressure field if in time reversal mode
if time_rev
    sensor_data = p(z1:z2, x1:x2, y1:y2);
end

% save the final pressure and velocity fields and cast variables back to
% double
if return_velocity
    field_data.p = double(p(z1:z2, x1:x2, y1:y2));
    field_data.ux = double(ux_r1(z1:z2, x1:x2, y1:y2));
    field_data.uy = double(uy_r2(z1:z2, x1:x2, y1:y2));
    field_data.uz = double(uz_r3(z1:z2, x1:x2, y1:y2));
    if use_sensor
        sensor_data.p = double(sensor_data.p);
        sensor_data.ux = double(sensor_data.ux);
        sensor_data.uy = double(sensor_data.uy);
        sensor_data.uz = double(sensor_data.uz);      
    end
else
    field_data = double(p(z1:z2, x1:x2, y1:y2));
    if use_sensor
        sensor_data = double(sensor_data);
    end
end

% reorder the sensor points if binary sensor mask was used for Cartesian
% sensor mask nearest neighbour interpolation
if use_sensor && reorder_data
    
    % update command line status
    disp('  reordering Cartesian measurement data...');
    
    if return_velocity
        
        % append the reordering data
        new_col_pos = length(sensor_data.p(1,:)) + 1;
        sensor_data.p(:, new_col_pos) = reorder_index;
        sensor_data.ux(:, new_col_pos) = reorder_index;
        sensor_data.uy(:, new_col_pos) = reorder_index;
        sensor_data.uz(:, new_col_pos) = reorder_index;        

        % reorder p0 based on the order_index
        sensor_data.p = sortrows(sensor_data.p, new_col_pos);
        sensor_data.ux = sortrows(sensor_data.ux, new_col_pos);
        sensor_data.uy = sortrows(sensor_data.uy, new_col_pos);
        sensor_data.uz = sortrows(sensor_data.uz, new_col_pos);

        % remove the reordering data
        sensor_data.p = sensor_data.p(:, 1:new_col_pos - 1);
        sensor_data.ux = sensor_data.ux(:, 1:new_col_pos - 1);
        sensor_data.uy = sensor_data.uy(:, 1:new_col_pos - 1);
        sensor_data.uz = sensor_data.uz(:, 1:new_col_pos - 1);   
        
    else
        
        % append the reordering data
        new_col_pos = length(sensor_data(1,:)) + 1;
        sensor_data(:, new_col_pos) = reorder_index;

        % reorder p0 based on the order_index
        sensor_data = sortrows(sensor_data, new_col_pos);

        % remove the reordering data
        sensor_data = sensor_data(:, 1:new_col_pos - 1);
        
    end
end

% filter the recorded time domain pressure signals if transducer parameters
% are given
kspaceFirstOrder_filterSensorData;

% return empty sensor data if not used
if ~use_sensor
    sensor_data = [];
end

% update command line status
disp(['  total computation time ' scaleTime(etime(clock, start_time))]);

% switch off log
if create_log
    diary off;
end

function planeplot(x_vec, y_vec, z_vec, data, data_title, plot_scale, prefix, color_map)
% Subfunction to produce a plot of a three-dimensional matrix through the
% three central planes
   
subplot(2, 2, 1), imagesc(x_vec, z_vec, squeeze(data(:, :, round(end/2))), plot_scale);
title([data_title 'x-z plane']);
axis image;

subplot(2, 2, 2), imagesc(y_vec, z_vec, squeeze(data(:, round(end/2), :)), plot_scale);
title('y-z plane');
axis image;
xlabel(['(All axes in ' prefix 'm)']);

subplot(2, 2, 3), imagesc(x_vec, y_vec, squeeze(data(round(end/2), :, :)).', plot_scale);
title('x-y plane');
axis image;
colormap(color_map); 
drawnow;