!鞫ｩ謫ｦ謖ｯ蜍墓婿遞句ｼ�
!  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