C 
C      PROGRAM    S Y N T A N
C      **********************
C
C      PROGRAM SYNTAN IS DESIGNED FOR THE COMPUTATION OF RAY
C      SYNTHETIC SEISMOGRAMS FROM THE IMPULSE SYNTHETIC SEISMOGRAMS
C      STORED IN THE FILE LU2 GENERATED BY THE PROGRAM ANRAY
C
C     **************************************************************
C
      CHARACTER*80 MTEXT,FILEIN,FILEOU,FILE1,FILE2,FILE3
      COMPLEX CAUX,CAX,CAY,CAZ,CAX1,CAY1,CAZ1,AUX,CSOUR,IMAG
      DIMENSION IDIST(2000),TIME(2000),A(2000),PH(2000),
     1SEIS(3001),TAS(2000),DIST(100),SMAX(100),ISEL(100),XX(100),
     2IEP(100),JS(20),AUX(3),p(3),pol(3),pol1(3)
      COMMON/SOUR/ROS
C
C**************************************************
C
      LIN=5
      LOU=6
      LU2=1
      LU4=2
      LU5=3
      FILEIN='syntan.dat'
      FILEOU='syntan.out'
      FILE1='lu2.in'
      FILE2='lu4.out'
      FILE3='lu5.out'
      WRITE(*,'(2A)') ' SPECIFY NAMES OF INPUT AND OUTPUT FILES',
     1' LIN, LOU, LU2, LU4, LU5: '
      READ(*,*) FILEIN,FILEOU,FILE1,FILE2,FILE3
      IF(FILE1.EQ.' ') GO TO 99
      IF(FILE2.EQ.' ') LU4=0
      IF(FILE3.EQ.' ') LU5=0
      OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LOU,FILE=FILEOU,FORM='FORMATTED')
      OPEN(LU2,FILE=FILE1,FORM='FORMATTED',STATUS='OLD')
      IF(LU4.NE.0)OPEN(LU4,FILE=FILE2,FORM='FORMATTED')
      IF(LU5.NE.0)OPEN(LU5,FILE=FILE3,FORM='UNFORMATTED')
C
C**************************************************
C
      IMAG=(0.,1.)
      PI=3.14159265
      SHF=0.
    1 READ(LIN,*)MPRINT
      WRITE(LOU,100)MPRINT
C
C     ***
C
      REWIND LU2
      IF(LU4.NE.0)REWIND LU5
      IF(LU4.NE.0)REWIND LU4
    2 READ(LU2,101)MTEXT
      READ(LU2,100)NDST,KSH,ILOC
      READ(LU2,102)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,ROS
      IF(NDST.EQ.1)THEN
        READ(LU2,102)XREC,YREC,XPRF,XATAN
      ELSE
        READ(LU2,102)(DIST(K),K=1,NDST)
      END IF
C
      WRITE(LOU,101)MTEXT
      WRITE(LOU,100)NDST,KSH,ILOC
      WRITE(LOU,102)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,ROS
      IF(NDST.EQ.1)WRITE(LOU,102)XREC,YREC,XPRF,XATAN
      IF(NDST.NE.1)WRITE(LOU,102)(DIST(K),K=1,NDST)
      IF(NDST.EQ.1)DIST(1)=XPRF
      READ(LIN,*)MCOMP,MRED,MSELEC,MEPIC,MSHIFT,KABS,MSOUR
      WRITE(LOU,100)MCOMP,MRED,MSELEC,MEPIC,MSHIFT,KABS,MSOUR
      READ(LIN,*)TMIN,DT,TMAX,FREQ,GAMA,PSI,VRED,SHIFT,AROT
      WRITE(LOU,102)TMIN,DT,TMAX,FREQ,GAMA,PSI,VRED,SHIFT,AROT
      CSRT=COS(AROT)
      SNRT=SIN(AROT)
      IF(MEPIC.EQ.0)GO TO 11
      READ(LIN,*)NEPIC,(IEP(I),I=1,NEPIC)
      WRITE(LOU,100)NEPIC,(IEP(I),I=1,NEPIC)
   11 IF(MSELEC.EQ.0)GO TO 12
      READ(LIN,*)NSELEC,(ISEL(I),I=1,NSELEC)
      WRITE(LOU,100)NSELEC,(ISEL(I),I=1,NSELEC)
   12 CONTINUE
      IF(KABS.EQ.0)GO TO 5
      READ(LIN,*) FREF,QRED
      WRITE(LOU,102)FREF,QRED
      IF(FREF.LT..00001)FREF=1.
      IF(QRED.LT..00001)QRED=1.
      RF=FREQ/FREF
C
    5 IF(MSOUR.NE.0)CALL SOURCE(LIN,LOU,0,0,MSOUR,P,POL,AMSOUR,PHSOUR)
C
      NCD=0
      I=0
    3 I=I+1
      IF(I.GT.2000)WRITE(LOU,104)
      IF(I.GT.2000)GO TO 99
      READ(LU2,116,END=7)NC,IDIST(I),TIME(I),CAX,CAY,CAZ,TAST
      IF(NC.LT.0)READ(LU2,117)CAX1,CAY1,CAZ1
      READ(LU2,118)(P(K),K=1,3)
      READ(LU2,118)(POL(K),K=1,3)
      IF(NC.LT.0)READ(LU2,118)(POL1(K),K=1,3)
      CAUX=CAX*CSRT+CAY*SNRT
      CAY=-CAX*SNRT+CAY*CSRT
      CAX=CAUX
      CAUX=CAX1*CSRT+CAY1*SNRT
      CAY1=-CAX1*SNRT+CAY1*CSRT
      CAX1=CAUX
      NC1=IABS(NC)
      NC1=NC1+NCD
      IF(MSELEC.EQ.0)GO TO 14
      DO 13 J=1,NSELEC
      IF(ISEL(J).EQ.NC1)GO TO 14
   13 CONTINUE
      GO TO 3
   14 IF(MEPIC.EQ.0)GO TO 16
      DO 15 J=1,NEPIC
      IF(IEP(J).EQ.IDIST(I))GO TO 16
   15 CONTINUE
      GO TO 3
   16 CONTINUE
C
      IF(MSOUR.NE.0)THEN
        CALL SOURCE(LIN,LOU,1,1,MSOUR,P,POL,AMDIR,PHDIR)
        CSOUR=AMDIR*CEXP(IMAG*PHDIR)
      END IF
      IF(MSOUR.EQ.0)CSOUR=(1.,0.)
      AUX(1)=CAX*CSOUR
      AUX(2)=CAY*CSOUR
      AUX(3)=CAZ*CSOUR
      IF(NC.LT.0)THEN
        IF(MSOUR.NE.0)THEN
          CALL SOURCE(LIN,LOU,1,1,MSOUR,P,POL1,AMDIR,PHDIR)
          CSOUR=AMDIR*CEXP(IMAG*PHDIR)
        END IF
        IF(MSOUR.EQ.0)CSOUR=(1.,0.)
        AUX(1)=AUX(1)+CAX1*CSOUR
        AUX(2)=AUX(2)+CAY1*CSOUR
        AUX(3)=AUX(3)+CAZ1*CSOUR
      END IF
      IF(MCOMP.EQ.0)K=3
      IF(MCOMP.EQ.1)K=1
      IF(MCOMP.EQ.2)K=2
      RW=REAL(AUX(K))
      CW=AIMAG(AUX(K))
      A(I)=SQRT(RW*RW+CW*CW)
      IF(ABS(RW).LT..000001.AND.ABS(CW).LT..000001)THEN
        PH(I)=0.
      ELSE
        PH(I)=ATAN2(CW,RW)
      END IF
      IF(MPRINT.GE.1)WRITE(LOU,114)
     1NC,IDIST(I),TIME(I),A(I),PH(I),P,TAST
      IF(KABS.NE.0)TAS(I)=TAST
      IF(KABS.EQ.2)PH(I)=PH(I)-2.*FREQ*TAST*ALOG(RF)
      GO TO 3
C
    7 ISUM=I-1
      IF(VRED.LT.0.0001)VRED=8.
      IF(FREQ.LT.0.0001)FREQ=4.
      IF(GAMA.LT.0.0001)GAMA=4.0
      IF(MSHIFT.EQ.0)SHF=0.
      IF(MSHIFT.EQ.1)SHF=.241506*GAMA/FREQ
      IF(MSHIFT.EQ.2)SHF=SHIFT
      OMEGA=2.*PI*FREQ
      OMRF=2.*PI*FREF
      OG=OMEGA/GAMA
      OGG=OG*OG
      TSH=0.45*GAMA/FREQ
C
C     CONSTRUCTION OF SYNTHETIC SEISMOGRAMS
C
      IF(NDST.GT.100)WRITE(LOU,105)
      IF(NDST.GT.100)GO TO 99
      NPTS=(TMAX-TMIN)/DT+1.5
      IF(NPTS.GT.3001)WRITE(LOU,106)
      IF(NPTS.GT.3001)GO TO 99
C
C
      MMD=NDST
      IF(MEPIC.NE.0)MMD=NEPIC
      IF(LU4.EQ.0)GO TO 6
      WRITE(LU5)NPTS,TMIN,DT,NDST,DIST(1),RSTEP
      WRITE(LU4,101)MTEXT
      WRITE(LU4,108)MMD,MRED,MCOMP,ILOC,VRED,RSTEP,XSOUR,DT
C
C     LOOP FOR RECEIVER POSITIONS
C
    6 CONTINUE
      DO 30 I=1,MMD
      JJ=I
      IF(MEPIC.NE.0)JJ=IEP(I)
      XX(I)=DIST(JJ)
      XXD=ABS(XX(I)-XSOUR)
      XXV=0.
      IF(MRED.NE.0)XXV=XXD/VRED
      AMAX=0.
      DO 20 J1=1,NPTS
   20 SEIS(J1)=0.
      DO 21 K=1,ISUM
      IF(IDIST(K).NE.JJ)GO TO 21
      TR=TIME(K)-XXV+SHF
      IF((TR+TSH).LT.TMIN)GO TO 21
      IF((TR-TSH).GT.TMAX)GO TO 21
      AMPL=A(K)
      PHASE=PH(K)
      IF(AMPL.LT.0.00005*AMAX)GO TO 21
      IF(AMPL.GT.AMAX)AMAX=AMPL
      IF1=IFIX((TR-TSH-TMIN)/DT+0.1)
      IF2=IFIX((TR+TSH-TMIN)/DT+0.1)
      IFF1=1
      IF(IF1.GT.0)IFF1=IF1
      IF1=IFF1
      IFF2=NPTS
      IF(IF2.LT.NPTS)IFF2=IF2
      IF2=IFF2
      IF(KABS.EQ.0)GO TO 27
      TAST=TAS(K)*QRED
      OMAS=OMEGA*(1.-OMEGA*TAST/(GAMA*GAMA))
      IF(OMAS.LE.0.)WRITE(LOU,112)
      IF(OMAS.LE.0.)GO TO 32
      AABS1=OMEGA*TAST/GAMA
      AABS1=.25*(AABS1*AABS1)
      AABS2=.5*OMAS*TAST
      IF(KABS.NE.2)GO TO 27
      AABS4=OMAS/OMRF
      AABS5=ALOG(AABS4)
      AABS3=(TAST/PI)*(1.+AABS5)
   27 DO 24 KK=IF1,IF2
      TT=TMIN+DT*FLOAT(KK-1)
      TD=TT-TR
      IF(KABS.EQ.2)TD=TD+AABS3
      AA=PSI-PHASE
      BB=-OGG*TD*TD
      IF(KABS.EQ.0)AA=AA+OMEGA*TD
      IF(KABS.NE.0)AA=AA+OMAS*TD
      IF(KABS.EQ.2)AA=AA-OMAS*TAST/PI
      IF(KABS.NE.0)BB=BB-AABS1-AABS2
      AA=AMPL*COS(AA)*EXP(BB)
   24 SEIS(KK)=SEIS(KK)+AA
   21 CONTINUE
      SMAX(I)=0.
      DO 25 KK=1,NPTS
      IF(SMAX(I).GT.ABS(SEIS(KK)))GO TO 25
      SMAX(I)=ABS(SEIS(KK))
   25 CONTINUE
C
C
      IF(LU4.NE.0)WRITE(LU5)(SEIS(J),J=1,NPTS)
   30 CONTINUE
      REWIND LU2
      IF(LU4.NE.0)REWIND LU5
C
C     END OF THE LOOP FOR RECEIVERS
C
      SMAXIM=0.
      DO 26 I=1,MMD
      IF(SMAXIM.GT.SMAX(I))GO TO 26
      SMAXIM=SMAX(I)
      XXX=XX(I)
   26 CONTINUE
      IF(LU4.NE.0)WRITE(LU4,113)XXX,SMAXIM
      WRITE(LOU,113)XXX,SMAXIM
      IF(LU4.EQ.0)GO TO 1
C
      IF(LU4.NE.0)READ(LU5)NPTS,TMIN,DT,NDST,DIST(1),RSTEP
      DO 31 I=1,MMD
      IF(LU4.NE.0)READ(LU5)(SEIS(J),J=1,NPTS)
C
C     NORMALIZATION OF SEISMOGRAMS
C
      IF1=0
      IF2=0
      SMAXI=SMAX(I)
      IF(SMAXI.LT..000001)GO TO 35
      DO 34 J=1,NPTS
      SEIS(J)=999.1*SEIS(J)/SMAXI
      IF(ABS(SEIS(J)).LT.1.)GO TO 33
      IF2=0
      IF(IF1.EQ.0.AND.J.EQ.1)IF1=1
      IF(IF1.EQ.0)IF1=J-1
      GO TO 34
   33 IF(IF1.EQ.0)GO TO 34
      IF(IF2.EQ.0)IF2=J
   34 CONTINUE
      TM=TMIN+FLOAT(IF1-1)*DT
      IF(IF2.EQ.0)IF2=NPTS
      NPS=IF2-IF1+1
   35 IF(IF1.EQ.0)NPS=0
      IF(IF1.EQ.0)TM=TMIN
C
      WRITE(LU4,109)XX(I),SMAX(I),TM,NPS
      WRITE(LOU,107)XX(I),SMAX(I),TM,NPS
      ISS=20
      IS=IF1-1
   41 IS1=IF2-IS
      IF(IS1.LT.20)ISS=IS1
      IF(IS.LT.0)IS=0
      DO 40 K=1,ISS
   40 JS(K)=SEIS(IS+K)
      WRITE(LU4,111)(JS(K),K=1,ISS)
      IF(MPRINT.EQ.2.AND.NPS.NE.0)
     1WRITE(LOU,110)(JS(K),K=1,ISS)
      IF(IS1.LE.20)GO TO 31
      IS=IS+20
      GO TO 41
   31 CONTINUE
      REWIND LU4
   32 CONTINUE
C
C
  100 FORMAT(26I3)
  101 FORMAT(A)
  102 FORMAT(8F10.5)
  104 FORMAT(1X,'NUMBER OF READINGS FROM LU2 GREATER THAN 2000.'/,
     11X,'COMPUTATION TERMINATED'/)
  105 FORMAT(1X,'NUMBER OF RECEIVER POSITIONS GREATER THAN 100.',/,
     11X,'COMPUTATION TERMINATED'/)
  106 FORMAT(1X,'NUMBER OF POINTS IN THE SYNTHETIC SEISMOGRAM
     1GREATER THAN 3001.'/1X,'COMPUTATION TERMINATED'/)
  107 FORMAT(1X,'EPICENTRAL DISTANCE =',F10.5,E15.9,1X,F10.5,I5)
  108 FORMAT(4I5,4F10.5)
  109 FORMAT(F10.5,E15.8,F10.5,I5)
  110 FORMAT(1X,20I4)
  111 FORMAT(20I4)
  112 FORMAT(1X,'FREQUENCY TOO HIGH OR GAMA TOO SMALL'/
     11X,'COMPUTATION TERMINATED'/)
  113 FORMAT(1X,'EPICENTRAL DISTANCE =',F10.5,1X,'SMAXIM =',E15.9)
  114 FORMAT(2I3,F10.5,E12.6,5F10.6)
  116 FORMAT(2I3,F12.7,6E12.6,F10.6)
  117 FORMAT(6E12.6)
  118 FORMAT(3F10.5)
C
   99 CONTINUE
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'source.for'
C     source.for
C
C=======================================================================
C