      PROGRAM sho                           ! 調和振動子
      IMPLICIT NONE
      REAL*8 dt, dt2, t, v, w2, x
      INTEGER counter, nshow, kbhit
      LOGICAL Ltmp1_
      CALL initial(x,v,w2,t,dt,dt2,nshow)
      counter = 0
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL update(x,v,w2,t,dt,dt2)
         counter = counter + 1
         IF(counter .EQ. nshow) THEN
            CALL output(x,v,t)
            counter = 0
         END IF
         Ltmp1_ = kbhit() .EQ. 0
      END DO
      END
C
      SUBROUTINE initial(x0,v,w2,t,dt,dt2,nshow)
      IMPLICIT NONE
      REAL*8 show_time, x0, v, w2, t, dt, dt2
      INTEGER nshow
      WRITE(*,'(A,$)') 'initial position = '
      READ(*,*) x0                          ! メートル
      t = 0
      v = 0                                 ! 初速度
*     もとの (角) 振動数
      WRITE(*,'(A,$)') 'ratio of k/m = '
      READ(*,*) w2
      WRITE(*,'(A,$)') 'time step = '
      READ(*,*) dt
      dt2 = 0.5D0*dt
      show_time = 0.1D0                     ! 出力の時間間隔
      nshow = int(show_time/dt)
      WRITE(*,'(A)') '         time     position     velocity'
      WRITE(*,*)                            ! 空白行
      END
C
      SUBROUTINE update(x,v,w2,t,dt,dt2)
*     オイラー-リチャードソンのアルゴリズム
      IMPLICIT NONE
      REAL*8 a, amid, vmid, xmid, x, v, w2, t, dt, dt2
      a = -w2*x
      vmid = v + a*dt2                      ! 中点での速度
      xmid = x + v*dt2                      ! 中点の位置
      amid = -w2*xmid                       ! 中点での加速度
      v = v + amid*dt
      x = x + vmid*dt
      t = t + dt
      END
C
      SUBROUTINE output(x,v,t)
*     - 上位桁に空白か負号を付けて数を表示
*     % 上位桁の零には`0'を表示
      IMPLICIT NONE
      REAL*8 x, v, t
      WRITE(*, '(3F13.4)') t, x, v
      END
