program force
      parameter( nmax=500 )
      real t(nmax+1), x(nmax+1), y(nmax+1), z(nmax+1)
      real dt, dx
      real x_0
      real k1, k2, k3, k4
      real j1, j2, j3, j4
      real tmin, tmax
      real alpha, omega, pi, delta, domega
      integer IW
      pi=3.14159
*-- 出力ファイルをダンプ (保存) するための命令 ---
*      CALL SWpSET( 'LDUMP', .TRUE. )
*      CALL SWpSET( 'LWAIT', .FALSE. )
*      CALL SWpSET( 'LWAIT1', .FALSE. )
*-- 係数の入力方式についての指定 ---
      write(*,*) "Which do you select?"
      write(*,*) "0: type, 1: pretype"
      read(*,*) IW
*-- 係数の決定 ---
      if(IW.eq.0)then
         write(*,*) "input the parameters"
         write(*,*) "omega(frequency)=, alpha(forcing F)=
     & , fo(forcing_amplitude)="
         read(*,*) omega, alpha, fo
      else
         omega=1.0
         alpha=1.5
         fo=1.0
      end if
      omega=omega**2
*-- 時間ステップの設定 ---
      tmin=0.0
      tmax=10.0
      dt=(tmax-tmin)/(nmax-1)
*-- 初期条件の設定 ---
      t(1)=0.0
      if(IW.eq.0)then
         write(*,*) "input the initialized numbers"
         write(*,*) "x_0(init position)=, u_0(init velocity)="
         read(*,*) x_0, y_0
         x(1)=x_0
         y(1)=y_0
         z(1)=x(1)
      else
         x(1)=0.0
         y(1)=1.0
         z(1)=x(1)
      end if
      domega=omega**2
*-- 実際の計算ここから ---
      do 10 i=1,nmax
*-- k の計算 (dx/dt=y で計算) ---
         t(i+1)=t(i)+dt
         k1=y(i)*dt
         j1=(-domega*x(i)+fo*cos(alpha*t(i)))*dt
         k2=(y(i)+0.5*j1)*dt
         j2=(-domega*(x(i)+0.5*k1)
     &      +fo*cos(alpha*(t(i)+0.5*dt)))*dt
         k3=(y(i)+0.5*j2)*dt
         j3=(-domega*(x(i)+0.5*k2)
     &      +fo*cos(alpha*(t(i)+0.5*dt)))*dt
         k4=(y(i)+j3)*dt
         j4=(-domega*(x(i)+k3)
     &      +fo*cos(alpha*(t(i)+dt)))*dt
*-- 時間ステップの計算 ---
         x(i+1)=x(i)+(1.0/6.0)*(k1+2.0*k2+2.0*k3+k4)
         y(i+1)=y(i)+(1.0/6.0)*(j1+2.0*j2+2.0*j3+j4)
         z(i+1)=(fo/(omega**2-alpha**2))*(cos(alpha*t(i+1))
     &          -cos(omega*t(i+1)))+sin(omega*t(i+1))
  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 USPFIT
      CALL GRSTRF

      CALL USSTTL( 'TIME', '', 'DEPARTMENT', '' )
      CALL USDAXS

      CALL UUSLNT(1)
      CALL UUSLNI(2)
      CALL UULIN(nmax+1, t, x)
      CALL USSPNT( nmax+1, t, z )
      CALL UUSLNT(2)
      CALL UUSLNI(3)
      CALL UULIN(nmax+1, t, z)
      CALL GRCLS

      stop
      end