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
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MWV=513)
      dimension XR(MWV),XR1(MWV)
      COMPLEX*16 FX(MWV),ft(MWV),ft1(MWV),ft2(MWV),ft3(MWV),im,dif
C
      dt=0.1
      nt=32
      dur=dt*nt
c      nfast=nt*2
c      dur2=dur*2
      write(*,*)'dt,nt',dt,nt
      write(*,*)'nfast,dur,dur2',nast,dur,dur2
c
c      rise=1
      vel=1.0
      write(*,*)'Enter Rise Time (sec)'
      read(*,*)rise
      write(*,*)'rise,disp',rise,disp
c
      xr(1)=0.0
      do i=2,nt
        t=dt*(i-1)
        if(t.le.rise)then
          xr(i)=vel
        else
          xr(i)=0.0
        end if
      end do
      write(*,*)'original xr'
      do i=1,nt
        fx(i)=dcmplx(xr(i),0.d0)
        write(*,*)i,xr(i)
      end do
c
      write(*,*)'Differenciate? (yes=1)'
      read(*,*)idif
C
      write(*,*)'fx'
      CALL WVFAST(Nt,FX,MWV,-1)
      df=1.0/dur
      im=(0.0,1.0)
      dom=df*dacos(-1.d0)*2.0
      dif=(1.d0,0.d0)
      do i=1,nt/2+1
        om=dom*(i-1)
        if(idif.eq.1)dif=im*om
        fx(i)=fx(i)*dt*dif
c        write(*,*)i,om,dif
c        write(*,*)i,cdabs(fx(i))
      end do
      do i=2,nt/2
        fx(nt+2-i)=dcmplx(dreal(fx(i)),-dimag(fx(i)))
c        write(*,*)i,cdabs(fx(i)),fx(i)
c        write(*,*)i,cdabs(fx(i))
      end do
      do i=1,nt
        write(*,*)i,cdabs(fx(i)),fx(i)
c        write(*,*)i,cdabs(fx(i))
      end do
c
      open(9,file='fft1_d',status='unknown')
      write(9,*)rise,disp,nt
      do i=1,nt+1
        f=df*(i-1)
        om=dom*(i-1)
        if(i.eq.1)then
          ft(i)=disp*(2*dur-rise-dt)/2.0
        else
          ot1=om*rise/2.0
          ot2=om*(dur-rise-dt)/2.0
          ot3=om*dt/2.0
         ft1(i)=im*(dsin(ot1)/ot1-cdexp(ot1*im))*rise*Disp/ot1/2.0
         ft2(i)=disp*(dur-rise-dt)*dsin(ot2)*cdexp(im*(ot2+om*rise))/ot2
         ft3(i)=im*(1.d0-dsin(ot3)/ot3*cdexp(ot3*im))
     *            *cdexp(im*om*(dur-dt))*dt*Disp/ot3/2.0
         ft(i)=ft1(i)+ft2(i)+ft3(i)
        end if
c        write(*,*)ft1(i)
c        write(*,*)ft2(i)
c        write(*,*)ft3(i)
c        write(*,9999)f,ft(i),fx(i)
c        write(*,9999)f,ft1(i),ft2(i),ft1(i)+ft2(i),fx(i)
        write(9,2000)f,fx(i)
c      stop
      end do

 2000 format(3e12.4)
c 9999 format('f,ft1,ft2,ft,fx',9e12.4)
 9999 format('f,ft,fx',9e12.4)
C
      write(*,*)'xr'
      CALL WVFAST(Nt,FX,MWV,1)
      DO I = 1,nt
         XR1(I)=dREAL(FX(I))/DUR
         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
      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

