Gravity Lab — Simulating the N-Body Problem in a Browser

Why the N-Body Problem?

The gravitational motion of two celestial bodies can be solved exactly with Kepler's laws. The moment a third body enters the picture, however, no analytical solution exists. This is the three-body problem — the same wall Newton ran into when trying to calculate the Moon's motion.

No analytical solution means you cannot use a single formula to compute an arbitrary future state. Instead, you approximate by integrating forces over very small time steps — numerical integration. The browser's RAF loop, which fires 60 times per second, serves as that integrator.

---

Body Types and the Mass System

The simulation features three types of celestial bodies. Mass differences determine gravitational strength, and each type is visually distinguished by size and color.

const BODY_DEF = {
  star: { mass: 500, radius: 16, color: '#fef08a', glow: 'rgba(254,240,138,0.95)' },
  planet: { mass: 80, radius: 9, color: '#60a5fa' },
  moon: { mass: 15, radius: 5, color: '#c4b5fd' },
};

A star weighs roughly six times a planet, and a planet about five times a moon. The actual Sun is nearly 1,000 times more massive than Jupiter, so these are heavily compressed ratios — but they produce the most interesting orbits within the confines of a browser window.

---

Gravity Calculation: F = G·m1·m2 / r² Every Frame, Every Pair

The core loop is straightforward. For N bodies, iterate over all pairs (i, j), compute the gravitational force, and accumulate velocity changes.

for (let i = 0; i < bodies.length; i++) {
  let vx = 0,
    vy = 0;
  for (let j = 0; j < bodies.length; j++) {
    if (i === j) continue;
    const b = bodies[j];
    const dx = b.x - bodies[i].x;
    const dy = b.y - bodies[i].y;
    const distSq = dx * dx + dy * dy;
    const dist = Math.sqrt(distSq);
    const angle = Math.atan2(dy, dx);
    const force = (G * b.mass) / distSq;
    vx += Math.cos(angle) * force;
    vy += Math.sin(angle) * force;
  }
  bodies[i].vx += vx;
  bodies[i].vy += vy;
}

Computational complexity is O(N²). Since body counts stay small (around 10 at most), this runs comfortably at 60 fps without further optimization.

---

Stable Orbits: Setting the Right Initial Velocity

Initial conditions are everything in a numerical simulation. To place a planet in a circular orbit around a star, the launch velocity must be calculated precisely.

The condition for a circular orbit is that centripetal acceleration equals gravitational acceleration:

// v_circular = sqrt(G * M / r)
const v = Math.sqrt((G * centralMass) / dist);
// Velocity direction is perpendicular to the radial direction
body.vx = -Math.sin(angle) * v + central.vx;
body.vy = Math.cos(angle) * v + central.vy;

The critical part is adding the central body's current velocity (central.vx, central.vy). In a binary star system where the center of mass itself is moving, omitting this causes the orbit to collapse immediately.

All three preset scenarios — solar system, binary star, and triple system — use this formula to compute initial velocities.

---

Preventing Energy Divergence: Minimum Distance Clamp

When two bodies get very close, r² approaches zero and the computed force blows up to infinity. In reality a collision would occur, but in numerical simulation the body can fly off the canvas in a single frame.

A minimum distance clamp prevents this:

const MIN_DIST_SQ = 25 * 25; // cap force below 25px separation
const distSq = Math.max(dx * dx + dy * dy, MIN_DIST_SQ);

Physically this can be interpreted as the bodies merging. Visually it looks like a smooth flyby.

---

Orbital Trail: Ring Buffer Approach

Displaying each body's recent path makes orbital shapes immediately readable. A ring buffer stores the last N positions per body.

interface Body {
  trail: Float32Array; // [x0, y0, x1, y1, ...]
  trailHead: number; // current write index
  trailLen: number; // number of stored entries
}

// Record position every frame
body.trail[body.trailHead * 2] = body.x;
body.trail[body.trailHead * 2 + 1] = body.y;
body.trailHead = (body.trailHead + 1) % TRAIL_MAX;
if (body.trailLen < TRAIL_MAX) body.trailLen++;

Using Float32Array keeps memory allocation stable and reduces GC pressure compared to plain arrays. During rendering, alpha values fade with age so older positions appear more transparent.

---

Canvas Performance — 4K Resolution Cap

The N-body simulation redraws every body, its glow effect, and its orbital trail on every frame. On 4K displays, the canvas buffer can exceed 8 million pixels, making each arc() and fillStyle assignment significantly more expensive.

A buffer area cap of 1920 × 1080 keeps the GPU workload constant regardless of display resolution:

const MAX_AREA = 1920 * 1080;
const area = w * h;
const scale = area > MAX_AREA ? Math.sqrt(MAX_AREA / area) : 1;
canvas.width = Math.round(w * scale);
canvas.height = Math.round(h * scale);

With 10 bodies, the O(N²) gravity loop itself is cheap — the real cost is rendering. Star glow effects use radialGradient fills, and each body's trail renders up to several hundred points per frame. The resolution cap reduces the per-pixel cost of all of these operations without any visible loss in visual quality.

---

Closing Thoughts

The first time I watched a planet trace an elliptical orbit after getting the numerical integrator working, I stared at the screen longer than expected. The realization that F = G·m1·m2 / r², a single line of code, produces all that motion still feels remarkable.

In the binary star scenario, two stars spiral around each other while a planet carves a complex rosette path between them — it is chaos, generated entirely by three lines of Newton's laws.