!---------------------------------------------
!   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_euler
  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(:)      ! 隗｣譫占ｧ｣
  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)

  ! 縺昴ｌ縺槭ｌ縺ｮ蛟､繧貞､夜Κ繝輔ぃ繧､繝ｫ縺ｸ譖ｸ縺榊�縺呎ｺ門ｙ
  open(10,file="data-euler.dat")
  write(10,*) '縲step謨ｰ縲',' 隗｣譫占ｧ｣縲','縲謨ｰ蛟､隗｣縲','縲邨ｶ蟇ｾ隱､蟾ｮ ','縲蠅怜ｹ�ｿよ焚縲','縲蠅怜ｹ�ｿよ焚隗｣譫仙､'

  ! 繝ｫ繝ｼ繝苓ｨ育ｮ�-------------------------------------------
  do loop = 0,step-1,1

     ! 谺｡縺ｮ繧ｹ繝�ャ繝励�U繧呈ｱゅａ繧�
     u_after = (1.0d0+im*pi*dt)*u_before

     ! u_after縺ｨ蜷後§譎る俣縺ｧ縺ｮ隗｣譫占ｧ｣繧呈ｱゅａ繧�
     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_euler