!---------------------------------------------
!   2011/08/25 福士 邦彦 作成
!   2011/10/18 福士 邦彦 更新
!   2012/11/08 増田 和孝 最終更新
!   後退スキームによるdU/dt=iwUの解
!   初期条件U(0)=1,他の条件はw=pi
!   求めるのはRe[U(t)]
!   解析解はU(t)=Re{exp(iwt)}
!
!----------------------------------------------
program scheme_koutai
  implicit none
  integer :: loop                            ! 作業用変数
  integer :: step                            ! ステップ数
  real(8) :: time                            ! 積分終了時間
  real(8) :: dt                              ! 時間の間隔
  real(8),parameter :: pi = 4*atan(1.0d0)    ! 振動数, 今はパイにしている
  complex(8),parameter :: im = (0.0d0, 1.0d0)          ! 虚数単位
  real(8),allocatable :: u_kaiseki(:)        ! 解析解
  complex(8),allocatable :: u_suuchi(:)      ! 数値解
  real(8),allocatable :: u_suuchi_real(:)    ! 数値解の実部
  real(8),allocatable :: error(:)            ! 絶対誤差
  real(8),allocatable :: worktime(:)         ! グラフ用時間変数
  real(8) :: zoufuku                         ! 増幅係数
  real(8) :: zoufuku1                        ! 増幅係数解析値
  integer :: iws
  complex(8) :: u_before                      ! ループ計算用変数(前のステップでのu)
  complex(8) :: u_after                       ! ループ計算用変数(次のステップでのu)


  ! 終了時刻(無次元)を入力
  write(*,*) 'time = '
  read(*,*) time

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

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


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

  ! 初期値を設定する
  u_kaiseki(0) = cmplx(1.0d0, 0.0d0)
  error(0) = 0.0d0
  worktime(0) = 0.0d0
  u_before = cmplx(1.0d0, 0.0d0)
  u_after = cmplx(1.0d0, 0.0d0)
  u_suuchi_real(0) = dble(u_before)

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

  ! それぞれの値を外部ファイルへ書き出す準備
  open(10,file="data-koutai.dat") 
  write(10,*) '　step数　',' 解析解　','　数値解　','　絶対誤差 ','　増幅係数　','　増幅係数解析値'

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

     ! 次のステップのUを求める
     u_after = ((1.0d0+im*pi*dt)/(1.0d0+(pi*dt)**2.0d0))*u_before 
     u_kaiseki(loop+1) = cos(pi*dt*(loop+1))


     ! 増幅係数を計算
     zoufuku = abs(u_after)/abs(u_before)   
     ! グラフ用時間変数を計算
     worktime(loop+1) = worktime(loop) + dt

     ! 実部を計算
     u_suuchi_real(loop+1) = dble(u_after)

     ! 絶対誤差を計算
     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

     ! u の更新
     u_before = u_after

  enddo
  close(10)

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

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

  call GROPN(iws)
  call GRFRM

  call USSTTL( 'time', '', 'U','')
  call USSPNT(step,real(worktime),real(u_kaiseki))
  call USSPNT(step,real(worktime),real(u_suuchi_real))
  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
  deallocate (u_kaiseki,u_suuchi_real,error,worktime)

end program scheme_koutai
