c
c *** Regularized Slip Velocity Function using Nakamura&Miyatake Function ***
c
      program miyakate
c
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(MOM=16400)
      complex*16 svAmp(MOM)
      dimension svWave(MOM),tm(MOM)

      dt=0.005d0
      nt=4096
      fmax=6
c	vm=20.837d0
c	slip=5.092d0
      vm=4.092d0
      slip=1.0d0
      tr=2.114d0
      ts=tr*1.5d0
      pai=3.141592654d0
      td=1.0d0/(pai*fmax)

      write(*,*)
      write(*,*)'dt,nt = ',dt,nt
      write(*,*)'fmax,vm,slip = ',fmax,vm,slip
      write(*,*)'td,tr,ts=',td,tr,ts
c
      do i=1,nt
        tm(i)=dt*dfloat(i-1)
      end do
      call vslip(tm,svWave,MOM,nt,slip,td,vm,Am,tr,ts,tb,b,c,ar,eps)
      write(*,*)'Am,tb,b,c,ar,eps =',Am,tb,b,c,ar,eps

      write(*,*)
      write(*,*)' Slip functions are stored in N&M_Slip.csv'
      write(*,*)' Fourier spectrum is stored in N&M_Amp.csv'
c
      open (10,file='N&M_Slip.csv',status='unknown')
        write(10,*)'time(s),velocity,displacement,acceleration'
	d0=0.0d0
	a0=0.0d0
	write(10,1000)0.d0,svWave(1),d0,a0
	do i=2,nt
        t=dfloat(i-1)*dt
	d1=d0+(svWave(i)+svWave(i-1))*dt/2.d0
	a1=(svWave(i)-svWave(i-1))/dt
	write(10,1000)t,svWave(i),d1,a1
	d0=d1
	end do
      close(10)
 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 Ramp: nfast > MOM !',nfast,MOM
c        write(10,*)'from Ramp: 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='N&M_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)
        svAmp(io)=dconjg(svAmp(io))*dt
c        write(12,1100)freq,cdabs(svAmp(io)),svAmp(io)
        write(12,1100)freq,cdabs(svAmp(io)),cdabs(svAmp(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 *** subroutine from Nakamura and Miyatake (2000) ***
c     (http://www.eri.u-tokyo.ac.jp/miyatake/SlipFunc-Prog.html)
c
      subroutine vslip(tm,rate,ndt,nt,slip,td,vm,Am,tr,ts,tb,b,c,ar,eps)
c
      IMPLICIT REAL*8(A-H,O-Z)
      dimension rate(ndt),tm(ndt)
      error = 1.0d-5
      Am = 2.d0 * Vm / td
      isw = 1
      h = td
      x = td*1.01d0
      call tslip(y1,slip,x ,Am,td,tr,ts,b,c,ar,eps)
    1 call tslip(y2,slip,x+h,Am,td,tr,ts,b,c,ar,eps)
      if(y1*y2.gt.0.d0) then
      x = x + h
      y1 = y2
      if(isw.eq.1) goto 1
      endif
      h = h / 2.d0
      isw = 2
      if(h.gt.error) goto 1
      tb = x
      do k=1,nt
      if(tm(k).le.tb) then
      v = Am * tm(k) * ( 1.d0 - 0.5d0 * tm(k) / td )
      endif
      if(tm(k).gt.tb.and.tm(k).le.tr) then
      v = b / sqrt( tm(k) - eps)
      endif
      if(tm(k).gt.tr.and.tm(k).le.ts) then
      v = c - ar * ( tm(k) - tr )
      endif
      if(tm(k).gt.ts) v = 0
      if(tm(k).lt. 0.d0) v = 0
      rate(k) = v
      enddo
      return
      end
c
      subroutine tslip(res,target,tb,Am,td,tr,ts,b,c,ar,eps)
c
      IMPLICIT REAL*8(A-H,O-Z)
      eps = tb + tb / 2.d0 * (1.d0-tb/td/2.d0)/(1.d0-tb/td)
      if(eps.gt.tb) eps = tb
      b = -2.d0 * Am * ( 1.d0 - tb / td ) * ( tb - eps ) ** 1.5d0
      c = b / dsqrt(tr-eps)
      ar = c / ( ts - tr )
      f1 = Am / 2.d0 * tb * tb - Am / 6.d0 * (tb**3) / td
      f2 = 2.d0 * b * ( dsqrt(tr-eps) - dsqrt(tb-eps) )
      f3 = c / 2.d0 * ( ts - tr )
      res = f1 + f2 + f3 - target
      return
      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

