c
c *** Regularized Yoffe Slip Velocity Function ***
c
      PROGRAM yoffe
c 
      IMPLICIT REAL*8(A-H,O-Z) 
      PARAMETER(MOM=16400)
      complex*16 yoffeAmp(MOM)
      dimension yoffeSV(MOM)
c
      Dmax=1.d0
c   * The shallow slip model of the Hutagawa Fault *
c      dt=0.01d0
c      nt=8192
c      TauS=1.2d0
c      TauR=2.8d0
c   * The slip model for Nakamura&Miyatake model *
      dt=0.005d0
      nt=4096
      TauS=0.02d0
      TauR=3.1d0
c
      write(*,*)
      write(*,*)'Regularized Yoffe Slip Velocity Function'
      write(*,*)
      write(*,*)' dt,nt=',dt,nt
      write(*,*)' Dmax,TauS,TauR=',Dmax,TauS,TauR
      if(2.d0*TauS.gt.TauR)then
         write(*,*)'2.d0*TauS < TauR, the STOP!'
         write(*,*)'TauS,2.d0*TauS,TauR=',TauS,2.d0*TauS,TauR
         stop
      end if    
      TauS2=TauS*TauS
      TauR2=TauR*TauR
      pai=3.141592654d0
      CK=2.d0/(pai*tauR*TauS2)
c
      write(*,*)
      write(*,*)' Slip functions are stored in Yoffe_Slip.csv'
      write(*,*)' Fourier spectrum of slip velocity is stored in Yoffe_A
     *mp.csv'
      open (11,file='Yoffe_Slip.csv',status='unknown')
      write(11,*)'time(s),velocity,displacement,acceleration'
      yoffeSV(1)=0.d0
      write(11,1000)0.d0,yoffeSV(1),0.d0,0.d0
c
      yoffeS=0.d0
      izero=0
      do 100 i=2,nt
c
        yoffeSV(i)=0.d0
        t=dfloat(i-1)*dt
        if(izero.eq.1)goto 10
c
        if(t.le.TauS)then
          C1=(0.5d0*t+0.25d0*TauR)*dsqrt(t*(TauR-t))
          C1=C1+(t*TauR-TauR2)*dasin(dsqrt(t/TauR))
          C1=C1-0.75d0*TauR2*datan(dsqrt((TauR-t)/t))
          C2=0.375d0*pai*TauR2
          yoffeSV(i)=C1+C2
        elseif(t.le.2.d0*TauS)then
          C1=(0.5d0*t+0.25d0*TauR)*dsqrt(t*(TauR-t))
          C1=C1+(t*TauR-TauR2)*dasin(dsqrt(t/TauR))
          C1=C1-0.75d0*TauR2*datan(dsqrt((TauR-t)/t))
          C2=0.375d0*pai*TauR2
          C3=(TauS-t-0.5d0*TauR)*dsqrt((t-TauS)*(TauR-t+TauS))
          C3=C3+2.d0*TauR*(TauR-t+TauS)*dasin(dsqrt((t-TauS)/TauR))
          C3=C3+1.5d0*TauR2*datan(dsqrt((TauR-t+TauS)/(t-TauS)))
          yoffeSV(i)=C1-C2+C3
        elseif(t.le.TauR)then
          C1=(0.5d0*t+0.25d0*TauR)*dsqrt(t*(TauR-t))
          C1=C1+(t*TauR-TauR2)*dasin(dsqrt(t/TauR))
          C1=C1-0.75d0*TauR2*datan(dsqrt((TauR-t)/t))
          C3=(TauS-t-0.5d0*TauR)*dsqrt((t-TauS)*(TauR-t+TauS))
          C3=C3+2.d0*TauR*(TauR-t+TauS)*dasin(dsqrt((t-TauS)/TauR))
          C3=C3+1.5d0*TauR2*datan(dsqrt((TauR-t+TauS)/(t-TauS)))
          Tau2S=2.d0*TauS
          C4=(-TauS+0.5d0*t+0.25d0*TauR)*dsqrt((t-Tau2S)*(TauR-t+Tau2S))
          C4=C4+TauR*(t-TauR-Tau2S)*dasin(dsqrt((t-Tau2S)/TauR))
          C4=C4-0.75d0*TauR2*datan(dsqrt((TauR-t+Tau2S)/(t-Tau2S)))
          yoffeSV(i)=C1+C3+C4
        elseif(t.lt.(TauR+TauS))then
          C3=(TauS-t-0.5d0*TauR)*dsqrt((t-TauS)*(TauR-t+TauS))
          C3=C3+2.d0*TauR*(TauR-t+TauS)*dasin(dsqrt((t-TauS)/TauR))
          C3=C3+1.5d0*TauR2*datan(dsqrt((TauR-t+TauS)/(t-TauS)))
          Tau2S=2.d0*TauS
          C4=(-TauS+0.5d0*t+0.25d0*TauR)*dsqrt((t-Tau2S)*(TauR-t+Tau2S))
          C4=C4+TauR*(t-TauR-Tau2S)*dasin(dsqrt((t-Tau2S)/TauR))
          C4=C4-0.75d0*TauR2*datan(dsqrt((TauR-t+Tau2S)/(t-Tau2S)))
          C5=0.5d0*pai*TauR*(t-TauR)
          yoffeSV(i)=C3+C4+C5
        elseif(t.le.(TauR+2.d0*TauS))then
          C4=(-TauS+0.5d0*t+0.25d0*TauR)*dsqrt((t-Tau2S)*(TauR-t+Tau2S))
          C4=C4+TauR*(t-TauR-Tau2S)*dasin(dsqrt((t-Tau2S)/TauR))
          C4=C4-0.75d0*TauR2*datan(dsqrt((TauR-t+Tau2S)/(t-Tau2S)))
          C6=0.5d0*pai*TauR*(2.d0*TauS-t+TauR)
          yoffeSV(i)=C4+C6
          if(C4+C6.lt.0.0005)then
            izero=1
            yoffeSV(i)=0.d0
          end if
        endif
c
        yoffeSV(i)=yoffeSV(i)*Dmax*CK
        if(i.eq.2)then
          yoffeS=yoffeS+yoffeSV(i)*dt/2.d0
          yoffeA=yoffeSV(i)/dt
        else
          yoffeS=yoffeS+(yoffeSV(i)+yoffeSV(i-1))*dt/2.d0
          yoffeA=(yoffeSV(i)-yoffeSV(i-1))/dt
        end if    
   10   write(11,1000)t,yoffeSV(i),yoffeS,yoffeA
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 Yoffe: nfast > MOM !',nfast,MOM
c        write(10,*)'from Yoffe: nfast > MOM !',nfast,MOM
      end if
      do i=1,nt
        yoffeAmp(i)=dCMPLX(yoffeSV(i),0.0d0)
      end do
      do i=nt+1,nfast
        yoffeAmp(I)=dCMPLX(0.d0,0.0d0)
      end do
c
      CALL dWVFAST(nfast,yoffeAmp,MOM,-1)

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

