C PROGRAM A N R A Y, VERSION 3.02 C C******************************************************************* C C PROGRAM ANRAY IS DESIGNED FOR RAY, TRAVEL TIME AND C AMPLITUDE COMPUTATIONS IN 3D GENERAL ANISOTROPIC AND ISOTROPIC C LATERALLY VARYING LAYERED MEDIA. THE PROGRAM MAKES POSSIBLE C COMPUTATION OF RAYS SPECIFIED BY INITIAL ANGLES AT THE SOURCE, C I.E., INITIAL-VALUE RAY TRACING, OR RAYS STARTING FROM THE C SOURCE AND TERMINATING ON A VERTICAL OR SURFACE PROFILE, I.E. C BOUNDARY-VALUE RAY TRACING. RAY AMPLITUDES CAN BE COMPUTED C ALONG RAYS. C C******************************************************************* C C character*80 mtext DIMENSION Y(18) COMMON /AUXI/ IANI(20),INTR,INT1,IOUT,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 /AUXX/ MMX(20),MMY(20),MMXY(20) 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 INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DIST/ DST(200),NDST,REPS,PROF(2),NDSTP,PREPS,LNDST, 1xprf,yprf 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),DINC,DTOLD,N,IREF,IND,IND1 COMMON /ZERO/ RNULL COMMON/VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP C C PI=3.14159265359 C C************************************************** C MODE=0 CALL SERV(MODE,LIN,LOU,LLU1,LLU2,LLU3) IF(MODE.EQ.0)LIN=5 IF(MODE.EQ.0)LOU=6 LOUT=LOU C C************************************************** C WRITE(LOU,777) 777 FORMAT(///,'***************************' 1,//,' PROGRAM A N R A Y 8 9 ',//, 2'***************************',//) NCODE=1 1 READ(LIN,115)MTEXT WRITE(LOU,115)MTEXT READ(LIN,106)LU1,LU2,INULL,iplot IF(INULL.EQ.0)INULL=4 RNULL=10.**(-INULL) WRITE(LOU,106)LU1,LU2,INULL,iplot IF(MODE.NE.0)LU1=LLU1 IF(MODE.NE.0)LU2=LLU2 IF(MODE.NE.0)LU3=LLU3 C C C SPECIFICATION OF THE MODEL C CALL MODEL(MTEXT,LIN) C C GENERATE FILE FOR PLOTTING VARIOUS CHARACTERISTIC SURFACES C if(iplot.ne.0)then call surfpl(lin,lu3) go to 99 end if C C SPECIFICATION OF SYNTHETIC SEISMOGRAMS C 2 READ(LIN,100)ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX, 1IPOL,ISPR,IRAYPL,IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,mori WRITE(LOU,102)ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX, 1IPOL,ISPR,IRAYPL,IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,mori IF(ICONT.EQ.(-1))GO TO 1 IF(ICONT.EQ.0)GO TO 99 C C IF(ITMAX.EQ.0)ITMAX=10 IVSP=0 IF(ILOC.EQ.0)ITPR=3 IF(ILOC.EQ.1)THEN IVSP=1 ITPR=43 MREG=1 END IF IF(ILOC.GT.1)THEN ITPR=ILOC+100 MREG=1 END IF C IF(MEP.EQ.0)THEN NDST=0 END IF C if(mep.eq.1)then ndst=1 read(lin,104)xrec,yrec write(lou,104)xrec,yrec go to 4 end if IF(MEP.LT.0)THEN NDST=-MEP READ(LIN,104)PROF(1),(DST(I),I=1,NDST),xprf,yprf WRITE(LOU,104)PROF(1),(DST(I),I=1,NDST),xprf,yprf IF(NDST.EQ.1)RSTEP=1. IF(NDST.EQ.1)DST(2)=DST(1)+1. IF(NDST.EQ.1)GO TO 4 RSTEP=(DST(NDST)-DST(1))/FLOAT(NDST-1) END IF C IF(MEP.GT.0)THEN NDST=MEP READ(LIN,104)PROF(1),RMIN,RSTEP,xprf,yprf WRITE(LOU,104)PROF(1),RMIN,RSTEP,xprf,yprf DO 13 I=1,MEP 13 DST(I)=RMIN+(I-1)*RSTEP if(ndst.eq.1)dst(2)=rmin+rstep END IF PROF(2)=PROF(1)+1. NDSTP=1 C IF(IVSP.EQ.1.AND.NDST.NE.0)THEN READ(LIN,104)XVSP,YVSP WRITE(LOU,104)XVSP,YVSP END IF C 4 READ(LIN,104)XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS WRITE(LOU,104)XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS C if(abs(xprf).lt..000001.and.abs(yprf).lt..000001)then xprf=xsour yprf=ysour end if if(mep.eq.1)then xe=xrec-xprf ye=yrec-yprf rprf=sqrt(xe*xe+ye*ye) xatan=atan2(ye,xe) prof(1)=xatan rmin=rprf dst(1)=rmin write(lou,104)rprf,xatan RSTEP=100. DST(2)=DST(1)+100. prof(2)=prof(1)+1. ndstp=1 end if C IF(IVSP.EQ.1.AND.NDST.NE.0)THEN XNRM=XVSP-XSOUR YNRM=YVSP-YSOUR PROF(1)=ATAN2(YNRM,XNRM) PROF(2)=PROF(1)+1. xprf=xsour yprf=ysour END IF IF(MCOD.EQ.0)THEN READ(LIN,104)AMIN,ASTEP,AMAX WRITE(LOU,104)AMIN,ASTEP,AMAX READ(LIN,104)BMIN,BSTEP,BMAX if(abs(bstep).lt..000001)then bmin=prof(1)-.3 bmax=prof(1)+.4 bstep=.6 end if WRITE(LOU,104)BMIN,BSTEP,BMAX NRAYA=IFIX(ABS(AMAX-AMIN)/ASTEP) END IF IF(MREG.EQ.0.AND.MDIM.NE.0) WRITE(LOU,'(/,A,/)') 1 ' COEFFICIENTS OF CONVERSION ARE APPLIED' IF(MREG.NE.0.AND.MDIM.NE.0) WRITE(LOU,'(/,A,/)') 1 ' COEFFICIENTS OF CONVERSION ARE *** NOT *** APPLIED' IF(DT.LT.0.0000001)DT=1. IF(AC.LT.0.0000001)AC=0.0001 TMAX=10000. IF(REPS.LT..00001)REPS=.05 IF(PREPS.LT..00001)PREPS=.05 IND=-1 NDER=1 CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,AMIN1,BMIN,PX,PY,PZ,XX,YY,ZZ,T, 1DT,AC) Y(1)=XSOUR Y(2)=YSOUR Y(3)=ZSOUR IF(IND.EQ.50)WRITE(LOU,111)IND IF(IND.EQ.50)GO TO 99 LAY=IND ISOUR=IND itype=3 CALL PARDIS(Y,0) VP=SQRT(A11) IF(IRHO.EQ.0)RO=1.7+.2*VP IF(IRHO.EQ.1)RO=RHO(IND) C C GENERATE FILE LU2 FOR SYNTHETIC SEISMOGRAM COMPUTATIONS C IF(LU2.NE.0.AND.NDST.NE.0)THEN WRITE(LU2,115)MTEXT KSH=2 WRITE(LU2,100)NDST,KSH,ILOC WRITE(LU2,104)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,RO if(mep.ne.1)WRITE(LU2,104)(DST(I),I=1,NDST) if(mep.eq.1)WRITE(LU2,104)xrec,yrec,rprf,xatan END IF C C LOOP FOR ELEMENTARY WAVES C 20 READ(LIN,100)KC,KREF,((CODE(I,K),K=1,2),I=1,KREF) WRITE(LOU,100)KC,KREF,((CODE(I,K),K=1,2),I=1,KREF) IF(KREF.EQ.0)GOTO 2 IF(MOUT.NE.0)WRITE(LOU,107) WRITE(LOU,103)NCODE,KC,KREF,((CODE(I,K),K=1,2),I=1,KREF) C IF(MCOD.NE.0)THEN READ(LIN,104)AMIN,ASTEP,AMAX WRITE(LOU,104)AMIN,ASTEP,AMAX READ(LIN,104)BMIN,BSTEP,BMAX if(abs(bstep).lt..000001)then bmin=prof(1)-.3 bmax=prof(1)+.4 bstep=.6 end if NRAYA=IFIX(ABS(AMAX-AMIN)/ASTEP) WRITE(LOU,104)BMIN,BSTEP,BMAX END IF C C GENERATE FILE LU1 FOR PLOTTING OF RAY DIAGRAMS, C TIME-DISTANCE AND AMPLITUDE-DISTANCE CURVES C IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 21 WRITE(LU1,100)ICONT,NDST,ILOC WRITE(LU1,104)RO NPN=2 APN=0. WRITE(LU1,100)NPN,NPN,NPN WRITE(LU1,101)APN,APN,APN,APN,APN WRITE(LU1,101)APN,APN,APN,APN,APN WRITE(LU1,104)Xprf,Yprf,0.0,PROF(1) WRITE(LU1,104)(DST(I),I=1,NDST) 21 CONTINUE C C C SEARCH FOR THE NUMBER OF THE ELEMENT OF THE RAY, STARTING FROM C WHICH THE WAVE DOES UNDERTAKE NEITHER REFLECTION NOR CONVERSION C ICOD=0 IF(IVSP.EQ.0)GO TO 35 DO 34 I=1,KREF ICOD=KREF-I+1 IF(ICOD.EQ.1) GO TO 34 IC1=CODE(ICOD,1) IC2=CODE(ICOD-1,1) IF((IC1-IC2).EQ.0)GO TO 35 IC1=CODE(ICOD,2) IC2=CODE(ICOD-1,2) IF((IC1-IC2).NE.0)GO TO 35 34 CONTINUE 35 CONTINUE IF(MOUT.NE.0)WRITE(LOU,108) C C CALL RECEIV(XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,ITMAX,AMIN,ASTEP, 1AMAX,BMIN,BSTEP,BMAX,MOUT,LU1,LU2,METHOD,ITPR,NCODE) IF(IND.EQ.20)WRITE(LOU,111)IND IF(IND.EQ.20)GO TO 99 IF(IND.EQ.14) WRITE(LOU,111) IND NCODE=NCODE+1 GOTO 20 C C END OF LOOP FOR ELEMENTARY WAVES C C 100 FORMAT(26I3) 101 FORMAT(5E15.5) 102 FORMAT(1H0,////,2X,26I3) 104 FORMAT(8F10.5) 106 FORMAT(17I5) 103 FORMAT(4X,I4,9X,36I3) 107 FORMAT(//2X,'INT.CODE',5X,'E X T E R N A L C O D E') 108 FORMAT(//) 111 FORMAT(/2X,'IND=',I5,/) 115 FORMAT(A) C 99 CONTINUE IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)ICONT,ICONT IF(LU1.NE.0)REWIND LU1 IF(LU2.NE.0)REWIND LU2 C STOP END