c subroutine RTCOEF calculates reflection/transmission coefficients
c for interface between two solid layers, based on the equations on 
c p. 150 of Aki and Richards.
c
c  Inputs:    vp1     =  P-wave velocity of layer 1 (top layer)
c  (real)     vs1     =  S-wave velocity of layer 1
c             den1    =  density of layer 1
c             vp2     =  P-wave velocity of layer 2 (bottom layer)
c             vs2     =  S-wave velocity of layer 2
c             den2    =  density of layer 2
c             hslow   =  horizontal slowness (ray parameter)
c  Returns:   rt(1)   =  down P to P up     (refl)
c  (complex)  rt(2)   =  down P to S up     (refl)
c             rt(3)   =  down P to P down   (tran)
c             rt(4)   =  down P to S down   (tran)
c             rt(5)   =  down S to P up     (refl)
c             rt(6)   =  down S to S up     (refl)
c             rt(7)   =  down S to P down   (tran)
c             rt(8)   =  down S to S down   (tran)
c             rt(9)   =    up P to P up     (tran)
c             rt(10)  =    up P to S up     (tran)
c             rt(11)  =    up P to P down   (refl)
c             rt(12)  =    up P to S down   (refl)
c             rt(13)  =    up S to P up     (tran)
c             rt(14)  =    up S to S up     (tran)
c             rt(15)  =    up S to P down   (refl)
c             rt(16)  =    up S to S down   (refl)
c
c NOTE:  All input variables are real.  
c        All output variables are complex!
c        Coefficients are not energy normalized.
c
      SUBROUTINE RTCOEF(vp1,vs1,den1,vp2,vs2,den2,hslow,rt)
      implicit complex (a-h,o-z)
      complex rt(16)
      real vp1,vs1,den1,vp2,vs2,den2,hslow
c
      alpha1=cmplx(vp1,0.)
      beta1=cmplx(vs1,0.)
      rho1=cmplx(den1,0.)
      alpha2=cmplx(vp2,0.)
      beta2=cmplx(vs2,0.)
      rho2=cmplx(den2,0.)
      p=cmplx(hslow,0.)
c
      cone=cmplx(1.,0.)
      ctwo=cmplx(2.,0.)
c
      si1=alpha1*p
      si2=alpha2*p          
      sj1=beta1*p
      sj2=beta2*p          
c
      ci1=csqrt(cone-si1**2)
      ci2=csqrt(cone-si2**2)
      cj1=csqrt(cone-sj1**2)
      cj2=csqrt(cone-sj2**2)         
c
      term1=(cone-ctwo*beta2*beta2*p*p)
      term2=(cone-ctwo*beta1*beta1*p*p)
      a=rho2*term1-rho1*term2
      b=rho2*term1+ctwo*rho1*beta1*beta1*p*p
      c=rho1*term2+ctwo*rho2*beta2*beta2*p*p
      d=ctwo*(rho2*beta2*beta2-rho1*beta1*beta1)
      E=b*ci1/alpha1+c*ci2/alpha2
      F=b*cj1/beta1+c*cj2/beta2
      G=a-d*ci1*cj2/(alpha1*beta2)
      H=a-d*ci2*cj1/(alpha2*beta1)
      DEN=E*F+G*H*p*p
c
c
      trm1=b*ci1/alpha1-c*ci2/alpha2          
      trm2=a+d*ci1*cj2/(alpha1*beta2)
      rt(1)=(trm1*F-trm2*H*p*p)/DEN              !refl down P to P up
c
      trm1=a*b+c*d*ci2*cj2/(alpha2*beta2)       
      rt(2)=(-ctwo*ci1*trm1*p)/(beta1*DEN)       !refl down P to S up
c 
      rt(3)=ctwo*rho1*ci1*F/(alpha2*DEN)         !trans down P to P down
c
      rt(4)=ctwo*rho1*ci1*H*p/(beta2*DEN)        !trans down P to S down
c
      trm1=a*b+c*d*ci2*cj2/(alpha2*beta2)       
      rt(5)=(-ctwo*cj1*trm1*p)/(alpha1*DEN)      !refl down S to P up
c
      trm1=b*cj1/beta1-c*cj2/beta2               
      trm2=a+d*ci2*cj1/(alpha2*beta1)
      rt(6)=-(trm1*E-trm2*G*p*p)/DEN             !refl down S to S up
c
      rt(7)=-ctwo*rho1*cj1*G*p/(alpha2*DEN)      !trans down S to P down 
c
      rt(8)=ctwo*rho1*cj1*E/(beta2*DEN)          !trans down S to S down
c
c
      trm1=b*ci1/alpha1-c*ci2/alpha2          
      trm2=a+d*ci2*cj1/(alpha2*beta1)
      rt(11)=-(trm1*F+trm2*G*p*p)/DEN            !refl up P to P down
c
      trm1=a*c+b*d*ci1*cj1/(alpha1*beta1)       
      rt(12)=(ctwo*ci2*trm1*p)/(beta2*DEN)       !refl up P to S down
c
      rt(9)=ctwo*rho2*ci2*F/(alpha1*DEN)         !trans up P to P up
c
      rt(10)=-ctwo*rho2*ci2*G*p/(beta1*DEN)      !trans up P to S up
c
      trm1=a*c+b*d*ci1*cj1/(alpha1*beta1)       
      rt(15)=(ctwo*cj2*trm1*p)/(alpha2*DEN)      !refl up S to P down
c
      trm1=b*cj1/beta1-c*cj2/beta2               
      trm2=a+d*ci1*cj2/(alpha1*beta2)
      rt(16)=(trm1*E+trm2*H*p*p)/DEN             !refl up S to S down
c                                               
      rt(13)=ctwo*rho2*cj2*H*p/(alpha1*DEN)      !trans up S to P up
c
      rt(14)=ctwo*rho2*cj2*E/(beta1*DEN)         !trans up S to S up
c
c
      return
      end