function [v, u, zz, zx] = raydep(p,om,d,ro,vp2,vs2)
%RAYDEP	Propagator matrix sol'n for P-SV waves, minor vector method
%   u   =  horizontal velocity AT TOP OF EACH LAYER
%   v   =  vertical velocity AT TOP OF EACH LAYER
%   sigzx = horizontal stress AT TOP OF EACH LAYER
%   sigzz = vertical stress AT TOP OF EACH LAYER
%   {  Normalized compliance = -k*v/(omega*sigzz)  }
% [v u sigzz sigzx] = raydep(p,omega,d,ro,vp2,vs2)
%   p     = slowness (s/m) of surface wave
%   omega = angular frequency (radians/sec) of surface wave
%   d     = thicknesses of the model layers (meters?)
%   rho   = density of the layer (kg/m^3) (= gm/cc * 1000)
%   vp2   = compressional velocity squared (m/s)^2
%   vs2   = shear velocity squared (m/s)^2

% W.C. Crawford 1-5-89
%   x = stress-displacement vector:
%		( vertical velocity,  horiz velocity, vert stress, horiz stress )
%	y = minor vector matrix
mu = ro .* vs2;
n = length(d);
ist = n;
ysav = 0;
psq = p*p;
r2 = 2*mu(ist)*p;
% R and S are the "Wavenumbers" of compress and shear waves in botlayer
% RoW and SoW are divided by ang freq
RoW = sqrt(psq- 1/vp2(ist));
SoW = sqrt(psq- 1/vs2(ist));
ym = zeros(ist,5);
i = ist;
y(3) =  RoW;
y(4) = -SoW;
y(1) = (RoW*SoW - psq)/ro(i);
y(2) = r2*y(1) + p;
y(5) = ro(i) - r2*(p + y(2));
ym(i,:) = y;
%*****PROPAGATE UP LAYERS*********
while i > 1
  i = i-1;
  ha = psq-1/vp2(i);
  [ca sa] = argdtray(om*d(i),ha);
  hb = psq - 1/vs2(i);
  [cb sb] = argdtray(om*d(i),hb);
  hbs = hb*sb;
  has = ha*sa;
  r1 = 1/ro(i);
  r2 = 2*mu(i)*p;
  b1 = r2*y(1) - y(2);
  g3 = ( y(5) + r2*(y(2)-b1) ) * r1;
  g1 = b1 + p*g3;
  g2 = ro(i)*y(1) - p*(g1+b1);
  e1 = cb*g2 - hbs*y(3);
  e2 = -sb*g2 + cb*y(3);
  e3 = cb*y(4) + hbs*g3;
  e4 = sb*y(4) + cb*g3;
  y(3) = ca*e2 - has*e4;
  y(4) = sa*e1 + ca*e3;
  g3 = ca*e4 - sa*e2;
  b1 = g1 - p*g3;
  y(1) = (ca*e1 + has*e3 + p*(g1+b1))*r1;
  y(2) = r2*y(1) - b1;
  y(5) = ro(i)*g3 - r2*(y(2) - b1);
  ym(i,:) = y;
end
de = y(5)/sqrt(y(1)*y(1) + y(2)*y(2));
ynorm = 1/y(3);
y(1:4) = [0  -ynorm  0  0];
%*****PROPAGATE BACK DOWN LAYERS*********
while i <= ist
  x(i,1) = -ym(i,2)*y(1) - ym(i,3)*y(2) + ym(i,1)*y(4);
  x(i,2) = -ym(i,4)*y(1) + ym(i,2)*y(2) - ym(i,1)*y(3);
  x(i,3) = -ym(i,5)*y(2) - ym(i,2)*y(3) - ym(i,4)*y(4);
  x(i,4) =  ym(i,5)*y(1) - ym(i,3)*y(3) + ym(i,2)*y(4);
  ls = i;
  if i >=2
    sum = abs( x(i,1) + i*x(i,2));
    pbsq = 1/vs2(i);
    if sum < 1e-4
        break
    end
  end
  ha = psq - 1/vp2(i);
  [ca sa] = argdtray(om*d(i),ha);
  hb = psq-1/vs2(i);
  [cb sb] = argdtray(om*d(i),hb);
  hbs = hb*sb;
  has = ha*sa;
  r2 = 2*p*mu(i);
  e2 = r2*y(2) - y(3);
  e3 = ro(i)*y(2) - p*e2;
  e4 = r2*y(1) - y(4);
  e1 = ro(i)*y(1) - p*e4;
  e6 = ca*e2 - sa*e1;
  e8 = cb*e4 - sb*e3;
  y(1) = (ca*e1 - has*e2+p*e8) / ro(i);
  y(2) = (cb*e3 - hbs*e4+p*e6) / ro(i);
  y(3) = r2*y(2) - e6;
  y(4) = r2*y(1) - e8;
  i = i+1;
end
if x(1,3) == 0
  error('vertical surface stress = 0 in DETRAY');
end
ist = ls;
v =	 x(:,1);
u =  x(:,2);
zz = x(:,3);
zx = x(:,4);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [c, s] = argdtray(wd,h)
%ARGDTRAY	Used by raydep.m to calculate surface wave vectors
%	[c s] = argdtray(wd, h)
%   wd = thickenss of layer times ang frequency
%   h = squared wavenumber/freq of shear or comp wave
%   c = 
%   s = 

% W.C. Crawford 5-5-93

hh = sqrt(abs(h)); % magnitude of wavenumber/freq
th=wd*hh;			% number of waves (or e-foldings) in layer in radians
if th >= 1.5e-14
    if h <= 0   % propagating wave
        c =  cos(th);
        s = -sin(th)/hh;
    else 			% evenescent wave
		if (th >= 88)
        	disp(sprintf('TH = %f, the freqs are possibly too high',th));
		end
    	d=exp(th);
    	c =  0.5*(d + 1/d);
    	s = -0.5*(d - 1/d)/hh;
	end
else
    c = 1;
    s = -wd;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
