program ABF
      parameter( nmax=500 )
      real t(nmax+1), x(nmax+1), y(nmax+1)  ! x=位置, y=速度
      real dt, dx
      real x_0
      real k1, k2, k3, k4
      real j1, j2, j3, j4
      real tmin, tmax
      real omega
      pi=3.14159
*-- 出力ファイルをダンプ (保存) するための命令 ---
*      CALL SWpSET( 'LDUMP', .TRUE. )
*      CALL SWpSET( 'LWAIT', .FALSE. )
*      CALL SWpSET( 'LWAIT1', .FALSE. )
*-- 係数の入力方式についての指定 ---
      write(*,*) "input the frequency."
      read(*,*) omega
*-- 時間ステップの設定 ---
      dt=0.004
*-- 初期条件の設定 ---
      x(1)=1.0
      y(1)=0.0
*-- 初期値を求めるため Runge-Kutta 法 (2 次) を用いて計算しておく ---
         k1=y(1)*dt
         j1=(-x(1)*omwga**2)*dt
         k2=(y(1)+0.5*k1)*dt
         j2=(-(x(1)+0.5*j1)*omwga**2)*dt
         k3=(y(1)+0.5*k2)*dt
         j3=(-(x(1)+0.5*j2)*omwga**2)*dt
         k4=(y(1)+k3)*dt
         j4=(-(x(1)+j3)*omwga**2)*dt
         x(2)=x(1)+(1.0/6.0)*(k1+2.0*k2+2.0*k3+k4)
         y(2)=y(1)+(1.0/6.0)*(j1+2.0*j2+2.0*j3+j4)

*-- 実際の計算ここから ---
*-- 時間ステップの計算 ---
      t(1)=0.0
      t(2)=dt
      do 10 i=2,nmax
         t(i+1)=t(i)+dt
         x(i+1)=x(i)+dt*(1.5*y(i)-0.5*y(i-1))
         y(i+1)=y(i)-dt*(1.5*x(i)-0.5*x(i-1))*omega**2
  10  continue

*-- グラフの出力装置の指定 (IWS=4 は GTK 指定) ----
      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
*      IWS=4
      READ (*,*) IWS

*-- 出力装置の呼び出し
      CALL GROPN( IWS )

      CALL GRFRM
      CALL USSPNT( nmax+1, t, x )

      CALL GRSWND(  0.0,  2.0,  -2.0,  2.0 )
      CALL GRSVPT(  0.2,  0.8,  0.2,  0.8 )
      CALL USPFIT
      CALL GRSTRF

      CALL USSTTL( 'VELOCITY', '', 'POSITION', '' )
      CALL USDAXS

      CALL UULIN( nmax+1, t, x )
      CALL GRCLS

      stop
      end