rem GENERAL PLANAR 3 BODY PROBLEM .... Rev 2.4
      rem A J Tooth / 17th Dec 2002 Updated 14th December 2004
      
      rem =====================================================
      rem Plots the path of 3 bodies interacting in free space.
      rem Restricted to PLANAR Motion.
      rem =====================================================
      
      on error if (err=17) then quit
      *FLOAT 64
      proc_fullscreen(xscreen%,yscreen%)  :  rem Set up use of Full Screen
      
      rem Choose colour scheme dependant on whether printing required
      print " To use a white screen, press -i- else press any key. "
      inv$=get$
      if inv$="i" or inv$="I" then
        colour 128,15 : colour 7,0 : colour 3,1 : colour 1,4
      else
        colour 128,0 : colour 7,7 : colour 3,3 : colour 1,1
      endif
      cls
      
      rem Display my Icon in compiled version
      proc_AJTicon(10,7*yscreen%/8)
      
      rem Enter Parametric Options
      proc_param
      
      rem Place Origin and draw Axes
      proc_axor
      *REFRESH OFF
      
      rem Sets initial and incremental values and ...
      rem Transfers Initial Positions and Velocities
      proc_setpar
      
      rem First Pass Indicator
      frst&=1
      
      rem Plots Initial Positions
      proc_plotbody(xAs,yAs,hor,1)
      proc_plotbody(xBs,yBs,hor,2)
      proc_plotbody(xCs,yCs,hor,3)
      *REFRESH
      
      rem MASTER PROCEDURE
      for par%=1 to lim%
        for tim%=0 to steeps%
          
          rem Computational Error Trapping
          on error local goto 1560
          
          rem Main calculations - Improved Euler Method
          r12=fn_r(xa,ya,xb,yb) : r23=fn_r(xb,yb,xc,yc) : r31=fn_r(xc,yc,xa,ya)
          rem Track Relative Distances
          proc_trac
          
          rem Quit the programme if approach gets too close.
          if (r12<0.001 or r23<0.001 or r31<0.001) then
            
            print tab(15,15);"COLLISION !!!!!"
            *REFRESH
            goto 1140
          endif
          
          sys "GetCurrentProcess" to hprocess%
          sys "SetPriorityClass", hprocess%, &100
          
          rem Track whether Total Energy remains reasonably constant
          proc_enertrac(frst&)
          frst&=0
          rem Reduce the time steps if approach is very close
          dt=fn_min(dtdef,r12,r23,r31)
          
          rem Calculate velocities resolved in x and y directions
          dua=fn_dvel(r12,r31,xa,xb,xc,dt,m2,m3)
          dva=fn_dvel(r12,r31,ya,yb,yc,dt,m2,m3)
          dub=fn_dvel(r12,r23,xb,xa,xc,dt,m1,m3)
          dvb=fn_dvel(r12,r23,yb,ya,yc,dt,m1,m3)
          duc=fn_dvel(r23,r31,xc,xb,xa,dt,m2,m1)
          dvc=fn_dvel(r23,r31,yc,yb,ya,dt,m2,m1)
          
          dxa=ua*dt : dya=va*dt : xa1=xa + dxa : ya1=ya + dya
          dxb=ub*dt : dyb=vb*dt : xb1=xb + dxb : yb1=yb + dyb
          dxc=uc*dt : dyc=vc*dt : xc1=xc + dxc : yc1=yc + dyc
          
          rem Second round calculations for Improved Euler Method
          r12=fn_r(xa1,ya1,xb1,yb1) : r23=fn_r(xb1,yb1,xc1,yc1) : r31=fn_r(xc1,yc1,xa1,ya1)
          
          dua1=fn_dvel(r12,r31,xa1,xb1,xc1,dt,m2,m3) : ua1=ua + ((dua+dua1)/2)
          dva1=fn_dvel(r12,r31,ya1,yb1,yc1,dt,m2,m3) : va1=va + ((dva+dva1)/2)
          dub1=fn_dvel(r12,r23,xb1,xa1,xc1,dt,m1,m3) : ub1=ub + ((dub+dub1)/2)
          dvb1=fn_dvel(r12,r23,yb1,ya1,yc1,dt,m1,m3) : vb1=vb + ((dvb+dvb1)/2)
          duc1=fn_dvel(r23,r31,xc1,xb1,xa1,dt,m2,m1) : uc1=uc + ((duc+duc1)/2)
          dvc1=fn_dvel(r23,r31,yc1,yb1,ya1,dt,m2,m1) : vc1=vc + ((dvc+dvc1)/2)
          
          dxa1=ua1*dt : dya1=va1*dt
          dxb1=ub1*dt : dyb1=vb1*dt
          dxc1=uc1*dt : dyc1=vc1*dt
          
          xa1=xa + ((dxa+dxa1)/2) : ya1=ya + ((dya+dya1)/2)
          xb1=xb + ((dxb+dxb1)/2) : yb1=yb + ((dyb+dyb1)/2)
          xc1=xc + ((dxc+dxc1)/2) : yc1=yc + ((dyc+dyc1)/2)
          
          rem Plot formatting section
          proc_plot(xa1,ya1,hor,1)
          proc_plot(xb1,yb1,hor,2)
          proc_plot(xc1,yc1,hor,3)
          *REFRESH
          
          rem Resets current values for x, y, u, v  at time T
          xa=xa1 : ya=ya1 : ua=ua1 : va=va1
          xb=xb1 : yb=yb1 : ub=ub1 : vb=vb1
          xc=xc1 : yc=yc1 : uc=uc1 : vc=vc1
          
          sys "GetCurrentProcess" to hprocess%
          sys "SetPriorityClass", hprocess%, &20
          
          rem Press spacebar to "crash" out the programme at any time
          if inkey(-99) then goto 1140
          
 1110   next tim%
      next par%
      
 1140 rem Repeat Options
      *REFRESH OFF
      repeat
        rep$=get$
      until rep$<>""
      cls : print tab(0,13);"Programme must now re-run or QUIT. ( r or q )"
      rep$=get$ : clg : cls
      if rep$="r" then run else quit
      
      rem END OF PROGRAMME
      end
      rem =====================================================================
      rem =====================================================================
      rem Various Functions and Procedures Section
      rem =====================================================================
      
      rem Main Function Definitions
      def fn_dvel(R2,R3,L1,L2,L3,T,M2,M3)
      local P2,P3
      
      P2=R2*R2*R2 : P3=R3*R3*R3
      vel=((M2*(L2-L1)/P2) + (M3*(L3-L1)/P3))*T
      =vel
      
      def fn_r(X1,Y1,X2,Y2)
      res=sqr(((X2-X1)*(X2-X1)) + ((Y2-Y1)*(Y2-Y1)))
      =res
      
      rem Calculates Total Energy - Kinetic and Potential Separately
      def fn_enerkin(U,V,M)
      = M*(U*U+V*V)/2
      
      def fn_enerpot(R1,R2,R3,M1,M2,M3)
      =-(M1*M2/R1)-(M2*M3/R2)-(M3*M1/R3)
      
      rem Reduces Time Steps if Approach is very close
      def fn_min(T,R1,R2,R3)
      den=(1/(R1*R1))+(1/(R2*R2))+(1/(R3*R3))
      if den>1.0 then DT=T/den else DT=T
      =DT
      rem =====================================================================
      
      rem Traps Computational Range Errors ONLY
 1560 if err=21 or err=23 or err=18 or err=20 or err=22 or err=24 then goto 1110
      else restore error
      endif
      goto 1110
      rem =====================================================================
      
      rem Plot Formatting Procedure
      def proc_plot(Xval,Yval,Hor,Col&)
      
      XP=Hor*Xval : YP=Hor*Yval
      sx%=sgn(XP) : sy%=sgn(YP)
      x%=sx%*int(abs(XP))
      y%=sy%*int(abs(YP))
      
      if abs(y%)<=1600 and abs(x%)<=2100 then
      gcol Col&
      plot 69,x%,y%
      endif
      endproc
      rem =====================================================================
      
      def proc_plotbody(Xval,Yval,Hor,Col&)
      local sx%,sy%,x%,y%,XP,YP
      
      gcol Col&
      XP=Hor*Xval : YP=Hor*Yval
      sx%=sgn(XP) : sy%=sgn(YP)
      x%=sx%*int(abs(XP))
      y%=sy%*int(abs(YP))
      
      circle fill x%,y%,10
      endproc
      rem =====================================================================
      
      def fn_dis(Num,SF%)
      
      result=(int(Num*SF%))/SF%
      =result
      rem =====================================================================
      
      def proc_trac
      
      if Cnt%=14 then
      r12dis=fn_dis(r12,10) :  r23dis=fn_dis(r23,10) : r31dis=fn_dis(r31,10)
      print tab(18,3);r12dis;
      print tab(18,4);r23dis;
      print tab(18,5);r31dis;
      endif
      endproc
      rem =====================================================================
      
      rem Track whether Total Energy remains reasonably constant
      def proc_enertrac(frst&)
      
      KinA=fn_enerkin(ua,va,m1) : KinB=fn_enerkin(ub,vb,m2) : KinC=fn_enerkin(uc,vc,m3)
      Potential=fn_enerpot(r12,r23,r31,m1,m2,m3)
      Kinetic= KinA + KinB + KinC
      E = Kinetic + Potential
      if frst&=1 then E_init=E
      del=100*(E_init-E)/E
      Cnt%+=1
      if Cnt%=15 then
      Edis=fn_dis(E,100)
      deldis=fn_dis(del,100)
      print tab(15,0);Edis
      print tab(17,7);deldis;" %     ";
      print tab(15,1);int(Kinetic);"   "; : print tab(15,2);int(Potential);"   ";
      Cnt%=0
      endif
      endproc
      rem =====================================================================
      
      rem Places axes, origin and preamble
      def proc_axor
      
      cls
      xmax%=xscreen% : ymax%=yscreen%
      
      origin xmax%,ymax% : gcol 7
      move -xmax%,0 : draw xmax%,0
      move 0,-ymax% : draw 0,ymax%
      
      print tab(0,0);"   Energy ="
      
      print tab(0,1);"  Kinetic ="
      print tab(0,2);"Potential ="
      print tab(0,3);"Distance A to B ="
      print tab(0,4);"Distance B to C ="
      print tab(0,5);"Distance C to A ="
      print tab(0,7);"Energy Error =";
      proc_AJTicon(10,7*yscreen%/8)
      endproc
      rem =====================================================================
      
      rem Parameter Entry
      def proc_param
      
      print " Just press -Enter- to select default parameter values!"
      print
      input " How many steps for the Time parameter T (2000)",steeps%
      if steeps%=0 then steeps%=2000
      rem Initial Positions
      input " What initial x,y coords for Body A (-3,0)",xAs$,tab(42,4),yAs$
      if xAs$="" then xAs=-3 else xAs=eval(xAs$)
      if yAs$="" then yAs=0 else yAs=eval(yAs$)
      input " What initial x,y coords for Body B (3,0) ",xBs$,tab(42,6),yBs$
      if xBs$="" then xBs=3 else xBs=eval(xBs$)
      if yBs$="" then yBs=0 else yBs=eval(yBs$)
      input " What initial x,y coords for Body C (0,3) ",xCs$,tab(42,8),yCs$
      if xCs$="" then xCs=0 else xCs=eval(xCs$)
      if yCs$="" then yCs=3 else yCs=eval(yCs$)
      rem Initial Velocities
      input " What intital u,v velocities for Body A (0,4) ",uAs$,tab(46,10),vAs$
      if uAs$="" then uAs=0 else uAs=eval(uAs$)
      if vAs$="" then vAs=4 else vAs=eval(vAs$)
      input " What intital u,v velocities for Body B (0,-4)",uBs$,tab(46,12),vBs$
      if uBs$="" then uBs=0 else uBs=eval(uBs$)
      if vBs$="" then vBs=-4 else vBs=eval(vBs$)
      input " What intital u,v velocities for Body C (0,0) ",uCs$,tab(46,14),vCs$
      if uCs$="" then uCs=0 else uCs=eval(uCs$)
      if vCs$="" then vCs=0 else vCs=eval(vCs$)
      
      input " Enter Mass of Body A (200)",m1
      if m1=0 then m1=200
      input " Enter Mass of Body B (210)",m2
      if m2=0 then m2=210
      input " Enter Mass of Body C (170)",m3
      if m3=0 then m3=170
      
      input " How many years (100)",lim%
      if lim%=0 then lim%=100
      input " Axis Scale  (250)",hor
      if hor=0 then hor=250
      ver=hor
      Cnt%=0
      endproc
      rem =====================================================================
      
      rem Set parameters
      def proc_setpar
      
      rem Sets initial and incremental values
      dtdef=1/steeps%
      
      rem Transfer Initial Positions and Velocities
      xa=xAs : ya=yAs
      xb=xBs : yb=yBs
      xc=xCs : yc=yCs
      
      ua=uAs : va=vAs
      ub=uBs : vb=vBs
      uc=uCs : vc=vCs
      endproc
      rem =====================================================================
      
      rem Set up use of Full Screen
      def proc_fullscreen(return xscreen%,return yscreen%)
      
      sys "GetSystemMetrics", 0 to xscreen%
      sys "GetSystemMetrics", 1 to yscreen%
      sys "SetWindowLong",@hwnd%,-16,&16000000
      sys "SetWindowPos",@hwnd%,-1,0,0,xscreen%,yscreen%,0
      vdu 23,22,xscreen%;yscreen%;8,16,16,1 : rem Set fullscreen mode
      mouse off : off : rem Turns off the Mouse Pointer and the Cursor
      endproc
      rem =====================================================================
      
      rem Displays my Icon in .exe version
      def proc_AJTicon(i%,j%)
      
      sys "GetModuleHandle", 0 to hm%
      sys "LoadImage", hm%, "BBCWrun", 1, 32, 32, 0 to hicon%
      w% = 32
      h% = 32
      sys "DrawIconEx", @memhdc%, i%, j%, hicon%, w%, h%, 0, 0, 3
      sys "InvalidateRect", @hwnd%, 0, 0
      endproc
      rem =====================================================================