!摩擦振動方程式
!  dU/dt = i\omega U - \kappa U
!  U(0) = 1.0
!  \omega = \pi, \kappa = 1/4
!振動方程式には leapfrog, 摩擦方程式には euler を用いる

      program leapeuler

!-----宣言文------------------------------------
      implicit none

      complex(8), allocatable :: U_a(:)            !U^{n+1}
      complex(8), allocatable :: U_b(:)            !U^{n-1}
      complex(8), allocatable :: U_c(:)            !U^{n}
      complex(8), allocatable :: Str(:)            !厳密解

      real(8), allocatable :: time(:)              !時間

      complex(8), parameter :: ei = (0.0d0,1.0d0)  !虚数単位
      real(8), parameter :: pi = 4.0d0*atan(1.0d0) !π
      real(8), parameter :: omega = pi             !振動数ω
      real(8), parameter :: kappa = 0.25d0         !摩擦係数κ
      real(8), parameter :: tmax = 5.0d0*pi        !時間終端
      real(8), parameter :: tmin = 0.0d0           !時間初期
      real(8) :: dt                                !時間間隔
      real, allocatable :: U_sing(:)               !数値解の単精度
      real, allocatable :: Str_sing(:)             !厳密解の単精度 
      real, allocatable :: time_sing(:)            !時間の単精度
      real, allocatable :: tang(:)                 !減衰解
      real, allocatable :: e_folding(:)            !e-folding time

      integer t_lupe                               !作業変数
      integer step                                 !分割数
      integer IWS                                  !dcl用作業変数

      write(*,*) "Enter in step number"

      read(*,*) step 

      ! 配列の動的振り付け
      allocate(U_a(0:step))
      allocate(U_b(0:step))
      allocate(U_c(0:step))
      allocate(time(0:step))
      allocate(U_sing(0:step))
      allocate(time_sing(0:step))
      allocate(Str(0:step))
      allocate(Str_sing(0:step))
      allocate(tang(0:step))
      allocate(e_folding(0:step))

      dt = (abs(tmax) - abs(tmin))/dble(step)

      write(*,*) "dt=", dt

!-----初期値代入---------------------------
      U_a(0) = (1.0d0,0.0d0)
      U_c(0) = (1.0d0,0.0d0)   
      U_b(0) = (1.0d0,0.0d0)
      time(0) = 0.0d0

      U_sing(0) = real(U_a(0))
      time_sing(0) = real(time(0))

     ! 厳密解の計算 u = u(0)exp((i\omega -\kappa)t)
      Str(0) = U_a(0)*exp((ei*omega -kappa)*time(0))
      Str_sing(0) = real(Str(0))

     ! 減衰解の計算 u = u(0)exp(-\kappa t)
      tang(0) = real(U_a(0))*exp(-kappa*time(0))
      e_folding(:) = 1.0e0/exp(1.0e0)

!------------数値積分(n=0)----------------------------------------
!     リープフロッグスキームを用いるため1段階数スキームを用いてu^{2}の値を求める.
!     用いるスキームはホインスキーム

      U_a(0) = U_b(0) + ei*omega*dt*U_b(0)-kappa*dt*U_b(0)
      U_c(0) = U_b(0) + 0.5d0*dt*(&
        & ei*omega*dt*U_b(0)-kappa*dt*U_b(0) &
        &  + ei*omega*dt*U_a(0)-kappa*dt*U_a(0))

!------------数値積分(n=1~N)------------------------
!      n=1 以降の計算
!      振動解の方にはリープフロッグスキーム,
!      減衰解の方にはオイラースキームを用いる

        do t_lupe = 1, step

         time(t_lupe) = dble(t_lupe)*dt

         ! リープフロッグスキーム
         U_a(t_lupe) = U_b(t_lupe-1) + 2.0d0*dt*(&
           & ei*omega*U_c(t_lupe-1) - kappa*U_b(t_lupe-1))

         U_b(t_lupe) = U_c(t_lupe-1)
         U_c(t_lupe) = U_a(t_lupe)

         ! 厳密解の計算
         Str(t_lupe) = U_a(0)*exp((ei*omega -kappa)*time(t_lupe))

         ! 単精度実数にする
         U_sing(t_lupe) = real(U_a(t_lupe))
         time_sing(t_lupe) = real(time(t_lupe))
         Str_sing(t_lupe) = real(Str(t_lupe))
         tang(t_lupe) = real(U_a(0))*exp(-kappa*time(t_lupe))

      end do

!---------dcl----------------------------------------------
      call sgpwsn
      read(*,*) IWS

      call gropn(IWS)
      call grfrm

      call usspnt(step+1,time_sing(0:step),U_sing(0:step))
      call usspnt(step+1,time_sing(0:step),Str_sing(0:step))
      call usspnt(step+1,time_sing(0:step),e_folding(0:step))
      call usspnt(step+1,time_sing(0:step),tang(0:step))
      call uspfit
      call grstrf
      call usdaxs

      ! 数値解(実線)
      call uulin(step+1,time_sing(0:step),U_sing(0:step))

      ! 厳密解(破線)
      call uuslnt(2)
      call uulin(step+1,time_sing(0:step),Str_sing(0:step))

      ! e-folding time(点線)
      call uuslnt(3)
      call uulin(step+1,time_sing(0:step),e_folding(0:step))

      ! 減衰解(1点鎖線)
      call uuslnt(4)
      call uulin(step+1,time_sing(0:step),tang(0:step))

      call grcls

      ! 配列の開放
      deallocate(U_a, U_b, U_c, time, U_sing,&
        &  time_sing, Str, Str_sing, tang, e_folding)

      stop
      end program leapeuler
