C 
C     PROGRAM SYNT - SYNTHETIC SEISMOGRAM COMPUTATION
C
C     **************************************************************
C
      EXTERNAL SIGNAL
C
      CHARACTER*80 FILEIN,FILEOU,FILE1,FILE2
      CHARACTER*80 TITLE,MPRINT,IPRINT,PSTEXT
      DIMENSION NDIS(100)
      COMPLEX S(2048),REF(2048),SS(2048)
      COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN
      COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3,
     1NPL4,IPR1,NDER,NINT
      COMMON/WIND/FLO,FLEFT,FRO,FRIGHT,FEXP,ILF,IRF,NFW
      COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4)
      COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT(
     14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE
     2,NPLOT,MCOMP
C
C
C**************************************************
C
      LIN=5
      LOU=6
      LU7=1
      LU8=2
      FILEIN='synfan.dat'
      FILEOU='synfan.out'
      FILE1='lu7.dat'
      FILE2='lu8.dat'
      WRITE(*,'(2A)') ' (SYNFAN) SPECIFY NAMES OF INPUT AND OUTPUT', 
     1' FILES LIN, LOU, LU7, LU8: '
      READ(*,*) FILEIN,FILEOU,FILE1,FILE2
      IF(FILE1.EQ.' ') GO TO 99
      IF(FILE2.EQ.' ') LU8=0
      OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LOU,FILE=FILEOU,FORM='FORMATTED')
      OPEN(LU7,FILE=FILE1,FORM='FORMATTED',STATUS='OLD')
      IF(LU8.NE.0)OPEN(LU8,FILE=FILE2,FORM='FORMATTED')
C
C**************************************************
C
C
      TITLE='SYNFAN'
      PSTEXT=' '
      IPR1=0
      IPR2=0
      READ(LIN,*)TITLE
      READ(LIN,*)IPR1,IPR2,PSTEXT
      WRITE(LOU,111)TITLE
      WRITE(LOU,103)IPR1,IPR2,PSTEXT
      REWIND LU7
      IF(LU8.NE.0)REWIND LU8
C
C     READ INPUT PARAMETERS
C
      NSIG=0
      NPTS=400
      NT=1024
      NWIN=0
      NFILT=0
      NPLOT=0
      NNPLOT=0
      NSTOP=0
      NSH=7
      NDER=0
      NINT=0
      READ(LIN,*)NSIG,NPTS,NT,NWIN,NFILT,NPLOT,NSTOP,NSH,NDER,NINT
      WRITE(LOU,101)NSIG,NPTS,NT,NWIN,NFILT,NPLOT,NSTOP,NSH,
     1NDER,NINT
      IF(NPLOT.NE.0)THEN
        IF(IPR1.LT.0)THEN
          CALL PLOTN(PSTEXT,0)
          IPR1=-IPR1
        END IF
        CALL PLOTS(LDUM1,LDUM2,7)
        READ(LIN,*)NNPLOT,(NDIS(I),I=1,NNPLOT)
        WRITE(LOU,101)NNPLOT,(NDIS(I),I=1,NNPLOT)
      END IF
      NT2=NT/2
      NTM=NT2+1
      READ(LU7,110)MPRINT
      READ(LU7,110)IPRINT
      READ(LU7,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,FL,DF
      READ(LU7,102)NDST,NFS,ILS,MCOMP,ILOC
      DT=1./(DF*FLOAT(NT))
      WRITE(LOU,111)MPRINT
      WRITE(LOU,111)IPRINT
      WRITE(LOU,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,FL,DF
      WRITE(LOU,101)NDST,NFS,ILS,MCOMP
      WRITE(LU8,110)MPRINT
      WRITE(LU8,110)IPRINT
      WRITE(LU8,110)TITLE
      WRITE(LU8,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,DT,DF
      WRITE(LU8,101)NDST,NT,MCOMP,ILOC
      IRS=ILS+NFS
      IRSS=IRS+1
      NI=10
      IF(NT.EQ.2048)NI=11
      IF(NT.EQ.512)NI=9
      IF(NT.EQ.256)NI=8
      IF(NT.EQ.128)NI=7
      READ(LIN,*)ARESP,JRESP
      READ(LIN,*)ASPECT,JSPECT
      READ(LIN,*)ASYNT,JSYNT
      WRITE(LOU,106)ARESP,JRESP
      WRITE(LOU,106)ASPECT,JSPECT
      WRITE(LOU,106)ASYNT,JSYNT
      IF(NPLOT.EQ.0)GO TO 7
      CALL PLOT(3.,5.,-3)
C
C     INPUT SIGNAL
C
    7 CALL SIGNAL(LIN,LOU)
C
C     SPECTRUM OF THE INPUT SIGNAL
C
      CALL FT(LIN,LOU)
      IF(NSTOP.EQ.1.AND.NPLOT.NE.0)CALL PLOT(0.,0.,999)
      IF(NSTOP.EQ.1)STOP
C
C     LOOP FOR RECEIVERS
C
      DO 40 IDIST=1,NDST
      NPR1=0
      NPR2=0
      NPR3=0
      READ(LU7,160)DST,AA
      WRITE(LOU,160)DST,AA
      NFF=2*NFS
      DO 49 I=1,NTM
      FF(I)=0.
      SS(I)=(0.,0.)
   49 REF(I)=(0.,0.)
      READ(LU7,161)(IS(I),I=1,NFF)
      DO 41 I=1,NFS
      A=AA*FLOAT(IS(2*(I-1)+1))/99999.1
      B=AA*FLOAT(IS(2*I))/99999.1
      FF(I+ILS)=SQRT(A*A+B*B)
   41 REF(I+ILS)=CMPLX(A,B)
      IF(NNPLOT.EQ.0)GO TO 42
      DO 43 I=1,NNPLOT
      IF(IDIST.EQ.NDIS(I))GO TO 45
   43 CONTINUE
      GO TO 42
   45 CONTINUE
      READ(LIN,*)NPR1,NPR2,NPR3,TMIN,TMAX,TLEN
      WRITE(LOU,170)NPR1,NPR2,NPR3,TMIN,TMAX,TLEN
      ASYNT(1)=TMIN
      ASYNT(2)=TMAX
      ASYNT(3)=TLEN
      CALL PLOT(3.,0.,-3)
      CALL SYMBOL(-0.25*TLEN,0.,0.05*TLEN,'X(KM)=',90.,6)
      CALL NUMBER(-0.25*TLEN,0.3*TLEN,0.05*TLEN,DST,90.,2)
   42 CONTINUE
C
C     PRINTING AND PLOTTING OF FREQUENCY RESPONSE
C
      IF(IPR2.LE.1.AND.NPR1.EQ.0)GO TO 50
      CALL REDUC(AREDUC,1,IRSS,FF,GG,IS)
      IF(IPR2.LE.1)GO TO 48
      WRITE(LOU,114)DST,DF,AREDUC,IRSS
      WRITE(LOU,107)(IS(I),I=1,IRSS)
   48 IF(NPR1.EQ.0)GO TO 50
      NOO=1
      NRED=NPR1
      IF(NPR1.EQ.1)CALL ZPLOT(0.,DF,ILS,IRSS,GG)
      IF(NPR1.EQ.2)CALL ZPLOT(0.,DF,ILS,IRSS,FF)
   50 CONTINUE
C
C     AMPLITUDE SPECTRUM OF SYNTHETIC SEISMOGRAM
C
      DO 18 I=1,NTM
   18 FF(I)=0.
      DO 11 I=ILF,IRF
      SS(I)=REF(I)*S(I)
      A=REAL(SS(I))
      B=AIMAG(SS(I))
   11 FF(I)=SQRT(A*A+B*B)
      DO 12 I=2,NT2
      IJ=NT-I+2
   12 SS(IJ)=CONJG(SS(I))
      IF(IPR2.LE.2.AND.NPR2.EQ.0)GO TO 60
      CALL REDUC(AREDUC,1,IRSS,FF,GG,IS)
      IF(IPR2.LE.2)GO TO 52
      WRITE(LOU,115)DST,DF,AREDUC,IRSS
      WRITE(LOU,107)(IS(I),I=1,IRSS)
   52 IF(NPR2.EQ.0)GO TO 60
      NOO=5
      NRED=NPR2
      IF(NPR2.EQ.1)CALL ZPLOT(0.,DF,ILS,IRSS,GG)
      IF(NPR2.EQ.2)CALL ZPLOT(0.,DF,ILS,IRSS,FF)
   60 CONTINUE
C
C     GOING BACK TO TIME DOMAIN. SYNTHETIC SEISMOGRAMS
C
      CALL FCOOLR(NI,SS,-1.0)
      DEL=1.0/FLOAT(NT)
      DO 20 I=1,NT
   20 FF(I)=DEL*REAL(SS(I))
      CALL REDUC(AREDUC,1,NT,FF,GG,IS)
      WRITE(LU8,162)DST,TO,AREDUC,NT
      WRITE(LU8,107)(IS(I),I=1,NT)
C
      IF(IPR2.EQ.0.AND.NPR3.EQ.0)GO TO 70
      IF(IPR2.EQ.0)GO TO 65
      WRITE(LOU,116)DST,TO,DT,AREDUC,NT
      WRITE(LOU,108)(IS(I),I=1,NT)
   65 IF(NPR3.EQ.0)GO TO 70
      NOO=6
      NRED=NPR3
      IF(AREDUC.EQ.0.)GO TO 66
      TSTART=TMIN-TO
      IPOM=IFIX(TSTART/DT)+1
   67 IF(IPOM.GT.0)GO TO 61
      IPOM=IPOM+NT
      GO TO 67
   61 IF(IPOM.LE.NT)GO TO 62
      IPOM=IPOM-NT
      GO TO 61
   62 CONTINUE
      IIAUX=NT-IPOM
      DO 63 J=1,NT
      IF(J.GE.IPOM)JAUX=J-IPOM+1
      IF(J.LT.IPOM)JAUX=J+IIAUX+1
   63 FF(JAUX)=GG(J)*AREDUC
      TLEN=ASYNT(2)-ASYNT(1)
      NFIN=IFIX(TLEN/DT)+1
      CALL REDUC(AREDUC,1,NFIN,FF,GG,IS)
   66 CONTINUE
      IF(NPR3.EQ.1)CALL ZPLOT(TMIN,DT,1,NFIN,GG)
      IF(NPR3.EQ.2)CALL ZPLOT(TMIN,DT,1,NFIN,FF)
   70 CONTINUE
   40 CONTINUE
C
C     END OF THE LOOP FOR RECEIVERS
C
   99 CONTINUE
      IF(NPLOT.NE.0)CALL PLOT(0.,0.,999)
C
  101 FORMAT(16I5)
  102 FORMAT(26I3)
  103 FORMAT(2I5,1X,A)
  106 FORMAT(6F10.5,4I5)
  107 FORMAT(20I4)
  108 FORMAT(1X,20I4)
  110 FORMAT(A)
  111 FORMAT(1X,A)
  114 FORMAT(/1X,'FREQUENCY RESPONSE-AMPLITUDES   X(KM)=',F10.5/,2X,
     1'DF=',F10.5,2X,'MAX=',E16.7,2X,'N=',I5)
  115 FORMAT(/1X,'AMPLITUDE SPECTRUM OF SYNTHETIC SEISMOGRAM   X(KM)=',
     1F10.5/2X,'DF=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5)
  116 FORMAT(/1X,'SYNTHETIC SEISMOGRAM   X(KM)=',F10.5/2X,'T0=',
     1F10.5,2X,'DT=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5)
  152 FORMAT(5F10.5,2E15.7)
  154 FORMAT(3E16.8)
  160 FORMAT(F10.3,E12.5)
  161 FORMAT(12I6)
  162 FORMAT(2F10.3,1E12.5,I5)
  170 FORMAT(3I5,3F10.5)
      REWIND LU7
      REWIND LU8
      STOP
      END
C
C
      SUBROUTINE BORD
      CHARACTER*80 TITLE
      COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT(
     14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE
     2,NPLOT,MCOMP
      IF(NOO.GT.1)GO TO 1
      RMIN=ARESP(1)
      RMAX=ARESP(2)
      GLR=ARESP(3)
      TRMIN=ARESP(4)
      TRMAX=ARESP(5)
      GLTR=ARESP(6)
      IR=JRESP(1)
      ITR=JRESP(2)
      NDECX=JRESP(3)
      NDECY=JRESP(4)
      GO TO 3
    1 IF(NOO.EQ.2.OR.NOO.GE.6)GO TO 2
      RMIN=ASPECT(1)
      RMAX=ASPECT(2)
      GLR=ASPECT(3)
      TRMIN=ASPECT(4)
      TRMAX=ASPECT(5)
      GLTR=ASPECT(6)
      IR=JSPECT(1)
      ITR=JSPECT(2)
      NDECX=JSPECT(3)
      NDECY=JSPECT(4)
      GO TO 3
    2 RMIN=ASYNT(1)
      RMAX=ASYNT(2)
      GLR=ASYNT(3)
      TRMIN=ASYNT(4)
      TRMAX=ASYNT(5)
      GLTR=ASYNT(6)
      IR=JSYNT(1)
      ITR=JSYNT(2)
      NDECX=JSYNT(3)
      NDECY=JSYNT(4)
    3 CONTINUE
      CALL SYMBOL(-0.12*GLR,0.,.07*GLR,TITLE,90.,80)
      CALL PLOT(0.,0.,3)
      IF(TRMAX.GE.0.)
     1CALL NUMBER(0.025*GLR,-0.2*GLR,0.05*GLR,TRMAX,90.,NDECY)
      IF(TRMAX.LT.0.)
     1CALL NUMBER(0.025*GLR,-0.25*GLR,0.05*GLR,TRMAX,90.,NDECY)
      CALL PLOT(0.,0.,3)
      ZR=(RMAX-RMIN)/FLOAT(IR)
      ZTR=(TRMAX-TRMIN)/FLOAT(ITR)
      A=GLR/FLOAT(IR)
      B=GLTR/FLOAT(ITR)
      DO 100 J=1,ITR
      F=B*FLOAT(J)
      TG=TRMAX-ZTR*FLOAT(J)
      CALL PLOT(F,0.,2)
      CALL PLOT(F,0.2,2)
      IF(TG.GE.0.)
     1CALL NUMBER(F+0.025*GLR,-0.2*GLR,0.05*GLR,TG,90.,NDECY)
      IF(TG.LT.0.)
     1CALL NUMBER(F+0.025*GLR,-0.25*GLR,0.05*GLR,TG,90.,NDECY)
      CALL PLOT(F,0.,3)
  100 CONTINUE
      CALL NUMBER(F+0.1*GLR,-.025*GLR,0.05*GLR,RMIN,90.,NDECX)
      IF(NOO.EQ.1)
     1CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR,'FREQUENCY RESPONSE',90.,18)
      IF(NOO.EQ.2)
     1CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR,'INPUT SIGNAL',90.,12)
      IF(NOO.EQ.3)CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR,
     1'AMPLITUDE SPECTRUM OF INPUT SIGNAL',90.,34)
      IF(NOO.EQ.4)CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR,
     1'AMPLITUDE SPECTRUM, WINDOW APPLIED',90.,34)
      IF(NOO.EQ.5)
     1CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR,'AMPLITUDE SPECTRUM',90.,18)
      IF(NOO.EQ.5)CALL SYMBOL(F+0.35*GLR,0.,0.05*GLR,
     1'OF SYNTHETIC SEISMOGRAM',90.,23)
      IF(NOO.EQ.6)THEN
        CALL SYMBOL
     1  (F+0.25*GLR,0.,0.05*GLR,'SYNTHETIC SEISMOGRAM',90.,20)
        IF(MCOMP.EQ.0)
     1  CALL SYMBOL(F+0.35*GLR,0.,0.05*GLR,'VERTICAL',90.,8)
        IF(MCOMP.EQ.1)
     1  CALL SYMBOL(F+0.35*GLR,0.,0.05*GLR,'X-COMPONENT',90.,11)
        IF(MCOMP.EQ.2)
     1  CALL SYMBOL(F+0.35*GLR,0.,0.05*GLR,'Y-COMPONENT',90.,11)
      END IF
      IF(NOO.EQ.7)CALL SYMBOL(F+0.25*GLR,0.,0.05*GLR,
     1'INPUT SIGNAL AFTER WINDOWING',90.,28)
      IF(NOO.EQ.7)NOO=2
      CALL PLOT(F,0.,3)
      DO 200 J=1,IR
      C=A*FLOAT(J)
      CALL PLOT(F,C,2)
      CALL PLOT(F-0.2,C,2)
      ZG=RMIN+ZR*FLOAT(J)
      CALL NUMBER(F+0.1*GLR,C-.08*GLR,0.05*GLR,ZG,90.,NDECX)
      CALL PLOT(F,C,3)
  200 CONTINUE
      IF(NOO.NE.2.AND.NOO.NE.6)
     1CALL SYMBOL(F+0.18*GLR,C-0.25*GLR,0.04*GLR,'F(HZ)',90.,5)
      IF(NOO.EQ.2.OR.NOO.EQ.6)
     1CALL SYMBOL(F+0.18*GLR,C-0.3*GLR,0.04*GLR,'TIME(S)',90.,7)
      CALL PLOT(F,C,3)
      DO 300 J=1,ITR
      D=B*FLOAT(ITR-J)
      CALL PLOT(D,C,2)
      CALL PLOT(D,C-0.2,2)
  300 CALL PLOT(D,C,2)
      IF(NRED.EQ.1)
     1CALL SYMBOL(D-0.05*GLR,0.,0.03*GLR,'REDUCTION FACTOR= ',90.,18)
      IF(NRED.EQ.1)
     1CALL NUMBER(D-0.05*GLR,0.7*GLR,0.03*GLR,AREDUC,90.,5)
      CALL PLOT(D,C,3)
      DO 400 J=1,IR
      E=A*FLOAT(IR-J)
      CALL PLOT(D,E,2)
      CALL PLOT(D+0.2,E,2)
  400 CALL PLOT(D,E,2)
      RETURN
      END
C
C
      SUBROUTINE ZPLOT(FMIN,DF,NMIN,NMAX,F)
      CHARACTER*80 TITLE
      DIMENSION F(1)
      COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT(
     14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE
     2,NPLOT,MCOMP
      CALL BORD
      GMR=GLR/(RMAX-RMIN)
      GMTR=GLTR/(TRMAX-TRMIN)
C
      DO 1 II=NMIN,NMAX
      VO=GMR*(FMIN-RMIN+DF*FLOAT(II-1))
      IF((VO+0.0001).GE.0.)GO TO 2
    1 CONTINUE
      GO TO 99
C
    2 IMAX1=0
      WO=GMTR*(TRMAX-F(II))
      VOS=VO
      WOS=WO
      IF(WO.GT.0.)GO TO 3
      IMAX1=1
      CALL PLOT(0.,VO,3)
      GO TO 4
    3 CALL PLOT(WO,VO,3)
    4 CONTINUE
      II=II+1
C
      DO 10 I=II,NMAX
      VO=GMR*(FMIN-RMIN+DF*FLOAT(I-1))
      WO=GMTR*(TRMAX-F(I))
      IF(VO.GT.(GLR+0.0001))GO TO 11
      IF(WO.LT.0.)GO TO 12
      IF(IMAX1.EQ.0)GO TO 13
      VV=VOS-WOS*(VO-VOS)/(WO-WOS)
      CALL PLOT(0.,VV,2)
      CALL PLOT(WO,VO,2)
      IMAX1=0
      GO TO 9
   13 CALL PLOT(WO,VO,2)
      IMAX1=0
      GO TO 9
   12 IF(IMAX1.EQ.0)GO TO 15
      CALL PLOT(0.,VO,2)
      IMAX1=1
      GO TO 9
   15 VV=VOS-WOS*(VO-VOS)/(WO-WOS)
      CALL PLOT(0.,VV,2)
      CALL PLOT(0.,VO,2)
      IMAX1=1
      GO TO 9
   11 IF(IMAX1.EQ.0)GO TO 16
      CALL PLOT(0.,RMAX,2)
      GO TO 99
   16 IF(WO.LT.0.)GO TO 18
   17 WW=WOS+(WO-WOS)*(GLR-VOS)/(VO-VOS)
      CALL PLOT(WW,GLR,2)
      GO TO 99
   18 VV=VOS-WOS*(VO-VOS)/(WO-WOS)
      IF(VV.GT.GLR)GO TO 17
      CALL PLOT(0.,VV,2)
      CALL PLOT(0.,GLR,2)
      GO TO 99
    9 VOS=VO
      WOS=WO
   10 CONTINUE
      CALL PLOT(GMTR*TRMAX,VOS,2)
      CALL PLOT(GMTR*TRMAX,GLR,2)
C
   99 CONTINUE
      CALL PLOT(GLTR+FLOAT(NSH),0.,-3)
      RETURN
      END
C
C
      SUBROUTINE SIGNAL(LIN,LOU)
      CHARACTER*80 TITLE
      COMPLEX S(2048),REF(2048)
      COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3,
     1NPL4,IPR,NDER,NINT
      COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN
      COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4)
      COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT(
     14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE
     2,NPLOT,MCOMP
C
      NPL1=0
      NPL2=0
      NPL3=0
      NPL4=0
      IF(NPLOT.EQ.0)GO TO 3
      READ(LIN,*) NPL1,NPL2,NPL3,NPL4
      WRITE(LOU,100)NPL1,NPL2,NPL3,NPL4
    3 READ(LIN,*)IAUX,AUX,TO
      WRITE(LOU,104)IAUX,AUX,TO
      DO 1 I=1,NPTS
      F(I)=0.
    1 GG(I)=0.
C
      IF(NSIG)30,20,30
   30 READ(LIN,*)(IS(I),I=1,NPTS)
      WRITE(LOU,120)NSIG
      WRITE(LOU,102)(IS(I),I=1,NPTS)
      A=10.**(-NSIG)
      DO 32 I=1,NPTS
      F(I)=A*FLOAT(IS(I))
      IS(I)=0
   32 CONTINUE
      GO TO 60
   20 CALL SFUN
   60 CONTINUE
C
      CALL REDUC(AREDUC,1,NPTS,F,GG,IS)
      IF(IPR.EQ.0)GO TO 61
      WRITE(LOU,121)TO,DT,AREDUC,NPTS
      WRITE(LOU,107)(IS(I),I=1,NPTS)
   61 IF(NPL1.EQ.0)GO TO 10
      NOO=2
      NRED=NPL1
      IF(NPL1.EQ.1)CALL ZPLOT(TO,DT,1,NPTS,GG)
      IF(NPL1.EQ.2)CALL ZPLOT(TO,DT,1,NPTS,F)
C
   10 CONTINUE
      RETURN
  100 FORMAT(16I5)
  102 FORMAT(2X,16I5)
  104 FORMAT(4I5,5F10.5)
  107 FORMAT(20I4)
  120 FORMAT(2X,'SIGNAL.   NSIG=',I3)
  121 FORMAT(/,1X,'INPUT SIGNAL   T0=',F10.5,2X,'DT=',F10.5,2X,
     1'MAX=',E16.8,2X,'N=',I5)
      END
C
C
      SUBROUTINE FT(LIN,LOU)
      CHARACTER*80 title
      COMPLEX S(2048),REF(2048),SAUX
      COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3,
     1NPL4,IPR,NDER,NINT
      COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN
      COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4)
      COMMON/WIND/FLO,FLEFT,FRO,FRIGHT,FEXP,ILF,IRF,NFW
      COMMON/PLOT1/ARESP(6),ASPECT(6),ASYNT(6),JRESP(4),JSPECT(4),JSYNT(
     14),NOO,RMIN,RMAX,GLR,TRMIN,TRMAX,GLTR,IR,ITR,NRED,NSH,AREDUC,TITLE
     2,NPLOT,MCOMP
C
      IF(NWIN.EQ.0)GO TO 1
      READ(LIN,*)FLO,FLEFT,FRIGHT,FRO,FEXP
      WRITE(LOU,101)FLO,FLEFT,FRIGHT,FRO,FEXP
    1 CONTINUE
C
C     PREPARE FOR FCOOLR
      DO 10 I=1,NPTS
   10 S(I)=CMPLX(F(I),0.0)
      IF(NT.EQ.NPTS)GO TO 11
      NPTS1=NPTS+1
      DO 12 I=NPTS1,NT
   12 S(I)=(0.0,0.0)
   11 CONTINUE
      CALL FCOOLR(NI,S,1.0)
      DO 21 I=1,NT
      A=REAL(S(I))
      B=AIMAG(S(I))
   21 FF(I)=SQRT(A*A+B*B)
C
C     DERIVATIVE (FOR NDER.NE.0) OR INTEGRAL (FOR NINT.NE.0)
C     OF THE INPUT SIGNAL
C
      IF(NDER.NE.0.OR.NINT.NE.0)THEN
      DO 22 I=2,NT
      AUX1=FLOAT(I-1)*DF
      IF(NDER.NE.0)SAUX=CMPLX(0.,-AUX1)
      IF(NINT.NE.0)SAUX=CMPLX(0.,1./AUX1)
      S(I)=S(I)*SAUX
   22 CONTINUE
      END IF
C
C     PLOTTING AMPLITUDE SPECTRUM OF THE INPUT SIGNAL
C
      CALL REDUC(AREDUC,1,IRS,FF,GG,IS)
      IF(IPR.EQ.0)GO TO 25
      WRITE(LOU,120)DF,AREDUC,IRS
      WRITE(LOU,107)(IS(I),I=1,IRS)
   25 CONTINUE
      IF(NPL2.EQ.0)GO TO 28
      NOO=3
      NRED=NPL2
      IF(NPL2.EQ.1)CALL ZPLOT(0.,DF,1,IRS,GG)
      IF(NPL2.EQ.2)CALL ZPLOT(0.,DF,1,IRS,FF)
   28 CONTINUE
C
C      WINDOWING OF AMPLITUDE SPECTRUM
C
      ILF=ILS+1
      IRF=IRS-1
      NFW=IRF-ILF+1
      IF(NWIN.EQ.0.AND.NFILT.EQ.0)GO TO 30
      IF(NWIN.NE.0)CALL FENSTR
      IF(NFILT.NE.0)CALL FILTER
      DO 31 I=1,IRS
      A=REAL(S(I))
      B=AIMAG(S(I))
   31 FF(I)=SQRT(A*A+B*B)
      CALL REDUC(AREDUC,1,IRS,FF,GG,IS)
      IF(IPR.EQ.0)GO TO 32
      WRITE(LOU,121)DF,AREDUC,IRS
      WRITE(LOU,107)(IS(I),I=1,IRS)
   32 CONTINUE
      IF(NPL3.EQ.0)GO TO 30
      NOO=4
      NRED=NPL3
      IF(NPL3.EQ.1)CALL ZPLOT(0.,DF,1,IRS,GG)
      IF(NPL3.EQ.2)CALL ZPLOT(0.,DF,1,IRS,FF)
   30 CONTINUE
C
C     INPUT SIGNAL, FREQUENCY WINDOW AND FILTER APPLIED
C
      DO 40 I=1,NT
   40 REF(I)=S(I)
      DEL=1./FLOAT(NT)
      CALL FCOOLR(NI,REF,-1.0)
      DO 41 I=1,NT
   41 FF(I)=DEL*REAL(REF(I))
      CALL REDUC(AREDUC,1,NPTS,FF,GG,IS)
      IF(IPR.EQ.0)GO TO 42
      WRITE(LOU,122)TO,DT,AREDUC,NPTS
      WRITE(LOU,107)(IS(I),I=1,NPTS)
   42 CONTINUE
      IF(NPL4.EQ.0)GO TO 50
      NRED=NPL4
      NOO=7
      IF(NPL4.EQ.1)CALL ZPLOT(TO,DT,1,NPTS,GG)
      IF(NPL4.EQ.2)CALL ZPLOT(TO,DT,1,NPTS,F)
   50 CONTINUE
      RETURN
C
  101 FORMAT(8F10.5)
  107 FORMAT(20I4)
  120 FORMAT(/1X,'AMPLITUDE SPECTRUM OF THE INPUT SIGNAL'/,
     12X,'DF=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5)
  121 FORMAT(/1X,'AMPLITUDE SPECTRUM OF INPUT SIGNAL,WINDOW APPLIED'/
     1,2X,'DF=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5)
  122 FORMAT(/1X,'INPUT SIGNAL, FREQUENCY WINDOW APPLIED'/2X,'T0=',
     1F10.5,2X,'DT=',F10.5,2X,'MAX=',E16.8,2X,'N=',I5)
      END
C
C
C
C
      SUBROUTINE FENSTR
      COMPLEX S(2048),REF(2048)
      COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN
      COMMON/WIND/FLO,FLEFT,FRO,FRIGHT,FEXP,ILF,IRF,NFW
C
      ILO=FLO/DF+1.1
      IRO=FRO/DF+1.1
      IF(IRO.GT.IRS)IRO=IRS
      IF(ILO.LT.ILS)ILO=ILS
      IF(IRO.GT.NTM)IRO=NTM
      ILEFT=FLEFT/DF+1.1
      IRIGHT=FRIGHT/DF+1.1
      IF(IRIGHT.GT.IRO)IRIGHT=IRO
      IF(ILEFT.LT.ILO)ILEFT=ILO
      IF(IRIGHT.GT.NTM)IRIGHT=NTM
      NLEFT=ILEFT-ILO
      NRIGHT=IRO-IRIGHT
      IF(NLEFT.GT.0)DLEFT=3.14159/FLOAT(NLEFT)
      IF(NRIGHT.GT.0)DRIGHT=3.14159/FLOAT(NRIGHT)
C
      DO 1 I=1,NTM
      J=NT+2-I
      FIF=I-IRIGHT
      FAF=I-ILO
      IF(I.GE.ILEFT.AND.I.LE.IRIGHT)GO TO 1
      IF(I.LE.ILO.OR.I.GE.IRO)S(I)=(0.0,0.0)
      IF(I.GT.ILO.AND.I.LT.ILEFT)
     1S(I)=S(I)*(0.5*(1.-COS(DLEFT*FAF)))**FEXP
      IF(I.GT.IRIGHT.AND.I.LT.IRO)
     1S(I)=S(I)*(0.5*(COS(DRIGHT*FIF)+1.))**FEXP
      IF(I.EQ.1.OR.I.EQ.NTM)GO TO 1
      S(J)=CONJG(S(I))
    1 CONTINUE
      ILF=ILO+1
      IRF=IRO-1
      NFW=IRF-ILF+1
      RETURN
      END
C
C
      SUBROUTINE SFUN
      COMMON/AUXIL/FF(2048),GG(2048),IS(2048),AUX(4),IAUX(4)
      COMMON/SIG/F(2048),TO,DT,NPTS,NO,NSIG,NPL1,NPL2,NPL3,
     1NPL4,IPR,NDER,NINT
C
      IF(IAUX(1).NE.1)GO TO 10
C
C     GABOR SIGNAL
C
      DO 1 I=1,NPTS
      T=TO-AUX(4)+DT*FLOAT(I-1)
      G=6.2831853*T*AUX(1)
      G1=G/AUX(2)
      G2=G1*G1
      F(I)=0.
      IF(G2.LT.20.)F(I)=EXP(-G2)*COS(G+AUX(3))
    1 CONTINUE
      RETURN
   10 CONTINUE
C
      IF(IAUX(1).NE.2)GO TO 20
C
C     BERLAGE SIGNAL
C
      DO 2 I=1,NPTS
      T=TO-AUX(4)+DT*FLOAT(I-1)
      G=6.2831853*T*AUX(1)
      G1=AUX(2)*T
      F(I)=0.
      IF(G1.LT.20..AND.T.GT.0.)F(I)=EXP(-G1)*SIN(G)*T**AUX(3)
    2 CONTINUE
      RETURN
   20 CONTINUE
      IF(IAUX(1).NE.3)GO TO 30
C
C     MULLER SIGNAL
C
      G1=FLOAT(IAUX(2))
      G2=G1/(G1+2)
      G3=G1+2
      G1=3.14159*G1
      G3=3.14159*G3
      DO 3 I=1,NPTS
      T=TO-AUX(2)+DT*FLOAT(I-1)
      F(I)=0.
      IF(T.GT.0.AND.T.LE.AUX(1))
     1F(I)=SIN(G1*T/AUX(1))-G2*SIN(G3*T/AUX(1))
    3 CONTINUE
      RETURN
   30 CONTINUE
      IF(IAUX(1).NE.4)GO TO 40
C
C     RICKER SIGNAL
C
      DO 4 I=1,NPTS
      T=TO-AUX(2)+DT*FLOAT(I-1)
      G1=AUX(1)*T
      G2=2.*G1
      G3=G1*G1
      F(I)=0.
      IF(G3.LT.20.)F(I)=(1.-0.5*G2*G2)*EXP(-G3)
    4 CONTINUE
      RETURN
   40 CONTINUE
      IF(IAUX(1).NE.5)GO TO 50
C
C     BOX-CAR FUNCTION
C
      DO 5 I=1,NPTS
      T=TO+DT*FLOAT(I-1)
      F(I)=0.
      IF(T.GE.AUX(1).AND.T.LE.AUX(2))F(I)=AUX(3)
    5 CONTINUE
      RETURN
   50 CONTINUE
      IF(IAUX(1).NE.6) GO TO 60
C
C     RAMP FUNCTION
C
      DO 6 I=1,NPTS
      F(I)=0.
      T=TO+DT*FLOAT(I-1)
      IF(T.LT.AUX(3).AND.T.GE.AUX(2))F(I)=AUX(4)
      IF(T.LE.AUX(2).AND.T.GE.AUX(1))
     1F(I)=AUX(4)*(T-AUX(1))/(AUX(2)-AUX(1))
    6 CONTINUE
      RETURN
   60 CONTINUE
      IF(IAUX(1).NE.7)GO TO 70
C
C     TRIANGLE FUNCTION
C
      DO 7 I=1,NPTS
      F(I)=0.
      T=TO+DT*FLOAT(I-1)
      IF(T.LE.AUX(2).AND.T.GE.AUX(1))
     1F(I)=AUX(4)*(T-AUX(1))/(AUX(2)-AUX(1))
      IF(T.LE.AUX(3).AND.T.GE.AUX(2))
     1F(I)=AUX(4)*(AUX(3)-T)/(AUX(3)-AUX(2))
    7 CONTINUE
      RETURN
   70 CONTINUE
      IF(IAUX(1).NE.8)GO TO 80
C
C     ONE SAMPLE FUNCTION
C
      DO 8 I=1,NPTS
      F(I)=0.
      IF(I.EQ.IAUX(2))F(I)=AUX(1)
    8 CONTINUE
      RETURN
   80 CONTINUE
      IF(IAUX(1).NE.9)GO TO 90
C
C     SQUARED SIN SIGNAL
C
      DO 9 I=1,NPTS
      T=TO-AUX(2)+DT*FLOAT(I-1)
      G1=3.1415926*T/AUX(1)
      F(I)=0.
      IF(T.GT.0.AND.T.LE.AUX(1))THEN
        G2=SIN(G1)
        F(I)=G2*G2
      END IF
    9 CONTINUE
      RETURN
   90 CONTINUE
C
      RETURN
      END
C
C
C       SUBROUTINE FCOOLR(K,D,SN)
C       FAST FOURIER TRANSFORM OF N = 2**K COMPLEX DATA POINTS
C       REPARTS HELD IN D(1,3,...2N-1), IMPARTS IN D(2,4,...2N).
C
        SUBROUTINE FCOOLR(K,D,SN)
        DIMENSION INU(15),D(4096)
        LX=2**K
        Q1=LX
        IL=LX
        SH=SN*6.28318530718/Q1
        DO 10 I=1,K
        IL=IL/2
10      INU(I)=IL
        NKK=1
        DO 40 LA=1,K
        NCK=NKK
        NKK=NKK+NKK
        LCK=LX/NCK
        L2K=LCK+LCK
        NW=0
        DO 40 ICK=1,NCK
        FNW=NW
        AA=SH*FNW
        W1=COS(AA)
        W2=SIN(AA)
        LS=L2K*(ICK-1)
        DO 20 I=2,LCK,2
        J1=I+LS
        J=J1-1
        JH=J+LCK
        JH1=JH+1
        Q1=D(JH)*W1-D(JH1)*W2
        Q2=D(JH)*W2+D(JH1)*W1
        D(JH)=D(J)-Q1
        D(JH1)=D(J1)-Q2
        D(J)=D(J)+Q1
20      D(J1)=D(J1)+Q2
        DO 29 I=2,K
        ID=INU(I)
        IL=ID+ID
        IF(NW-ID-IL*(NW/IL)) 40,30,30
30      NW=NW-ID
29      CONTINUE
40      NW=NW+ID
        NW=0
        DO 6 J=1,LX
        IF(NW-J) 8,7,7
7       JJ=NW+NW+1
        J1=JJ+1
        JH1=J+J
        JH=JH1-1
        Q1=D(JJ)
        D(JJ)=D(JH)
        D(JH)=Q1
        Q1=D(J1)
        D(J1)=D(JH1)
        D(JH1)=Q1
8       DO 9 I=1,K
        ID=INU(I)
        IL=ID+ID
        IF(NW-ID-IL*(NW/IL)) 6,5,5
5       NW=NW-ID
9       CONTINUE
6       NW=NW+ID
        RETURN
        END
C
C
      SUBROUTINE REDUC(AREDUC,N1,N2,F,G,IS)
      DIMENSION F(1),G(1),IS(1)
      DO 1 I=N1,N2
      IS(I)=0
    1 G(I)=0.
      AREDUC=0.
      DO 2 I=N1,N2
    2 IF(ABS(F(I)).GT.AREDUC)AREDUC=ABS(F(I))
      IF(AREDUC.EQ.0.)GO TO 5
      B=1./AREDUC
      C=999.1*B
      GO TO 6
    5 B=0.
      C=0.
    6 CONTINUE
      DO 3 I=N1,N2
      G(I)=B*F(I)
    3 IS(I)=IFIX(C*F(I))
      RETURN
      END
C
C
      SUBROUTINE FILTER
      COMPLEX S(2048),REF(2048)
      COMMON/SPECTR/S,REF,DF,NI,NT,NT2,NTM,ILS,IRS,NFS,NFILT,NWIN
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'calcops.for'
C     calcops.for
C
C=======================================================================
C