C     PROGRAM A N R A Y P L
C     *********************
C
C     THE PROGRAM ANRAYPL IS DESIGNED FOR PLOTTING OF RAY DIAGRAMS,
C     TRAVEL TIMES AND AMPLITUDES OF SEISMIC BODY WAVES FROM THE
C     FILE LU1 GENERATED BY PROGRAM ANRAY89
C
C     ****************************************************************
C
      character*80 text
      complex ax,ay,az,aux,caux,csour,imag
      DIMENSION A(30),B(30),C(30),D(30),X1(30),NPNT(20),NR(40)
      DIMENSION X(500),T(500),XZ(500),AX(2,500),AY(2,500),AZ(2,500),
     1Y(500),Z(500),INDI(500),DST(200),
     2aux(3),am(3),ph(3),p(3),pol(3),pol1(3)
      COMMON/SOUR/ROS
C
C     ***
      mode=0
      call serv(mode,lin,lou,llu,lu2,lu3)
      if(mode.eq.0)call plots(ldum1,ldum2,7)
      if(mode.eq.0)lin=5
      if(mode.eq.0)lou=6
C     ***
C
      imag=(0.,1.)
      XB = 0.
      IRUN = 0
      READ(lin,106)ISHFTx,ishfty,IPRINT,lu
      if(mode.ne.0)lu=llu
      WRITE(lou,106)ISHFTx,ishfty,IPRINT,lu
C
      REWIND LU
      IF(ISHFTx.EQ.0)ISHFTx=8
      IF(ISHFTy.EQ.0)ISHFTy=8
      xSHIFT=FLOAT(ISHFTx)
      ySHIFT=FLOAT(ISHFTy)
      CALL PLOT(3.,3.,-3)
c
  200 READ(LU,101)ICONT,NDST,iloc
      WRITE(lou,101)ICONT,NDST,iloc
      IF(ICONT.EQ.0)WRITE(lou,107)LU
      IF(ICONT.EQ.0)GO TO 99
      read(lu,102)ros
      write(lou,102)ros
      READ(lin,100)TEXT
      WRITE(lou,100)TEXT
      READ(lin,106)NTICX,NTICY,ntich,NTICT,NTICA,IPROF,NRAY,IBOUND,
     1IRED,IRS,NDX,NDY,ndh
      WRITE(lou,106)NTICX,NTICY,ntich,NTICT,NTICA,IPROF,NRAY,IBOUND,
     1IRED,IRS,NDX,NDY,ndh
C
      IF(NTICX.EQ.0)GO TO 99
      IF(IPROF.EQ.0)IPROF=1
      XB = 0.
      IF(IBOUND.EQ.0) IBOUND =100
      READ(lin,102)XMIN,XMAX,XLEN,DTICX,SC,arot
      WRITE(lou,102)XMIN,XMAX,XLEN,DTICX,SC,arot
      csrt=cos(arot)
      snrt=sin(arot)
      IF(ABS(SC).LT..00001)SC=1.
      IF(ABS(XMAX-XMIN).LT..00001)GO TO 32
      XMER = XLEN/(XMAX-XMIN)
      GLH = 1.5
      READ(lin,102)YMIN,YMAX,YLEN,DTICY,hmin,hmax,hlen,dtich
      WRITE(lou,102)YMIN,YMAX,YLEN,DTICY,hmin,hmax,hlen,dtich
      if(nticy.eq.0)go to 40
      YMER = YLEN/(YMAX-YMIN)
C
C     PLOTTING OF BORDER FOR VERTICAL RAY DIAGRAM
C
      IF(IRUN.EQ.1)CALL PLOT(xSHIFT/3.,0.,-3)
      IRUN=1
      CALL BORDER(XLEN,DTICX,YLEN,DTICY,SC,TEXT,1,XMIN,XMAX,YMIN,YMAX,
     1NTICX,NTICY,NDX,NDY)
      TX=.5*(XLEN-6.3*SC)
      CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14)
      TX=.5*(YLEN-4.95*SC)
      U=-(1.6+.4*NDX)*SC
      CALL SYMBOL(U,TX,.45*SC,'DEPTH IN KM',90.,11)
C
   40 if(ntich.eq.0)go to 32
      hMER = hLEN/(hMAX-hMIN)
C
C     PLOTTING OF BORDER FOR HORIZONTAL RAY DIAGRAM
C
      IF(IRUN.EQ.1.and.nticy.ne.0)CALL PLOT(0.,yLEN+ySHIFT,-3)
      IF(IRUN.EQ.1.and.nticy.eq.0)CALL PLOT(XLEN+xSHIFT,0.,-3)
      IRUN=1
      CALL BORDER(XLEN,DTICX,hLEN,DTICh,SC,TEXT,1,XMIN,XMAX,hMIN,hMAX,
     1NTICX,NTICh,NDX,NDh)
      call plot(0.,.5*hlen,3)
      call plot(xlen,.5*hlen,2)
      call plot(0.,0.,3)
      TX=.5*(XLEN-6.3*SC)
      CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14)
      TX=.5*(hLEN-8.1*SC)
      U=-(1.6+.4*NDX)*SC
      CALL SYMBOL(U,TX,.45*SC,'TRANSV.DIST. IN KM',90.,18)
      if(nticy.ne.0)call plot(0.,-ylen-yshift,-3)
C
C     PLOTTING OF INTERFACES
C
   32 READ(LU,101)NINT,(NPNT(I),I=1,NINT)
      IF(IPRINT.EQ.3)WRITE(lou,106)NINT,(NPNT(I),I = 1,NINT)
      IB=IBOUND
      IF(IBOUND.LT.0)IB=-IBOUND
      RD = (XMAX-XMIN)/FLOAT(IB-1)
      YM = YMAX
      DO 5 I = 1,NINT
      N = NPNT(I)-1
      READ(LU,113)(A(J),B(J),C(J),D(J),X1(J),J = 1,N)
      IF(IPRINT.EQ.3)WRITE(lou,105)(A(J),B(J),C(J),D(J),X1(J),
     1J=1,N)
      IF(NTICY.EQ.0)GOTO 5
      MMM=1
      IF(IBOUND.LT.0)YM=YM-.06/YMER
      IF(IBOUND.LT.0)MMM=3
      DO 4 M=1,MMM
      RB=XMIN
      IF(IBOUND.LT.0)YM=YM+.03/YMER
      IPL=0
      DO 4 J = 1,IB
      IF(RB.GT.(XMAX+0.01))GOTO 4
      DO 1 K = 1,N
      IF(K.EQ.N) K1 = K
      IF(RB.GE.X1(K))GOTO 1
      K1 = K-1
      GOTO 2
  1   CONTINUE
  2   X2 = RB-X1(K1)
      X33=X1(K1+1)
      IF(K1.EQ.N)X33=XMAX
      IF(XMAX.LT.X33)X33=XMAX
      X3=X33-RB
      IF(X3.LT.RD)RB=X33
      IF(X3.LT.RD)X2=X33-X1(K1)
      HXX=YM-(((D(K1)*X2+C(K1))*X2+B(K1))*X2+A(K1)-YMIN)
      IF(HXX.LT.YMIN.OR.HXX.GT.YMAX)IPL=0
      IF(HXX.LT.YMIN.OR.HXX.GT.YMAX)GO TO 4
      IF(IPL.NE.0)GO TO 3
      IPL=1
      XSV=(RB-XMIN)*XMER
      YSV = HXX*YMER
      CALL PLOT(XSV,YSV,3)
      GOTO 4
  3   XSV = (RB-XMIN)*XMER
      YSV = HXX*YMER
      CALL PLOT(XSV,YSV,2)
  4   RB = RB+RD
      IF(IBOUND.LT.0)YM=YM-0.03/YMER
  5   CONTINUE
   20 READ(LU,102)X0,Y0,Z0,fi
      READ(LU,102)(DST(I),I=1,NDST)
      IF(IPRINT.EQ.3)WRITE(lou,102)X0,Y0,Z0,fi
      IF(IPRINT.EQ.3)WRITE(lou,102)(DST(I),I=1,NDST)
      IF(NTICY.EQ.0)GOTO 21
      CSF=COS(FI)
      SNF=SIN(FI)
C
C     PLOTTING OF RAYS
C
  26  IF(NRAY.NE.0) READ(lin,101)(NR(I),I = 1,NRAY)
      IF(NRAY.NE.0) WRITE(lou,104)(NR(I),I = 1,NRAY)
      I = 1
      K = 1
  21  READ(LU,101)N,IND
      IPL = 0
      IF(IPRINT.GE.2)WRITE(lou,106)N,IND
      IF(N.EQ.0)GOTO 24
      READ(LU,111)(X(J),Y(J),Z(J),J = 1,N)
      IF(IPRINT.GE.2)WRITE(lou,102)(X(J),Y(J),Z(J),J = 1,N)
      IF(NTICY.EQ.0.and.ntich.eq.0)GO TO 21
      if(nticy.eq.0) go to 41
      IF(NRAY.EQ.0)GOTO 25
      IF(NR(K).EQ.I)GOTO 22
  25  DO 10 J = 1,N
      XX=(X(J)-X0)*CSF+(Y(J)-Y0)*SNF
      YY=Z(J)
      IF(XX.LT.XMIN.OR.XX.GT.XMAX)GOTO 8
      IF(YY.LT.YMIN.OR.YY.GT.YMAX)GOTO 9
      IF(IPL.EQ.1)GOTO 6
      IF(J.LE.IRS)GO TO 10
      IPL = 1
      INDEX = 3
      GO TO 7
  6   INDEX = 2
  7   XNEW = (XX-XMIN)*XMER
      YNEW = (YMAX-(YY-YMIN))*YMER
      CALL PLOT(XNEW,YNEW,INDEX)
      GOTO 10
  8   IF (IPL.EQ.0) GOTO 10
      IF(XX.LT.XMIN) XXX = XMIN
      IF(XX.GT.XMAX) XXX = XMAX
      XOL=(X(J-1)-X0)*CSF+(Y(J-1)-Y0)*SNF
      TG=(YY-Z(J-1))/(XX-XOL)
      YY = TG*(XXX-XX)+YY
      XX = XXX
      IPL=0
      GOTO 6
  9   IF (IPL.EQ.0) GOTO 10
      IF(YY.LT.YMIN) YYY=YMIN
      IF(YY.GT.YMAX) YYY=YMAX
      XOL=(X(J-1)-X0)*CSF+(Y(J-1)-Y0)*SNF
      TG=(XX-XOL)/(YY-Z(J-1))
      XX = TG*(YYY-YY)+XX
      YY = YYY
      IPL=0
      GOTO 6
  10  CONTINUE
c
  41  if(ntich.eq.0)go to 23
      if(nticy.eq.0)shf=0.
      if(nticy.ne.0)shf=ylen+yshift-hmin*hmer
      ipl=0
      DO 42 J = 1,N
      XX=(X(J)-X0)*CSF+(Y(J)-Y0)*SNF
      hh=-(x(j)-x0)*snf+(y(j)-y0)*csf
      IF(XX.LT.XMIN.OR.XX.GT.XMAX)GOTO 45
      IF(hh.LT.hMIN.OR.hh.GT.hMAX)GOTO 46
      IF(IPL.EQ.1)GOTO 43
      IF(J.LE.IRS)GO TO 42
      IPL = 1
      INDEX = 3
      GO TO 44
  43  INDEX = 2
  44  XNEW = (XX-XMIN)*XMER
      hNEW = (hMAX-(hh-hMIN))*hMER
      CALL PLOT(XNEW,hNEW+shf,INDEX)
      GOTO 42
  45  IF (IPL.EQ.0) GOTO 42
      IF(XX.LT.XMIN) XXX = XMIN
      IF(XX.GT.XMAX) XXX = XMAX
      XOL=(X(J-1)-X0)*CSF+(Y(J-1)-Y0)*SNF
      hol=-(x(j-1)-x0)*snf+(y(j-1)-y0)*csf
      TG=(hh-hol)/(XX-XOL)
      hh = TG*(XXX-XX)+hh
      XX = XXX
      IPL=0
      GOTO 43
  46  IF (IPL.EQ.0) GOTO 42
      IF(hh.LT.hMIN) hhh=hMIN
      IF(hh.GT.hMAX) hhh=hMAX
      XOL=(X(J-1)-X0)*CSF+(Y(J-1)-Y0)*SNF
      hol=-(x(j-1)-x0)*snf+(y(j-1)-y0)*csf
      TG=(XX-XOL)/(hh-hol)
      XX = TG*(hhh-hh)+XX
      hh = hhh
      IPL=0
      GOTO 43
  42  CONTINUE
      GOTO 23
  22  IF(K.GE.NRAY)GOTO 23
      K = K+1
  23  I = I+1
      GOTO 21
  24  CONTINUE
C
C     PLOTTING OF TIME-DISTANCE CURVE
C
  11  READ(LU,101)NS
      IF(IPRINT.GE.1)WRITE(lou,101)NS
      NWTYP=NS
      IF(NS.LT.0)NS=-NS
      IF(NS.NE.0)then
        do 300 i=1,ns
        READ(LU,112)INDI(I),X(I),z(i),T(I),AX(1,I),AY(1,I),AZ(1,I)
        if(nwtyp.lt.0)READ(LU,115)AX(2,I),AY(2,I),AZ(2,I)
        read(lu,115)(p(k),k=1,3)
        read(lu,115)(pol(k),k=1,3)
        if(nwtyp.lt.0)read(lu,115)(pol1(k),k=1,3)
  300   continue
      end if
      IF(NS.NE.0.AND.IPRINT.GE.1)WRITE(lou,108)(INDI(I),X(I),z(i),
     1T(I),(AX(j,I),AY(j,I),AZ(j,I),j=1,2),I=1,NS)
      NSS = NS
  35  CONTINUE
      IF(NTICT.EQ.0)GOTO 15
      READ(lin,102)TMIN,TMAX,TLEN,DTICT,VRED
      WRITE(lou,102)TMIN,TMAX,TLEN,DTICT,VRED
      shf=0.
      if(nticy.ne.0)shf=shf+ylen+yshift
      if(ntich.ne.0)shf=shf+hlen+yshift
      IF(IRUN.EQ.1)CALL PLOT(0.,SHF,-3)
      IRUN=1
      IF(NSS.EQ.0)GOTO 30
      IEX=0
      if(iloc.eq.1)then
        xmer=ylen/(ymax-ymin)
        dticx=dticy
        xmin=ymin
        xmax=ymax
        xlen=ylen
      end if
      YMER = TLEN/(TMAX-TMIN)
      IF(ABS(VRED).LT..00001) VRED = 6.
      NAUX=3
      CALL BORDER(XLEN,DTICX,TLEN,DTICT,SC,TEXT,0,XMIN,XMAX,TMIN,
     1TMAX,NTICX,NTICT,NDX,NDY)
      TX=.5*(XLEN-6.3*SC)
      if(iloc.ne.1)
     1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14)
      if(iloc.eq.1)
     1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DEPTH IN KM',0.,11)
      TX=.5*(TLEN-8.1*SC)
      U=-(1.6+.4*NDX)*SC
      IF(IRED.EQ.0)
     1CALL SYMBOL(U,TX,.45*SC,'TRAVEL TIME IN SEC',90.,18)
      IF(IRED.EQ.0)GO TO 27
      CALL SYMBOL(U,TX,.45*SC,'T-D/ ',90.,5)
      TX=TX+1.8*SC
      CALL NUMBER(U,TX,.45*SC,VRED,90.,2)
      TX=TX+2.7*SC
      CALL SYMBOL(U,TX,.45*SC,'(IN SEC)',90.,8)
  27  INDEX=3
  12  DO 13 I = 1,NS
      XNEW = X(I)
      IF(IEX.EQ.1)XNEW=XZ(I)
      IF(XNEW.LT.XMIN.OR.XNEW.GT.XMAX)GOTO 13
      YNEW=T(I)
      IF(IRED.NE.0)YNEW=T(I)-ABS(XNEW)/VRED
      XNEW = (XNEW-XMIN)*XMER
      IF(YNEW.LT.TMIN.OR.YNEW.GT.TMAX)GOTO 13
      YNEW = (YNEW-TMIN)*YMER
      CALL SYMBOL(XNEW,YNEW,.15,char(index),0.,-1)
  13  CONTINUE
      IF(IEX.EQ.0)GOTO 30
      NS = NSS
      GOTO 14
  30  READ(lin,101)NEXP
      WRITE(lou,104)NEXP
      IF(NEXP.EQ.0)GOTO 14
      NS = NEXP
      READ(lin,102)(XZ(I),T(I),I=1,NS)
      WRITE(lou,102)(XZ(I),T(I),I=1,NS)
      IF(NSS.EQ.0)GO TO 15
      IEX = 1
      INDEX=4
      GOTO 12
   14 CONTINUE
      CALL PLOT(0.,-SHF,-3)
C
C     PLOTTING OF AMPLITUDE-DISTANCE CURVE
C
   15 IF(NTICA.EQ.0)GO TO 200
      IRUN1=0
      ALBACK=0.
   19 READ(lin,109)AMIN,AMAX,ALEN,DTICA,ICOMP,msour
      WRITE(lou,109)AMIN,AMAX,ALEN,DTICA,ICOMP,msour
      if(msour.ne.0)call source(lin,lou,0,0,msour,p,pol,amsour,phsour)
      IF(ALEN.LT..00001)CALL PLOT(xshift+xlen,-ALBACK,-3)
      IF(ALEN.LT..00001)GO TO 200
   38 IF(NSS.EQ.0)GO TO 36
      IF(IRUN.EQ.1.AND.IRUN1.EQ.0)CALL PLOT(XLEN+xSHIFT,0.,-3)
      IF(IRUN1.EQ.1)CALL PLOT(0.,ALEN+ySHIFT,-3)
      IF(IRUN1.EQ.1)ALBACK=ALBACK+ALEN+ySHIFT
      IRUN1=1
      IRUN=1
      IEX=0
      YMER = ALEN/(AMAX-AMIN)
      NAUX=2
      CALL BORDER(XLEN,DTICX,ALEN,DTICA,SC,TEXT,0,XMIN,XMAX,AMIN,
     1AMAX,NTICX,NTICA,NDX,NDY)
      TX=.5*(XLEN-6.3*SC)
      if(iloc.ne.1)
     1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14)
      if(iloc.eq.1)
     1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DEPTH IN KM',0.,11)
      TX=.5*(ALEN-7.65*SC)
      U=-(1.6+.4*NDX)*SC
      CALL SYMBOL(U,TX,.45*SC,'A M P L I T U D E',90.,17)
      IF(ICOMP.EQ.0)
     1CALL SYMBOL(.45*SC,ALEN+SC,.45*SC,'VERTICAL',0.,8)
      IF(ICOMP.EQ.1)
     1CALL SYMBOL(.45*SC,ALEN+SC,.45*SC,'RADIAL',0.,6)
      IF(ICOMP.EQ.2)
     1CALL SYMBOL(.45*SC,ALEN+SC,.45*SC,'TRANSVERSE',0.,10)
  28  INDEX=3
  16  DO 17 I=1,NS
      if(msour.ne.0)then
        call source(lin,lou,1,3,msour,p,pol,amsour,phsour)
        csour=amsour*cexp(imag*phsour)
      end if
      if(msour.eq.0)csour=(1.,0.)
      aux(1)=ax(1,i)*csour
      aux(2)=ay(1,i)*csour
      aux(3)=az(1,i)*csour
      if(nwtyp.lt.0)then
        if(msour.ne.0)then
          call source(lin,lou,1,1,msour,p,pol1,amsour,phsour)
          csour=amsour*cexp(imag*phsour)
        end if
        if(msour.eq.0)csour=(1.,0.)
        aux(1)=aux(1)+ax(2,i)*csour
        aux(2)=aux(2)+ay(2,i)*csour
        aux(3)=aux(3)+az(2,i)*csour
      end if
      caux=aux(1)*csrt+aux(2)*snrt
      aux(2)=-aux(1)*snrt+aux(2)*csrt
      aux(1)=caux
      do 37 k=1,3
      rw=real(aux(k))
      cw=aimag(aux(k))
      am(k)=sqrt(rw*rw+cw*cw)
      if(abs(cw).le..000001.and.abs(rw).le.000001)ph(k)=0.
      if(abs(cw).gt..000001.or.abs(rw).gt.000001)ph(k)=atan2(cw,rw)
   37 continue
      if(iprint.eq.2)write(lou,114)(am(k),ph(k),k=1,3)
      XNEW = X(I)
      if(icomp.eq.0)ynew=am(3)
      if(icomp.eq.1)ynew=am(1)
      if(icomp.eq.2)ynew=am(2)
      IF(IEX.EQ.0.OR.NEXP.EQ.0)GO TO 31
      XNEW=XZ(I)
      YNEW=T(I)
   31 IF(XNEW.LE.XMIN.OR.XNEW.GE.XMAX)GO TO 17
      XNEW=(XNEW-XMIN)*XMER
      IF(ABS(YNEW).LT.1E-10)GO TO 17
      YNEW = ALOG10(ABS(YNEW))
      IF(YNEW.LT.AMIN.OR.YNEW.GT.AMAX)GOTO 17
      YNEW = (YNEW-AMIN)*YMER
      CALL SYMBOL(XNEW,YNEW,.15,char(index),0.,-1)
   17 CONTINUE
   36 IF(IEX.EQ.1)GOTO 18
      READ(lin,101)NEXP
      WRITE(lou,104)NEXP
      IF(NEXP.EQ.0)GOTO 18
      NS = NEXP
      READ(lin,102) (XZ(I),T(I),I=1,NS)
      WRITE(lou,102)(XZ(I),T(I),I=1,NS)
      IF(NSS.EQ.0)GO TO 18
      IEX=1
      INDEX=4
      GOTO 16
   18 IF(IEX.EQ.0)GO TO 29
      NS=NSS
   29 CONTINUE
C
C
      GO TO 19
C
C
  100 FORMAT(A)
  101 FORMAT(26I3)
  102 FORMAT(8F10.5)
  103 FORMAT(2X,I2,3(2F7.2,E10.3))
  104 FORMAT(2X,26I3)
  105 FORMAT(5E12.6,I5)
  106 FORMAT(16I5)
  107 FORMAT(2X,'END OF FILE',I3)
  108 FORMAT(I5,3F10.5,6E15.9,25x,6e15.9)
  109 FORMAT(4F10.5,4I5)
  110 FORMAT(F10.5,E15.9)
  111 FORMAT(3E15.5)
  112 FORMAT(I5,9E15.5)
  113 FORMAT(5E15.5)
  114 format(2x,3(e15.5,f10.5))
  115 format(6e15.5)
C
   99 CONTINUE
C
      CALL PLOT(0.,0.,999)
      STOP
      END