function [sensor_data p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor, varargin)
%KSPACEFIRSTORDER2D     2D time-domain simulation of wave propagation.
%
% DESCRIPTION:
%       kspaceFirstOrder2D simulates the time-domain propagation of linear
%       compressional waves through a two-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 a heterogeneous sound speed
%       and density and power law absorption. 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.
%  
%       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 the same size as the
%       grid defined by kgrid) 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 2 by N matrix corresponding
%       to the x and z positions, respectively. 
%
%       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.  
%
%       The code may also be used for time-reversal image reconstruction by
%       assigning the boundary data to the input field
%       sensor.time_reversal_boundary_data. This data is then enforced as a
%       time varying Dirichlet boundary condition over the 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 then 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 then the boundary data must be given as an array
%       ordered using MATLAB's standard column-wise linear matrix indexing.
% 
%       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 = kspaceFirstOrder2D(kgrid, medium, source, sensor)
%       sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, ...) 
%       [sensor_data p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor)
%       [sensor_data p_final] = kspaceFirstOrder2D(kgrid, medium, source, sensor, ...) 
%
% INPUTS:
% The fields that are required for a photoacoustic forward simulation are marked with a *. 
%
%       kgrid*              - k-Wave 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 [dB/(MHz^y cm)]
%       medium.alpha_coeff  - power law absorption coefficient 
%                             [dB/(MHz^y cm)] 
%
%       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 
%
%       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.time_reversal_adapt_thresh - adaptive boundary condition
%                             threshold 
%       sensor.directivity_angle - matrix of directivity angles (direction
%                             of maximum response) for each sensor element
%                             defined in sensor.mask. The angles are in
%                             radians where 0 = max sensitivity in z
%                             direction (up/down) and pi/2 or -pi/2 = max
%                             sensitivity in x direction (left/right) 
%       sensor.directivity_size - equivalent element size (the larger the
%                             element size the more directional the
%                             response)  
%
% 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. 
%
% 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 = 'linear'). 
%       'CreateLog' - Boolean controlling whether the command line output
%                     is saved using the diary function with a date and
%                     time stamped filename. 
%       '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 fft2 and ifft2 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. 
%       'MovieName' - Name of the movie produced when 'RecordMovie' is set
%                     to true (default = 'date-time-kspaceFirstOrder2D').
%       'MovieArgs' - Settings for movie2avi. Parameters must be given as
%                     {param, value, ...} pairs within a cell array
%                     (default = {}).
%       'PlotFreq'  - The number of iterations which must pass before the
%                     simulation plot is updated (default = 10).
%       'PlotLayout'- Boolean controlling whether a four panel plot of the
%                     initial simulation layout is produced (initial
%                     pressure, sensor mask, sound speed, density)
%                     (default = false).
%       'PlotScale' - [min, max] values used to control the scaling for
%                     imagesc (visualisation) and im2frame (movie capture)
%                     (default = [-1 1]).
%       'PlotSim'   - Boolean controlling whether the simulation iterations
%                     are progressively plotted (default = true).
%       'PMLAlpha'  - attenuation in Nepers per m of the absorption within
%                     the perfectly matched layer (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 pixels. By
%                     default, the PML is added evenly to all sides of the
%                     grid, however, both PMLSize and PMLAlpha can be given
%                     as two element arrays to specify the x 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 = 20).
%       'RecordMovie' - Boolean controlling whether the displayed image
%                     frames are captured using im2frame and stored as a
%                     movie using movie2avi (default = false).
%       'Smooth'    - Boolean controlling whether the p0, c, and rho
%                     distributions are smoothed. Smooth can also be given
%                     as a 3 element array to control the smoothing of p0,
%                     c, and rho, respectively (default = true).
%
% OUTPUTS:
%       sensor_data - array of pressure time-series recorded at the sensor
%                     positions given by sensor_mask
%       p_final     - final pressure field
%
% ABOUT:
%       author      - Bradley Treeby and Ben Cox
%       date        - 25th February 2009
%       last update - 23rd 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
%
% See also fft2, ifft2, im2frame, imagesc, kspaceFirstOrder1D,
% kspaceFirstOrder3D, makeGrid, makeTime, movie2avi, 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/>.

% start the timer
tic;

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

% minimum number of input variables
NUM_REQ_INPUT_VARIABLES = 4;

% optional input defaults
CARTESIAN_INTERP_DEF = 'linear';
CREATE_LOG_DEF = false;
DATA_CAST_DEF = 'off';
DISPLAY_MASK_DEF = 'sensor';
LINEAR_ATTEN_DEF = 1.01;
LINEAR_ATTEN_TOLERANCE = 0.1;
MOVIE_ARGS_DEF = {};
MOVIE_NAME_DEF = [makeDateString() '-kspaceFirstOrder2D'];
PLOT_FREQ_DEF = 10;
PLOT_LAYOUT_DEF = false;
PLOT_SCALE_DEF = [-1 1];
PLOT_SIM_DEF = true;
PML_ALPHA_DEF = 4;
PML_INSIDE_DEF = true;
PML_SIZE_DEF = 20;
RECORD_MOVIE_DEF = false;
SMOOTH_DEF = true;

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

% set additional literals
COLOR_MAP = getColorMap;
DT_WARNING_CFL = 0.5;  
LOG_NAME = ['k-Wave-Log-' makeDateString()];
PLOT_SCALE_WARNING = 10;
     
% =========================================================================
% EXTRACT AND CHECK INPUTS
% =========================================================================

% check for release B.0.1 inputs where sensor_data = kspaceFirstOrder2D(p0,
% kgrid, c, rho, t_array, sensor_mask, ...) instead of sensor_data =
% kspaceFirstOrder2D(kgrid, medium, source, sensor, ...)
if ~isstruct(kgrid)
    
    % display warning message
    disp('WARNING: input usage deprecated, please see documentation.');
    disp('In future releases this usage will no longer be functional.');
    
    % extract old inputs
    p0_old = kgrid;
    kgrid_old = medium;
    c_old = source;
    rho_old = sensor;
    t_array_old = varargin{1};
    sensor_mask_old = varargin{2};
    clear kgrid medium source sensor;
    
    % reassign inputs
    kgrid = kgrid_old;
    kgrid.t_array = t_array_old;
    medium.sound_speed = c_old;
    medium.density = rho_old;
    if all(size(p0_old) == [kgrid_old.Nz, kgrid_old.Nx])
        source.p0 = p0_old;
    else
        source = [];
        sensor.time_reversal_boundary_data = p0_old;
    end
    sensor.mask = sensor_mask_old;
    
    % cleanup old inputs
    clear *_old;
    
    % repeat the function call with the correct inputs
    if nargin > 6
        sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor, varargin{3:end});
    else
        sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor);
    end
    return
end 

% check medium fields
checkFieldNames(medium, {'sound_speed', 'density', 'alpha_coeff', 'alpha_power'});
enforceFields(medium, {'sound_speed', 'density'});

% check medium absorption inputs
if isfield(medium, 'alpha_coeff') || isfield(medium, 'alpha_power')
    
    % if one absorption parameter is given, enforce the other
    enforceFields(medium, {'alpha_coeff', 'alpha_power'});
    
    % assign the appropriate equation of state
    if strcmp(medium.alpha_power, 'stokes')
        equation_of_state = 'stokes';
        medium.alpha_power = 2;
    else
        equation_of_state = 'absorbing';
        if medium.alpha_power == 1
            medium.alpha_power = LINEAR_ATTEN_DEF;
        end
        if abs(medium.alpha_power - 1) < LINEAR_ATTEN_TOLERANCE
            display('WARNING: Power law exponents close to 1 may increase numerical dispersion');
        end        
    end
else
    equation_of_state = 'linear';
end 

% check sensor fields
checkFieldNames(sensor, {'mask', 'directivity_pattern', 'directivity_angle', 'directivity_size',...
    'time_reversal_boundary_data', 'time_reversal_adapt_thresh'});
enforceFields(sensor, {'mask'});

% check for sensor directivity input
if ~isfield(sensor, 'directivity_angle')
    compute_directivity = false;
else
    compute_directivity = true;
end

% check for time reversal inputs
if isfield(sensor, 'time_reversal_boundary_data')
    if isfield(sensor, 'time_reversal_adapt_thresh')
        time_rev = 2;
    else
        time_rev = 1;
    end
else
    time_rev = 0;
end

% check source input for empty initial pressure
if isfield(source, 'p0')
    if isempty(source.p0) || ~sum(source.p0(:))
        source = rmfield(source, 'p0');
    end
end

% check source inputs
if ~isstruct(source)
    % allow an invalid or empty source input if computing time reversal    
    if ~time_rev
        error('Invalid source input');
    end
else
    checkFieldNames(source, {'p0', 'p', 'p_mask'});
end

% check for time varying source input and set source flag
if isfield(source, 'p') || isfield(source, 'p_mask')
    % if one time varying source parameter is given, enforce the other
    enforceFields(source, {'p', 'p_mask'});
    p_source = true;
else
    p_source = false;
end

% check kgrid for t_array existance and stability
if strcmp(kgrid.t_array, 'auto')
    if ~time_rev
        % create the time array
        [kgrid.t_array dt] = makeTime(kgrid, medium.sound_speed);
    else
        % throw error requesting for t_array
        error('t_array must be given explicitly in time reversal mode');
    end
else
    % assign dt
    dt = kgrid.t_array(2) - kgrid.t_array(1);
    
    % check the time steps are increasing
    if dt <= 0
        error('kgrid.t_array must be monotonically increasing');
    end
    
    % check the time array is evenly spaced
    if (kgrid.t_array(2:end) - kgrid.t_array(1:end-1)) ~= dt
        error('kgrid.t_array must be evenly spaced');
    end
end

% shorten commonly used field names
t_array = kgrid.t_array;
c = medium.sound_speed;
rho0 = medium.density;
   
% assign the default input parameters
cartesian_interp = CARTESIAN_INTERP_DEF;
create_log = CREATE_LOG_DEF;
data_cast = DATA_CAST_DEF;
display_mask = DISPLAY_MASK_DEF;
movie_args = MOVIE_ARGS_DEF;
movie_name = MOVIE_NAME_DEF;
plot_freq = PLOT_FREQ_DEF;
plot_layout = PLOT_LAYOUT_DEF;
plot_scale = PLOT_SCALE_DEF;
plot_sim = PLOT_SIM_DEF;
PML_inside = PML_INSIDE_DEF;
PML_x_alpha = PML_ALPHA_DEF;
PML_z_alpha = PML_ALPHA_DEF;
PML_x_size = PML_SIZE_DEF;
PML_z_size = PML_SIZE_DEF;
record_movie = RECORD_MOVIE_DEF;
smooth_c = SMOOTH_DEF;
smooth_p0 = SMOOTH_DEF;
smooth_rho = SMOOTH_DEF;

% replace defaults with user defined values if provided and check inputs    
if nargin < NUM_REQ_INPUT_VARIABLES
    error('Not enough input parameters');
elseif rem(nargin - NUM_REQ_INPUT_VARIABLES, 2)
    error('Optional input parameters must be given as param, value pairs');
elseif ~isempty(varargin)
    for input_index = 1:2:length(varargin)
        switch varargin{input_index}           
            case 'CartInterp'
                cartesian_interp = varargin{input_index + 1}; 
                if ~(strcmp(cartesian_interp, 'linear') || strcmp(cartesian_interp, 'nearest'))
                    error('Optional input CartInterp must be set to linear or nearest');
                end
            case 'CreateLog'
                create_log = varargin{input_index + 1}; 
                if ~islogical(create_log)
                    error('Optional input CreateLog must be Boolean');
                end
            case 'DataCast'
                data_cast = varargin{input_index + 1};
                if ~ischar(data_cast)
                    error('Optional input DataCast must be a string');
                end
            case 'DisplayMask'
                display_mask = varargin{input_index + 1};
                if ~all(size(display_mask) == size(kgrid.k))
                    error('Optional input DisplayMask must be the same size as the k-space grid');
                end
            case 'MovieArgs'
                movie_args = varargin{input_index + 1};  
            case 'MovieName'
                movie_name = varargin{input_index + 1};
                if ~ischar(movie_name)
                    error('Optional input MovieName must be a string');
                end                
            case 'PlotFreq'
                plot_freq = varargin{input_index + 1}; 
                if ~(numel(plot_freq) == 1 && isnumeric(plot_freq))
                    error('Optional input PlotFreq must be a single numerical value');
                end
            case 'PlotLayout'
                plot_layout = varargin{input_index + 1}; 
                if ~islogical(plot_layout)
                    error('Optional input PlotLayout must be Boolean');
                end                
            case 'PlotScale'
                plot_scale = varargin{input_index + 1};
                if ~(numel(plot_scale) == 2 && isnumeric(plot_scale))
                    error('Optional input PlotScale must be a 2 element numerical array');
                end
            case 'PlotSim'
                plot_sim = varargin{input_index + 1};
                if ~islogical(plot_sim)
                    error('Optional input PlotSim must be Boolean');
                end                
            case 'PMLAlpha'
                if length(varargin{input_index + 1}) > 2
                    error('Optional input PMLAlpha must be a 1 or 2 element numerical array');
                end                
                PML_x_alpha = varargin{input_index + 1}(1);
                PML_z_alpha = varargin{input_index + 1}(end);
            case 'PMLInside'
                PML_inside = varargin{input_index + 1};   
                if ~islogical(PML_inside)
                    error('Optional input PMLInside must be Boolean');
                end
            case 'PMLSize'
                if length(varargin{input_index + 1}) > 2
                    error('Optional input PMLSize must be a 1 or 2 element numerical array');
                end
                PML_x_size = varargin{input_index + 1}(1);
                PML_z_size = varargin{input_index + 1}(end);                
            case 'RecordMovie'
                record_movie = varargin{input_index + 1};    
                if ~islogical(record_movie)
                    error('Optional input RecordMovie must be Boolean');
                end
            case 'Smooth'
                if length(varargin{input_index + 1}) > 3 || ~islogical(varargin{input_index + 1})
                    error('Optional input Smooth must be a 1, 2 or 3 element Boolean array');
                end
                smooth_p0 = varargin{input_index + 1}(1);
                smooth_c = varargin{input_index + 1}(ceil((end + 1)/2));
                smooth_rho = varargin{input_index + 1}(end);   
            case 'TimeRev'
                % do nothing but allow the input name to exist for B.0.1
                % release compatability
                disp('WARNING: Optional input TimeRev has been deprecated');
            otherwise
                error(['Unknown optional input ' varargin{input_index}]);
        end
    end
end

% cleanup unused variables
clear *_DEF NUM_REQ_INPUT_VARIABLES;

% switch off layout plot in time reversal mode
plot_layout = plot_layout && ~time_rev;

% force visualisation if record_movie is true
if record_movie
    plot_sim = true;
end

% ensure p0 smoothing is switched off if p0 is empty
if ~isfield(source, 'p0')
    smooth_p0 = false;
end

% start log if required
if create_log
    diary([LOG_NAME '.txt']);
end

% enforce GPUmat compatability by using wrapper for GPUsingle
if strcmp(data_cast, 'GPUsingle');
    data_cast = 'kWaveGPUsingle';
end

% update command line status
if ~time_rev
    disp('Running k-space simulation...'); 
else
    disp('Running k-space time reversal...');
end

% check plot scaling
if ~time_rev && plot_sim
    
    % find the maximum input pressure amplitude
    if isfield(source, 'p0') && isfield(source, 'p')
        max_val = max([source.p0(:); source.p(:)]);
    elseif isfield(source, 'p0')
        max_val = max(source.p0(:));
    elseif isfield(source, 'p')
        max_val = max(source.p(:));
    else
        max_val = 0;
    end
    
    % check the plot scaling
    if max_val > PLOT_SCALE_WARNING*plot_scale(2) || PLOT_SCALE_WARNING*max_val < plot_scale(2)
        disp('  WARNING: visualisation plot scale may not be optimal for given source');
    end
    
    clear max_val;    
end

% check kgrid.t_array for stability given medium properties
if (numel(medium.sound_speed) > 1 || numel(medium.density) > 1) &&...
        (dt > DT_WARNING_CFL*max([kgrid.dz, kgrid.dx])/max(medium.sound_speed(:)))
    disp('  WARNING: time step may be too large for a stable simulation');
end

% =========================================================================
% CHECK AND PREPARE SENSOR MASK
% =========================================================================

% switch off Cartesian reorder flag
reorder_data = false;

% set usage to binary sensor mask
binary_sensor_mask = true;

% check if sensor mask is a binary grid or a set of interpolation points
if all(size(sensor.mask) == size(kgrid.k))
    
    % check the grid is binary
    if sum(sensor.mask(:)) ~= numel(sensor.mask) - sum(sensor.mask(:) == 0)
        error('sensor.mask must be a binary grid (numeric values must be 0 or 1)');
    end

else
   
    % extract Cartesian data from sensor mask
    sensor_x = sensor.mask(1, :);
    sensor_z = sensor.mask(2, :);
    
    % compute an equivalent sensor mask using nearest neighbour
    % interpolation, if time_rev = false and cartesian_interp = 'linear'
    % then this is only used for display, if time_rev = true or
    % cartesian_interp = 'nearest' this grid is used as the sensor.mask 
    [sensor.mask, order_index, reorder_index] = cart2grid(kgrid, sensor.mask);
   
    if ~time_rev && strcmp(cartesian_interp, 'nearest')
        % use the interpolated binary sensor mask but switch on Cartesian
        % reorder flag 
        reorder_data = true;
        
        % check if any duplicate points have been discarded
        num_discarded_points = length(sensor_x) - sum(sensor.mask(:));
        if num_discarded_points ~= 0
            disp(['  WARNING: ' num2str(num_discarded_points) ' duplicated sensor points discarded (nearest neighbour interpolation)']);
        end
    else
        % use the Cartesian points instead of the binary sensor mask
        binary_sensor_mask = false;
    end
        
    % reorder the sensor data in the order of the binary sensor.mask 
    if time_rev
        
        % append the reordering data
        new_col_pos = length(sensor.time_reversal_boundary_data(1,:)) + 1;
        sensor.time_reversal_boundary_data(:, new_col_pos) = order_index;

        % reorder p0 based on the order_index
        sensor.time_reversal_boundary_data = sortrows(sensor.time_reversal_boundary_data, new_col_pos);
        
        % remove the reordering data
        sensor.time_reversal_boundary_data = sensor.time_reversal_boundary_data(:, 1:new_col_pos - 1);
        
    end
end

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

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.z_size)); %#ok<*ASGLU>
disp(['  input grid size: ' num2str(kgrid.Nx) ' by ' num2str(kgrid.Nz) ' pixels (' num2str(kgrid.x_size*scale) ' by ' num2str(kgrid.z_size*scale) prefix 'm)']);
disp(['  maximum supported frequency: ' scaleSI( kgrid.k_max * min(c(:)) / (2*pi) ) 'Hz']);

% =========================================================================
% SMOOTH AND ENLARGE INPUT GRIDS
% =========================================================================
        
% smooth p0 if required, restoring the 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.Nz + 2*PML_z_size, kgrid.dz);
       
    % enlarge the grid of sound speed by extending the edge values into the
    % expanded grid 
    if numDim(c) == 2
        disp('  expanding sound speed grid...');
        c = expandMatrix(c, [PML_x_size, PML_z_size]);
    end

    % enlarge the grid of density by extending the edge values into the
    % expanded grid
    if numDim(rho0) == 2
        disp('  expanding density grid...');
        rho0 = expandMatrix(rho0, [PML_x_size, PML_z_size]);
    end
    
    % create indexes to allow the source input to be placed into the larger
    % simulation grid
    x1 = (PML_x_size + 1);
    x2 = kgrid.Nx - PML_x_size;
    z1 = (PML_z_size + 1);
    z2 = kgrid.Nz - PML_z_size;
    
    % enlarge the sensor mask
    sensor_mask_new = zeros(kgrid.Nz, kgrid.Nx);
    sensor_mask_new(z1:z2, x1:x2) = sensor.mask;
    sensor.mask = sensor_mask_new;
    
    % enlarge the display mask if given
    if ~strcmp(display_mask, 'sensor')
        display_mask_new = zeros(kgrid.Nz, kgrid.Nx);
        display_mask_new(z1:z2, x1:x2) = display_mask;
        display_mask = display_mask_new;
        clear display_mask_new;
    end
       
    % enlarge the source mask if given
    if p_source
        source_p_mask = zeros(kgrid.Nz, kgrid.Nx);
        source_p_mask(z1:z2, x1:x2) = source.p_mask;
        source.p_mask = source_p_mask;
        clear source_p_mask;
    end

    % enlarge the directivity angle if given
    if compute_directivity
        sensor_directivity_angle = zeros(kgrid.Nz, kgrid.Nx);
        sensor_directivity_angle(z1:z2, x1:x2) = sensor.directivity_angle;
        sensor.directivity_angle = sensor_directivity_angle;
        clear sensor_directivity_angle;
    end
    
    % clean up unused variables
    clear sensor_mask_new;
    
    % update command line status
    disp(['  computational grid size: ' num2str(kgrid.Nx) ' by ' num2str(kgrid.Nz) ' pixels']);
else
    % create indexes to place the source input exactly into the simulation
    % grid
    x1 = 1;
    x2 = kgrid.Nx;
    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) == 2
    disp('  smoothing sound speed distribution...');  
    c = smooth(kgrid, c);
end
    
% smooth rho0 distribution if required
if smooth_rho && numDim(rho0) == 2
    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;
z_sg = kgrid.z + kgrid.dz/2;

% interpolate the values of the density at the staggered grid locations
% where r1 = (x + dx/2, z) and r2 = (x, z + dz/2), replace values outside
% of the interpolation range with their original values
if numDim(rho0) == 2
    % rho0 is heterogeneous  
    rho0_r1 = interp2(kgrid.x, kgrid.z, rho0, x_sg, kgrid.z, '*linear');
    rho0_r1(isnan(rho0_r1)) = rho0(isnan(rho0_r1));
    rho0_r2 = interp2(kgrid.x, kgrid.z, rho0, kgrid.x, z_sg, '*linear');
    rho0_r2(isnan(rho0_r2)) = rho0(isnan(rho0_r2));
else
    % rho0 is homogeneous
    rho0_r1 = rho0;
    rho0_r2 = 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;
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);

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);

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_z = exp(-az_sg*dt/2);
abc_x_alt = exp(-ax*dt);
abc_z_alt = exp(-az*dt);

% define the modified first order k-space derivative operators and
% multiply by the staggered grid shift operators
ddx_k_shift = 1i*kgrid.kx.*sinc(c_max*dt*kgrid.k/2) .* exp(1i*kgrid.kx*kgrid.dx/2);
ddx_k_shift_min = 1i*kgrid.kx.*sinc(c_max*dt*kgrid.k/2) .* exp(-1i*kgrid.kx*kgrid.dx/2);
ddz_k_shift = 1i*kgrid.kz.*sinc(c_max*dt*kgrid.k/2) .* exp(1i*kgrid.kz*kgrid.dz/2);
ddz_k_shift_min = 1i*kgrid.kz.*sinc(c_max*dt*kgrid.k/2) .* exp(-1i*kgrid.kz*kgrid.dz/2);

% pre-shift variables used as frequency domain multipliers
ddx_k_shift = ifftshift(ddx_k_shift);
ddx_k_shift_min = ifftshift(ddx_k_shift_min);
ddz_k_shift = ifftshift(ddz_k_shift);
ddz_k_shift_min = ifftshift(ddz_k_shift_min);

% cleanup unused variables
clear ax* az* x0_min x0_max x_sg z0_min z0_max z_sg PML*;

% =========================================================================
% PREPARE ABSORPTION VARIABLES
% =========================================================================

% define the lossy derivative operators and proportionality coefficients
if strcmp(equation_of_state, 'absorbing')
            
    % convert the absorption coefficient to nepers.(rad/s)^-y.m^-1
    medium.alpha_coeff = db2neper(medium.alpha_coeff, medium.alpha_power);

    % compute the absorbing fractional Laplacian operator and coefficient
    absorb_nabla1 = (kgrid.k.*sinc(c_max*dt*kgrid.k/2)).^(medium.alpha_power-2); 
    absorb_nabla1(isinf(absorb_nabla1)) = 0;
    absorb_nabla1 = ifftshift(absorb_nabla1);
    absorb_tau = -2*medium.alpha_coeff*c_max^(medium.alpha_power - 1);
    absorb_param = absorb_tau.*absorb_nabla1.*rho0;
    
    % compute the dispersive fractional Laplacian operator and coefficient
    absorb_nabla2 = (kgrid.k.*sinc(c_max*dt*kgrid.k/2)).^(medium.alpha_power-1);  
    absorb_nabla2(isinf(absorb_nabla2)) = 0;
    absorb_nabla2 = ifftshift(absorb_nabla2);            
    absorb_eta = 2*medium.alpha_coeff*c_max^(medium.alpha_power)*tan(pi*medium.alpha_power/2);
    dispers_param = -absorb_eta.*absorb_nabla2;
    
    % cleanup unused variables
    clear absorb_nabla* absorb_tau absorb_eta;
    
elseif strcmp(equation_of_state, 'stokes')
    
    % convert the absorption coefficient to nepers.(rad/s)^-y.m^-1
    medium.alpha_coeff = db2neper(medium.alpha_coeff, medium.alpha_power);
    
    % compute the proportionality coefficient
    absorb_tau = -2*medium.alpha_coeff*c_max;
    absorb_param = absorb_tau.*rho0;
    
    % cleanup unused variables
    clear absorb_tau;
    
end

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

% create mask indices
sensor_mask_ind  = find(sensor.mask ~= 0);

% initialise threshold variable used for the adaptive bc
threshold_mask = 1;

% create storage and scaling variables
switch time_rev
    case 0
        
        % preallocate storage variables
        if binary_sensor_mask
            sensor_data = zeros(sum(sensor.mask(:)), length(t_array));
        else
            sensor_data = zeros(length(sensor_x), length(t_array));
        end
               
    case 2
        
        % preallocate the sensor variable
        p_sensor = zeros(1, kgrid.Nx*kgrid.Nz);
        
        % preallocate the threshold variable used for visualisation
        threshold_mask = zeros(kgrid.Nz, kgrid.Nx);

        % precompute a threshold index which defines whether there is any
        % boundary data above the threshold for that time step
        p_sensor_threshold_index  = abs(sensor.time_reversal_boundary_data) >= sensor.time_reversal_adapt_thresh;
        p_sensor_threshold_index  = any(p_sensor_threshold_index);
        
        % reverse the order of the threshold index as the boundary data is
        % later reversed
        p_sensor_threshold_index = fliplr(p_sensor_threshold_index);        
        
        % set pressure values below the threshold to zero
        sensor.time_reversal_boundary_data(abs(sensor.time_reversal_boundary_data) < sensor.time_reversal_adapt_thresh) = 0;     
        
end

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

% create the velocity and pressure variables
p = zeros(kgrid.Nz, kgrid.Nx);
ux_r1 = zeros(kgrid.Nz, kgrid.Nx);
uz_r2 = zeros(kgrid.Nz, kgrid.Nx);
    
% add p0 to p
if isfield(source, 'p0')
    p(z1:z2, x1:x2) = source.p0;
end    

% add p_source(t = t1) to p and create source variables
if p_source
          
    % create an empty matrix to place the source terms into
    ps = zeros(kgrid.Nz, kgrid.Nx);
    ps_index = find(source.p_mask ~= 0);
    
    % compute the scaling parameter
    if numel(c) == 1
        source_scale_factor = 2*c*dt./kgrid.dx;
    else
        source_scale_factor = 2.*c(ps_index).*dt./kgrid.dx;
    end
    
    % extract the initial source pressure and add to p
    ps(ps_index) = source.p(:, 1).*source_scale_factor;
    p = p + ps;
    
end

% compute rho(t = t1) using a linear (adiabatic) equation of state
rhox = 0.5*p./(c.^2);
rhoz = 0.5*p./(c.^2);
    
% 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;
    index_end = length(t_array);
end

% prepare data storage
if ~time_rev
    
    % precomputate the triangulation points used in gridDataFast
    if ~binary_sensor_mask && ~time_rev
        [zi del_tri] = gridDataFast(kgrid.z, kgrid.x, p, sensor_z, sensor_x);
    end    
    
    % store p(t = t1)
    if binary_sensor_mask
        sensor_data(:, 1) = p(sensor_mask_ind);
    else
        sensor_data(:, 1) = gridDataFast(kgrid.z, kgrid.x, p, sensor_z, sensor_x, del_tri);
    end
    
end

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

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

% plot the simulation layout
if plot_layout
    figure;
    subplot(2, 2, 1), imagesc(kgrid.x(1,:)*scale, kgrid.z(:,1)*scale, p, plot_scale);
    axis image;
    colormap(COLOR_MAP);
    title('Initial Pressure');
    subplot(2, 2, 2), imagesc(kgrid.x(1,:)*scale, kgrid.z(:,1)*scale, sensor.mask, [-1 1]);
    axis image;
    title('Sensor Mask');
    subplot(2, 2, 3), imagesc(kgrid.x(1,:)*scale, kgrid.z(:,1)*scale, c, [0.5*min(c(:)), max(c(:))/0.5]);
    axis image;
    title('Sound Speed');
    subplot(2, 2, 4), imagesc(kgrid.x(1,:)*scale, kgrid.z(:,1)*scale, rho0, [0.5*min(rho0(:)), max(rho0(:))/0.5]);
    axis image;
    title('Density');   
    xlabel(['(All axes in ' prefix 'm)']);
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
    
    % define frame index
    frame_index = 1;
    
    %  save the first frame if in forward mode
    if ~time_rev
    
        % scale the image data from 1 -> length(COLOR_MAP) using the
        % plot_scale used for the simulation animation
        p0_mov = double(p);
        p0_mov(p0_mov > plot_scale(2)) = plot_scale(2);
        p0_mov(p0_mov < plot_scale(1)) = plot_scale(1);
        
        % add display mask onto plot, in adaptive time reversal mode only
        % display the active sensor elements
        if strcmp(display_mask, 'sensor')
            p0_mov(threshold_mask.*sensor.mask == 1) = plot_scale(2);
        else
            p0_mov(display_mask ~= 0) = plot_scale(2);
        end        

        % apply the scaling
        p0_mov = (length(COLOR_MAP) - 1).*(p0_mov + abs(plot_scale(1)))./(plot_scale(2) + abs(plot_scale(1))) + 1;

        % save the movie frame
        movie_frames(frame_index) = im2frame(p0_mov, COLOR_MAP);

        % update frame index
        frame_index  = frame_index  + 1;

        % cleanup unused variables
        clear p0_mov;
    end
end

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

% cast variables used in time loop
if ~strcmp(data_cast, 'off')
    
    % update command line status
    disp(['  casting variables to ' data_cast ' type...']);

    % create list of variable names to cast
    cast_variables = {'ddx_k_shift', 'ddx_k_shift_min', 'ddz_k_shift', 'ddz_k_shift_min',...
        'abc_x', 'abc_z', 'abc_x_alt', 'abc_z_alt',...
        'ux_r1', 'uz_r2', 'p',...
        'rho0', 'rho0_r1', 'rho0_r2', 'rhox', 'rhoz',...
        'sensor_mask_ind', 'c', 'dt'};
    
    % additional variables only used in time reversal
    if ~time_rev
        cast_variables = [cast_variables, {'sensor_data'}];
    elseif time_rev
        cast_variables = [cast_variables, {'sensor.time_reversal_boundary_data'}];
    elseif time_rev == 2
        cast_variables = [cast_variables, {'sensor.time_reversal_boundary_data', 'threshold_mask', 'p_sensor', 'p_sensor_threshold_index'}];
    end
    
    % additional variables only used if the medium is absorbing
    if strcmp(equation_of_state, 'absorbing')
        cast_variables = [cast_variables, {'absorb_param', 'dispers_param'}];
    elseif strcmp(equation_of_state, 'stokes')
        cast_variables = [cast_variables, {'absorb_param'}];
    end
    
    % additional variables only used if there is a time varying source term
    if p_source
        cast_variables = [cast_variables, {'ps', 'ps_index', 'source.p', 'source_scale_factor'}];
    end    
    
    % cast variables
    for cast_index = 1:length(cast_variables)
        eval([cast_variables{cast_index} ' = ' data_cast '(' cast_variables{cast_index} ');']);
    end
end

% =========================================================================
% 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 = fft2(p);

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

    % enforce time reversal bounday condition
    switch time_rev
        
        % enforce a time-varying Dirichlet boundary condition
        case 1
            
            % load pressure value on boundary
            p(sensor_mask_ind) = sensor.time_reversal_boundary_data(:, t_index);
            
            % update p_k
            p_k = fft2(p);
            
            % compute rhox and rhoz using an adiabatic equation of state
            rhox_mod = 0.5*p./(c.^2);
            rhoz_mod = 0.5*p./(c.^2);
            rhox(sensor_mask_ind) = rhox_mod(sensor_mask_ind);
            rhoz(sensor_mask_ind) = rhoz_mod(sensor_mask_ind);
            
        % enforce a thresholded Dirichlet boundary condition
        case 2
            
            % first check if there are any values above the threshold
            if p_sensor_threshold_index(t_index)

                % place the boundary pressure values into a larger grid
                p_sensor(sensor_mask_ind) = sensor.time_reversal_boundary_data(:, t_index);
                p_sensor_index = find(p_sensor ~= 0);
                
                % apply adaptive boundary condition
                p(p_sensor_index) = p_sensor(p_sensor_index);
                
                % update p_k
                p_k = fft2(p);
                
                % compute rhox and rhoz using an adiabatic equation of state
                rhox_mod = 0.5*p./(c.^2);
                rhoz_mod = 0.5*p./(c.^2);
                rhox(p_sensor_index) = rhox_mod(p_sensor_index);
                rhoz(p_sensor_index) = rhoz_mod(p_sensor_index);
                
                % update threshold mask used for visualisation
                threshold_mask(:) = 0;
                threshold_mask(p_sensor_index) = 1; %#ok<AGROW>

            else
                threshold_mask(:) = 0;
            end
    end

    % calculate ux and uz at t + dt/2 using dp/dx and dp/dz at t
    ux_r1 = abc_x .* (  abc_x.*ux_r1 - dt./rho0_r1 .* real(ifft2(ddx_k_shift .* p_k))  );
    uz_r2 = abc_z .* (  abc_z.*uz_r2 - dt./rho0_r2 .* real(ifft2(ddz_k_shift .* p_k))  );

    % calculate dux/dx and duz/dz at t + dt/2
    duxdx_k = ddx_k_shift_min .* fft2(ux_r1);
    duzdz_k = ddz_k_shift_min .* fft2(uz_r2);     

    % calculate rhox and rhoz at t + dt
    rhox = abc_x_alt .* (  rhox - dt.*rho0 .* real(ifft2(duxdx_k))  );
    rhoz = abc_z_alt .* (  rhoz - dt.*rho0 .* real(ifft2(duzdz_k))  );     
    
    switch equation_of_state
        case 'linear'
            % calculate p using an adiabatic equation of state
            p = c.^2.*(rhox + rhoz);
        case 'absorbing'
            % calculate p using an absorbing equation of state          
            p = c.^2.*(  (rhox + rhoz) + real(ifft2( absorb_param.*(duxdx_k + duzdz_k) + dispers_param.*fft2(rhox + rhoz) ))  );
        case 'stokes'
            % calculate p using stokes' equation of state
            p = c.^2.*(  (rhox + rhoz) + absorb_param .* real(ifft2(duxdx_k + duzdz_k))  );
    end

    % add in the source term
    if p_source
        % update the pressure term with the source pressure
        ps(ps_index) = source.p(:, t_index).*source_scale_factor;
        p = p + ps;
        
        % update the acoustic density values using a linear equation of
        % state
        rhox = rhox + 0.5*ps./(c.^2);
        rhoz = rhoz + 0.5*ps./(c.^2);
    end

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

    % extract required data
    if ~time_rev
        if ~compute_directivity
            if binary_sensor_mask
                sensor_data(:, t_index) = p(sensor_mask_ind);
            else
                sensor_data(:, t_index) = gridDataFast(kgrid.z, kgrid.x, p, sensor_z, sensor_x, del_tri);
            end
        else
           sensor_data(:, t_index) = directionalResponse(kgrid, sensor, p_k); 
        end
    end

    % plot data if required
    if plot_sim && rem(t_index, plot_freq) == 0

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

        % ensure p is cast as a CPU variable
        p_plot = double(p);

        % add display mask onto plot, in adaptive time reversal mode only
        % display the active sensor elements
        if strcmp(display_mask, 'sensor')
            p_plot(double(threshold_mask).*sensor.mask == 1) = plot_scale(2);
        else
            p_plot(display_mask ~= 0) = plot_scale(2);
        end

        % update plot
        imagesc(kgrid.x(1,:)*scale, kgrid.z(:,1)*scale, p_plot, plot_scale);
        colormap(COLOR_MAP);
        ylabel(['z-position [' prefix 'm]']);
        xlabel(['x-position [' prefix 'm]']);
        axis image;
        drawnow;

        % save movie frame if required
        if record_movie

            % scale the image data from 1 -> length(COLOR_MAP) using the
            % plot_scale used for the simulation animation
            p_plot(p_plot > plot_scale(2)) = plot_scale(2);
            p_plot(p_plot < plot_scale(1)) = plot_scale(1);
            p_plot = (length(COLOR_MAP) - 1).*(p_plot + abs(plot_scale(1)))./(plot_scale(2) + abs(plot_scale(1))) + 1;

            % save the movie frame
            movie_frames(frame_index) = im2frame(p_plot, COLOR_MAP);

            % update frame index
            frame_index  = frame_index  + 1;

        end
    end
end

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

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

% cast pressure variable back to double if required and save the final
% pressure to return to user
if ~strcmp(data_cast, 'off')
    sensor_data = double(sensor_data);
    p_final = double(p(z1:z2, x1:x2));
else
    p_final = p(z1:z2, x1:x2);
end

% reorder the sensor points if a binary sensor mask was used for Cartesian
% sensor mask nearest neighbour interpolation
if reorder_data
    
    % update command line status
    disp('  reordering Cartesian measurement data...');
    
    % 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

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

% store the final movie frame and then save the movie
if record_movie
        
    % update command line status
    disp('  saving movie file...');    
    
    % scale the image data from 1 -> length(COLOR_MAP) using the
    % plot_scale used for the simulation animation
    p0_mov = double(p);
    p0_mov(p0_mov > plot_scale(2)) = plot_scale(2);
    p0_mov(p0_mov < plot_scale(1)) = plot_scale(1);

    % add display mask onto plot, in adaptive time reversal mode only
    % display the active sensor elements
    if strcmp(display_mask, 'sensor')
        p0_mov(threshold_mask.*sensor.mask == 1) = plot_scale(2);
    else
        p0_mov(display_mask ~= 0) = plot_scale(2);
    end       
    
    % apply the scaling
    p0_mov = (length(COLOR_MAP) - 1).*(p0_mov + abs(plot_scale(1)))./(plot_scale(2) + abs(plot_scale(1))) + 1;

    % save the movie frame
    movie_frames(frame_index) = im2frame(p0_mov, COLOR_MAP);
          
    % check if user provided a compression setting
    search_found = false;
    if ~isempty(movie_args)
        search_index = 1;
        while ~search_found && search_index <= length(movie_args)
            if strcmp(movie_args{search_index}, 'Compression')
                search_found = true;
            end
            search_index = search_index + 1;
        end
    end
    
    if search_found
        % use compression setting provided by user
        movie2avi(movie_frames, movie_name, movie_args{:});
    else
        % check if on UNIX or 64 bit computer
        if sum(strfind(computer, '64'))
            movie2avi(movie_frames, movie_name, 'Compression', MOVIE_COMP_64B, movie_args{:});
        elseif sum(strfind(computer, 'GLNX'))
            movie2avi(movie_frames, movie_name, 'Compression', MOVIE_COMP_LNX, movie_args{:});
        else
            movie2avi(movie_frames, movie_name, 'Compression', MOVIE_COMP_WIN, movie_args{:});
        end
    end
end

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

% switch off log
if create_log
    diary off;
end