C
C C ************************************************************* C SUBROUTINE RAY2(X0,Z0,T0,FI0,X,Z,T,FI,TMAX,DT,AC) C C RAY TRACING BY THE RUNGE-KUTTA METHOD C INTEGER CODE DIMENSION Y(4),DERY(4),DIN(4),PRM(5),AUX(8,4),VEL(6) COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20), 1BRD(2),III(30,20),NPNT(20),NINT COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT, 1NTR,IMET COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1 COMMON/COD/CODE(50),KREF,KC EXTERNAL FCT,OUT C C Y(1)=X0 Y(2)=Z0 IREFR=0 KRE=KREF IF(KC.EQ.0)KRE=0 N=0 IREF=1 IOUT=0 isour=ind C C DETERMINATION OF INITIAL CONDITIONS FOR THE RUNGE-KUTTA PROCEDURE C C LAY=ISOUR IF(ISOUR.NE.IABS(CODE(1)))IND=14 IF(ISOUR.NE.IABS(CODE(1)))RETURN ITYPE=ISIGN(1,CODE(1)) Y(3)=0. Y(4)=0. CALL VELOC(Y,VEL) V=VEL(1) Y(3)=COS(FI0)/V Y(4)=SIN(FI0)/V IND=0 IND1=0 PRM(1)=T0 PRM(2)=TMAX PRM(3)=DT PRM(4)=AC T=PRM(1) DO 21 I=1,4 21 DIN(I)=.25 9 DO 10 I=1,4 10 DERY(I)=DIN(I) C C COMPUTATION OF THE RAY C CALL RKGS(PRM,Y,DERY,4,IHLF,FCT,OUT,AUX) IF(IHLF.EQ.11)IND=5 IF(IHLF.EQ.12)IND=6 IF(IHLF.EQ.13)IND=7 IF(IND.GE.5.AND.IND.LE.7)RETURN IF(IND.EQ.20)RETURN IF(IND.EQ.21)RETURN IF(ABS(PRM(5)).GT..0001) GO TO 11 C C INTERUPTION OF THE COMPUTATION OF THE RAY DUE TO (T.GE.PRM(2)) C IF PRM(2).EQ.TMAX - TERMINATION OF THE COMPUTATION C IF PRM(2).NE.TMAX - THE COMPUTATION CONTINUES FURTHER ON WITH C THE REGULAR TIME STEP ALONG THE RAY C T=AY(1,N) IF(T.GE.TMAX)GO TO 15 GO TO 14 11 IF(ABS(PRM(5)-1.1).GT..0001) GO TO 13 C C LOOP FOR THE ITERATIVE DETERMINATION OF THE POINT OF INCIDENCE C (NOT USED IN THIS VERSION) PRM(1)=AY(1,N) PRM(3)=PRM(3)/(2.**(IHLF+1)) DO 12 I=1,4 12 Y(I)=AY(I+1,N) T=AY(1,N) N=N-1 GO TO 9 13 X=Y(1) Z=Y(2) T=AY(1,N) IF(ABS(PRM(5)-2.).GT..0001.AND.IREFR.EQ.1)IND1=-IND1 IF(ABS(PRM(5)-2.).GT..0001)GO TO 20 C C INTEGRATION FROM THE POINT OF REFLECTION/TRANSMISSION C TO THE CLOSEST POINT OF THE RAY CORRESPONDING TO THE REGULAR C TIME STEP ALONG THE RAY C PRM(1)=AY(1,N) I=INT((PRM(1)-T0)/DT) PRM(2)=FLOAT(I+1)*DT+T0 GO TO 9 14 PRM(1)=PRM(2) PRM(2)=TMAX PRM(3)=DT N=N-1 GO TO 9 15 IND=12 X=Y(1) Z=Y(2) IF(IREFR.EQ.1)IND1=-IND1 C 20 V=AY(6,N) CS=Y(3)*V SN=Y(4)*V FI=ATAN2(SN,CS) RETURN END C C ************************************************************* C SUBROUTINE OUT(X,Y,DERY,IHLF,NDIM,PRMT) C C EXTERNAL ROUTINE IN THE RUNGE-KUTTA RAY TRACING. C PERFORMS VARIOUS COMPUTATIONS AT EACH POINT OF THE RAY C INTEGER CODE DIMENSION Y(4),DERY(4),PRMT(5),YOLD(2),YINT(2),VEL(6),Y1(2) COMMON/COD/CODE(50),KREF,KC COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20), 1BRD(2),III(30,20),NPNT(20),NINT COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1 COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT, 1NTR,IMET COMMON/DEN/RHO1(19),RHO2(19),NRHO common/vsp/xvsp,icod,ivsp C C INCD=0 INSPL=0 N=N+1 NRR=N NTR=1 IF(N.GT.400)GO TO 50 C C NORMALIZATION OF SLOWNESS VECTOR C YAN=Y(3) Y(3)=101. CALL VELOC(Y,VEL) Y(3)=YAN IF(N.EQ.1)GO TO 67 VN=VEL(1) VRT=VN*VN*(Y(3)*Y(3)+Y(4)*Y(4)) VRT=1./SQRT(VRT) Y(3)=Y(3)*VRT Y(4)=Y(4)*VRT C C CHECK OF THE POSITION OF THE POINT OF RAY WITH RESPECT TO INTERFAC C INTR=LAY NL=NPNT(INTR)-1 ICH=1 IF(Y(1).LT.BRD(1).OR.Y(1).GT.BRD(2))GO TO 7 65 ICH=0 IND=0 DO 2 I=1,NL J=I IF(Y(1).LT.X1(I+1,INTR)) GOTO 3 2 CONTINUE 3 IF(IREF.EQ.1)GO TO 27 IF(KINT(IREF-1).EQ.(N-1))GO TO 21 27 A2=A1(J,INTR) B2=B1(J,INTR) C2=C1(J,INTR) D2=D1(J,INTR) X2=X1(J,INTR) AUX=Y(1)-X2 ZINT=((D2*AUX+C2)*AUX+B2)*AUX+A2 ZINT1=ZINT IF(Y(2).LT.ZINT) GO TO 9 INTR=LAY+1 NL=NPNT(INTR)-1 DO 4 I=1,NL J=I IF(Y(1).LT.X1(I+1,INTR)) GO TO 5 4 CONTINUE 5 A2=A1(J,INTR) B2=B1(J,INTR) C2=C1(J,INTR) D2=D1(J,INTR) X2=X1(J,INTR) AUX=Y(1)-X2 ZINT=((D2*AUX+C2)*AUX+B2)*AUX+A2 IF(ZINT1.GT.ZINT) GO TO 54 IF(Y(2).GT.ZINT) GO TO 9 21 IF(Y(1).LE.BRD(1).OR.Y(1).GE.BRD(2))GO TO 7 c c CHECK WHETHER THE RAY DID NOT CROSS THE BOREHOLE if(icod.eq.0)go to 67 if(icod.gt.iref)go to 67 aux=(ay(2,n-1)-xvsp)*(y(1)-xvsp) if(aux.gt.0.)go to 67 c c THE RAY CROSSED THE BOREHOLE c 69 ind=22 avert=xvsp go to 33 C C THE RAY DID NOT CROSS ANY INTERFACE C 67 AY(1,N)=X DO 6 I=1,4 6 AY(I+1,N)=Y(I) RETURN C C THE RAY CROSSED ONE OF THE VERTICAL BOUNDARIES C 7 IF(Y(1).LE.BRD(1)) IND=1 IF(Y(1).GE.BRD(2)) IND=2 avert=brd(ind) 33 AUX=Y(1)-AY(2,N-1) IF(ABS(AUX).LT..00001)AUX=0. IF(ABS(AUX).GE..00001)AUX=(avert-AY(2,N-1))/AUX X=AY(1,N-1)+AUX*(X-AY(1,N-1)) Y(1)=avert DO 8 I=2,4 8 Y(I)=AY(I+1,N-1)+AUX*(Y(I)-AY(I+1,N-1)) IF(ICH.EQ.1)GO TO 65 AY(1,N)=X DO 64 I=1,4 64 AY(I+1,N)=Y(I) PRMT(5)=3. DS(2,IREF)=Y(3) DS(3,IREF)=Y(4) KINT(IREF)=N GO TO 49 C C THE RAY CROSSED AN INTERFACE C C THE DETERMINATION OF THE POINT OF INCIDENCE C 9 YOLD(1)=AY(2,N-1) YOLD(2)=AY(3,N-1) IRR=IREF CALL ROOT(YOLD,Y,YINT,ANYX,ANYZ,RC,J,IROOT) c c CHECK WHETHER THE RAY DID NOT CROSS THE BOREHOLE if(icod.eq.0)go to 12 if(icod.gt.iref)go to 12 aux=(yold(1)-xvsp)*(yint(1)-xvsp) if(aux.gt.0.)go to 12 c c THE RAY CROSSED THE BOREHOLE c go to 69 12 ind1=intr IF(IROOT.EQ.100)GO TO 30 C C THE LOOP FOR THE ITERATIVE DETERMINATION OF THE POINT OF INCIDENCE C (NOT USED IN THIS VERSION) IF(IOUT.NE.0)GO TO 11 10 N=N-1 Y1(1)=YINT(1) Y1(2)=YINT(2) PRMT(5)=1.1 IOUT=1 RETURN 11 IAC=IAC+1 IF(IAC.EQ.10)GO TO 30 AUX1=Y1(1)-YINT(1) AUX2=Y1(2)-YINT(2) AUX=AUX1*AUX1+AUX2*AUX2 IF(AUX.GT..00001) GO TO 10 C C DETERMINATION OF PARAMETERS OF INCIDENT WAVE C AT THE POINT OF INCIDENCE C 30 IAC=0 IOUT=0 AY(2,N)=YINT(1) AY(3,N)=YINT(2) XDIF=Y(1)-AY(2,N-1) YDIF=Y(2)-AY(3,N-1) IF(ABS(XDIF).LT.ABS(YDIF))GO TO 60 IF(ABS(XDIF).LT..001)GO TO 61 AUX=(AY(2,N)-AY(2,N-1))/XDIF GO TO 62 60 IF(ABS(YDIF).LT..001)GO TO61 AUX=(AY(3,N)-AY(3,N-1))/YDIF GO TO 62 61 AUX=0. 62 CONTINUE AY(1,N)=AY(1,N-1)+AUX*(X-AY(1,N-1)) X=AY(1,N) AY5=(X-AY(1,N-1))/AY(6,N-1) AY4=AY(4,N-1)-AY(7,N-1)*AY5 AY5=AY(5,N-1)-AY(8,N-1)*AY5 Y(1)=AY(2,N) Y(2)=AY(3,N) C C NORMALIZATION OF SLOWNESS VECTOR AT THE POINT OF INCIDENCE C CALL VELOC(Y,VEL) VN=VEL(1) VRT=AY4*AY4+AY5*AY5 VRT=VN*VN*VRT VRT=1./SQRT(VRT) Y(3)=AY4*VRT Y(4)=AY5*VRT AY(4,N)=Y(3) AY(5,N)=Y(4) C IF(IREF.EQ.1)GO TO 49 NNN=KINT(IREF-1) IF(NNN.LE.0)GO TO 49 IF(ABS(AY(2,NNN)-YINT(1)).GT..0001.OR.INTR.NE.INT1)GO TO 49 IRR=IRR-1 NTR=2 IF(AINDEX.GT.1)GO TO 28 NTR=3 GO TO 26 C C 49 ITYPE1=ITYPE ITYPE=1 CALL VELOC(Y,VEL) DS(4,IREF)=VEL(1) DS(8,IREF)=RHO1(LAY)+RHO2(LAY)*VEL(1) ITYPE=-1 CALL VELOC(Y,VEL) DS(5,IREF)=VEL(1) ITYPE=ITYPE1 VVV=DS(4,IREF) IF(ITYPE.LT.0)VVV=DS(5,IREF) DS(1,IREF)=0. DS(6,IREF)=0. DS(7,IREF)=0. DS(9,IREF)=0. DS(10,IREF)=0. DS(11,IREF)=0. NTR=4 IF(IND.EQ.1.OR.IND.EQ.2.or.ind.eq.22)GO TO 32 LAY1=LAY PX=Y(3) PZ=Y(4) YAN=Y(3) AUX=ANYX*PX+ANYZ*PZ IF(AUX.LE.0)GO TO 13 ANYX=-ANYX ANYZ=-ANYZ AUX=-AUX RC=-RC 13 AUXX=(ANYX*PZ-ANYZ*PX)*VVV AUXZ=AUX*VVV DS(1,IREF)=RC DS(2,IREF)=AUXX DS(3,IREF)=AUXZ NTR=5 IF(KRE.EQ.1)GO TO 24 IF(KRE.LT.1)GO TO 31 C C MULTIPLY REFLECTED WAVE C NTR=6 IF((IREF+1).GT.KRE.AND.INTR.EQ.INT1)GO TO 26 IF((IREF+1).GT.KRE) GO TO 24 NCD=CODE(IREF+1)-CODE(IREF) NTR=7 IF(KC.GT.0.AND.IREF.EQ.1.AND.INTR.EQ.LAY)GO TO 26 NTR=8 IF(KC.LT.0.AND.IREF.EQ.1.AND.INTR.NE.LAY)GO TO 26 IF(IREF.EQ.1)GO TO 46 NTR=9 IF(INTR.EQ.INT1.AND.NCD.NE.0)GO TO 26 IF(INTR.NE.INT1.OR.NCD.NE.0)GO TO 46 IREFR=1 KINT(IREF)=0 IREF=IREF+1 IRR=IREF DS(1,IREF)=DS(1,IREF-1) DS(2,IREF)=DS(2,IREF-1) DS(3,IREF)=DS(3,IREF-1) DS(4,IREF)=DS(4,IREF-1) DS(5,IREF)=DS(5,IREF-1) DS(6,IREF)=0. DS(7,IREF)=0. DS(8,IREF)=DS(8,IREF-1) DS(9,IREF)=0. DS(10,IREF)=0. DS(11,IREF)=0. DO 53 I=1,9 53 DS(I,IREF-1)=0. NCD=CODE(IREF+1)-CODE(IREF) NTR=10 IF((IREF+1).GT.KRE) GO TO 24 46 INT1=INTR J11=J I11=INTR L11=LAY1 IF(NCD.NE.0) GO TO 22 C C REFLECTION OF UNCONVERTED WAVE C IIII=III(J,INTR) NTR=11 IF(IIII.EQ.(-2))GO TO 48 NTR=12 IF(IIII.GT.0)GO TO 23 Y(3)=101. CALL VELOC(Y,VEL) Y(4)=PZ-2.*AUX*ANYZ Y(3)=PX-2.*AUX*ANYX IREF=IREF+1 V1=VEL(1) V2=V1 AINDEX=1. GO TO 20 C C REFRACTED WAVE C 31 NTR=13 IF(INTR.EQ.LAY.AND.LAY.EQ.1) GO TO 24 NTR=14 IF(INTR.GT.LAY.AND.INTR.EQ.NINT) GO TO 24 NCD=1 INT1=INTR GO TO 14 C C REFRACTION OR REFLECTION OF CONVERTED WAVE C 22 NCD=IABS(CODE(IREF+1))-IABS(CODE(IREF)) IIII=III(J,INTR) NTR=15 IF(NCD.EQ.0.AND.IIII.EQ.(-2)) GO TO 48 NTR=16 IF(NCD.EQ.0.AND.IIII.GT.0)GO TO 23 14 Y(3)=101. CALL VELOC(Y,VEL) Y(3)=YAN V1=VEL(1) IF(KRE.GT.1)ITYPE=ISIGN(1,CODE(IREF)) 45 IREF=IREF+1 C C CHECK WHETHER A TRANSMISSION TAKES PLACE AT AN INTERFACE C WHICH COINCIDES WITH ANOTHER INTERFACE C IF(NCD.EQ.0.AND.INSPL.EQ.0)GO TO 16 IF(NCD.EQ.0.AND.INSPL.EQ.1)GO TO 68 IF(INTR.EQ.LAY) GO TO 15 NTR=17 IF(NCD.LT.0)GO TO 26 NTR=18 IF(KRE.LE.1.AND.INTR.EQ.NINT)GO TO 24 NTR=19 IF(INTR.EQ.NINT) GO TO 26 LAY=LAY+1 GO TO 29 68 INTR=INTRA J=JA INCD=1 GO TO 37 15 NTR=20 IF(NCD.GT.0.AND.KRE.GT.1)GO TO 26 NTR=21 IF(KRE.LE.1.AND.LAY.EQ.1)GO TO 24 NTR=22 IF(LAY.EQ.1)GO TO 26 LAY=LAY-1 29 INTRA=INTR JA=J IF(INTR.EQ.LAY1)GO TO 36 NC=NPNT(INTR+1) DO 34 I=1,NC JJ1=I II1=III(I,INTR+1) IF(J.EQ.II1)GO TO 35 34 CONTINUE GO TO 16 35 INTR=INTR+1 J=JJ1 GO TO 37 36 II1=III(J,INTR) IF(II1.LE.0)GO TO 16 INTR=INTR-1 J=II1 C C TRANSMISSION AT AN INTERFACE WHICH COINCIDES C WITH ANOTHER INTERFACE C 37 N=N+1 N1=N+1 NTR=23 IF(N1.Ge.400) GO TO 50 INSPL=1 NDR=8 IF(NDER.EQ.1)NDR=11 DO 38 I=N,N1 DO 38 L=1,NDR 38 AY(L,I)=AY(L,N-1) N=N1 LAY1=LAY IND1=INTR KINT(IREF)=-1 DO 47 I=1,11 47 DS(I,IREF)=0. NTR=24 IF(KRE.GT.1.AND.(CODE(IREF)*ITYPE).LT.0)GO TO 66 NTR=25 IF(KRE.GT.1.AND.(IREF+1).GT.KRE)GO TO 24 INT1=INTR IF(KRE.GT.1)NCD=IABS(CODE(IREF+1))-IABS(CODE(IREF)) NTR=26 IF(NCD.EQ.0.AND.III(J,INTR).GT.0)GO TO 23 GO TO 45 C C DETERMINATION OF THE ANGLE OF THE RAY OF GENERATED WAVE C AT THE POINT OF INCIDENCE C 16 LAY2=LAY IF(KRE.GT.1)ITYPE=ISIGN(1,CODE(IREF)) CALL VELOC(Y,VEL) V2=VEL(1) IF(INCD.EQ.1)NCD=0 INCD=0 INSPL=0 AINDEX=V2/V1 AUX1=1./(V2*V2)-1./(V1*V1)+AUX*AUX IF(AUX1.GT.0.) GO TO 17 NTR=27 GO TO 28 17 IF(NCD.EQ.0)AUXZ=-SQRT(AUX1) IF(NCD.NE.0)AUXZ=SQRT(AUX1) AUX1=AUXZ+AUX AUXZ=AUXZ*V2 Y(3)=PX-ANYX*AUX1 Y(4)=PZ-ANYZ*AUX1 C C C DETERMINATION OF PARAMETERS OF GENERATED WAVE C AT THE POINT OF INCIDENCE C 20 PRMT(5)=2. KINT(IRR)=NRR DS(1,IRR)=RC DS(10,IRR)=AINDEX*AUXX DS(11,IRR)=-AUXZ IF(NCD.NE.0)GO TO 52 C C DETERMINATION OF VELOCITIES AND DENSITIES C ON THE OPPOSITE SIDE OF A REFLECTOR C JA=J11 INTRA=I11 LAY2=L11 IF(INTRA.EQ.L11)GO TO 43 40 LAY2=LAY2+1 NTR=28 IF(INTRA.EQ.NINT.AND.NDER.NE.0)GO TO 51 IF(INTRA.EQ.NINT)RETURN NC=NPNT(INTRA+1) DO 41 I=1,NC JJ1=I II1=III(I,INTRA+1) IF(JA.EQ.II1)GO TO 42 41 CONTINUE GO TO 44 42 INTRA=INTRA+1 JA=JJ1 GO TO 40 43 LAY2=LAY2-1 II1=III(JA,INTRA) IF(II1.LE.0)GO TO 44 INTRA=INTRA-1 JA=II1 GO TO 43 44 IF(LAY2.NE.0)GO TO 52 DS(6,IRR)=0. DS(7,IRR)=0. DS(9,IRR)=0. RETURN 52 ITYPE1=ITYPE ITYPE=1 L11=LAY LAY=LAY2 CALL VELOC(Y,VEL) DS(6,IRR)=VEL(1) DS(9,IRR)=RHO1(LAY2)+RHO2(LAY2)*VEL(1) ITYPE=-1 CALL VELOC(Y,VEL) DS(7,IRR)=VEL(1) ITYPE=ITYPE1 LAY=L11 RETURN C C CHECK WHETHER THE RAY DOES NOT TERMINATE AT A FICTITIOUS INTERFACE C 24 IND=INTR+100 IF(IREF.EQ.1.AND.KC.GT.0.AND.INTR.EQ.INT1)GO TO 26 IF(IREF.EQ.1.AND.KC.LT.0.AND.INTR.NE.INT1)GO TO 26 IF(INTR.EQ.1)IND=3 IF(INTR.EQ.NINT)IND=4 NTR=29 IF(III(J,INTR).EQ.(-2))GO TO 48 IF(III(J,INTR).LE.0)GO TO 39 IEL=J IDS=INTR 63 II1=III(IEL,IDS) NTR=30 IF(III(II1,IDS-1).EQ.(-2)) GO TO 48 IF(III(II1,IDS-1).LE.0)GO TO 39 IEL=III(II1,IDS-1) IDS=IDS-1 GO TO 63 C C TERMINATION OF COMPUTATION C 48 IND=16 GO TO 25 51 IND=15 GO TO 39 23 IND=17 GO TO 39 66 IND=18 GO TO 39 26 IND=8 GO TO 39 54 IND=20 GO TO 39 28 IND=9 39 IREF=IRR 25 PRMT(5)=3. KINT(IRR)=NRR GO TO 32 50 IND=13 PRMT(5)=3. 32 IF(IND.EQ.9)RETURN YAN=Y(3) Y(3)=101. CALL VELOC(Y,VEL) Y(3)=YAN RETURN END C C ***************************************************************** C SUBROUTINE FCT(X,Y,DERY) C C COMPUTATION OF THE RIGHT-HAND SIDES C OF THE RAY TRACING SYSTEM EQUATIONS C DIMENSION VEL(6),Y(4),DERY(4) COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT, 1NTR,IMET C CALL VELOC(Y,VEL) VV=VEL(1) V1=1./VV VV=VV*VV DERY(1)=VV*Y(3) DERY(2)=VV*Y(4) DERY(3)=-V1*VEL(2) DERY(4)=-V1*VEL(3) RETURN END C C **************************************************************** C SUBROUTINE ROOT(YOLD,YNEW,YINT,ANYX,ANYZ,RC,J,IROOT) C C DETERMINATION OF THE POINT OF INCIDENCE AT AN INTERFACE C (INTERSECTION OF THE RAY WITH AN INTERFACE) C C DIMENSION YDIF(3),YOLD(2),YNEW(2),YINT(2) COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20), 1BRD(2),III(30,20),NPNT(20),NINT COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT, 1NTR,IMET C SH(X9)=(EXP(X9)-EXP(-X9))/2. CH(X9)=(EXP(X9)+EXP(-X9))/2. ARCCOS(X9)=ATAN(SQRT(1.-X9*X9)/X9) ARCCH(X9)=ALOG(X9+SQRT(X9*X9-1.)) ARCSH(X9)=ALOG(X9+SQRT(X9*X9+1.)) C C IROOT=100 C C DETERMINATION OF THE ELEMENT OF THE INTERFACE CONTAINING C THE POINT OF INCIDENCE C LRT=0 30 NC=NPNT(INTR) XA=YOLD(1) XB=YNEW(1) AUX=XB-XA AUXAB=AUX X2=X1(J,INTR) A2=A1(J,INTR) B2=B1(J,INTR) C2=C1(J,INTR) D2=D1(J,INTR) AUX1=XB-X2 IF(ABS(AUX).GT..00001)GO TO 1 YINT(1)=YOLD(1) GO TO 2 1 PP=(YNEW(2)-YOLD(2))/AUX QQ=YOLD(2)-PP*XA IF(AUX.LT.0.)GO TO 21 IF(XA.GE.X2)GO TO 2 X3=XB A3=((D2*AUX1+C2)*AUX1+B2)*AUX1+A2 Y2=PP*X2+QQ Y3=PP*X3+QQ 20 AUX1=(Y3-A3)*(Y2-A2) IF(AUX1.LE.0.)GO TO 2 J=J-1 X3=X2 X2=X1(J,INTR) A3=A2 A2=A1(J,INTR) IF(XA.GE.X2)GO TO 2 Y3=Y2 Y2=PP*X2+QQ GO TO 20 21 X3=X1(J+1,INTR) IF(XA.LE.X3)GO TO 2 A3=A1(J+1,INTR) X2=XB A2=((D2*AUX1+C2)*AUX1+B2)*AUX1+A2 Y2=PP*X2+QQ Y3=PP*X3+QQ 22 AUX1=(Y3-A3)*(Y2-A2) IF(AUX1.LE.0.)GO TO 2 J=J+1 X2=X3 X3=X1(J+1,INTR) A2=A3 IF(XA.LE.X3)GO TO 2 A3=A1(J+1,INTR) Y2=Y3 Y3=PP*X3+QQ GO TO 22 C 2 X2=X1(J,INTR) A2=A1(J,INTR) B2=B1(J,INTR) C2=C1(J,INTR) D2=D1(J,INTR) IF(ABS(AUX).LE..00001)GO TO 17 XX=XA-X2 C C DETERMINATION OF COEFFICIENTS OF CUBIC EQUATION C D*X**3+C*X**2+B*X+A=0 C ROOTS ARE LOOKED FOR IN INTERVAL (0,1) C D=D2*AUX*AUX*AUX C=(3*D2*XX+C2)*AUX*AUX B=((3.*D2*XX+2.*C2)*XX+B2-PP)*AUX A=((D2*XX+C2)*XX+B2)*XX+A2-PP*XA-QQ C IF(ABS(D).LT..000001)GO TO 10 C C TRANSFORMATION OF CUBIC EQUATION INTO FORM Y**3+3*P*Y+Q=0 C SUBSTITUTING Y=X+C/(3*D) C AUX1=C/(3.*D) Q=AUX1*AUX1*AUX1-.5*(B*AUX1-A)/D P=B/(3.*D)-AUX1*AUX1 DISKR=Q*Q+P*P*P IF(Q.EQ.0.) GO TO 8 IF(P.EQ.0.) GO TO 7 R=SIGN(1.,Q)*SQRT(ABS(P)) AX=Q/(R*R*R) IF(P.GT.0.) GO TO 6 IF(DISKR)3,3,5 C C P.LT.0.AND.DISKR.LE.0 C 3 D=ARCCOS(AX)/3. XR=-2.*R*COS(D)-AUX1 XR1=2.*R*COS(1.047198-D)-AUX1 XR2=2.*R*COS(1.047198+D)-AUX1 25 NRT=0 IF(XR.LT.0..OR.XR.GT.1.)GO TO 18 RR=XA+XR*AUX-YOLD(1) NRT=NRT+1 YDIF(NRT)=RR 18 IF(XR1.LT.0..OR.XR1.GT.1.)GO TO 19 RR=XA+XR1*AUX-YOLD(1) NRT=NRT+1 YDIF(NRT)=RR 19 IF(XR2.LT.0..OR.XR2.GT.1.)GO TO 4 RR=XA+XR2*AUX-YOLD(1) NRT=NRT+1 YDIF(NRT)=RR 4 IF(NRT.EQ.1)GO TO 26 IF(NRT.NE.0)GO TO 40 IF(ABS(XR).LE.ABS(XR1))GO TO 41 XR=XR1 41 IF(ABS(XR).LE.ABS(XR2))GO TO 42 XR=XR2 42 YDIF(1)=XA+XR*AUX-YOLD(1) GO TO 26 40 Y1=YDIF(1) Y2=YDIF(2) IF(ABS(Y1).LT..00001.AND.INTR.EQ.INT1)Y1=1000000. IF(ABS(Y2).LT..00001.AND.INTR.EQ.INT1)Y2=1000000. IF(ABS(Y2).LT.ABS(Y1))YDIF(1)=Y2 IF(ABS(Y2).LT.ABS(Y1))Y1=Y2 IF(NRT.EQ.2)GO TO 26 Y3=YDIF(3) IF(ABS(Y3).LT..00001.AND.INTR.EQ.INT1)Y3=1000000. IF(ABS(Y3).LT.ABS(Y1))YDIF(1)=Y3 26 YINT(1)=YDIF(1)+YOLD(1) GO TO 17 C C P.LT.0..AND.DISKR.GT.0. C 5 XR=-2.*R*CH(ARCCH(AX)/3.)-AUX1 GO TO 15 C C P.GT.0 C 6 XR=-2.*R*SH(ARCSH(AX)/3.)-AUX1 GO TO 15 C C P.EQ.0. C 7 XR=-SIGN(1.,Q)*EXP(ALOG(2.*ABS(Q))/3.)-AUX1 GO TO 15 C C Q.EQ.0 C 8 XR=-AUX1 IF(P)9,15,15 C C Q.EQ.0..AND.P.LT.0. C 9 XR1=SQRT(-3.*P)-AUX1 XR2=-XR1-2.*AUX1 GO TO 25 C C REDUCTION OF CUBIC EQUATION TO QUADRATIC EQUATION C 10 IF(ABS(C).LT..000001) GO TO 14 DISKR=B*B-4.*C*A P=-B/(2.*C) IF(DISKR)11,11,12 C 11 XR=P GO TO 15 C 12 Q=SQRT(DISKR)/(2.*C) XR=P+Q XR1=P-Q 13 NRT=2 IF(XR.GE.0..AND.XR.LE.1.)GO TO 23 NRT=NRT-1 XR=XR1 23 IF(XR1.GE.0..AND.XR1.LE.1.)GO TO 24 NRT=NRT-1 24 IF(NRT.EQ.1)GO TO 15 RR=XA+XR*AUX-YOLD(1) IF(ABS(RR).LT..00001.AND.INTR.EQ.INT1)RR=1000000. RR1=XA+XR1*AUX-YOLD(1) IF(ABS(RR1).LT..00001.AND.INTR.EQ.INT1)RR1=1000000. IF(ABS(RR1).LT.ABS(RR))XR=XR1 GO TO 15 C C REDUCTION OF CUBIC EQUATION TO LINEAR EQUATION C 14 XR=-A/B 15 YINT(1)=XA+XR*(XB-XA) GO TO 17 27 J=J+1 GO TO 30 28 J=J-1 GO TO 30 17 IF(YINT(1).GE.YOLD(1).AND.YINT(1).LE.YNEW(1))GO TO 16 IF(YINT(1).GE.YNEW(1).AND.YINT(1).LE.YOLD(1))GO TO 16 YINT(1)=.5*(YOLD(1)+YNEW(1)) 16 AUX=YINT(1)-X2 YINT(2)=((D2*AUX+C2)*AUX+B2)*AUX+A2 AUX1=((3.*D2*AUX+2.*C2)*AUX+B2) AUX2=1./(AUX1*AUX1+1.) AUX2=SQRT(AUX2) ANYX=-AUX1*AUX2 ANYZ=AUX2 IF(LRT.NE.1)GO TO 31 C C CHECK WHETHER THE RAY DID NOT CROSS THE INTERFACE TWICE C XAA=YINT(1)-YOLD(1) XBB=YINT(2)-YOLD(2) IF(ABS(XBB).LT..0001.OR.ABS(XAA).LT..0001)GO TO 35 XAA=YNEW(1)-YOLD(1) XBB=YNEW(2)-YOLD(2) AUX=XAA*ANYX+XBB*ANYZ IF(INTR.EQ.LAY.AND.AUX.LT.0.)GO TO 31 IF(INTR.NE.LAY.AND.AUX.GT.0.)GO TO 31 IF(YINT(1).GT.YOLD(1))GO TO 28 IF(YINT(1).LT.YOLD(1))GO TO 27 31 IF(INTR.EQ.LAY)GO TO 32 IF(III(J,INTR).LE.0)GO TO 35 LRT=1 J=III(J,INTR) INTR=INTR-1 GO TO 30 32 NC1=NPNT(INTR+1)-1 DO 33 I=1,NC1 I1=I IF(J.EQ.III(I,INTR+1))GO TO 34 33 CONTINUE GO TO 35 34 J=I1 LRT=1 INTR=INTR+1 GO TO 30 35 RC=(6.*D2*AUX+2.*C2) RC=RC*ANYZ*ANYZ*ANYZ RETURN END C C ************************************************************* C SUBROUTINE AMPL(FJ,AX,AAY,AZ,PHX,PHY,PHZ) C C COMPUTATION OF AMPLITUDES C INTEGER CODE COMMON/COD/CODE(50),KREF,KC COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1 COMMON/SOUR/ROS,VPS,VSS COMMON/SH/KSH,NSH C PH=0. PHSH=0. Q=1. QSH=1. V=1. AL=1. AX=0. AAY=0. AZ=0. PHX=0. PHY=0. PHZ=0. IREF1=IREF-1 IF(IREF1.EQ.0)IC2=CODE(1) IF(IREF1.EQ.0)GO TO 10 C C LOOP FOR INTERFACES C DO 5 I=1,IREF1 I1=KINT(I) IF(I1.LE.0)GO TO 5 V1=AY(6,I1) PP=DS(2,I)/V1 P=ABS(PP) VP1=DS(4,I) VS1=DS(5,I) RO1=DS(8,I) VP2=DS(6,I) VS2=DS(7,I) RO2=DS(9,I) IF(KREF.LE.1)IC1=CODE(1) IF(KREF.LE.1)IC2=CODE(1) IF(KREF.LE.1)GO TO 11 IC1=CODE(I) II=I 12 II=II+1 IF(KINT(II).LE.0)GO TO 12 IC2=CODE(II) IF((IABS(IC2)-IABS(IC1)).EQ.0)GO TO 2 11 IF(IC1.LT.0)GO TO 1 AV=RO1*VP1 IF(IC2.GT.0)NC=3 IF(IC2.GT.0)AV=(RO2*VP2)/AV IF(IC2.LT.0)NC=4 IF(IC2.LT.0)AV=(RO2*VS2)/AV GO TO 4 1 AV=RO1*VS1 IF(IC2.GT.0)NC=7 IF(IC2.GT.0)AV=(RO2*VP2)/AV IF(IC2.LT.0)NC=8 IF(NSH.NE.0)NCY=10 IF(IC2.LT.0)AV=(RO2*VS2)/AV GO TO 4 2 IF(ABS(VP2).LT..00001)RO2=0. IF(IC1.LT.0)GO TO 3 IF(IC2.GT.0)NC=1 IF(IC2.GT.0)AV=1. IF(IC2.LT.0)NC=2 IF(IC2.LT.0)AV=VS1/VP1 GO TO 4 3 IF(IC2.GT.0)NC=5 IF(IC2.GT.0)AV=VP1/VS1 IF(IC2.LT.0)NC=6 IF(NSH.NE.0)NCY=9 IF(IC2.LT.0)AV=1. IF(ABS(RO2).LT..00001)NC=NC-2 IF(ABS(RO2).LT..00001.AND.NSH.NE.0)NCY=9 GO TO 4 4 V=V*AV ND=0 IF(PP.LT.0.)ND=1 IF(KSH.EQ.0)GO TO 15 CALL COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,NCY,ND,RY,PHSY) QSH=QSH*RY PHSH=PHSH+PHSY IF(NSH.EQ.1)GO TO 16 15 CALL COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,NC,ND,R,PHS) Q=Q*R PH=PH+PHS 16 SN1=DS(3,I) SN2=DS(11,I) AL=AL*ABS(SN1/SN2) 5 CONTINUE C C END OF LOOP FOR INTERFACES C 10 V0=AY(6,1) V=V*ROS*V0 RO1=DS(8,IREF) I=KINT(IREF) VF=AY(6,I) V=V/(RO1*VF) V=SQRT(V) AL=AL*FJ FJ=SQRT(AL) U=(V*Q)/FJ C C TX=AY(4,I)*VF TZ=-AY(5,I)*VF IF(IND.NE.3)GO TO 8 IF(MREG.EQ.1)GO TO 8 VP1=DS(4,IREF) VS1=DS(5,IREF) RO2=0. P=DS(2,IREF)/VF IF(NSH.EQ.0)GO TO 14 CALL COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,10,ND,RY,PHSY) USH=(V*QSH)/FJ AAY=USH*RY PHY=PHSY+PHSH+3.14159265 IF(KSH.EQ.1)RETURN 14 TG=DS(2,IREF) TH=DS(3,IREF) ND=0 IF(P.LT.0.)ND=1 P=ABS(P) IF(IC2.LT.0)GO TO 6 CALL COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,5,ND,RX,PHX) CALL COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,6,ND,RZ,PHZ) GO TO 7 6 CONTINUE CALL COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,7,ND,RX,PHX) CALL COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,8,ND,RZ,PHZ) 7 AX=U*RX AZ=U*RZ PHX=PHX+PH PHZ=PHZ+PH AUX1=TX*TH-TZ*TG IF(ABS(AUX1).GT..00001)GO TO 13 PHX=PHX+3.14159265 PHZ=PHZ+3.14159265 GO TO 9 13 AUX2=TG*TX+TH*TZ AUX3=AX*COS(PHX) AUX4=AZ*COS(PHZ) AUX5=AX*SIN(PHX) AUX6=AZ*SIN(PHZ) AR1=(AUX3*AUX2+AUX4*AUX1) AI1=(AUX5*AUX2+AUX6*AUX1) AR2=(AUX4*AUX2-AUX3*AUX1) AI2=(AUX6*AUX2-AUX5*AUX1) AX=SQRT(AR1*AR1+AI1*AI1) PHX=ATAN2(AI1,AR1) AZ=SQRT(AR2*AR2+AI2*AI2) PHZ=ATAN2(AI2,AR2) GO TO 9 8 IF(IC2.LT.0)GO TO 17 AX=U*TX AZ=-U*TZ GO TO 19 17 IF(KSH.EQ.1)GO TO 18 AX=-U*TZ AZ=-U*TX IF(NSH.EQ.0)GO TO 19 18 AAY=(V*QSH)/FJ PHY=PHSH IF(AAY.LT.0.)PHY=PHSH-3.14159 IF(AAY.LT.0.)AAY=-AAY IF(KSH.EQ.1)RETURN 19 PHX=PH PHZ=PH IF(AX.LT.0.)PHX=PH-3.14159 IF(AX.LT.0.)AX=-AX IF(AZ.LT.0.)PHZ=PH-3.14159 IF(AZ.LT.0.)AZ=-AZ 9 RETURN END C C *********************************************************** C SUBROUTINE JPAR(Y,AMULT,KMAH) C C DYNAMIC RAY TRACING C INTEGER CODE DIMENSION Y(15) COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1 COMMON/COD/CODE(50),KREF,KC COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,KPRINT,MPRINT, 1NTR,IMET C IRF=0 AMULT=1. KMAH=0 NDIM=6 IF(IMET.EQ.1)NDIM=15 IREF1=IREF X=AY(1,1) Y(1)=1. Y(2)=0. Y(3)=0. Y(4)=1. DO 10 I=5,NDIM 10 Y(I)=0. 2 I1=KINT(IREF1) IF(I1.LE.0)IREF1=IREF1-1 IF(I1.LE.0)GO TO 2 8 IRF=IRF+1 I1=KINT(IRF) IF(I1.LE.0)GO TO 8 1 I2=I1+1 IRF1=IRF 3 IRF1=IRF1+1 IF(IRF1.GE.IREF1)GO TO 5 I3=KINT(IRF1) IF(I3.LT.0)I2=I2+2 IF(I3.LT.0)GO TO 3 5 CONTINUE CALL RK(X,Y,IRF,NDIM,KMAH) IF(N.EQ.KINT(IREF1))RETURN D11=.5*DS(1,IRF) IF(KREF.LE.1)GO TO 6 6 V1=AY(6,I1) VX=AY(7,I1) VY=AY(8,I1) CS=AY(5,I1) SN=AY(4,I1) AKAPA1=VX*CS-VY*SN VS1=(VX*SN+VY*CS)*V1 V2=AY(6,I2) VX=AY(7,I2) VY=AY(8,I2) CS=AY(5,I2) SN=AY(4,I2) AKAPA2=VX*CS-VY*SN VS2=(VX*SN+VY*CS)*V2 CS1=DS(2,IRF) SN1=-DS(3,IRF) SN2=-DS(11,IRF) S1=2.*(AKAPA1*SN1-AKAPA2*SN2)*CS1/V1 S1=S1+2.*D11*(SN1/V1-SN2/V2) S1=S1+(VS1-VS2)*(CS1*CS1/(V1*V1)) Y(2)=(Y(2)*SN1-Y(1)*S1/SN1)/SN2 Y(1)=Y(1)*SN2/SN1 AMULT=AMULT*SN1/SN2 Y(4)=(Y(4)*SN1-Y(3)*S1/SN1)/SN2 Y(3)=Y(3)*SN2/SN1 7 IRF=IRF+1 I1=KINT(IRF) IF(I1.LE.0)GO TO 7 GO TO 1 C END C C *********************************************************** C SUBROUTINE RK(X,Y,IRF,NDIM,KMAH) C C MODIFIED EULER'S METHOD TO SOLVE A SYSTEM OF LINEAR C ORDINARY DIFFERENTIAL EQUATIONS OF FIRST ORDER C DIMENSION Y(15),DERY(15),Y1(15),Y2(15) COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1 C I1=KINT(IRF) N=N+1 1 H=AY(1,N+1)-AY(1,N) X=X+H CALL FCTA(Y,DERY,NDIM) DO 2 I=1,NDIM Y1(I)=DERY(I) 2 Y2(I)=Y(I)+H*Y1(I) N=N+1 CALL FCTA(Y2,DERY,NDIM) AUX=Y(3) DO 3 I=1,NDIM 3 Y(I)=Y(I)+.5*H*(Y1(I)+DERY(I)) IF(AUX*Y(3).LT.0.)KMAH=KMAH+1 IF(N.EQ.I1) GO TO 4 GO TO 1 4 RETURN END C C *********************************************************** C SUBROUTINE FCTA(Y,DERY,NDIM) C C COMPUTATION OF THE RIGHT-HAND SIDES OF THE DYNAMIC C RAY TRACING SYSTEM EQUATIONS C DIMENSION Y(NDIM),DERY(NDIM) COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1 COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT, 1NTR,IMET C V=AY(6,N) VV=V*V CS=AY(5,N)*V SN=AY(4,N)*V VXX=AY(9,N) VXY=AY(10,N) VYY=AY(11,N) VN=(AY(7,N)*SN-AY(8,N)*CS)**2 V22=VXX*CS*CS-2.*VXY*CS*SN+VYY*SN*SN DERY(1)=VV*Y(2) DERY(2)=-(V22*Y(1))/V DERY(3)=VV*Y(4) DERY(4)=-(V22*Y(3))/V DERY(5)=0. IF(ABS(AY(12,N)).GT..1)DERY(5)=1./AY(12,N) DERY(6)=VV IF(IMET.NE.1)RETURN DERY(7)=VV*Y(1)*Y(1) DERY(8)=VV*Y(1)*Y(3) DERY(9)=VV*Y(3)*Y(3) DERY(10)=VV*Y(2)*Y(2) DERY(11)=VV*Y(2)*Y(4) DERY(12)=VV*Y(4)*Y(4) DERY(13)=VN*Y(1)*Y(1) DERY(14)=VN*Y(1)*Y(3) DERY(15)=VN*Y(3)*Y(3) RETURN END C C *********************************************************** C SUBROUTINE COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,ND,RMOD,RPH) C C THE ROUTINE COEF8 IS DESIGNED FOR THE COMPUTATION OF REFLECTION C AND TRANSMISSION COEFFICIENTS AT A PLANE INTERFACE BETWEEN TWO C HOMOGENEOUS SOLID HALFSPACES OR AT A FREE SURFACE OF A HOMOGENEOUS C SOLID HALFSPACE. C COMPLEX B(4),RR,C1,C2,C3,C4,H1,H2,H3,H4,H5,H6,H,HB,HC DIMENSION PRMT(4),D(4),DD(4) IF(NCODE.GE.9)GO TO 300 IF(RO2.LT.0.000001)GO TO 150 PRMT(1)=VP1 PRMT(2)=VS1 PRMT(3)=VP2 PRMT(4)=VS2 A1=VP1*VS1 A2=VP2*VS2 A3=VP1*RO1 A4=VP2*RO2 A5=VS1*RO1 A6=VS2*RO2 Q=2.*(A6*VS2-A5*VS1) PP=P*P QP=Q*PP X=RO2-QP Y=RO1+QP Z=RO2-RO1-QP G1=A1*A2*PP*Z*Z G2=A2*X*X G3=A1*Y*Y G4=A4*A5 G5=A3*A6 G6=Q*Q*PP DO 21 I=1,4 DD(I)=P*PRMT(I) 21 D(I)=SQRT(ABS(1.-DD(I)*DD(I))) IF(DD(1).LE.1..AND.DD(2).LE.1..AND.DD(3).LE.1..AND.DD(4).LE.1.) 1GO TO 100 C C COMPLEX COEFFICIENTS DO 22 I=1,4 IF(DD(I).GT.1.)GO TO 23 B(I)=CMPLX(D(I),0.) GO TO 22 23 B(I)= CMPLX(0.,D(I)) 22 CONTINUE C1=B(1)*B(2) C2=B(3)*B(4) C3=B(1)*B(4) C4=B(2)*B(3) H1=G1 H2=G2*C1 H3=G3*C2 H4=G4*C3 H5=G5*C4 H6=G6*C1*C2 H=1./(H1+H2+H3+H4+H5+H6) HB=2.*H HC=HB*P GO TO (1,2,3,4,5,6,7,8),NCODE 1 RR=H*(H2+H4+H6-H1-H3-H5) GO TO 26 2 RR=VP1*B(1)*HC*(Q*Y*C2+A2*X*Z) IF(ND.NE.0)RR=-RR GO TO 26 3 RR=A3*B(1)*HB*(VS2*B(2)*X+VS1*B(4)*Y) GO TO 26 4 RR=-A3*B(1)*HC*(Q*C4-VS1*VP2*Z) IF(ND.NE.0)RR=-RR GO TO 26 5 RR=-VS1*B(2)*HC*(Q*Y*C2+A2*X*Z) IF(ND.NE.0)RR=-RR GO TO 26 6 RR=H*(H2+H5+H6-H1-H3-H4) GO TO 26 7 RR=A5*B(2)*HC*(Q*C3-VP1*VS2*Z) IF(ND.NE.0)RR=-RR GO TO 26 8 RR=A5*B(2)*HB*(VP1*B(3)*Y+VP2*B(1)*X) GO TO 26 C REAL COEFFICIENTS 100 E1=D(1)*D(2) E2=D(3)*D(4) E3=D(1)*D(4) E4=D(2)*D(3) S1=G1 S2=G2*E1 S3=G3*E2 S4=G4*E3 S5=G5*E4 S6=G6*E1*E2 S=1./(S1+S2+S3+S4+S5+S6) SB=2.*S SC=SB*P GO TO (101,102,103,104,105,106,107,108),NCODE 101 R=S*(S2+S4+S6-S1-S3-S5) GO TO 250 102 R=VP1*D(1)*SC*(Q*Y*E2+A2*X*Z) IF(ND.NE.0)R=-R GO TO 250 103 R=A3*D(1)*SB*(VS2*D(2)*X+VS1*D(4)*Y) GO TO 250 104 R=-A3*D(1)*SC*(Q*E4-VS1*VP2*Z) IF(ND.NE.0)R=-R GO TO 250 105 R=-VS1*D(2)*SC*(Q*Y*E2+A2*X*Z) IF(ND.NE.0)R=-R GO TO 250 106 R=S*(S2+S5+S6-S1-S3-S4) GO TO 250 107 R=A5*D(2)*SC*(Q*E3-VP1*VS2*Z) IF(ND.NE.0)R=-R GO TO 250 108 R=A5*D(2)*SB*(VP1*D(3)*Y+VP2*D(1)*X) GO TO 250 C C EARTHS SURFACE,COMPLEX COEFFICIENTS AND COEFFICIENTS OF CONVERSION 150 A1=VS1*P A2=A1*A1 A3=2.*A2 A4=2.*A1 A5=A4+A4 A6=1.-A3 A7=2.*A6 A8=2.*A3*VS1/VP1 A9=A6*A6 DD(1)=P*VP1 DD(2)=P*VS1 DO 151 I=1,2 151 D(I)=SQRT(ABS(1.-DD(I)*DD(I))) IF(DD(1).LE.1..AND.DD(2).LE.1.)GO TO 200 DO 154 I=1,2 IF(DD(I).GT.1.)GO TO 155 B(I)=CMPLX(D(I),0.) GO TO 154 155 B(I)= CMPLX(0.,D(I)) 154 CONTINUE H1=B(1)*B(2) H2=H1*A8 H=1./(A9+H2) GO TO (161,162,163,164,165,166,167,168),NCODE 161 RR=(-A9+H2)*H GO TO 26 162 RR=-A5*B(1)*H*A6 IF(ND.NE.0)RR=-RR GO TO 26 163 RR=A5*B(2)*H*A6*VS1/VP1 IF(ND.NE.0)RR=-RR GO TO 26 164 RR=-(A9-H2)*H GO TO 26 165 RR=A5*H1*H IF(ND.NE.0)RR=-RR GO TO 26 166 RR=-A7*B(1)*H GO TO 26 167 RR=A7*B(2)*H GO TO 26 168 RR=A5*VS1*H1*H/VP1 IF(ND.NE.0)RR=-RR 26 Z2=REAL(RR) Z3=AIMAG(RR) IF(Z2.EQ.0..AND.Z3.EQ.0.)GO TO 157 RMOD=SQRT(Z2*Z2+Z3*Z3) RPH=ATAN2(Z3,Z2) RETURN 157 RMOD=0. RPH=0. RETURN C C EARTHS SURFACE,REAL COEFFICIENTS AND COEFFICIENTS OF CONVERSION 200 S1=D(1)*D(2) S2=A8*S1 S=1./(A9+S2) GO TO (201,202,203,204,205,206,207,208),NCODE 201 R=(-A9+S2)*S GO TO 250 202 R=-A5*D(1)*S*A6 IF(ND.NE.0)R=-R GO TO 250 203 R=A5*D(2)*S*A6*VS1/VP1 IF(ND.NE.0)R=-R GO TO 250 204 R=(S2-A9)*S GO TO 250 205 R=A5*S1*S IF(ND.NE.0)R=-R GO TO 250 206 R=-A7*D(1)*S GO TO 250 207 R=A7*D(2)*S GO TO 250 208 R=A5*VS1*S1*S/VP1 IF(ND.NE.0)R=-R 250 IF(R.LT.0.)GO TO 251 RMOD=R RPH=0. RETURN 251 RMOD=-R RPH=-3.14159 RETURN C C SH COEFFICIENTS OF REFLECTION, TRANSMISSION AND CONVERSION C 300 if(ro2.lt..000001)go to 302 A1=P*VS1 A2=P*VS2 A3=RO1*VS1 A4=RO2*VS2 A5=SQRT(ABS(1.-A1*A1)) A6=SQRT(ABS(1.-A2*A2)) C1=A5 IF(A2.LE.1.)C2=A6 IF(A2.GT.1.)C2=CMPLX(0.,A6) C1=C1*A3 C2=C2*A4 H=1./(C1+C2) IF(NCODE.EQ.10)GO TO 301 RR=H*(C1-C2) GO TO 26 301 RR=2.*C1*H GO TO 26 302 if(ncode.eq.10)go to 303 RMOD=1. RPH=0. RETURN 303 RMOD=2. RPH=0. RETURN END C C ************************************************************* C SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) DIMENSION Y(1),DERY(1),AUX(8,1),A(4),B(4),C(4),PRMT(1) DO 1 I=1,NDIM 1 AUX(8,I)=.0666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) IF(H*(XEND-X))38,37,2 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0.0 IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GO TO 10 15 IF(ITEST)16,16,20 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GO TO 9 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GO TO 9 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GO TO 18 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GO TO 4 36 IHLF=11 CALL FCT(X,Y,DERY) GO TO 39 37 IHLF=12 GO TO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END C