function sensor_data = kspaceFirstOrder1D(p0, kgrid, c, rho, t_array, sensor_mask, varargin) 
%KSPACEFIRSTORDER1D     1D time-domain simulation of wave propagation.
%
% DESCRIPTION:
%       kspaceFirstOrder1D simulates the time-domain propagation of linear
%       compressional waves through a one-dimensional homogeneous or
%       heterogeneous acoustic medium defined by c and rho given the
%       initial pressure distribution p0. The size and discretisation of
%       the acoustic domain are defined by the k-space grid structure
%       kgrid. At each time step the pressure at the positions defined by
%       sensor_mask are recorded and stored.
%
%       The computation is based on a first-order k-space model which
%       allows a heterogeneous sound speed and density. An absorbing
%       boundary condition (in this case a perfectly matched layer) 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 (PML). 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 p0)
%       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 1 by N matrix
%
%       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
%       setting the optional input 'TimeRev' to true. This enforces the
%       pressure given by p0 as a time varying Dirichlet boundary condition
%       over the sensor mask. In this mode, the input pressure p0 must be
%       indexed as p0(sensor position, time). If sensor_mask is given as a
%       set of Cartesian coordinates then p0 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 p0 must be given as an
%       array ordered using MATLAB's standard column-wise linear matrix
%       indexing.
%
%       To run the inverse crime, the sensor_data returned when 'TimeRev'
%       is set to false can be used without modification as p0 when
%       'TimeRev' is set to true.
%
% USAGE:
%       sensor_data = kspaceFirstOrder1D(p0, kgrid, c, rho, t_array, sensor_mask)
%       sensor_data = kspaceFirstOrder1D(p0, kgrid, c, rho, t_array, sensor_mask, ...) 
%
% INPUTS:
%       p0          - map of the initial pressure within the medium over
%                     the discretisation given by kgrid (or the time
%                     varying pressure across the sensor mask if 'TimeRev'
%                     is set to true)  
%       kgrid       - kspace grid structure returned by makeGrid.m
%       c / rho     - maps of the sound speed and density over the medium
%                     discretisation given by kgrid. For homogeneous media,
%                     c and rho may be given as a single value.
%       t_array     - evenly spaced array of time values [s] (t_array can
%                     alternatively be set to 'auto' to automatically
%                     generate the array using makeTime)
%       sensor_mask - binary grid or a set of Cartesian points where the
%                     pressure is recorded at each time step
%
% OPTIONAL INPUTS:
%       Optional 'string', value pairs that may be used to modify the
%       default computational settings.
%
%       'AdaptThresh' - Adaptive boundary condition threshold used when
%                     'TimeRev' is set to 3 (default = 0.005).
%       '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').     
%       '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 fft and ifft 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').
%       '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 plot scaling
%                     (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. To
%                     remove the PML, set the appropriate PMLAlpha to zero
%                     rather than forcing the PML to be of zero size
%                     (default = 20). 
%       'MovieName' - Name of the movie produced when 'RecordMovie' is set
%                     to true (default = 'date-time-kspacemovie').
%       'MovieArgs' - Settings for movie2avi. Parameters must be given as
%                     {param, value, ...} pairs within a cell array
%                     (default = {}).
%       'RecordMovie' - Boolean controlling whether the displayed image
%                     frames are captured using getframe 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).
%       'TimeRev'   - Boolean controlling whether the code is used in time
%                     reversal mode. If set to true (or 1), the time
%                     reversal is computed by enforcing the pressure values
%                     given by p0 over the sensor surface at each time step
%                     (conventional time reversal). If set to 2, the time
%                     reversal is computed by introducing the pressure
%                     values given by p0 over the sensor surface as a
%                     source term. If set to 3, the time reversal is
%                     computed by using an adaptive boundary condition with
%                     the threshold set by 'AdaptThresh' (default = false).
%
% OUTPUTS:
%       sensor_data - array of pressure time-series recorded at the sensor
%                     positions given by sensor_mask
%
% ABOUT:
%       author      - Bradley Treeby and Ben Cox
%       date        - 22nd April 2009
%       last update - 20th July 2009
%       
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
%
% See also fft, ifft, getframe, kspaceFirstOrder2D, kspaceFirstOrder3D,
% makeGrid, makeTime, movie2avi, smooth, unmaskSensorData 

% KNOWN ISSUES:
%       - status bar does not always update correctly in time reversal mode
%         for earlier versions of matlab
%       - progress bar does not close using 'close all' if ctrl+c forced
%         break is used to exit the simulation
%       - movie_frames storage variable is not preallocated
%       - indeo3 and indeo5 video compression fails on some PCs

% start the timer
tic;
time = clock;

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

% minimum number of input variables
NUM_REQ_INPUT_VARIABLES = 6;

% optional input defaults
ADAPTIVE_BC_THRESHOLD_DEF = 0.005;
CARTESIAN_INTERP_DEF = 'linear';
DATA_CAST_DEF = 'off';
MOVIE_ARGS_DEF = {};
MOVIE_NAME_DEF = [date '-' num2str(time(4)) '-' num2str(time(5)) '-kspaceFirstOrder1D'];
PLOT_FREQ_DEF = 10;
PLOT_LAYOUT_DEF = false;
PLOT_SCALE_DEF = [-1.1 1.1];
PLOT_SIM_DEF = true;
PML_SIZE_DEF = 20;
PML_ALPHA_DEF = 4;
PML_INSIDE_DEF = true;
RECORD_MOVIE_DEF = false;
SMOOTH_DEF = true;
TIME_REV_DEF = false;

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

% =========================================================================
% EXTRACT OPTIONAL INPUTS
% =========================================================================

% assign default input parameters
adaptive_bc_threshold = ADAPTIVE_BC_THRESHOLD_DEF;
cartesian_interp = CARTESIAN_INTERP_DEF;
data_cast = DATA_CAST_DEF;
movie_args = MOVIE_ARGS_DEF;
movie_name = MOVIE_NAME_DEF;
plot_freq = PLOT_FREQ_DEF;
plot_scale = PLOT_SCALE_DEF;
plot_sim = PLOT_SIM_DEF;
plot_layout = PLOT_LAYOUT_DEF;
PML_inside = PML_INSIDE_DEF;
PML_x_alpha = PML_ALPHA_DEF;
PML_x_size = PML_SIZE_DEF;
record_movie = RECORD_MOVIE_DEF;
smooth_c = SMOOTH_DEF;
smooth_p0 = SMOOTH_DEF;
smooth_rho = SMOOTH_DEF;
time_rev = TIME_REV_DEF;

% replace with user defined values if provided
if nargin < NUM_REQ_INPUT_VARIABLES
    error('Not enough input parameters');
elseif rem(nargin, 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 'AdaptThresh'
                adaptive_bc_threshold = varargin{input_index + 1};
            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 'DataCast'
                data_cast = varargin{input_index + 1};     
            case 'MovieArgs'
                movie_args = varargin{input_index + 1};  
            case 'MovieName'
                movie_name = varargin{input_index + 1};                 
            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'             
                PML_x_alpha = varargin{input_index + 1}(1);            
            case 'PMLInside'
                PML_inside = varargin{input_index + 1};
                if ~islogical(PML_inside)
                    error('Optional input PMLInside must be Boolean');
                end                
            case 'PMLSize'
                PML_x_size = varargin{input_index + 1}(1);     
            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'
                time_rev = varargin{input_index + 1};
              
            otherwise
                error('Unknown optional input');
        end
    end
end

% cleanup unused variables
clear *_DEF NUM_REQ_INPUT_VARIABLES time;

% 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

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

% check scaling input
if ~time_rev && plot_sim
    if max(p0(:)) > 10*plot_scale(2) || 10*max(p0(:)) < plot_scale(2)
        disp('  WARNING: visualisation plot scale may not be optimal for given p0');
    end
end

% =========================================================================
% SETUP TIME VARIABLE
% =========================================================================

% automatically create a suitable time array if required
if strcmp(t_array, 'auto')
    if ~time_rev
        % create the time array
        [t_array dt] = makeTime(kgrid, max(c(:)));
    else
        % throw error requesting for t_array
        error('t_array must be given explicitly in time reversal mode');
    end
else
    % extract the time step from the input data
    dt = t_array(2) - t_array(1); 
    
    % check for stability
    if (numel(c) > 1 || numel(rho) > 1) && (dt > 0.5*kgrid.dx/max(c(:)))
        disp('  WARNING: time step may be too large for a stable simulation');
    end
end

% setup the time index variable
if ~time_rev
    index_start = 1;
    index_step = 1;
    index_end = length(t_array);  
else
    index_start = length(t_array);
    index_step = -1;
    index_end = 1;    
end

% =========================================================================
% CHECK SENSOR MASK INPUT
% =========================================================================

% 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 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;
    
    % compute an equivalent sensor mask using nearest neighbour
    % interpolation, 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_x);
   
    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 p0 input data in the order of the binary sensor_mask 
    if time_rev
        
        % append the reordering data
        new_col_pos = length(p0(1,:)) + 1;
        p0(:, new_col_pos) = order_index;

        % reorder p0 based on the order_index
        p0 = sortrows(p0, new_col_pos);
        
        % remove the reordering data
        p0 = p0(:, 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))]);
disp(['  input grid size: ' num2str(kgrid.Nx) ' pixels (' scaleSI(kgrid.x_size) 'm)']);

% =========================================================================
% PREPARE COMPUTATIONAL GRIDS
% =========================================================================

% smooth p0 distribution if required restoring the maximum magnitude
if smooth_p0 && ~time_rev
    disp('  smoothing p0 distribution...');  
    p0 = smooth(p0, kgrid, 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);
    
    % 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;
    
    % enlarge the grid of sound speed by extending the edge values into the
    % expanded grid 
    if numel(c) > 1
        disp('  expanding sound speed grid...');
        c_new = ones(kgrid.Nx, 1);
        c_new(x1:x2) = c;  
        c_new(1:x1-1) = c(1);
        c_new(x2+1:end) = c(end); 
        c = c_new;
    end
    
    if numel(rho) > 1
        disp('  expanding density grid...');
        rho_new = ones(kgrid.Nx, 1);
        rho_new(x1:x2) = rho;  
        rho_new(1:x1-1) = rho(1);
        rho_new(x2+1:end) = rho(end); 
        rho = rho_new;
    end

    % enlarge the sensor mask
    sensor_mask_new = zeros(kgrid.Nx, 1);
    sensor_mask_new(x1:x2) = sensor_mask;
    sensor_mask = sensor_mask_new;
    
    % clean up unused variables
    clear c_new rho_new sensor_mask_new;
    
    % update command line status
    disp(['  computation grid size: ' num2str(kgrid.Nx) ' pixels']);
else
    % create indexes to place the source input exactly into the simulation
    % grid
    x1 = 1;
    x2 = kgrid.Nx;
end

% select reference sound speed based on heterogeneity map
c0 = max(c);

% smooth c distribution if required
if smooth_c && (sum(size(c)) > 2)
    disp('  smoothing c distribution...');     
    c = smooth(c, kgrid);
end
    
% smooth rho distribution if required
if smooth_rho && (sum(size(rho)) > 2)
    disp('  smoothing rho distribution...');     
    rho = smooth(rho, kgrid);
end

% =========================================================================
% PREPARE STAGGERED DENSITY AND PML GRIDS
% =========================================================================

% create the staggered grid
x_sg = kgrid.x + kgrid.dx/2;

% interpolate the values of the density at the staggered grid locations
if length(rho) > 1
    % rho is heterogeneous
    rho_sg = interp1(kgrid.x, rho, x_sg, '*linear');
    
    % set values outside of the interpolation range to original values 
    rho_sg(isnan(rho_sg)) = rho(isnan(rho_sg)); 
else
    % rho is homogeneous
    rho_sg = rho;
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;

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

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

% precompute shift and absorbing boundary condition operators
x_shift = exp(1i*kgrid.kx*kgrid.dx/2);
x_shift_min = exp(-1i*kgrid.kx*kgrid.dx/2);
abc_x = exp(-ax_sg*dt/2);
abc_x_alt = exp(-ax*dt);

% clean up unused variables
clear ax* x0_min x0_max x_sg PML*;

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

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

% 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
        
        % extract the sound speed at sensor location
        if length(c) > 1
            c_sens = c(sensor_mask ~= 0);
        else
            c_sens = c;
        end

        % assign source term scale values
        time_rev_x_scale = (2*c_sens*dt)/(2*kgrid.dx);
        
        % cleanup unused variables
        clear c_sens;
        
    case 3
        
        % preallocate the sensor variable
        p_sensor = zeros(1, kgrid.Nx);

        % precompute the threshold index
        p_sensor_threshold_index  = abs(p0) >= adaptive_bc_threshold;
        p_sensor_threshold_index  = sum(p_sensor_threshold_index);
        
        % set pressure values below threshold to zero
        p0(abs(p0) < adaptive_bc_threshold) = 0;        

end

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

% set the initial pressure and velocity
p = zeros(kgrid.Nx, 1);
u = zeros(kgrid.Nx, 1);

% define the modified first order k-space derivative operator
ddx_k = 1i*kgrid.kx.*sinc(c0*dt*kgrid.k/2);

% set initial pressure if in forward mode
if ~time_rev
    p(x1:x2) = p0;
end

% pre-shift variables used as frequency domain multipliers
ddx_k = ifftshift(ddx_k);
x_shift = ifftshift(x_shift);
x_shift_min = ifftshift(x_shift_min);

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

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

% plot the simulation layout
if plot_layout
    figure;
    p0_plot = zeros(kgrid.Nx, 1);
    p0_plot(x1:x2) = p0;
    subplot(2, 2, 1), plot(kgrid.x*scale, p0_plot);
    axis tight;
    title('Initial Pressure');
    subplot(2, 2, 2), plot(kgrid.x*scale, sensor_mask);
    axis tight;
    title('Sensor Mask');
    subplot(2, 2, 3), plot(kgrid.x*scale, c);
    axis tight;
    title('Sound Speed');
    subplot(2, 2, 4), plot(kgrid.x*scale, rho);
    axis tight;
    title('Density');
    xlabel(['(All horizontal axes in ' prefix 'm)']);
    clear p0_plot;
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;
    
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...']);

    % cast computation variables to data_cast type
    eval(['x_shift = ' data_cast '(x_shift);']);
    eval(['x_shift_min = ' data_cast '(x_shift_min);']);
    eval(['u = ' data_cast '(u);']);
    eval(['p = ' data_cast '(p);']);
    eval(['ddx_k = ' data_cast '(ddx_k);']);
    eval(['abc_x = ' data_cast '(abc_x);']);
    eval(['abc_x_alt = ' data_cast '(abc_x_alt);']);
    eval(['dt = ' data_cast '(dt);']);
    eval(['rho = ' data_cast '(rho);']);
    eval(['rho_sg = ' data_cast '(rho_sg);']);
    eval(['c = ' data_cast '(c);']);

    % cast storage variables to data_cast type
    eval(['sensor_mask_ind  = ' data_cast '(sensor_mask_ind);']);
    if ~time_rev
        eval(['sensor_data = ' data_cast '(sensor_data);']);
    else
        eval(['p0  = ' data_cast '(p0);']);  
    end
    
    if time_rev == 2
        eval(['time_rev_x_scale = ' data_cast '(time_rev_x_scale);']);           
    end
    
    if time_rev == 3 
        eval(['p_sensor  = ' data_cast '(p_sensor);']);
        eval(['p_sensor_threshold_index  = ' data_cast '(p_sensor_threshold_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;

% set adaptive boundary condition loop flag
if time_rev == 3
    null_boundary_condition = true;
    disp('  skipping time steps with no boundary input...');
else
    null_boundary_condition = false;
end

for t_index = index_start:index_step:index_end
    
    % enforce time reversal bounday condition
    switch time_rev
        case 1
            p(sensor_mask_ind) = p0(:, t_index);
        case 2
            p(sensor_mask_ind) = p(sensor_mask ~= 0) + time_rev_x_scale.*p0(:, t_index);    
        case 3
            
            % check if any values above the threshold were found
            if p_sensor_threshold_index(t_index)

                % update the boundary condition parameter
                null_boundary_condition = false;

                % place the boundary pressure values into a larger grid
                p_sensor(sensor_mask_ind) = p0(:, t_index);
                p_sensor_index = find(p_sensor ~= 0);
                
                % apply adaptive boundary condition
                p(p_sensor_index) = p_sensor(p_sensor_index);
                
            end
    end
    
    % skip loop in time reversal mode 
    if ~null_boundary_condition

        % calculate dp(t+dt)/dx from p(t+dt)
        dpdx = real(ifft(  ddx_k .* fft(p) .* x_shift));

        % calculate u(t+dt/2) from u(t-dt/2) and dp(t)/dx
        u = real(abc_x .* (  abc_x.*u - dt./rho_sg.*dpdx));

        % calculate du(t+dt/2)/dx from u(t+dt/2)
        dudx = real(ifft(  ddx_k .* fft(u) .* x_shift_min  ));  

        % calculate p(t+dt) from p(t) and du(t+dt/2)/dx 
        p = real(abc_x_alt .* (  p - dt*rho.*c.^2.*dudx)); 
         
        % extract required data
        if ~time_rev
            p_array = reshape(p, 1, []);
            if binary_sensor_mask
                sensor_data(:, t_index) = p_array(sensor_mask_ind);
            else
                sensor_data(:, t_index) = interp1(kgrid.x, p, sensor_x);
            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;
            
            % update plot
            plot(kgrid.x*scale, double(p));
            xlabel(['x-position [' prefix 'm]']);
            set(gca,'YLim',plot_scale);
                        
            % 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
            
            % force plot update
            drawnow;
        end    
    end
end

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

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

% cast pressure variable back to double if required
if ~strcmp(data_cast, 'off')
    sensor_data = double(sensor_data);
end

% reorder the sensor points if 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

% save movie if required
if record_movie
    
    % update command line status
    disp('  saving movie file...');
    
    % 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)]);