%
% PSvResponse
%
   function y = PSvResponse(wavetype, frequencies, model)
%
% Frequency Domain Solution to the Surface
% Response of a Layered Halfspace to
% P wave incident from Below
%
%  wavetype = [incidence angle in degrees, 'P' or 'S']
%
%  frequencies = vector of frequencies 
%               
%  model = [h;a;b;rho;qp;qs]
%
%     where
%
%       h = vector of layer thicknesses (terminated by zero)
%       a = vector of P velocities
%       b = vector of S velocities
%       rho = vector of densities
%       qp = vector of Q-sub-P's
%       qs = vector of Q-sub-S's
%
% unpack the model
%
    h = model(1,:);
    a = model(2,:);
    b = model(3,:);
    rho = model(4,:);
    qp = model(5,:);
    qs = model(6,:);
    N = length(a);
%
% unpack wavetype
%
    if wavetype(2)=='P'
         p=sin(wavetype(1)*pi/180)/a(N);
    elseif wavetype(2)=='S'
        p=sin(wavetype(1)*pi/180)/b(N);
    else
        error('wavetype(2) must be P or S')
    end
%
% Set up storage 
%
    Hp = zeros(size(frequencies));
    Hs = zeros(size(frequencies));
    Vp = zeros(size(frequencies));
    Vs = zeros(size(frequencies));
%
% Compute layer quantities
%
    gs = ones(1,N)-.5*i*ones(1,N)./qs;
    gp = ones(1,N)-.5*i*ones(1,N)./qp;

    etas = ((ones(1,N)./(b(1:N).*gs(1:N)).^2)-p^2*ones(1,N)).^.5;
    etap = ((ones(1,N)./(a(1:N).*gp(1:N)).^2)-p^2*ones(1,N)).^.5;
    epss = (2*rho.*etas).^(-.5);
    epsp = (2*rho.*etap ).^(-.5);

    Mup = [-i*etap(1)*epsp(1), p*epss(1); p*epsp(1), -i*etas(1)*epss(1)];
    Mdp = [i*etap(1)*epsp(1), p*epss(1); p*epsp(1), i*etas(1)*epss(1)];
    
    Nup = rho(1)*[(2*b(1)^2*p^2-1)*epsp(1),-2*i*b(1)^2*p*etas(1)*epss(1); ...
    -2*i*b(1)^2*p*etap(1)*epsp(1),(2*b(1)^2*p^2-1)*epss(1)];

    Ndp = rho(1)*[(2*b(1)^2*p^2-1)*epsp(1),2*i*b(1)^2*p*etas(1)*epss(1); ...
    2*i*b(1)^2*p*etap(1)*epsp(1),(2*b(1)^2*p^2-1)*epss(1)];

    Rfree = -inv(Ndp)*Nup;
    Mufree = Mup;
    Mdfree = Mdp;

for k = 1: N-1

    loc = 2*k-1:2*k;
    Mum = Mup;
    Mdm = Mdp;
    Num = Nup;
    Ndm = Ndp;
  
    Mup = [-i*etap(k+1)*epsp(k+1), p*epss(k+1); p*epsp(k+1), -i*etas(k+1)*epss(k+1)];
    Mdp = [i*etap(k+1)*epsp(k+1), p*epss(k+1); p*epsp(k+1), i*etas(k+1)*epss(k+1)];

    Nup = rho(k+1)*[(2*b(k+1)^2*p^2-1)*epsp(k+1), ...
    -2*i*b(k+1)^2*p*etas(k+1)*epss(k+1); ...
    -2*i*b(k+1)^2*p*etap(k+1)*epsp(k+1),(2*b(k+1)^2*p^2-1)*epss(k+1)];

    Ndp = rho(k+1)*[(2*b(k+1)^2*p^2-1)*epsp(k+1), ...
     2*i*b(k+1)^2*p*etas(k+1)*epss(k+1); ...
     2*i*b(k+1)^2*p*etap(k+1)*epsp(k+1),(2*b(k+1)^2*p^2-1)*epss(k+1)];

    denom = inv(conj(Mum')*Ndp-conj(Num')*Mdp);

    Rd(loc,1:2) = -(conj(Mdm')*Ndp-conj(Ndm')*Mdp)*denom;
    Ru(loc,1:2) = -denom*(conj(Mum')*Nup-conj(Num')*Mup);
    Td(loc,1:2) = i*denom; 
    Tu(loc,1:2) = conj(Td(loc,1:2)');
  
end

    Rhalf = -inv(Ndp)*Nup;
    Muhalf = Mup;
    Mdhalf = Mdp;
    nf = 0;
    matid = eye(2,2);
%
% Do the recursion (following Kennett, 1983)
%
    for f = frequencies
        nf = nf+1;
        Rdtot = zeros(2,2);
        Tutot = matid;
  
        for k = N-1:-1:1
            loc = 2*k-1:2*k;
            Reverb = inv(matid-Rdtot*Ru(loc,1:2));
            
            Tutot = Tu(loc,1:2)*Reverb*Tutot;   
            Rdtot = Rd(loc,1:2)+Tu(loc,1:2)*Reverb*Rdtot*Td(loc,1:2);
            
            Ed = diag([exp(i*2*pi*f*etap(k)*h(k)),exp(i*2*pi*f*etas(k)*h(k))]);
            Rdtot = Ed*Rdtot*Ed;
            Tutot = Ed*Tutot;
        end
%
        surface = (Mufree+Mdfree*Rfree)*inv(matid-Rdtot*Rfree)*Tutot;
        basement = Muhalf+Mdhalf*Rhalf;
        
        Vp(nf) = i*(2*rho(N)/etap(N))^.5*surface(1,1);
        Hs(nf) = i*(2*rho(N)/etas(N))^.5*surface(2,2);
        Vs(nf) = (2*rho(N)/etas(N))^.5*surface(1,2);
        Hp(nf) = (2*rho(N)/etap(N))^.5*surface(2,1);
    end
    
y = [Vp; Hp; Vs; Hs];

