C
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) 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.000001)then th=pi end if if(a1.gt.0.000001)then th=q*q*q th=p/sqrt(th) th=acos(th) end if a1=-2.*sqrt(q) a2=-r/3. a3=th/3. v(1)=a1*cos(a3)+a2 a3=(th+2.*pi)/3. v(2)=a1*cos(a3)+a2 a3=(th+4.*pi)/3. v(3)=a1*cos(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(4),NINT, 1 XINTA 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) 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.400) 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,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 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) ns=n2-n1 if(ns.eq.0)ns=1 if(ipol.eq.2)ns=1 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,ns 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.ne.0)write(lou,'(a,/,6f10.5)') ' x1,x2,x3,tx,ty,tz', 1ay(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 IF(IPOL.EQ.0.AND.J.NE.N2) GO TO 1100 if(ipol.eq.1.and.(j.ne.n1.and.j.ne.n2))go to 1100 P1=AY(5,J) P2=AY(6,J) P3=AY(7,J) XL=SQRT(P1*P1+P2*P2) if(xl.lt..000001)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 if(ipol.ne.0)write(lou,'(a,/,4f10.5/6f10.5)') 1' 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,ns 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 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.ne.0)write(lou,'(a,/,6f10.5)') ' x1,x2,x3,gx,gy,gz', 1ay(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 C