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=513)
      REAL XR(MWV),XR1(MWV)
      COMPLEX*8 FX(MWV),ft(MWV),im
C
      dt=0.1
      nt=32
      dur=dt*nt
      nfast=nt*2
      dur2=dur*2
      write(*,*)'dt,nt',dt,nt
      write(*,*)'nfast,dur,dur2',nast,dur,dur2
c
      rise=1
      disp=1.0
      write(*,*)'Enter Rise Time (sec)'
      read(*,*)rise
      write(*,*)'rise,disp',rise,disp
c
      do i=1,nt
        t=dt*(i-1)
        if(t.le.rise)then
          xr(i)=disp*t/rise
        else
          xr(i)=disp
        end if
      end do
      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

      open(9,file='fft1_d',status='unknown')
      write(9,*)rise,disp,nt
      df=1.0/dur2
      dom=df*acos(-1.0)*2
      im=(0.0,1.0)
      do i=1,nt+1
        f=df*(i-1)
        om=dom*(i-1)
      om=0.1
      i=2
        if(i.eq.1)then
          ft(i)=disp*(2*dur-rise)/2.0
        else
          ot=om*rise/2.0
          ft(i)=(im*sin(ot)*cexp(ot*im)/ot
     *          +cexp(im*om*dur))*disp/(im*om)
        end if
        write(*,9999)f,ft(i),fx(i)
        write(9,2000)f,fx(i)
      stop

      end do

 2000 format(3e12.4)
 9999 format('f,ft,fx',5e12.4)
C
      write(*,*)'xr'
      CALL WVFAST(NFAST,FX,MWV,1)
      DO I = 1,nt
         XR1(I)=REAL(FX(I))/DUR2
         write(*,*)i,xr1(i)
         write(9,2000)dt*(i-1),xr(i),xr1(i)
      end do
      close(9)
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

