rem N Body Problem ...... Rev 4.4
      rem A J Tooth // May 2005
      
      rem ========================================================
      rem Global variables start with a CAPITAL.
      rem Local variables are lower-case.
      rem ========================================================
      rem Simulates N asteroids interacting under gravity, near to
      rem a massive planet, which also moves under the net
      rem gravitational force of the asteroids.
      rem ========================================================
      rem on error if (err=17) then quit
      
      rem Setup
      proc_setup
      
      rem Initial Net View
      proc_astroview
      
      Trk%=0
      colour 7 : print tab(1,0);"Energy Error"
      print tab(1,5);"Dt= ";
      
      repeat
        rem Calculate Resultant Force from viewpoint of each asteroid in turn
        proc_Force
        
        rem Main Iterative Movement Calculation
        proc_calc
        
        rem Energy Tracking
        proc_energy
        
        rem Next view
        proc_astroview
        
        Trk%+=1 : if Trk%=100 then Trk%=0
        a$=inkey$(5)
        mouse x,y,b&
      until (a$<>"" or b&=4 or Efl%=1)
      
      a$=inkey$(20)
      repeat
        a$=inkey$(1)
        mouse x,y,b&
      until (a$<>"" or b&=4)
      quit
      
      end
      rem End of Program ===================================================
      
      rem Setup fullscreen
      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
      off : mouse off
      endproc
      rem ==================================================================
      
      rem Setup
      def proc_setup
      local a%,b&,x,y,tst%,i%
      
      himem=lomem + 20000000
      *FLOAT 64
      dim rgb% 3, re% 0, gr% 0, bl% 0
      dim pt$(3)
      pt$()="short","medium","long"
      
      proc_fullscreen(xscreen%,yscreen%)
      
      *REFRESH OFF
      Num%=50
      proc_back("How many Asteroids?",100,1250,500,100,200,0,200)
      proc_box(4,700,1250,200,str$(Num%),50,50,0)
      repeat
        mouse x,y,b&
        tst&=asc(inkey$(1))
        if tst&=140 or tst&=138 then Num%+=1
        if tst&=141 or tst&=139 then Num%-=1
        if Num%<1 then Num%=1
        if Num%>1000 then Num%=1000
        if tst&=140 or tst&=141 or tst&=138 or tst&=139 then proc_box(4,700,1250,200,str$(Num%),50,50,0)
        tst%=0
        *REFRESH
      until b&=4
      a$=inkey$(20)
      
      proc_back("Choose a default time interval",100,1000,700,100,200,0,200)
      i%=0
      proc_box(4,900,1000,250,pt$(i%),50,50,0)
      repeat
        mouse x,y,b&
        tst&=asc(inkey$(1))
        if tst&=140 or tst&=138 then i%+=1
        if tst&=141 or tst&=139 then i%-=1
        if i%<0 then i%=2
        if i%>2 then i%=0
        if tst&=140 or tst&=141 or tst&=138 or tst&=139 then proc_box(4,900,1000,250,pt$(i%),50,50,0)
        tst%=0
        *REFRESH
      until b&=4
      
      t$=left$(pt$(i%),1)
      case t$ of
        when "s": DT=0.001
        when "m": DT=0.005
        when "l": DT=0.01
      endcase
      G=10.0 : Flg%=0 : Efl%=0
      cls
      *REFRESH OFF
      
      dim Ver(Num%,2), Vel(Num%,4), Acc(Num%,2), Ms(Num%), dVer(Num%,4,2)
      dim Dis(Num%,Num%)
      
      Xlim%=xscreen% : Ylim%=yscreen%
      origin Xlim%,Ylim%
      
      rem Masses, and initial positions and velocities
      Ms(1)=1000000.0
      Ver(1,1)=0.0
      Ver(2,1)=0.0
      Vel(1,1)=0.0
      Vel(2,1)=0.0
      Vero1=0.0
      Vero2=0.0
      for a%=2 to Num%
        Ms(a%)=1000.0*rnd(5)
        
        Ver(a%,1)=1.0*(rnd(Xlim%) - Xlim%/2)
        Ver(a%,2)=1.0*(rnd(Ylim%) - Ylim%/2)
        
        Vel(a%,1)=100.0*(rnd(11)-6)
        Vel(a%,2)=100.0*(rnd(11)-6)
      next a%
      
      endproc
      rem ==================================================================
      
      rem Initial Net View
      def proc_astroview
      local a%
      
      rem Draw the vertices themselves
      for a%=1 to Num%
        if a%=1 then gcol 0 : circle fill 2*Vero1,2*Vero2,20
        if Acc(a%,0)<0.5 then
          !rgb%=100000*a%
          ?re%=?rgb%
          ?gr%=rgb%?(1)
          ?bl%=rgb%?(2)
          colour 10,?re%,?gr%,?bl% : gcol 10
          plot 2*Ver(a%,1),2*Ver(a%,2)
        endif
        if a%=1 then circle fill 2*Ver(a%,1),2*Ver(a%,2),20 : Vero1=Ver(a%,1) : Vero2=Ver(a%,2)
      next a%
      *REFRESH
      endproc
      rem ==================================================================
      
      rem Calculate Resultant Force experienced by each asteroid in turn
      def proc_Force
      local i%,j%,cumFx,cumFy,Rx,Ry,D,Cub,R,Tadj
      
      Coll%=0 : rem Set to 1 if an asteroid collides with the planet.
      Dt=DT   : rem Default maximum value for time slot.
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &100
      
      for j%=1 to Num%
        cumFx=0.0 : cumFy=0.0 : Tadj=1.0
        if Acc(j%,0)<0.5 then
          for i%=1 to Num%
            if (i%<>j% and Acc(i%,0)<0.5) then
              Rx=Ver(i%,1) - Ver(j%,1)
              Ry=Ver(i%,2) - Ver(j%,2)
              D=sqr(Rx*Rx + Ry*Ry)
              if D<0.01 then D=0.01
              Dis(j%,i%)=D
              if (j%=1 and D<10.0) then Acc(i%,0)=1.0 : Coll%=1
              rem Increase Tadj whenever the subject asteroid
              rem gets too close to any of the object asteroids
              if Tadj<1000/D then Tadj=1000/D
              Cub=D*D*D
              cumFx+=Ms(i%)*Rx/Cub
              cumFy+=Ms(i%)*Ry/Cub
            endif
          next i%
        endif
        rem Dt is set at the minimum value reached
        if DT/Tadj<Dt then Dt=DT/Tadj
        
        cumFx=G*cumFx : cumFy=G*cumFy
        Acc(j%,1)=cumFx : Acc(j%,2)=cumFy
      next j%
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &20
      endproc
      rem ==================================================================
      
      rem Main Iterative Movement Calculation
      def proc_calc
      local j%,k%
      
      if (Trk% mod 20 =0) then colour 3 : print tab(5,5);Dt;"            "
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &100
      
      rem Add the mass of each new asteroid which hits the planet
      rem to the mass of the planet, but only once.
      for k%=2 to Num%
        if (Acc(k%,0)>0.5 and Acc(k%,0)<1.5) then Ms(1)+=Ms(k%) : Acc(k%,0)=2.0
      next k%
      
      rem Solve the equations of motion for each asteroid
      rem AND the planet.
      for j%=1 to Num%
        if Acc(j%,0)<0.5 then
          rem Velocity
          Vel(j%,3)=Vel(j%,1) +  Dt*Acc(j%,1)
          Vel(j%,4)=Vel(j%,2) +  Dt*Acc(j%,2)
          
          rem Adjustements
          dVer(j%,1,1)=Dt*Vel(j%,1) : dVer(j%,1,2)=Dt*Vel(j%,2) : rem C1
          dVer(j%,2,1)=(Dt/2)*(Vel(j%,1) + Vel(j%,3))/2         : rem C2
          dVer(j%,2,2)=(Dt/2)*(Vel(j%,2) + Vel(j%,4))/2
          
          rem Reassign values
          Ver(j%,1)=Ver(j%,1) + (dVer(j%,1,1) + 2*dVer(j%,2,1))/2
          Ver(j%,2)=Ver(j%,2) + (dVer(j%,1,2) + 2*dVer(j%,2,2))/2
          Vel(j%,1)=Vel(j%,3) : Vel(j%,2)=Vel(j%,4)
        endif
      next j%
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &20
      
      endproc
      rem ==================================================================
      
      rem Energy Tracking
      def proc_energy
      local j%,k%,vx,vy,vsq,PK,PE,PEj
      
      rem Potential Energy
      PE=0.0
      for j%=1 to Num%
        PEj=0.0
        for k%=1 to Num%
          if (Acc(k%,0)<0.5 and j%<>k%) then PEj+=Ms(k%)/Dis(j%,k%)
        next k%
        if Acc(j%,0)<0.5 then PEj=PEj*Ms(j%)*G : PE+=PEj
      next j%
      rem Potential energy value must be halved, since it has been exactly double-counted.
      rem It also has to be negated.
      PE=-PE/2
      
      rem Kinetic Energy
      PK=0.0
      for j%=1 to Num%
        if Acc(j%,0)<0.5 then
          vx=Vel(j%,1) : vy=Vel(j%,2)
          vsq=vx*vx + vy*vy
          PK+=Ms(j%)*vsq/2
        endif
      next j%
      
      rem Combine, compare and report
      Energy=PK + PE
      if Trk% mod 50 =0 then print tab(10,1);"       "
      if Coll%=1 then Flg%=0 : colour 1 : print tab(10,1);"RESET!"
      if Flg%=0 then InitEnergy=Energy : Flg%=1
      NetEnergy=Energy - InitEnergy
      colour 7 : print tab(1,0);"Energy Error"
      Dif=(int(1000*NetEnergy/InitEnergy))/10
      colour 3: print tab(1,1);Dif;" %   ";
      if abs(Dif)<=10.0 then colour 2: print tab(1,2);"within limits     "
      if abs(Dif)>10.0 and abs(Dif)<=30.0 then colour 5 : print tab(1,2);"losing accuracy"
      if abs(Dif)>30.0 then Efl%=1 : colour 1 : print tab(1,2);"ACCURACY LOST : PROGRAM STOPPED !!"
      *REFRESH
      endproc
      rem ==================================================================
      
      rem Prints a coloured button to the screen
      def proc_box(Mul%,Xpos%,Ypos%,Ws%,Msg$,re%,gr%,bl%)
      local h%
      vdu 5
      for h%=0 to 30
        colour 8,Mul%*h%*re%/30,Mul%*h%*gr%/30,Mul%*h%*bl%/30
        gcol 8 : rectangle fill (Xpos%+2*h%),(Ypos%+2*h%),(Ws%-4*h%),(150-4*h%)
      next h%
      gcol 0 : move (Xpos%+60),(Ypos%+90) : print;Msg$;
      vdu 4
      endproc
      rem ========================================================================
      
      rem Prints a text-message backdrop
      def proc_back(Msg$,Xs%,Ys%,Ws%,Hs%,Rr%,Gg%,Bb%)
      local h%,Rf%,Gf%,Bf%
      
      vdu 5
      for h%=0 to 30
        Rf%=Rr%*h%/30 : Gf%=Gg%*h%/30 : Bf%=Bb%*h%/30
        colour 9,Rf%,Gf%,Bf% : gcol 9
        rectangle fill Xs%+h%,Ys%+h%,Ws%-2*h%,Hs%-2*h%
      next h%
      gcol 0 : move Xs%+50,Ys%+65 : print;Msg$;
      vdu 4
      endproc
      rem ========================================================================