10.1 Rigid Bodies & Collision

Everything so far has been about a single object moving through space: its position, its rotation, the vectors that describe it. This chapter is about what happens when many of those objects share the same space and start bumping into each other — a ball landing on the floor, a crate sliding into a wall, a stack of boxes staying stacked instead of falling through each other. That is rigid body physics, and it rests entirely on two things you already have: the calculus chapter's semi-implicit Euler stepping and fixed timestep, and the geometry chapter's vectors, dot products, and AABBs (axis-aligned bounding boxes). In fact, the geometry chapter's tip box about circles and boxes already used the exact words "broad phase" and "narrow phase" without explaining them. This chapter is that promise, paid off in full.

Every idea below follows the usual shape: a small runnable C++ program, its real printed output, then a plain explanation of what happened. Sections 1 through 11 build a tiny physics engine completely from scratch in C++, reusing the Vec2 struct and add/sub/scale/dot/length functions from the vectors chapter, so you can see exactly what is happening under the hood. Section 12 then switches to C# and shows how Unity's Rigidbody and Collider components do this same work for you, so you can recognize every concept inside the tool you will actually use day to day.

1. What is a rigid body?

A rigid body is an object whose shape never changes — the distance between any two points on it stays fixed no matter how it moves or spins. A crate, a ball, a car chassis: all rigid bodies. Cloth, rope, and skin that squashes are not — those are soft bodies, a different and harder topic this chapter does not cover. Restricting ourselves to shapes that never deform is exactly what makes rigid body physics tractable: instead of tracking every point on the object separately, you only need to track a handful of numbers for the whole object at once.

Those numbers are the rigid body's state:

• Position — where the body's reference point (usually its center of mass) is in the world, a vector.
• Rotation (also called orientation) — which way the body is currently turned. In 2D this is a single angle; in 3D, as the rotation chapter covered, engines store it as a quaternion to avoid gimbal lock.
• Linear velocity — how fast the position is changing, a vector, in units per second.
• Angular velocity — how fast the rotation is changing. In 2D a single number (radians per second); in 3D a vector — an axis to spin around plus a speed.
• Mass — how much the body resists a push. A bowling ball needs more force than a tennis ball to reach the same speed.

Here is that state as a C++ struct, the foundation every later section builds on:

#include <iostream>

struct Vec2 { float x, y; };

struct RigidBody {
    Vec2 position;
    float rotation;         // 2D: a single angle, in radians
    Vec2 velocity;          // linear velocity, units per second
    float angularVelocity;  // 2D: radians per second
    float mass;
    float invMass;          // 1 / mass; 0 means "infinite mass", i.e. never moves
};

RigidBody makeBody(Vec2 pos, float mass) {
    RigidBody b;
    b.position = pos;
    b.rotation = 0.0f;
    b.velocity = {0.0f, 0.0f};
    b.angularVelocity = 0.0f;
    b.mass = mass;
    b.invMass = (mass > 0.0f) ? 1.0f / mass : 0.0f;
    return b;
}

int main() {
    RigidBody box   = makeBody({0.0f, 5.0f}, 2.0f);
    RigidBody floor = makeBody({0.0f, 0.0f}, 0.0f);   // mass 0 = static, never moves

    std::cout << "box:   mass=" << box.mass   << " invMass=" << box.invMass   << "\n";
    std::cout << "floor: mass=" << floor.mass << " invMass=" << floor.invMass << "\n";
}

Output:

box:   mass=2 invMass=0.5
floor: mass=0 invMass=0

Notice invMass (inverse mass, 1/mass) instead of storing mass alone. Every later section divides by "how much do these two bodies resist moving, combined" — and a wall or floor should never move no matter how hard something hits it. Representing "infinite mass" directly is impossible (you cannot divide by infinity), but representing "zero inverse mass" is trivial, and multiplying anything by zero cleanly gives zero movement. This one trick is why invMass, not mass, shows up in every formula for the rest of this chapter. A body with invMass = 0 is called static; anything else is dynamic.

2. Integrating motion: semi-implicit Euler, applied to a body

Recall the calculus chapter: a physics engine never solves the exact equations of motion. It approximates, one small fixed step at a time, using semi-implicit (symplectic) Euler — update velocity first from the current acceleration, then update position using the new velocity, not the old one. That one-line ordering is what keeps energy bounded instead of leaking away or exploding, and it is exactly what every real physics engine does every single step. Here it is applied to a RigidBody, stepping under gravity:

#include <iostream>

struct Vec2 { float x, y; };

struct RigidBody {
    Vec2 position;
    Vec2 velocity;
    float invMass;
};

void integrate(RigidBody& body, Vec2 gravity, float dt) {
    if (body.invMass <= 0.0f) return;         // static bodies never move

    // 1. velocity FIRST, using the current acceleration (here: just gravity)
    body.velocity.x += gravity.x * dt;
    body.velocity.y += gravity.y * dt;

    // 2. position using the NEW velocity -- this is what makes it semi-implicit
    body.position.x += body.velocity.x * dt;
    body.position.y += body.velocity.y * dt;
}

int main() {
    RigidBody ball;
    ball.position = {0.0f, 0.0f};
    ball.velocity = {3.0f, 0.0f};
    ball.invMass  = 1.0f;

    Vec2 gravity = {0.0f, -10.0f};
    float dt = 0.1f;

    for (int tick = 1; tick <= 4; tick++) {
        integrate(ball, gravity, dt);
        std::cout << "tick " << tick << ": pos=(" << ball.position.x << ", " << ball.position.y
                   << ") vel=(" << ball.velocity.x << ", " << ball.velocity.y << ")\n";
    }
}

Output:

tick 1: pos=(0.3, -0.1) vel=(3, -1)
tick 2: pos=(0.6, -0.3) vel=(3, -2)
tick 3: pos=(0.9, -0.6) vel=(3, -3)
tick 4: pos=(1.2, -1) vel=(3, -4)

If step 2 instead used the old velocity to move position (explicit Euler, the unstable variant the calculus chapter warned about), the y positions after the same four ticks would be 0, -0.1, -0.3, -0.6 instead of -0.1, -0.3, -0.6, -1.0 — a full step behind, because it is always moving position with a velocity that has not "seen" this step's gravity yet. Same amount of arithmetic, same input numbers, meaningfully different (and less stable) trajectory. This is exactly the one-line ordering the game loop chapter's FixedUpdate example already used, without naming it — updating velocityY before using it to move transform.position.

Explicit Euler: watch the energy blow up

Why insist on the "velocity first" ordering? Because the other order — explicit (forward) Euler, which moves position with the old velocity and only then updates velocity — does not merely lag by a step; on anything that oscillates it silently adds energy every step until the simulation explodes. The clearest place to see it is a spring: a mass pulled back toward the origin with acceleration a = -k*x, which in the real world swings forever at a constant amplitude. Run both integrators on that exact spring and print the total energy (0.5*(v*v + k*x*x), which physically must stay constant):

#include <cstdio>
#include <cmath>

int main() {
    float k = 1.0f, dt = 0.2f;      // spring constant, fixed step
    float ex = 1.0f, ev = 0.0f;     // explicit Euler: position, velocity
    float sx = 1.0f, sv = 0.0f;     // semi-implicit Euler: position, velocity

    printf("step | explicit x   energy | semi-impl x   energy\n");
    for (int step = 0; step <= 40; step++) {
        if (step % 8 == 0) {
            float eE = 0.5f * (ev * ev + k * ex * ex);
            float sE = 0.5f * (sv * sv + k * sx * sx);
            printf("%3d  | %9.4f  %.4f | %9.4f  %.4f\n",
                   step, ex, eE, sx, sE);
        }
        // explicit (forward) Euler: move x with OLD v, then v with OLD x
        float exNew = ex + ev * dt;
        float evNew = ev - k * ex * dt;
        ex = exNew; ev = evNew;
        // semi-implicit: v FIRST (current x), then x with the NEW v
        sv = sv - k * sx * dt;
        sx = sx + sv * dt;
    }
}

Output:

step | explicit x   energy | semi-impl x   energy
  0  |    1.0000  0.5000 |    1.0000  0.5000
  8  |   -0.0098  0.6843 |   -0.1323  0.5133
 16  |   -1.3684  0.9365 |   -0.9916  0.4936
 24  |    0.0402  1.2817 |    0.1955  0.5196
 32  |    1.8719  1.7540 |    0.9791  0.4875
 40  |   -0.0917  2.4005 |   -0.2580  0.5256

The spring should hold a constant energy of 0.5 forever. Semi-implicit Euler (right two columns) does exactly that — its energy wobbles a hair around 0.5 (0.49 to 0.53) and never drifts, so the mass keeps swinging between about -1 and +1 for as long as you run it. Explicit Euler (left two columns) instead climbs without stopping: 0.50, 0.68, 0.94, 1.28, 1.75, 2.40 — nearly five times the starting energy after only 40 steps, its swing growing wider every cycle. Give it a few hundred more steps and the numbers run off to infinity. This is not a rounding artifact you can shrink away by using double; it is a directional bias baked into the ordering. It is the reason no shipping engine uses explicit Euler for dynamics, and why the calculus chapter's one-line reordering matters so much.

Verlet integration: the velocity you never store

A third integrator is worth knowing because a whole family of effects — cloth, rope, hair, and position-based dynamics — is built on it: Verlet integration. Its trick is to not store velocity at all. Instead it keeps the previous position and infers motion from the gap between "where I was" and "where I am":

#include <cstdio>

int main() {
    float g = -10.0f, dt = 0.1f;   // gravity, fixed step
    float y     = 0.0f;            // current position (start at rest at 0)
    float yPrev = 0.0f;            // previous position; at rest it equals y

    printf("step | y (Verlet) | implied v = (y - yPrev)/dt\n");
    for (int step = 1; step <= 5; step++) {
        // position Verlet: next = 2*current - previous + a*dt*dt
        float yNext = 2.0f * y - yPrev + g * dt * dt;
        float impliedV = (yNext - y) / dt;   // velocity is never stored -- it is inferred
        yPrev = y;
        y = yNext;
        printf("%3d  | %8.2f  | %8.1f\n", step, y, impliedV);
    }
}

Output:

step | y (Verlet) | implied v = (y - yPrev)/dt
  1  |    -0.10  |     -1.0
  2  |    -0.30  |     -2.0
  3  |    -0.60  |     -3.0
  4  |    -1.00  |     -4.0
  5  |    -1.50  |     -5.0

Those y values — -0.10, -0.30, -0.60, -1.00, -1.50 — are the exact same fall semi-implicit Euler produced two listings above, reached without ever storing a velocity variable. Verlet is symplectic too (it will not blow the spring up), and it has one property that makes it the default for cloth and rope: because there is no separate velocity to contradict the position, you can grab a point and move it — pin a flag's corner, snap a rope back to its maximum length — and the next step's (current - previous) automatically turns that change into a sensible velocity. Correcting positions directly and letting the integrator recover velocity for free is the whole idea behind position-based dynamics.

Why physics insists on a fixed step

The game loop chapter showed Update's deltaTime jittering frame to frame, and the calculus chapter built a full fixed-timestep-plus-accumulator loop to fix it. Collision detection is exactly why physics cannot skip that machinery. Three concrete reasons stack on top of each other:

• Stability. A larger, uneven dt makes any Euler-family integrator's error grow faster — a single lag spike can make one step wildly wrong.
• Determinism. The exact same starting state, stepped with the exact same fixed dt every time, always produces the exact same result — critical for replays, and for keeping two players' simulations agreeing in a networked game.
• Tunneling risk. A fast object's travel-per-step is speed * dt. A bigger, unpredictable dt means bigger, unpredictable jumps — and section 10 of this very chapter shows exactly how a big enough jump can skip clean over a thin wall undetected. A fixed step keeps that jump size fixed and known, which is the only thing that makes the fix in section 10 possible at all.

So every body in the world gets integrated by the same fixed dt inside the same physics step — which sets up the rest of this chapter perfectly: once every body's new position for this step is known, the engine has to figure out which of them now overlap.

The accumulator: one wobbly frame time, whole fixed steps

"Use a fixed dt" runs straight into a wall: real frames do not arrive on a fixed schedule. One frame takes 12 ms, the next 22 ms, the next 9 ms, depending on what the GPU and the rest of the game are doing that instant. If you simply passed each frame's real duration to integrate, physics would step by a different dt every frame — and Euler's answer depends on dt, so the same motion would land in a different place depending on how the frames happened to fall:

#include <iostream>

// integrate a falling body for a total of 1.0s, split into `steps` equal parts
float fallFor(float totalTime, int steps) {
    float y = 0.0f, v = 0.0f, g = -10.0f;
    float dt = totalTime / steps;
    for (int i = 0; i < steps; i++) {
        v += g * dt;      // semi-implicit Euler
        y += v * dt;
    }
    return y;
}

int main() {
    std::cout << "1.0s of gravity, stepped different ways:\n";
    std::cout << "  1 step  of 1.00s : y = " << fallFor(1.0f, 1) << "\n";
    std::cout << "  4 steps of 0.25s : y = " << fallFor(1.0f, 4) << "\n";
    std::cout << "  10 steps of 0.10s: y = " << fallFor(1.0f, 10) << "\n";
}

Output:

1.0s of gravity, stepped different ways:
  1 step  of 1.00s : y = -10
  4 steps of 0.25s : y = -6.25
  10 steps of 0.10s: y = -5.5

One second of the same gravity lands the body at -10, or -6.25, or -5.5, purely depending on how the second was chopped up (smaller steps are closer to the true answer of -5.0, but that is not the point here). That spread is the non-determinism: hand physics whatever dt the frame happened to take, and the result is at the mercy of frame timing — impossible to reproduce for a replay, and guaranteed to drift two networked players apart. The fix is the accumulator: add each frame's real elapsed time to a running total, then spend that total in fixed-size chunks, carrying whatever is left over into the next frame.

#include <iostream>

int main() {
    const float FIXED = 0.25f;   // physics ALWAYS steps by this, never anything else
    float acc = 0.0f;
    float frameTimes[] = {0.30f, 0.10f, 0.28f, 0.40f};   // jittery real frame durations

    for (int f = 0; f < 4; f++) {
        acc += frameTimes[f];
        int stepsThisFrame = 0;
        while (acc >= FIXED) {
            acc -= FIXED;        // integrate(FIXED) would run here
            stepsThisFrame++;
        }
        std::cout << "  frame dt=" << frameTimes[f]
                  << " -> ran " << stepsThisFrame
                  << " fixed step(s), leftover=" << acc << "\n";
    }
}

Output:

  frame dt=0.3 -> ran 1 fixed step(s), leftover=0.05
  frame dt=0.1 -> ran 0 fixed step(s), leftover=0.15
  frame dt=0.28 -> ran 1 fixed step(s), leftover=0.18
  frame dt=0.4 -> ran 2 fixed step(s), leftover=0.08

The frame times jitter wildly — 0.30, 0.10, 0.28, 0.40 — but every physics step is exactly 0.25. A long frame (0.40) runs two steps to catch up; a short one (0.10) runs zero steps and just banks its time in the accumulator for later. Physics never once sees a variable dt, so it stays fully deterministic no matter how the frames stutter. This is exactly what Unity's FixedUpdate does: it is called zero, one, or several times per rendered frame so physics always advances in Time.fixedDeltaTime-sized chunks. The one wrinkle to know: the render frame usually lands between two fixed steps (that leftover sitting in the accumulator), so smooth engines interpolate the visual position between the last two physics states — otherwise a fast object visibly stutters even though the simulation underneath is perfect.

3. Broad phase: don't test everything against everything

With n bodies in the world, there are n * (n-1) / 2 possible pairs that might be touching. Testing every single pair with an exact shape test is wasteful: 100 bodies means 4,950 pairs checked every physics step, and almost all of them are nowhere near each other. Broad phase collision detection's only job is to throw out the obviously-not-touching pairs fast and cheap, using the AABBs from the geometry chapter, leaving a much smaller list of candidate pairs for an exact test.

Here are four circles, using the Circle/AABB structs from the geometry chapter, with a naive broad phase that checks every pair's AABB overlap:

#include <iostream>

struct Vec2 { float x, y; };
struct AABB { Vec2 min, max; };

AABB circleAABB(Vec2 center, float radius) {
    return { {center.x - radius, center.y - radius},
             {center.x + radius, center.y + radius} };
}

bool aabbOverlap(AABB a, AABB b) {
    if (a.max.x < b.min.x || a.min.x > b.max.x) return false;
    if (a.max.y < b.min.y || a.min.y > b.max.y) return false;
    return true;
}

int main() {
    struct Named { const char* name; Vec2 center; float radius; };
    Named bodies[] = {
        {"A", {0.0f, 0.0f}, 1.0f},
        {"B", {1.7f, 0.0f}, 1.0f},
        {"C", {5.0f, 5.0f}, 1.0f},
        {"D", {0.0f, 3.0f}, 1.0f},
    };

    int pairsChecked = 0, candidates = 0;
    for (int i = 0; i < 4; i++) {
        for (int j = i + 1; j < 4; j++) {
            AABB boxI = circleAABB(bodies[i].center, bodies[i].radius);
            AABB boxJ = circleAABB(bodies[j].center, bodies[j].radius);
            pairsChecked++;
            if (aabbOverlap(boxI, boxJ)) {
                std::cout << bodies[i].name << "-" << bodies[j].name << ": CANDIDATE\n";
                candidates++;
            }
        }
    }
    std::cout << "checked " << pairsChecked << " pairs, found " << candidates << " candidate(s)\n";
}

Output:

A-B: CANDIDATE
checked 6 pairs, found 1 candidate(s)

Six pairs checked, one survives. That one pair — A and B — is the only one narrow phase (section 4) needs to bother with. This naive version is still O(n^2) in how many pairs it checks (it just makes each check cheap), which is fine for a handful of objects but starts to add up once a scene has hundreds of bodies. Two common ways to avoid checking every pair at all:

Uniform grid: sort bodies into cells

Divide the world into fixed-size square cells and drop each body into whichever cell its center falls in (using floor(position / cellSize)). Only bodies sharing a cell (or a neighboring one, for objects near a cell edge) ever get compared — bodies in far-apart cells are never even considered:

#include <iostream>
#include <cmath>

struct Vec2 { float x, y; };

int cellOf(float coord, float cellSize) {
    return (int)std::floor(coord / cellSize);
}

int main() {
    struct Named { const char* name; Vec2 center; };
    Named bodies[] = {
        {"A", {0.0f, 0.0f}},
        {"B", {1.7f, 0.0f}},
        {"C", {5.0f, 5.0f}},
        {"D", {0.0f, 3.0f}},
    };
    float cellSize = 2.0f;

    for (auto& b : bodies) {
        int cx = cellOf(b.center.x, cellSize);
        int cy = cellOf(b.center.y, cellSize);
        std::cout << b.name << ": cell (" << cx << ", " << cy << ")\n";
    }
}

Output:

A: cell (0, 0)
B: cell (0, 0)
C: cell (2, 2)
D: cell (0, 1)

A and B land in the same cell; C and D each land alone. A real implementation groups bodies by cell (typically with a hash map keyed on the cell coordinates) and only tests pairs that share a cell or a neighboring one — the same one candidate pair the AABB sweep found, reached a different way. A grid is simple and fast when objects are roughly evenly spread out, but wastes memory on mostly-empty cells if objects cluster tightly in one area or spread across a huge, sparse world.

Sweep and prune: sort along one axis

Sweep and prune (sometimes called the "sort and sweep" method) takes a different angle: sort every body's AABB by its minimum x coordinate, then walk the sorted list left to right, keeping track of which AABBs are currently "open" (their span has started but not yet ended). Two AABBs can only possibly overlap if they are both open at the same time:

Engines like Box2D favor sweep and prune because, frame to frame, an object's neighbors rarely change — the sorted order from last frame is already almost correct, so re-sorting is nearly free (a handful of swaps, not a full sort) most steps. Whichever method a broad phase uses, its job is always the same: turn "all possible pairs" into "candidate pairs," fast, so the exact (and much more expensive) test in the next section only has to run a handful of times, not thousands.

4. Narrow phase: the exact test

Narrow phase collision detection takes the small list of candidate pairs broad phase produced and runs an exact, shape-specific test on each one: do these two shapes actually overlap, and if so, by how much? For two circles this is simple — the exact test the geometry chapter's point-in-circle idea generalizes directly to:

#include <iostream>
#include <cmath>

struct Vec2 { float x, y; };

Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
float length(Vec2 v) { return std::sqrt(v.x * v.x + v.y * v.y); }

bool circlesOverlap(Vec2 centerA, float rA, Vec2 centerB, float rB) {
    float dist = length(sub(centerB, centerA));
    return dist < (rA + rB);
}

int main() {
    Vec2 A = {0.0f, 0.0f}, B = {1.7f, 0.0f};
    std::cout << "A-B exact test: " << (circlesOverlap(A, 1.0f, B, 1.0f) ? "COLLIDING" : "clear") << "\n";
}

Output:

A-B exact test: COLLIDING

Broad phase already told us A and B were the only pair worth checking, so narrow phase only has one test to run — not the six a naive all-pairs check would have needed. "Two centers closer together than the sum of their radii" is the entire exact test for circles, which is why circles (and spheres in 3D) are the cheapest shape a physics engine can collide.

Circle versus box: clamp to the nearest point

Between "two circles" and "two polygons" sits a case common enough to earn its own shortcut: a circle against an axis-aligned box. The exact test is a tidy trick — find the point on the box closest to the circle's center by clamping the center into the box's range on each axis, then it is just a circle-versus-point check from there:

#include <iostream>
#include <cmath>

struct Vec2 { float x, y; };
struct AABB { Vec2 min, max; };

float clampf(float v, float lo, float hi) {
    return v < lo ? lo : (v > hi ? hi : v);
}

int main() {
    AABB box = { {0.0f, 0.0f}, {3.0f, 2.0f} };
    Vec2 center = {3.5f, 1.0f};
    float radius = 1.0f;

    // closest point on the box to the circle center = clamp center into the box
    Vec2 closest = { clampf(center.x, box.min.x, box.max.x),
                     clampf(center.y, box.min.y, box.max.y) };
    Vec2 d = { center.x - closest.x, center.y - closest.y };
    float dist = std::sqrt(d.x * d.x + d.y * d.y);

    std::cout << "closest point on box = (" << closest.x << ", " << closest.y << ")\n";
    std::cout << "distance to center   = " << dist << "\n";
    if (dist < radius) {
        Vec2 normal = { d.x / dist, d.y / dist };   // from box toward circle
        float penetration = radius - dist;
        std::cout << "COLLIDING  normal=(" << normal.x << ", " << normal.y
                  << ")  penetration=" << penetration << "\n";
    } else {
        std::cout << "clear\n";
    }
}

Output:

closest point on box = (3, 1)
distance to center   = 0.5
COLLIDING  normal=(1, 0)  penetration=0.5

The circle sits at x=3.5, just past the box's right edge at x=3. Clamping the center into the box's [0,3] x [0,2] range snaps it to (3, 1) — the nearest point on the box — which is 0.5 away, inside the radius of 1, so they collide by 0.5 with a normal pointing straight out along +x. The one case this simple version does not handle is the center being inside the box (then the clamp returns the center itself, dist is 0, and you cannot divide to get a normal); a full implementation detects that and pushes out along the nearest face instead. Engines keep hand-written shortcuts like this for circle-vs-circle, circle-vs-box, and box-vs-box precisely because they are this cheap.

Boxes and polygons: the Separating Axis Theorem

Circles are the easy case. For boxes and other convex polygons (a shape is convex if a straight line between any two points inside it never leaves the shape), the classic exact test is the Separating Axis Theorem (SAT): two convex shapes are not colliding if and only if you can find at least one axis where projecting both shapes onto it leaves a gap between them. Try every candidate axis (for polygons, that means every edge's perpendicular normal); if every single one shows overlap, the shapes are colliding, and the axis with the smallest overlap tells you the collision normal.

SAT in code: projecting two boxes onto an axis

Here is that theorem as a runnable test. For each candidate axis (every edge normal of both shapes), project both polygons onto it — reduce each to the [min, max] shadow it casts on that axis — and check whether the shadows overlap. The moment one axis shows a gap, stop: the shapes are apart. If every axis overlaps, they collide, and the axis with the smallest overlap is the collision normal, its overlap the penetration depth:

#include <iostream>
#include <cmath>
#include <vector>

struct Vec2 { float x, y; };
float dot(Vec2 a, Vec2 b) { return a.x * b.x + a.y * b.y; }

// project every vertex of `poly` onto `axis`, return the [min,max] shadow
void projectOnto(const std::vector<Vec2>& poly, Vec2 axis, float& mn, float& mx) {
    mn = mx = dot(poly[0], axis);
    for (const Vec2& v : poly) {
        float p = dot(v, axis);
        if (p < mn) mn = p;
        if (p > mx) mx = p;
    }
}

bool satOverlap(const std::vector<Vec2>& A, const std::vector<Vec2>& B,
                Vec2& normal, float& depth) {
    depth = 1e30f;
    for (int poly = 0; poly < 2; poly++) {               // axes from BOTH shapes
        const std::vector<Vec2>& P = (poly == 0) ? A : B;
        for (size_t i = 0; i < P.size(); i++) {
            Vec2 a = P[i], b = P[(i + 1) % P.size()];
            Vec2 edge = { b.x - a.x, b.y - a.y };
            Vec2 axis = { -edge.y, edge.x };             // edge normal (perpendicular)
            float len = std::sqrt(axis.x * axis.x + axis.y * axis.y);
            axis = { axis.x / len, axis.y / len };
            float minA, maxA, minB, maxB;
            projectOnto(A, axis, minA, maxA);
            projectOnto(B, axis, minB, maxB);
            float overlap = std::min(maxA, maxB) - std::max(minA, minB);
            if (overlap <= 0.0f) return false;           // a gap -> separated, done
            if (overlap < depth) { depth = overlap; normal = axis; }
        }
    }
    return true;
}

void test(const char* name, std::vector<Vec2> A, std::vector<Vec2> B) {
    Vec2 n; float d;
    if (satOverlap(A, B, n, d))
        std::cout << name << ": COLLIDING  normal=(" << n.x << ", " << n.y
                  << ")  depth=" << d << "\n";
    else
        std::cout << name << ": separated (a gap axis was found)\n";
}

int main() {
    std::vector<Vec2> A = {{0,0},{2,0},{2,2},{0,2}};                 // axis-aligned square
    std::vector<Vec2> Bhit  = {{2.8f,0},{3.8f,1},{2.8f,2},{1.8f,1}}; // diamond, overlaps A
    std::vector<Vec2> Bmiss = {{4.2f,0},{5.2f,1},{4.2f,2},{3.2f,1}}; // diamond, clears A
    test("overlapping", A, Bhit);
    test("separated  ", A, Bmiss);
}

Output:

overlapping: COLLIDING  normal=(-1, 0)  depth=0.2
separated  : separated (a gap axis was found)

The second box is a diamond (a square turned 45 degrees), so this is genuine SAT, not a disguised AABB check — projecting onto the diamond's slanted edge normals is real work. In the overlapping case, square A spans x in [0, 2] and the diamond's leftmost tip sits at x=1.8, so their shadows on the x-axis overlap by just 0.2 — smaller than the overlap on any other axis tried — so that axis wins as the collision normal with a penetration of 0.2. (The sign of the normal is whichever edge produced it; a real engine flips it to always point from A toward B, so here it would report (1, 0).) In the separated case the diamond has been slid right until its tip clears A's edge, and the very first axis showing a gap makes satOverlap return early — that early-out is what makes SAT fast in the common not-touching case.

General convex shapes: GJK

SAT works cleanly for polygons because "try every edge normal" is a short, known list of axes. For a general convex shape — a rounded capsule, an arbitrary convex mesh — there is no small fixed list of axes to try. The standard answer is the GJK algorithm (Gilbert-Johnson-Keerthi, named after its three authors), which sidesteps the axis-search problem entirely. GJK works on the Minkowski difference of the two shapes (informally: "shape A minus shape B," built by subtracting every point of B from every point of A) — a beautiful fact about this construction is that the two original shapes overlap if and only if this combined shape contains the origin point (0,0).

GJK never has to build the whole Minkowski difference shape (which can be huge). Instead it uses a support function — "give me the point of this shape furthest in direction D" — to grow a small triangle (a simplex) step by step, each time picking a new direction that moves the simplex closer to enclosing the origin. If the simplex ever encloses the origin, the shapes overlap; if the search direction ever points away from the origin with no way to get closer, they do not. The full algorithm is genuinely intricate to implement correctly, which is exactly why almost nobody hand-writes it — Unity's PhysX backend, Unreal's Chaos physics, and virtually every general-purpose physics engine ships a battle-tested GJK (often paired with a related algorithm, EPA, to recover penetration depth once GJK confirms overlap). What matters for you is recognizing the name and knowing what problem it solves: the exact overlap test for general convex shapes, when a short list of axes like SAT's is not available.

5. Building a contact: point, normal, penetration depth

A yes/no "colliding" answer is not enough to actually resolve a collision — the engine needs to know exactly where the shapes touch and which way to push them apart. Narrow phase's real output, once it confirms an overlap, is a small struct called a contact (or manifold when a single pair produces more than one contact point, like a box resting flat on a floor):

• Point — roughly where, in world space, the two shapes are touching.
• Normal — a unit vector (length 1) pointing from body A's surface toward body B, perpendicular to the surfaces at the contact.
• Penetration depth — how far the shapes currently overlap along that normal; a positive number means they are stuck inside each other by that much.

For two circles, every one of those three pieces falls straight out of the exact test from section 4 — the center-to-center line already is the normal direction:

#include <iostream>
#include <cmath>

struct Vec2 { float x, y; };

Vec2 add(Vec2 a, Vec2 b) { return { a.x + b.x, a.y + b.y }; }
Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
Vec2 scale(Vec2 a, float s) { return { a.x * s, a.y * s }; }
float length(Vec2 v) { return std::sqrt(v.x * v.x + v.y * v.y); }

struct Contact {
    Vec2 point;
    Vec2 normal;       // points from A toward B
    float penetration;
};

Contact makeCircleContact(Vec2 posA, float rA, Vec2 posB, float rB) {
    Vec2 delta = sub(posB, posA);
    float dist = length(delta);
    Vec2 normal = scale(delta, 1.0f / dist);

    Vec2 pointOnA = add(posA, scale(normal, rA));
    Vec2 pointOnB = sub(posB, scale(normal, rB));
    Vec2 midpoint = scale(add(pointOnA, pointOnB), 0.5f);

    float penetration = (rA + rB) - dist;
    return { midpoint, normal, penetration };
}

int main() {
    Contact c = makeCircleContact({0.0f, 0.0f}, 1.0f, {1.7f, 0.0f}, 1.0f);
    std::cout << "point=(" << c.point.x << ", " << c.point.y << ")\n";
    std::cout << "normal=(" << c.normal.x << ", " << c.normal.y << ")\n";
    std::cout << "penetration=" << c.penetration << "\n";
}

Output:

point=(0.85, 0)
normal=(1, 0)
penetration=0.3

A (center at the origin) and B (center 1.7 units to the right) overlap by 0.3 units, the contact sits at the midpoint of that overlap, and the normal points straight from A toward B. Every remaining section in this chapter works entirely off this one small struct — it is the hand-off point between "did they collide" (sections 1-4) and "what do we do about it" (sections 6 onward).

6. Resolving penetration: positional correction

Left alone, two overlapping bodies would just sit stuck inside each other. Positional correction directly moves the two bodies apart along the contact normal, enough to remove the penetration. Which body moves how much follows the same invMass logic from section 1: a body with more inverse mass (a lighter, easier-to-push object) gets pushed further; a static body (invMass = 0) never moves at all, and the other body absorbs the entire correction:

#include <iostream>

struct Vec2 { float x, y; };
struct RigidBody { Vec2 position; float invMass; };
struct Contact { Vec2 point; Vec2 normal; float penetration; };

Vec2 add(Vec2 a, Vec2 b) { return { a.x + b.x, a.y + b.y }; }
Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
Vec2 scale(Vec2 a, float s) { return { a.x * s, a.y * s }; }

void correctPenetration(RigidBody& a, RigidBody& b, Contact c) {
    float totalInvMass = a.invMass + b.invMass;
    if (totalInvMass <= 0.0f) return;   // both static -- nothing can move

    float shareA = a.invMass / totalInvMass;
    float shareB = b.invMass / totalInvMass;

    Vec2 correction = scale(c.normal, c.penetration);
    a.position = sub(a.position, scale(correction, shareA));
    b.position = add(b.position, scale(correction, shareB));
}

int main() {
    RigidBody a = { {0.0f, 0.0f}, 0.5f };   // mass 2 -> invMass 0.5
    RigidBody b = { {1.7f, 0.0f}, 1.0f };   // mass 1 -> invMass 1.0
    Contact c = { {0.85f, 0.0f}, {1.0f, 0.0f}, 0.3f };

    correctPenetration(a, b, c);

    std::cout << "a.position=(" << a.position.x << ", " << a.position.y << ")\n";
    std::cout << "b.position=(" << b.position.x << ", " << b.position.y << ")\n";
}

Output:

a.position=(-0.1, 0)
b.position=(1.9, 0)

Body A (twice as heavy, half the inverse mass) only moves 0.1 units left; body B (lighter, full inverse mass) moves 0.2 units right — together they cover the full 0.3 penetration, split in proportion to how "easy" each one is to push. Check the new distance between them: 1.9 - (-0.1) = 2.0, exactly the sum of their radii — the overlap is completely gone.

Baumgarte: push out through the velocity solver instead

Editing positions directly, as above, works but fights the velocity solver coming up in section 7 — you move a body one way, then the impulse step moves it another. Most engines instead fold penetration correction into the velocity solve using Baumgarte stabilization: turn the leftover overlap into a small extra separating velocity, a bias, and let the normal impulse deliver it. The bias is (beta / dt) * (penetration - slop), where beta is a small fraction (typically 0.1 to 0.2) setting how aggressively to close the gap:

#include <iostream>

int main() {
    float penetration = 0.30f;   // how deep the bodies overlap
    float slop        = 0.01f;   // permitted overlap left uncorrected
    float beta        = 0.20f;   // Baumgarte factor (fraction closed per step)
    float dt          = 1.0f / 60.0f;

    float correctable = penetration - slop;
    if (correctable < 0.0f) correctable = 0.0f;
    float biasVel = (beta / dt) * correctable;   // added to the target separating speed

    float pushThisStep = biasVel * dt;           // how far it moves them apart this step
    std::cout << "penetration      = " << penetration << "\n";
    std::cout << "bias velocity    = " << biasVel << " units/s\n";
    std::cout << "push-out this dt = " << pushThisStep << "\n";
    std::cout << "= beta*(pen-slop)= " << beta * correctable << "  (same number)\n";
}

Output:

penetration      = 0.3
bias velocity    = 3.48 units/s
push-out this dt = 0.058
= beta*(pen-slop)= 0.058  (same number)

The punchline is in the last two lines: pushing the bodies apart with a bias velocity of 3.48 for one step moves them 0.058 apart — which is exactly beta * (penetration - slop) = 0.2 * 0.29. The /dt inside the bias and the *dt of one integration step cancel, so beta is simply "what fraction of the overlap to remove per step" — the same 20%-per-step idea as the split positional correction above, just delivered as a velocity the existing impulse solver already knows how to apply. That is why it is popular: no separate position pass, one solver does both jobs.

7. Resolving velocity: the impulse and restitution

Fixing position stops bodies from overlapping, but it says nothing about how they should be moving after the hit — a ball should bounce, not just teleport out of the floor and keep falling through it next step. That is the job of an impulse: an instantaneous change in velocity, applied directly (not gradually, the way a force applied over time works), specifically along the contact normal.

The key ingredient is relative velocity along the normal — how fast the two bodies are approaching each other, measured along the direction they are about to separate on. If that number is already zero or positive, they are separating (or just touching) on their own and no impulse is needed. If it is negative, they are approaching, and the impulse needs to cancel that approach and, depending on restitution (the bounciness coefficient, from 0 = objects stick, "totally inelastic," to 1 = a perfectly elastic bounce with no energy lost), send them apart again:

The impulse magnitude formula (a standard result, derived from requiring momentum to be conserved and the post-collision separating speed to equal restitution * approaching speed) is:

#include <iostream>

struct Vec2 { float x, y; };
struct RigidBody { Vec2 velocity; float invMass; };
struct Contact { Vec2 normal; };

Vec2 add(Vec2 a, Vec2 b) { return { a.x + b.x, a.y + b.y }; }
Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
Vec2 scale(Vec2 a, float s) { return { a.x * s, a.y * s }; }
float dot(Vec2 a, Vec2 b) { return a.x * b.x + a.y * b.y; }

float resolveVelocity(RigidBody& a, RigidBody& b, Contact c, float restitution) {
    Vec2 relVel = sub(b.velocity, a.velocity);
    float velAlongNormal = dot(relVel, c.normal);

    if (velAlongNormal >= 0.0f) return 0.0f;   // already separating -- nothing to do

    float totalInvMass = a.invMass + b.invMass;
    float j = -(1.0f + restitution) * velAlongNormal / totalInvMass;

    Vec2 impulse = scale(c.normal, j);
    a.velocity = sub(a.velocity, scale(impulse, a.invMass));
    b.velocity = add(b.velocity, scale(impulse, b.invMass));
    return j;
}

int main() {
    RigidBody a = { {4.0f, 0.0f}, 0.5f };    // mass 2
    RigidBody b = { {-2.0f, 0.0f}, 1.0f };   // mass 1
    Contact c = { {1.0f, 0.0f} };

    float momentumBefore = (1.0f / a.invMass) * a.velocity.x + (1.0f / b.invMass) * b.velocity.x;
    float j = resolveVelocity(a, b, c, 0.5f);
    float momentumAfter = (1.0f / a.invMass) * a.velocity.x + (1.0f / b.invMass) * b.velocity.x;

    std::cout << "impulse j=" << j << "\n";
    std::cout << "a.velocity=(" << a.velocity.x << ", " << a.velocity.y << ")\n";
    std::cout << "b.velocity=(" << b.velocity.x << ", " << b.velocity.y << ")\n";
    std::cout << "momentum before=" << momentumBefore << " after=" << momentumAfter << "\n";
}

Output:

impulse j=6
a.velocity=(1, 0)
b.velocity=(4, 0)
momentum before=6 after=6

Body A (mass 2) was moving right at 4, body B (mass 1) was moving left at 2 — they were closing the gap at 6 units per second. After the impulse, A slows to 1 and B speeds up to 4: they now separate at 3 units per second, exactly half the approach speed, matching the restitution of 0.5. Momentum (mass * velocity, summed over both bodies) reads 6 both before and after — a real physical law falling out of the formula for free, not something the code enforces directly, which is a strong sign the formula is correct.

8. Friction: impulses along the surface

Section 7's impulse only touches velocity along the normal — it says nothing about sliding sideways. Friction is a second impulse, applied along the surface (the tangent direction, perpendicular to the normal), that resists that sideways sliding. Real friction follows an approximate rule called Coulomb friction: the maximum sideways force available is proportional to how hard the surfaces are being pressed together, capped by a friction coefficient (commonly written mu, the Greek letter mu). Concretely: compute the impulse that would fully cancel the sideways relative velocity, then clamp its size to mu * normalImpulse — if the full stop would need more than that, the surfaces slide (kinetic friction); if not, they grip completely:

#include <iostream>
#include <cmath>

struct Vec2 { float x, y; };
struct RigidBody { Vec2 velocity; float invMass; };
struct Contact { Vec2 normal; };

Vec2 add(Vec2 a, Vec2 b) { return { a.x + b.x, a.y + b.y }; }
Vec2 sub(Vec2 a, Vec2 b) { return { a.x - b.x, a.y - b.y }; }
Vec2 scale(Vec2 a, float s) { return { a.x * s, a.y * s }; }
float dot(Vec2 a, Vec2 b) { return a.x * b.x + a.y * b.y; }
float length(Vec2 v) { return std::sqrt(v.x * v.x + v.y * v.y); }

void applyFriction(RigidBody& a, RigidBody& b, Contact c, float normalImpulse, float friction) {
    Vec2 relVel = sub(b.velocity, a.velocity);

    // build the tangent: the part of relVel that is NOT along the normal
    float velAlongNormal = dot(relVel, c.normal);
    Vec2 tangentRaw = sub(relVel, scale(c.normal, velAlongNormal));
    float tangentLen = length(tangentRaw);
    if (tangentLen < 0.0001f) return;   // no sideways motion to resist
    Vec2 tangent = scale(tangentRaw, 1.0f / tangentLen);

    float totalInvMass = a.invMass + b.invMass;
    float velAlongTangent = dot(relVel, tangent);
    float jt = -velAlongTangent / totalInvMass;   // impulse to fully stop sliding

    float maxFriction = friction * normalImpulse;   // Coulomb's clamp
    if (jt > maxFriction) jt = maxFriction;
    if (jt < -maxFriction) jt = -maxFriction;

    Vec2 frictionImpulse = scale(tangent, jt);
    a.velocity = sub(a.velocity, scale(frictionImpulse, a.invMass));
    b.velocity = add(b.velocity, scale(frictionImpulse, b.invMass));
}

int main() {
    RigidBody floorBody = { {0.0f, 0.0f}, 0.0f };    // static
    RigidBody box       = { {3.0f, -2.0f}, 1.0f };   // sliding in and down
    Contact c = { {0.0f, 1.0f} };                    // normal points up

    // step 1: resolve the normal impulse (section 7's formula, restitution 0.2)
    Vec2 relVel = sub(box.velocity, floorBody.velocity);
    float velAlongNormal = dot(relVel, c.normal);
    float totalInvMass = floorBody.invMass + box.invMass;
    float j = -(1.0f + 0.2f) * velAlongNormal / totalInvMass;
    box.velocity = add(box.velocity, scale(c.normal, j * box.invMass));
    std::cout << "after normal impulse: box.velocity=(" << box.velocity.x << ", " << box.velocity.y << ")\n";

    // step 2: friction, clamped by mu = 0.5
    applyFriction(floorBody, box, c, j, 0.5f);
    std::cout << "after friction:       box.velocity=(" << box.velocity.x << ", " << box.velocity.y << ")\n";
}

Output:

after normal impulse: box.velocity=(3, 0.4)
after friction:       box.velocity=(1.8, 0.4)

The box hit the floor moving right at 3 and down at 2. The normal impulse (restitution 0.2) kills the downward motion and gives a small upward bounce of 0.4, but leaves the full sideways speed of 3 untouched. Friction then tries to fully cancel that sideways speed — that would need an impulse of magnitude 3 — but the Coulomb clamp only allows 0.5 * 2.4 = 1.2, so the box keeps sliding, just slower: 3 drops to 1.8. If the friction coefficient had been high enough to cover the full 3, the box would have stopped sliding completely in one step instead — the clamp is exactly the boundary between "grips" and "slides."

9. Stacking and why solvers iterate

One contact resolves cleanly in a single pass, as sections 6 and 7 just showed. Real scenes rarely have just one contact — a crate resting on the floor with another crate resting on top of it has two contacts sharing one body (the middle crate), and fixing one can quietly break the other. Watch it happen: a static floor, box2 resting on it (sunk in by 0.3), box3 resting on box2 (also sunk in by 0.3), each box exactly 1 unit tall:

#include <iostream>

struct Box { float bottom; float invMass; };   // 1D vertical stack; height is always 1

void resolveContact(Box& lower, Box& upper) {
    float lowerTop = lower.bottom + 1.0f;
    float penetration = lowerTop - upper.bottom;
    if (penetration <= 0.0f) return;

    float totalInv = lower.invMass + upper.invMass;
    if (totalInv <= 0.0f) return;

    float shareLower = lower.invMass / totalInv;
    float shareUpper = upper.invMass / totalInv;
    lower.bottom -= penetration * shareLower;
    upper.bottom += penetration * shareUpper;
}

int main() {
    Box floorBox = { 0.0f, 0.0f };   // static: top always at 1.0
    Box box2     = { 0.7f, 1.0f };   // sunk 0.3 into the floor
    Box box3     = { 1.4f, 1.0f };   // sunk 0.3 into box2

    for (int iter = 1; iter <= 3; iter++) {
        resolveContact(floorBox, box2);   // contact A: floor - box2
        resolveContact(box2, box3);       // contact B: box2 - box3

        float penA = (floorBox.bottom + 1.0f) - box2.bottom;
        float penB = (box2.bottom + 1.0f) - box3.bottom;
        std::cout << "iteration " << iter << ": penetrationA=" << penA
                   << " penetrationB=" << penB << "\n";
    }
}

Output:

iteration 1: penetrationA=0.3 penetrationB=0
iteration 2: penetrationA=0.15 penetrationB=0
iteration 3: penetrationA=0.075 penetrationB=0

Trace what happens in iteration 1: fixing contact A (floor-box2) pushes box2 fully up by 0.3, which shoves box2 deeper into box3 — the box2-box3 penetration jumps from 0.3 to 0.6. Fixing contact B then splits that 0.6 gap 50/50 between box2 and box3, which pulls box2 back down by 0.3 — right back to where it started, re-opening contact A. After one full pass, contact A is exactly as broken as when the step began. Only by processing both contacts again (iteration 2) does contact A improve — the penetration halves, to 0.15, then halves again to 0.075 in iteration 3, and would keep halving forever without ever quite reaching zero in a finite number of steps.

This is why physics engines are called iterative solvers: a single pass through every contact (this pattern is called sequential impulse when applied to velocities, and is a form of the Gauss-Seidel method mathematically) rarely satisfies every contact simultaneously the moment two or more contacts share a body. Running several passes — a handful of iterations, commonly somewhere around 4 to 10 — converges the whole system toward a state that is close enough to "every contact satisfied at once" to look and feel solid. More iterations mean a stiffer, more accurate stack at more CPU cost per step; fewer iterations are cheaper but let stacks visibly sink, jitter, or feel slightly springy. This is a real, tunable trade-off, not a bug to be fixed once and forgotten — which is exactly why engines like Unity expose the iteration count as a setting rather than hard-coding it.

10. Continuous collision detection: stopping tunneling

Every test so far is discrete: it looks at where a body is at the start and end of a step and checks those two snapshots for overlap. That works fine for anything moving slower than its own size per step. It fails completely for anything fast and thin — a bullet, a thrown knife, a speeding car versus a chain-link fence — because the object can be entirely on one side of a thin wall at the start of a step and entirely on the other side at the end, having passed straight through the middle without either snapshot ever catching it mid-flight. This failure is called tunneling.

#include <iostream>

struct Wall { float xStart, xEnd; };

bool discreteHitsWall(float x, Wall w) {
    return x >= w.xStart && x <= w.xEnd;
}

int main() {
    Wall wall = { 5.0f, 5.1f };   // a thin wall, only 0.1 units thick
    float x0 = 4.5f;              // position at the START of this step
    float speed = 80.0f;
    float dt = 0.02f;
    float x1 = x0 + speed * dt;   // position at the END of this step

    std::cout << "x0=" << x0 << " x1=" << x1 << "\n";
    std::cout << "discrete check at x0: " << (discreteHitsWall(x0, wall) ? "HIT" : "clear") << "\n";
    std::cout << "discrete check at x1: " << (discreteHitsWall(x1, wall) ? "HIT" : "clear") << "\n";
}

Output:

x0=4.5 x1=6.1
discrete check at x0: clear
discrete check at x1: clear

The ball travels 1.6 units this step (80 * 0.02) but the wall is only 0.1 units thick — the ball leaps clean over the wall's entire span between one sample and the next, and both discrete checks report "clear." The ball tunneled straight through and nothing noticed.

Continuous collision detection (CCD) fixes this by testing the entire path the object swept through this step, not just its two endpoints — a swept test. The simplest version treats the movement as a line segment from the start position to the end position and finds exactly where that segment first crosses the wall:

#include <iostream>

struct Wall { float xStart, xEnd; };

bool sweptHitsWall(float x0, float x1, Wall w, float& hitT, float& hitX) {
    if (x0 >= w.xStart) return false;   // already at/past the wall's near face
    if (x1 < w.xStart) return false;    // never reached the wall at all

    hitT = (w.xStart - x0) / (x1 - x0);
    hitX = x0 + hitT * (x1 - x0);
    return true;
}

int main() {
    Wall wall = { 5.0f, 5.1f };
    float x0 = 4.5f, x1 = 6.1f;

    float hitT, hitX;
    if (sweptHitsWall(x0, x1, wall, hitT, hitX)) {
        std::cout << "swept test: HIT at t=" << hitT << ", x=" << hitX << "\n";
    } else {
        std::cout << "swept test: clear\n";
    }
}

Output:

swept test: HIT at t=0.3125, x=5

The swept test catches exactly what the two discrete snapshots missed: the ball crosses the wall's near face 31.25% of the way through this step's movement, at position x = 5.0. A real engine responds by stopping the object right there (or resolving a full collision at that point) instead of letting it continue to x = 6.1. Real CCD implementations go further than this simple line test — conservative advancement expands the swept segment by the moving shape's own radius so a fast sphere, not just a fast point, is caught correctly, and some engines instead use speculative contacts, which predict an upcoming collision one step early and start resolving it slightly before the shapes actually touch. The core idea stays the same either way: look at the whole path swept this step, not just its two ends.

11. Resting contacts: jitter, restitution, and sleeping

A box dropped on the floor should, after a bounce or two, just sit there — dead still, costing nothing. Getting a body to truly come to rest is its own small problem, because the naive loop from the sections above never quite lets it. Gravity pulls the resting box down a hair every step, the contact catches it, and section 7's restitution faithfully bounces that hair back up. The result is a body that never stops twitching:

#include <iostream>

// one step of a ball above a floor at y=0: semi-implicit Euler, then a bounce.
// `threshold`: approach speeds slower than this get restitution 0 (no bounce).
void step(float& y, float& v, float g, float dt, float e, float threshold) {
    v += g * dt;
    y += v * dt;
    if (y < 0.0f) {                 // hit the floor this step
        y = 0.0f;
        if (v < 0.0f) {
            float useE = (-v < threshold) ? 0.0f : e;   // slow contact -> no bounce
            v = -useE * v;
        }
    }
}

void run(const char* label, float threshold) {
    float y = 0.0f, v = 0.0f, g = -10.0f, dt = 1.0f / 60.0f, e = 0.5f;
    std::cout << label << "\n";
    for (int i = 1; i <= 6; i++) {
        step(y, v, g, dt, e, threshold);
        std::cout << "  step " << i << ": y=" << y << "  v=" << v << "\n";
    }
}

int main() {
    run("no threshold (restitution always on) -> never settles:", 0.0f);
    run("with threshold 1.0 -> settles to rest:", 1.0f);
}

Output:

no threshold (restitution always on) -> never settles:
  step 1: y=0  v=0.0833333
  step 2: y=0  v=0.0416667
  step 3: y=0  v=0.0625
  step 4: y=0  v=0.0520833
  step 5: y=0  v=0.0572917
  step 6: y=0  v=0.0546875
with threshold 1.0 -> settles to rest:
  step 1: y=0  v=0
  step 2: y=0  v=0
  step 3: y=0  v=0
  step 4: y=0  v=0
  step 5: y=0  v=0
  step 6: y=0  v=0

With restitution always on (top), the velocity never reaches zero — it locks into a permanent buzz around +0.05, because every step gravity feeds in a little downward speed and restitution kicks it straight back up. In a real scene that residual upward velocity lifts the box a hair off the floor each step: the visible resting jitter of a stack that will not hold still. The fix (bottom) is a restitution threshold (also called a bounce threshold): if the approach speed is below some small value, treat restitution as zero, so a slow contact just stops instead of bouncing. One line, and the box settles on the very first step and stays put. Every engine does this; Unity exposes it as Physics.bounceThreshold.

Sleeping: switching a settled body off

The threshold gets a body to rest, but a resting body still costs CPU: every step it is integrated, tested for collisions, and run through the solver, only to end up exactly where it was. In a scene with hundreds of settled objects — a warehouse of crates, a pile of rubble — that is enormous waste, and worse, the tiny numerical noise left in the solver keeps nudging those bodies so they never sit perfectly still. The answer is sleeping: once a body's speed has stayed below a threshold for long enough, mark it asleep and stop simulating it entirely until something touches it and wakes it up.

#include <iostream>
#include <cmath>

struct Body {
    float x, v;
    float idleTime;   // how long we have been below the sleep threshold
    bool  asleep;
};

int main() {
    Body b = {0.0f, 5.0f, 0.0f, false};
    float dt = 1.0f / 60.0f;
    float sleepSpeed = 0.10f;   // below this counts as "still"
    float sleepDelay = 0.05f;   // must stay still this long to fall asleep

    for (int i = 1; i <= 14; i++) {
        if (!b.asleep) {
            b.v *= 0.6f;                 // friction bleeds speed off each step
            b.x += b.v * dt;
            if (std::fabs(b.v) < sleepSpeed) b.idleTime += dt;
            else                             b.idleTime = 0.0f;
            if (b.idleTime >= sleepDelay) b.asleep = true;   // go to sleep
        }
        // an asleep body skips integration entirely -- zero CPU until woken
        std::cout << "step " << i << ": x=" << b.x << "  v=" << b.v
                  << (b.asleep ? "  [ASLEEP]" : "") << "\n";
    }
}

Output:

step 1: x=0.05  v=3
step 2: x=0.08  v=1.8
step 3: x=0.098  v=1.08
step 4: x=0.1088  v=0.648
step 5: x=0.11528  v=0.3888
step 6: x=0.119168  v=0.23328
step 7: x=0.121501  v=0.139968
step 8: x=0.1229  v=0.0839808
step 9: x=0.12374  v=0.0503885
step 10: x=0.124244  v=0.0302331  [ASLEEP]
step 11: x=0.124244  v=0.0302331  [ASLEEP]
step 12: x=0.124244  v=0.0302331  [ASLEEP]
step 13: x=0.124244  v=0.0302331  [ASLEEP]
step 14: x=0.124244  v=0.0302331  [ASLEEP]

The body slows under friction until, at step 10, its speed has been under 0.1 for the required 0.05 seconds; it falls asleep, and from then on its position is frozen and the loop does no work on it at all. A real engine goes one step further and sleeps islands together — a stack of boxes all resting on each other must sleep and wake as a group, since waking the bottom one has to wake the ones above it. This is also a famous source of "bugs that are not bugs": a crate sitting on a pressure plate that a designer expects to trigger, or an object that should react to a slow nudge, does nothing — because it is asleep. In Unity the knobs are Rigidbody.sleepThreshold, the query rb.IsSleeping(), and the escape hatch rb.WakeUp() to force a body awake when your own code has changed something the physics engine cannot see.

12. Using the engine's physics: Rigidbody and Collider in Unity

Every idea in sections 1 through 11 is exactly what Unity's built-in physics (backed by PhysX) is doing for you, every fixed step, without you writing a single line of broad phase, narrow phase, or impulse-resolution code. Recognizing the concept inside the component is the whole point of this section — you will spend far more time configuring Rigidbody and Collider than writing your own solver.

using UnityEngine;

public class BouncyBall : MonoBehaviour
{
    void Start()
    {
        Rigidbody rb = GetComponent<Rigidbody>();

        rb.mass = 2f;                  // section 1's RigidBody.mass
        rb.linearVelocity = new Vector3(3f, 0f, 0f);   // section 1's RigidBody.velocity
        rb.useGravity = true;          // section 2's integrate(), running every FixedUpdate

        // section 10: switch on CCD only for objects at real risk of tunneling
        rb.collisionDetectionMode = CollisionDetectionMode.ContinuousDynamic;
    }
}

Rigidbody is section 1's RigidBody struct, field for field: mass, linearVelocity, angularVelocity, position/rotation (the latter stored as a quaternion, exactly as the rotation chapter described). Every FixedUpdate tick — the same fixed timestep from the game loop chapter — Unity runs section 2's semi-implicit Euler step on it internally, using Time.fixedDeltaTime, for exactly the stability and determinism reasons section 2 explained.

A Collider (BoxCollider, SphereCollider, CapsuleCollider, MeshCollider for convex meshes) is the shape PhysX's broad phase and narrow phase actually operate on — PhysX runs its own broad phase over every collider's world-space bounds (sections 3's job) and its own narrow phase (SAT-like tests for simple shapes, GJK-family algorithms for convex meshes, exactly section 4) to find contacts, entirely off-screen from your code:

using UnityEngine;

public class ContactLogger : MonoBehaviour
{
    void OnCollisionEnter(Collision collision)
    {
        // section 5's Contact struct, handed to you directly by Unity
        foreach (ContactPoint contact in collision.contacts)
        {
            Debug.Log("point=" + contact.point + " normal=" + contact.normal
                     + " separation=" + contact.separation);   // separation < 0 means penetrating
        }
    }
}

Every field from section 5's Contact struct is right there in ContactPoint: point, normal, and separation (Unity's name for penetration depth, with the sign flipped — negative means overlapping). Restitution and friction — sections 7 and 8 — live on a Physic Material asset assigned to the collider, with a combine mode (average, minimum, maximum, or multiply) controlling how two different bodies' materials blend when they touch:

using UnityEngine;

public class BouncyMaterialSetup : MonoBehaviour
{
    void Start()
    {
        PhysicsMaterial mat = new PhysicsMaterial();
        mat.bounciness = 0.5f;                              // section 7's restitution
        mat.dynamicFriction = 0.5f;                          // section 8's friction (sliding)
        mat.staticFriction = 0.6f;                           // section 8's friction (starting from rest)
        mat.frictionCombine = PhysicsMaterialCombine.Average;

        GetComponent<Collider>().material = mat;
    }
}

Section 9's iteration count is a real, exposed setting too — Project Settings > Physics > Solver Iterations (position iterations) and Solver Velocity Iterations control exactly how many passes PhysX's internal solver runs per step, the same convergence-through-repetition idea from section 9's stacked-box trace, just running on every contact in the scene at once rather than two contacts by hand.

13. Glossary

• Rigid body — an object whose shape never deforms; tracked by position, rotation, velocity, angular velocity, and mass.
• Static body — a rigid body with infinite mass (invMass = 0) that never moves, like a floor or wall.
• Dynamic body — a rigid body that moves and responds to forces and collisions.
• Inverse mass (invMass) — 1/mass; used so "infinite mass" (never moves) can be represented as 0 instead of an undefined division.
• Semi-implicit (symplectic) Euler — integrating velocity first, then position from the new velocity; the stable default from the calculus chapter, applied here to a rigid body.
• Explicit (forward) Euler — integrating position from the old velocity before updating velocity; not symplectic, gains energy on oscillating systems, and never used for game dynamics.
• Verlet integration — an integrator that stores the previous position instead of a velocity and infers motion from the two; symplectic, and the basis of cloth, rope, and position-based dynamics.
• Fixed timestep — advancing physics by the same dt every step, for stability, determinism, and predictable tunneling risk.
• Accumulator — a running total of elapsed frame time spent in fixed-size physics steps, carrying the remainder forward, so physics always sees the same dt however the frame rate jitters.
• Broad phase — a cheap first pass (AABB checks, a uniform grid, or sweep and prune) that throws out obviously-not-touching pairs before an exact test runs.
• Narrow phase — the exact, shape-specific overlap test run only on the candidate pairs broad phase produced.
• Uniform grid — a broad-phase method that buckets bodies into fixed-size cells and only compares bodies sharing a cell.
• Sweep and prune — a broad-phase method that sorts AABBs along one axis and only compares bodies whose sorted spans overlap.
• Separating Axis Theorem (SAT) — an exact test for convex polygons: shapes are not colliding if any candidate axis shows a gap between their projections.
• Closest-point test — the circle-versus-box narrow-phase check: clamp the circle's center into the box to find the nearest point, then compare that distance to the radius.
• GJK algorithm — an exact overlap test for general convex shapes, based on searching the Minkowski difference for the origin.
• Minkowski difference — the shape formed by subtracting every point of one shape from every point of another; contains the origin exactly when the two shapes overlap.
• Support function — a function returning the furthest point of a shape in a given direction; the core operation GJK is built from.
• Contact — the output of narrow phase: a point, a normal, and a penetration depth describing how two shapes touch.
• Penetration depth — how far two overlapping shapes currently overlap, measured along the contact normal.
• Positional correction — directly moving two overlapping bodies apart along the contact normal to remove penetration.
• Slop — a small amount of penetration real engines deliberately allow uncorrected, to avoid jitter from over-aggressive correction.
• Baumgarte stabilization — removing penetration by adding a small bias velocity (beta/dt * (penetration - slop)) into the velocity solver rather than moving positions directly.
• Split impulse / position projection — applying the penetration-correction bias on a throwaway pseudo-velocity so it never adds real momentum, avoiding the bounciness of plain Baumgarte.
• Impulse — an instantaneous change in velocity, applied directly rather than accumulated gradually like a force.
• Restitution — the bounciness coefficient (0 = no bounce, 1 = perfectly elastic) controlling how much separating speed an impulse produces.
• Restitution (bounce) threshold — a minimum approach speed below which restitution is treated as zero, so slow contacts stop dead instead of buzzing.
• Resting jitter — the permanent twitch of a body that should be at rest, caused by restitution bouncing back the tiny per-step velocity gravity keeps adding.
• Sleeping — deactivating a body whose speed has stayed below a threshold, skipping its integration and solving until it is touched and woken.
• Coulomb friction — the model where the maximum available sideways (frictional) impulse is proportional to the normal impulse, via a friction coefficient.
• Tangent — the direction along a contact surface, perpendicular to the contact normal; where friction impulses act.
• Sequential impulse / Gauss-Seidel solver — resolving contacts one at a time, in passes, rather than all at once; the standard approach for real-time physics.
• Solver iteration — one full pass over every contact in a scene; more iterations converge a stack of contacts closer to a fully consistent state.
• Continuous collision detection (CCD) — testing the entire path an object swept through a step, not just its start and end position.
• Tunneling — a fast object passing completely through a thin shape because a discrete check never sampled it mid-flight.
• Discrete collision detection — checking only a body's position at the start and end of a step, the default and cheaper mode.
• Rigidbody / Collider (Unity) — Unity's components for a rigid body's physical state and its collision shape, backed internally by PhysX.
• Physic Material — a Unity asset holding a collider's restitution (bounciness) and friction values, plus how they combine with another collider's material.

14. Exercises

AABBs:
E: min=(1, 1)   max=(3, 3)
F: min=(2.2, 1) max=(4.2, 3)
G: min=(9, 9)   max=(11, 11)

E-F: overlap (x: 1<=4.2 and 3>=2.2; y: identical range) -> CANDIDATE
E-G: no overlap (E.max.x=3 < G.min.x=9)
F-G: no overlap (F.max.x=4.2 < G.min.x=9)

narrow phase E-F:
distance = 1.2, sum of radii = 2  ->  1.2 < 2, colliding
penetration = 2 - 1.2 = 0.8
normal = (1, 0)
point on E = (3, 2), point on F = (2.2, 2), contact = midpoint = (2.6, 2)

invMassA = 1/3, invMassB = 1
relVel = velB - velA = (-4,0) - (2,0) = (-6, 0)
velocityAlongNormal = -6

j = -(1 + 1) * (-6) / (1/3 + 1) = 12 / (4/3) = 9

impulse = (9, 0)
velA_new = (2,0) - (9,0)*(1/3) = (2-3, 0) = (-1, 0)
velB_new = (-4,0) + (9,0)*1     = (-4+9, 0) = (5, 0)

momentum before = 3*2 + 1*(-4) = 6 - 4 = 2
momentum after  = 3*(-1) + 1*5 = -3 + 5 = 2   -- matches

travel this step = 150 * 0.02 = 3.0
x0 = 9.0, x1 = 9.0 + 3.0 = 12.0

discrete check at x0=9.0: 10.0 <= 9.0 <= 10.2 ? NO
discrete check at x1=12.0: 10.0 <= 12.0 <= 10.2 ? NO
-> tunneling: the discrete check misses the wall completely

swept test:
hitT = (10.0 - 9.0) / (12.0 - 9.0) = 1.0 / 3.0 = 0.3333
hitX = 9.0 + 0.3333 * 3.0 = 10.0

closest point = clamp center (6,4) into [2,5] x [1,3]
              = (5, 3)   -- a corner of the box

d    = center - closest = (6-5, 4-3) = (1, 1)
dist = sqrt(1*1 + 1*1)  = sqrt(2) = 1.41421
1.41421 < 1.5  ->  COLLIDING

penetration = 1.5 - 1.41421 = 0.08579
normal      = d / dist = (0.7071, 0.7071)

That is the full pipeline a physics engine runs every fixed step: integrate every body's position and rotation forward with semi-implicit Euler, narrow a world full of possible pairs down to a handful of real candidates with broad phase, confirm and measure each one exactly with narrow phase, turn a confirmed overlap into a contact, and then resolve that contact twice over — once to fix position (pushing penetration out) and once to fix velocity (impulses, restitution, and friction) — repeating the whole resolve step a few times so that stacks and piles of contacts settle into something consistent instead of fighting each other, and finally letting any body that has come fully to rest sleep so a quiet scene costs almost nothing. Continuous collision detection is the one addition needed to stop that pipeline from missing anything fast and thin. Unity's Rigidbody and Collider run every one of these steps for you, which is exactly why understanding them here turns "physics feels broken" into a short, specific checklist instead of a guess.