C
C *** grfft.f is a FORTRAN code to obtain seismograms in the time domain,
C     which are computed from grpoint.f and grfault.f
C
      PROGRAM grfft
C
      PARAMETER(MWV=4100)
      CHARACTER INFILE1*20
      REAL TIME(MWV)
      REAL XR1(MWV),XI1(MWV),YR1(MWV),YI1(MWV),ZR1(MWV),ZI1(MWV)
      COMPLEX*8 FX1(MWV),FY1(MWV),FZ1(MWV)
C
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) 
     *'"grfft" - To obtain seismograms for grflt3.f and grpoint.f'
      WRITE(*,*)'           at a single observation point.'
      WRITE(*,*)
C
C * Read data *
c
      WRITE(*,1010) 
 1010 FORMAT('Choose the input datafile for FFT (Enter Number 1 to 6):'/ 
     *' 1) grfault.dxyz (x,y,z-comp. of disp. from seismic sources)'/
     *' 2) grfault.dqvw (u,v,w-comp. of disp. from seismic sources)'/
     *' 3) grpoint.dxyz (x,y,z-comp. of disp. from point sources)'/
     *' 4) grpoint.sxyz (x,y,z-comp. of stress from point sources)'/
     *' 5) grpoint.dqvw (u,v,w-comp. of disp.  from point sources)'/
     *' 6) grpoint.sqvw (u,v,w-comp. of stress from point sources)')
      READ(*,*) IDF
C    
      IF(IDF.EQ.1)THEN
        INFILE1='grfault.dxyz'
      ELSE 
        IF(IDF.EQ.2)THEN
          INFILE1='grfault.dqvw'
        ELSE 
          IF(IDF.EQ.3)THEN
            INFILE1='grpoint.dxyz'
          ELSE
            IF(IDF.EQ.4)THEN
              INFILE1='grpoint.sxyz'
            ELSE 
              IF(IDF.EQ.5)THEN
                INFILE1='grpoint.dqvw'
              ELSE 
                IF(IDF.EQ.6)THEN
                  INFILE1='grpoint.sqvw'
                END IF
              END IF
            END IF
          END IF
        END IF
      END IF
c
      OPEN(11,FILE=INFILE1,STATUS='OLD')
c
      READ(11,1100)DOM,IOM,nom,omi
 1100 FORMAT(1X,F10.0,2I10,f10.0)
      WRITE(*,*)
      WRITE(*,*)' DOM,IOM,nom,omi =',DOM,IOM,nom,omi
      NDATA=nom+1
      WRITE(*,*)' New NDATA including OM=0 =',NDATA
      READ(11,1110)NOB
 1110 FORMAT(1X,I10///)
      WRITE(*,*)
      WRITE(*,*)' Total number of receivers =',NOB
      WRITE(*,1120) 
 1120 FORMAT(1X,' Enter the receiver number to output : ')
      READ(*,*)NOB1
      WRITE(*,*)
      WRITE(*,*)' XR,XI,YR,YI,ZR,ZI='
      DO 150 J=1,NOB1
      DO 150 I = IOM+1,NDATA
	 READ(11,1130) XR1(I),XI1(I),YR1(I),YI1(I),ZR1(I),ZI1(I)
         IF(J.NE.NOB1)GOTO 150
c	 WRITE(*,1130) XR1(I),XI1(I),YR1(I),YI1(I),ZR1(I),ZI1(I)
  150 CONTINUE
 1130 FORMAT(1X,6E12.5)
C
      CLOSE(11)
C     
C   * change of signs *
C
      WRITE(*,*)
      WRITE(*,*)' Change signs of amplitudes? (Yes=1 or No=0) :'
      READ(*,*)MYES
C 
      IF(MYES.EQ.1)THEN
C
        WRITE(*,*)'   Change signs of X? (Yes=1 or No=0) :'
        READ(*,*)MY1
C
         IF(MY1.EQ.1)THEN
           DO 180 I=1,NDATA
             XR1(I)=-XR1(I)
             XI1(I)=-XI1(I)
  180      CONTINUE
         END IF
C
        WRITE(*,*)'   Change signs of Y? (Yes=1 or No=0) :'
        READ(*,*)MY2
C
         IF(MY2.EQ.1)THEN
           DO 182 I=1,NDATA
             YR1(I)=-YR1(I)
             YI1(I)=-YI1(I)
  182      CONTINUE
         END IF
C
        WRITE(*,*)'   Change signs of Z? (Yes=1 or No=0) :'
        READ(*,*)MY3
C
         IF(MY3.EQ.1)THEN
           DO 184 I=1,NDATA
             ZR1(I)=-ZR1(I)
             ZI1(I)=-ZI1(I)
  184      CONTINUE
         END IF
C
      END IF
C     
C * DIFFERENTIATION *
C
      WRITE(*,*)
      WRITE(*,*)
     * ' You can transfer the displacements to vel. or acc.;'
      WRITE(*,*)
     * '   Do not transfer                              (Enter 0),'
      WRITE(*,*)
     * '   Transfer from displacements to velocities    (Enter 1),'
      WRITE(*,*)
     * '   Transfer from displacements to accelerations (Enter 2):'
      READ(*,*)NYES
C     
C * BAND PASS *
C
      WRITE(*,*)
      WRITE(*,*)' Band Pass Filter? (Yes=1 or No=0) :'
      READ(*,*)NBD
C
      OML1=0.0
      OMH2=99999.9
      IF(NBD.EQ.1)THEN
        WRITE(*,*)'   Data for TRAPEZOIDAL BAND PASS FILTER:'
        WRITE(*,*)'   INITIAL AND FINAL OMs =',DOM*IOM,DOM*nom
        WRITE(*,*)
        WRITE(*,*)
     *  '   Enter OML1,OML2,OMH1 and OMH2 (OML1<OML2<OMH1<OMH2): '
        READ(*,*)OML1,OML2,OMH1,OMH2
      END IF
C
C * Write data to a complex array *
C
      FX1(1) = CMPLX(0.0,0.0)
      FY1(1) = CMPLX(0.0,0.0)
      FZ1(1) = CMPLX(0.0,0.0)
      DO 200 I = 2,NDATA
           OMEGA=FLOAT(I-1)*DOM
           IF(OMEGA.LE.OML1)GOTO 200
           IF(OMEGA.GE.OMH2)GOTO 200
           AA=1.0
           IF(NBD.EQ.1)THEN
             IF(OMEGA.LT.OML2)AA=(OMEGA-OML1)/(OML2-OML1)
             IF(OMEGA.GT.OMH1)AA=(OMH2-OMEGA)/(OMH2-OMH1)
           END IF
	   FX1(I) = AA*CMPLX(XR1(I),XI1(I))
	   FY1(I) = AA*CMPLX(YR1(I),YI1(I))
	   FZ1(I) = AA*CMPLX(ZR1(I),ZI1(I))
           IF(NYES.EQ.1)THEN
             FX1(I)=FX1(I)*OMEGA*(0.0,1.0)
             FY1(I)=FY1(I)*OMEGA*(0.0,1.0)
             FZ1(I)=FZ1(I)*OMEGA*(0.0,1.0)
           END IF
           IF(NYES.EQ.2)THEN
             FX1(I)=-FX1(I)*OMEGA*OMEGA
             FY1(I)=-FY1(I)*OMEGA*OMEGA
             FZ1(I)=-FZ1(I)*OMEGA*OMEGA
           END IF
  200 CONTINUE
c
c * Output Fourier Spectrum data *
c
      write(*,*)
      write(*,*)'The Fourier Spectrum data are stored in "amp_dat"'
      open(10,file='amp_dat',status='unknown')
      DO 250 I = 2,NDATA
        OMEGA=FLOAT(I-1)*DOM
        write(10,1500)omega,FX1(I),FY1(I),FZ1(I)
  250 CONTINUE
 1500 format(f8.4,6e12.4)
      close(10)
C
C * Find closest (but larger) 2**J value to the input data dimension *
C
      WRITE(*,*)
      DO 300 J = 1,999
           N = 2**J
           IF(nom.LE.N)THEN
                NUM = N
                GO TO 20
           END IF
  300 CONTINUE
C
   20 IF(NUM.GT.MWV)THEN
        WRITE(*,*)' NUM > MWV !!!'
        STOP
      END IF

      WRITE(*,*)' N and J of N=2**J =',NUM,J
      WRITE(*,*)' Increase J for smoother seismograms? (Yes=1 or No=0)'
      READ(*,*)IYES
      IF(IYES.EQ.1)THEN
        J=J+1
        NUM=NUM*2
        GO TO 20
      END IF
C
C * Pad the remaining spots (up to N=2**J) of the input array with zeroes *
C
      IF(NUM*2.GT.MWV)THEN
      WRITE(*,*)
      WRITE(*,*)' NUM*2.GT.MWV !!!!!'
      WRITE(*,*)' Increase MWV in the parameter statement!'
      STOP
      ENDIF
C
      DO 350 J = NDATA+1,NUM+1
           FX1(J) = CMPLX(0.0,0.0)
           FY1(J) = CMPLX(0.0,0.0)
           FZ1(J) = CMPLX(0.0,0.0)
  350 CONTINUE
C
C * Define some variables (First freq., duration of seismogram, number
C     of IFT data points, total number of output points, delta time
C     of seismogram)
C
      PI = ACOS(-1.0)
c      FOM=0.0
      DUR = (2.0*PI)/DOM
      NFAST = NUM*2
      DT = DUR/FLOAT(NFAST)
C
C * INVERSE FOURIER TRANSFORM  *
C
      DO 400 I=2,NFAST/2
           FX1(NFAST+2-I)=CMPLX(REAL(FX1(I)),-AIMAG(FX1(I)))
           FY1(NFAST+2-I)=CMPLX(REAL(FY1(I)),-AIMAG(FY1(I)))
           FZ1(NFAST+2-I)=CMPLX(REAL(FZ1(I)),-AIMAG(FZ1(I)))
  400 CONTINUE
C
      CALL WVFAST(NFAST,FX1,MWV,1)
      CALL WVFAST(NFAST,FY1,MWV,1)
      CALL WVFAST(NFAST,FZ1,MWV,1)
C
C * WAVE DATA *
C
      DO 450 I = 1,NFAST
         TIME(I) = FLOAT(I-1)*DT
         a=exp(omi*time(i))
         XR1(I)=REAL(FX1(I))/DUR*a
         YR1(I)=REAL(FY1(I))/DUR*a
         ZR1(I)=REAL(FZ1(I))/DUR*a
  450 CONTINUE
C
      WRITE(*,*)
      WRITE(*,1550) 
 1550 FORMAT(1X,'Change from (m) to (cm)? (yes=1)    : ')
      READ(*,*)ICM
      ACM=1.0
      IF(ICM.EQ.1)ACM=100.0
C
      DO 462 I=1,NFAST+1
        XR1(I)=ACM*XR1(I)
        YR1(I)=ACM*YR1(I)
        ZR1(I)=ACM*ZR1(I)
  462 CONTINUE
C
C * MAX AND MIN AMPLITUDE *
C
      XMAX=0.0 
      YMAX=0.0 
      ZMAX=0.0 
      XMIN=0.0 
      YMIN=0.0 
      ZMIN=0.0 
C
      DO 550 I=1,NFAST
        IF(XR1(I).GT.XMAX)XMAX=XR1(I)
        IF(XR1(I).LT.XMIN)XMIN=XR1(I)
        IF(YR1(I).GT.YMAX)YMAX=YR1(I)
        IF(YR1(I).LT.YMIN)YMIN=YR1(I)
        IF(ZR1(I).GT.ZMAX)ZMAX=ZR1(I)
        IF(ZR1(I).LT.ZMIN)ZMIN=ZR1(I)
  550 CONTINUE
C
      AMAX=XMAX
      AMIN=XMIN
      IF(YMAX.GT.AMAX)AMAX=YMAX
      IF(ZMAX.GT.AMAX)AMAX=ZMAX
      IF(YMIN.LT.AMIN)AMIN=YMIN
      IF(ZMIN.LT.AMIN)AMIN=ZMIN
C
      AAMAX=AMAX
      IF(ABS(AMIN).GT.AAMAX)AAMAX=ABS(AMIN)
C
      WRITE(*,*)
      WRITE(*,*)' Nmber of data points    =',NFAST
      WRITE(*,*)' Duration and Delta Time =',DUR,DT
      WRITE(*,*)' Initial and final times =',0.0,DUR-DT
      WRITE(*,*)
      WRITE(*,*)' Min X, Y and Z          =',XMIN,YMIN,ZMIN
      WRITE(*,*)' Max X, Y abd Z          =',XMAX,YMAX,ZMAX
      WRITE(*,*)' Min and Max amplitudes  =',AMIN,AMAX
      WRITE(*,*)' Absolute Max amplitude  =',AAMAX
C      
C * Output data file  *
C
      WRITE(*,*)
      WRITE(*,*)'Output data are stored in "wave_dat" and "amp_dat"'
      WRITE(*,*)
     * 'You can plot the waveform using "wave.m" and "amp.m" for Matlab'
      WRITE(*,*)
      write(*,*)'Note: When omi > 0, "amp_dat" are damped amplitudes!'
      WRITE(*,*)
      OPEN(11,FILE='wave_dat',STATUS='unknown')
c      WRITE(11,*)NFAST,DT,AMIN,AMAX
      do i=1,NFAST
      write(11,2000)TIME(I),XR1(I),YR1(I),ZR1(I)
      end do
      CLOSE(11)
 2000 FORMAT(F10.5,3E12.4)
c
      STOP
      END
C
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE WVFAST(N,A,MOM,IND)
C
C * N   : NUMBER OF FOURIER COEFF. OR WAVE POINTS (=2**I)
C * A(I): FOURIER COEFFICIENTS ( N/2+1 ) OR WAVE DATA (REAL PART)
C    ( FOURIER SPECTRUM = FOURIER COEFFICIENTS * DURATION/2.0 )
C * MOM : N IN MAIN DIMMENSION
C * IND : = -1 : TRANS.,   = 1 : INVERSE TRANS.
C
      COMPLEX*8 A(MOM),TEMP,THETA
C
C * FFT and inverse FFT *
C
      PAI=3.1415926540d0
      J=1
      DO 140 I=1,N
        IF(I.GE.J) GO TO 110
        TEMP=A(J)
        A(J)=A(I)
        A(I)=TEMP
  110   M=N/2
  120   IF(J.LE.M) GO TO 130
        J=J-M
        M=M/2
        IF(M.GE.2) GO TO 120
  130   J=J+M
  140 CONTINUE
      KMAX=1
  150 IF(KMAX.GE.N) RETURN
      ISTEP=KMAX*2
      DO 170 K=1,KMAX
        THETA=CMPLX(0.0,PAI*FLOAT(IND*(K-1))/FLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*CEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END

