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
