SUBROUTINE RAYA (X0,Y0,Z0,T0,FI0,PSI0,PX,PY,PZ,XX,YY,ZZ,T,DT,AC) C C INITIAL-VALUE RAY TRACING BY THE RUNGE-KUTTA METHOD C DIMENSION Y(18),DEP(6),PRM(5),DERY(18),AUX(8,18),DIN(18),VSQ(3) dimension pn(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,IOUT,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 /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),brd(4),NINT COMPLEX PS COMMON /RAY/ AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax, 1 PS(3,7,20),IS(8,20),DINC,DTOLD,N,IREF,IND,IND1 COMMON /ZERO/ RNULL 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 EXTERNAL FCT,OUT C IAC=0 kmah=0 Y(1)=X0 Y(2)=Y0 Y(3)=Z0 IREFR=0 KRE=KREF IF(KC.EQ.0) KRE=0 N=0 IREF=1 IOUT=0 IF(IND.GT.0) GOTO 6 C C FOR IND=-1 DETERMINATION OF THE NUMBER OF THE LAYER IN WHICH THE C INITIAL POINT OF THE RAY IS SITIUATED C IF(Y(1).LT.BRD(1).OR.Y(1).GT.BRD(2).OR.Y(2).LT.BRD(3).OR.Y(2).GT. 1 BRD(4)) GOTO 4 INTR=1 1 CALL DISC (Y,DEP) ZINT=DEP(1) IF(ABS(Y(3)).GT.0.00001) GOTO 2 ISOUR=1 C C REDEFINITION OF Z-COORDINATES FOR A SOURCE ON THE EARTH SURFACE C Z0=ZINT GOTO 5 2 IF(Y(3).LT.ZINT.AND.INTR.EQ.1) GOTO 4 IF(Y(3).LT.ZINT) GOTO 5 IF(ABS(Y(3)-ZINT).LT.RNULL.AND.INTR.EQ.NINT) GOTO 5 IF(INTR.EQ.NINT) GOTO 4 ISOUR=INTR INTR=INTR+1 GOTO 1 4 WRITE(lou,'(A,/,6F10.5,/,4F10.5)')' Y,BRD',Y,BRD IND=50 RETURN C C DETERMINATION OF INITIAL CONDITIONS FOR THE RUNGE-KUTTA PROCEDURE C 5 IF(IND.GE.0) GOTO 6 IND=ISOUR RETURN 6 LAY=ISOUR is(3,1)=lay INT1=ISOUR IF(ISOUR.NE.CODE(1,1)) IND=14 IF(ISOUR.NE.CODE(1,1)) RETURN ITYPE=CODE(1,2) IF(IANI(ISOUR).EQ.0.AND.ITYPE.EQ.2.and.kc.ne.0) IND=38 IF(IANI(ISOUR).EQ.0.AND.ITYPE.EQ.2.and.kc.ne.0) RETURN CALL PARDIS(Y,0) C C DETERMINATION OF INITIAL VALUES FOR RAY TRACING C AND DYNAMIC RAY TRACING C csp=cos(psi0) snp=sin(psi0) csf=cos(fi0) snf=sin(fi0) if(mori.eq.0)then Y(4)=CSP*CSF Y(5)=SNP*CSF Y(6)=SNF else y(4)=csf*csp y(5)=snp y(6)=snf*csp end if do 8 k=1,3 pn(k)=y(k+3) 8 continue if(nder.gt.1)then do 3 k=7,12 y(k)=0. 3 continue if(mori.eq.0)then y(13)=-snp*csf y(14)=csp*csf y(15)=0. y(16)=-csp*snf y(17)=-snp*snf y(18)=csf else y(13)=-snf*csp y(14)=0. y(15)=csf*csp y(16)=-csf*snp y(17)=csp y(18)=-snf*snp end if end if IF(IANI(ISOUR).ne.0)then C C SOURCE LOCATED IN AN ANISOTROPIC LAYER C CALL INIT(pn,VSQ) IF(IPRINT.GT.2)WRITE(lou,'(a,3F14.6)')' V1,V2,V3=', VSQ VP=AMAX1(VSQ(1),VSQ(2),VSQ(3)) VS1=AMIN1(VSQ(1),VSQ(2),VSQ(3)) VS2=VSQ(1)+VSQ(2)+VSQ(3)-VP-VS1 IF(ITYPE.EQ.3)V=SQRT(VP) IF(ITYPE.EQ.1)V=SQRT(VS1) IF(ITYPE.EQ.2)V=SQRT(VS2) do 7 i=4,6 y(i)=y(i)/v 7 continue if(nder.gt.1)then nder=1 call fct(0.,y,dery) nder=2 vg=sqrt(dery(1)*dery(1)+dery(2)*dery(2)+dery(3)*dery(3)) if(mori.eq.0)then auxp=-dery(1)*snp*csf+dery(2)*csp*csf auxf=-dery(1)*csp*snf-dery(2)*snp*snf+dery(3)*csf y(13)=y(13)-auxp*csp*csf/v y(14)=y(14)-auxp*snp*csf/v y(15)=y(15)-auxp*snf/v y(16)=y(16)-auxf*csp*csf/v y(17)=y(17)-auxf*snp*csf/v y(18)=y(18)-auxf*snf/v else auxp=-dery(1)*snf*csp+dery(3)*csf*csp auxf=-dery(1)*csf*snp+dery(2)*csp-dery(3)*snf*snp y(13)=y(13)-auxp*csf*csp/v y(14)=y(14)-auxp*snp/v y(15)=y(15)-auxp*snf*csp/v y(16)=y(16)-auxf*csf*csp/v y(17)=y(17)-auxf*snp/v y(18)=y(18)-auxf*snf*csp/v end if do 11 i=13,18 y(i)=y(i)/v 11 continue C C DETERMINATION OF KMAH INDEX AT SOURCE IN ANISTROPIC MEDIUM C call fct(0.,y,dery) if(mori.eq.0)then el=-(-dery(7)*snp*csf+dery(8)*csp*csf)/vg em=-(-dery(10)*snp*csf+dery(11)*csp*csf)/vg en=-(-dery(10)*csp*snf-dery(11)*snp*snf+dery(12)*csf)/vg else el=-(-dery(7)*snf*csp+dery(9)*csf*csp)/vg em=-(-dery(10)*snf*csp+dery(12)*csf*csp)/vg en=-(-dery(10)*csf*snp+dery(11)*csp-dery(12)*snf*snp)/vg end if ee=y(13)*y(13)+y(14)*y(14)+y(15)*y(15) ff=y(13)*y(16)+y(14)*y(17)+y(15)*y(18) gg=y(16)*y(16)+y(17)*y(17)+y(18)*y(18) egf=(ee*gg-ff*ff)/v/v elnm=el*en-em*em be=el*gg+en*ee-2.*em*ff if(abs(egf).lt..000001)egf=sign(1.,egf)*.00001 gk=elnm/egf ak=.5*be/egf if(gk.gt.0.)then if(ak.lt.0.)kmah=0 if(ak.gt.0.)kmah=2 end if if(gk.lt.0.)kmah=1 end if end if C C SOURCE LOCATED IN AN ISOTROPIC LAYER C IF(IANI(ISOUR).eq.0)then IF(ITYPE.EQ.3)V=SQRT(A11) IF(ITYPE.NE.3)V=SQRT(A44) kdim=6 if(nder.gt.1)kdim=18 do 9 i=4,kdim y(i)=y(i)/v 9 continue end if C IND=0 IND1=0 PRM(1)=T0 PRM(2)=TMAX PRM(3)=DT IF(ITYPE.NE.3) PRM(3)=DT*1.7 PRM(4)=AC T=PRM(1) 20 DO 10 I=1,3 auxx=y(4)*y(4)+y(5)*y(5)+y(6)*y(6) auxx=sqrt(auxx) DIN(I)=auxx din(i+3)=prm(3)/auxx 10 CONTINUE kdim=6 if(nder.gt.1)kdim=18 do 25 i=4,kdim din(i)=0. 25 continue DO 30 I=1,kdim DERY(I)=DIN(I) 30 CONTINUE C C COMPUTATION OF THE RAY C CALL RKGS(PRM,Y,DERY,kdim,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(ABS(PRM(5)).GT.0.0001) GOTO 35 IF(IND.eq.12) GOTO 70 GOTO 60 35 IF(ABS(PRM(5)-1.1).GT.0.0001) GOTO 50 C C LOOP FOR THE ITERATIVE DETERMINATION OF THE POINT OF INCIDENCE C PRM(1)=AY(1,N) DO 40 I=1,6 Y(I)=Yold(I) 40 CONTINUE if(nder.gt.1)then do 37 i=1,12 y(i+6)=yold(i+6) 37 continue end if T=AY(1,N) N=N-1 GOTO 20 50 XX=Y(1) YY=Y(2) ZZ=Y(3) T=AY(1,N) IF(ABS(PRM(5)-2.).GT.0.0001) GOTO 80 C C INTEGRATION FROM THE POINT OF REFLECTION/TRANSMISSION TO THE CLOSEST C POINT OF THE RAY CORRESPONDING TO REGULAR TIME STEP C PRM(1)=AY(1,N) I=INT((PRM(1)-T0)/DT) PRM(2)=FLOAT(I+1)*DT+T0 PRM(3)=DT GOTO 20 60 PRM(1)=PRM(2) PRM(2)=TMAX PRM(3)=DT IF(ITYPE.NE.3) PRM(3)=1.7*DT N=N-1 GOTO 20 70 CONTINUE XX=Y(1) YY=Y(2) ZZ=Y(3) 80 continue c if(kmah.ne.0)ind1=ind1+50 IF(IREFR.EQ.1) IND1=-IND1 PX=Y(4) PY=Y(5) PZ=Y(6) if(nder.gt.1)then do 90 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) 90 continue end if RETURN END