      program ex053b
      implicit real*8 (a-h,o-z)
c
c Set default parameters
      n = 1000           ! mesh points
      gam = 0.1d0        ! gamma, Friction (Masatsu)
      omg  = 2.0d0       ! omega, Restoring Force
      ti = 0.0d0         ! Initial Value of t
      tf = 100.0d0       ! Final   Value of t
      xi = 5.0d0         ! Initial  Value of x
      vi = 0.0d0         ! Initial  Value of v
C
c Parameter inputs
      write(*,810) '# gam, omg =?'
      read(*,*)       gam, omg
c
      imethod=4
      write(*,810) '# Method = ?'
      write(*,810) '# 1: Euler'
      write(*,810) '# 2: Euler (Improved)'
      write(*,810) '# 3: Euler (Modified)'
      write(*,810) '# 4: Heun'
      read(*,*) imethod
 810  format(a)                  ! format(a): $BJ8;zNs$r6uGr$r3+$1$:$K=q$/(B
      open(16,file='ex053b.dat') ! open(unit,file='filename')
c
      write( *,820) '# gam, omg, imethod =',gam,omg,imethod
      write(16,820) '# gam, omg, imethod =',gam,omg,imethod
c
 820  format(a,2(1x,f15.7),1x,i2) ! $BJ8;z(B 2$B2s(B($B6uGr(B $B<B?t$r(B15$B7e$G(B) $B6uGr(B $B@0?t$r(B2$B7e$G(B
c
      dt = (tf-ti)/n     ! Time Mesh
      t=ti               ! Initial Condition
      x=xi               ! Initial Condition
      v=vi               ! Initial Condition
c
         call exact(xexa,t,xi,vi,gam,omg)
         write( *,800) t,x,v,xexa,abs(x-xexa)
         write(16,800) t,x,v,xexa,abs(x-xexa)
      do i = 1,n
c...
c Select one of the following methods.
            if(imethod.eq.1) then
               call euler(x,v,t,dt,gam,omg)
            else if(imethod.eq.2) then
               call impeul(x,v,t,dt,gam,omg)
            else if(imethod.eq.3) then
               call modeul(x,v,t,dt,gam,omg)
            else if(imethod.eq.4) then
               call heun(x,v,t,dt,gam,omg)
            endif
c...Exact Solution
               call exact(xexa,t,xi,vi,gam,omg)
c...
         write( *,800) t,x,v,xexa,abs(x-xexa)
         write(16,800) t,x,v,xexa,abs(x-xexa)
      end do
c
      stop
 800  format(5(1x,f15.7))
      end
c *********************************************************************
      subroutine impeul(x,v,t,dt,gam,omg)
c *********************************************************************
      implicit real*8 (a-h,o-z)
c 
      v0 = v
      f0 = f(v,x,t,gam,omg)
c
      t1 = t + dt
      x1 = x + dt*v0
      v1 = v + dt*f0
      f1 = f(v1,x1,t1,gam,omg)
c
      dx = ( v0 + v1 )*dt/2	! Improved Euler
      dv = ( f0 + f1 )*dt/2	! Improved Euler
c
      t = t + dt
      x = x + dx
      v = v + dv
c
      end
c *********************************************************************
      subroutine euler(x,v,t,dt,gam,omg)
c *********************************************************************
      implicit real*8 (a-h,o-z)
c
      dx = v                * dt ! Euler
      dv = f(v,x,t,gam,omg) * dt ! Euler
      t = t + dt
      x = x + dx
      v = v + dv

      end
c *********************************************************************
      subroutine modeul(x,v,t,dt,gam,omg)
c *********************************************************************
      implicit real*8 (a-h,o-z)
c 
      v0 = v
      f0 = f(v,x,t,gam,omg)
c
      t1 = t + dt   /2
      x1 = x + dt*v0/2
      v1 = v + dt*f0/2
      f1 = f(v1,x1,t1,gam,omg)
c
      dx = v1*dt	! Modified Euler
      dv = f1*dt	! Modified Euler
c
      t = t + dt
      x = x + dx
      v = v + dv
c
      end
c *********************************************************************
      subroutine heun(x,v,t,dt,gam,omg)
c *********************************************************************
*     second order runge-kutta method
*        z(t+h) = z(t) + h*{ c1*f[z(t)]+c2*f(z(t)+c3*h*f(z(t))) }
*     where c1=1-c2, c3=1.0/2/c2, 0<c2<1
*     before call, dmat must be called
*
*     for practical use
*        z1 = z0+c3*h*f0
*        z2 = z0+h*c1*f0+h*c2*f1
*           = z1+h*(1-c2)*f0+h*c2*f1
*
*     note:c2 must be 0< c2 <=1
*       c2=0.5 => improved Euler method i
*       c2=1   => modified Euler method ii
*       c2=0.75=> Heun method
*
************************************************************************
      implicit real*8(a-h,o-z)
c 
      c2=0.75d0      ! Heun
c     c2=0.5d0       ! improved Euler
c     c2=1.0d0       ! modified Euler
      c1=1.0d0-c2
      c3=1.0d0/2/c2
c 
      v0 = v
      f0 = f(v,x,t,gam,omg)
c
      t1 = t + dt   *c3
      x1 = x + dt*v0*c3
      v1 = v + dt*f0*c3
      f1 = f(v1,x1,t1,gam,omg)
c
      dx = ( c1*v0 + c2*v1 )*dt  ! Heun
      dv = ( c1*f0 + c2*f1 )*dt  ! Heun
      t = t + dt
      x = x + dx
      v = v + dv

      end
c *********************************************************************
      subroutine exact(x,t,xi,vi,gam,omg)
c *********************************************************************
      implicit real*8(a-h,o-z)
c
      pi=4*atan(1.0d0)
c
c Special Solution with External force
      omgf=1.0d0
      ff=1.0d0
      ek=gam/2
      ampf=ff/sqrt((omg**2-omgf**2)**2+4*ek*ek*omgf*omgf)
      deltaf=atan2(2*omgf*ek,omg**2-omgf**2)
      if(deltaf.lt.0) deltaf=deltaf+pi
c...
      x=ampf*cos(omgf*t-deltaf)
      x0f=ampf*cos(deltaf)
      v0f=ampf*sin(deltaf)*omgf
c
      if(omg.gt.ek) then
         omg1=sqrt(omg**2-ek**2)
         amp1=xi-x0f
         amp2=(ek*(xi-x0f)+(vi-v0f))/omg1
         x=x+exp(-ek*t)*(amp1*cos(omg1*t)+amp2*sin(omg1*t))
      else
         ek0=sqrt(ek**2-omg**2)
         ek1=ek+ek0
         ek2=ek-ek0
         amp1=(-ek2*(xi-x0f)-(vi-v0f))/2/ek0
         amp2=( ek1*(xi-x0f)+(vi-v0f))/2/ek0
         x=x+amp1*exp(-ek1*t)+amp2*exp(-ek2*t)
      endif

      end
c **********************************************************************
      function f(v,x,t,gam,omg)
c **********************************************************************
      implicit real*8 (a-h,o-z)
c Force
      f = -gam*v-omg*omg*x+cos(t)
      end
