C
```C
C     *********************************************************
C
SUBROUTINE SPLIN(X,FX,NMIN,NMAX)
C
DIMENSION A(200),B(200),H(200),F(200),X(200),FX(200)
C
C     CUBIC SPLINE INTERPOLATION WITH ZERO CURVATURES AT
C     THE END POINTS
C
IF((NMAX-NMIN).EQ.1)GO TO 4
NMIN1=NMIN+1
NMAX1=NMAX-1
H(NMIN1)=X(NMIN1)-X(NMIN)
D2=(FX(NMIN1)-FX(NMIN))/H(NMIN1)
DO 1 I=NMIN1,NMAX1
H(I+1)=X(I+1)-X(I)
D1=D2
D2=(FX(I+1)-FX(I))/H(I+1)
B(I)=H(I)+H(I+1)
FX(I)=3.*(D2-D1)/B(I)
A(I)=H(I)/B(I)
1 B(I)=H(I+1)/B(I)
4 FX(NMIN)=0.
FX(NMAX)=0.
IF((NMAX-NMIN).EQ.1)RETURN
H(NMIN)=0.
F(NMIN)=0.
DO 2 I=NMIN1,NMAX1
XPOM=2.+A(I)*H(I-1)
H(I)=-B(I)/XPOM
2 F(I)=(FX(I)-A(I)*F(I-1))/XPOM
DO 3 I=NMIN,NMAX1
J=NMAX1-(I-NMIN)
3 FX(J)=H(J)*FX(J+1)+F(J)
RETURN
END
C
C     *********************************************************
C
SUBROUTINE STRESS(P,XN,U,TAU)
C
C  ROUTINE FOR THE COMPUTATION OF NORMAL STRESSES
C
IMPLICIT COMPLEX (P,U,T,E)
REAL N1,N2,N3
DIMENSION XN(3),P(3),U(3),TAU(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,
COMMON/DENS/RHO(20)
C
P1=P(1)
P2=P(2)
P3=P(3)
U1=U(1)
U2=U(2)
U3=U(3)
N1=XN(1)
N2=XN(2)
N3=XN(3)
IF(IANI(LAY).EQ.0) GO TO 1000
C
C  ANISOTROPIC CASE
C
RO=1.7+0.2*SQRT(A11)
IF(IRHO.NE.0) RO=RHO(LAY)
C
C  DISPLACEMENT TENSOR
C
E11=U1*P1
E12=0.5*(U1*P2+U2*P1)
E13=0.5*(U1*P3+U3*P1)
E22=U2*P2
E23=0.5*(U2*P3+U3*P2)
E33=U3*P3
C
C  STRESS TENSOR
C
P11=A11*E11+A12*E22+A13*E33+2.*A14*E23+2.*A15*E13+2.*A16*E12
P22=A12*E11+A22*E22+A23*E33+2.*A24*E23+2.*A25*E13+2.*A26*E12
P33=A13*E11+A23*E22+A33*E33+2.*A34*E23+2.*A35*E13+2.*A36*E12
P23=A14*E11+A24*E22+A34*E33+2.*A44*E23+2.*A45*E13+2.*A46*E12
P13=A15*E11+A25*E22+A35*E33+2.*A45*E23+2.*A55*E13+2.*A56*E12
P12=A16*E11+A26*E22+A36*E33+2.*A46*E23+2.*A56*E13+2.*A66*E12
C
C  NORMAL STRESS
C
TAU(1)=P11*N1+P12*N2+P13*N3
TAU(2)=P12*N1+P22*N2+P23*N3
TAU(3)=P13*N1+P23*N2+P33*N3
TAU(1)=TAU(1)*RO
TAU(2)=TAU(2)*RO
TAU(3)=TAU(3)*RO
RETURN
C
C  ISOTROPIC CASE
C
1000 A12=A11-2.*A44
RO=1.7+0.2*SQRT(A11)
IF(IRHO.NE.0) RO=RHO(LAY)
IF(ICOEF.EQ.0) GOTO 100
WRITE(LOUT,'(A,4F12.7)') 'STRESS: A11,A44,A12,RO',A11,A44,A12,RO
WRITE(LOUT,'(A,/,3(2F12.7,/))') 'STRESS: PSTR',P
WRITE(LOUT,'(A,/,3(2F12.7,/))') 'STRESS: USTR',U
100 TETA=U1*P1+U2*P2+U3*P3
IF(ICOEF.GT.0)
1WRITE(LOUT,'(A,/,5F12.7)') 'STRESS: UN,TETA',XN,TETA
P11=A12*TETA+2.*A44*P1*U1
P12=A44*(P1*U2+P2*U1)
P13=A44*(P1*U3+P3*U1)
P22=A12*TETA+2.*A44*P2*U2
P23=A44*(P2*U3+P3*U2)
P33=A12*TETA+2.*A44*P3*U3
IF(ICOEF.GT.0)
1WRITE(LOUT,'(A,/,6(2F12.7,/))')'STRESS: P11,P12,P13,P22,P23,P33'
1,P11,P12,P13,P22,P23,P33
TAU(1)=P11*N1+P12*N2+P13*N3
TAU(2)=P12*N1+P22*N2+P23*N3
TAU(3)=P13*N1+P23*N2+P33*N3
TAU(1)=TAU(1)*RO
TAU(2)=TAU(2)*RO
TAU(3)=TAU(3)*RO
RETURN
END
C
C     *********************************************************
C
SUBROUTINE SURFPL(LIN,LU3)
C
DIMENSION XX(500),YY(500),ZZ(500)
COMMON /AUXI/  IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT,
1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU,
INTEGER CODE
COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
COMMON /DIST/ XDST(200),NDSTX,XREPS,DST(2),NDST,REPS,LNDST,
1XPRF,YPRF
COMPLEX PS
COMMON /RAY/   AY(28,2000),DS(20,20),KINT(20),HHH(3,3),TMAX,
1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
COMMON /RAY2/  DRY(3,2000)
C
MORI=0
NDST=1
REPS=.001
TMAX=.01
IND=1
NPAR=1
LAY=1
X0=1.
Y0=1.
Z0=1.
DDELTA=.05
AZIM=0.
WRITE(LOU,100)NPAR,LAY
WRITE(LOU,102)X0,Y0,Z0,DDELTA,AZIM
DST(1)=AZIM
DST(2)=AZIM+.1
XPRF=X0
YPRF=Y0
KREF=1
ITYPE=0
CODE(1,1)=LAY
ISOUR=LAY
1 CONTINUE
INUM=0
DELTA=0.
ITYPE=ITYPE+1
CODE(1,2)=ITYPE
IF(NPAR.LE.2)NDSTX=0
IF(NPAR.EQ.3)THEN
NDSTX=1
NDER=2
MDIM=2
END IF
2 CONTINUE
INUM=INUM+1
IF(NPAR.LE.2)THEN
CALL PROFIL(X0,Y0,Z0,0.,DELTA,PAZM,RANG,
1  X,Y,Z,T,.1,.0001,AZIM,1.,AZIM+.1,10,3,1,0,12,0)
XX(INUM)=AY(5,1)
YY(INUM)=AY(6,1)
ZZ(INUM)=AY(7,1)
END IF
IF(NPAR.EQ.3)THEN
CALL PROFIL(X0,Y0,Z0,0.,DELTA,PAZM,RANG,
1  X,Y,Z,T,.1,.0001,AZIM-.5,1.,AZIM+.6,10,3,1,0,12,0)
IF(IND.NE.0)THEN
XX(INUM)=DRY(1,1)
YY(INUM)=DRY(2,1)
ZZ(INUM)=DRY(3,1)
END IF
END IF
DELTA=DELTA+DDELTA
IF(INUM.LT.500.AND.DELTA.LT.6.3)GO TO 2
WRITE(LU3,100)ITYPE,NPAR,INUM
DELTA=0.
DO 3 I=1,INUM
WRITE(LU3,101)DELTA,XX(I),YY(I),ZZ(I)
DELTA=DELTA+DDELTA
3 CONTINUE
IF(ITYPE.LT.3)GO TO 1
ITYPE=0
WRITE(LU3,100)ITYPE
C
100 FORMAT(16I5)
101 FORMAT(4E15.5)
102 FORMAT(8F10.5)
RETURN
END
C
C     *********************************************************
C
SUBROUTINE TRANSL(Y,P,PNEW,UN,ITRANS,IASW)
C
C  ROUTINE FOR THE COMPUTATION OF THE TRANSFORMATION OF THE SLOWNESS
C  VECTOR AT AN INTERFACE
C
SAVE A,AI,B,BI,C,CI,CD,CDI,XCOE,AT,BT,DETB
SAVE XKO,PN,Y1,Z,NDER1,LAY1,ISG,IMAG
DOUBLE PRECISION XCOE,COE
DIMENSION A(3,3),AI(3,3),B(3,3),BI(3,3),C(3,3),CI(3,3),CD(3,3),
*         CDI(3,3),AC(3,3),ACI(3,3),XH1(3),XH2(3),
*         P(3),PNEW(3),UN(3),Y(18),PN(3),Y1(18),DERY(18),XCOE(7),
*         COE(7),XSI(3),ISG(6),IR(3),IT(3),ICR(3),ICT(3),XKO(7)
COMPLEX Z(6),Z0,CSQ,CSI(3),XPR,IMAG
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,
INTEGER CODE
COMMON /COD/  CODE(50,2),KREF,KC,ITYPE
COMMON /DENS/ RHO(20)
COMPLEX PS
COMMON /RAY/   AY(28,2000),DS(20,20),KINT(20),HHH(3,3),tmax,
1               PS(3,7,20),IS(8,20),N,IREF,IND,IND1
COMMON /ZERO/ RNULL
C
C  IASW :      SWITCH INDICATING FROM WHICH PROGRAM TRANSL IS CALLED
C              IASW=0   CALLED FROM AMPL
C              IASW=1   CALLED FROM OUT
C
IMAG=(0.,1.)
IRF=IREF-1
C
C     FILLING DS FIELD
C
IF(IASW.GT.0) THEN
DS(1,IRF)=UN(1)
DS(2,IRF)=UN(2)
DS(3,IRF)=UN(3)
IF(IRHO.EQ.0)DS(8,IRF)=0.2*SQRT(A11)+1.7
KINT(IRF)=N
IS(1,IRF)=LAY
IS(2,IRF)=ITRANS
END IF
CALL PARDIS(Y,0)
IF(ITRANS.EQ.0)IS(7,IRF)=LAY
IF(ITRANS.NE.0)IS(8,IRF)=LAY
DO 1 K=1,3
IT(K)=0
IR(K)=0
ICT(K)=0
ICR(K)=0
ISG(K)=0
ISG(K+3)=0
1  CONTINUE
KDIM=6
IF(NDER.GT.1)KDIM=18
DO 2 K=1,KDIM
2  Y1(K)=Y(K)
PUN=P(1)*UN(1)+P(2)*UN(2)+P(3)*UN(3)
PU1=PUN*UN(1)
PU2=PUN*UN(2)
PU3=PUN*UN(3)
PN(1)=P(1)-PU1
PN(2)=P(2)-PU2
PN(3)=P(3)-PU3
IF(IASW.EQ.0) GO TO 9
NDER1=NDER
NDER=1
LAY1=LAY
LAY=IS(3,IRF)
CALL FCT(X,Y1,DERY)
IF(IND.EQ.10)RETURN
DUN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3)
LAY=LAY1
CALL FCT(X,Y1,DERY)
IF(IND.EQ.10)RETURN
NDER=NDER1
DS(4,IRF)=PN(1)
DS(5,IRF)=PN(2)
DS(6,IRF)=PN(3)
DS(7,IRF)=abs(dun)
9  IF(KC.NE.0) ITYPE=CODE(IREF,2)
IF(IREF.GT.KREF) ITYPE=CODE(KREF,2)
IF(IANI(LAY).EQ.0) GOTO 230
C
C  TRANSMISSION OR REFLECTION FOR THE ANISOTROPIC CASE
C
IF(MTRNS.GT.0) THEN
WRITE(LOUT,999)UN,PN,(Y(K),K=1,KDIM),P
999    FORMAT(' UN,PN,Y,P',/,6(E12.5,2X))
WRITE(LOUT,'(A)')' ITYPE,IREF,LAY,IANI(LAY),N,KINT(IREF-1)'
WRITE(LOUT,'(6I5)')ITYPE,IREF,LAY,IANI(LAY),N,KINT(IREF-1)
WRITE(LOUT,'(A,/,10F10.5,/,11F10.5)')' A(I,J,K,L)=',A11,A12,A13,
1  A14,A15,A16,A22,A23,A24,A25,A26,A33,A34,A35,A36,A44,A45,A46,A55,
2  A56,A66
END IF
CALL CHRM1(A,UN,UN)
AT=A(1,1)+A(2,2)+A(3,3)
CALL INV(A,AI,DET)
XCOE(7)=DET
CALL CHRM1(C,PN,PN)
CALL INV(C,CI,DET)
DO 10 J=1,3
DO 10 K=1,3
CD(J,K)=C(J,K)
IF(J.EQ.K)CD(J,K)=CD(J,K)-1
10 CONTINUE
CALL INV(CD,CDI,DET)
XCOE(1)=DET
CALL CHRM1(B,PN,UN)
CALL INV(B,BI,DET)
BT=B(1,1)+B(2,2)+B(3,3)
DETB=DET
DO 20 J=1,3
DO 20 K=1,3
AC(J,K)=A(J,K)+C(J,K)
20  CONTINUE
CALL INV(AC,ACI,DET)
DO 22 J=1,3
XH1(J)=0.0
DO 21 M=1,3
21  XH1(J)=AI(M,J)*B(J,M)+XH1(J)
22  CONTINUE
XCOE(6)=XH1(1)+XH1(2)+XH1(3)
XCOE(6)=2.*XCOE(6)
DO 30 J=1,3
XH1(J)=0.0
DO 29 M=1,3
29  XH1(J)=AI(M,J)*CD(J,M)+XH1(J)
30  CONTINUE
DO 40 J=1,3
XH2(J)=0.0
DO 39 M=1,3
39  XH2(J)=BI(M,J)*A(J,M)+XH2(J)
40  CONTINUE
XCOE(5)=0.0
DO 60 K=1,3
XCOE(5)=XH1(K)+4.*XH2(K)+XCOE(5)
60  CONTINUE
DO 70 K=1,3
DO 70 J=1,3
ACI(J,K)=ACI(J,K)-AI(J,K)-CI(J,K)
70  CONTINUE
DO 80 J=1,3
XH1(J)=0.0
DO 79 M=1,3
79  XH1(J)=ACI(M,J)*B(J,M)+XH1(J)
80  CONTINUE
DO 90 J=1,3
XH2(J)=0.0
DO 89 M=1,3
89  XH2(J)=A(M,J)*B(J,M)+XH2(J)
90  CONTINUE
XCOE(4)=0.0
DO 100 J=1,3
XCOE(4)=XH1(J)+XH2(J)+XCOE(4)
100  CONTINUE
XCOE(4)=4.*DETB+XCOE(4)-AT*BT
XCOE(4)=2.*XCOE(4)
DO 110 J=1,3
XH1(J)=0.0
DO 109 M=1,3
109  XH1(J)=CDI(M,J)*A(J,M)+XH1(J)
110  CONTINUE
DO 120 J=1,3
XH2(J)=0.0
DO 119 M=1,3
119  XH2(J)=BI(M,J)*CD(J,M)+XH2(J)
120  CONTINUE
XCOE(3)=0.0
DO 130 J=1,3
XCOE(3)=XH1(J)+4.*XH2(J)+XCOE(3)
130  CONTINUE
DO 140 J=1,3
XH1(J)=0.0
DO 139 M=1,3
139  XH1(J)=CDI(M,J)*B(J,M)+XH1(J)
140  CONTINUE
XCOE(2)=XH1(1)+XH1(2)+XH1(3)
XCOE(2)=2.*XCOE(2)
DO 141 K=1,7
M=8-K
141 XKO(M)=XCOE(K)
CALL POLRT(XCOE,COE,6,Z,IER)
iww=0
142 IF(MTRNS.GT.0)then
WRITE(LOUT,1000)XKO
WRITE(LOUT,1010) Z
1000   FORMAT(' COEF. OF POLYNOMIAL',/,4(E12.5,4X))
1010   FORMAT(' ROOTS OF POLYNOMIAL',/,4(E12.5,4X))
end if
IF(IER.EQ.131) THEN
WRITE(LOUT,1030)
1030   FORMAT(/' NOT ALL ZEROS FOUND , PROGRAM TERMINATES'//)
STOP
END IF
IF(IER.EQ.130) THEN
WRITE(LOUT,1040)
1040   FORMAT(/' DEGREE OF POLYNOMAL IS REDUCED %%%%%%%%%%%%%%%%'//)
END IF
DO 150 J=1,6
ZI=AIMAG(Z(J))
IF(ABS(ZI).GT..01)THEN
ISG(J)=0
GOTO 165
else
z(j)=z(j)-zi*IMAG
END IF
C
C CHECK OF REAL-VALUED SLOWNESS VECTORS OF GENERATED WAVES
C
Z0=Z(J)
XPR=XKO(7)+Z0*(XKO(6)+Z0*(XKO(5)+Z0*(XKO(4)+Z0*(XKO(3)+Z0*(XKO(2)+
*Z0*XKO(1))))))
IF(MTRNS.NE.0.OR.(IWW.EQ.1.AND.CABS(XPR).GT..00001))THEN
WRITE(LOUT,1041) XPR,Z0
1041   FORMAT(' PRECISSION OF ROOTS',/,4(E12.5,4X))
END IF
DO 160 K=1,3
Y1(K+3)=PN(K)+REAL(Z(J))*UN(K)
160 CONTINUE
NDER1=NDER
NDER=1
CALL FCT(X,Y1,DERY)
IF(IND.EQ.10)RETURN
NDER=NDER1
AAA=Y1(4)*DERY(1)+Y1(5)*DERY(2)+Y1(6)*DERY(3)
IF(MTRNS.GT.0)WRITE(LOUT,2000)AAA,Y1(4),Y1(5),Y1(6)
2000 FORMAT(1X,'GROUP*SLOWNESS, SLOWNESS'/1X,4E15.5)
IF(ABS(AAA-1.).GT.1.0E-02)THEN
IND=10
RETURN
END IF
SG=UN(1)*DERY(1)+UN(2)*DERY(2)+UN(3)*DERY(3)
IF(MTRNS.GT.0)WRITE(LOUT,1042)SG
1042 FORMAT(' DIRECTION OF RAY VECTOR (UN(I)*DERY(I),I=1,3',F10.5)
IF(MTRNS.GT.0) WRITE(LOUT,1043)(DERY(i),i=1,3)
1043 FORMAT(' RAY VECTOR (DERY)'/3F10.5)
IF(SG.GE.0.)ISG(J)=1
IF(SG.LT.0.)ISG(J)=-1
165 CONTINUE
150 CONTINUE
MCT=0
MCR=0
C
C CHECK OF COMPLEX-VALUED SLOWNESS VECTORS OF GENERATED WAVES
C
DO 179 K=1,6
IF(ISG(K).NE.0) GOTO 179
DIR=AIMAG(Z(K))*UN(3)
IF(UN(3).GT.0.0) GOTO 175
C
C  INCIDENT WAVE STRIKES THE INTERFACE FROM ABOVE
C
IF(DIR.LT.0.0)THEN
MCR=MCR+1
ICR(MCR)=K
END IF
IF(DIR.GT.0.0)THEN
MCT=MCT+1
ICT(MCT)=K
END IF
GOTO 179
C
C  INCIDENT WAVE STRIKES THE INTERFACE FROM BELOW
C
175 IF(DIR.GT.0.0)THEN
MCR=MCR+1
ICR(MCR)=K
END IF
IF(DIR.LT.0.0)THEN
MCT=MCT+1
ICT(MCT)=K
END IF
179 CONTINUE
MT=0
DO 180 K=1,6
IF(ISG(K).EQ.(-1))THEN
MT=MT+1
IT(MT)=K
END IF
180 CONTINUE
MR=0
DO 190 K=1,6
IF(ISG(K).EQ.1)THEN
MR=MR+1
IR(MR)=K
END IF
190 CONTINUE
C
IF(MTRNS.EQ.0) GOTO 191
WRITE(LOUT,1021)ISG
1021 FORMAT(' INDICATIONS,ISG(6),MCR ICR(3),MCT ICT(3),MR IR(3)
1MT IT(3) ',/,6I5)
WRITE(LOUT,1020)MCR,ICR
WRITE(LOUT,1020)MCT,ICT
WRITE(LOUT,1020)MR,IR
WRITE(LOUT,1020)MT,IT
1020  FORMAT(14I5)
191  IS(5,IRF)=MCR
IS(6,IRF)=MCT
IF(ITRANS.EQ.0)GO TO 200
C
C     TRANSMISSION
C
IF(MT.EQ.3) THEN
Z1=REAL(Z(IT(1)))
Z2=REAL(Z(IT(2)))
Z3=REAL(Z(IT(3)))
XSI(3)=AMAX1(Z1,Z2,Z3)
XSI(1)=AMIN1(Z1,Z2,Z3)
XSI(2)=Z1+Z2+Z3-XSI(3)-XSI(1)
END IF
IF(MT.EQ.2) THEN
Z1=REAL(Z(IT(1)))
Z2=REAL(Z(IT(2)))
XSI(2)=AMAX1(Z1,Z2)
XSI(1)=Z1+Z2-XSI(2)
CSI(1)=Z(ICT(1))
END IF
IF(MT.EQ.1) THEN
XSI(1)=REAL(Z(IT(1)))
IF(CABS(Z(ICT(1))).GT.CABS(Z(ICT(2)))) THEN
CSI(1)=Z(ICT(2))
CSI(2)=Z(ICT(1))
ELSE
CSI(1)=Z(ICT(1))
CSI(2)=Z(ICT(2))
END IF
END IF
IF(MT.EQ.0) THEN
CSI(1)=CMPLX(999.,999.)
CSI(3)=CMPLX(0.,0.)
DO 189 K=1,3
IF(CABS(Z(ICT(K))).LT.CABS(CSI(1))) CSI(1)=Z(ICT(K))
IF(CABS(Z(ICT(K))).GT.CABS(CSI(3))) CSI(3)=Z(ICT(K))
189 CONTINUE
CSI(2)=Z(ICT(1))+Z(ICT(2))+Z(ICT(3))-CSI(1)-CSI(3)
END IF
IF(MDIM.LT.1) GOTO 196
C
C     COMPUTATION OF TRANSMITTED SLOWNESS VECTORS FOR EVALUATING
C     AMPLITUDES IN ROUTINE AMPL
C
DO 192 K=1,MT
DO 193 L=1,3
PS(L,K+3,IRF)=PN(L)+XSI(K)*UN(L)
193 CONTINUE
192 CONTINUE
DO 194 K=1,MCT
I=ICT(K)
M=MT+K
DO 195 L=1,3
PS(L,M+3,IRF)=PN(L)+CSI(K)*UN(L)
195 CONTINUE
194 CONTINUE
IF(IASW.EQ.0)RETURN
C
C     OVERCRITICAL INCIDENCE
C
196 IF(ITYPE.GT.MT) IND=9
IF(IND.EQ.9) RETURN
C
GOTO 210
C
C     REFLECTION
C
200 IF(MR.EQ.3) THEN
Z1=REAL(Z(IR(1)))
Z2=REAL(Z(IR(2)))
Z3=REAL(Z(IR(3)))
XSI(3)=AMIN1(Z1,Z2,Z3)
XSI(1)=AMAX1(Z1,Z2,Z3)
XSI(2)=Z1+Z2+Z3-XSI(3)-XSI(1)
END IF
IF(MR.EQ.2) THEN
Z1=REAL(Z(IR(1)))
Z2=REAL(Z(IR(2)))
XSI(2)=AMIN1(Z1,Z2)
XSI(1)=Z1+Z2-XSI(2)
CSI(1)=Z(ICR(1))
END IF
IF(MR.EQ.1) THEN
XSI(1)=REAL(Z(IR(1)))
IF(CABS(Z(ICR(1))).GT.CABS(Z(ICR(2)))) THEN
CSI(1)=Z(ICR(2))
CSI(2)=Z(ICR(1))
ELSE
CSI(1)=Z(ICR(1))
CSI(2)=Z(ICR(2))
END IF
END IF
IF(MR.EQ.0) THEN
CSI(1)=CMPLX(0.0)
CSI(3)=CMPLX(999.,999.)
DO 211 K=1,3
IF(CABS(Z(ICR(K))).GT.CABS(CSI(1))) CSI(3)=Z(ICR(K))
IF(CABS(Z(ICR(K))).LT.CABS(CSI(3))) CSI(1)=Z(ICR(K))
211 CONTINUE
CSI(2)=Z(ICR(1))+Z(ICR(2))+Z(ICR(3))-CSI(1)-CSI(3)
END IF
C
C     COMPUTATION OF REFLECTED SLOWNESS VECTORS FOR EVALUATING
C     AMPLITUDES IN ROUTINE AMPL
C
IF(MDIM.LT.1) GOTO 209
DO 201 K=1,MR
DO 202 L=1,3
PS(L,K,IRF)=PN(L)+XSI(K)*UN(L)
202 CONTINUE
201 CONTINUE
DO 203 K=1,MCR
I=ICR(K)
M=MR+K
DO 204 L=1,3
PS(L,M,IRF)=PN(L)+CSI(K)*UN(L)
204 CONTINUE
203 CONTINUE
IF(IASW.EQ.0) RETURN
C
C     OVERCRITICAL INCIDENCE
C
209  IF(ITYPE.GT.MR) IND=9
IF(IND.EQ.9) RETURN
C
210 IF(MTRNS.GT.0) WRITE(LOUT,1111) XSI,Z1,Z2,Z3,ITYPE
1111  FORMAT(' XSI,Z1,Z2,Z3,ITYPE',/,6E14.6,I10)
XSSI=XSI(ITYPE)
PU1=XSSI*UN(1)
PU2=XSSI*UN(2)
PU3=XSSI*UN(3)
PNEW(1)=PN(1)+PU1
PNEW(2)=PN(2)+PU2
PNEW(3)=PN(3)+PU3
Y(4)=PNEW(1)
Y(5)=PNEW(2)
Y(6)=PNEW(3)
NDER1=NDER
NDER=1
CALL FCT(X,Y,DERY)
IF(IND.EQ.10)RETURN
NDER=NDER1
DUN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3)
RHO2=0.2*SQRT(A11)+1.7
IF(IRHO.NE.0) RHO2=RHO(LAY)
DS(10,IRF)=abs(dun)
DS(11,IRF)=RHO2
IF(MTRNS.GT.0)WRITE(LOUT,1100)PNEW
1100 FORMAT(' PNEW',/,3F12.8)
RETURN
C
C  COMPUTATION OF THE TRANSFORMATION OF THE SLOWNESS VECTOR
C  IN THE ISOTROPIC CASE
C
230 IF(ITYPE.EQ.3)CI2=1./A11
IF(ITYPE.NE.3)CI2=1./A44
PP=0.0
DO 250 K=1,3
PP=PP+P(K)*P(K)
250 CONTINUE
PUN2=PUN*PUN
ROO=CI2-PP+PUN2
ROT1=1./A11-PP+PUN2
ROT2=1./A44-PP+PUN2
IF(MTRNS.GT.0) WRITE(LOUT,1065) ROO,ROT1,ROT2
1065 FORMAT('ROO,ROT1,ROT2',3F10.5)
C
C  ROT1.LT.ROT2 FOR FOR REALISTIC MATERIALS
C
IF(ITRANS.EQ.0.AND.MDIM.GE.1) THEN
IF(ROT1.GT.0.AND.ROT2.GT.0) THEN
MR=3
MCR=0
END IF
IF(ROT1.LE.0.AND.ROT2.GT.0) THEN
MR=2
MCR=1
END IF
IF(ROT1.LE.0.AND.ROT2.LE.0) THEN
MR=0
MCR=3
END IF
IS(5,IRF)=MCR
IF(MTRNS.GT.0) WRITE(LOUT,1066)MR,MCR
1066   FORMAT(' MR,MCR',/,2I5)
DO 251 K=1,MR
IF(K+MCR.EQ.1) RO=ROT1
IF(K+MCR.GT.1) RO=ROT2
SQU=SQRT(RO)
IF(MTRNS.GT.0) WRITE(LOUT,1062) SQU
1062   FORMAT(' SQU',/,2F12.6)
IS(5,IRF)=MCR
J=MR+1-K
DO 252 L=1,3
PS(L,J,IRF)=P(L)-(PUN-SQU)*UN(L)
252   CONTINUE
251   CONTINUE
DO 253 K=1,MCR
M=MR+K
IF(M.EQ.3) CSQ=CSQRT(CMPLX(ROT1,0.0))
IF(M.NE.3) CSQ=CSQRT(CMPLX(ROT2,0.0))
IF(MTRNS.GT.0) WRITE(LOUT,1064)CSQ
1064   FORMAT(' CSQ',2F10.5)
DO 254 L=1,3
PS(L,M,IRF)=P(L)-(PUN-CSQ)*UN(L)
254   CONTINUE
253   CONTINUE
IF(IASW.EQ.0) RETURN
END IF
IF(ITRANS.EQ.1.AND.MDIM.GE.1) THEN
IF(ROT1.GT.0.AND.ROT2.GT.0) THEN
MT=3
MCT=0
END IF
IF(ROT1.LE.0.AND.ROT2.GT.0) THEN
MT=2
MCT=1
END IF
IF(ROT1.LE.0.AND.ROT2.LE.0) THEN
MCT=3
MT=0
END IF
IS(6,IRF)=MCT
IF(MTRNS.GT.0) WRITE(LOUT,1061)MT,MCT
1061   FORMAT(' MT,MCT',/,2I5)
DO 255 K=1,MT
IF(K+MCT.EQ.1) RO=ROT1
IF(K+MCT.NE.1) RO=ROT2
SQU=SQRT(RO)
IF(MTRNS.GT.0) WRITE(LOUT,1062)SQU
J=MT+1-K
DO 256 L=1,3
PS(L,J+3,IRF)=P(L)-(PUN+SQU)*UN(L)
256   CONTINUE
255   CONTINUE
DO 257 K=1,MCT
M=MT+K
IF(M.EQ.3)CSQ=CSQRT(CMPLX(ROT1,0.0))
IF(M.NE.3)CSQ=CSQRT(CMPLX(ROT2,0.0))
IF(MTRNS.GT.0) WRITE(LOUT,1064)CSQ
DO 258 L=1,3
PS(L,M+3,IRF)=P(L)-(PUN+CSQ)*UN(L)
258   CONTINUE
257   CONTINUE
IF(IASW.EQ.0) RETURN
END IF
IF(ROO.LE.0.0) GOTO 300
SQU=SQRT(ROO)
IF(ITRANS.EQ.1) GOTO 280
C
C  REFLECTION
C
PU1=PUN-SQU
DO 260 K=1,3
PNEW(K)=P(K)-PU1*UN(K)
260 CONTINUE
GOTO 275
C
C  TRANSMISSION
C
280 PU1=PUN+SQU
DO 270 K=1,3
PNEW(K)=P(K)-PU1*UN(K)
270 CONTINUE
275 IF(MTRNS.GT.0) WRITE(LOUT,1063)PNEW
1063 FORMAT(' PNEW',/,3F12.6)
PP=0.0
DO 290 K=1,3
PP=PP+PNEW(K)*PNEW(K)
290 CONTINUE
RHO2=0.2*SQRT(A11)+1.7
IF(IRHO.NE.0) RHO2=RHO(LAY)
PUNG=PNEW(1)*UN(1)+PNEW(2)*UN(2)+PNEW(3)*UN(3)
DS(10,IRF)=ABS(PUNG)/PP
DS(11,IRF)=RHO2
RETURN
300 IND=9
NTR=102
RETURN
END
C                                                                 ```