rem SPRING PENDULUM Simulation .... Rev 1.1
      rem A J Tooth / 28th November 2003
      
      rem +++++++++++++++++++++++++++++++++++++++++++++++++
      rem This program models the non-linear motion of
      rem a pendulum attached to a spring.
      rem +++++++++++++++++++++++++++++++++++++++++++++++++
      
      on error if (err=17) then quit
      proc_fullscreen(22)  :  rem Set up use of Full Screen
      *FONT Arial,10,B
   60 rem Choose colour scheme dependant on whether printing required
      cls : print "Press -i- to use a WHITE background"
      inv$=get$
      if inv$="i" or inv$="I" then
        colour 128,15 : colour 7,0 : colour 3,1 : colour 6,0
      else
        colour 128,0 : colour 7,7 : colour 3,3 : colour 6,6
      endif
      clg
      
      rem Format of Main Differential Functions
      pdot$="-D*G*(SIN(P))/(E+X)"
      
      xdot$="(((E+X)*Q*Q)+(G*COS(P))-((K*(C+X)/M)))*D"
      
      rem Enter Parametric Options
      print " Suggested trial values appear in brackets. Just press -Enter- to use the defaults."
      print:input " How many steps for the Time parameter T (15000)",stps%
      if stps%=0 then stps%=15000
      print:input " What value for initial angle (degrees) from vertical (45)",P0S
      if P0S=0 then P0S=45
      print:input " What initial value for angular velocity of pendulum, if any (0)",Q0S
      print:input " What value for the mass of the bob (0.1)",m
      if m=0 then m=0.1
      print:print " Unstretched rod length is 1m. "
      input " How much further would a 1kg weight stretch the spring hanging still. (0.2)",l
      if l=0 then l=0.2
      print:input " What is the initial extension (0)",X0
      print:input " What is the initial radial speed (0) ",Xdot0
      print:input " How many seconds (30) ",lim%
      if lim%=0 then lim%=30
      
      scal=600 : g=9.81
      L=1+l : lam=g/l : rem Calculate spring modulus
      clg : cls
      xmax%=1024 : ymax%=768
      rem Place Origin and draw Axes
      setx%=xmax% : sety%=(3/2)*ymax%
      origin setx%,sety% : gcol 3
      move -xmax%,0 : draw xmax%,0
      move 0,0 : draw 0,ymax%/2
      *FONT Castellar,18
      colour 2
      print tab(5,5);"The Spring - Pendulum"
      rem Curve Drawing Section
      gcol 6
      rem Sets initial and incremental values
      dtim=1/stps%
      P0=rad(P0S) : Q0=Q0S
      x0=X0 : xdot0=Xdot0
      rem Plots Initial Value
      proc_plot(P0,scal,(L+x0))
      
      for par%= 0 to lim%
        for tim%=0 to stps%
          
          sys "GetCurrentProcess" to hprocess%
          sys "SetPriorityClass", hprocess%, &100
          rem Main calculations - Improved Euler Method
          
          dQ=fn_dfuncp(pdot$,P0,dtim,x0,L,g)
          Q1=Q0 + dQ
          dP=Q0*dtim : P1=P0 + dP
          
          dxdot=fn_dfuncx(xdot$,Q0,P0,dtim,x0,lam,l,L,g,m)
          xdot1=xdot0 + dxdot
          dx=xdot0*dtim : x1=x0 + dx
          
          dQ2=fn_dfuncp(pdot$,P1,dtim,x1,L,g)
          dP2=Q1*dtim
          Q1=Q0 + ((dQ+dQ2)/2) : P1=P0 + ((dP+dP2)/2)
          
          dxdot2=fn_dfuncx(xdot$,Q1,P1,dtim,x1,lam,l,L,g,m)
          dx2=xdot1*dtim
          xdot1=xdot0 + ((dxdot+dxdot2)/2) : x1=x0 + ((dx+dx2)/2)
          
          rem Plot formatting section
          proc_plot(P1,scal,(L+x1))
          sys "GetCurrentProcess" to hprocess%
          sys "SetPriorityClass", hprocess%, &20
          
          rem Resets current values for P,Q,x and xdot at time T
          Q0=Q1 : P0=P1 : x0=x1 : xdot0=xdot1
  750   next tim%
      next par%
      
      repeat until inkey(0)=-1
      rem Repeat Options
      rep$=get$
      clg
      print tab(10,10);" Re-run, or QUIT the programme? ( r or q )"
      rep$=get$ : clg
      if rep$="r" or rep$="R" then flag%=1 else flag%=2
      on flag% goto 60,860
  860 quit
      :
      rem END OF PROGRAMME
      end
      rem =====================================================================
      
      rem Various Functions and Procedures Section
      
      rem Main Calculation for x and y
      def fn_dfuncp(main1$,P,D,X,E,G)
      =eval(main1$)
      
      def fn_dfuncx(main2$,Q,P,D,X,K,C,E,G,M)
      =eval(main2$)
      
      rem Plot Formatting Procedure
      def proc_plot(Pval,Scal,Lgth)
      local XP,YP,sx%,sy%,x%,y%,xl%,yl%,Rod
      Rod=Scal*Lgth
      XP=-Rod*sin(Pval) : YP=-Rod*cos(Pval)
      sx%=sgn(XP) : sy%=sgn(YP)
      x%=sx%*int(abs(XP))
      y%=sy%*int(abs(YP))
      
      if abs(y%)>1600 then goto 1240
      if abs(x%)>2100 then goto 1240
      rem Prints pendulum shaft, and locus of end-point
      xl%= int(0.6*x%) : yl%= int(0.6*y%)
      move 0,0 : gcol 1 : plot 5,xl%,yl% : gcol 6 : move x%,y% : draw x%,y%
      sys "UpdateWindow",@hwnd%
      rem Un-prints pendulum shaft, to simulate movement
      move 0,0 : gcol 0 : plot 5,xl%,yl%
 1240 endproc
      
      rem Set up use of Full Screen
      
      def proc_fullscreen(N%)
      sys "GetSystemMetrics", 0 to xscreen%
      sys "GetSystemMetrics", 1 to yscreen%
      sys "SetWindowLong",@hwnd%,-16,&16000000
      sys "SetWindowPos",@hwnd%,-1,0,0,xscreen%,yscreen%,0
      mode N%
      mouse off : off : rem Turns off the Mouse Pointer and the Cursor
      endproc
      
      def proc_enertrac
      rem Energy Tracking
      T=m*((xdot0*xdot0) + ((L+x0)*(L+x0)*Q0*Q0))/2
      V=-m*g*((L+x0)*cos(P0)) + lam*(x0*l + (x0*x0/2))
      Energy=T+V : print tab(5,5);"ENERGY = ";Energy
      endproc