C 
C Subroutine file 'coef.for' - reflection/transmission coefficients
C
C Date: 1989, December 8
C Coded by: Ludek Klimes
C
C.......................................................................
C
C This file consists of:
C     COEFSH..Subroutine designed to evaluate SH R/T coefficients, see
C             C.R.T.5.9.7.
C             COEFSH
C     COEF50..Subroutine designed to evaluate P-SV R/T coefficients, see
C             C.R.T.5.9.7.
C             COEF50
C
C
C=======================================================================
C
C     
C
      SUBROUTINE COEFSH(P,VS1,RO1,VS2,RO2,NCODE,RMOD,RPH)
C
C     ******************************************************************
C
C     The routine COEFSH is designed for the computation of reflection
C     and transmission coefficients at a plane interface between two
C     homogeneous solid halfspaces or at a free surface of a homogeneous
C     solid halfspace.
C
C     I n p u t   p a r a m e t e r s
C           P...Ray parameter
C           VS1,RO1...Parameters of the first halfspace
C           VS2,RO2...Parameters of second halfspace. For the free
C                    surface take RO2=0. And arbitrary
C                    value of VS2
C           NCODE...Code of the computed coefficient
C                    S1S1...NCODE=1
C                    S1S2...NCODE=2
C
C     O u t p u t   p a r a m e t e r s
C           RMOD,RPH...Modul and argument of the coefficient
C
C     N o t e s
C     1/ Time factor of incident wave...EXP(-i*OMEGA*T)
C     2/ Formulae are taken from Cerveny,Molotkov and Psencik,
C        Ray Method in Seismology, page 34
C
C     ******************************************************************
C
      COMPLEX A,B
      X= 1.-P*P*VS1*VS1
      Y= 1.-P*P*VS2*VS2
      C= RO1*VS1*SQRT(ABS(X))
      D= RO2*VS2*SQRT(ABS(Y))
      A= CMPLX(C,0.)
      IF(X.LT.0.) A= CMPLX(0.,C)
      B= CMPLX(D,0.)
      IF(Y.LT.0.) B= CMPLX(0.,D)
      GO TO (1,2), NCODE
    1 A= (A-B)/(A+B)
      GO TO 3
    2 A= (A+A)/(A+B)
    3 RMOD= SQRT(REAL(A)*REAL(A)+AIMAG(A)*AIMAG(A))
      RPH= ATAN2(AIMAG(A),REAL(A))
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE COEF50(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,ND,RMOD,RPH)
C
C     ******************************************************************
C
C     The routine COEF50 is designed to evaluate the reflection and
C     transmission coefficients at a plane interface between two
C     homogeneous solid halfspaces or at a free surface of a homogeneous
C     solid halfspace.
C
C     The kinds of individual coefficients are specified by the
C     following numbers
C     a/ Interface between two solid halfspaces
C     P1P1...1       P1S1...2       P1P2...3       P1S2...4
C     S1P1...5       S1S1...6       S1P2...7       S1S2...8
C     b/ Free surface (for RO2.LT.0.00001)
C     PP.....1       PX.....5       PX,PZ...X- and Z- components of the
C     PS.....2       PZ.....6       coef.of conversion,incident P wave
C     SP.....3       SX.....7       SX,SZ...X- and Z- components of the
C     SS.....4       SZ.....8       coef.of conversion,incident S wave
C
C     I n p u t   p a r a m e t e r s
C           P...Ray parameter
C           VP1,VS1,RO1...Parameters of the first halfspace
C           VP2,VS2,RO2...Parameters of second halfspace. For the free
C                    surface take RO2.LT.0.000001,eg.RO2=0., and
C                    arbitrary values of VP2 and VS2
C           NCODE...Code of the computed coefficient
C           ND...=0  when the interface is situated at the right-hand
C                    side of the incident ray (X against P)
C                =1  when the interface is situated at the left-hand
C                    side of the incident ray (X along P)
C
C     O u t p u t   p a r a m e t e r s
C           RMOD,RPH...Modul and argument of the coefficient
C
C     N o t e s
C     1/ Positive P...In the direction of propagation
C     2/ Positive S...To the left from P
C     3/ Positive X...To the right from Z (to the right from P)
C     4/ Positive Z...In the direction of P
C     5/ Time factor of incident wave...EXP(-i*OMEGA*T)
C     6/ Formulae are taken from Cerveny and Ravindra, Theory of Seismic
C        Head Waves,Pages 63-67. Due to the note 2, the signs at certain
C        coefficients are opposite (and time factor is changed, L.K.)
C        (signs of conversion coefficients are opposite to Cerveny,
C        Molotkov and Psencik, Ray Method in Seismology, page 35)
C
C     ******************************************************************
C
      COMPLEX B(4),RR,C1,C2,C3,C4,H1,H2,H3,H4,H5,H6,H,HB,HC
      DIMENSION PRMT(4),D(4),DD(4)
      PRMT(1)=VP1
      PRMT(2)=VS1
      PRMT(3)=VP2
      PRMT(4)=VS2
      IF(RO2.LT.0.000001)GO TO 150
      A1=VP1*VS1
      A2=VP2*VS2
      A3=VP1*RO1
      A4=VP2*RO2
      A5=VS1*RO1
      A6=VS2*RO2
      Q=2.*(A6*VS2-A5*VS1)
      PP=P*P
      QP=Q*PP
      X=RO2-QP
      Y=RO1+QP
      Z=RO2-RO1-QP
      G1=A1*A2*PP*Z*Z
      G2=A2*X*X
      G3=A1*Y*Y
      G4=A4*A5
      G5=A3*A6
      G6=Q*Q*PP
      DO 21 I=1,4
      DD(I)=P*PRMT(I)
   21 D(I)=SQRT(ABS(1.-DD(I)*DD(I)))
      IF(DD(1).LE.1..AND.DD(2).LE.1..AND.DD(3).LE.1..AND.DD(4).LE.1.)
     1GO TO 100
C
C     Complex coefficients
      DO 22 I=1,4
      IF(DD(I).GT.1.)GO TO 23
      B(I)=CMPLX(D(I),0.)
      GO TO 22
   23 B(I)= CMPLX(0.,D(I))
   22 CONTINUE
      C1=B(1)*B(2)
      C2=B(3)*B(4)
      C3=B(1)*B(4)
      C4=B(2)*B(3)
      H1=CMPLX(G1,0.)
      H2=G2*C1
      H3=G3*C2
      H4=G4*C3
      H5=G5*C4
      H6=G6*C1*C2
      H=1./(H1+H2+H3+H4+H5+H6)
      HB=2.*H
      HC=HB*P
      GO TO (1,2,3,4,5,6,7,8),NCODE
    1 RR=H*(H2+H4+H6-H1-H3-H5)
      GO TO 26
    2 RR=VP1*B(1)*HC*(Q*Y*C2+A2*X*Z)
      IF(ND.NE.0)RR=-RR
      GO TO 26
    3 RR=A3*B(1)*HB*(VS2*B(2)*X+VS1*B(4)*Y)
      GO TO 26
    4 RR=-A3*B(1)*HC*(Q*C4-VS1*VP2*Z)
      IF(ND.NE.0)RR=-RR
      GO TO 26
    5 RR=-VS1*B(2)*HC*(Q*Y*C2+A2*X*Z)
      IF(ND.NE.0)RR=-RR
      GO TO 26
    6 RR=H*(H2+H5+H6-H1-H3-H4)
      GO TO 26
    7 RR=A5*B(2)*HC*(Q*C3-VP1*VS2*Z)
      IF(ND.NE.0)RR=-RR
      GO TO 26
    8 RR=A5*B(2)*HB*(VP1*B(3)*Y+VP2*B(1)*X)
      GO TO 26
C     Real coefficients
  100 E1=D(1)*D(2)
      E2=D(3)*D(4)
      E3=D(1)*D(4)
      E4=D(2)*D(3)
      S1=G1
      S2=G2*E1
      S3=G3*E2
      S4=G4*E3
      S5=G5*E4
      S6=G6*E1*E2
      S=1./(S1+S2+S3+S4+S5+S6)
      SB=2.*S
      SC=SB*P
      GO TO (101,102,103,104,105,106,107,108),NCODE
  101 R=S*(S2+S4+S6-S1-S3-S5)
      GO TO 250
  102 R=VP1*D(1)*SC*(Q*Y*E2+A2*X*Z)
      IF(ND.NE.0)R=-R
      GO TO 250
  103 R=A3*D(1)*SB*(VS2*D(2)*X+VS1*D(4)*Y)
      GO TO 250
  104 R=-A3*D(1)*SC*(Q*E4-VS1*VP2*Z)
      IF(ND.NE.0)R=-R
      GO TO 250
  105 R=-VS1*D(2)*SC*(Q*Y*E2+A2*X*Z)
      IF(ND.NE.0)R=-R
      GO TO 250
  106 R=S*(S2+S5+S6-S1-S3-S4)
      GO TO 250
  107 R=A5*D(2)*SC*(Q*E3-VP1*VS2*Z)
      IF(ND.NE.0)R=-R
      GO TO 250
  108 R=A5*D(2)*SB*(VP1*D(3)*Y+VP2*D(1)*X)
      GO TO 250
C
C     Earths surface,complex coefficients and coefficients of conversion
  150 A1=VS1*P
      A2=A1*A1
      A3=2.*A2
      A4=2.*A1
      A5=A4+A4
      A6=1.-A3
      A7=2.*A6
      A8=2.*A3*VS1/VP1
      A9=A6*A6
      DD(1)=P*VP1
      DD(2)=P*VS1
      DO 151 I=1,2
  151 D(I)=SQRT(ABS(1.-DD(I)*DD(I)))
      IF(DD(1).LE.1..AND.DD(2).LE.1.)GO TO 200
      DO 154 I=1,2
      IF(DD(I).GT.1.)GO TO 155
      B(I)=CMPLX(D(I),0.)
      GO TO 154
  155 B(I)= CMPLX(0.,D(I))
  154 CONTINUE
      H1=B(1)*B(2)
      H2=H1*A8
      H=1./(A9+H2)
      GO TO (161,162,163,164,165,166,167,168),NCODE
  161 RR=(-A9+H2)*H
      GO TO 26
  162 RR=-A5*B(1)*H*A6
      IF(ND.NE.0)RR=-RR
      GO TO 26
  163 RR=A5*B(2)*H*A6*VS1/VP1
      IF(ND.NE.0)RR=-RR
      GO TO 26
  164 RR=-(A9-H2)*H
      GO TO 26
  165 RR=-A5*H1*H
      IF(ND.NE.0)RR=-RR
      GO TO 26
  166 RR=A7*B(1)*H
      GO TO 26
  167 RR=-A7*B(2)*H
      GO TO 26
  168 RR=-A5*VS1*H1*H/VP1
      IF(ND.NE.0)RR=-RR
   26 Z2=REAL(RR)
      Z3=AIMAG(RR)
      IF(Z2.EQ.0..AND.Z3.EQ.0.)GO TO 157
      RMOD=SQRT(Z2*Z2+Z3*Z3)
      RPH=ATAN2(Z3,Z2)
      RETURN
  157 RMOD=0.
      RPH=0.
      RETURN
C
C     Earths surface,real coefficients and coefficients of conversion
  200 S1=D(1)*D(2)
      S2=A8*S1
      S=1./(A9+S2)
      GO TO (201,202,203,204,205,206,207,208),NCODE
  201 R=(-A9+S2)*S
      GO TO 250
  202 R=-A5*D(1)*S*A6
      IF(ND.NE.0)R=-R
      GO TO 250
  203 R=A5*D(2)*S*A6*VS1/VP1
      IF(ND.NE.0)R=-R
      GO TO 250
  204 R=(S2-A9)*S
      GO TO 250
  205 R=-A5*S1*S
      IF(ND.NE.0)R=-R
      GO TO 250
  206 R=A7*D(1)*S
      GO TO 250
  207 R=-A7*D(2)*S
      GO TO 250
  208 R=-A5*VS1*S1*S/VP1
      IF(ND.NE.0)R=-R
  250 IF(R.LT.0.)GO TO 251
      RMOD=R
      RPH=0.
      RETURN
  251 RMOD=-R
      RPH=3.14159
      RETURN
      END
C
C=======================================================================
C