% DESCRIPTION:
%       subscript to create absorption and storage variables
%
% ABOUT:
%       author      - Bradley Treeby
%       date        - 26th November 2010
%       last update - 26th November 2010
%       
% 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/>.

% =========================================================================
% CREATE ABSORPTION VARIABLES
% =========================================================================

% define the lossy derivative operators and proportionality coefficients
if strcmp(equation_of_state, 'absorbing')
            
    % make sure the operators are positive and real
    medium.alpha_coeff = abs(real(medium.alpha_coeff));
    medium.alpha_power = abs(real(medium.alpha_power));
    
    % 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
    if ~(isfield(medium, 'alpha_mode') && strcmp(medium.alpha_mode, 'no_absorption'))
        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;
    else
        absorb_param = 0;
    end
       
    % compute the dispersive fractional Laplacian operator and coefficient
    if ~(isfield(medium, 'alpha_mode') && strcmp(medium.alpha_mode, 'no_dispersion'))
        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;
    else
        dispers_param = 0;
    end
        
    % filter the absorption parameters
    if isfield(medium, 'alpha_filter');
                
        % update command line status
        disp('  filtering absorption variables...');        
        
        % frequency shift the absorption parameters
        absorb_param = fftshift(absorb_param);
        dispers_param = fftshift(dispers_param);
                        
        % apply the filter
        absorb_param = absorb_param.*medium.alpha_filter;
        dispers_param = dispers_param.*medium.alpha_filter;

        % shift the parameters back
        absorb_param = ifftshift(absorb_param);
        dispers_param = ifftshift(dispers_param); 
           
    end    

    % modify the sign of the absorption operators
    if isfield(medium, 'alpha_sign')
       if numel(medium.alpha_sign) == 2
           % if two parameters are given, apply the first to the absorption
           % parameter and the second to the disperion parameters
           absorb_param = sign(medium.alpha_sign(1))*absorb_param;
           dispers_param = sign(medium.alpha_sign(2))*dispers_param;
       else
           error('medium.alpha_sign must be given as a 2 element array controlling absorption and dispersion, respectively.');
       end
    end
        
    % cleanup unused variables
    clear absorb_nabla* absorb_tau absorb_eta;
    
elseif strcmp(equation_of_state, 'stokes')
    
    % make sure the operators are positive
    medium.alpha_coeff = abs(medium.alpha_coeff);
    medium.alpha_power = abs(medium.alpha_power);
    
    % 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
% =========================================================================

if use_sensor

    % 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

        % time reversal is off    
        case 0

            % 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
                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
                switch kgrid.dim
                    case 1
                        sensor_data.ux = sensor_data.p;
                    case 2
                        sensor_data.ux = sensor_data.p;
                        sensor_data.uz = sensor_data.p;
                    case 3
                        sensor_data.ux = sensor_data.p;
                        sensor_data.uy = sensor_data.p;
                        sensor_data.uz = sensor_data.p;  
                end
            else
                if binary_sensor_mask
                    sensor_data = zeros(sum(sensor.mask(:)), length(t_array));
                else
                    sensor_data = zeros(length(sensor_x), length(t_array));
                end
            end                      

        % thresholded time reversal    
        case 2

            % preallocate the sensor variable
            if kgrid.dim == 1
                p_sensor = zeros(kgrid.Nx, 1);
            else
                p_sensor = zeros(1, numel(kgrid.k));
            end

            if kgrid.dim == 2
                % preallocate the threshold variable used for visualisation
                threshold_mask = zeros(kgrid.Nz, kgrid.Nx);
            end

            % 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 threshold to zero
            sensor.time_reversal_boundary_data(abs(sensor.time_reversal_boundary_data) < sensor.time_reversal_adapt_thresh) = 0;        

    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