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=30)
      REAL XR(MWV)
      COMPLEX*8 FX(MWV)
C
      dt=0.2
      nt=4
      dur=dt*nt
      nfast=nt*2
      dur2=dur*2
      write(*,*)'dt,nt',dt,nt
      write(*,*)'nfast,dur,dur2',nast,dur,dur2
c
      xr(1)=0.0
      xr(2)=1.2
      xr(3)=2.3
      xr(4)=0.0
c      xr(5)=3.0
c      xr(6)=2.0
c      xr(7)=1.0
c      xr(8)=0.0
c
      write(*,*)'original xr'
      do i=1,nt
        fx(i)=cmplx(xr(i),0.0)
        write(*,*)i,xr(i)
      end do
C
      write(*,*)'fx'
      CALL WVFAST(NFAST,FX,MWV,-1)
      do i=1,nfast
        fx(i)=fx(i)*dt
        write(*,*)i,cabs(fx(i)),fx(i)
      end do
C
      write(*,*)'xr'
      CALL WVFAST(NFAST,FX,MWV,1)
      DO I = 1,NFAST
         XR(I)=REAL(FX(I))/DUR2
         write(*,*)i,xr(i)
      end do
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

