Build a Rigid-Body Physics Engine From Scratch in JavaScript

1. Rigid Bodies and Integration

A rigid body stores position, orientation, linear velocity, angular velocity, and mass. We advance the simulation using semi-implicit Euler integration (symplectic Euler) — more stable than explicit Euler for oscillatory systems:

class Body {
  constructor(mass, x, y) {
    this.pos  = { x, y };
    this.vel  = { x: 0, y: 0 };
    this.force = { x: 0, y: 0 };
    this.invMass = mass > 0 ? 1 / mass : 0; // 0 = static body
    this.angle  = 0;
    this.omega  = 0; // angular velocity [rad/s]
    this.torque = 0;
    this.invI   = mass > 0 ? 12 / (mass * 1) : 0; // I=m·L²/12 for unit square
  }

  integrate(dt) {
    // Semi-implicit Euler: update velocity first, then position
    this.vel.x += this.invMass * this.force.x * dt;
    this.vel.y += this.invMass * this.force.y * dt;
    this.pos.x += this.vel.x * dt;
    this.pos.y += this.vel.y * dt;
    this.omega  += this.invI * this.torque * dt;
    this.angle  += this.omega * dt;
    this.force.x = this.force.y = this.torque = 0; // reset forces
  }
}

2. Broad Phase: AABB Overlap Test

Testing every pair of N bodies is O(N²). The broad phase quickly eliminates non-overlapping pairs using Axis-Aligned Bounding Boxes (AABBs):

function aabbOverlap(a, b) {
  // AABB: {minX, maxX, minY, maxY}
  return a.maxX >= b.minX && b.maxX >= a.minX
      && a.maxY >= b.minY && b.maxY >= a.minY;
}

// Sweep-and-prune broad phase: sort by minX, sweep right
function broadPhase(bodies) {
  bodies.sort((a, b) => a.aabb.minX - b.aabb.minX);
  const pairs = [];
  for (let i = 0; i < bodies.length; i++)
    for (let j = i + 1; j < bodies.length; j++) {
      if (bodies[j].aabb.minX > bodies[i].aabb.maxX) break; // sorted → no more overlaps
      if (aabbOverlap(bodies[i].aabb, bodies[j].aabb))
        pairs.push([bodies[i], bodies[j]]);
    }
  return pairs;
}

3. Narrow Phase: Separating Axis Theorem (SAT)

For convex shapes, the Separating Axis Theorem states: two convex shapes do not overlap if and only if there exists a separating axis — a direction along which their projected intervals don't overlap. The candidate axes are the face normals of both shapes.

// Returns null if no collision, or {normal, depth} if collision
function satCollision(bodyA, bodyB) {
  const axes = [...getFaceNormals(bodyA), ...getFaceNormals(bodyB)];
  let minDepth = Infinity, minAxis = null;

  for (const axis of axes) {
    const [minA, maxA] = project(bodyA.vertices, axis);
    const [minB, maxB] = project(bodyB.vertices, axis);
    const overlap = Math.min(maxA, maxB) - Math.max(minA, minB);
    if (overlap <= 0) return null; // separating axis found → no collision
    if (overlap < minDepth) { minDepth = overlap; minAxis = axis; }
  }
  return { normal: minAxis, depth: minDepth };
}

function project(vertices, axis) {
  const dots = vertices.map(v => v.x * axis.x + v.y * axis.y);
  return [Math.min(...dots), Math.max(...dots)];
}

4. Impulse-Based Collision Response

When a collision is detected, we apply instantaneous velocity changes (impulses) to both bodies such that they no longer penetrate and their relative velocity along the collision normal satisfies the coefficient of restitution:

function resolveCollision(a, b, contact) {
  const { normal: n, depth, pointA: rA, pointB: rB } = contact;
  const e = 0.5; // coefficient of restitution

  const rAxn = rA.x * n.y - rA.y * n.x; // 2D cross product
  const rBxn = rB.x * n.y - rB.y * n.x;

  const vRelN =
    (b.vel.x - a.vel.x) * n.x + (b.vel.y - a.vel.y) * n.y
    + b.omega * rBxn - a.omega * rAxn;

  if (vRelN > 0) return; // bodies already separating

  const denom = a.invMass + b.invMass + rAxn**2 * a.invI + rBxn**2 * b.invI;
  const j = -(1 + e) * vRelN / denom;

  a.vel.x -= j * a.invMass * n.x;
  a.vel.y -= j * a.invMass * n.y;
  b.vel.x += j * b.invMass * n.x;
  b.vel.y += j * b.invMass * n.y;
  a.omega -= j * a.invI * rAxn;
  b.omega += j * b.invI * rBxn;

  // Positional correction (Baumgarte stabilization) — prevent sinking
  const correction = Math.max(depth - 0.01, 0) * 0.2 / (a.invMass + b.invMass);
  a.pos.x -= a.invMass * correction * n.x;
  a.pos.y -= a.invMass * correction * n.y;
  b.pos.x += b.invMass * correction * n.x;
  b.pos.y += b.invMass * correction * n.y;
}

5. Friction

After applying the normal impulse, apply a tangential impulse bounded by the friction cone (Coulomb's law: |j_t| ≤ μ · |j_n|):

function resolveFriction(a, b, contact, jNormal) {
  const mu = 0.4; // friction coefficient
  const { normal: n, pointA: rA, pointB: rB } = contact;
  const t = { x: -n.y, y: n.x }; // 2D tangent = rotate normal 90°

  const rAxT = rA.x * t.y - rA.y * t.x;
  const rBxT = rB.x * t.y - rB.y * t.x;
  const vRelT =
    (b.vel.x - a.vel.x) * t.x + (b.vel.y - a.vel.y) * t.y
    + b.omega * rBxT - a.omega * rAxT;

  const denom = a.invMass + b.invMass + rAxT**2*a.invI + rBxT**2*b.invI;
  let   jt = -vRelT / denom;

  // Clamp to friction cone
  jt = Math.max(-mu * jNormal, Math.min(mu * jNormal, jt));

  a.vel.x -= jt * a.invMass * t.x;
  a.vel.y -= jt * a.invMass * t.y;
  b.vel.x += jt * b.invMass * t.x;
  b.vel.y += jt * b.invMass * t.y;
  a.omega -= jt * a.invI * rAxT;
  b.omega += jt * b.invI * rBxT;
}

6. Angular Dynamics — Moment of Inertia

The moment of inertia I depends on shape. Common 2D values for a body of mass m:

7. Sequential Impulse Constraints (PGS)

Stacking boxes requires solving multiple simultaneous contacts. The Projected Gauss-Seidel (PGS) method iterates over all contacts multiple times, applying small impulses each pass, until convergence:

function solveContacts(contacts, iterations = 10) {
  for (let iter = 0; iter < iterations; iter++) {
    for (const contact of contacts) {
      resolveCollision(contact.a, contact.b, contact);
      resolveFriction(contact.a, contact.b, contact, contact.jNormal);
    }
  }
}

// Warm starting: keep accumulated impulses from previous frame
// → dramatically improves convergence speed and stacking stability

8. Substep Loop and Continuous Collision Detection

Fast-moving objects can tunnel through thin walls in a single timestep. Two solutions:

• Substepping: divide each frame (1/60 s) into N substeps of dt/N. More accurate, costs N× more CPU. Typical: 2–4 substeps.
• CCD (Continuous Collision Detection): sweep the body's AABB over the whole timestep, find TOI (time of impact), integrate to TOI, resolve, then integrate the remainder. Used in Cannon-es via body.ccdSpeedThreshold.

function update(dt) {
  const SUBSTEPS = 4;
  const subDt = dt / SUBSTEPS;
  for (let s = 0; s < SUBSTEPS; s++) {
    for (const body of bodies) body.applyForce({x:0, y:-9.8 * (1/body.invMass)});
    const pairs = broadPhase(bodies);
    const contacts = pairs.flatMap(([a,b]) => satCollision(a,b) ? [{a,b,...satCollision(a,b)}] : []);
    solveContacts(contacts);
    for (const body of bodies) body.integrate(subDt);
  }
}