! 4 段 4 位の Runge-Kutta 公式を用いて振動方程式を数値的に解く
! du/dt = iwu
! u(0) = 1, \omega = \pi
! u = exp(i \pi t) 

program runge

!暗黙の型宣言禁止

implicit none


!変数定義

real(8) , allocatable :: Ur(:) !U(t) の実部の before
real(8) , allocatable :: Ui(:) !U(t) の虚部の before
real(8) , allocatable :: time(:)   !時間配列
real(8) , parameter :: omega = 4.0d0*atan(1.0d0)   !振動数
real(8) :: dt   !時間間隔
real(8) :: kei_a_r   !k1
real(8) :: kei_a_i   !k1
real(8) :: kei_b_r   !k2
real(8) :: kei_b_i   !k2
real(8) :: kei_c_r   !k3
real(8) :: kei_c_i   !k3
real(8) :: kei_d_r   !k4
real(8) :: kei_d_i   !k4
real(8) :: zouhuku
real(8) :: lambda
real(8) , allocatable :: true(:)
real(8) , allocatable :: dep(:)
integer :: unit ! 装置番号
integer :: nmax !ステップ数
integer :: it   !do 作業用変数
integer :: ins

nmax = 1000
dt = 0.010d0
!配列割付
allocate (Ur(1:nmax+1), Ui(1:nmax+1), time(1:nmax+1), true(1:nmax+1), dep(1:nmax+1) )

!値を書き出す
unit = 10
open(unit, file = 'out.d')

!dt = 0.10d0
!配列初期化
do it = 1, nmax+1
   Ur(it) = 0.0d0
   Ui(it) = 0.0d0
   time(it) = 0.0d0
   true(it) = 0.0d0
   dep(it) = 0.0d0
end do
do it = 1, nmax
   time(it+1) = time(it) + dt
end do

!初期値生成
!U(0)=1 から, U(1) を求める
Ur(1) = 1.0d0
Ui(1) = 0.0d0

do it = 1, nmax
   ! k1 を求める
   kei_a_r = - dt * omega * Ui(it)
   kei_a_i = dt * omega * Ur(it)
   
   ! k2 をもとめる
   kei_b_r = -dt * omega * ( Ui(it) + 0.50d0 * kei_a_i )
   kei_b_i = dt * omega * ( Ur(it) + 0.50d0 * kei_a_r )

   ! k3 をもとめる
   kei_c_r = -dt * omega * ( Ui(it) + 0.50d0 * kei_b_i )
   kei_c_i = dt * omega * (Ur(it) + 0.50d0 * kei_b_r )
   ! k4 を求める
   kei_d_r = -dt * omega * ( Ui(it) + kei_c_i )
   kei_d_i = dt * omega * ( Ur(it) + kei_c_r)

   ! U をもとめる
   Ur(it+1) = Ur(it) + (1.0d0 / 6.0d0) *( kei_a_r + 2.0d0 * kei_b_r + 2.0d0 * kei_c_r + kei_d_r)
   Ui(it+1) = Ui(it) + (1.0d0 / 6.0d0) *( kei_a_i + 2.0d0 * kei_b_i + 2.0d0 * kei_c_i + kei_d_i)
   time(it+1) = time(it) + dt

   ! 増幅係数

   write(unit,*) Ur(it+1)

   ! 解析解
   true(it) = cos(omega*time(it))

   ! 解析解と数値解の差
   dep(it) = abs(true(it) - dep(it))

end do

!増幅係数の理論値
zouhuku = sqrt(1.0d0 + (1.0d0 / 576.0d0) * (omega * dt)**8.0d0 - (1.0d0/72.0d0) * (omega*dt)**6.0d0)
write(*,*) "zouhuku=", zouhuku 

!output ファイルクローズ
close(unit)

!DCL 開始
WRITE(*,*) ' WORKSTATION ID (I)  ? ;'  
    CALL SGPWSN  
    READ (*,*) ins  
  
  
    CALL GROPN( ins )  
!    CALL GROPN( 2 )  
    CALL GRFRM  
  
 
    CALL USSTTL( 'time', '', 'Ur',  '')  
    CALL USGRPH( nmax+1, real(time), real(Ur) )  
  
   CALL GRCLS


deallocate (Ur, Ui, time)

end program runge
