How to build n-body gravity algorithm with source code

This is essentially the solution to 3D-n-body-gravity within the theoretical framework of Newtonian-Planck gravity. This computer algorithm demonstrates the gravitational effects that the various planetary bodies have on one another. The code is from Visual Basic 6, but should be easily translated into any other language. The most important lines of code are bold-yellow, other code is in white, remarks in teal.
     
 

'GRate is a time-scale variable that is a square of TimeRate!
100
GRate = TimeRate^2
'set both rates to = 1 for simplest solution.

For Ss3 = 0 to Planets 'zero is the sun's number.
TraVX = 0: TraVY = 0: TraVZ = 0 'planet velocity Traversed XYZ
For Ss2 = 0 to Planets
If Ss3 = Ss2 Then GoTo 200 'no gravity on planet itself.
Ss2 is the effective planet, Ss3 is is the effected planet.

'PoX PoY PoZ are real positions in km for XYZ axes.
'Preset the positions yourself prior to this routine.
Xx = PoX(Ss2) - PoX(Ss3)'Get the distance between the bodies.
Yy = PoY(Ss2) - PoY(Ss3)
Zz = PoZ(Ss2) - PoZ(Ss3)
'3d pythagorus, Dist = distance between bodies.
Dist = ((Xx * Xx) + (Yy * Yy) + (Zz * Zz)) ^ (0.5)

'Radbody = radius of body. Exclude close encounters here:
If Dist > RadBody Then
  '3d Newton's equation * time-scale variable:
  Gfor = ((Mass(Ss2)*Gee) / (Dist^2)) * GRate
  'GRate = time-scale.
  'Gfor = virtual gravity force.

  Zp = Zz / Dist 'Zp is Z-proportion from Zz distance, ditto Xx & Yy.
  Yp = Yy / Dist
  Xp = Xx / Dist
  GravX = Gfor * Xp 'Gravity distributed proportionally.
  GravY = Gfor * Yp
  GravZ = Gfor * Zp
  'GraVX GraVY GraVZ are individual shifts in velocity caused by
  'gravity-force for each isolated interaction in quantum time.
End If

TraVX = TraVX + GraVX
TraVY = TraVY + GraVY
TraVZ = TraVZ + GraVZ
'TraVX TraVY TraVZ are accumulated velocities traversing.
'for each planet in quantum time.

200
Next Ss2

'TotalX TotalY TotalZ are combined shifts in velocity.
'for a single step in time from all planets.
TotalX(Ss3) = TraVX
TotalY(Ss3) = TraVY
TotalZ(Ss3) = TraVZ
Next Ss3

'The first pair of for...next loops are now complete.
'They have determined what the changes to the planetary velocities will be.

'Now this is where the velocities are actually changed for each time unit.
'MoX MoY MoZ are motion variables for all time, so do not set them to zero.
'Starting motion/velocities not here included, define those yourself prior to this routine.
For Ss4 = 0 to Planets
MoX(Ss4) = MoX(Ss4) + TotalX(Ss4)
MoY(Ss4) = MoY(Ss4) + TotalY(Ss4)
MoZ(Ss4) = MoZ(Ss4) + TotalZ(Ss4)
Next Ss4

'The velocity has changed, now we move the positions of each planet.
For Ss5 = 0 to Planets 'positions moved by motion variables.
PoX(Ss5) = PoX(Ss5) + MoX(Ss5)
PoY(Ss5) = PoY(Ss5) + MoY(Ss5)
PoZ(Ss5) = PoZ(Ss5) + MoZ(Ss5)
Next Ss5

'now put the actual dots on the screen:
For Ss6 = 0 To Planets
'Zoom is distance scale,
'balX balY balZ balance position on screen. Figure these yourself.
PsX = (PoX(Ss6) / Zoom) + balX
PsY = (PoY(Ss6) / Zoom) + balY
'PsX PsY PsZ are screen positions in pixels.
balX, balY, balZ are your screen position variables balanced to center-screen.
PSet (PsX, PsY), PlanetColor 'draw pixel: top-view, main screen.

'zxaxisZ will also require your zoom and balance adjustments
zxaxisZ = (PoZ(Ss6) / Zoom) + balZ
'zxaxisY and zxaxisX are similar to PsX & PsY
ViewXZ.PSet (zxaxisX, zxaxisZ), PlanetColor 'draw side-view (XZ-axis).
ViewZY.PSet (zyaxisZ, zyaxisY), PlanetColor 'draw side-view (ZY-axis).

Next Ss6

Goto 100 'not required if using an event-timer that loops itself.

  
     

By breaking it down like this, the computation yields the same results regardless of which planet takes which number in the sequence. Its just a matter of not trying to do all the steps at once. That is why we complete the loops for determining the potential motions for all the planetary interactions (Ss2 & Ss3), before we actually move any of them (Ss5). Using a single large formula can also often increase the rounding errors, which is why accurate coding uses numerous simple steps.

Importantly, GRate is the most interesting variable.
Consider the first line of code in bold: GRate = TimeRate^2

So if we triple the TimeRate, then we must multiply the computational force of gravity (GFor) by 9. It should be realized that this is the reverse consequence of the inverse-of-the-square law of Newton. Using the GRate variable yields an exquisitely close approximation which gets more accurate as the size of the jumps in time get smaller. But it also allows one to speed up the process in order to get good averages for many orbits. When changing between time-rates in mid-scenario it is still ok, but not entirely accurate, depending on the amount of change.

Using quantum jumps of 2.5 hours in the algorithm OGS12 (orbit-gravity-sim-12.exe) an error-margin for the distance from Mercury to the Sun was attained that totaled just 177mm over 20 orbits. The orbit did fluctuate +-8m between these orbits. Smaller quantum jumps in time than 15000 seconds improves this error-margin.

Perihelion Precession requires a finer time quantum. 15 virtual seconds per calculation for Mercury was an optimal balance between accuracy and the speed of the computational process itself. 150 virtual seconds per iteration gives a good average over a 237-year sample.

Of course Mass(Ss2) is the mass of each planet whereby the variable Ss2 is the ordinal number for each planet. If you are not used to computer code, it would be easy to be confused by two different types of brackets. (Ss2) is not a mathematical use, instead it designates the planet's number, where normally Mass(0) is the mass of the Sun, and Mass(1) holds the mass of Mercury, etc.
     
 

The txt file following contains the Horizon Ephemeris positional data for the major planets as well as my re-evaluation of the planetary velocities as they should be as a result of the given orbital duration averages.

exact-planet-velocity-position.txt
(11 mb)

The data in this file is for 3 separate temporal points in the years 1773, 1900 and 1940:

Jan/02/1900 @ 00:08 - Earth-Moon at Perihelion
Jan/24/1940 @ 00:00 - Jupiter quite near Perihelion
Dec/18/1773 @ 12h53 - Jupiter at Perihelion

Your questions can be asked at this forum: cosmology.africamotion.net

See the section on Uranus for proof; that has been repeated within these computations, like thousands of times. Horizon Ephemeris velocity data has a 7%+ error.

Also see: Sorting Horizons for more details on that 7%+ error.

  
     

To calculate the Perihelion Precession, the angular shift in arc-seconds requires a polygon-orbit of about 13 million calculations for an error margin of 0.1 arc-seconds per orbit. This is because there are 1.296 million arc-seconds in any orbit.

The rate of such computation will vary considerably for different computers when using more efficient algorithms like OGS13 and higher. These algorithms can operate with turbo-time, which does not visibly demonstrate the change in velocity and distances for each step in the orbit. Turbo-time is faster to compute, but less-detailed, giving orbital statistics only after each orbit is complete.

So you have to wait for the computer to crunch those 13 million or more calculations per orbit. My bucket-of-bolts from 2014 cranked along at about 7 million per hour for 9-body gravity. Of course it would be ideal if the quanta were at Planck-time, but computers are far too slow for a perfect polygon-orbit which would have about 7 x 10^50 sides to it.

There is quite a great deal more coding which measures the various details of the orbits, like orbit-count, aphelion, and duration, etc. That has been here omitted for the sake of simplicity in explanation, so obviously you'll need to add these in yourselves. You'll also need to add variable declarations and place the visual objects like the 'pictureboxes', which I have simply called 'viewXZ' (last lines of code). This shows the side-view with the Z-axis.

You need to also note the starting positions and velocities for the planets. To prevent the entire solar-system wandering off the side of the screen it has to be realized that the Sun has a sum of counter-balancing velocities proportional to the masses of the planets.

For Jupiter, the Sun's velocity is about 1000th the velocity of Jupiter in precisely the opposite direction to Jupiter. This is because the Sun is 1000 times the mass of Jupiter. Momentums between the two must be perfectly equal and opposite to each other.

I have been using the Horizon Ephemeris for starting positions, but it is strongly advised to read the analysis of that process on the next page before starting to place planets in your algorithm. While the OGS algorithms are precise algebraically, finding the exact positions is not quite plain sailing. Often those parameters are based on computations that use statistical Keplerian processes, or less reliable 2D formulae; and combinations of both of these.

This is the crux of the matter. The given velocities fromHorizon Ephemeris certainly have an error of more than 7%. You will soon see this if you try and apply their given velocities with the N-body utility. This is why I have provided all my conclusions of the improved planetary velocities in the file: exact-planet-velocity-position.txt (see details earlier on this page). These were painstakingly reverse engineered from their given positions and durations of their orbits. NASA take note: Your bloody velocities are wrong! Dammit! Corroborated independently by Bernard Burchell too!

     
   
     

 

Sections of this Article by web-page



Visit homepage :

n-body gravity from www.flight-light-and-spin.com

 

Visit forum:

Cosmology Forum