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=2050,mobv=12,mobh=122)
      CHARACTER INFILE1*20
      REAL XR1(MWV),XI1(MWV),YR1(MWV),YI1(MWV),ZR1(MWV),ZI1(MWV)
      COMPLEX FX1(MWV),FY1(MWV),FZ1(MWV)
      complex gx(mobh,mobv,mwv),gy(mobh,mobv,mwv),gz(mobh,mobv,mwv)
C
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*) 
     *'"grfft8" - To obtain seismograms for NS,EW,UD components'
      WRITE(*,*)
C
C * Read data from 'grflt8.in' *
C
      PI = ACOS(-1.0)
      write(*,*)'Data from "grflt8.in"'

      OPEN (7,FILE='grflt8.in',STATUS='OLD')
      read(7,1100)dt,nt
      a=nt
 1    a=a/2
      if(a.gt.1.d0)then
        goto 1
      else
        if(a.ne.1.d0)then
          write(*,*)
          write(*,*)'NT is not a power of 2 !!'
          write(10,*)'NT is not a power of 2 !!'
          stop
        end if
      end if
      read(7,1110)period
 1100 format(/f8.0,i8)
 1110 format(f8.0)
      duration=dt*nt
      dom=pi*2.0/duration
      iom=1
      nom=duration/period
      if(nom.gt.mwv)then
        write(*,*)'nom > mwv !'
        stop
      end if
      write(*,*)
      write(*,*)'Delta Time (s) and Number of Times =',dt,nt
      write(*,*)'Minimum Period to be computed (s)  =',period
      write(*,*)'Number of Omega and Delta Omega    =',nom,dom
C
      READ(7,1200)NL
c      write(*,*)'nl=',nl
      DO 100 IL=1,NL
        READ(7,*)
  100 CONTINUE
 1200 FORMAT(/I10)
C
      READ(7,1300)NF
c      write(*,*)'nf=',nf
 1300 FORMAT(/2I10)
      DO 250 I=1,NF
        READ(7,*)
        READ(7,*)
        READ(7,*)
        READ(7,*)
        READ(7,*)
  250 CONTINUE
C
      READ(7,*)
      READ(7,*)
      READ(7,*)
      READ(7,*)
      READ(7,*)
      READ(7,*)
      READ(7,*)
C
      READ(7,1700)ISIGN
      IF(ISIGN.EQ.1)WRITE(*,2700)
      IF(ISIGN.EQ.0)WRITE(*,2710)
 1700 FORMAT(/I10)
 2700 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE CHANGED')
 2710 FORMAT(/' SIGNS OF IMAGINARY PARTS FOR RESULTS ARE NOT CHANGED')
c
      read(7,1500)
 1500 format(//)
      read(7,1510)nobh
      read(7,1510)nobv
 1510 format(i10)
      if(nobh.gt.mobh)then
        write(*,*)'nobh > mobh !'
        stop
      end if
      if(nobv.gt.mobv)then
        write(*,*)'nobv > mobv !'
        stop
      end if
      write(*,*)
      write(*,*)'Number of Receivers along V and H axes =',NOBV,NOBH

      close(7)
c
c *** Read Green'f Functions ***
c
      open(12,file='grflt.dxyz',form='unformatted',status='old')

      do 200 i=1,nom
      do 200 j=1,nobv
      do 200 k=1,nobh
      read(12)gx(k,j,i),gy(k,j,i),gz(k,j,i)
c      write(*,*)gx(k,j,i),gy(k,j,i),gz(k,j,i)
 200  continue

      close(12)

      if(nobv.gt.mobv)then
        write(*,*)'nobv > mobv !'
        stop
      end if
      if(nobh.gt.mobh)then
        write(*,*)'nobh > mobh !'
        stop
      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 *** Change Unit ***
c
      WRITE(*,*)
      WRITE(*,*)'Unit is changed fro M to CM' 
      ACM=100.0
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
      DUR = (2.0*PI)/DOM
      ndata=nom+1
      num=nt/2
      NFAST = NUM*2
c
      WRITE(*,*)
      WRITE(*,*)'Output data are stored in "wave_dat"'
      WRITE(*,*)'  1st line: iV and iH on the Receiver Plane'
      WRITE(*,*)'  2nd line: NS, EW, and UD components at t=0,'
      write(*,*)'  3rd line: NS, EW, and UD components at t=dt,'
      write(*,*)'           ..........................'
      write(*,*)' NTth line: NS, EW, and UD comp. at t=duration-dt'
      write(*,*)'  (Those are Repeated until iV=nV and iH=nH)'
      write(*,*)'  (Data at t=duration are same as those at t=0)'
c
c *** Big Loop for Observation Points ***
c
      OPEN(11,FILE='wave_dat',STATUS='unknown')
      do 400 iv=1,nobv
      do 400 ih=1,nobh

c   * Values at om=0 are assumed to be zero *
      xr1(1)=0.0
      xi1(1)=0.0
      yr1(1)=0.0
      yi1(1)=0.0
      zr1(1)=0.0
      zi1(1)=0.0

      do 410 i=1,nom
        xr1(i+1)=real(gx(ih,iv,i))*acm
        xi1(i+1)=imag(gx(ih,iv,i))*acm
        yr1(i+1)=real(gy(ih,iv,i))*acm
        yi1(i+1)=imag(gy(ih,iv,i))*acm
        zr1(i+1)=real(gz(ih,iv,i))*acm
        zi1(i+1)=imag(gz(ih,iv,i))*acm
        if(isign.eq.1)then
          xi1(i+1)=-xi1(i+1)
          yi1(i+1)=-yi1(i+1)
          zi1(i+1)=-zi1(i+1)
        end if
 410  continue
C     
C   * change of signs for Z components *
C
      DO 420 I=1,NDATA
        ZR1(I)=-ZR1(I)
        ZI1(I)=-ZI1(I)
 420  CONTINUE
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 440 I = 2,NDATA
           OMEGA=FLOAT(I-1)*DOM
c           IF(OMEGA.LE.OML1)GOTO 200
c           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
  440 CONTINUE
C
C * Pad the remaining spots (up to N=2**J) of the input array with zeroes *
C
      DO 450 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)
  450 CONTINUE
C
C * INVERSE FOURIER TRANSFORM  *
C
      DO 460 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)))
  460 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 470 I = 1,NFAST
         XR1(I)=REAL(FX1(I))/DUR
         YR1(I)=REAL(FY1(I))/DUR
         ZR1(I)=REAL(FZ1(I))/DUR
  470 CONTINUE
C      
C * Output data file  *
C
      write(11,4010)iv,ih
 4010 format(2i5)
      do i=1,NFAST
      write(11,4000)XR1(I),YR1(I),ZR1(I)
      end do
 4000 FORMAT(3E12.4)
c
 400  continue
      CLOSE(11)
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

















