C
      SUBROUTINE COEF(STU,UC,ITRANS)
C
C  ROUTINE FOR THE COMPUTATION OF REFLECTION/TRANSMISSION COEFFICIENTS
C  FOR GENERAL ISOTROPIC/ANISOTROPIC MEDIA
C
C  INPUT PARAMETERS:
C  STU(1-3)... CARTESIAN COMPONENTS OF THE DISPLACEMENT VECTOR
C  STU(4-6)... CARTESIAN COMPONETS OF THE TRACTION VECTOR CALCULATED
C              IN THIS ROUTINE
C  ITRANS... ITRANS=O: REFLECTION
C            ITRANS=1: TRANSMISSION
C            ITRANS=999: SURFACE CONVERSION
C
C  OUTPUT PARAMETERS:
C  UC(1-3)...  RAY-CENTERED COMPONENTS OF THE DISPLACEMENT VECTOR
C
      COMPLEX A(6,6),C(6,6),COA(3,3),COC(3,3),AA(6),TAU(3),P(3),U(3),
     1CP1,CP2,CP3,DETA,DETC,UC(3),STU(6),CAB
      DIMENSION Y(18),UN(3)
      COMMON /AUXI/IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      INTEGER CODE
      COMMON /COD/CODE(50,2),KREF,KC,ITYPE
      COMPLEX PS
      COMMON /RAY/   AY(28,2000),DS(20,50),KINT(50),HHH(3,3),TMAX,
     1               PS(3,7,50),IS(8,50),N,IREF,IND,IND1
      COMMON /ZERO/ RNULL
C
      INDI=0
      IF(ITRANS.EQ.999) THEN
        ITRANS=0
        INDI=1
      END IF
      I1=1
      I2=3
      LAY1=IS(1,IREF)
      LAY=IS(7,IREF)
      DO 910 K=1,6
      IF(K.GT.3) GOTO 900
      UN(K)=DS(K,IREF)
  900 Y(K)=AY(K+1,N)
  910 CONTINUE
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,/,3(3F14.7,/),3I5)')'COEF: Y(1)-Y(6),UN,IREF,LAY,
     2ITYPE',(Y(K),K=1,6),UN,IREF,LAY,ITYPE
      CALL PARDIS (Y,0)
C
C  DISPLACEMENT AND STRESS VECTOR OF INCIDENT WAVE
C
      BSTU=SQRT(REAL(STU(1)*CONJG(STU(1))+STU(2)*CONJG(STU(2))
     1               +STU(3)*CONJG(STU(3))))
      BST=STU(1)*HHH(ITYPE,1)+STU(2)*HHH(ITYPE,2)
     1               +STU(3)*HHH(ITYPE,3)
      BSTU=SIGN(1.,BST)*BSTU
      DO 950 K=1,3
      P(K)=PS(K,7,IREF)
      STU(K)=STU(K)/BSTU
      U(K)=STU(K)
  950 CONTINUE
      CALL STRESS (P,UN,U,TAU)
      DO 960 K=1,3
      STU(K+3)=TAU(K)
  960 CONTINUE
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,/,6(2F14.7,/))') 'COEF: INCIDENT U,TAU',U,TAU
 1000 IF(IANI(LAY).EQ.0) GOTO 3000
C
C  COMPUTATION OF STRESS-DISPLACEMENT-VECTORS FOR GENERATED WAVES
C
C  ANISOTROPIC CASE
C
 1010 IF(I1.GT.3)CALL PARDIS(Y,0)
      DO 1020 K=I1,I2
      DO 1030 L=1,3
      P(L)=PS(L,K,IREF)
 1030 CONTINUE
C
C  EVALUATION OF DISPLACEMENT VECTOR (GENERALLY COMPLEX-VALUED)
C
      CALL DISPL(P,U,K)
      CALL STRESS(P,UN,U,TAU)
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7,/))')'COEF: P,U,TAU'
     1,P,U,TAU
      DO 1060 L=1,3
      A(K,L)=U(L)
      A(K,L+3)=TAU(L)
 1060 CONTINUE
 1020 CONTINUE
      IF(INTR.EQ.1.AND.INDI.EQ.1) GOTO 5000
      LAY=IS(8,IREF)
      IF(LAY.EQ.0) THEN
         LAY=IS(7,IREF)
         GO TO 5000
      END IF
      IF(I2.EQ.6) GOTO 4000
      I1=4
      I2=6
      IF(IANI(LAY).NE.0) GOTO 1010
C
C  ISOTROPIC CASE
C
 3000 IF(I1.GT.3) CALL PARDIS(Y,0)
      DO 3010 L=1,3
      P(L)=PS(L,I1,IREF)
 3010 CONTINUE
      CP1=P(1)
      CP2=P(2)
      IVERT=0
      IF(CABS(CP1).LT..0000001.AND.CABS(CP2).LT..0000001)IVERT=1
      IF(IVERT.EQ.1)THEN
        U(1)=1.
        U(2)=0.
        U(3)=0.
        GO TO 3020
      END IF
      CP3=-(CP1*CP1+CP2*CP2)
      CP1=CP1*P(3)
      CP2=CP2*P(3)
      CAB=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3)))
      U(1)=CP1/CAB
      U(2)=CP2/CAB
      U(3)=CP3/CAB
 3020 CONTINUE
      CALL STRESS(P,UN,U,TAU)
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7/))')'COEF: P,U,TAU'
     1,P,U,TAU
      DO 3030 L=1,3
      A(I1,L)=U(L)
      A(I1,L+3)=TAU(L)
 3030 CONTINUE
      IF(IVERT.EQ.1)THEN
        U(1)=0.
        U(2)=1.
        IF(REAL(P(3)).LT.0.)U(2)=-1.
        U(3)=0.
        GO TO 3040
      END IF
      CP1=-P(2)
      CP2=P(1)
      CP3=CMPLX(0.,0.)
      CAB=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3)))
      U(1)=CP1/CAB
      U(2)=CP2/CAB
      U(3)=CP3/CAB
 3040 CONTINUE
      CALL STRESS (P,UN,U,TAU)
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7/))')'COEF: P,U,TAU'
     1,P,U,TAU
      DO 3050 L=1,3
      A(I1+1,L)=U(L)
      A(I1+1,L+3)=TAU(L)
 3050 CONTINUE
      DO 3210 L=1,3
      P(L)=PS(L,I2,IREF)
 3210 CONTINUE
      CP1=P(1)
      CP2=P(2)
      CP3=P(3)
      CAB=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3)))
      U(1)=CP1/CAB
      U(2)=CP2/CAB
      U(3)=CP3/CAB
      CALL STRESS (P,UN,U,TAU)
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7/))')'COEF: P,U,TAU'
     1,P,U,TAU
      DO 3220 L=1,3
      A(I2,L)=U(L)
      A(I2,L+3)=TAU(L)
3220  CONTINUE
      IF(INTR.EQ.1.AND.INDI.EQ.1) GOTO 5000
      LAY=IS(8,IREF)
      IF(LAY.EQ.0) THEN
         LAY=IS(7,IREF)
         GO TO 5000
      END IF
      IF(I2.EQ.6)GOTO 4000
      I1=4
      I2=6
      GOTO 1000
C
C  REVERSE OF SIGNS FOR REFLECTIONS
C
 4000 DO 4010 L=1,3
      DO 4010 K=1,6
      A(L,K)=-A(L,K)
 4010 CONTINUE
      IF(ICOEF.EQ.0) GOTO 4012
      WRITE(LOUT,'(//,A)') 'COEF: MATRIX A'
      DO 4011 K=1,6
      WRITE(LOUT,'(6(12F10.5))') (A(L,K),L=1,6)
 4011 CONTINUE
C
C  SOLUTION OF SYSTEM OF LINEAR EQUATIONS (KRAMER'S METHOD)
C
 4012 CALL DET6(A,DETA)
      DO 4020 K=1,3
      UC(K)=CMPLX(0.,0.)
 4020 CONTINUE
      IF(CABS(DETA).LT.0.0000001) THEN
        WRITE(LOUT,'(A)')
     1  'COEF: MATRIX FOR R/T COEF. SINGULAR: AMPLITUDE NOT EVALUATED'
        IND=11
        RETURN
      END IF
      IF(INDI.EQ.1) THEN
         L1=1
         L2=3
         GO TO 4041
      END IF
      IF(IANI(LAY1).NE.0.OR.ITYPE.EQ.3)THEN
         L1=ITYPE
         IF(ITRANS.NE.0)L1=L1+3
         L2=L1
      END IF
      IF(IANI(LAY1).EQ.0.AND.ITYPE.NE.3)THEN
         L1=1
         IF(ITRANS.NE.0)L1=L1+3
         L2=L1+1
      END IF
 4041 CONTINUE
      IF(ICOEF.NE.0) WRITE(LOUT,'(A,2I5)')'COEF: L1,L2',L1,L2
      DO 4044 L=L1,L2
      DO 4043 J=1,6
      DO 4042 K=1,6
      C(J,K)=A(J,K)
      IF(L.EQ.J)C(J,K)=STU(K)
 4042 CONTINUE
 4043 CONTINUE
      CALL DET6(C,DETC)
      AA(L)=DETC/DETA
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,2F12.6)') 'COEF: AA=',AA(L)
 4044 CONTINUE
      IF(INDI.EQ.1)GO TO 5045
C
      DO 4100 L=L1,L2
      LL=L
      IF(ITRANS.NE.0)LL=L-3
      UC(LL)=AA(L)*BSTU
 4100 CONTINUE
      RETURN
C
C  INTERACTION WITH THE FREE SURFACE
C
 5000 DO 5010 L=1,3
      DO 5010 K=1,3
      A(L,K)=-A(L,K)
      COA(L,K)=-A(L,K+3)
 5010 CONTINUE
      IF(ICOEF.EQ.0)GO TO 5012
      WRITE(LOUT,'(//,A)')'COEF: MATRIX COA'
      WRITE(LOUT,'(3(6F10.5/))')((COA(L,K),L=1,3),K=1,3)
C
C  SOLUTION OF SYSTEM OF 3 LINEAR EQUATIONS (CRAMER'S METHOD)
C
 5012 CALL DET3(COA,DETA)
      IF(CABS(DETA).LT.0.0000001)THEN
        WRITE(LOUT,'(A)')
     1  'COEF: MATRIX FOR R/T COEF. SINGULAR: AMPLITUDE NOT EVALUATED'
        IND=11
        RETURN
      END IF
      DO 5044 L=1,3
      DO 5043 J=1,3
      DO 5042 K=1,3
      COC(J,K)=COA(J,K)
      IF(L.EQ.J)COC(J,K)=STU(K+3)
 5042 CONTINUE
 5043 CONTINUE
      CALL DET3(COC,DETC)
      AA(L)=DETC/DETA
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,2F12.6)')'COEF: AA=',AA(L)
 5044 CONTINUE
      IF(INDI.NE.1)GO TO 5200
C
C     CONVERSION AT AN INTERFACE
C
 5045 CONTINUE
      IF(MREG.EQ.2) THEN
C
C     CALCULATION OF CONVERSION FOR PRESSURE
C
        CP1=STU(1)
        CP2=STU(2)
        CP3=STU(3)
        ARE=REAL(CP1)
        AIM=AIMAG(CP1)
        APHI=ATAN2(AIM,ARE)
        ARE=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3)))
        UC(1)=ARE*(1.+AA(3))*CMPLX(COS(APHI),SIN(APHI))
        UC(2)=(0.,0.)
        UC(3)=(0.,0.)
        IF(ICOEF.GT.0)WRITE(LOUT,'(A,4F10.5)') 'COEF: UC(1),ARE,APHI',
     1  UC(1),ARE,APHI
        GO TO 5120
      END IF
C
      DO 5100 K=1,3
      U(K)=CMPLX(0.,0.)
      DO 5050 J=1,3
      U(K)=-AA(J)*A(J,K)+U(K)
 5050 CONTINUE
      UC(K)=U(K)+STU(K)
 5100 CONTINUE
 5120 CONTINUE
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,3(2F12.5/))') 'COEF:
     1 NORMALIZED CONVERSION-VECTOR',UC
      DO 5150 K=1,3
      UC(K)=UC(K)*BSTU
 5150 CONTINUE
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,3(2F12.5/))') 'COEF:
     1 CONVERSION-VECTOR',UC
      RETURN
C
C     REFLECTION FROM THE FREE SURFACE
C
 5200 L1=1
      IF(ITYPE.EQ.3)L1=3
      L2=2
      IF(ITYPE.EQ.3)L2=3
      DO 5250 K=1,3
      UC(K)=CMPLX(0.,0.)
 5250 CONTINUE
      DO 5300 L=L1,L2
      UC(L)=AA(L)*BSTU
 5300 CONTINUE
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE DET6(A,D6)
      COMPLEX A(6,6),B(5,5),D6,D
      D6=(0.,0.)
      DO 10 J=1,6
      IF(A(1,J).NE.(0.,0.)) THEN
      CALL SET(A,B,1,J,6)
      CALL DET5(B,D)
      C=(-1)**(J+1)
      D6=D6+CMPLX(C,0.)*D*A(1,J)
      ENDIF
  10  CONTINUE
      RETURN
      END
C
C     *********************************************************
C
      SUBROUTINE DET5(A,D5)
      COMPLEX A(5,5),B(4,4),D5,D
      D5=(0.,0.)
      DO 10 J=1,5
      IF(A(1,J).NE.(0.,0.)) THEN
      CALL SET(A,B,1,J,5)
      CALL DET4(B,D)
      C=(-1)**(J+1)
      D5=D5+CMPLX(C,0.)*D*A(1,J)
      ENDIF
  10  CONTINUE
      RETURN
      END
C
C     ***********************************************************
C
      SUBROUTINE DET4(A,D4)
      COMPLEX A(4,4),B(3,3),D4,D
      D4=(0.,0.)
      DO 10 J=1,4
      IF(A(1,J).NE.(0.,0.)) THEN
      CALL SET(A,B,1,J,4)
      CALL DET3(B,D)
      C=(-1)**(J+1)
      D4=D4+CMPLX(C,0.)*D*A(1,J)
      ENDIF
  10  CONTINUE
      RETURN
      END
C
C     **********************************************************
C
      SUBROUTINE DET3(A,D3)
      COMPLEX A(3,3),D3
      D3=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))
     1+  A(1,2)*(A(2,3)*A(3,1)-A(2,1)*A(3,3))
     2+  A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1))
      RETURN
      END
C
C     ************************************************************
C
      SUBROUTINE SET(A,B,I,J,K)
      COMPLEX A(K,K),B(K-1,K-1)
      IFL=0
      JFL=0
      DO 20 I1=1,K
      IF(I1.EQ.I) THEN
      IFL=1
      GOTO 20
      ENDIF
      DO 10 J1=1,K
      IF(J1.EQ.J) THEN
      JFL=1
      GOTO 10
      ENDIF
      IF(IFL.EQ.0.AND.JFL.EQ.0) THEN
      B(I1,J1)=A(I1,J1)
      ELSE IF(IFL.EQ.0.AND.JFL.EQ.1) THEN
      B(I1,J1-1)=A(I1,J1)
      ELSE IF(IFL.EQ.1.AND.JFL.EQ.0) THEN
      B(I1-1,J1)=A(I1,J1)
      ELSE
      B(I1-1,J1-1)=A(I1,J1)
      ENDIF
  10  CONTINUE
      JFL=0
  20  CONTINUE
      RETURN
      END
C
C     **********************************************************
C
      SUBROUTINE DISC(Y,DEP)
C
C     DETERMINATION OF DEPTH OF 3D INTERFACES AND ITS DERIVATIVES
C     FOR BICUBIC POLYNOMIAL APPROXIMATION
C
      SAVE B,DX,DY,LLAY,IBB,IU,IL,JU,JL
      DIMENSION Y(18),DEP(6),B(16,2),DX(2),DY(2)
C
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      COMMON /AUXX/  MMX(20),MMY(20),MMXY(20)
      COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT,
     1   XINTA
      COMMON/ZCOEF/ A02(1000),A20(1000),A22(1000)
C
      IBB=0
      LB=2
      IF(INTR.EQ.LAY)LB=1
      MX=NX(INTR)
      MY=NY(INTR)
      DO 1 I=2,MX
      K=MMX(INTR)+I-1
      IF(Y(1).LT.SX(K))GO TO 2
    1 CONTINUE
    2 I1=K
      DO 3 I=2,MY
      K=MMY(INTR)+I-1
      IF(Y(2).LT.SY(K))GO TO 4
    3 CONTINUE
    4 J1=K
      IF(MAUX.EQ.0) GOTO 8
      IF(LB.EQ.2) GOTO 5
      IF(I1.EQ.IU.AND.J1.EQ.JU.AND.LLAY.EQ.LAY) GOTO 10
      GOTO 8
    5 IF(I1.EQ.IL.AND.J1.EQ.JL.AND.LLAY.EQ.LAY) GOTO 10
      IL=I1
      JL=J1
      GOTO 9
    8 IU=I1
      JU=J1
    9 LLAY=LAY
      I=I1-MMX(INTR)
      J=J1-MMY(INTR)
      DX(LB)=SX(I1-1)
      DY(LB)=SY(J1-1)
      MM=MMXY(INTR)-1
      K=MM+(I-1)*MY+J
      B20=A20(K)
      B02=A02(K)
      B22=A22(K)
      B00=Z(K)
      K=MM+I*MY+J
      C20=A20(K)
      C02=A02(K)
      C22=A22(K)
      C00=Z(K)
      K=MM+(I-1)*MY+J+1
      D20=A20(K)
      D02=A02(K)
      D22=A22(K)
      D00=Z(K)
      K=MM+I*MY+J+1
      E20=A20(K)
      E02=A02(K)
      E22=A22(K)
      E00=Z(K)
      HX=SX(I1)-DX(LB)
      HY=SY(J1)-DY(LB)
      XA=3.*HX
      YA=3.*HY
      XB=HX/3.
      YB=HY/3.
      D32=(E22-D22)/XA
      D30=(E20-D20)/XA
      B30=(C20-B20)/XA
      B32=(C22-B22)/XA
      D12=(E02-D02)/HX-XB*(E22+2.*D22)
      D10=(E00-D00)/HX-XB*(E20+2.*D20)
      B10=(C00-B00)/HX-XB*(C20+2.*B20)
      B12=(C02-B02)/HX-XB*(C22+2.*B22)
      B03=(D02-B02)/YA
      B13=(D12-B12)/YA
      B23=(D22-B22)/YA
      B33=(D32-B32)/YA
      B01=(D00-B00)/HY-YB*(D02+2.*B02)
      B11=(D10-B10)/HY-YB*(D12+2.*B12)
      B21=(D20-B20)/HY-YB*(D22+2.*B22)
      B31=(D30-B30)/HY-YB*(D32+2.*B32)
      MAUX=1
      B(1,LB)=B00
      B(2,LB)=B01
      B(3,LB)=B02
      B(4,LB)=B03
      B(5,LB)=B10
      B(6,LB)=B11
      B(7,LB)=B12
      B(8,LB)=B13
      B(9,LB)=B20
      B(10,LB)=B21
      B(11,LB)=B22
      B(12,LB)=B23
      B(13,LB)=B30
      B(14,LB)=B31
      B(15,LB)=B32
      B(16,LB)=B33
      IBB=1
   10 AX=Y(1)-DX(LB)
      AZ=Y(2)-DY(LB)
      IF(IBB.EQ.1) GOTO 11
      B00=B(1,LB)
      B01=B(2,LB)
      B02=B(3,LB)
      B03=B(4,LB)
      B10=B(5,LB)
      B11=B(6,LB)
      B12=B(7,LB)
      B13=B(8,LB)
      B20=B(9,LB)
      B21=B(10,LB)
      B22=B(11,LB)
      B23=B(12,LB)
      B30=B(13,LB)
      B31=B(14,LB)
      B32=B(15,LB)
      B33=B(16,LB)
   11 AUX1=((B33*AZ+B32)*AZ+B31)*AZ+B30
      AUX2=((B23*AZ+B22)*AZ+B21)*AZ+B20
      AUX3=((B13*AZ+B12)*AZ+B11)*AZ+B10
      AUX4=((B03*AZ+B02)*AZ+B01)*AZ+B00
      DEP(1)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4
      IF(NDER.EQ.0) RETURN
      DEP(2)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3
      IF(NDER.EQ.1)GO TO 7
      DEP(4)=6.*AUX1*AX+2.*AUX2
    7 AUX1=(3.*B33*AZ+2.*B32)*AZ+B31
      AUX2=(3.*B23*AZ+2.*B22)*AZ+B21
      AUX3=(3.*B13*AZ+2.*B12)*AZ+B11
      AUX4=(3.*B03*AZ+2.*B02)*AZ+B01
      DEP(3)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4
      IF(NDER.EQ.1) RETURN
      DEP(5)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3
      AUX1=3.*B33*AZ+B32
      AUX2=3.*B23*AZ+B22
      AUX3=3.*B13*AZ+B12
      AUX4=3.*B03*AZ+B02
      DEP(6)=2.*(((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4)
      RETURN
      END
C
C     **********************************************************
C
      SUBROUTINE DISPL(P,U,I1)
C
C  ROUTINE FOR THE COMPUTATION OF DISPLACEMENT VECTOR OF GENERATED
C  WAVES AT AN INTERFACE. IT WORKS EVEN FOR A COMPLEX REFRACTION VECTOR.
C  ROUTINES CHRM AND POLAR WORK ONLY FOR REAL REFRACTION VECTORS
C
C  INPUT PARAMETERS:
C  P(1-3)... SLOWNESS VECTOR
C  I1... TYPE OF THE WAVE WITH THE SLOWNESS VECTOR P
C        I1=1... REFLECTED S WAVE (SV COMPONENT IN ISOTROPIC CASE)
C        I1=2... REFLECTED S WAVE (SH COMPONENT IN ISOTROPIC CASE)
C        I1=3... REFLECTED P WAVE
C        I1=4... TRANSMITTED S WAVE (SV COMPONENT IN ISOTROPIC CASE)
C        I1=5... TRANSMITTED S WAVE (SH COMPONENT IN ISOTROPIC CASE)
C        I1=6... TRANSMITTED P WAVE
C
C  OUTPUT PARAMETERS:
C  U(1-3)... CORRESPONDING DISPLACEMENT VECTOR
C
      complex uxn(3,3)
      COMPLEX U,U1,U2,U3,P,P1,P2,P3,P2P3,P1P2,P1P3,P1P1,P2P2,P3P3,C11,
     1  C12,C13,C22,C23,C33,C11N,C22N,C33N,C23SQ,C13SQ,C12SQ,CD11,CD12,
     2  CD13,CD22,CD23,CD33,CZ01,CZ02,CZ03,CDTR
      DIMENSION U(3),P(3),xx(3)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      COMPLEX PS
      COMMON /RAY/   AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax,
     1               PS(3,7,50),IS(8,50),N,IREF,IND,IND1
C
C
C  THIS PART CORRESPONDS TO CHRM ,BUT FOR COMPLEX REFRACTION VECTORS
C
 900  P1=P(1)
      P2=P(2)
      P3=P(3)
 910  continue
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: START P',
     1p1,p2,p3
      P2P3=P2*P3
      P1P2=P1*P2
      P1P3=P1*P3
      P1P1=P1*P1
      P2P2=P2*P2
      P3P3=P3*P3
      C11=P1P1*A11+P2P2*A66+P3P3*A55
     1+2.*P2P3*A56+2.*P1P3*A15+2.*P1P2*A16
      C22=P1P1*A66+P2P2*A22+P3P3*A44
     1+2.*P2P3*A24+2.*P1P3*A46+2.*P1P2*A26
      C33=P1P1*A55+P2P2*A44+P3P3*A33
     1+2.*P2P3*A34+2.*P1P3*A35+2.*P1P2*A45
      C23=P1P1*A56+P2P2*A24+P3P3*A34
     1   +P2P3*A2344+P1P3*A3645+P1P2*A2546
      C13=P1P1*A15+P2P2*A46+P3P3*A35
     1   +P2P3*A3645+P1P3*A1355+P1P2*A1456
      C12=P1P1*A16+P2P2*A26+P3P3*A45
     1   +P2P3*A2546+P1P3*A1456+P1P2*A1266
      C11N=C11-1.
      C22N=C22-1.
      C33N=C33-1.
      C23SQ=C23*C23
      C13SQ=C13*C13
      C12SQ=C12*C12
      CD11=C22N*C33N-C23SQ
      CD22=C11N*C33N-C13SQ
      CD33=C11N*C22N-C12SQ
      CD12=C13*C23-C12*C33N
      CD13=C12*C23-C13*C22N
      CD23=C12*C13-C23*C11N
      CDTR=CD11+CD22+CD33
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,4(/4F12.6))')'DISPL: CDIJ'
     1,CD11,CD12,CD13,CD22,CD23,CD33,CDTR
c     if(abs(cdtr).lt..000001)write(lout,'(A)')'DISPL: SHEAR WAVE
c    1 SINGULARITY IN CALCULATION OF R/T COEFFICIENT'
C
C  THIS PART CORRESPONDS TO ROUTINE POLAR BUT FOR COMPLEX POLARIZATION
C  VECTORS
C
      U1=CD11
      U2=CD12
      U3=CD13
      uxn(1,1)=u1
      uxn(1,2)=u2
      uxn(1,3)=u3
      XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3))
      xx(1)=xn
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN
      U1=CD12
      U2=CD22
      U3=CD23
      uxn(2,1)=u1
      uxn(2,2)=u2
      uxn(2,3)=u3
      XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3))
      xx(2)=xn
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN
      U1=CD13
      U2=CD23
      U3=CD33
      uxn(3,1)=u1
      uxn(3,2)=u2
      uxn(3,3)=u3
      XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3))
      xx(3)=xn
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN
      xn=0.
      do 920 l=1,3
        if(xn.ge.xx(l))go to 920
        xn=xx(l)
        nx=l
  920 continue
      xn=1./sqrt(xn)
      u1=uxn(nx,1)*xn
      u2=uxn(nx,2)*xn
      u3=uxn(nx,3)*xn
      U(1)=U1
      U(2)=U2
      U(3)=U3
      IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: U',U
C
C  CHECK OF PRECISION
C
      CZ01=C11N*u1+C12*u2+C13*u3
      CZ02=C12*u1+C22N*u2+C23*u3
      CZ03=C13*u1+C23*u2+C33N*u3
      IF(ICOEF.GT.0)
     1WRITE(LOUT,'(A,3(/2F14.7))')'DISPL: PRECISSION OF DISPL'
     1,CZ01,CZ02,CZ03
      RETURN
      END
C
C     **********************************************************
C
      SUBROUTINE DMAT
C
C     EVALUATES ELEMENTS OF THE MATRIX DIJ
C
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMMON/DJK/D11,D12,D13,D22,D23,D33,DTR
      COMMON/GAM/G11,G12,G13,G22,G23,G33
      COMMON/AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      COMPLEX PS
      COMMON /RAY/   AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax,
     1               PS(3,7,50),IS(8,50),N,IREF,IND,IND1
C
      G1=G11-1.
      G2=G22-1.
      G3=G33-1.
      D11=G2*G3-G23*G23
      D22=G1*G3-G13*G13
      D33=G1*G2-G12*G12
      D12=G13*G23-G3*G12
      D13=G12*G23-G2*G13
      D23=G12*G13-G1*G23
      DTR=D11+D22+D33
      IF(ABS(DTR).LT.0.0001)THEN
        write(lout,'(A)')'DMAT: SHEAR WAVE SINGULARITY'
        IND=10
      END IF
      RETURN
      END
C
C     **********************************************************
C
      SUBROUTINE DDMAT(DG,g)
C
C     EVALUATES DERIVATIVES OF ELEMENTS OF THE MATRIX DIJ
C
      DIMENSION DG(3,3)
      COMMON/DDJK/DD11,DD12,DD13,DD22,DD23,DD33,DDTR
      COMMON/GAM/G11,G12,G13,G22,G23,G33
C
      DG11=DG(1,1)
      DG22=DG(2,2)
      DG33=DG(3,3)
      DG12=DG(1,2)
      DG13=DG(1,3)
      DG23=DG(2,3)
      G1=G11-1.
      G2=G22-1.
      G3=G33-1.
      DD11=G3*DG22+G2*DG33-2.*G23*DG23+g*(2.-g22-g33)
      DD22=G3*DG11+G1*DG33-2.*G13*DG13+g*(2.-g33-g11)
      DD33=G2*DG11+G1*DG22-2.*G12*DG12+g*(2.-g11-g22)
      DD12=G23*DG13+G13*DG23-DG12*G3-G12*DG33+g*g12
      DD13=G23*DG12+G12*DG23-DG13*G2-G13*DG22+g*g13
      DD23=G13*DG12+G12*DG13-DG23*G1-G23*DG11+g*g23
      DDTR=DD11+DD22+DD33
      RETURN
      END
C
C     **********************************************************
C
      FUNCTION TR(G)
C
C     TRACE OF PRODUCT OF MATRICES G AND D
C
      DIMENSION G(3,3)
      COMMON/DJK/ D11,D12,D13,D22,D23,D33,DTR
C
      TR=G(1,1)*D11+G(2,2)*D22+G(3,3)*D33+
     1    2.*(G(1,2)*D12+G(1,3)*D13+G(2,3)*D23)
      RETURN
      END
C
C     **********************************************************
C
      FUNCTION TRD(G)
C
C     TRACE OF PRODUCT OF MATRICES G AND DD
C
      DIMENSION G(3,3)
      COMMON/DDJK/ DD11,DD12,DD13,DD22,DD23,DD33,DDTR
C
      TRD=G(1,1)*DD11+G(2,2)*DD22+G(3,3)*DD33+
     1    2.*(G(1,2)*DD12+G(1,3)*DD13+G(2,3)*DD23)
      RETURN
      END
C
C     **********************************************************
C
      SUBROUTINE FCT(X,Y,DERY)
C
C  COMPUTATION OF THE RIGHT HAND SIDES OF THE RAY TRACING EQUATIONS
C
      save g,gpp,gpx,gxx,gx,gp,gx1,gx2,gx3,gp1,gp2,gp3
      DIMENSION DERY(18),Y(18),AUX(2),AUXX(3,2),VX(3)
      DIMENSION G(3,3),GPP(3,3),GPX(3,3),GXX(3,3),GX(3),GP(3)
      DIMENSION GX1(3,3),GX2(3,3),GX3(3,3),GP1(3,3),GP2(3,3),GP3(3,3)
      COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,
     1               A34,A35,A36,A44,A45,A46,A55,A56,A66,
     1               DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23,
     1               DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44,
     1               DXA45,DXA46,DXA55,DXA56,DXA66,
     1               DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23,
     1               DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44,
     1               DYA45,DYA46,DYA55,DYA56,DYA66,
     1               DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23,
     1               DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44,
     1               DZA45,DZA46,DZA55,DZA56,DZA66,
     1               A2546,A1266,A1355,A1456,A3645,A2344
      COMMON /APROX1/ E(21),EX(21),EY(21),EZ(21),EXX(21),EXY(21),
     1    EXZ(21),EYY(21),EYZ(21),EZZ(21)
      COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
     1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT,
     2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori
      INTEGER CODE
      COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
      COMPLEX PS
      COMMON /RAY/   AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax,
     1               PS(3,7,50),IS(8,50),N,IREF,IND,IND1
      COMMON /DJK/   D11,D12,D13,D22,D23,D33,DTR
      COMMON /DDJK/  DD11,DD12,DD13,DD22,DD23,DD33,DDTR
      COMMON /ZERO/ RNULL
C
      IF(IANI(LAY).NE.0) GOTO 20
      CALL PARDIS(Y,0)
      IF(ITYPE.LT.3)ITT=16
      IF(ITYPE.EQ.3)ITT=1
C
C  COMPUTATION OF RIGHT HAND SIDES OF RAY TRACING EQUATIONS
C  IN AN ISOTROPIC LAYER
C
C  P-WAVE FOR ITT=1, S-WAVE FOR ITT=16
C
      V=E(ITT)
      VX(1)=EX(ITT)
      VX(2)=EY(ITT)
      VX(3)=EZ(ITT)
      DO 1 I=1,3
        DERY(I)=V*Y(3+I)
        DERY(3+I)=-.5*VX(I)/V
    1 CONTINUE
      IF(NDER.EQ.1)RETURN
      VXX=EXX(ITT)
      VXY=EXY(ITT)
      VXZ=EXZ(ITT)
      VYY=EYY(ITT)
      VYZ=EYZ(ITT)
      VZZ=EZZ(ITT)
      DO 3 J=1,2
        AUX(J)=0.
        DO 2 I=1,3
          II=3+I+3*J
          AUX(J)=AUX(J)+VX(I)*Y(II)
    2   CONTINUE
        JJ=4+3*J
        AUXX(1,J)=VXX*Y(JJ)+VXY*Y(JJ+1)+VXZ*Y(JJ+2)
        AUXX(2,J)=VXY*Y(JJ)+VYY*Y(JJ+1)+VYZ*Y(JJ+2)
        AUXX(3,J)=VXZ*Y(JJ)+VYZ*Y(JJ+1)+VZZ*Y(JJ+2)
    3 CONTINUE
      DO 4 I=1,3
        DERY(I+6)=AUX(1)*Y(I+3)+V*Y(I+12)
        DERY(I+9)=AUX(2)*Y(I+3)+V*Y(I+15)
        DERY(I+12)=.5*(VX(I)*AUX(1)-V*AUXX(I,1))/V/V
        DERY(I+15)=.5*(VX(I)*AUX(2)-V*AUXX(I,2))/V/V
    4 CONTINUE
      RETURN
C
C
C  COMPUTATION OF RIGHT HAND SIDES OF RAY TRACING EQUATIONS
C  IN AN ANISOTROPIC LAYER
C
C  DETERMINATION OF PARAMETERS OF THE MEDIUM
C
  20  CALL PARDIS(Y,0)
      CALL CHRM2(Y,G,1)
      CALL DMAT
      IF(IND.EQ.10)RETURN
      CALL PCHRM(Y,GP1,1,1)
      CALL PCHRM(Y,GP2,2,1)
      CALL PCHRM(Y,GP3,3,1)
      CALL CHRM2(Y,GX1,2)
      CALL CHRM2(Y,GX2,3)
      CALL CHRM2(Y,GX3,4)
      GP(1)=TR(GP1)/DTR
      GP(2)=TR(GP2)/DTR
      GP(3)=TR(GP3)/DTR
      GX(1)=TR(GX1)/DTR
      GX(2)=TR(GX2)/DTR
      GX(3)=TR(GX3)/DTR
      DO 5 I=1,3
        DERY(I)=.5*GP(I)
        DERY(I+3)=-.5*GX(I)
    5 CONTINUE
      IF(NDER.EQ.1)RETURN
      CALL CHRM2(Y,G,1)
      CALL DDMAT(GP1,gp(1))
      GPP(1,1)=TRD(GP1)-GP(1)*DDTR
      GPP(1,2)=TRD(GP2)-GP(2)*DDTR
      GPP(1,3)=TRD(GP3)-GP(3)*DDTR
      GPX(1,1)=TRD(GX1)-GX(1)*DDTR
      GPX(1,2)=TRD(GX2)-GX(2)*DDTR
      GPX(1,3)=TRD(GX3)-GX(3)*DDTR
      CALL DDMAT(GP2,gp(2))
      GPP(2,1)=TRD(GP1)-GP(1)*DDTR
      GPP(2,2)=TRD(GP2)-GP(2)*DDTR
      GPP(2,3)=TRD(GP3)-GP(3)*DDTR
      GPX(2,1)=TRD(GX1)-GX(1)*DDTR
      GPX(2,2)=TRD(GX2)-GX(2)*DDTR
      GPX(2,3)=TRD(GX3)-GX(3)*DDTR
      CALL DDMAT(GP3,gp(3))
      GPP(3,1)=TRD(GP1)-GP(1)*DDTR
      GPP(3,2)=TRD(GP2)-GP(2)*DDTR
      GPP(3,3)=TRD(GP3)-GP(3)*DDTR
      GPX(3,1)=TRD(GX1)-GX(1)*DDTR
      GPX(3,2)=TRD(GX2)-GX(2)*DDTR
      GPX(3,3)=TRD(GX3)-GX(3)*DDTR
      CALL DDMAT(GX1,gx(1))
      GXX(1,1)=TRD(GX1)-GX(1)*DDTR
      GXX(1,2)=TRD(GX2)-GX(2)*DDTR
      GXX(1,3)=TRD(GX3)-GX(3)*DDTR
      CALL DDMAT(GX2,gx(2))
      GXX(2,1)=TRD(GX1)-GX(1)*DDTR
      GXX(2,2)=TRD(GX2)-GX(2)*DDTR
      GXX(2,3)=TRD(GX3)-GX(3)*DDTR
      CALL DDMAT(GX3,gx(3))
      GXX(3,1)=TRD(GX1)-GX(1)*DDTR
      GXX(3,2)=TRD(GX2)-GX(2)*DDTR
      GXX(3,3)=TRD(GX3)-GX(3)*DDTR
C
      DO 11 L=1,3
      DO 11 M=L,3
      CALL PPCHRM(G,L,M,1)
      AUX1=TR(G)
      GPP(M,L)=(GPP(M,L)+AUX1)/DTR
      IF(L.NE.M)GPP(L,M)=(GPP(L,M)+AUX1)/DTR
   11 CONTINUE
      DO 12 L=1,3
      CALL PCHRM(Y,G,L,2)
      GPX(L,1)=(GPX(L,1)+TR(G))/DTR
      CALL PCHRM(Y,G,L,3)
      GPX(L,2)=(GPX(L,2)+TR(G))/DTR
      CALL PCHRM(Y,G,L,4)
      GPX(L,3)=(GPX(L,3)+TR(G))/DTR
   12 CONTINUE
      CALL CHRM2(Y,G,5)
      GXX(1,1)=(GXX(1,1)+TR(G))/DTR
      CALL CHRM2(Y,G,8)
      GXX(2,2)=(GXX(2,2)+TR(G))/DTR
      CALL CHRM2(Y,G,10)
      GXX(3,3)=(GXX(3,3)+TR(G))/DTR
      CALL CHRM2(Y,G,6)
      AUX1=(GXX(1,2)+TR(G))/DTR
      GXX(1,2)=AUX1
      GXX(2,1)=AUX1
      CALL CHRM2(Y,G,7)
      AUX1=(GXX(1,3)+TR(G))/DTR
      GXX(1,3)=AUX1
      GXX(3,1)=AUX1
      CALL CHRM2(Y,G,9)
      AUX1=(GXX(2,3)+TR(G))/DTR
      GXX(2,3)=AUX1
      GXX(3,2)=AUX1
      CALL CHRM2(Y,G,1)
C
      DO 13 I=1,3
        DERY(I+6)=0.
        DERY(I+9)=0.
        DERY(I+12)=0.
        DERY(I+15)=0.
        DO 13 K=1,3
          DERY(I+6)=DERY(I+6)+.5*(GPX(I,K)*Y(K+6)+GPP(I,K)*Y(K+12))
          DERY(I+9)=DERY(I+9)+.5*(GPX(I,K)*Y(K+9)+GPP(I,K)*Y(K+15))
          DERY(I+12)=DERY(I+12)-.5*(GXX(I,K)*Y(K+6)+GPX(K,I)*Y(K+12))
          DERY(I+15)=DERY(I+15)-.5*(GXX(I,K)*Y(K+9)+GPX(K,I)*Y(K+15))
   13 CONTINUE
      RETURN
      END
C