% DESCRIPTION:
%       subscript to check input structures and optional input parameters
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 21st December 2010
%       last update - 11th February 2011
%       
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009, 2010, 2011 Bradley Treeby and Ben Cox

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
% 
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details. 
% 
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% =========================================================================
% CHECK INPUT COMPATABILITY
% =========================================================================

% set exit flag
early_exit = false;

% 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 ~strcmp(class(kgrid), 'kWaveGrid') && ~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 ((numDim(kgrid.k) == 1) && all(size(p0_old) == [kgrid_old.Nx, 1])) || ...
            ((numDim(kgrid.k) == 2) && all(size(p0_old) == [kgrid_old.Nz, kgrid_old.Nx])) || ...
            ((numDim(kgrid.k) == 3) && all(size(p0_old) == [kgrid_old.Nz, kgrid_old.Nx, kgrid_old.Ny]))
        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 = eval([MFILE '(kgrid, medium, source, sensor, varargin{3:end})']);
    else
        sensor_data = eval([MFILE '(kgrid, medium, source, sensor)']);
    end
    
    % set exit flag
    early_exit = true;
    return
end 

% =========================================================================
% CHECK MEDIUM STRUCTURE INPUTS
% =========================================================================

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

% allow the density field to be blank if the medium is homogeneous
if ~isfield(medium, 'density') && numel(medium.sound_speed) == 1
    user_medium_density_input = false;
    medium.density = 1;
else
    enforceFields(medium, {'density'});
    user_medium_density_input = true;
end

% 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')
        % debugging option to ensure the encoded absorption matches stokes
        % equation for y = 2
        equation_of_state = 'stokes';
        medium.alpha_power = 2;
    else
        % check y is within 0 to 3
        if medium.alpha_power > 3 || medium.alpha_power < 0
            error('medium.alpha_power must be within 0 and 3');
        end
        
        % display warning if y is close to 1 and the dispersion term has
        % not been set to zero
        if ~(isfield(medium, 'alpha_mode') && strcmp(medium.alpha_mode, 'no_dispersion'))
            if medium.alpha_power == 1
                error('The power law dispersion term in the equation of state is not valid for medium.alpha_power = 1. This error can be avoided by choosing a power law exponent close to, but not exactly, 1. If modelling acoustic absorption for medium.alpha_power = 1 is important and modelling dispersion is not critical, this error can also be avoided by setting medium.alpha_mode to ''no_dispersion''.');
            end
        end
        equation_of_state = 'absorbing';
    end
    
    % check the absorption mode input
    if isfield(medium, 'alpha_mode')
        if ~ischar(medium.alpha_mode) || (~strcmp(medium.alpha_mode, 'no_absorption') && ~strcmp(medium.alpha_mode, 'no_dispersion'))
            error('medium.alpha_mode must be set to ''no_absorption'' or ''no_dispersion''');
        end
    end    
    
    % check the absorption filter input
    if isfield(medium, 'alpha_filter') && ~all(size(medium.alpha_filter) == size(kgrid.k))
        error('medium.alpha_filter must be the same size as the simulation grid');
    end
    
    % check the absorption sign input
    if isfield(medium, 'alpha_sign') && (~isnumeric(medium.alpha_sign) || numel(medium.alpha_sign) > 2)
        error('medium.alpha_sign must be given as a 2 element numerical array')
    end    
else
    equation_of_state = 'lossless';
end

if isfield(medium, 'BonA')
    nonlinear = true;
else
    nonlinear = false;
end 

% =========================================================================
% CHECK SENSOR STRUCTURE INPUTS
% =========================================================================

% check sensor fields
if isempty(sensor)
    use_sensor = false;
    time_rev = false;
    compute_directivity = false;
else
    use_sensor = true;
    
    % check if sensor is a transducer, otherwise check inputs
    if ~strcmp(class(sensor), 'kWaveTransducer')
        if kgrid.dim == 2

            % check field names including the directivity inputs
            checkFieldNames(sensor, {'mask', 'directivity_pattern', 'directivity_angle', 'directivity_size',...
                'time_reversal_boundary_data', 'frequency_response'});

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

        else
            % check field names without directivity inputs (not supported in 1 or 3D)
            checkFieldNames(sensor, {'mask', 'time_reversal_boundary_data', 'frequency_response'});
        end
        enforceFields(sensor, {'mask'});

        % check for time reversal inputs
        if isfield(sensor, 'time_reversal_boundary_data')
            time_rev = true;
        else
            time_rev = false;
        end
    else

        % check for time reversal inputs
        if isfield(sensor, 'time_reversal_boundary_data')
            % don't allow time reversal with a transducer just yet
            error('Transducer inputs not yet supported with time reversal');
        else
            time_rev = 0;
        end
        
        if kgrid.dim == 2
            % switch off directivity flag if a transducer is given
            compute_directivity = false;
        end
    end
end

% check for directivity inputs with time reversal
if kgrid.dim == 2 && compute_directivity && time_rev
    disp('WARNING: sensor directivity fields are not used for time reversal');
end

% =========================================================================
% CHECK SOURCE STRUCTURE INPUTS
% =========================================================================

% 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
    switch kgrid.dim
        case 1
            checkFieldNames(source, {'p0', 'p', 'p_mask', 'ux', 'u_mask'});
        case 2
            checkFieldNames(source, {'p0', 'p', 'p_mask', 'ux', 'uz', 'u_mask'});
        case 3
            checkFieldNames(source, {'p0', 'p', 'p_mask', 'ux', 'uy', 'uz', 'u_mask'});
    end
end

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

% check for time varying pressure source input
if isfield(source, 'p')
    
    % set source flag
    p_source = true;
    
    % if source.p_mask is not given, check if the sensor is a transducer
    if ~isfield(source, 'p_mask') && strcmp(class(sensor), 'kWaveTransducer')

        % if more than one time series is given, check the number matches
        % the number of transducer elements
        if length(source.p(:,1)) == 1
            % create indexing variable to assign source data
            p_element_index = 1;
        elseif length(source.p(:,1)) ~= sensor.num_elements
            error('The number of time series in source.p must match the number of sensor elements in the transducer');
        else
            % create indexing variable to assign source data
            p_element_index = sensor.element_voxel_index;
        end
        
        % create an indexing variable corresponding to the source elements
        ps_index = find(sensor.mask ~= 0);        
    else
        % force p_mask to be given
        enforceFields(source, {'p_mask'});
        
        % if more than one time series is given, check the number matches
        % the number of source elements
        if length(source.p(:,1)) == 1        
            % create indexing variable to assign source data
            p_element_index = 1;
        elseif length(source.p(:,1)) ~= sum(source.p_mask(:))
            error('The number of time series in source.p must match the number of source elements in source.p_mask');
        else
            % create indexing variable to assign source data
            p_element_index = 1:length(source.p(:,1));
        end
        
        % create an indexing variable corresponding to the source elements
        ps_index = find(source.p_mask ~= 0);        
    end
else
    % set source flag
    p_source = false;
end

% add check for length(ux) == length(uy) etc

% check for time varying velocity source input and set source flag
if isfield(source, 'ux') || isfield(source, 'uy') || isfield(source, 'uz') || isfield(source, 'u_mask')
    
    % if one time varying source parameter is given, enforce all of them
    switch kgrid.dim
        case 2
            enforceFields(source, {'ux', 'uz'});
        case 3
            enforceFields(source, {'ux', 'uy', 'uz'});
    end    
    
    % set source flag
    u_source = true;
    
    % if source.u_mask is not given, check if the sensor is a transducer
    if ~isfield(source, 'u_mask') && strcmp(class(sensor), 'kWaveTransducer')
        
        % if more than one time series is given, check the number matches
        % the number of transducer elements
        if length(source.ux(:,1)) == 1
            % create indexing variable to assign source data
            u_element_index = 1;
        elseif length(source.ux(:,1)) ~= sensor.num_elements
            error('The number of time series in source.ux (etc) must match the number of sensor elements in the transducer');
        else
            % create indexing variable to assign source data
            u_element_index = sensor.element_voxel_index;
        end
        
        % create an indexing variable corresponding to the source elements
        us_index = find(sensor.mask ~= 0);           
        
    else
        
        % force u_mask to be given
        enforceFields(source, {'u_mask'});
        
        % if more than one time series is given, check the number matches
        % the number of source elements
        if length(source.ux(:,1)) == 1        
            % create indexing variable to assign source data
            u_element_index = 1;
        elseif length(source.ux(:,1)) ~= sum(source.u_mask(:))
            error('The number of time series in source.ux (etc) must match the number of source elements in source.u_mask');
        else
            % create indexing variable to assign source data
            u_element_index = 1:length(source.ux(:,1));
        end
        
        % create an indexing variable corresponding to the source elements
        us_index = find(source.u_mask ~= 0);          
        
    end
else
    % set source flag
    u_source = false;
end

% =========================================================================
% CHECK KGRID STRUCTURE INPUTS
% =========================================================================

% 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('kgrid.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
    
    % 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, kgrid.dy])/max(medium.sound_speed(:)))
        disp('  WARNING: time step may be too large for a stable simulation');
    end    
end

% shorten commonly used field names
t_array = kgrid.t_array;
c = medium.sound_speed;
rho0 = medium.density;

% =========================================================================
% CHECK OPTIONAL INPUTS
% =========================================================================

% 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;
log_scale_comp_factor = LOG_SCALE_COMPRESSION_FACTOR_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_scale_log = LOG_SCALE_DEF;
plot_sim = PLOT_SIM_DEF;
plot_PML = PLOT_PML_DEF;
PML_inside = PML_INSIDE_DEF;
record_movie = RECORD_MOVIE_DEF;
return_velocity = RETURN_VELOCITY_DEF;
smooth_c = SMOOTH_C0_DEF;
smooth_p0 = SMOOTH_P0_DEF;
smooth_rho0 = SMOOTH_RHO0_DEF;
use_kspace = USE_KSPACE_DEF;
use_sg = USE_SG_DEF;

% assign the default input parameters that vary for different dimensions
switch kgrid.dim
    case 1
        PML_x_alpha = PML_ALPHA_DEF;
        PML_x_size = PML_SIZE_DEF;        
    case 2
        PML_x_alpha = PML_ALPHA_DEF;
        PML_z_alpha = PML_ALPHA_DEF;
        PML_x_size = PML_SIZE_DEF;
        PML_z_size = PML_SIZE_DEF;
        mesh_plot = MESH_PLOT_DEF;
        movie_type = MOVIE_TYPE_DEF;
    case 3
        PML_x_alpha = PML_ALPHA_DEF;
        PML_y_alpha = PML_ALPHA_DEF;
        PML_z_alpha = PML_ALPHA_DEF;
        PML_x_size = PML_SIZE_DEF;
        PML_y_size = PML_SIZE_DEF;
        PML_z_size = PML_SIZE_DEF;
end

% 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 ~(strcmp(display_mask, 'off') || all(size(display_mask) == size(kgrid.k)))
                    error('Optional input ''DisplayMask'' must be the same size as the k-space grid or set to ''off''');
                end
            case 'LogScale'
                plot_scale_log = varargin{input_index + 1};
                if numel(plot_freq) == 1 && isnumeric(plot_scale_log) && plot_scale_log > 0
                    log_scale_comp_factor = plot_scale_log;
                    plot_scale_log = true;
                elseif ~islogical(plot_scale_log)
                    error('Optional input ''LogScale'' must be Boolean or a single numerical value > 0');
                end              
            case 'MeshPlot'
                if kgrid.dim == 2
                    mesh_plot = varargin{input_index + 1};
                    if ~islogical(mesh_plot)
                        error('Optional input ''MeshPlot'' must be Boolean');
                    end
                else
                    error('Optional input ''MeshPlot'' only supported in 2D');
                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 'MovieType'
                if kgrid.dim == 2
                    movie_type = varargin{input_index + 1};
                    if ~(strcmp(movie_type, 'frame') || strcmp(movie_type, 'image'))
                        error('Optional input ''MovieType'' must be set to ''frame'' or ''image''');
                    end
                else
                    error('Optional input ''MovieType'' only supported in 2D');
                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 'PlotPML'
                plot_PML = varargin{input_index + 1};
                if ~islogical(plot_PML)
                    error('Optional input ''PlotPML'' must be Boolean');
                end                 
            case 'PlotScale'
                plot_scale = varargin{input_index + 1};
                if ~strcmp(plot_scale, 'auto') && (~(numel(plot_scale) == 2 && isnumeric(plot_scale)))
                    error('Optional input ''PlotScale'' must be a 2 element numerical array or set to ''auto''');
                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}) > kgrid.dim
                    if kgrid.dim > 1
                        error(['Optional input ''PMLAlpha'' must be a 1 or ' kgrid.dim ' element numerical array']);
                    else
                        error('Optional input ''PMLAlpha'' must be a single numerical value');
                    end                    
                end
                switch kgrid.dim
                    case 1
                        PML_x_alpha = varargin{input_index + 1}(1); 
                    case 2
                        PML_x_alpha = varargin{input_index + 1}(1);
                        PML_z_alpha = varargin{input_index + 1}(end);                        
                    case 3
                        PML_x_alpha = varargin{input_index + 1}(1);
                        PML_y_alpha = varargin{input_index + 1}(ceil((end + 1)/2));
                        PML_z_alpha = varargin{input_index + 1}(end);
                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}) > kgrid.dim
                    if kgrid.dim > 1
                        error(['Optional input ''PMLSize'' must be a 1 or ' kgrid.dim ' element numerical array']);
                    else
                        error('Optional input ''PMLSize'' must be a single numerical value');
                    end
                end
                switch kgrid.dim
                    case 1
                        PML_x_size = varargin{input_index + 1}(1);
                    case 2
                        PML_x_size = varargin{input_index + 1}(1);
                        PML_z_size = varargin{input_index + 1}(end);
                    case 3
                        PML_x_size = varargin{input_index + 1}(1);
                        PML_y_size = varargin{input_index + 1}(ceil((end + 1)/2));
                        PML_z_size = varargin{input_index + 1}(end);
                end
            case 'RecordMovie'
                record_movie = varargin{input_index + 1};    
                if ~islogical(record_movie)
                    error('Optional input ''RecordMovie'' must be Boolean');
                end
            case 'ReturnVelocity'
                return_velocity = varargin{input_index + 1};
                if ~islogical(return_velocity)
                    error('Optional input ''ReturnVelocity'' 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_rho0 = 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');
            case 'UsekSpace'
                use_kspace = varargin{input_index + 1}; 
                if ~islogical(use_kspace)
                    error('Optional input ''UsekSpace'' must be Boolean');
                end
            case 'UseSG'
                use_sg = varargin{input_index + 1}; 
                if ~islogical(use_sg)
                    error('Optional input ''UseSG'' must be Boolean');
                end                   
            otherwise
                error(['Unknown optional input ' varargin{input_index}]);
        end
    end
end

% enforce density input if velocity sources or output are being used
if ~user_medium_density_input && (u_source || return_velocity)
    error('medium.density must be explicitly defined if velocity inputs or outputs are used, even in homogeneous media');
end

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

% check for automatic plot scaling
if strcmp(plot_scale, 'auto') 
    plot_scale_auto = true;
else
    plot_scale_auto = false;
end

% check for log plot scaling
if plot_scale_log && ~plot_scale_auto
    % store the linear and log plot scales
    alt_plot_scale_lin = plot_scale;
    alt_plot_scale_log = log10(abs(plot_scale) + log_scale_comp_factor) - log10(log_scale_comp_factor);
    alt_plot_scale_log(1) = -alt_plot_scale_log(1);
end

% 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

% ensure default display mask is switched off if sensor input is empty
if ~use_sensor && strcmp(display_mask, 'default')
    display_mask = 'off';
end

% switch off default display mask if using mesh plot
if kgrid.dim == 2 && mesh_plot && strcmp(display_mask, 'default')
    display_mask = 'off';
end

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

% enforce GPUmat compatability by using wrapper for GPUsingle and GPUdouble
if strcmp(data_cast, 'GPUsingle');
    data_cast = 'kWaveGPUsingle';
elseif strcmp(data_cast, 'GPUdouble');
    data_cast = 'kWaveGPUdouble';
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 p0 is given
if isfield(source, 'p0') && ~time_rev && plot_sim && ~plot_scale_auto
    
    % find the maximum input pressure amplitude
    if isfield(source, 'p')
        max_val = max([source.p0(:); source.p(:)]);
    else
        max_val = max(source.p0(:));
    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

% cleanup unused variables
clear *_DEF NUM_REQ_INPUT_VARIABLES user_medium_density_input;

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

if use_sensor

    % 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 (kgrid.dim == 3 && numDim(sensor.mask) == 3) || (kgrid.dim ~= 3 && 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
        switch kgrid.dim
            case 1
                sensor_x = sensor.mask;
            case 2
                sensor_x = sensor.mask(1, :);
                sensor_z = sensor.mask(2, :);
            case 3
                sensor_x = sensor.mask(1, :);
                sensor_y = sensor.mask(2, :);
                sensor_z = sensor.mask(3, :);    
        end

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

            % check if directivity is being used
            if kgrid.dim == 2 && compute_directivity
                error('sensor directivity fields are only compatible with binary sensor masks or ''CartInterp'' set to ''nearest''');
            end        
        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(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
end