c
c *** Regularized Slip Velocity Function using Exponential Function ***
c
      PROGRAM ExpSlipV
c 
      IMPLICIT REAL*8(A-H,O-Z) 
      PARAMETER(MOM=16400)
      complex*16 svAmp(MOM),examp,im
      dimension svWave(MOM)
c
      dt=0.01d0
      nt=8192
c      Tau=0.5d0
      Tau=1.d0
      pai=3.141592654d0
      Tau2=Tau*Tau
c
      write(*,*)
      write(*,*)'Slip Velocity Function using Exponential Function'
      write(*,*)
      write(*,*)' dt,nt=',dt,nt
      write(*,*)' Tau(s) =',Tau
c
      write(*,*)
      write(*,*)' Slip functions are stored in Exp_Slip.csv'
      write(*,*)' Fourier spectrum is stored in Exp_Amp.csv'
      open (11,file='Exp_Slip.csv',status='unknown')
      write(11,*)'time(s),velocity,displacement,acceleration'
      svWave(1)=0.d0
      svS=0.d0
      svA=1/Tau/Tau
      write(11,1000)0.d0,svWave(1),svS,svA
c
      svS=0.d0
      do 100 i=2,nt
        t=dfloat(i-1)*dt
        svWave(i)=t/(Tau2*dexp(t/Tau))
        svS=svS+(svWave(i)+svWave(i-1))*dt/2.d0
        svA=(svWave(i)-svWave(i-1))/dt
        write(11,1000)t,svWave(i),svS,svA
  100 continue
      close(11)
 1000 format(f8.4,3(',',f15.7))
c
c   * Fourier Transform of the Slip-Rate Dunction *
c
      dur=dt*dfloat(nt)
      dom=2.d0*pai/dur
c      iom=1 
c      nom=int(duration/period)
      nom=nt/2
      EOM=DOM*NOM
      WRITE(*,*)
      WRITE(*,*)' DUR,DOM,NOM,EOM=',DUR,DOM,NOM,EOM
      WRITE(*,*)' Highest Frequency (hz)=',EOM/6.2831853070d0
c
      nfast=nt
      if(nfast.ge.MOM)then
        write(*,*)'from Triangle: nfast > MOM !',nfast,MOM
c        write(10,*)'from Triangle: nfast > MOM !',nfast,MOM
      end if
      do i=1,nt
        svAmp(i)=dCMPLX(svWave(i),0.0d0)
      end do
      do i=nt+1,nfast
        svAmp(I)=dCMPLX(0.d0,0.0d0)
      end do
c
      CALL dWVFAST(nfast,svAmp,MOM,-1)

      open (12,file='Ex_Amp.csv',status='unknown')
c      write(12,*)'freq(Hz),abs(velocity),real,imag'
      write(12,*)'freq(Hz),abs(vel),abs(acc),theory(vel),theory(acc)'
      im=dCMPLX(0.d0,1.0d0)
      do io=2,nom
        om=DOM*DFLOAT(io-1)
        freq=om/(2.d0*pai)
        svAmp(io)=dconjg(svAmp(io))*dt
        examp=1.d0/((1-im*om*Tau)*(1-im*om*Tau))
        GA=cdabs(examp)
c        write(12,1100)freq,cdabs(svAmp(io)),svAmp(io)
        write(12,1100)freq,cdabs(svAmp(io)),cdabs(svAmp(io))*om,GA,GA*om
      end do
      close(12)
 1100 format(f8.4,4(',',e12.4))
c
      write(*,*)
      write(*,*)'Hit any key to stop'
      read(*,*)
c  
      stop
      end
c
c
C *** FOURIER TRANSFORM OR INVERSE FOURIER TRANSFORM
C
      SUBROUTINE dWVFAST(N,A,MOM,IND)
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                
                
                    

