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 grfftsp
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MWV=73400)
      CHARACTER INFILE1*20
      dimension TIME(MWV)
      dimension XR1(MWV),XI1(MWV),YR1(MWV),YI1(MWV),ZR1(MWV),ZI1(MWV)
      COMPLEX*16 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,NDATA
	 READ(11,1130) XR1(I),XI1(I),YR1(I),YI1(I),ZR1(I),ZI1(I)
         IF(J.NE.NOB1)GOTO 150
	 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 velocities to the accelerations;'
      WRITE(*,*)
     * '   Do not transfer (Output Velocities)       (Enter 0):'
      WRITE(*,*)
     * '   Transfer from velocities to accelerations (Enter 1):'
      READ(*,*)NYES
C     
C * BAND PASS *
C
      WRITE(*,*)
      WRITE(*,*)' Band Pass Filter? (Yes=1 or No=0) :'
      READ(*,*)NBD
C
      pi2=datan(1.d0)*8.d0
      OML1=0.d0
      OMH2=99999.9d0
      IF(NBD.EQ.1)THEN
        f1=DOM*(IOM-1)/pi2
        f2=DOM*(NDATA-1)/pi2
        WRITE(*,*)'   Data for TRAPEZOIDAL BAND PASS FILTER:'
        WRITE(*,*)'   INITIAL AND FINAL Frequencies =',f1,f2
        WRITE(*,*)
        WRITE(*,*)
     *  '   Enter fL1,fL2,fH1 and fH2 (fL1<fL2<fH1<fH2: Hz): '
        READ(*,*)fL1,fL2,fH1,fH2
        OML1=fL1*pi2
        OML2=fL2*pi2
        OMH1=fH1*pi2
        OMH2=fH2*pi2
      END IF
C
C * Write data to a complex array *
C
      DO 200 I = 1,NDATA
           OMEGA=dFLOAT(I-1)*DOM
           IF(OMEGA.lt.OML1)GOTO 200
           IF(OMEGA.GE.OMH2)GOTO 200
           AA=1.d0
           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*dCMPLX(XR1(I),XI1(I))
           FY1(I) = AA*dCMPLX(YR1(I),YI1(I))
           FZ1(I) = AA*dCMPLX(ZR1(I),ZI1(I))
           IF(NYES.EQ.1)THEN
             FX1(I)=FX1(I)*OMEGA*(0.d0,1.d0)
             FY1(I)=FY1(I)*OMEGA*(0.d0,1.d0)
             FZ1(I)=FZ1(I)*OMEGA*(0.d0,1.d0)
           END IF
c        write(*,1500)omega,FX1(I),FY1(I),FZ1(I)
  200 CONTINUE
c
c * Output Fourier Spectrum data *
c
      write(*,*)
      write(*,*)'The Fourier Spectrum data are stored in "amp.csv"'
c      write(*,*)'The amplitude of "amp.csv" is *100'
      open(10,file='amp.csv',status='unknown')
        write(10,*)'frequency(Hz), X(NS), Y(EW), Z(UD)'
        DO 250 I = 1,NDATA
          OM=dFLOAT(I-1)*DOM
          f=om/pi2
c          x=cdabs(FX1(I))*100.d0
c          y=cdabs(FY1(I))*100.d0
c          z=cdabs(FZ1(I))*100.d0
          x=cdabs(FX1(I))
          y=cdabs(FY1(I))
          z=cdabs(FZ1(I))
          write(10,1500)f,x,y,z
  250   CONTINUE
      close(10)
 1500 format(f8.4,6(',',e12.4))
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(*,*)
      WRITE(*,*)' Increase NT 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) = dCMPLX(0.0,0.0)
           FY1(J) = dCMPLX(0.0,0.0)
           FZ1(J) = dCMPLX(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 = dACOS(-1.d0)
c      FOM=0.0
      DUR = (2.d0*PI)/DOM
      NFAST = NUM*2
      DT = DUR/dFLOAT(NFAST)
C
C * INVERSE FOURIER TRANSFORM  *
C
      DO 400 I=2,NFAST/2
           FX1(NFAST+2-I)=dCMPLX(dREAL(FX1(I)),-dIMAG(FX1(I)))
           FY1(NFAST+2-I)=dCMPLX(dREAL(FY1(I)),-dIMAG(FY1(I)))
           FZ1(NFAST+2-I)=dCMPLX(dREAL(FZ1(I)),-dIMAG(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) = dFLOAT(I-1)*DT
         a=exp(omi*time(i))
         XR1(I)=dREAL(FX1(I))/DUR*a
         YR1(I)=dREAL(FY1(I))/DUR*a
         ZR1(I)=dREAL(FZ1(I))/DUR*a
  450 CONTINUE
C
      WRITE(*,*)
      WRITE(*,1550) 
 1550 FORMAT(1X,'Change from (m) to (cm)? (yes=1)    : ')
      READ(*,*)ICM
      ACM=1.d0
      IF(ICM.EQ.1)ACM=100.d0
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.d0
      YMAX=0.d0
      ZMAX=0.d0
      XMIN=0.d0
      YMIN=0.d0
      ZMIN=0.d0
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.csv"'
      WRITE(*,*)'You can plot the waveform using "wave.m" for Matlab'
      WRITE(*,*)
      OPEN(11,FILE='wave.csv',STATUS='unknown')
        write(11,*)'time(s), X(NS), Y(EW), Z(UD)'
        do i=1,NFAST
c         write(11,2000)TIME(I),XR1(I)-xr1(1),YR1(I)-yr1(1),ZR1(I)-zr1(1)
         write(11,2000)TIME(I),XR1(I),YR1(I),ZR1(I)
        end do
      CLOSE(11)
 2000 FORMAT(F12.6,3(',',E12.4))
c
      STOP
      END
c
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
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 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=dCMPLX(0.d0,PAI*dFLOAT(IND*(K-1))/dFLOAT(KMAX))
        DO 160 I=K,N,ISTEP
          J=I+KMAX
          TEMP=A(J)*CdEXP(THETA)
          A(J)=A(I)-TEMP
          A(I)=A(I)+TEMP
  160   CONTINUE
  170 CONTINUE
      KMAX=ISTEP
      GO TO 150
C
      END

