N-Body Gravity
Newton's law is a sentence. The consequences take forever to unfold. Two bodies: exactly solvable, orbits closed, elegant. Three bodies: Poincaré proved in 1889 that no general analytic solution exists — the phase space is chaotic almost everywhere. Yet in 2000, Chenciner and Montgomery found a needle in that haystack: three equal masses tracing a perfect figure-8, periodic, choreographic, improbable. It runs below. Drag the time scale up and watch the trails carve the orbit. Then add your own body and watch what happens to the mission parameters.
The equations
Each body i experiences the gravitational attraction of every other body j:
Fij = G mi mj / rij² · r̂ij
ai = Σj≠i G mj (rj − ri) / (|rj − ri|² + ε²)3/2
The softening parameter ε = 0.03 prevents the acceleration from diverging at close approach — a standard regularisation used in cosmological simulations. It introduces a small bias to the potential (equivalent to extending point masses into smooth spheres of radius ε), but keeps the integration stable across encounters.
The figure-8
The Chenciner–Montgomery choreography is a periodic orbit of the equal-mass 3-body problem in which all three bodies follow the same figure-8 curve, separated in phase by T/3. The initial conditions (in normalised units with G = 1, m = 1) are:
r₁ = ( 0.97000436, −0.24308753), ṙ₁ = (0.46620369, 0.43236573)
r₂ = (−0.97000436, 0.24308753), ṙ₂ = (0.46620369, 0.43236573)
r₃ = (0, 0), ṙ₃ = (−0.93240737, −0.86473146)
The period is T ≈ 6.3259. Unlike generic 3-body trajectories — which are chaotic almost everywhere in phase space — the figure-8 is orbitally stable to small perturbations, though you can break it decisively by adding a fourth body.
Conservation laws
With no external forces, total linear momentum, total angular momentum L, and total energy E = KE + PE
are conserved exactly by the continuous equations. The RK4 integrator conserves them only approximately;
the |ΔE / E₀| readout tracks energy drift over time. At the default settings
(Δt = 0.002, 5 substeps/frame) drift over a full figure-8 period stays below 10⁻³ %. Increase the
time scale and watch the budget tighten.
Binary star preset
Two stars of mass m = 2 sit at x = ±0.5 in a circular mutual orbit. Their orbital speed follows from balancing gravity against centripetal acceleration: v = √(Gm / 4r) = 1.0. A low-mass planet orbits much further out at r = 2.2, where the binary looks like a single point mass of M = 4, so its circular speed is v ≈ √(GM/r) ≈ 1.35. The planet's orbit is stable while r ≫ binary separation — drag its starting distance inward and watch stability collapse.
Numerical note
Integration is 4th-order Runge–Kutta with Δt = 0.002 per substep. Each animation frame advances max(1, round(5 × time_scale)) substeps. Trail points are stored in a circular buffer of 700 entries per body; the trail length slider controls how many are drawn. The centre-of-mass crosshair stays fixed for closed systems — any drift signals numerical error, not physics.