program force
      parameter( nmax=12600 )
      real t(nmax+1), x(nmax+1), y(nmax+1)
      real u(nmax+1), v(nmax+1)
      real dt
      real k1, k2, k3, k4
      real j1, j2, j3, j4
      real l1, l2, l3, l4
      real m1, m2, m3, m4
      real tmin, tmax
      real x_0, y_0, u_0, v_0
      real f, g, l
      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(*,*) "f=, l=, g="
         read(*,*) f, l, g
      else
         f=0.015
         l=10000.0
         g=9.8
         a=g/l
      end if
*-- 時間ステップの設定 ---
      dt=0.1
*-- 初期条件の設定 ---
      t(1)=0.0
      if(IW.eq.0)then
         write(*,*) "input the initialized numbers"
         write(*,*) "x_0, y_0(init position)=, u_0, v_0(init velocity)="
         read(*,*) x_0, y_0, u_0, v_0
         x(1)=x_0
         y(1)=y_0
         u(1)=u_0
         v(1)=v_0
      else
         x(1)=0.0
         y(1)=0.0
         u(1)=0.1
         v(1)=0.1
      end if
*-- 実際の計算ここから ---
      do 10 i=1,nmax
*-- k の計算 (dx/dt=y で計算) ---
         t(i+1)=t(i)+dt
         k1=u(i)*dt
         j1=v(i)*dt
         l1=(f*v(i)-a*x(i))*dt
         m1=-(f*u(i)+a*y(i))*dt
         k2=(u(i)+0.5*l1)*dt
         j2=(v(i)+0.5*m1)*dt
         l2=(f*(v(i)+0.5*m1)-a*(x(i)+0.5*k1))*dt
         m2=-(f*(u(i)+0.5*l1)+a*(y(i)+0.5*j1))*dt
         k3=(u(i)+0.5*l2)*dt
         j3=(v(i)+0.5*m2)*dt
         l3=(f*(v(i)+0.5*m2)-a*(x(i)+0.5*k2))*dt
         m3=-(f*(u(i)+0.5*l2)+a*(y(i)+0.5*j2))*dt
         k4=(u(i)+l3)*dt
         j4=(v(i)+m3)*dt
         l4=(f*(v(i)+m3)-a*(x(i)+k3))*dt
         m4=-(f*(u(i)+l3)+a*(y(i)+j3))*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)
         u(i+1)=u(i)+(1.0/6.0)*(l1+2.0*l2+2.0*l3+l4)
         v(i+1)=v(i)+(1.0/6.0)*(m1+2.0*m2+2.0*m3+m4)
  10  continue

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

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

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

      CALL GRSWND(  -5.0,  5.0,  -5.0,  5.0 )
      CALL GRSVPT(  0.2,  0.8,  0.2,  0.8 )
      CALL USPFIT
      CALL GRSTRF

      CALL USSTTL( 'X-AXIS', '', 'Y-AXIS', '' )
      CALL USDAXS

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

      stop
      end