How to Simulate Gravity in JavaScript (In About 40 Lines)

the whole thing is one equation

Newton worked this out in the 1680s and it has stubbornly refused to get more complicated since:

F = G · m₁ · m₂ / r²

Two things pull on each other. The pull gets stronger with mass, and weaker with the square of the distance. That's it. That's the physics. Everything else in a gravity simulator is bookkeeping.

I find that genuinely lovely, and slightly annoying. You spend a week fighting a shader and then discover the actual universe fits on one line.

step one: every pair pulls on every other pair

With two bodies you work out one force. With three you work out three. With ten you work out forty-five, because every body pulls on every other body and you'd be double-counting if you weren't careful.

The standard trick is to start the inner loop one past the outer one, so each pair is visited exactly once:

for (let i = 0; i < bodies.length; i++) {
  for (let j = i + 1; j < bodies.length; j++) {
    const b1 = bodies[i];
    const b2 = bodies[j];
    const dx = b2.x - b1.x;
    const dy = b2.y - b1.y;
    const distSq = dx * dx + dy * dy;
    ...
  }
}

Note distSq. Distance squared, and no square root yet. The force equation wants r² anyway, and Math.sqrt is the expensive bit. You'll need the actual distance in a moment to work out direction, but there's no sense computing it before you know you need it.

Then Newton's third law does you a favour. Every pair gives you two forces for the price of one, equal and opposite:

const force = (G * b1.mass * b2.mass) / (distSq + EPSILON * EPSILON);
const dist = Math.sqrt(distSq);
const fx = force * (dx / dist);
const fy = force * (dy / dist);

forces[i].fx += fx; forces[i].fy += fy;
forces[j].fx -= fx; forces[j].fy -= fy;

Dividing dx by dist turns the offset into a unit direction, so multiplying by force splits the pull into its x and y parts. And += on one body, -= on the other. Same strength, opposite way. Half your work, done for you, by a man who died in 1727.

That EPSILON in the denominator is doing something important and I've written about why a gravity simulator needs a softening parameter separately, because it deserves more than a footnote. Short version: without it, two bodies that get very close divide by nearly zero and leave the solar system.

step two: force becomes movement

Once every body has its total force, the rest is the most famous equation in physics rearranged slightly. F = ma, so a = F/m:

const ax = forces[i].fx / b.mass;
const ay = forces[i].fy / b.mass;

b.vx += ax * dt;
b.vy += ay * dt;
b.x += b.vx * dt;
b.y += b.vy * dt;

Acceleration changes velocity, velocity changes position. Four lines and you have orbits.

Except (and this is the part that took me an unreasonably long time to appreciate) the order of those four lines matters enormously. Velocity is updated first, and then position is moved using the brand-new velocity rather than the old one. Swap those two pairs around and your orbits slowly spiral outward until everything drifts off screen. Same maths, same number of lines, completely different universe. That one's worth its own post too.

step three: pick your units and stop apologising

Here's where I think a lot of tutorials lose people. They reach for real values (G is 6.674 × 10⁻¹¹, Earth is 5.97 × 10²⁴ kg) and then everything is in metres, one pixel is a continent, and nothing visible happens for four hundred years.

I set G = 1 and made up masses. A planet is 500, a star is 3000, a black hole is 10000. None of those are kilograms. They're just numbers that produce nice orbits at a size that fits on a screen.

This bothered me for about a day, until I noticed the maths doesn't care. Gravity is scale-invariant in the ways that matter here: change the units and you change the numbers, not the shapes. The orbit is still an ellipse. You still get the same beautiful mess with three bodies. And when you want a circular orbit you can still work out the speed exactly:

// A body at distance 250 from a 15000-mass sun, with G = 1
bodies.push(new Body(250, 0, 0, Math.sqrt(15000 / 250), 300, 8, '#4dabf7'));

That's the real orbital velocity formula, v = √(GM/r), doing honest work in a made-up unit system. Very freeing, once you let yourself.

Though "circular" deserves an asterisk I didn't give it originally. That speed is exactly right for a circular orbit around a fixed sun with nothing else in the system. My shipped preset has neither: the sun is free to move and recoils, and there are three planets tugging on each other. Run it and that body's distance actually swings between about 141 and 296 rather than sitting at 250. The formula is right; the situation I dropped it into isn't the one the formula assumes. Pin the sun and remove the siblings and it holds to within a couple of units over thousands of frames.

the bit I got wrong

I'd like to end on a clean note but that would be dishonest, and there's a genuine mistake sitting in my collision handling that I only spotted while writing this.

When two bodies touch, I merge or bounce them, and then:

if (distSq < (b1.radius + b2.radius) ** 2) {
  handleCollision(i, j);
  return; // Exit current step, arrays modified
}

The reasoning was sound. handleCollision rebuilds the bodies array, so carrying on through a loop that's indexing into the old one is asking for trouble. Bailing out is the safe move.

The problem is how far it bails. That return leaves updatePhysics entirely (not the inner loop, the whole function) so on any frame where a collision happens, the force-application loop underneath never runs at all. Nothing moves. The entire universe holds perfectly still for one frame while two rocks sort themselves out.

At four bodies you will never see it. Load the galaxy preset, though (82 bodies constantly bumping into each other) and you're dropping frames of motion regularly. It's not a crash, it's not an error, it's just a universe with a slight stutter that I'd never have found by looking at it.

Before fixing it I wanted to know how bad it actually was, so I stubbed out the canvas, ran updatePhysics six hundred times headlessly, and counted frames where a body that should have moved didn't. Over three collision modes on the galaxy preset: 161 stalled frames. Elastic was the worst at 97 out of 600, which surprised me. I'd assumed elastic was fine, because the two colliding bodies do get pushed apart on the way out of the function, so something visibly moves. The other eighty bodies just quietly didn't.

That's the argument for measuring before fixing, by the way.

The fix is to stop resolving collisions in the middle of the loop. Queue them instead (holding references to the bodies rather than their indices, since the array is about to be rebuilt underneath them) finish the force pass, move everything, and drain the queue at the end:

if (distSq < (b1.radius + b2.radius) ** 2) {
  collisions.push([b1, b2]);
  continue; // No mutual gravity for an overlapping pair.
}

One wrinkle worth mentioning: a body can end up in two collisions in the same frame, which happens constantly in a crowded galaxy. If you don't guard against it, the second pair merrily merges a body the first pair already consumed, and you've duplicated mass out of nothing. A Set of already-removed bodies, checked before each resolution, handles it.

Same test after the change: zero stalled frames across all six preset-and-mode combinations, collisions still merging, no NaN anywhere. I ran it against the old code first to make sure the test could actually detect the problem. A green test that was never red isn't evidence of anything.

where to go next

If you want to see all of this actually running, the gravity sandbox is free and needs no signup. Click to place bodies, drag to launch them. And if you're curious why three bodies is where the maths stops being polite, that's the three-body problem, which I can demonstrate very convincingly because my own figure-eight preset doesn't work.