c
c *** Regularized Slip Velocity Function using multiple triangles (Hisada, 2001) ***
c
      PROGRAM triangle
c 
      IMPLICIT REAL*8(A-H,O-Z) 
      PARAMETER(MOM=16400)
      complex*16 triAmp(MOM)
      dimension triSV(MOM)
c
      Dmax=1.d0
c   * Data 1: Shallow Slip model for Yoffe (TauS=1.2 s, TauR=2.8 s) *
c      dt=0.01d0
c      nt=8192
c      fmax=0.55d0
c      Tr=1.162d0
c      Ar=0.95d0
c      Nv=8
c   * Data 2: Deep Slip model for Yoffe (TauS=0.02, TauR=3.1) *
      dt=0.005d0
      nt=4096
      fmax=10.0d0
      Tr=1.633d0
      Ar=0.95d0
      Nv=8
c
      TauMin=1.d0/fmax
      TauMax=Tr**(Nv-1)/fmax
c
      write(*,*)
      write(*,*)'Regularized Triangle Slip Velocity Function'
      write(*,*)
      write(*,*)' dt,nt=',dt,nt
      write(*,*)' Dmax,fmax,Tr,Ar=',Dmax,fmax,Tr,Ar
      write(*,*)' TauMin,TauMax=',TauMin,TauMax
      pai=3.141592654d0
c
      write(*,*)
      write(*,*)' Slip functions are stored in Tri_Slip.csv'
      write(*,*)' Fourier spectrum of slip velocity is stored in Tri_Amp
     *.csv'
      open (11,file='Tri_Slip.csv',status='unknown')
      write(11,*)'time(s),velocity,displacement,acceleration'
      triSV(1)=0.d0
      write(11,1000)0.d0,triSV(1),0.d0,0.d0
c
      triS=0.d0
      do 100 i=2,nt
c
        t=dfloat(i-1)*dt
        triSV(i)=0.d0
c
        ArSum=0.d0
        do 110 j=1,Nv
          Tauj=TauMin*(Tr**(j-1))
          Arj=Ar**(j-1)
          ArSum=ArSum+Arj
          Tau1j=Tauj/Tr
          Tau2j=Tauj-Tau1j
          Vj=2.d0/Tauj
          A1j=Vj/Tau1j
          A2j=Vj/Tau2j
          if(t.le.Tau1j)then
            sVj=A1j*t
          elseif(t.le.Tauj)then
            sVj=A2j*(Tauj-t)
          else
            sVj=0.d0
          end if
          triSV(i)=triSV(i)+sVj*Arj
  110   continue
        triSV(i)=triSV(i)*Dmax/ArSum
c
        if(i.eq.2)then
          triS=triS+triSV(i)*dt/2.d0
          triA=triSV(i)/dt
        else
          triS=triS+(triSV(i)+triSV(i-1))*dt/2.d0
          triA=(triSV(i)-triSV(i-1))/dt
        end if    
        write(11,1000)t,triSV(i),triS,triA
c
  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
        triAmp(i)=dCMPLX(triSV(i),0.0d0)
      end do
      do i=nt+1,nfast
        triAmp(I)=dCMPLX(0.d0,0.0d0)
      end do
c
      CALL dWVFAST(nfast,triAmp,MOM,-1)

      open (12,file='Tri_Amp.csv',status='unknown')
c      write(12,*)'freq(Hz),abs(velocity),real,imag'
      write(12,*)'freq(Hz),abs(vel),abs(acc)'
      write(12,1100)0.d0,1.d0,0.d0
      do io=2,nom
        om=DOM*DFLOAT(io-1)
        freq=om/(2.d0*pai)
        triAmp(io)=dconjg(triAmp(io))*dt
c        write(12,1100)freq,cdabs(triAmp(io)),triAmp(io)
        write(12,1100)freq,cdabs(triAmp(io)),cdabs(triAmp(io))*om
      end do
      close(12)
 1100 format(f8.4,3(',',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                
                
                    

