! ホインの 3 次式を用いて振動方程式を数値的に解く
! du/dt = iwu
! u(0) = 1, \omega = \pi
! u = exp(i \pi t) 


program heun

!暗黙の型宣言禁止

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) :: keisuA   !
real(8) :: keisuB   !
real(8) :: zouhuku
real(8) :: lambda
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) )
unit = 10
open(unit, file = 'lambda.d')

!配列初期化
do it = 1, nmax+1
   Ur(it) = 0.0d0
   Ui(it) = 0.0d0
   time(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) + (1.0d0 / 3.0d0) * kei_a_i )
   kei_b_i = dt * omega * ( Ur(it) + (1.0d0 / 3.0d0) * kei_a_r )

   ! k3 をもとめる
   kei_c_r = -dt * omega * ( Ui(it) + (2.0d0 / 3.0d0) * kei_b_i )
   kei_c_i = dt * omega * (Ur(it) + (2.0d0 / 3.0d0) * kei_b_r )

   ! U をもとめる
   Ur(it+1) = Ur(it) + 0.250d0 *( kei_a_r + 3.0d0 * kei_c_r )
   Ui(it+1) = Ui(it) + 0.250d0 *( kei_a_i + 3.0d0 * kei_c_i )
   time(it+1) = time(it) + dt

   ! 増幅係数

   lambda = abs( Ur(it+1) ) / abs( Ur(it) )
   write(unit,*) lambda

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 
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 heun
