%----------------------------------------------------------
% input:
%   theta   incidence angle (degrees)
%   model:  array containing P and S velocities
%           denisties, and Q values 
%   dt:     time between samples
%   NT:     number of time sample
% ouput:
%   time:   array of time values
%   RF_filt receiver function
%   
function [time, RF_filt] = RecFun(theta, model, dt, NT)
% Computes receiver function for a flat layered structure  
% Calls PSvResponse
%
% time parameters
%

    tminplot = -5;
    tmaxplot = 20;
%     tmaxplot = 10;

    tmin = 0;
    tmax = tmin + (NT - 1) * dt;

% frequencies
    [fim, fmin, fmax, df, NF, freq] = get_freqs(tmaxplot, dt, NT);
 
% filter parameters
    npoles = 3;
    f0 = 0.25;

% set type of wave
    wavetype = [theta,'P'];

% calculate response
    response = PSvResponse(wavetype, freq, model);
    
% extract response
    Vp = response(1,:);
    Hp = response(2,:);
    
% filter
    [time, RF_filt] = do_filter(Hp, Vp, f0, fim, freq, npoles, NF, tminplot, tmaxplot, tmax, dt);

% normalize
    RF_filt = RF_filt ./ max(RF_filt);
 
%----------------------------------------------------------
% calculates offset and delay time
% given P and S wave velocities
% input:
%   theta     incidence angle (degrees)
%   a1        P wave velocity (top layer)
%   b1        S wave velocity (top layer)
% returns:
%   offset    distance
%   delay     time difference between P and S (converted)
%----------------------------------------------------------
function [offset, delay] = get_offset(theta, a1, b1, h1)
    offset = h1 * sin(pi * theta/180.0);  
    path = sqrt(offset^2 + h1^2);
    delay = path/b1 - path/a1;
% end function
%----------------------------------------------------------
% calculates frequencies
% input:
%   tmaxplot  length of time window 
%   dt        time between samples
%   NT        number of time samples
% returns:
%   fmin    minimum frequency
%   fmax    maximum frequency
%   df      spacing between frequencies
%   NF      number of frequencies
%   freq    complex array of frequencies
%----------------------------------------------------------
function [fim, fmin, fmax, df, NF, freq] = get_freqs(tmaxplot, dt, NT)
% set minimum
    fim = 0.5/tmaxplot;  
    fmin = 0;
    fmax = 0.5/dt;
    df = 1/(dt*NT);
    
% number of frequencies

    NF = floor((fmax-fmin)/df+1.5);
    freq = fmin:df:(NF-1)*df;
    freq = freq+i*fim*ones(size(freq));   

% end function
%----------------------------------------------------------
% calculates frequencies
% input:
%   HP       
%   Vp          P wave velocity
%   f0          minimum frequency
%   freq        frequency bins
%   npoles      number of poles in filter
%   NF          number of frequencies
%   tminplot    minimum time
%   tmax        maximum time
%   dt          time spacing
% returns:
%   time       time axis
%   RF_filt    filtered receiver function
%----------------------------------------------------------
function [time, RF_filt] = do_filter(Hp, Vp, f0, fim, freq, npoles, NF, tminplot, tmaxplot, tmax, dt)
    
%    Filter = ones(size(freq))./ ...
%    sqrt(ones(size(freq))+(freq./(ones(size(freq)).*f0)).^(2*npoles*ones(size(freq))));
    alpha=1.0;
    Filter=exp(-(pi.*freq./alpha).^2); % alpha is Gaussian filter
    
    Shift = exp(-i*2*pi*freq*tminplot);
    ratio = (Hp./Vp).*Filter.*Shift;
    decon = (1./dt)*ifft([conj(ratio) ratio(NF-1:-1:2)]);
    
    time = (tminplot:dt:min(tmaxplot,tmax));
    RF_filt = exp(2.*pi*fim*(time-tminplot)).*real(decon(1:length(time)));

    % end function
%----------------------------------------------------------