Lesson 12: Build a Physics Engine (+ pymunk)

🏗️ What Is in a Physics Engine?

A physics engine is a small set of data types and one function that is called over and over:

Body: position, velocity, angle, spin, and the inverse mass and inertia (0 for static bodies).
Shape: the collision geometry a body carries (in our engine, a circle or an axis-aligned box).
Contact: one touching pair found this step, with its normal, depth, point and the impulses applied so far.
Joint: a rule that ties two bodies together (our pendulum rod).
World: owns all of these, plus gravity and damping, and advances them with step(dt).

The order inside step matters. Our engine uses the order below; each stage only needs the results of the stage before it.

graph LR A["1. Forces gravity, damping change velocities"] --> B["2. Detect find contacts at current positions"] B --> C["3. Solve impulses on contacts and joints, 8 passes"] C --> D["4. Move positions += velocity × dt"] D -->|next fixed step| A

Stages 1 and 4 are semi-implicit Euler, split around the solver: velocities are updated first, the solver corrects them so bodies don't push into each other, and only then do positions move with the corrected velocities. That is why the solver works on velocities: by the time anything moves, the contacts have already had their say.

⏱️ Fixed Steps and Damping

In Velocity, Acceleration & Timesteps you built the fixed-step accumulator. An engine depends on it: contact solving and joints are tuned for one step size, and a step that is sometimes 1/144 s and sometimes 1/20 s makes piles jitter and rods stretch. The loop from that lesson, with one addition:

accumulator += min(clock.tick(60) / 1000, MAX_STEPS * STEP)   # never more than 8 steps' worth
while accumulator >= STEP:
    world.step(STEP)
    accumulator -= STEP

The min(...) caps how much time one frame can add. Without it, one slow frame (a window drag, a hiccup) adds so much time that the next frame has to run many steps, which makes it slow too, and the game can fall further and further behind. Capping loses a little simulated time after a hiccup instead.

A fixed step also makes runs repeatable: the same starting state and the same inputs, stepped with the same STEP, produce the same result each time you run the same program on the same machine. With a variable dt, every run takes slightly different steps, and small differences grow.

Damping must use dt. "Multiply the velocity by 0.99 every step" sounds harmless, but the effect per second depends on how many steps you run: at 60 steps per second it keeps 0.99⁶⁰ ≈ 55% of the velocity after one second, at 120 steps per second only 0.99¹²⁰ ≈ 30%, and applied inside 10 sub-steps of a 60 FPS frame it keeps 0.99⁶⁰⁰ ≈ 0.2%. Express damping as a rate per second and use the exponential form from Interpolation & Easing:

lin = math.exp(-self.linear_damping * dt)        # linear_damping is "per second"
body.vel *= lin

Now one second of simulation keeps exactly exp(-linear_damping) of the velocity, however you slice the second. The lab's tests check this with steps of 1/60, 1/120 and 1/240 s.

🧮 Integrators Compared

An integrator turns "acceleration now" into "position and velocity one step later". A mass on a spring is the classic test, because the exact answer is known (it swings forever with constant energy) and weak integrators visibly fail. Here are four, all given the same function accel(x):

Explicit Euler moves with the old velocity, then updates the velocity. One force evaluation per step.
Semi-implicit Euler updates the velocity first, then moves with the new one. Also one evaluation; it is the "Euler" in Velocity, Acceleration & Timesteps and in most game engines.
Velocity Verlet moves using the start-of-step acceleration, then evaluates the acceleration again at the new position and updates the velocity with the average of the two.
RK4 (fourth-order Runge–Kutta) samples the slope four times across the step and combines them with weights 1, 2, 2, 1.

This complete program runs all four on the spring for 10 seconds at 30 steps per second and prints the error and the energy change:

"""Integrator Shoot-out: Advanced Lesson 12 "try it" lab (solution).

A mass on a spring (period 1 s) is simulated for 10 seconds with four
integrators at 30 steps per second. The exact answer is known, so we can
measure each method's error and how much energy it gains or loses.
No window: the results print as a table.
"""
import math

OMEGA = 2 * math.pi            # rad/s: the spring swings once per second
X0, V0 = 100.0, 0.0            # start 100 px from rest, not moving
DT = 1 / 30
SECONDS = 10


def spring_accel(x):
    """Hooke's law with mass 1: a = -omega^2 * x (px/s^2)."""
    return -OMEGA * OMEGA * x


def explicit_euler(x, v, dt, accel):
    """Moves with the OLD velocity. Gains energy on a spring."""
    a = accel(x)
    return x + v * dt, v + a * dt


def semi_implicit_euler(x, v, dt, accel):
    """Updates velocity first, then moves with the NEW velocity (what most game engines use)."""
    v = v + accel(x) * dt
    return x + v * dt, v


def velocity_verlet(x, v, dt, accel):
    """Uses the acceleration at the start AND at the end of the step."""
    a0 = accel(x)
    x = x + v * dt + 0.5 * a0 * dt * dt
    a1 = accel(x)                      # re-evaluate the force at the new position
    return x, v + 0.5 * (a0 + a1) * dt


def rk4(x, v, dt, accel):
    """Classic 4th-order Runge-Kutta: four samples of the slope, weighted 1-2-2-1."""
    k1x, k1v = v, accel(x)
    k2x, k2v = v + 0.5 * dt * k1v, accel(x + 0.5 * dt * k1x)
    k3x, k3v = v + 0.5 * dt * k2v, accel(x + 0.5 * dt * k2x)
    k4x, k4v = v + dt * k3v, accel(x + dt * k3x)
    x_new = x + dt / 6 * (k1x + 2 * k2x + 2 * k3x + k4x)
    v_new = v + dt / 6 * (k1v + 2 * k2v + 2 * k3v + k4v)
    return x_new, v_new


def energy(x, v):
    return 0.5 * v * v + 0.5 * OMEGA * OMEGA * x * x


def simulate(method, dt=DT, seconds=SECONDS):
    x, v = X0, V0
    for _ in range(round(seconds / dt)):
        x, v = method(x, v, dt, spring_accel)
    return x, v


def main():
    exact_x = X0 * math.cos(OMEGA * SECONDS)
    e0 = energy(X0, V0)
    print(f"Spring, period 1 s, {SECONDS} s at dt = 1/{round(1 / DT)} s. Exact x = {exact_x:.2f} px")
    print(f"{'method':<22}{'x after 10 s':>14}{'error (px)':>12}{'energy change':>15}")
    for method in (explicit_euler, semi_implicit_euler, velocity_verlet, rk4):
        x, v = simulate(method)
        change = (energy(x, v) - e0) / e0 * 100
        print(f"{method.__name__:<22}{x:>14.2f}{abs(x - exact_x):>12.3f}{change:>14.2f}%")


if __name__ == "__main__":
    main()

Its output, which you can reproduce exactly (the arithmetic has no randomness):

Spring, period 1 s, 10 s at dt = 1/30 s. Exact x = 100.00 px
method                  x after 10 s  error (px)  energy change
explicit_euler              39151.23   39051.231   39200205.62%
semi_implicit_euler            98.12       1.878         -2.38%
velocity_verlet                99.33       0.665         -0.01%
rk4                            99.98       0.018         -0.03%

Explicit Euler gains energy every step until the spring explodes. Semi-implicit Euler is off by a couple of pixels, but its energy only wobbles around the true value instead of drifting away, which is what keeps game physics from blowing up. Verlet is several times more accurate for one extra evaluation, and RK4 is far more accurate still, for four evaluations per step.

Step:

So why don't engines use RK4? Because contacts don't fit its assumptions. RK4 needs to evaluate the force at several points inside the step, but a collision is not a smooth force: it is an impulse that changes the velocity all at once when the solver runs. Semi-implicit Euler has exactly one "change velocity" moment per step, which is where the solver slots in. Verlet and RK4 are the better choice for smooth forces without contacts, such as orbits, springs, cloth or a projectile preview line.

✅ Growth Mindset: Measure, Don't Memorize

You may have read that one integrator is simply "the best". The table shows why that question has no single answer: it depends on the forces, the step size and what you can afford per step. You don't need to remember a ranking. You now have a short program that measures one; change the step, the force or the run time, and let the numbers answer.

🧱 Bodies, Shapes and Contacts

Our engine's bodies are dynamic circles, and static circles or axis-aligned boxes. That small set covers balls, pegs, floors, walls and ledges, and every pair of shapes that can meet has a test:

def collide(a, b):
    if isinstance(a.shape, Circle) and isinstance(b.shape, Circle):
        return circle_circle(a, b)
    if isinstance(a.shape, Circle) and isinstance(b.shape, Box):
        return circle_box(a, b)
    if isinstance(a.shape, Box) and isinstance(b.shape, Circle):
        return circle_box(b, a)                         # keep the circle as A
    return None                                         # box-box: both static here

A dispatcher like this should cover every combination your bodies can produce, even if the answer for some is "never collides". Dynamic boxes would need the SAT code from SAT & Rotational Collisions, so this engine refuses them up front: Body raises ValueError if you give a Box a mass.

Static bodies get inv_mass = inv_inertia = 0. Every formula already multiplies by the inverse mass, so a static body needs no special cases: impulses and pushes simply multiply to zero. (Setting only inv_mass and forgetting inv_inertia is a classic bug: the anchor of a pendulum then starts to spin.)

Circle against box clamps the circle's center into the box to find the closest point. If the center is outside the box, the normal points from the center to that point; if the center is already inside, the circle leaves through the nearest face:

def circle_box(a, b):
    """A is a circle, B an axis-aligned box. The normal points from the circle into the box."""
    lo, hi = b.pos - b.shape.half, b.pos + b.shape.half
    c, r = a.pos, a.shape.radius
    closest = Vector2(max(lo.x, min(c.x, hi.x)), max(lo.y, min(c.y, hi.y)))
    delta = closest - c
    if delta.length_squared() > 0:                      # center outside the box
        dist = delta.length()
        if dist >= r:
            return None
        return Contact(a, b, delta / dist, r - dist, closest)
    # center inside the box: leave through the nearest face
    faces = [(c.x - lo.x, Vector2(1, 0)), (hi.x - c.x, Vector2(-1, 0)),
             (c.y - lo.y, Vector2(0, 1)), (hi.y - c.y, Vector2(0, -1))]
    gap, normal = min(faces, key=lambda f: f[0])
    return Contact(a, b, normal, r + gap, c - normal * gap)

🔧 The Solver: Accumulated Impulses and Joints

In SAT & Rotational Collisions each solver pass applied a fresh impulse if the contact was still closing. The engine uses a refinement that production solvers rely on: each contact remembers the total impulse applied so far this step (c.jn), and every pass adjusts that total, clamped so it never becomes negative:

dj = (c.target - vn) / effective_mass(a, b, ra, rb, n)
new_jn = max(c.jn + dj, 0.0)                     # contacts push, never pull
apply_impulse(a, b, n * (new_jn - c.jn), ra, rb)
c.jn = new_jn

Why bother? When one contact's pass overshoots, a later pass can take some of that push back, as long as the total stays a push. With fresh impulses only, an overshoot can never be undone, and piles buzz. The friction total c.jt is clamped the same way, to ±μ times c.jn.

The target speed. Before the passes, each contact gets a target separating speed: the bounce -e * vn for real impacts (none below 30 px/s, so resting bodies stay still), or a small push-out speed BETA / dt * (depth - SLOP) when bodies overlap, whichever is larger. Feeding a little of the overlap back into the velocity each step (the Baumgarte method) gently separates sunken bodies without a separate position pass.

Spin makes balls roll. The effective mass includes the spin terms from last lesson:

def effective_mass(a, b, ra, rb, axis):
    """1 / (the mass the impulse 'feels' along axis), with the spin terms."""
    ran, rbn = cross(ra, axis), cross(rb, axis)
    return a.inv_mass + b.inv_mass + ran * ran * a.inv_inertia + rbn * rbn * b.inv_inertia

For a ball on the floor, the contact point is straight below the center, so along the normal the spin term is 0, but along the tangent it equals r²/I. That term is why friction turns sliding into rolling: the lab's test slides a ball at 300 px/s and checks that after one second its spin times its radius matches its speed, the condition for rolling without slipping. Leave the spin terms out and friction pushes on the wrong mass, so the ball's spin and speed never match up.

Joints are velocity constraints too. A pendulum rod says "the distance stays L". Correcting only the positions each step leaves the velocity still pointing outward, so the rod keeps fighting it and stretches like rubber. Our DistanceJoint instead removes the part of the relative velocity along the rod, plus a small Baumgarte bias for any length error:

def solve(self):
    if self.k == 0:
        return
    rel = (self.b.vel - self.a.vel).dot(self.axis)
    impulse = self.axis * ((self.bias - rel) / self.k)
    self.a.vel -= impulse * self.a.inv_mass
    self.b.vel += impulse * self.b.inv_mass

The anchor is a static body, so a.inv_mass is 0 and only the bob moves. In the lab the rod stays within a fraction of a pixel of its length, with balls bouncing off the bob.

✅ Growth Mindset: Engines Are Built One Test at a Time

Nobody writes a physics engine in one go, and when yours misbehaves it is rarely clear which part is at fault. That is not a sign you aren't ready for this; it is the nature of systems where everything affects everything. The way through is small scenes: one ball on a floor, one sliding ball, one pendulum. Each of the lab's tests is one of those scenes. When a big scene breaks, rebuild it one body at a time until it breaks again.

🐍 The Same Scene in pymunk

pymunk is a Python physics library built on a C engine derived from Chipmunk2D. Install it with pip install pymunk (this course was tested with pymunk 7.3). It uses the same ideas you just built, bodies, shapes, a space, constraints and a fixed-step solver, and adds everything a small engine leaves out.

Here is the mini engine's scene rebuilt in pymunk, plus a stack of four rotating boxes, which our engine can't do:

"""The Same Scene in pymunk: Advanced Lesson 12 comparison lab (solution).

The mini engine's scene rebuilt with pymunk (a Python physics library on
top of a C engine derived from Chipmunk2D): floor, walls, ledge, pegs, a
pendulum and bouncing balls, plus a stack of rotating boxes that our
own engine cannot do.
Click to drop a ball, B drops a box, R resets.
"""
import random

import pygame
import pymunk


WIDTH, HEIGHT = 900, 600
STEP = 1 / 120
MAX_STEPS = 8
GREY = (90, 96, 110)
TEXT_COLOR = (230, 230, 230)


def add_static_box(space, center, size):
    (cx, cy), (w, h) = center, size
    corners = [(cx - w / 2, cy - h / 2), (cx + w / 2, cy - h / 2), (cx + w / 2, cy + h / 2), (cx - w / 2, cy + h / 2)]
    shape = pymunk.Poly(space.static_body, corners)      # static shapes live on space.static_body
    shape.friction = 0.6
    shape.elasticity = 0.3
    space.add(shape)
    return shape


def add_ball(space, rng, pos):
    r = rng.randint(10, 22)
    mass = r * r / 100
    body = pymunk.Body(mass, pymunk.moment_for_circle(mass, 0, r))
    body.position = pos
    shape = pymunk.Circle(body, r)
    shape.friction = 0.6
    shape.elasticity = 0.4
    shape.color = (rng.randint(90, 255), rng.randint(90, 255), rng.randint(90, 255), 255)
    space.add(body, shape)
    return body


def add_box(space, pos, size=(50, 36)):
    mass = size[0] * size[1] / 1000
    body = pymunk.Body(mass, pymunk.moment_for_box(mass, size))
    body.position = pos
    shape = pymunk.Poly.create_box(body, size)
    shape.friction = 0.7
    shape.elasticity = 0.1
    shape.color = (210, 150, 90, 255)
    space.add(body, shape)
    return body


def make_space(rng):
    space = pymunk.Space()
    space.gravity = (0, 980)                     # y grows down, just like pygame
    space.damping = 0.95                         # keep 95% of velocity per second
    add_static_box(space, (WIDTH / 2, HEIGHT - 20), (WIDTH, 40))        # floor
    add_static_box(space, (-20, HEIGHT / 2), (40, HEIGHT))              # left wall
    add_static_box(space, (WIDTH + 20, HEIGHT / 2), (40, HEIGHT))       # right wall
    add_static_box(space, (200, 260), (300, 20))                        # ledge
    for i in range(5):
        peg = pymunk.Circle(space.static_body, 8, offset=(375 + i * 60, 400))
        peg.friction, peg.elasticity = 0.6, 0.3
        space.add(peg)

    bob = pymunk.Body(6, pymunk.moment_for_circle(6, 0, 24))
    bob.position = (840, 90)
    bob_shape = pymunk.Circle(bob, 24)
    bob_shape.elasticity, bob_shape.friction = 0.2, 0.6
    bob_shape.color = (240, 120, 90, 255)
    rod = pymunk.PinJoint(space.static_body, bob, (700, 90), (0, 0))   # keeps the distance fixed
    space.add(bob, bob_shape, rod)

    for i in range(8):
        add_ball(space, rng, (120 + i * 70, 60 + (i % 3) * 40))
    stack = [add_box(space, (760, HEIGHT - 40 - 18 - level * 36)) for level in range(4)]
    return space, stack                                                 # a stack our engine can't do


def draw(screen, space):
    for shape in space.shapes:
        color = getattr(shape, "color", None) or GREY
        color = color[:3]
        if isinstance(shape, pymunk.Circle):
            center = shape.body.local_to_world(shape.offset)
            pygame.draw.circle(screen, color, center, shape.radius)
            spoke = shape.body.local_to_world(shape.offset + pymunk.Vec2d(shape.radius, 0))
            pygame.draw.line(screen, (30, 30, 30), center, spoke, 2)
        elif isinstance(shape, pymunk.Poly):
            points = [shape.body.local_to_world(v) for v in shape.get_vertices()]
            pygame.draw.polygon(screen, color, points)
            pygame.draw.polygon(screen, (20, 20, 20), points, 2)
    for joint in space.constraints:
        a = joint.a.local_to_world(joint.anchor_a)
        b = joint.b.local_to_world(joint.anchor_b)
        pygame.draw.line(screen, (220, 220, 220), a, b, 2)


def main():
    pygame.init()
    screen = pygame.display.set_mode((WIDTH, HEIGHT))
    pygame.display.set_caption("The Same Scene in pymunk")
    clock = pygame.time.Clock()
    font = pygame.font.Font(None, 24)
    rng = random.Random(12)
    space, stack = make_space(rng)
    accumulator = 0.0
    steps = 0

    running = True
    while running:
        accumulator += min(clock.tick(60) / 1000, MAX_STEPS * STEP)
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                running = False
            elif event.type == pygame.MOUSEBUTTONDOWN and event.button == 1:
                add_ball(space, rng, event.pos)
            elif event.type == pygame.KEYDOWN and event.key == pygame.K_b:
                add_box(space, (rng.randint(420, 700), 40))
            elif event.type == pygame.KEYDOWN and event.key == pygame.K_r:
                space, stack = make_space(rng)

        while accumulator >= STEP:
            space.step(STEP)                     # pymunk wants a fixed step too
            accumulator -= STEP
            steps += 1

        screen.fill((24, 26, 34))
        draw(screen, space)
        hud = f"bodies {len(space.bodies)} | [click] ball  [B] box  [R] reset"
        screen.blit(font.render(hud, True, TEXT_COLOR), (10, 10))
        pygame.display.flip()

    pygame.quit()
    print(f"Ran {steps} pymunk steps with {len(space.bodies)} bodies")
    print(f"Box stack standing: {all(abs(box.angle) < 0.2 for box in stack)}")


if __name__ == "__main__":
    main()

Things to notice:

No "up" in pymunk. Gravity is just a vector. With pygame's y-down screen, (0, 980) pulls things down and you can draw positions directly, no flipping.
Mass and inertia are separate. pymunk.moment_for_circle(mass, 0, r) and moment_for_box compute the moment of inertia, with the formulas from SAT & Rotational Collisions. Pass float("inf") and the body can never spin.
Static shapes attach to space.static_body, the same idea as our inv_mass = 0.
space.damping = 0.95 means each body keeps 95% of its velocity per second, a per-second rate like ours.
Still a fixed step. pymunk's documentation recommends stepping with a constant dt, so the accumulator stays.

	Our mini engine	pymunk
Shapes	dynamic circles; static circles and boxes	circles, segments and convex polygons, all can be dynamic
Stacking rotating boxes	not supported	works (the stack in the program above)
Joints	one distance joint	pin, slide, pivot, groove, damped spring, motors and more
Resting bodies	always simulated	can sleep (`space.sleep_time_threshold`)
Game events	read `world.contacts` yourself	collision callbacks (`space.on_collision`) and point, segment and shape queries
You can read every line	yes, about 300 lines with drawing and the game loop	no, it is a compiled library

So which should you use? For a game whose physics is a few bouncing, rolling things, your own engine is small, debuggable and does exactly what you want. For piles, ragdolls, vehicles or anything with many joints, use a library, and spend your time on the game. Either way, having built one means you can read the library's documentation and know what every knob does.

🏋️ Practice Exercise: Mini Physics Engine

Objective: finish a small rigid-body engine in which balls bounce off pegs and ledges, roll along the floor and knock into a swinging pendulum, running on a fixed step.

Time: about 40 minutes, then 20 minutes for the pymunk version. Starter files: mini_engine_starter.py and pymunk_scene_starter.py (your instructor has them). The world, bodies, drawing and most of the solver are in place; right now balls fall through the floor, the rod does nothing and the timing depends on the frame rate. The numbered comments match the steps below.

Write circle_box: clamp the center into the box, return None if the closest point is out of reach, otherwise a Contact. Handle a center inside the box with the nearest face. Balls now land on the floor and ledge. (≈ 10 min)
Add the spin terms to effective_mass. Watch the spokes: balls that land moving sideways now roll. (≈ 4 min)
Write DistanceJoint.solve. The pendulum now swings on a rigid rod. (≈ 6 min)
Clamp the accumulated friction to ±limit. (≈ 3 min)
Replace the 0.99 damping with math.exp(-damping * dt). (≈ 3 min)
Replace the one-step-per-frame call with the fixed-step accumulator. (≈ 5 min)
Open pymunk_scene_starter.py and finish its four TODOs (gravity, moment of inertia, pin joint, accumulator). Compare it with your engine, then press B to drop boxes. (≈ 20 min)

You are done when:

balls bounce off the pegs, rest on the ledge and roll along the floor (the spokes turn as they move);
the pendulum swings on a rod that doesn't stretch: on exit the program prints a Pendulum length error below 1 px;
the scene runs at the same speed whatever your frame rate;
the pymunk version shows the same scene, and its box stack stands until you knock it over.

💡 Hint

For circle_box, closest is Vector2(max(lo.x, min(c.x, hi.x)), max(lo.y, min(c.y, hi.y))); if it equals the center, the center is inside the box. For the joint, the formula mirrors a contact without the clamp: a rod can push and pull, so its impulse may be negative. If balls sink slowly into the floor, check that c.target includes the BETA / dt push-out term.

✅ Example Solution

The lab file your instructor hands out also contains a few lines marked lab runtime and and frame_budget() in the loop condition. They let the checker run the program for a fixed number of frames; they do nothing when you run it yourself. The pymunk solution is the program in The Same Scene in pymunk.

"""Mini Physics Engine: Advanced Lesson 12 practice exercise (solution).

A small rigid-body engine: dynamic circles that roll, static boxes and
pegs, a pendulum on a distance joint, a fixed-step accumulator and an
iterative impulse solver with friction. Click to drop a ball, R resets.
"""
import math
import random
from dataclasses import dataclass

import pygame
from pygame import Vector2


WIDTH, HEIGHT = 900, 600
STEP = 1 / 120                  # fixed physics step (s)
MAX_STEPS = 8                   # never run more than this per frame (avoids the spiral of death)
ITERATIONS = 8                  # solver passes per step
BETA = 0.2                      # share of the overlap fixed per step (Baumgarte)
SLOP = 0.5                      # px of overlap we allow before correcting
RESTING_SPEED = 30              # px/s: slower impacts do not bounce
MAX_BALLS = 40
TEXT_COLOR = (230, 230, 230)


def cross(a, b):
    return a.x * b.y - a.y * b.x


def spin_velocity(omega, r):
    return Vector2(-omega * r.y, omega * r.x)


@dataclass
class Circle:
    radius: float


@dataclass
class Box:
    half: Vector2               # half width, half height (axis-aligned, static only)


class Body:
    def __init__(self, shape, pos, mass=0.0, restitution=0.3, friction=0.6, color=(200, 200, 200)):
        self.shape = shape
        self.pos = Vector2(pos)
        self.vel = Vector2()
        self.angle = 0.0
        self.omega = 0.0
        self.restitution = restitution
        self.friction = friction
        self.color = color
        if mass > 0:
            if not isinstance(shape, Circle):
                raise ValueError("this engine only moves circles; boxes must be static")
            self.inv_mass = 1 / mass
            self.inv_inertia = 1 / (0.5 * mass * shape.radius ** 2)   # solid disc
        else:
            self.inv_mass = self.inv_inertia = 0.0                   # static: infinite mass

    @property
    def is_static(self):
        return self.inv_mass == 0


class Contact:
    def __init__(self, a, b, normal, depth, point):
        self.a, self.b = a, b
        self.normal = normal            # unit vector from A to B
        self.depth = depth
        self.point = point
        self.jn = 0.0                   # accumulated normal impulse
        self.jt = 0.0                   # accumulated friction impulse


def circle_circle(a, b):
    delta = b.pos - a.pos
    radii = a.shape.radius + b.shape.radius
    dist_sq = delta.length_squared()
    if dist_sq >= radii * radii:
        return None
    dist = math.sqrt(dist_sq)
    normal = delta / dist if dist > 0 else Vector2(0, 1)
    return Contact(a, b, normal, radii - dist, a.pos + normal * a.shape.radius)


def circle_box(a, b):
    """A is a circle, B an axis-aligned box. The normal points from the circle into the box."""
    lo, hi = b.pos - b.shape.half, b.pos + b.shape.half
    c, r = a.pos, a.shape.radius
    closest = Vector2(max(lo.x, min(c.x, hi.x)), max(lo.y, min(c.y, hi.y)))
    delta = closest - c
    if delta.length_squared() > 0:                      # center outside the box
        dist = delta.length()
        if dist >= r:
            return None
        return Contact(a, b, delta / dist, r - dist, closest)
    # center inside the box: leave through the nearest face
    faces = [(c.x - lo.x, Vector2(1, 0)), (hi.x - c.x, Vector2(-1, 0)),
             (c.y - lo.y, Vector2(0, 1)), (hi.y - c.y, Vector2(0, -1))]
    gap, normal = min(faces, key=lambda f: f[0])
    return Contact(a, b, normal, r + gap, c - normal * gap)


def collide(a, b):
    if isinstance(a.shape, Circle) and isinstance(b.shape, Circle):
        return circle_circle(a, b)
    if isinstance(a.shape, Circle) and isinstance(b.shape, Box):
        return circle_box(a, b)
    if isinstance(a.shape, Box) and isinstance(b.shape, Circle):
        return circle_box(b, a)                         # keep the circle as A
    return None                                         # box-box: both static here


def apply_impulse(a, b, impulse, ra, rb):
    a.vel -= impulse * a.inv_mass
    a.omega -= cross(ra, impulse) * a.inv_inertia
    b.vel += impulse * b.inv_mass
    b.omega += cross(rb, impulse) * b.inv_inertia


def effective_mass(a, b, ra, rb, axis):
    """1 / (the mass the impulse 'feels' along axis), with the spin terms."""
    ran, rbn = cross(ra, axis), cross(rb, axis)
    return a.inv_mass + b.inv_mass + ran * ran * a.inv_inertia + rbn * rbn * b.inv_inertia


class DistanceJoint:
    """Keeps two body centers exactly `length` apart (a rigid rod)."""

    def __init__(self, a, b):
        self.a, self.b = a, b
        self.length = (b.pos - a.pos).length()

    def prepare(self, dt):
        delta = self.b.pos - self.a.pos
        dist = delta.length()
        self.axis = delta / dist if dist > 0 else Vector2(1, 0)
        self.bias = -BETA / dt * (dist - self.length)   # pull back toward the right length
        self.k = self.a.inv_mass + self.b.inv_mass

    def solve(self):
        if self.k == 0:
            return
        rel = (self.b.vel - self.a.vel).dot(self.axis)
        impulse = self.axis * ((self.bias - rel) / self.k)
        self.a.vel -= impulse * self.a.inv_mass
        self.b.vel += impulse * self.b.inv_mass


class World:
    def __init__(self, gravity=(0, 980), linear_damping=0.05, angular_damping=0.3):
        self.bodies = []
        self.joints = []
        self.gravity = Vector2(gravity)
        self.linear_damping = linear_damping            # per second
        self.angular_damping = angular_damping
        self.contacts = []

    def add(self, body):
        self.bodies.append(body)
        return body

    def find_contacts(self):
        found = []
        for i, a in enumerate(self.bodies):
            for b in self.bodies[i + 1:]:
                if a.is_static and b.is_static:
                    continue
                contact = collide(a, b)
                if contact is not None:
                    found.append(contact)
        return found

    def solve_contact(self, c, dt):
        a, b, n = c.a, c.b, c.normal
        t = Vector2(-n.y, n.x)
        ra, rb = c.point - a.pos, c.point - b.pos
        rel = (b.vel + spin_velocity(b.omega, rb)) - (a.vel + spin_velocity(a.omega, ra))
        vn = rel.dot(n)
        # normal impulse, accumulated and never negative (contacts push, never pull)
        dj = (c.target - vn) / effective_mass(a, b, ra, rb, n)
        new_jn = max(c.jn + dj, 0.0)
        apply_impulse(a, b, n * (new_jn - c.jn), ra, rb)
        c.jn = new_jn
        # friction impulse, clamped by Coulomb's law to mu * normal impulse
        rel = (b.vel + spin_velocity(b.omega, rb)) - (a.vel + spin_velocity(a.omega, ra))
        djt = -rel.dot(t) / effective_mass(a, b, ra, rb, t)
        limit = math.sqrt(a.friction * b.friction) * c.jn
        new_jt = max(-limit, min(limit, c.jt + djt))
        apply_impulse(a, b, t * (new_jt - c.jt), ra, rb)
        c.jt = new_jt

    def step(self, dt):
        """One fixed step: forces -> contacts and joints -> solve velocities -> move."""
        lin = math.exp(-self.linear_damping * dt)        # dt-scaled damping
        ang = math.exp(-self.angular_damping * dt)
        for body in self.bodies:
            if not body.is_static:
                body.vel += self.gravity * dt
                body.vel *= lin
                body.omega *= ang
        self.contacts = self.find_contacts()
        for c in self.contacts:
            ra, rb = c.point - c.a.pos, c.point - c.b.pos
            rel = (c.b.vel + spin_velocity(c.b.omega, rb)) - (c.a.vel + spin_velocity(c.a.omega, ra))
            vn = rel.dot(c.normal)
            e = min(c.a.restitution, c.b.restitution)
            bounce = -e * vn if vn < -RESTING_SPEED else 0.0
            c.target = max(bounce, BETA / dt * max(c.depth - SLOP, 0.0))
        for joint in self.joints:
            joint.prepare(dt)
        for _ in range(ITERATIONS):
            for c in self.contacts:
                self.solve_contact(c, dt)
            for joint in self.joints:
                joint.solve()
        for body in self.bodies:                         # semi-implicit Euler: new velocity moves the body
            if not body.is_static:
                body.pos += body.vel * dt
                body.angle += body.omega * dt


def ball(rng, pos):
    r = rng.randint(10, 22)
    color = (rng.randint(90, 255), rng.randint(90, 255), rng.randint(90, 255))
    return Body(Circle(r), pos, mass=r * r / 100, restitution=0.4, friction=0.6, color=color)


def make_world(rng):
    world = World()
    grey = (90, 96, 110)
    world.add(Body(Box(Vector2(WIDTH / 2, 20)), (WIDTH / 2, HEIGHT - 20), color=grey))     # floor
    world.add(Body(Box(Vector2(20, HEIGHT / 2)), (-20, HEIGHT / 2), color=grey))           # left wall
    world.add(Body(Box(Vector2(20, HEIGHT / 2)), (WIDTH + 20, HEIGHT / 2), color=grey))    # right wall
    world.add(Body(Box(Vector2(150, 10)), (200, 260), color=grey))                           # ledge
    for i in range(5):
        world.add(Body(Circle(8), (375 + i * 60, 400), color=grey))                          # pegs
    anchor = world.add(Body(Circle(6), (700, 90), color=grey))
    bob = world.add(Body(Circle(24), (840, 90), mass=6, restitution=0.2, color=(240, 120, 90)))
    world.joints.append(DistanceJoint(anchor, bob))
    for i in range(8):
        world.add(ball(rng, (120 + i * 70, 60 + (i % 3) * 40)))
    return world


def draw(screen, world):
    for body in world.bodies:
        if isinstance(body.shape, Box):
            rect = pygame.FRect(body.pos - body.shape.half, body.shape.half * 2)
            pygame.draw.rect(screen, body.color, rect)
        else:
            r = body.shape.radius
            pygame.draw.circle(screen, body.color, body.pos, r)
            spoke = body.pos + Vector2(r, 0).rotate_rad(body.angle)
            pygame.draw.line(screen, (30, 30, 30), body.pos, spoke, 2)     # shows the spin
    for joint in world.joints:
        pygame.draw.line(screen, (220, 220, 220), joint.a.pos, joint.b.pos, 2)


def main():
    pygame.init()
    screen = pygame.display.set_mode((WIDTH, HEIGHT))
    pygame.display.set_caption("Mini Physics Engine")
    clock = pygame.time.Clock()
    font = pygame.font.Font(None, 24)
    rng = random.Random(12)
    world = make_world(rng)
    accumulator = 0.0
    steps = 0

    running = True
    while running:
        accumulator += min(clock.tick(60) / 1000, MAX_STEPS * STEP)
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                running = False
            elif event.type == pygame.MOUSEBUTTONDOWN and event.button == 1:
                if sum(1 for b in world.bodies if not b.is_static) < MAX_BALLS:
                    world.add(ball(rng, event.pos))
            elif event.type == pygame.KEYDOWN and event.key == pygame.K_r:
                world = make_world(rng)

        while accumulator >= STEP:
            world.step(STEP)
            accumulator -= STEP
            steps += 1

        screen.fill((24, 26, 34))
        draw(screen, world)
        moving = sum(1 for b in world.bodies if not b.is_static)
        hud = f"bodies {moving} | contacts {len(world.contacts)} | [click] drop ball  [R] reset"
        screen.blit(font.render(hud, True, TEXT_COLOR), (10, 10))
        pygame.display.flip()

    pygame.quit()
    joint = world.joints[0]
    stretch = abs((joint.b.pos - joint.a.pos).length() - joint.length)
    print(f"Ran {steps} fixed steps of {STEP:.4f} s")
    print(f"Pendulum length error: {stretch:.2f} px")


if __name__ == "__main__":
    main()

📓 Learning Journal

Take five minutes to write in your learning journal. Jot down:

Key concepts you learned today
Techniques that clicked (and the ones that haven't, yet)
Questions or confusion to bring to the next session
Ideas to try in your own game
Progress and feelings: how did this lesson go for you?

✍️ This lesson's prompts:

Walk through one call of World.step in your own words, stage by stage. What would go wrong if positions moved before the solver ran?
You measured four integrators. Which would you choose for a physics puzzle game, for a space game with orbits, and for a trajectory preview line? Why?
After building an engine and using pymunk, where would you draw the line between "write it myself" and "use a library" in your own projects?

📝 Summary

You assembled a physics engine from parts you already understood. A World owns bodies, shapes, contacts and joints and advances them with a fixed step: forces and exponential damping change velocities, contacts are found, a solver adjusts accumulated impulses over several passes, and finally positions move. You measured four integrators on a spring and saw why engines favor semi-implicit Euler, gave static bodies zero inverse mass and inertia, made balls roll through the spin terms in the effective mass, and kept a pendulum rigid with a velocity-level joint. Then you rebuilt the scene in pymunk and compared what each approach gives you.

🎓 Key Takeaways

An engine is data (bodies, shapes, contacts, joints) plus a step(dt) with a fixed order: forces, detect, solve, move.
Run physics on a fixed step with a capped accumulator; damp with exp(-k * dt), never a per-step constant.
Semi-implicit Euler: one force evaluation and a single velocity update that the solver can correct. Verlet and RK4 are more accurate for smooth forces.
Static bodies are just zero inverse mass and zero inverse inertia.
Accumulated, clamped impulses let later solver passes correct earlier ones; joints must act on velocities, not only positions.
pymunk offers the same concepts plus polygons, stacking, many joints, sleeping and callbacks; it still wants a fixed step.

🔭 Looking Ahead

That completes the physics module. Next, in moderngl Foundations, you move drawing onto the graphics card with shaders, and learn the numpy you need to feed it data.

❓ Common Questions

How many solver iterations should I use?

Enough that your scenes look calm, and no more. More passes make piles and chains stiffer and cost proportionally more time per step. Our engine uses 8; pymunk's space.iterations defaults to 10. Try 2 and 20 in the lab and watch the pendulum and the pile.

Should I interpolate between physics steps when drawing?

If you see stutter because frames and steps don't line up, yes: keep the previous position of each body and draw at prev + (pos - prev) * (accumulator / STEP). At 120 steps per second and 60 frames per second the lab doesn't need it, but a 30-step engine on a 144 Hz display would.

Can I hold Shift to spawn heavier balls?

Yes, but check the modifier state, not a key constant: pygame has K_LSHIFT and K_RSHIFT for the keys, and pygame.key.get_mods() & pygame.KMOD_SHIFT for "either Shift is held". There is no pygame.K_SHIFT.

Are fixed-step results identical on every computer?

Don't count on it. The same program with the same inputs and step size repeats exactly on the same machine and Python build. Across different CPUs, operating systems or library versions, floating-point results can differ in the last bits, and physics amplifies those differences. Games that need identical results everywhere (lockstep multiplayer, for example) use fixed-point math or send state rather than inputs.

Why does pymunk ask for a moment of inertia at all?

Because only you know how your body's mass is distributed. A disc, a ring and a box of the same mass spin up differently. pymunk's moment_for_* helpers compute it for its standard shapes, and float("inf") gives a body that never rotates, handy for a player character.

🎯 Quick Quiz

Question 1: Why does the engine damp with vel *= math.exp(-linear_damping * dt) instead of vel *= 0.99 every step?

Question 2: How does the mini engine make a body static?

Question 3: Which integrator moves with the start-of-step acceleration, then re-evaluates the force at the new position and averages the two accelerations?

Question 4: Why does each contact keep an accumulated normal impulse c.jn, clamped so it never drops below 0?

Question 5: The mini engine uses gravity (0, 980) px/s² on pygame's y-down screen. What should space.gravity be for the same scene in pymunk?

🌟 Going Further

Sleeping bodies: give each body a timer that grows while its speed and spin stay tiny, skip it in step after half a second, and wake it when a contact with a moving body appears.
Springs: add a SpringJoint that applies Hooke's law as a force in stage 1, with damping. Then compare it with pymunk's DampedSpring.
Broad phase: reuse the spatial hash from Spatial Hashing & Object Pools so find_contacts only tests nearby pairs, and measure the difference with 200 balls.
Collision callbacks in pymunk: give balls and pegs collision_type values and use space.on_collision(...) to play a sound or count hits.
Read more: the pymunk documentation (tutorials, API reference and examples) and Glenn Fiedler's classic article Fix Your Timestep!.

🎯 Learning Objectives