!-------------------------------------------
!  2011/09/29 作成
!  2011/10/18 最終更新
!  リープフロッグスキームによるdU/dt=iwUの解
!  初期条件はU(0)=1,他の条件はw=pi
!  求めるのはRe[U(t)]
!  解析解はU(t)=Re{exp(iwt)}
!
!--------------------------------------------

program scheme_leapfrog
  implicit none
  integer :: loop                          ! 作業用変数
  integer :: step                          ! ステップ数
  real(4) :: time                          ! 積分終了時間
  real(4) :: dt                            ! 時間の間隔
  real(4),parameter :: pi=4*atan(1.0e0)    ! 振動数,今はパイにしている
  complex(4) :: im = (0.0e0, 1.0e0)        ! 虚数単位
  real(4),allocatable :: u_kaiseki(:)      ! 解析解
  complex(4),allocatable :: u_suuchi(:)    ! 数値解
  real(4),allocatable :: u_suuchi_real(:)  ! 数値解の実部
  real(4),allocatable :: error(:)          ! 絶対誤差
  real(4),allocatable :: worktime(:)       ! 作業用時間変数
  real(4) :: zoufuku                       ! 増幅係数
  real(4) :: zoufuku1                      ! 増幅係数解析値
  integer :: iws
  integer :: ins

  ! 積分の終わりの時間を入力
  write(*,*) 'time = '
  read(*,*) time

  ! ステップ数を入力
  write(*,*) 'step = '
  read(*,*) step

  ! 時間の刻み幅を入力
  dt = time/real(step)

  ! 配列の形を決める
  allocate (u_kaiseki(0:step),u_suuchi(0:step),u_suuchi_real(0:step),error(0:step),worktime(0:step))

  ! 初期値を設定する
  u_kaiseki(0) = cmplx(1.0e0, 0.0e0)
  u_suuchi(0) = cmplx(1.0e0, 0.0e0)
  u_suuchi_real(0) = real(u_suuchi(0))
  error(0) = 0.0e0
  worktime(0) = 0.0e0

  ! 2番目(配列の(1)に相当するもの)を計算する
  ! eulerスキームの場合
  ! u_suuchii(1) = (1.0e0 + im*pi*dt)*u_suuchi(0)           

  ! 後退スキームの場合
  ! u_suuchi(1) = (1.0e0 + im*pi*dt)/(1 + (pi*dt)**2)*u_suuchi(0)

  !　台形スキームの場合
  u_suuchi(1) = (1.0e0 - 0.25e0*(pi*dt)**2.0e0 + im*pi*dt)/&
&                (1.0e0 + 0.25e0*(pi*dt)**2.0e0)

  ! 1ステップ目の数値解の実部をとる
  u_suuchi_real(1) = real(u_suuchi(1))

  ! 1ステップ目の解析解を計算
  u_kaiseki(1) = cos(pi*dt)

  ! 1ステップ目の作業時間の計算
  worktime(1) = worktime(0) + dt

  ! 1ステップ目の増幅係数の計算
  zoufuku = abs(u_suuchi(1))/abs(u_suuchi(0))

  ! 1ステップ目の絶対誤差の計算
  error(1) = abs(u_kaiseki(1)) - abs(u_suuchi(1))

  ! 増幅係数解析値の計算
  zoufuku1 = abs(sqrt(1.0e0 - (pi*dt)**2.0e0) + im*pi*dt)

  ! それぞれの値をファイルに書き出す
  open(10,file='data-leap.dat')
  write(10,*) '回数, 解析解, 数値解, 絶対誤差, 増幅係数, 増幅係数解析値'
  write(10,*) '0',u_kaiseki(0),u_suuchi_real(0)
  write(10,*) '1',u_kaiseki(1),u_suuchi_real(1),error(1),zoufuku,zoufuku1

  ! ループ計算-------------------------------------------------
  do loop = 1,step-1,1

     ! 次のステップのUを求める
     u_suuchi(loop+1) = u_suuchi(loop-1) + 2.0e0*im*pi*dt*u_suuchi(loop)

     ! u_suuchi(loop+1)と同じ時間での解析解を求める
     u_kaiseki(loop+1) = cos(pi*dt*(loop+1))

     ! 増幅係数を計算
     zoufuku = abs(u_suuchi(loop+1))/abs(u_suuchi(loop))

     ! グラフ用時間変数を計算
     worktime(loop+1) = worktime(loop)+dt

     ! 実部を計算
     u_suuchi_real(loop+1) = real(u_suuchi(loop+1))

     ! 絶対誤差を計算
     error(loop+1) = u_suuchi_real(loop+1) - u_kaiseki(loop+1)


     ! 書き出し
    write(10,*) loop+1,u_kaiseki(loop+1),u_suuchi_real(loop+1),error(loop+1),zoufuku,zoufuku1

  enddo
  close(10)

!------------------------------------------------------

  write(*,*) 'WORKSTATION ID(I) ?'
  call SGPWSN
  read(*,*) iws

  call GROPN(iws)
  call GRFRM

  call USSTTL( 'time', 'second', 'U','y-unit')
  call USSPNT(step,real(worktime),real(u_suuchi_real))
  call USSPNT(step,real(worktime),real(u_kaiseki))
  call USPFIT
  call GRSTRF

  call usdaxs

  call UULIN(step,real(worktime),real(u_suuchi_real))
  
  call UUSLNT(2)
  call UULIN(step,real(worktime),real(u_kaiseki))

  call GRCLS

end program scheme_leapfrog
