rem Rocket II ...... Rev 2.0
      rem A J Tooth // May 2007
      
      rem ===================================================================
      on error if (err=17) then quit
      *FLOAT 64
      install @lib$+"MyUtils.bbc"
      rem ===================================================================
      
      rem Initial Setup
      proc_setup(Xs%,Ys%)
      
      rem Draw the axes
      proc_axes(Xs%,Ys%)
      
      rem Parameters
      proc_param
      
      for b=1.0 to 5.0 step 0.5
        rem Launch the Rocket
        proc_Launch(b)
      next b
      
      proc_event(a$,b&)
      if (b&=4 or a$=" ") then run else quit
      quit
      rem End of Program =================================================
      rem ================================================================
      
      rem Initial Setup
      def proc_setup(return xs%,return ys%)
      
      proc_fullscreen(xscreen%,yscreen%)
      
      xs%=80 : ys%=500 : rem Origin placement
      
      origin xs%,ys%
      
      endproc
      rem ================================================================
      
      rem Draw axes
      def proc_axes(Xs%,Ys%)
      
      gcol 1
      move -Xs%,0 : draw 2*xscreen%-Xs%,0
      move 0,-Ys% : draw 0,2*yscreen%-Ys%
      endproc
      rem ================================================================
      
      rem Parameters
      def proc_param
      
      rem Universal variables
      X0=0.0     : rem Initial x-position
      Y0=0.0     : rem Initial y-position
      g=9.81     : rem Acceleration due to gravity [m/s/s]
      Lam=0.01   : rem Strength of drag
      M=0.5      : rem Mass of Rocket [kg]
      m0=3.0     : rem Initial Mass of Fuel [kg]
      b=1.0      : rem Burn-rate of fuel [kg/s]
      u=120.0    : rem Exhaust velocity of gases [m/s]
      dt=0.005   : rem Time step [s]
      Tlm=15.0    : rem Time limit [s]
      Xd=250.0    : rem X-scale [m]
      Yd=250.0    : rem Y-scale [m]
      
      endproc
      rem ================================================================
      
      rem Launch the Rocket
      def proc_Launch(b)
      
      x0=X0 : y0=Y0
      xdot0=0.01 : ydot0=0.09
      Tim=0.0
      
      repeat
        
        com=1.0
        proc_common(xdot0,ydot0,Tim,Mtr&,com)
        
        cx1=dt*xdot0*com
        cy1=dt*(ydot0*com - g)
        
        xdot11=xdot0+cx1 : ydot11=ydot0+cy1
        
        Tim+=dt
        proc_common(xdot11,ydot11,Tim,Mtr&,com)
        
        cx2=dt*xdot11*com
        cy2=dt*(ydot11*com - g)
        
        xdot1=xdot0 + (cx1+cx2)/2
        ydot1=ydot0 + (cy1+cy2)/2
        
        x1=x0 + (xdot0+xdot1)*dt/2
        y1=y0 + (ydot0+ydot1)*dt/2
        
        rem Plot next point
        proc_pltt(Mtr&,x1,y1,Xd,Yd,Xs%,Ys%)
        
        x0=x1 : y0=y1 : xdot0=xdot1 : ydot0=ydot1
      until Tim>Tlm
      endproc
      rem ================================================================
      
      rem Function for common factor
      def proc_common(x,y,T,return C&,return res)
      
      if m0>b*T then
        m=m0-b*T : C&=1
      else
        m=0.0 : C&=0
      endif
      res=(b*u*C& - Lam*(x*x + y*y))/((M+m)*sqr(x*x + y*y))
      endproc
      rem ================================================================
      
      rem Plot a point
      def proc_pltt(Mtr&,x,y,Xm,Ys,Sx%,Sy%)
      local Xp%,Yp%
      
      Xp%=int(x*(2*xscreen%-Sx%)/Xm)
      Yp%=int(y*(2*yscreen%-Sy%)/Ys)
      
      case Mtr& of
        when 1: gcol 3
        when 0: gcol 4
      endcase
      plot Xp%,Yp%
      
      endproc
      rem ================================================================