!---------------------------------------------
!   2011/08/25 作成
!   2011/10/18 最終更新
!   台形スキームによるdU/dt=iwUの解
!   初期条件U(0)=1,他の条件はw=pi
!   求めるのはRe[U(t)]
!   解析解はU(t)=Re{exp(iwt)}
!
!----------------------------------------------
program scheme_daikei
  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

  ! 積分終了時間を入力
  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

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

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

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

     ! 次のステップのUを求める
     u_suuchi(loop+1) = ((1.0e0-0.25e0*(pi*dt)**2.0e0+im*pi*dt)/&
&                          (1.0e0+0.25e0*(pi*dt)**2.0e0))*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( 'worktime', 'second', 'U','y-unit')
  call USSPNT(step,worktime,u_suuchi_real)
  call USSPNT(step,worktime,u_kaiseki)
  call USPFIT
  call GRSTRF

  call usdaxs

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

  call GRCLS

end program scheme_daikei
