SUBROUTINE INIT(PN,V) C C ROUTINE FOR THE CALCULATION OF THE PHASE VELOCITIES FOR THE C SPECIFIED UNIT NORMAL TO THE PHASE FRONT C SOLUTION OF THE EIGENVALUE PROBLEM C DIMENSION C(3,3),PN(3),V(3) DOUBLE PRECISION A1,A2,A3,R,S,T,Q,P,TH 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 C CALL CHRM1(C,PN,PN) C C COMPUTATION OF THE ROOTS OF THE CHARACTERISTIC EQUATION WITH C CARDANS FORMULA C C COMPUTATION OF THE COEFFICIENTS OF THE CHARACTERISTIC EQUATION C X**3+R*X**2+S*X+T=0 C C11=C(1,1) C12=C(1,2) C13=C(1,3) C22=C(2,2) C23=C(2,3) C33=C(3,3) A1=C22*C33-C23*C23 A2=C11*C33-C13*C13 A3=C11*C22-C12*C12 R=-(C11+C22+C33) S=A1+A2+A3 A2=C12*C33-C13*C23 A3=C12*C23-C13*C22 T=-C11*A1+C12*A2-C13*A3 C C SOLUTION OF CUBIC EQUATION FOLLOWING THE ALGORITHM C FROM NUMERICAL RECIPES C Q=(R*R-3.*S)/9. P=(2.*R*R*R-9.*R*S+27.*T)/54. A1=Q*Q*Q-P*P PI=3.14159265 IF(A1.LE.0.00001)THEN TH=PI END IF IF(A1.GT.0.00001)THEN TH=Q*Q*Q TH=P/DSQRT(TH) IF(TH.LE.(-1.))TH=-1. IF(TH.GE.(1.))TH=1. TH=DACOS(TH) END IF A1=-2.*DSQRT(Q) A2=-R/3. A3=TH/3. V(1)=A1*DCOS(A3)+A2 A3=(TH+2.*PI)/3. V(2)=A1*DCOS(A3)+A2 A3=(TH+4.*PI)/3. V(3)=A1*DCOS(A3)+A2 X=V(1) AUX1=X**3+R*X**2+S*X+T X=V(2) AUX2=X**3+R*X**2+S*X+T X=V(3) AUX3=X**3+R*X**2+S*X+T IF(IPRINT.GT.2)WRITE(LOU,'(A/1X,3E15.6)')'INIT: PRECISION OF SOLUT *IONS OF CUBIC EQUATION ', AUX1,AUX2,AUX3 RETURN END C C ********************************************************* C SUBROUTINE INV(A,AINV,DET) C C ROUTINE INV DETERMINESS AN ADJOINT MATRIX OF A REAL 3*3 MATRIX C DIMENSION A(3,3),AINV(3,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 COMMON /ZERO/ RNULL C DET=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))-A(2,1)*(A(1,2)*A(3,3) 1-A(3,2)*A(1,3))+A(3,1)*(A(1,2)*A(2,3)-A(2,2)*A(1,3)) INULL=0 IF(ABS(DET).GT.RNULL) GOTO 100 DO 10 K=1,3 DO 10 J=1,3 IF(ABS(A(J,K)).GT.RNULL) INULL=1 AINV(J,K)=0.0 10 CONTINUE IF(INULL.EQ.0) THEN IF(IPRINT.GT.2)WRITE(LOUT,1000) A ELSE IF(IPRINT.GT.2)WRITE(LOUT,1010) END IF 100 AINV(1,1)=(A(2,2)*A(3,3)-A(3,2)*A(2,3)) AINV(2,1)=-(A(2,1)*A(3,3)-A(2,3)*A(3,1)) AINV(3,1)=(A(2,1)*A(3,2)-A(2,2)*A(3,1)) AINV(1,2)=-(A(1,2)*A(3,3)-A(1,3)*A(3,2)) AINV(2,2)=(A(1,1)*A(3,3)-A(1,3)*A(3,1)) AINV(3,2)=-(A(1,1)*A(3,2)-A(1,2)*A(3,1)) AINV(1,3)=(A(1,2)*A(2,3)-A(1,3)*A(2,2)) AINV(2,3)=-(A(1,1)*A(2,3)-A(1,3)*A(2,1)) AINV(3,3)=(A(1,1)*A(2,2)-A(1,2)*A(2,1)) IF(IPRINT.LE.2)RETURN WRITE(LOUT,1030)A WRITE(LOUT,1031)AINV WRITE(LOUT,1032)DET DO 30 K=1,3 E1=0.0 E2=0.0 E3=0.0 DO 20 J=1,3 E1=E1+A(J,K)*AINV(1,J) E2=E2+A(J,K)*AINV(2,J) E3=E3+A(J,K)*AINV(3,J) 20 CONTINUE WRITE(LOUT,1033)E1,E2,E3 30 CONTINUE C 1000 FORMAT(1x,' INV: ALL ELEMENTS OF MATRIX A ZERO',/,3(3F12.6,/)) 1010 FORMAT(1x,' INV :DETERMINANT OF MATRIX A ZERO') 1030 FORMAT(1X,' INV: A'/3(3F16.8/)) 1031 FORMAT(1X,' INV: AINV'/3(3F16.8/)) 1032 FORMAT(1X,' INV: DET ',F16.8) 1033 FORMAT(1X,' INV: ',3F16.8) RETURN END C C ********************************************************* C SUBROUTINE OUT(X,Y,DERY,IHLF,NDIM,PRMT) C DIMENSION Y(18),DERY(18),UN(3),POLD(3),PNEW(3),YA(3),YB(3),DEP(6), *PRMT(5),DN(3,3),DNG(3,2),XK(3),XKG(3),YK(3),YKG(3),TM(2), *XDAUX(3),YDAUX(3) 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 /DENS/ RHO(20) COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT, 1 XINTA 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) COMMON /ZERO/ RNULL COMMON /VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP COMMON /DYN/XDYN(3,3),ydyn(3,3) common/appr/ xold,xnew,yold(18),dold(18),ynew(18),dnew(18) COMMON /KM/KMAH,SPROLD,TSTOLD C N=N+1 NTR=20 IF(IND.EQ.10)GO TO 25 IF(N.GE.2000) GOTO 152 IF(X.GE.TMAX)GO TO 155 C C TEST OF SATISFACTION OF THE EIKONAL EQUATION AND OF C THE DERIVATIVES OF THE EIKONAL EQUATION C IF(IPREC.EQ.1.AND.NDER.GT.1)THEN TEST0=DERY(1)*Y(4)+DERY(2)*Y(5)+DERY(3)*Y(6) IF(ABS(TEST0).GT.RNULL)THEN DO 80 I=4,6 Y(I)=Y(I)/TEST0 80 CONTINUE END IF V2=1./(Y(4)*Y(4)+Y(5)*Y(5)+Y(6)*Y(6)) TEST3=Y(7)*Y(4)+Y(8)*Y(5)+Y(9)*Y(6) IF(ABS(TEST3).GT.RNULL)THEN DO 81 I=4,6 Y(I+3)=Y(I+3)-V2*TEST3*Y(I) 81 CONTINUE END IF ALP1=DERY(4)*Y(7)+DERY(5)*Y(8)+DERY(6)*Y(9) ALP2=DERY(1)*Y(13)+DERY(2)*Y(14)+DERY(3)*Y(15) IF(ABS(ALP1-ALP2).GT.RNULL)THEN IF(ABS(ALP1).GT..000001)THEN ALP=ALP2/ALP1 DO 82 I=7,9 Y(I)=Y(I)*ALP 82 CONTINUE ELSE DO 86 I=4,6 DERY(I+3)=DERY(I+3)-V2*ALP2*Y(I) 86 CONTINUE END IF END IF TEST4=Y(10)*Y(4)+Y(11)*Y(5)+Y(12)*Y(6) IF(ABS(TEST4).GT.RNULL)THEN DO 83 I=4,6 Y(I+6)=Y(I+6)-V2*TEST4*Y(I) 83 CONTINUE END IF ALP1=DERY(4)*Y(10)+DERY(5)*Y(11)+DERY(6)*Y(12) ALP2=DERY(1)*Y(16)+DERY(2)*Y(17)+DERY(3)*Y(18) IF(ABS(ALP1-ALP2).GT.RNULL)THEN IF(ABS(ALP1).GT..000001)THEN ALP=ALP2/ALP1 DO 84 I=10,12 Y(I)=Y(I)*ALP 84 CONTINUE ELSE DO 87 I=4,6 DERY(I+6)=DERY(I+6)-V2*ALP2*Y(I) 87 CONTINUE END IF END IF END IF IF(IPREC.EQ.2.AND.NDER.GT.1)THEN TEST0=DERY(1)*Y(4)+DERY(2)*Y(5)+DERY(3)*Y(6) TEST1=DERY(1)*Y(13)+DERY(2)*Y(14)+DERY(3)*Y(15) 1 -DERY(4)*Y(7)-DERY(5)*Y(8)-DERY(6)*Y(9) TEST2=DERY(1)*Y(16)+DERY(2)*Y(17)+DERY(3)*Y(18) 1 -DERY(4)*Y(10)-DERY(5)*Y(11)-DERY(6)*Y(12) TEST3=Y(7)*Y(4)+Y(8)*Y(5)+Y(9)*Y(6) TEST4=Y(10)*Y(4)+Y(11)*Y(5)+Y(12)*Y(6) IF(ABS(TEST0-1.).GT.RNULL)WRITE(LOU,100)TEST0 IF(ABS(TEST1).GT.RNULL.OR.ABS(TEST2).GT.RNULL) 1 WRITE(LOU,101)TEST1,TEST2 IF(ABS(TEST3).GT.RNULL.OR.ABS(TEST4).GT.RNULL) 1 WRITE(LOU,102)TEST3,TEST4 100 FORMAT(1X,'P_I*V_I=',E15.5) 101 FORMAT(1X,'DER.EIKON.EQ.=',2E15.5) 102 FORMAT(1X,'P_I*Q_{IJ}=',2E15.5) END IF C C CHECK THE POSITION OF THE RAY WITH RESPECT TO BOUNDARIES OF THE MODEL C NTR=21 IF(N.EQ.1)GO TO 13 IF(IREF.GT.1.AND.(N-1).EQ.KINT(IREF-1)) GOTO 13 NTR=22 IF(Y(1).LT.BRD(1).OR.Y(1).GT.BRD(2).OR.Y(2).LT.BRD(3).OR.Y(2). 1GT.BRD(4)) GOTO 20 C C CHECK WETHER THE RAY CROSSED THE BOREHOLE C IF(ICOD.EQ.0)GO TO 4 IF(ICOD.GT.IREF)GO TO 4 FXYZA=(AY(2,N-1)-XVSP)*XNRM+(AY(3,N-1)-YVSP)*YNRM FXYZB=(Y(1)-XVSP)*XNRM+(Y(2)-YVSP)*YNRM IF(FXYZA*FXYZB.GT.0.)GO TO 4 C C THE RAY CROSSED THE VERTICAL PLANE PERPENDICULAR TO THE LINE C SOURCE-BOREHOLE AND CONTAINING THE BOREHOLE C DETERMINATION OF THE POINT OF INTERSECTION C XNEW=X KDIM=6 IF(NDER.GT.1)KDIM=18 DO 41 I=1,KDIM YNEW(I)=Y(I) DNEW(I)=DERY(I) 41 CONTINUE XA=XOLD XB=XNEW C C DETERMINATION OF THE POINT OF INTERSECTION WITH THE VERTICAL PLANE C PERPENDICULAR TO THE LINE SOURCE-BOREHOLE C CALL ROOT(XA,FXYZA,XB,FXYZB,X,PRMT,-1) IF(IAC.GE.10)GO TO 153 CALL APPROX(X,Y,DERY,KDIM) IND=43 4 CONTINUE C C CHECK WETHER THE RAY CROSSED AN INTERFACE C NTR=23 INTR=LAY+1 CALL DISC(Y,DEP) ZPL=DEP(1) INTR=LAY CALL DISC(Y,DEP) ZPU=DEP(1) NTR=24 IF(ZPL.LE.ZPU) THEN C C TWO NEIGHBOURHOOD INTERFACES CROSS EACH OTHER C WRITE(LOU,'(A)') ' MAUX, LAY, IND,IND1,X,Y, Z,ZU,ZL' WRITE(LOU,'(4I5,5F12.5,/)')MAUX,LAY,IND,IND1,Y(1),Y(2),Y(3),ZPU, 1ZPL GOTO 150 END IF NTR=25 IS(3,IREF)=LAY IF(Y(3).LT.ZPU.OR.Y(3).GT.ZPL) GOTO 30 NTR=27 C C THE RAY DID NOT CROSS AN INTERFACE C C IF(IND.EQ.1.OR.IND.EQ.2)GO TO 148 IF(IND.EQ.30.OR.IND.EQ.31)GO TO 148 13 CONTINUE AY(1,N)=X XOLD=X DO 15 I=1,6 AY(I+1,N)=Y(I) YOLD(I)=Y(I) DOLD(I)=DERY(I) 15 CONTINUE DRY(1,N)=DERY(1) DRY(2,N)=DERY(2) DRY(3,N)=DERY(3) IF(NDER.GT.1)THEN C C DETERMINATION OF KMAH INDEX C IF(N.GT.1)THEN SPR1=Y(8)*Y(12)-Y(9)*Y(11) SPR2=Y(9)*Y(10)-Y(7)*Y(12) SPR3=Y(7)*Y(11)-Y(8)*Y(10) SPR=SPR1*DERY(1)+SPR2*DERY(2)+SPR3*DERY(3) IF(SPR*SPROLD.LT.0..AND.(X-TSTOLD).GT.0.)KMAH=KMAH+1 SPROLD=SPR TSTOLD=X END IF DO 14 I=1,12 YOLD(i+6)=y(i+6) DOLD(I+6)=DERY(I+6) 14 CONTINUE END IF NTR=99 IF(IND.EQ.10)GO TO 25 IF(IND.EQ.43)GO TO 25 C C WRITING IN AY FIELD C CALL PARDIS (Y,1) RETURN C C THE RAY CROSSED ONE OF THE VERTICAL BOUNDARIES OF THE MODEL C 20 IF(Y(1).LT.BRD(1)) IND=1 IF(Y(1).GT.BRD(2)) IND=2 IF(Y(2).LT.BRD(3)) IND=3 IF(Y(2).GT.BRD(4)) IND=4 IF(IND.EQ.3.OR.IND.EQ.4) GOTO 23 AUX=(BRD(IND)-AY(2,N-1))/(Y(1)-AY(2,N-1)) TR=X X=AY(1,N-1)+AUX*(X-AY(1,N-1)) K=6 IF(NDER.GT.1)K=18 DO 21 I=1,K 21 Y(I)=Y(I)+DERY(I)*(X-TR) Y(1)=BRD(IND) AY(1,N)=X DO 3 I=1,6 3 AY(I+1,N)=Y(I) GO TO 4 23 AUX=(BRD(IND)-AY(3,N-1))/(Y(2)-AY(3,N-1)) TR=X X=AY(1,N-1)+AUX*(X-AY(1,N-1)) K=6 IF(NDER.GT.1)K=18 DO 16 I=1,K 16 Y(I)=Y(I)+DERY(I)*(X-TR) Y(2)=BRD(IND) AY(1,N)=X DO 17 I=1,6 17 AY(I+1,N)=Y(I) IF(IND.EQ.3) IND=30 IF(IND.EQ.4) IND=31 GO TO 4 C C THE RAY CROSSED AN INTERFACE C 30 CONTINUE IF(IREF.LE.1.OR.KC.GT.0) GOTO 39 IF(N-KINT(IREF-1).EQ.2) THEN IND=9 IND1=LAY GO TO 25 END IF C C WHICH INTERFACE WAS CROSSED C 39 INTR=LAY IF(Y(3).GE.ZPL) INTR=LAY+1 XNEW=X KDIM=6 IF(NDER.GT.1)KDIM=18 DO 40 I=1,KDIM YNEW(I)=Y(I) DNEW(I)=DERY(I) 40 CONTINUE DO 35 I=1,3 YA(I)=YOLD(I) YB(I)=YNEW(I) 35 CONTINUE CALL DISC(YA,DEP) FXYZA=DEP(1)-YA(3) CALL DISC(YB,DEP) FXYZB=DEP(1)-YB(3) C C DETERMINATION OF THE POINT OF INTERSECTION WITH THE INTERFACE C XA=XOLD XB=XNEW CALL ROOT(XA,FXYZA,XB,FXYZB,XINT,PRMT,1) IF(IAC.GE.10)GO TO 153 NTR=10 C C THE QUANTITIES OF RAY TRACING AND DYNAMIC RAY TRACING C AT THE POINT OF INCIDENCE WILL BE DETERMINED C IND1=INTR CALL APPROX(XINT,Y,DERY,KDIM) AY(1,N)=XINT DO 63 I=1,6 AY(I+1,N)=y(I) 63 CONTINUE CALL PARDIS(Y,1) CALL FCT(X,Y,DERY) DO 64 I=1,3 DRY(I,N)=DERY(I) 64 CONTINUE IF(NDER.GT.1)THEN IF(N.GT.1)THEN SPR1=Y(8)*Y(12)-Y(9)*Y(11) SPR2=Y(9)*Y(10)-Y(7)*Y(12) SPR3=Y(7)*Y(11)-Y(8)*Y(10) SPR=SPR1*DERY(1)+SPR2*DERY(2)+SPR3*DERY(3) IF(SPR*SPROLD.LT.0..AND.(X-TSTOLD).GT.0.)KMAH=KMAH+1 SPROLD=SPR TSTOLD=X END IF DO 62 I=1,3 XDYN(I,1)=Y(I+6) XDYN(I,2)=Y(I+9) XDYN(I,3)=DERY(I) YDYN(I,1)=Y(I+12) YDYN(I,2)=Y(I+15) YDYN(I,3)=DERY(I+3) XK(I)=DERY(I) YK(I)=DERY(I+3) 62 CONTINUE END IF C C DETERMINATION OF THE UNIT NORMAL TO THE INTERFACE AT THE POINT OF C INCIDENCE C CALL DISC (Y,DEP) DFX=DEP(2) DFY=DEP(3) ROO=1+DFX*DFX+DFY*DFY UN3=ROO**(-0.5) UN(3)=-UN3 UN(2)=UN3*DFY UN(1)=UN3*DFX C C SCALAR PRODUCT OF GROUP VELOCITY WITH UNIT NORMAL VECTOR C PN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3) C C ORIENTATION OF THE UNIT NORMAL SO THAT IT POINTS TO THE MEDIUM IN C WHICH THE ENERGY OF INCIDENT WAVE PROPAGATES C IF(PN.LT.0) GOTO 120 UN(1)=-UN(1) UN(2)=-UN(2) UN(3)=-UN(3) UN3=-UN3 120 CONTINUE IF(NDER.GT.1)THEN XIINC=Y(4)*UN(1)+Y(5)*UN(2)+Y(6)*UN(3) UN33=UN3*UN3*UN3 DFXX=DEP(4) DFXY=DEP(5) DFYY=DEP(6) TM(1)=0. TM(2)=0. DO 126 J=1,2 AUX=0. DO 125 I=1,3 TM(J)=TM(J)+UN(I)*XDYN(I,J) AUX=AUX+UN(I)*XK(I) 125 CONTINUE TM(J)=-TM(J)/AUX 126 CONTINUE C C DETERMINATION OF DERIVATIVES OF THE UNIT NORMAL TO INTERFACE C AUX1=1.+DFX*DFX AUX2=1.+DFY*DFY AUX3=DFX*DFY DN(1,1)=(DFXX*AUX2-AUX3*DFXY)*UN33 DN(1,2)=(DFXY*AUX2-AUX3*DFYY)*UN33 DN(1,3)=0. DN(2,1)=(DFXY*AUX1-AUX3*DFXX)*UN33 DN(2,2)=(DFYY*AUX1-AUX3*DFXY)*UN33 DN(2,3)=0. DN(3,1)=(DFX*DFXX+DFY*DFXY)*UN33 DN(3,2)=(DFX*DFXY+DFY*DFYY)*UN33 DN(3,3)=0. DO 129 L=1,3 DO 128 M=1,2 DNG(L,M)=0. DO 127 K=1,3 DNG(L,M)=DNG(L,M)+DN(L,K)*(XDYN(K,M)+XK(K)*TM(M)) 127 CONTINUE 128 CONTINUE 129 CONTINUE END IF NTR=20 IF(KRE.EQ.1) GOTO 24 NTR=30 IF(KRE.EQ.0) GOTO 130 C C MULTIPLE REFLECTED WAVE C IF((IREF+1).GT.KRE.AND.INTR.EQ.INT1) GOTO 151 C C TERMINATION OF THE RAY AT AN INNER INTERFACE OR AT THE FREE SURFACE C OR AT THE BOTTOM OF THE MODEL C IF((IREF+1).GT.KRE) IRR=IREF NTR=35 IF((IREF+1).GT.KRE) GOTO 148 C C NCD : INDICATOR WHETHER A REFLECTION OR TRANSMISSION TAKES PLACE C WITH RESPECT TO THE CODE OF THE WAVE C NCD=CODE(IREF+1,1)-CODE(IREF,1) C C NCD1 : INDICATOR FOR THE TYPE OF THE WAVE IN THE NEW LAYER C NCD1=CODE(IREF+1,2)-CODE(IREF,2) NTR=40 IF(KC.GT.0.AND.IREF.EQ.1.AND.INTR.EQ.LAY) GOTO 151 NTR=50 IF(KC.LT.0.AND.IREF.EQ.1.AND.INTR.NE.LAY) GOTO 151 NTR=60 IF(IREF.EQ.1) GOTO 170 NTR=70 C IF(INTR.EQ.INT1.AND.NCD.NE.0.OR.NCD1.NE.0) GOTO 151 IF(INTR.EQ.INT1.AND.NCD.NE.0)GOTO 151 NTR=75 IF(INTR.EQ.INT1.AND.NCD1.NE.0)GOTO 151 C C REFLECTION OR TRANSMISSION OF THE RAY C NTR=80 IF(INTR.NE.INT1.OR.NCD.NE.0.OR.NCD1.NE.0) GOTO 170 C C COMPOUND ELEMENT OF THE RAY C IREFR=1 KINT(IREF)=0 IREF=IREF+1 NCD=CODE(IREF+1,1)-CODE(IREF,1) NCD1=CODE(IREF+1,2)-CODE(IREF,2) C C TERMINATION AT AN INNER INTERFACE C IRR=IREF NTR=90 IF((IREF+1).GT.KRE) GOTO 24 NTR=100 GOTO 170 C C REFRACTED WAVE C C ORDINARY TERMINATION AT UPPER BOUNDARY C 130 NTR=110 IF(INTR.EQ.LAY.AND.LAY.EQ.1) GOTO 148 NCD=1 NCD1=0 C C SPECIFICATION OF THE LAYER OF THE GENERATED WAVE C 170 INT1=INTR IF(IRHO.NE.0)DS(8,IREF)=RHO(LAY) IRR=IREF IREF=IREF+1 NTR=120 IF(NCD.EQ.0) GOTO 200 NTR=130 IF(INTR.EQ.LAY) GOTO 190 C C TRANSMISSION AT THE LOWER INTERFACE C NTR=140 IF(NCD.LT.0) GOTO 151 C C ORDINARY TERMINATION AT LOWER BOUNDARY C NTR=150 IF(KRE.LE.1.AND.INTR.EQ.NINT) GOTO 148 NTR=160 IF(INTR.EQ.NINT) GOTO 151 LAY=LAY+1 NTR=170 GOTO 200 C C TRANSMISSION AT THE UPPER INTERFACE C 190 NTR=180 IF(NCD.GT.0.AND.KRE.GT.1) GOTO 151 NTR=190 IF(KRE.LE.1.AND.LAY.EQ.1) GOTO 24 NTR=200 IF(LAY.EQ.1) GOTO 151 LAY=LAY-1 C C TRANSFORMATION OF THE SLOWNESS VECTOR C 200 IF(INTR.EQ.NINT.AND.MDIM.NE.0) GOTO 154 DO 210 I=1,3 POLD(I)=Y(3+I) 210 CONTINUE ITRANS=0 C C CHECK WHETHER A REFLECTION OR TRANSMISSION TAKES PLACE C IF(KC.EQ.0) GOTO 217 IF(CODE(IREF,1).EQ.CODE(IREF-1,1)) GOTO 218 217 ITRANS=1 218 CALL TRANSL (Y,POLD,PNEW,UN,ITRANS,1) IF(IND.EQ.9.OR.IND.EQ.10) GOTO 25 DO 220 I=1,3 Y(I+3)=PNEW(I) 220 CONTINUE IF(NDER.GT.1)THEN CALL FCT(X,Y,DERY) DO 225 I=1,3 XKG(I)=DERY(I) YKG(I)=DERY(I+3) 225 CONTINUE XI=PNEW(1)*UN(1)+PNEW(2)*UN(2)+PNEW(3)*UN(3) DXTN=XKG(1)*UN(1)+XKG(2)*UN(2)+XKG(3)*UN(3) DO 228 M=1,2 AUX=XKG(1)*DNG(1,M)+XKG(2)*DNG(2,M)+XKG(3)*DNG(3,M) DO 224 I=1,3 XDAUX(I)=XDYN(I,M) YDAUX(I)=YDYN(I,M) 224 CONTINUE DO 227 I=1,3 AUXX=(DNG(I,M)-UN(I)*AUX/DXTN)*(XI-XIINC) DO 226 K=1,3 AUXX=AUXX-UN(I)*(XKG(K)*(YDaux(K)+YK(K)*TM(M))- 1 YKG(K)*(XDaux(K)+XK(K)*TM(M)))/DXTN 226 CONTINUE XDYN(I,M)=XDYN(I,M)+(XK(I)-XKG(I))*TM(M) YDYN(I,M)=YDYN(I,M)+(YK(I)-YKG(I))*TM(M)+AUXX 227 CONTINUE 228 CONTINUE DO 229 I=1,3 Y(I+6)=XDYN(I,1) Y(I+9)=XDYN(I,2) Y(I+12)=YDYN(I,1) Y(I+15)=YDYN(I,2) 229 CONTINUE END IF PRMT(5)=2. RETURN C C ORDINARY TERMINATION OF THE RAY C 148 IF(IND.EQ.1.OR.IND.EQ.2.OR.IND.EQ.30.OR.IND.EQ.31) GOTO 25 24 IND=INTR+100 IF(KRE.LE.1)IRR=IREF IF(INTR.EQ.1) IND=3 IF(INTR.EQ.NINT) IND=4 KINT(IRR)=N IF(IND.NE.3) GOTO 25 C C COMPUTATION OF REFLECTED WAVES AT THE EARTH SURFACE FOR COEFFICIENTS C OF CONVERSION C IF(MREG.EQ.0.AND.MDIM.NE.0) THEN POLD(1)=Y(4) POLD(2)=Y(5) POLD(3)=Y(6) IREF=IREF+1 CALL TRANSL(Y,POLD,PNEW,UN,0,1) IREF=IREF-1 END IF GOTO 25 150 IND=20 GOTO 25 151 IND=8 GOTO 25 152 IND=13 GOTO 25 153 IND=39 GOTO 25 154 IND=15 GOTO 25 155 IND=12 25 PRMT(5)=3.0 IF(IND.NE.3.AND.IND.NE.43) RETURN ITYPE=CODE(IREF,2) CALL PARDIS(Y,2) RETURN END C C ********************************************************* C SUBROUTINE POLAR (N1,N2,NN,I) C C ROUTINE COMPUTING POLARIZATION VECTORS C C IF(IPOL.EQ.0) POLARIZATION VECTORS ARE NOT PRINTED C IF(IPOL.EQ.1) POLARIZATION VECTORS ARE PRINTED ONLY C FOR THE POINTS OF INCIDENCE OF A RAY AT INTERFACES C IF(IPOL.EQ.2) POLARIZATION VECTORS ARE PRINTED FOR ALL C COMPUTED POINTS ALONG THE RAY C N1,N2 - SUCCESSIVE NUMBERS OF THE FIRST AND THE LAST POINT C OF AN ELEMENT OF A RAY (IN CASE IPOL.EQ.0, N1=N2 AND THIS C NUMBER CORRESPONDS TO THE POINT OF INCIDENCE OF A RAY AT C AN INTERFACE; NEGATIVE N2 - INDICATION OF COMPOSED ELEMENT C OF THE RAY C NN - TOTAL NUMBER OF POINTS ALONG THE RAY C I - NUMBER OF THE ELEMENT OF THE RAY C C CALLED FROM ROUTINES: AMPL C ROUTINES CALLED: CHRM C DIMENSION PXN(3,3),Y(6),B1(3),B2(3),XX(3),E(21),E1(3),E2(3),T(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,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,MORI INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DJK/ D11,D12,D13,D22,D23,D33,DTR COMMON /GAM/ C11,C12,C13,C22,C23,C33 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 EQUIVALENCE(E(1),A11),(E(2),A12),(E(3),A13),(E(4),A14),(E(5),A15) 1 ,(E(6),A16),(E(7),A22),(E(8),A23),(E(9),A24) 2 ,(E(10),A25),(E(11),A26),(E(12),A33),(E(13),A34),(E(14),A35) 2 ,(E(15),A36),(E(16),A44),(E(17),A45),(E(18),A46),(E(19),A55) 2 ,(E(20),A56),(E(21),A66) C IF(N2.LT.0)THEN LAY=IS(3,I-1) N2=-N2 ELSE LAY=IS(3,I) END IF IF(IPOL.GT.0) 1WRITE(LOU,'(A,6I5)')' LAY,I,N1,N2,ITYPE,IANI',LAY,I,N1,N2,ITYPE, 2IANI(LAY) IF(IANI(LAY).NE.0) GO TO 2000 C C ISOTROPIC CASE C IF(NN.EQ.N2) GO TO 1005 IF(ITYPE.NE.3) GO TO 1020 C C P-WAVE POLARIZATION VECTOR C 1005 DO 1010 J=N1,N2 V=1./SQRT(AY(5,J)*AY(5,J)+AY(6,J)*AY(6,J)+AY(7,J)*AY(7,J)) T(1)=V*AY(5,J) T(2)=V*AY(6,J) T(3)=V*AY(7,J) IF((IPOL.EQ.1.AND.J.EQ.N2).OR.IPOL.EQ.2) 1WRITE(LOU,'(A,/,6F10.5)') ' X1,X2,X3,TX,TY,TZ', 2AY(2,J),AY(3,J),AY(4,J),T 1010 CONTINUE DO 1015 K=1,3 HHH(1,K)=0. HHH(2,K)=0. HHH(3,K)=T(K) 1015 CONTINUE IF(N2.EQ.NN) GO TO 1020 RETURN C C S-WAVE POLARIZATION VECTORS C 1020 CONTINUE DO 1100 J=N1,N2 V=SQRT(AY(12,J)) V2=2.*V M1=5 M2=6 M3=7 DVS1=AY(13,J)/V2 DVS2=AY(14,J)/V2 1030 P1=AY(M1,J) P2=AY(M2,J) P3=AY(M3,J) XL=SQRT(P1*P1+P2*P2) IF(XL.LT.RNULL)THEN IF(M1.EQ.4)THEN M1=6 M2=7 M3=5 DVS1=AY(14,J)/V2 DVS2=AY(15,J)/V2 GO TO 1030 END IF IF(M1.EQ.5)THEN M1=7 M2=5 M3=6 DVS1=AY(15,J)/V2 DVS2=AY(13,J)/V2 GO TO 1030 END IF END IF C C COMPUTATION OF ANGLE XI SPECIFYING POLARIZATION VECTORS C IF(IPOL.EQ.2)WRITE(LOU,'(A,/,6F12.6)') 1' P1,P2,P3,XL,D1,D2',P1,P2,P3,XL,DVS1,DVS2 IF(J.EQ.N1)THEN X1=0. TIME=0. XO=(P3*(DVS1*P2-DVS2*P1))/(XL*XL) XI=0. GO TO 1040 END IF TIME=AY(1,J)-AY(1,J-1) X1=(P3*(DVS1*P2-DVS2*P1))/(XL*XL) XI=.5*(X1+XO)*TIME+XI 1040 IF(IPOL.EQ.2)WRITE(LOU,'(A,/,6F12.6)') 1' XI0,XI1,TSTEP,XISTEP',XO,X1,TIME,XI XO=X1 P1=AY(5,J) P2=AY(6,J) P3=AY(7,J) XL=SQRT(P1*P1+P2*P2) IF(XL.LT..0001)THEN B1(1)=1. B1(2)=0. B1(3)=0. B2(1)=0. B2(2)=1. IF(P3.LT.0.)B2(2)=-1. B2(3)=0. ELSE U=SQRT(P1*P1+P2*P2+P3*P3) P13=P1*P3/XL P23=P2*P3/XL B1(1)=P13/U B1(2)=P23/U B1(3)=-xl/U B2(1)=-P2/XL B2(2)=P1/XL B2(3)=0. END IF SN=SIN(XI) CS=COS(XI) DO 1050 K=1,3 E1(K)=CS*B1(K)-SN*B2(K) E2(K)=SN*B1(K)+CS*B2(K) 1050 CONTINUE AY(17,J)=E1(1) AY(18,J)=E1(2) AY(19,J)=E1(3) AY(20,J)=E2(1) AY(21,J)=E2(2) AY(22,J)=E2(3) IF((IPOL.EQ.1.AND.J.EQ.N2).OR.IPOL.EQ.2) 1WRITE(LOU,'(A,/,4F10.5/6F10.5)') 2' TIME,X1,X2,X3/ E1X,E1Y,E1Z,E2X,E2Y,E2Z',(AY(K,J),K=1,4),E1,E2 1100 CONTINUE DO 1110 K=1,3 HHH(1,K)=E1(K) HHH(2,K)=E2(K) IF(N2.EQ.NN) GO TO 1110 HHH(3,K)=0. 1110 CONTINUE RETURN C C ANISOTROPIC CASE C C C COMPUTATION OF POLARIZATION VECTORS C 2000 DO 2300 J=N1,N2 DO 2110 K=1,21 E(K)=AY(K+7,J) 2110 CONTINUE DO 2120 K=1,6 Y(K)=AY(K+1,J) 2120 CONTINUE A2546=A25+A46 A1266=A12+A66 A1355=A13+A55 A1456=A14+A56 A3645=A36+A45 A2344=A23+A44 CALL CHRM(Y) 2125 IF(ABS(DTR).LT..000001)GO TO 2135 PXN(1,1)=D11 PXN(1,2)=D12 PXN(1,3)=D13 XX(1)=D11*D11+D12*D12+D13*D13 PXN(2,1)=D12 PXN(2,2)=D22 PXN(2,3)=D23 XX(2)=D12*D12+D22*D22+D23*D23 PXN(3,1)=D13 PXN(3,2)=D23 PXN(3,3)=D33 XX(3)=D13*D13+D23*D23+D33*D33 XN=0. DO 2130 L=1,3 IF(XN.GE.XX(L))GO TO 2130 XN=XX(L) NX=L 2130 CONTINUE XN=SQRT(XN) IF(XN.GT.RNULL)GO TO 2140 2135 CONTINUE WRITE(LOU,'(A,A,5I5)') 1' XN.LT.RNULL IN POLAR - SHEAR-WAVE SINGULARITY ???', 2' LAY,N,N1,N2',LAY,J,J,N1,N2 IND=10 RETURN 2140 P1=PXN(NX,1)/XN P2=PXN(NX,2)/XN P3=PXN(NX,3)/XN IF(J.NE.N1)THEN XN=P1OLD*P1+P2OLD*P2+P3OLD*P3 IF(XN.LT.0.)THEN P1=-P1 P2=-P2 P3=-P3 ENDIF ENDIF P1OLD=P1 P2OLD=P2 P3OLD=P3 C C CHECK OF PRECISION C Z01=(C11-1.)*P1+C12*P2+C13*P3 Z02=C12*P1+(C22-1.)*P2+C23*P3 Z03=C13*P1+C23*P2+(C33-1.)*P3 Z01=ABS(Z01) Z02=ABS(Z02) Z03=ABS(Z03) IF(Z01.GT.0.01.OR.Z02.GT.0.01.OR.Z03.GT.0.01) THEN WRITE(LOU,'(A,4I5)') ' CHRISTOFFEL EQUATION IS SATISFIED 1 WITH PRECISION LESS THAN 0.01: LAY,ITYPE,RAY ELEMENT,NPOINT ', 2 LAY,ITYPE,I,J END IF IF((IPOL.EQ.1.AND.J.EQ.N2).OR.IPOL.EQ.2) 1WRITE(LOU,'(A,/,6F10.5)') 2' X1,X2,X3,GX,GY,GZ',AY(2,J),AY(3,J),AY(4,J),P1,P2,P3 IF(J.EQ.N2)THEN HHH(ITYPE,1)=P1 HHH(ITYPE,2)=P2 HHH(ITYPE,3)=P3 END IF 2300 CONTINUE RETURN END C C ********************************************************* C SUBROUTINE POLRT(XCOF,COF,M,ZERO,IER) DIMENSION XCOF(7),COF(7) COMPLEX ZERO(6) DOUBLE PRECISION XO,YO,X,Y,XPR,YPR,UX,UY,V,YT,XT,U,XT2,YT2,SUMSQ, / DX,DY,TEMP,ALPHA,xcof,cof C IFIT=0 N=M IER=0 IF(XCOF(N+1)) 10,25,10 10 IF(N) 15,15,32 C 15 IER=1 20 RETURN C 25 IER=4 GO TO 20 C 30 IER=2 GO TO 20 32 IF(N-36) 35,35,30 35 NX=N NXX=N+1 N2=1 KJ1=N+1 DO 40 L=1,KJ1 MT=KJ1-L+1 40 COF(MT)=XCOF(L) C 45 XO=.00500101 YO=0.01000101 C IN=0 50 X=XO C XO=-10.0*YO YO=-10.0*X C X=XO Y=YO IN=IN+1 GO TO 59 55 IFIT=1 XPR=X YPR=Y C 59 ICT=0 60 UX=0.0 UY=0.0 V =0.0 YT=0.0 XT=1.0 U=COF(N+1) IF(U) 65,130,65 65 DO 70 I=1,N L=N-I+1 TEMP=COF(L) XT2=X*XT-Y*YT YT2=X*YT+Y*XT U=U+TEMP*XT2 V=V+TEMP*YT2 FI=I UX=UX+FI*XT*TEMP UY=UY-FI*YT*TEMP XT=XT2 70 YT=YT2 SUMSQ=UX*UX+UY*UY IF(SUMSQ) 75,110,75 75 DX=(V*UY-U*UX)/SUMSQ X=X+DX DY=-(U*UY+V*UX)/SUMSQ Y=Y+DY 78 IF(DABS(DY)+DABS(DX)-1.0D-05) 100,80,80 C 80 ICT=ICT+1 IF(ICT-500) 60,85,85 85 IF(IFIT) 100,90,100 90 IF(N-5) 50,95,95 C 95 IER=3 GO TO 20 100 DO 105 L=1,NXX MT=KJ1-L+1 TEMP=XCOF(MT) XCOF(MT)=COF(L) 105 COF(L)=TEMP ITEMP=N N=NX NX=ITEMP IF(IFIT) 120,55,120 110 IF(IFIT) 115,50,115 115 X=XPR Y=YPR 120 IFIT=0 if(dabs(x)-1.0D-5)121,121,122 121 x=0. 122 IF(DABS(Y)-1.0D-3*DABS(X))135,125,125 125 ALPHA=X+X SUMSQ=X*X+Y*Y N=N-2 GO TO 140 130 X=0.0 NX=NX-1 NXX=NXX-1 135 Y=0.0 SUMSQ=0.0 ALPHA=X N=N-1 140 COF(2)=COF(2)+ALPHA*COF(1) 145 DO 150 L=2,N 150 COF(L+1)=COF(L+1)+ALPHA*COF(L)-SUMSQ*COF(L-1) 155 ZERO(N2)=CMPLX(X,Y) N2=N2+1 IF(SUMSQ) 160,165,160 160 Y=-Y SUMSQ=0.0 GO TO 155 165 IF(N) 20,20,45 END