Love the presets! The Broucke one is my favourite

Love this.

It's not obvious from the UI but you can enter small mass changes and watch things slowly fall apart. E.g. 1.0001 will work even though the UI displays 1.0 after you hit enter.

As the default simulation played out beautifully on-load, I immediately started to question: "hang on, I thought there wasn't a solution for 3 bodies, but this looks stable".

Before I could complete the thought, it fell apart magnificently :)

Discussed before here:

https://news.ycombinator.com/item?id=45967079 (245 and 112 comments)

Wow, this looks really nice on a 240 Hz display.

Here is a thought; Instead of using F=ma, use the equations of motion from GR:

================================================================================ SCHWARZSCHILD METRIC AND GEODESIC EQUATIONS OF MOTION (SUMMARY) ================================================================================

I. THE SCHWARZSCHILD METRIC (g_uv)

The spacetime geometry is defined by the *line element*, ds^2, which relates coordinate changes (dt, dr, d(phi), etc.) to physical distance or proper time: ds^2 = g_uv * dx^u * dx^v

For the Schwarzschild vacuum solution, the line element in the equatorial plane (theta = pi/2) is: ds^2 = -(1 - r_s / r) * c^2 * dt^2 + (1 - r_s / r)^(-1) * dr^2 + r^2 * d(phi)^2

The corresponding non-zero metric components (g_uv) are: g_tt = -(1 - r_s / r) * c^2 g_rr = 1 / (1 - r_s / r) g_phiphi = r^2

Where: r_s = 2 G * M / c^2 (Schwarzschild Radius)

The Lagrangian L for the geodesic path is constructed directly from the metric: L = (1/2) * [ g_tt * (dt/d(lambda))^2 + g_rr * (dr/d(lambda))^2 + g_phiphi (d(phi)/d(lambda))^2 ]

--------------------------------------------------------------------------------

II. CONSERVATION LAWS (FROM EULER-LAGRANGE EQUATIONS)

A. TIME EOM (Conserved Energy E) Since the metric is time-independent, the quantity conjugate to t is conserved: *Specific Energy (E)*.

EQUATION (1): Time Evolution d(t)/d(lambda) = E / ( c^2 * (1 - r_s / r) )

B. PHI EOM (Conserved Angular Momentum L_z) Since the metric is symmetric with respect to phi, the quantity conjugate to phi is conserved: *Specific Angular Momentum (L_z)*.

EQUATION (2): Angular Evolution d(phi)/d(lambda) = L_z / r^2

--------------------------------------------------------------------------------

III. RADIAL EQUATION OF MOTION (FROM THE METRIC CONSTRAINT)

The radial EOM is derived by imposing the metric normalization condition (g_uv * u^u * u^v = epsilon).

A. MASSIVE PARTICLES (Mass m > 0) The proper time (tau) is the affine parameter (lambda=tau), and the normalization is epsilon = c^2. The final EOM is: (dr/d(tau))^2 = E^2/c^2 - V_eff^2

EQUATION (3M): Radial EOM (Massive) (dr/d(tau))^2 = E^2/c^2 - c^2 * (1 - r_s/r) * ( 1 + L_z^2 / (c^2 * r^2) )

B. MASSLESS PARTICLES (Mass m = 0) The normalization is epsilon = 0. The final EOM is: (dr/d(lambda))^2 = E^2 - V_eff^2

EQUATION (3P): Radial EOM (Massless / Photon) (dr/d(lambda))^2 = E^2 - (1 - r_s/r) * L_z^2 / r^2

--------------------------------------------------------------------------------

IV. SUMMARY OF GEODESIC EQUATIONS OF MOTION (EOM)

The motion of any particle (massive or massless) in the Schwarzschild spacetime is determined by the following three coupled first-order differential equations:

A. TIME EVOLUTION: d(t)/d(lambda) = E / ( c^2 * (1 - r_s / r) )

B. ANGULAR EVOLUTION: d(phi)/d(lambda) = L_z / r^2

C. RADIAL EVOLUTION (Specific): 1. Massive Particle (using d(tau)): (dr/d(tau))^2 = E^2/c^2 - c^2 * (1 - r_s/r) * ( 1 + L_z^2 / (c^2 * r^2) )

2. Massless Particle (using d(lambda)): (dr/d(lambda))^2 = E^2 - (1 - r_s/r) * L_z^2 / r^2

================================================================================

This also holds for a non-rotating black hole.