3D N-Body Orbit Bench

What is modelled

Point masses attracting each other with Newton's law of gravitation in three dimensions. The units are dimensionless "code units" with G = 1 by default, so only ratios matter: a star of mass 1000 with a planet at distance 5 has a period of 2 pi sqrt(r^3 / (G (M + m))) = 2.2 time units. Nothing here is scaled to a real solar system.

Equations

• Acceleration: a_i = G sum_j m_j (r_j - r_i) / (|r_j - r_i|^2 + eps^2)^(3/2); eps is a small optional softening length that stops two bodies from producing infinite forces. Set it to 0 for pure point masses.
• Integrator: kick-drift-kick leapfrog (velocity Verlet). It is symplectic, so energy does not drift away secularly - it wobbles within a small band. The time step shrinks automatically during close encounters.
• Total energy E = sum m v^2/2 - G sum_(i<j) m_i m_j / r_ij and angular momentum L = sum m r x v are conserved by the true dynamics; graph 2 shows how well the code keeps them.
• Selected body: osculating elements about the heaviest body with mu = G (M + m): orbital energy v^2/2 - mu/r, a = -mu / (2 energy), e = sqrt(1 + 2 energy h^2 / mu^2), i = acos(h_z / |h|) with h = r x v, and the Kepler period T = 2 pi sqrt(a^3 / mu).
• Figure-eight: three equal masses on one planar curve (Chenciner and Montgomery, 2000; initial data from Simo), period 6.3259 in these units, momentum zero.

Checks run in the page

• A circular orbit (r = 5, M = 1000) is integrated until the angle has advanced by 2 pi; the measured period is compared with Kepler's third law (within 0.1 %).
• The figure-eight is integrated for exactly one period: bodies must return to their starting positions, energy must stay within 1e-5 and the z-component of total angular momentum must stay 0.
• A 28 degree tilted ellipse is integrated for one Kepler period; the osculating a, e and inclination computed from the state vectors must stay constant (a two-body orbit is exactly Keplerian).

Honest simplifications

• Body sizes are exaggerated (slider) so that planets are visible; collisions use the true size of the point-mass model, so a merge can look like a near miss on screen. Mergers are perfectly inelastic: momentum is conserved, kinetic energy is lost, so the energy reference resets.
• Only Newtonian point-mass gravity: no tides, no relativity, no radiation pressure, no non-spherical shapes, no spin-orbit coupling. Planet spin and the ring are cosmetic.
• Many-body orbits are chaotic: from three bodies on, tiny differences grow exponentially, so the long-run trajectories of the cluster preset depend on step size and are not predictions.
• The Blender models are artistic; their surface textures are procedural and not real photographs of any body.

What to try

• Run the figure-eight, then launch a mass 0.05 rock on a wide orbit: the choreography breaks.
• In the solar-type system, switch to Edge-on to see that the orbits are tilted by different angles; tap a planet to read its inclination.
• Launch a body with speed ratio 0.6: the launch point becomes the apoapsis of an ellipse with e = 1 - 0.6^2 = 0.64 (tap the new body to compare).
• Switch softening to 0 and release the cluster: watch the energy graph jump when bodies merge.