N‑BODY

# Overview

The cursor is a star. Thousands of small bodies orbit it, and each body leaves a short light trail. Slow bodies are blue, and fast bodies near the star are magenta and white. Behind them is a fixed starfield. A press of the button sends a shockwave that throws the bodies away from the star.

The effect is plain canvas 2D with no physics library. The state is four Float32Array with one value for each body: x, y, vx, vy.

The name N-BODY is not exact. In a real N-body simulation, each body attracts each other body, and the cost grows with N². This module does not do that. Only the star attracts the bodies. The bodies have no mass, do not attract each other, and do not collide. In physics terms, this is a central-force problem with N test bodies. The cost grows with N, and that is the reason that thousands of bodies are possible on the CPU.

The element count in the statistics (this.count) is N, the number of bodies. init sets it one time: 3,200, or 1,400 when the viewport is narrower than 800 px. A body that leaves the screen is not deleted. It gets a new position, so N stays constant.

# Interaction

• Move. The star is at the smoothed pointer position (pointer.sx, pointer.sy), which follows the cursor with a time constant of 0.1 s. The bodies do not move with the star. When the star moves, near bodies get new, more eccentric orbits, and far bodies fall slowly to the new centre.
• Move fast. The module does not read pointer.speed, pointer.vx, or pointer.vy. Speed has an effect only through the star position: bodies that the star goes near get a large, short acceleration and fly away on long orbits.
• Press edge. On the first frame with the button down, nova is set to 1. It decreases by 1.6 for each second, so it is zero after 0.625 s. While it is more than zero, all bodies get a strong outward acceleration, the star glow is larger, and a white ring expands from the star.
• Hold the button. The sign of the gravity changes and its strength is multiplied by 1.5. The star is a repulsor, and its glow gets an amber core. On release, the bodies fall to the star again.
• Idle. After 2.5 s with no input, the shared pointer goes into autopilot and moves on a slow path. The star follows it.

update reads the reduced flag. With reduced motion, the time step is multiplied by 0.6. The orbits, the decrease of nova, and the trail fade all run at 60% of the normal rate. The number of bodies and the forces do not change.

There is no URL parameter. The only mobile branch is the number of bodies.

# Model

Softened gravity from one centre

For each body, the code computes the vector from the body to the star, (dx, dy), and then:

d2 = dx² + dy² + soft          = r² + ε²
d  = √d2
a  = GM / d2
v += (dx / d, dy / d) · a · hdt

(dx / d, dy / d) is not a unit vector, because d is the softened distance and not r. Its length is r / d. The acceleration that the code applies is thus:

acc   = GM · (dx, dy) / (r² + ε²)^(3/2)
|acc| = GM · r / (r² + ε²)^(3/2)

This is the field of a mass that is spread out over a radius of approximately ε. Far from the star it is the Newton law GM / r². Near the star it has a maximum of approximately 940 px/s² at r = ε / √2 (21 px), and it is zero at the centre. With no softening, the acceleration goes to infinity at the centre, and one time step there throws a body off the screen.

Other accelerations

hold  :  a = −1.5 · GM / d2
nova  :  a = a − novaF · 3.5e7 / (d2 + 4000)
centre:  v += (cx − x, cy − y) · 0.05 · hdt

• The nova term is a second softened inverse-square repulsion. Its constant, 3.5e7 px³/s², is 16 times GM. Because d2 already contains soft, the denominator is r² + 4900. The maximum is approximately 4,300 px/s² at 35 px. The term acts on all bodies on the same frame. The white ring is only a drawing, not a force front.
• The centre term is a weak spring to the centre of the viewport: 0.05 px/s² for each px of distance. It is not gravity. It keeps the swarm on the screen. For a star at the centre, spring and gravity are equal at approximately 350 px (GM / r² = 0.05 · r). Farther out, the spring is the stronger force.

After these terms, the speed is limited to 1,200 px/s.

Integration

sdt = min(dt, 1/30) · (reduced ? 0.6 : 1)
hdt = sdt / sub                    sub = 2

two times:
v = v + acc(x) · hdt
x = x + v · hdt

Explicit Euler moves the position with the old velocity, and each orbit then gains energy and becomes a spiral that goes out. The semi-implicit method is symplectic: the energy error stays bounded, and orbits stay closed, at the same cost. It is not leapfrog, because there is no half step between velocity and position.

The two sub-steps make the step 1/120 s at 60 Hz. The error is largest near the star, where the acceleration changes quickly across a short distance. The star position is constant during the two sub-steps.

Start conditions

Each body starts at a random angle a and at a radius r = 80 + random · 0.55 · min(w, h) from the star. The radius is uniform, not the area, so the disc is more dense near the star.

v  = √(GM / r) · (0.85 + 0.3 · random)
vx = −sin(a) · v
vy =  cos(a) · v

√(GM / r) is the speed at which gravity is equal to the centripetal acceleration v² / r: 166 px/s at 80 px, 86 px/s at 300 px. The direction (−sin a, cos a) is at 90° to the radius. All bodies turn in the same direction, clockwise on the screen because the y axis points down. A factor of 1 gives a circle, and other factors give ellipses. The escape speed is √2 times the circular speed, and 1.15 is below that, so all start orbits are closed.

Too far and too near

A body that is more than 200 px out of the viewport gets a new state. It starts 10 px out of a random edge with a speed v of 60 to 180 px/s: 0.4 · v in the direction of the star, plus v at 90° to that direction with a random sign. These bodies thus turn in the two directions, and they go around the star and not into it.

There is no limit for a body that goes near the star. It is not deleted and it does not merge. The softening limits its acceleration, the speed limit bounds its velocity, and it goes through the centre and out on the other side.

Colour by speed

The speed gives an index from 0 to 63: min(63, floor(speed / 620 · 63)). The index selects a colour from a table that goes from deep blue through cyan and magenta to white. The same index sets the alpha (0.45 to 1) and the size (1.2 to 2.8 px). Orbit speed decreases with distance, so the colour shows the field: a white core, a cyan ring, a blue outer region.

# Frame pipeline

1. Time stepreads dt, reducedwrites sdt, hdt

   dt is limited to 1/30 s and multiplied by 0.6 for reduced motion. hdt is one half of the result.

2. Nova statereads pointer.down, wasDown, novawrites nova, wasDown, novaF, repel

   A press edge sets nova to 1. novaF keeps the value for this frame, and nova then decreases by 1.6 · sdt.

3. Integrate, two sub-stepsreads pointer.sx, pointer.sy, x, y, vx, vywrites x, y, vx, vy

   For each body: softened gravity or repulsion, nova term, centre spring, speed limit, then the position step.

4. Fade the trailsreads trail canvaswrites trail canvas

   A full-canvas destination-out fill with alpha 1 − 0.82^(60 · sdt) removes that part of the alpha of each pixel. The blend mode is then set to lighter.

5. Respawn and sortreads x, y, vx, vywrites buckets, state of respawned bodies

   A body more than 200 px out of the viewport gets a new state and is not drawn on this frame. Each other body goes into one of 64 buckets by its speed.

6. Draw the bodiesreads buckets, lut, x, ywrites trail canvas

   For each bucket that is not empty: set colour, alpha, and size one time, then draw one square for each body.

7. Draw the starreads pointer.sx, pointer.sy, pointer.down, novaFwrites trail canvas

   A radial gradient of radius 70 + 260 · novaF px, plus 30 px while the button is down. While novaF is more than zero, a white ring of radius (1 − novaF) · 520 px is drawn. The blend mode then goes back to source-over.

The starfield canvas is not in this list. It is drawn in resize only: a dark fill, one radial gradient for the nebula, and 700 stars of 0.8, 1.4, or 2.2 px. resize also clears the trail canvas. It does not move the bodies.

# Code

Start orbits

The direction vector (−sin a, cos a) is the position vector (cos a, sin a) turned by 90°.

// ring around the cursor with orbital velocity
const a = Math.random() * Math.PI * 2, r = 80 + Math.random() * Math.min(w, h) * 0.55;
this.x[i] = p.sx + Math.cos(a) * r; this.y[i] = p.sy + Math.sin(a) * r;
const v = Math.sqrt(this.GM / r) * (0.85 + Math.random() * 0.3);
this.vx[i] = -Math.sin(a) * v; this.vy[i] = Math.cos(a) * v;

Respawn from an edge

The + 1 in d prevents a division by zero. The tangential part (−dy / d, dx / d) is 2.5 times the radial part, so the new body goes into an orbit and does not fall straight in.

// enter from a random screen edge, drifting inward
const side = (Math.random() * 4) | 0, s = Math.random();
if (side === 0) { this.x[i] = -10; this.y[i] = s * h; } else if (side === 1) { this.x[i] = w + 10; this.y[i] = s * h; }
else if (side === 2) { this.x[i] = s * w; this.y[i] = -10; } else { this.x[i] = s * w; this.y[i] = h + 10; }
const dx = p.sx - this.x[i], dy = p.sy - this.y[i], d = Math.hypot(dx, dy) + 1;
const v = 60 + Math.random() * 120, side_ = Math.random() < 0.5 ? 1 : -1;
this.vx[i] = (dx / d) * v * 0.4 + (-dy / d) * v * side_; this.vy[i] = (dy / d) * v * 0.4 + (dx / d) * v * side_;

Press edge and hold

wasDown makes the difference between the press edge and the hold. novaF is read before nova decreases, so the first frame of the shockwave has the full value of 1.

if (p.down && !this.wasDown) this.nova = 1;   // press edge → shockwave
this.wasDown = p.down;
const novaF = this.nova; this.nova = Math.max(0, this.nova - sdt * 1.6);
const repel = p.down ? 1 : 0;

The force and the step

The velocity lines come before the position line: that order makes the method semi-implicit. There is no loop across pairs of bodies. The only position that a body reads, apart from its own, is sx, sy.

for (let s = 0; s < sub; s++) {
  for (let i = 0; i < N; i++) {
    let dx = sx - x[i], dy = sy - y[i];
    const d2 = dx * dx + dy * dy + soft, d = Math.sqrt(d2);
    let a = GM / d2;
    // hold: the star flips to a repulsor; the press itself sends a shockwave
    if (repel) a = -a * 1.5;
    if (novaF > 0) a -= novaF * 3.5e7 / (d2 + 4000);
    vx[i] += dx / d * a * hdt; vy[i] += dy / d * a * hdt;
    // gentle pull to screen centre keeps the swarm on screen
    vx[i] += (cx - x[i]) * 0.05 * hdt; vy[i] += (cy - y[i]) * 0.05 * hdt;
    // speed cap
    const sp2 = vx[i] * vx[i] + vy[i] * vy[i];
    if (sp2 > 1.44e6) { const k = 1200 / Math.sqrt(sp2); vx[i] *= k; vy[i] *= k; }
    x[i] += vx[i] * hdt; y[i] += vy[i] * hdt;
  }
}

Trail fade

In destination-out, the colour of the fill is not used. Only its alpha is.

// fade trails
ctx.globalCompositeOperation = 'destination-out';
ctx.fillStyle = `rgba(0,0,0,${1 - Math.pow(0.82, sdt * 60)})`;
ctx.fillRect(0, 0, w, h);
ctx.globalCompositeOperation = 'lighter';

Draw by speed bucket

The bodies are sorted into 64 arrays first, so fillStyle and globalAlpha are set 64 times at most, not 3,200 times. The arrays are made one time and emptied with length = 0. The respawn test is in the same loop.

// bodies, batched by colour bucket
const lut = this.lut, buckets = this.buckets || (this.buckets = Array.from({ length: 64 }, () => []));
for (let b = 0; b < 64; b++) buckets[b].length = 0;
const margin = 200;
for (let i = 0; i < N; i++) {
  if (x[i] < -margin || x[i] > w + margin || y[i] < -margin || y[i] > h + margin) { this.spawn(i, false); continue; }
  const sp = Math.sqrt(vx[i] * vx[i] + vy[i] * vy[i]);
  const b = Math.min(63, (sp / 620 * 63) | 0);
  buckets[b].push(i);
}
for (let b = 0; b < 64; b++) {
  const arr = buckets[b]; if (!arr.length) continue;
  ctx.fillStyle = lut[b];
  ctx.globalAlpha = 0.45 + b / 63 * 0.55;
  const sz = 1.2 + b / 63 * 1.6;
  for (let k = 0; k < arr.length; k++) { const i = arr[k]; ctx.fillRect(x[i] - sz / 2, y[i] - sz / 2, sz, sz); }
}

# Tuned constants

# Performance

All the work is on the CPU, and it is linear in the number of bodies. There is no pair loop.

• Scales with N. The force loop runs 2 × 3,200 = 6,400 times for each frame, with one square root each time. The draw loop adds 3,200 square roots for the colour index and 3,200 fillRect calls, with 64 state changes at most.
• Scales with pixels. The destination-out fill touches each pixel of the trail canvas on each frame. The star gradient covers a square of 140 px, and as much as 720 px at the start of a nova. The browser also blends the transparent trail canvas on the starfield canvas.

Why two canvases: the usual way to fade trails is to fill the canvas with the background colour at a low alpha. That makes the canvas opaque, so a background picture must be drawn again on each frame. Here the background is a gradient and 700 stars. destination-out decreases the alpha of the trail pixels and adds no colour. The trail canvas fades to transparent, and the starfield below it is drawn only in resize. The price is one more full-size canvas.

Limits that the module applies:

• The device pixel ratio is limited to 2. The value is set in init. resize does not read its dpr argument.
• The starfield canvas is opaque (alpha: false). Only the trail canvas has an alpha channel.
• 1,400 bodies and not 3,200 when the viewport is narrower than 800 px at init. A later resize does not change N.
• The frame loop makes no new arrays.
• There is no resolution limit other than the dpr limit, and no fallback path.

# References

• Wikipedia, N-body simulation, section Softening. The source of the force law: r² + ε² replaces r². The page also describes the real N-body problem with pair forces, which this module does not compute.
• Wikipedia, Plummer model. The mass distribution whose field is the exact form in the code.
• Wikipedia, Semi-implicit Euler method. The integration method, and why it keeps the energy bounded.
• Wikipedia, Circular orbit. The derivation of the start speed √(GM / r).
• MDN, globalCompositeOperation. The destination-out mode for the fade and the lighter mode for the glow.