program main
!
! 振動方程式の数値計算
! 20101104 川畑拓也
!
! dU/dt + iwU = 0
! 厳密解は U = e^{-iwt}
! 簡単のため w = 
! 解くべき式は
! Re[dU/dt + iwU = 0]
! Im[dU/dt + iwU = 0]
! 最終的に求める解は Re[U(t)] である

implicit none
real(4)    :: Urn1 ! 次のステップの U の実部
real(4)    :: Uin1 ! 次のステップの U の虚部
real(4)    :: Urn  ! 今のステップの U の実部
real(4)    :: Uin  ! 今のステップの U の虚部
real(4)    :: U    ! 厳密解
real(4)    :: dt   ! 時間のインターバル
real(4)   , parameter :: w = 4 * atan(1.0e0) ! 振動数. 今は パイに設定
integer :: n    ! 時間のステップ数
integer :: i    ! do 作業用整数
integer :: k
integer :: l
real(4) , allocatable  :: a(:)
real(4) , allocatable :: b(:)
integer :: ins

! 時間の間隔を入力
write(*,*) 'interval dt?'
read(*,*) dt
! 分割数を入力
write(*,*) 'n = '
read(*,*) n

allocate ( a(1:n), b(1:n) )

! 初期値の設定
Urn = 1.0e0
Uin = 0.0e0
Urn1= 0.0e0
Uin1= 0.0e0


! 実際に計算
do i = 1, n, 1
   n = real(n)
   k = real(i)
   Urn1 = Urn + w*dt*Uin
   Uin1 = Uin - w*dt*Urn
   Urn = Urn1
   Uin = Uin1
   U = cos( w*k*dt)
   a(i) = i * dt
   b(i) = Urn
 !  write(*,*) 'n = ', i, 'Urn = ', Urn, 'U = ', U
end do

WRITE(*,*) ' WORKSTATION ID (I)  ? ;'  
    CALL SGPWSN  
!    READ (*,*) ins  


! PS 形式で図を出力  
  
!    CALL GROPN( ins )  
    CALL GROPN( 2 )  
    CALL GRFRM  
  
 
    CALL USSTTL( 'X-TITLE', 'x-unit', 'Y-TITLE', 'y-unit' )  
    CALL USGRPH( n, a, b )  
  
   CALL GRCLS
end program main
