program shindou_ab
 ! 振動方程式をアダムスバッシュフォース法で解いてみる

!!!!!!!!!!!!!!!型宣言!!!!!!!!!!!!!!!
 implicit none
 real(8) U, Ure_n, Uim_n, Ure_np1, Uim_np1, Ure_nm1, Uim_nm1! real&image,n-1&n&n+1の計6つの U を用いる 
 real(8) delta_t, omega, genmitukai, time, pi, lambda 
 integer n
 
!!!!!!!!!!!!!!!定数設定!!!!!!!!!!!!!!!
 pi=4.0d0*atan(1.0d0)
 omega = pi
 Ure_n = 1.0d0 !初期条件
 Uim_n = 0.0d0 !初期条件
 delta_t = 0.01d0 

 open(50,file='shindou_abd001.out')

!!!!!!!!!!!!!!!計算!!!!!!!!!!!!!!!
!最初はオイラースキームで2項目まで計算
 Ure_np1 = Ure_n - ( omega * delta_t * Uim_n )
 Uim_np1 = Uim_n + ( omega * delta_t * Ure_n )

 Ure_nm1 = Ure_n
 Ure_n = Ure_np1
 Uim_nm1 = Uim_n
 Uim_n = Uim_np1
!ここからABスキーム
 do n = 2,400
  Ure_np1 = Ure_n - ( (1.5d0) * omega * delta_t * Uim_n ) - ( (0.5d0) * omega * delta_t * Uim_nm1)
  Uim_np1 = Uim_n + ( (1.5d0) * omega * delta_t * Ure_n ) - ( (0.5d0) * omega * delta_t * Ure_nm1)

  lambda = Ure_n / Ure_nm1

  Ure_nm1 = Ure_n
  Uim_nm1 = Uim_n
  Ure_n = Ure_np1
  Uim_n = Uim_np1

  genmitukai = cos( - omega * n * delta_t)
  time = n * delta_t
  write(50,*) time, Ure_n, genmitukai, lambda
 enddo

 close(50)

end program shindou_ab
