% DESCRIPTION:
%       subscript to create storage variables
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 1st August 2011
%       last update - 7th November 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/>.

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

if use_sensor

    % check sensor mask based on the Cartesian interpolation setting
    if ~binary_sensor_mask && strcmp(cartesian_interp, 'nearest')
        
        % extract the data using the binary sensor mask created in
        % inputChecking, but switch on Cartesian reorder flag so that the
        % final data is returned in the correct order (not in time
        % reversal mode).
        binary_sensor_mask = true;
        if ~time_rev
            reorder_data = true;
        end

        % check if any duplicate points have been discarded in the
        % conversion from a Cartesian to binary mask
        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        

    end
    
    % create mask indices (this works for both normal sensor and transducer
    % inputs)
    sensor_mask_ind = find(sensor.mask ~= 0);
    
    % convert the data type depending on the number of indices (this saves
    % memory)
    eval(['sensor_mask_ind = ' index_data_type '(sensor_mask_ind);']); 

    % create storage and scaling variables
    if ~time_rev
        
        % preallocate storage variables, if return_velocity is true, the
        % outputs are assigned as structure fields, otherwise the
        % appropriate pressure data is assigned directly to sensor_data
        if return_velocity
            
            % pre-allocate the sensor_data variable based on the number of
            % binary or Cartesian sensor points
            if binary_sensor_mask
                sensor_data.p = zeros(sum(sensor.mask(:)), length(t_array));
            else
                sensor_data.p = zeros(length(sensor_x), length(t_array));
            end
            
            % pre-allocate the velocity fields based on the number of
            % dimensions in the simulation
            switch kgrid.dim
                case 1
                    sensor_data.ux = sensor_data.p;
                case 2
                    sensor_data.ux = sensor_data.p;
                    sensor_data.uy = sensor_data.p;
                case 3
                    sensor_data.ux = sensor_data.p;
                    sensor_data.uy = sensor_data.p;
                    sensor_data.uz = sensor_data.p;  
            end
        
        % return velocity is false    
        else
            
            % if streaming data to disk, reduce to the size of the
            % sensor_data matrix based on the value of stream_to_disk
            if kgrid.dim == 3 && stream_to_disk
                num_time_points = stream_to_disk;
                
                % initialise the file index variable
                stream_data_index = 1;
            else
                num_time_points = length(t_array);
            end
            
            % binary sensor mask
            if binary_sensor_mask
                if transducer_sensor
                    if transducer_receive_elevation_focus
                        % if there is elevation focusing, a buffer is
                        % needed to store a short time history at each
                        % sensor point before averaging
                        sensor_data_buffer_size = max(sensor.elevation_beamforming_delays) + 1;
                        if sensor_data_buffer_size > 1
                            sensor_data_buffer = zeros(sum(sensor.mask(:)), sensor_data_buffer_size);
                        else
                            clear sensor_data_buffer sensor_data_buffer_size;
                            transducer_receive_elevation_focus = false;
                        end
                   end
                    
                    % the grid points can be summed on the fly and so the
                    % sensor is the size of the number of active elements 
                    sensor_data = zeros(sensor.number_active_elements, num_time_points);
                    
                else
                    if store_time_series
                        % store time series
                        sensor_data = zeros(sum(sensor.mask(:)), num_time_points);
                    else
                        % store only statistics at each sensor point
                        sensor_data.p_rms = zeros(sum(sensor.mask(:)), 1);
                        sensor_data.p_max = zeros(sum(sensor.mask(:)), 1);
                    end
                end
                
            % Cartesian sensor mask    
            else
                if store_time_series
                    % store time series
                    sensor_data = zeros(length(sensor_x), num_time_points);
                else
                    % store only statistics at each sensor point
                    sensor_data.p_rms = zeros(length(sensor_x), 1);
                    sensor_data.p_max = zeros(length(sensor_x), 1);
                end
                
            end
            
            clear num_time_points;
        end       
    end
end

% =========================================================================
% PRECOMPUTE DATA STORAGE CASE
% =========================================================================

if time_rev || ~use_sensor
    % do not store any data in time reversal mode
    extract_data_case = 0;
elseif ~(kgrid.dim == 2 && compute_directivity)
    if ~return_velocity && ~binary_sensor_mask
        % return velocity = false
        % binary sensor mask = false
        extract_data_case = 1;
        
    elseif ~return_velocity && binary_sensor_mask
        % return velocity = false
        % binary sensor mask = true      
        extract_data_case = 2;

    elseif return_velocity && ~binary_sensor_mask 
        % return velocity = true
        % binary sensor mask = false     
        extract_data_case = 3;

    elseif return_velocity && binary_sensor_mask 
        % return velocity = true
        % binary sensor mask = true    
        extract_data_case = 4;

    else
        error('Unknown data output combination...');
    end
elseif ~return_velocity && binary_sensor_mask 
    % compute directivity = true (only supported in 2D)
    % return velocity = false (must be false if directivity = true)
    % binary sensor mask = true (must be true if directivity = true)
    extract_data_case = 5;      
else
    error('Unknown data output combination...');    
end