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, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori 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(300),yy(300),zz(300) COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori 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,400),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,400) c dst(1)=0. mori=1 ndst=1 reps=.001 tmax=.01 ind=1 isour=1 read(lin,100)npar read(lin,102)x0,y0,z0,ddelta if(abs(x0).lt..0001.or.abs(y0).lt..0001.or.abs(z0).lt..0001)then x0=1. y0=1. z0=1. end if if(abs(ddelta).lt..000001)ddelta=.05 write(lou,100)npar write(lou,102)x0,y0,z0,ddelta xprf=x0 yprf=y0 kref=1 itype=0 code(1,1)=1 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,.2,.0001,0.,.6,.4,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,.004,.0001,-.5,1.,.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.300.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(4f10.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, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DENS/ RHO(20) COMPLEX PS COMMON /RAY/ AY(28,400),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 IF(ITRANS.EQ.0)IS(5,IRF)=MCR IF(ITRANS.GT.0)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 ROT1.LT.ROT2 FOR PHYSICAL SUFFCIENT ELAST. COEF. IF(ITRANS.EQ.0.AND.MDIM.GE.1) THEN IF(ROT1.GT.0.AND.ROT2.GT.0) MR=3 IF(ROT1.GT.0.AND.ROT2.GT.0) MCR=0 IF(ROT1.LE.0.AND.ROT2.GT.0) MR=2 IF(ROT1.LE.0.AND.ROT2.GT.0) MCR=1 IF(ROT1.LE.0.AND.ROT2.LE.0) MCR=3 IF(ROT1.LE.0.AND.ROT2.LE.0) MR=0 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) MT=3 IF(ROT1.GT.0.AND.ROT2.GT.0) MCT=0 IF(ROT1.LE.0.AND.ROT2.GT.0) MT=2 IF(ROT1.LE.0.AND.ROT2.GT.0) MCT=1 IF(ROT1.LE.0.AND.ROT2.LE.0) MCT=3 IF(ROT1.LE.0.AND.ROT2.LE.0) MT=0 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