diff --git a/changes/physics-jolt.md b/changes/physics-jolt.md
new file mode 100644
index 00000000..b90d59ec
--- /dev/null
+++ b/changes/physics-jolt.md
@@ -0,0 +1,8 @@
+bump: minor
+type: feat
+**`ludic.physics`: Jolt Physics as a package.** Shapes (box, sphere, capsule, cylinder, a
+boulder's dome, a heightfield, a mesh), still, kinematic and dynamic bodies that can be removed
+for real, rays, a foot's landing height, a push out of what a body stands in, overlaps, buoyancy
+against a `PhysWater` port with its current, and hard contacts and splashes as facts - stepped at a
+fixed 1/60 s, deterministic across machines. Jolt v5.6.0 is built here from its pinned tag
+(`native/build.sh`) through the phase-15 route.
diff --git a/packages/README.md b/packages/README.md
index fa65680d..444d193b 100644
--- a/packages/README.md
+++ b/packages/README.md
@@ -29,6 +29,7 @@ section. The rules are in [ludic.base](ludic.base/README.md).
 | [ludic.minimap](ludic.minimap/README.md) | a map's logic: the explored fog over a grid and its land share, the player's marks, and the view (world to map and back, zoom, pan, the camp, you, focus) |
 | [ludic.needs](ludic.needs/README.md) | a body's warmth, food, water and energy, and the countdown to a collapse |
 | [ludic.npc](ludic.npc/README.md) | the other people in a place: a routine by the hour, walking round what is in the way, facing a player who comes near, lines as data, a guest's copy |
+| [ludic.physics](ludic.physics/README.md) | rigid bodies, removable still shapes, queries and buoyancy over Jolt Physics (a native library, phase 16) |
 | [ludic.photo](ludic.photo/README.md) | a camera's photographs: what is in the frame, a grade on size, framing, light and the moment, the roll, its worth, the best of each subject |
 | [ludic.save](ludic.save/README.md) | versioned save files: a migration chain the game declares, torn writes told apart, a backup, a newer file refused and read-only |
 | [ludic.settings](ludic.settings/README.md) | a game's settings as data: one store, a fact per change, ranges, a safe set |
diff --git a/packages/ludic.physics/README.md b/packages/ludic.physics/README.md
new file mode 100644
index 00000000..642128f7
--- /dev/null
+++ b/packages/ludic.physics/README.md
@@ -0,0 +1,63 @@
+# ludic.physics
+
+Rigid bodies, still shapes that can be removed, queries and buoyancy, over
+[Jolt Physics](https://github.com/jrouwe/JoltPhysics) (MIT) built here from a pinned tag
+(phase 16). Uses `ludic.base` and nothing else.
+
+```ludic
+import "ludic.physics"
+```
+
+## The rules it keeps
+
+- **A package with a native library** (phase 15): `native/build.sh` fetches Jolt v5.6.0, checks
+  its SHA-256, and builds it with the shim `native/shim/jph_shim.cpp` into `lib/<target>/`
+  (`libludicjolt.dylib`, `ludicjolt.dll`). The shim keeps the rules in
+  [packages/README.md](../README.md#native-libraries): no callback into Ludic, contacts recorded
+  during the step and drained after it, Jolt's job threads inside C.
+- **Deterministic across machines**: `JPH_CROSS_PLATFORM_DETERMINISTIC` and no floating point
+  contraction, so the Mac and the PC agree in co-op and a run repeats to the bit (a test holds it).
+- **Fixed steps**: `phys_tick` runs as many 1/60 s steps as the frame's time holds, at most four -
+  frame rate never changes the rules, and a stall is not replayed.
+- **It saves nothing.** A thing at rest is saved by the Thing's owner; the world is rebuilt on
+  load and on a map swap (`phys_close`, `phys_open`).
+- **A shape and the world are handles it keeps**; a game names a shape and a body by an int.
+- **Removal is real**: `phys_remove` takes a body out of the world - no parking it far away.
+
+## The API
+
+| | |
+| --- | --- |
+| `phys_open(st, max_bodies, threads) -> bool`, `phys_close(st)`, `phys_is_open(st)` | a world, and everything in it let go |
+| `phys_box`, `phys_sphere`, `phys_capsule`, `phys_cylinder`, `phys_dome(r, h)`, `phys_offset(shape, x, y, z, yaw)` | shapes, each an int (`-1` refused); a dome is a boulder standing on its origin |
+| `phys_heightfield(st, h, n, ox, oz, cell)`, `phys_mesh(st, v, nv, tri, nt)` | the ground (`PHYS_HOLE` is none) and a dock or cabin |
+| `phys_ground_add`, `phys_static_add`, `phys_kinematic_add`, `phys_body_add(..., mass)` | bodies: the ground, still things, what the game moves (the player), what falls and floats |
+| `phys_settle(st)` | after adding many still things, before the first query |
+| `phys_remove(st, id)`, `phys_count(st)` | gone for good |
+| `phys_pose(st, id) -> PhysPose`, `phys_yaw_of(p)` | position, rotation, velocity |
+| `phys_place`, `phys_impulse`, `phys_set_velocity`, `phys_set_kinematic_target`, `phys_awake` | moving a body |
+| `phys_ray(st, o, d, mask) -> PhysHit` | the first thing along a ray |
+| `phys_top_at(st, x, z, r, from_y, depth, mask) -> float` | what a foot would land on (`CharacterGround.top_at`), `PHYS_NONE` if nothing |
+| `phys_push(st, x, z, r, y0, y1, mask) -> PhysPush` | the move that stands an upright body clear (`CharacterGround.push`) |
+| `phys_overlap(st, x, y, z, r, mask) -> []int` | the bodies a ball touches |
+| `phys_float(st, id, buoyancy, linear_drag, angular_drag)`, `phys_sink`, `phys_floating` | a body buoyed each step against `PhysWater` and carried by its current |
+| `phys_facts(st)` | `PhysFact`: `PHYS_HIT` (two bodies met at `ph_hard` m/s or more), `PHYS_SPLASH` (a floating body reached water) |
+| `phys_tick`, `phys_advance(st, dt) -> steps`, `phys_step`, `phys_system()` | the system, in `PH_SIMULATE` |
+| `phys_config(st, step, sub, hard)`, `phys_gravity(st, y)` | the numbers |
+
+A query's mask is `PHYS_M_GROUND`, `PHYS_M_STATIC`, `PHYS_M_MOVING`, or `PHYS_M_SOLID` /
+`PHYS_M_ALL` - the body's own questions leave the ground out, because `hk_ground` stays the pure
+ground query.
+
+## The port
+
+`PhysWater { height(x, z), current_x(x, z), current_z(x, z) }` - the surface under a point
+(`PHYS_NONE` where there is no water) and the current there. Unbound, there is no water. The game
+binds it to the same function its water is drawn with, so a boat never floats on a different wave
+from the one on screen.
+
+## Tests
+
+`ludic test packages/ludic.physics` runs the real library: a box onto a heightfield and a ray to
+it, still things removed, a boulder's top, a push out of a post, a crate floating level and
+drifting, hard and soft landings, a kinematic shove, fixed steps, and a pile run twice to the bit.
diff --git a/packages/ludic.physics/bodies.ludic b/packages/ludic.physics/bodies.ludic
new file mode 100644
index 00000000..6942331e
--- /dev/null
+++ b/packages/ludic.physics/bodies.ludic
@@ -0,0 +1,70 @@
+# bodies.ludic - bodies by id: still, kinematic or dynamic; removed for real; moved and pushed
+
+export property PhysPose {
+  x: float = 0.0
+  y: float = 0.0
+  z: float = 0.0
+  qx: float = 0.0
+  qy: float = 0.0
+  qz: float = 0.0
+  qw: float = 1.0
+  vx: float = 0.0
+  vy: float = 0.0
+  vz: float = 0.0
+}
+
+function ph_add(physics_st: mut PhysicsState, shape: int, x: float, y: float, z: float, yaw: float, motion: int, layer: int, mass: float) -> int {
+  let s = ph_shape(physics_st, shape)
+  if physics_st.ph_world == null or s == null { return -1 }
+  return jph_body_add(physics_st.ph_world, s, x, y, z, yaw, motion, layer, mass, 0.6, 0.1)
+}
+
+# the ground: a heightfield or a mesh at the world's origin
+export function phys_ground_add(physics_st: mut PhysicsState, shape: int) -> int { return ph_add(physics_st, shape, 0.0, 0.0, 0.0, 0.0, 0, 0, 0.0) }
+# a thing that never moves: a trunk, a boulder, a tent, a bear box
+export function phys_static_add(physics_st: mut PhysicsState, shape: int, x: float, y: float, z: float, yaw: float) -> int { return ph_add(physics_st, shape, x, y, z, yaw, 0, 1, 0.0) }
+# a thing moved by the game that pushes what it meets: the player's body
+export function phys_kinematic_add(physics_st: mut PhysicsState, shape: int, x: float, y: float, z: float, yaw: float) -> int { return ph_add(physics_st, shape, x, y, z, yaw, 1, 2, 0.0) }
+# a thing that falls, rolls and floats; mass in kilograms (0: from the shape at 1000 kg/m3)
+export function phys_body_add(physics_st: mut PhysicsState, shape: int, x: float, y: float, z: float, yaw: float, mass: float) -> int { return ph_add(physics_st, shape, x, y, z, yaw, 2, 2, mass) }
+
+# after many still things are added at once, before the first query
+export function phys_settle(physics_st: mut PhysicsState) -> void { if physics_st.ph_world != null { jph_world_optimize(physics_st.ph_world) } }
+
+export function phys_count(physics_st: PhysicsState) -> int {
+  if physics_st.ph_world == null { return 0 }
+  return jph_world_bodies(physics_st.ph_world)
+}
+
+# gone from the world, for good: what a struck tent needs
+export function phys_remove(physics_st: mut PhysicsState, id: int) -> void {
+  if physics_st.ph_world == null or id < 0 { return }
+  phys_sink(physics_st, id)
+  jph_body_remove(physics_st.ph_world, id)
+}
+
+# where a body is, how it is turned and how fast it goes; null when there is no such body
+export function phys_pose(physics_st: mut PhysicsState, id: int) -> PhysPose {
+  if physics_st.ph_world == null { return null }
+  ph_buffers(physics_st)
+  if jph_body_read(physics_st.ph_world, id, physics_st.ph_vals) != 0 { return null }
+  let v = physics_st.ph_vals
+  let p = new PhysPose
+  p.x = v[0];  p.y = v[1];  p.z = v[2]
+  p.qx = v[3];  p.qy = v[4];  p.qz = v[5];  p.qw = v[6]
+  p.vx = v[7];  p.vy = v[8];  p.vz = v[9]
+  return p
+}
+
+# the heading a pose's rotation gives, in radians about y
+export function phys_yaw_of(p: PhysPose) -> float { return Math.atan2(2.0 * (p.qw * p.qy + p.qx * p.qz), 1.0 - 2.0 * (p.qy * p.qy + p.qx * p.qx)) }
+
+export function phys_place(physics_st: mut PhysicsState, id: int, x: float, y: float, z: float, yaw: float) -> void { if physics_st.ph_world != null { jph_body_place(physics_st.ph_world, id, x, y, z, yaw) } }
+export function phys_impulse(physics_st: mut PhysicsState, id: int, x: float, y: float, z: float) -> void { if physics_st.ph_world != null { jph_body_impulse(physics_st.ph_world, id, x, y, z) } }
+export function phys_set_velocity(physics_st: mut PhysicsState, id: int, x: float, y: float, z: float) -> void { if physics_st.ph_world != null { jph_body_velocity(physics_st.ph_world, id, x, y, z) } }
+export function phys_awake(physics_st: PhysicsState, id: int) -> bool { return physics_st.ph_world != null and jph_body_awake(physics_st.ph_world, id) == 1 }
+
+# a kinematic body driven to where it should be one step from now, pushing what it meets
+export function phys_set_kinematic_target(physics_st: mut PhysicsState, id: int, x: float, y: float, z: float, yaw: float) -> void {
+  if physics_st.ph_world != null { jph_body_target(physics_st.ph_world, id, x, y, z, yaw, physics_st.ph_step) }
+}
diff --git a/packages/ludic.physics/index.ludic b/packages/ludic.physics/index.ludic
new file mode 100644
index 00000000..8f5aada1
--- /dev/null
+++ b/packages/ludic.physics/index.ludic
@@ -0,0 +1,12 @@
+# ludic.physics - a physics world over Jolt: shapes as handles it keeps, bodies by id, queries,
+# water that floats a body, and the contacts that matter as facts. Fixed 1/60 steps.
+module ludic_physics uses ludic_base
+numbers float
+import "ludic.base"
+import "native.ludic"
+import "state.ludic"
+import "shapes.ludic"
+import "bodies.ludic"
+import "queries.ludic"
+import "water.ludic"
+import "system.ludic"
diff --git a/packages/ludic.physics/lib/macos-arm64/libludicjolt.dylib b/packages/ludic.physics/lib/macos-arm64/libludicjolt.dylib
new file mode 100755
index 00000000..d6266968
--- /dev/null
+++ b/packages/ludic.physics/lib/macos-arm64/libludicjolt.dylib
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8cc67c97be16787659f781337f0c80551aab3663b6a216e96ac43bcb3eecda47
+size 1099872
diff --git a/packages/ludic.physics/native.ludic b/packages/ludic.physics/native.ludic
new file mode 100644
index 00000000..2b5a5ec7
--- /dev/null
+++ b/packages/ludic.physics/native.ludic
@@ -0,0 +1,32 @@
+# native.ludic - the shim's symbols (native/shim/jph_shim.cpp). A world and a shape are handles
+# the package keeps; a body is an int id. None of this is exported.
+extern function jph_world_new(max_bodies: int, threads: int) -> pointer = "jph_world_new"
+extern function jph_world_free(w: pointer) -> void = "jph_world_free"
+extern function jph_world_gravity(w: pointer, y: float) -> void = "jph_world_gravity"
+extern function jph_world_step(w: pointer, dt: float, sub: int) -> int = "jph_world_step"
+extern function jph_world_optimize(w: pointer) -> void = "jph_world_optimize"
+extern function jph_world_bodies(w: pointer) -> int = "jph_world_bodies"
+extern function jph_contacts_drain(w: pointer, ids: pointer, vals: pointer, cap: int) -> int = "jph_contacts_drain"
+extern function jph_contacts_lost(w: pointer) -> int = "jph_contacts_lost"
+extern function jph_shape_box(hx: float, hy: float, hz: float) -> pointer = "jph_shape_box"
+extern function jph_shape_sphere(r: float) -> pointer = "jph_shape_sphere"
+extern function jph_shape_capsule(half_h: float, r: float) -> pointer = "jph_shape_capsule"
+extern function jph_shape_cylinder(half_h: float, r: float) -> pointer = "jph_shape_cylinder"
+extern function jph_shape_dome(r: float, h: float) -> pointer = "jph_shape_dome"
+extern function jph_shape_offset(s: pointer, x: float, y: float, z: float, yaw: float) -> pointer = "jph_shape_offset"
+extern function jph_shape_heightfield(h: pointer, n: int, ox: float, oz: float, cell: float) -> pointer = "jph_shape_heightfield"
+extern function jph_shape_mesh(v: pointer, nv: int, tri: pointer, nt: int) -> pointer = "jph_shape_mesh"
+extern function jph_shape_free(s: pointer) -> void = "jph_shape_free"
+extern function jph_body_add(w: pointer, s: pointer, x: float, y: float, z: float, yaw: float, motion: int, layer: int, mass: float, friction: float, bounce: float) -> int = "jph_body_add"
+extern function jph_body_remove(w: pointer, id: int) -> void = "jph_body_remove"
+extern function jph_body_read(w: pointer, id: int, out: pointer) -> int = "jph_body_read"
+extern function jph_body_place(w: pointer, id: int, x: float, y: float, z: float, yaw: float) -> void = "jph_body_place"
+extern function jph_body_impulse(w: pointer, id: int, x: float, y: float, z: float) -> void = "jph_body_impulse"
+extern function jph_body_velocity(w: pointer, id: int, x: float, y: float, z: float) -> void = "jph_body_velocity"
+extern function jph_body_target(w: pointer, id: int, x: float, y: float, z: float, yaw: float, dt: float) -> void = "jph_body_target"
+extern function jph_body_awake(w: pointer, id: int) -> int = "jph_body_awake"
+extern function jph_body_buoyancy(w: pointer, id: int, sy: float, nx: float, ny: float, nz: float, buoy: float, lin: float, ang: float, cx: float, cz: float, dt: float) -> int = "jph_body_buoyancy"
+extern function jph_ray(w: pointer, ox: float, oy: float, oz: float, dx: float, dy: float, dz: float, mask: int, out: pointer) -> int = "jph_ray"
+extern function jph_top_at(w: pointer, x: float, z: float, r: float, from_y: float, depth: float, mask: int) -> float = "jph_top_at"
+extern function jph_overlap(w: pointer, x: float, y: float, z: float, r: float, mask: int, ids: pointer, cap: int) -> int = "jph_overlap"
+extern function jph_push(w: pointer, x: float, z: float, r: float, y0: float, y1: float, mask: int, out: pointer) -> int = "jph_push"
diff --git a/packages/ludic.physics/native/build.sh b/packages/ludic.physics/native/build.sh
new file mode 100755
index 00000000..191b443e
--- /dev/null
+++ b/packages/ludic.physics/native/build.sh
@@ -0,0 +1,28 @@
+#!/bin/sh
+# builds lib/<target>/ for ludic.physics: Jolt Physics at a pinned tag, and the shim over it.
+# Cross-platform determinism is on (the Mac and the PC must agree in co-op), and floating point
+# contraction off, which it needs. Objects go to build/native; about a minute on eight cores.
+set -eu
+PKG="$(cd "$(dirname "$0")/.." && pwd)"
+. "$PKG/../../tools/native/lib.sh"
+JOLT_TAG=v5.6.0
+JOLT_SHA=6e069ee0172478cc78182047aac87e5310ba14a67a53348ae14cc37801fd3f8e
+SRC="$PKG/build/src/jolt-$JOLT_TAG"
+native_fetch "$SRC" "https://github.com/jrouwe/JoltPhysics/archive/refs/tags/$JOLT_TAG.tar.gz" "$JOLT_SHA"
+
+OBJ="$PKG/build/native/$(native_target)"
+mkdir -p "$OBJ"
+CXX="$(native_cxx)"
+DEFS="-DNDEBUG -DJPH_CROSS_PLATFORM_DETERMINISTIC -DJPH_OBJECT_LAYER_BITS=16"
+FLAGS="$(native_cflags) -std=c++17 -ffp-contract=off -fno-exceptions -fno-rtti $DEFS -I$SRC"
+case "$(native_target)" in windows-*) FLAGS="$FLAGS -D_CRT_SECURE_NO_WARNINGS" ;; esac
+JOBS="$(getconf _NPROCESSORS_ONLN 2>/dev/null || echo 4)"
+
+# one object per Jolt source, named by its path so two files called the same do not collide
+( cd "$SRC" && find Jolt -name '*.cpp' ) | while read -r f; do
+  o="$OBJ/$(echo "$f" | tr '/' '_' | sed 's/\.cpp$/.o/')"
+  [ "$o" -nt "$SRC/$f" ] || echo "$f $o"
+done | xargs -P "$JOBS" -n 2 sh -c '"$0" '"$FLAGS"' -c "'"$SRC"'/$1" -o "$2"' "$CXX"
+"$CXX" $FLAGS -I"$PKG/native/shim" -c "$PKG/native/shim/jph_shim.cpp" -o "$OBJ/jph_shim.o"
+
+native_link "$PKG" ludicjolt "$OBJ"/*.o
diff --git a/packages/ludic.physics/native/shim/jph_bodies.inl b/packages/ludic.physics/native/shim/jph_bodies.inl
new file mode 100644
index 00000000..2ed464b9
--- /dev/null
+++ b/packages/ludic.physics/native/shim/jph_bodies.inl
@@ -0,0 +1,64 @@
+// jph_bodies.inl - bodies by id: added, removed, moved, pushed, read back
+
+// motion 0 still, 1 kinematic, 2 dynamic; layer 0 ground, 1 static, 2 moving; mass 0 = from shape
+JPH_SHIM int jph_body_add(void *w, void *shape, float x, float y, float z, float yaw, int motion, int layer,
+                          float mass, float friction, float bounce) {
+  World *ww = static_cast<World *>(w);
+  EMotionType mt = motion == 2 ? EMotionType::Dynamic : (motion == 1 ? EMotionType::Kinematic : EMotionType::Static);
+  BodyCreationSettings st(static_cast<Shape *>(shape), RVec3(x, y, z), yaw_q(yaw), mt, ObjectLayer(layer));
+  st.mFriction = friction;
+  st.mRestitution = bounce;
+  if (mt == EMotionType::Dynamic && mass > 0.0f) {
+    st.mOverrideMassProperties = EOverrideMassProperties::CalculateInertia;
+    st.mMassPropertiesOverride.mMass = mass;
+  }
+  EActivation act = mt == EMotionType::Static ? EActivation::DontActivate : EActivation::Activate;
+  BodyID id = ww->sys.GetBodyInterface().CreateAndAddBody(st, act);
+  if (id.IsInvalid()) return -1;
+  return int(id.GetIndexAndSequenceNumber());
+}
+
+JPH_SHIM void jph_body_remove(void *w, int id) {
+  BodyInterface &bi = static_cast<World *>(w)->sys.GetBodyInterface();
+  bi.RemoveBody(bid(id));
+  bi.DestroyBody(bid(id));
+}
+
+// position, rotation (x, y, z, w) and velocity into out[10]
+JPH_SHIM int jph_body_read(void *w, int id, float *out) {
+  BodyInterface &bi = static_cast<World *>(w)->sys.GetBodyInterface();
+  if (!bi.IsAdded(bid(id))) return -1;
+  RVec3 p; Quat q;
+  bi.GetPositionAndRotation(bid(id), p, q);
+  Vec3 v = bi.GetLinearVelocity(bid(id));
+  float vals[10] = { float(p.GetX()), float(p.GetY()), float(p.GetZ()), q.GetX(), q.GetY(), q.GetZ(), q.GetW(),
+                     v.GetX(), v.GetY(), v.GetZ() };
+  for (int i = 0; i < 10; i++) out[i] = vals[i];
+  return 0;
+}
+
+JPH_SHIM void jph_body_place(void *w, int id, float x, float y, float z, float yaw) {
+  static_cast<World *>(w)->sys.GetBodyInterface().SetPositionAndRotation(bid(id), RVec3(x, y, z), yaw_q(yaw), EActivation::Activate);
+}
+JPH_SHIM void jph_body_impulse(void *w, int id, float x, float y, float z) {
+  static_cast<World *>(w)->sys.GetBodyInterface().AddImpulse(bid(id), Vec3(x, y, z));
+}
+JPH_SHIM void jph_body_velocity(void *w, int id, float x, float y, float z) {
+  static_cast<World *>(w)->sys.GetBodyInterface().SetLinearVelocity(bid(id), Vec3(x, y, z));
+}
+// a kinematic body driven to where it should be after dt, pushing what it meets
+JPH_SHIM void jph_body_target(void *w, int id, float x, float y, float z, float yaw, float dt) {
+  static_cast<World *>(w)->sys.GetBodyInterface().MoveKinematic(bid(id), RVec3(x, y, z), yaw_q(yaw), dt);
+}
+JPH_SHIM int jph_body_awake(void *w, int id) { return static_cast<World *>(w)->sys.GetBodyInterface().IsActive(bid(id)) ? 1 : 0; }
+
+// water under a body: the surface's height and normal, how buoyant, the drags and the current
+JPH_SHIM int jph_body_buoyancy(void *w, int id, float sy, float nx, float ny, float nz, float buoy, float lin, float ang,
+                               float cx, float cz, float dt) {
+  World *ww = static_cast<World *>(w);
+  BodyInterface &bi = ww->sys.GetBodyInterface();
+  RVec3 p = bi.GetPosition(bid(id));
+  bool in = bi.ApplyBuoyancyImpulse(bid(id), RVec3(p.GetX(), sy, p.GetZ()), Vec3(nx, ny, nz), buoy, lin, ang,
+                                     Vec3(cx, 0.0f, cz), ww->sys.GetGravity(), dt);
+  return in ? 1 : 0;
+}
diff --git a/packages/ludic.physics/native/shim/jph_queries.inl b/packages/ludic.physics/native/shim/jph_queries.inl
new file mode 100644
index 00000000..76642b22
--- /dev/null
+++ b/packages/ludic.physics/native/shim/jph_queries.inl
@@ -0,0 +1,81 @@
+// jph_queries.inl - what is there: a ray, the top a body could step onto, overlaps, a push out
+
+// a ray from o along d (its length the reach) against the layers in mask; the body id or -1, and
+// out[7] the fraction, the point and the normal
+JPH_SHIM int jph_ray(void *w, float ox, float oy, float oz, float dx, float dy, float dz, int mask, float *out) {
+  World *ww = static_cast<World *>(w);
+  RRayCast ray{ RVec3(ox, oy, oz), Vec3(dx, dy, dz) };
+  RayCastResult hit;
+  MaskFilter f(mask);
+  if (!ww->sys.GetNarrowPhaseQuery().CastRay(ray, hit, {}, f)) return -1;
+  RVec3 p = ray.GetPointOnRay(hit.mFraction);
+  Vec3 n = Vec3::sAxisY();
+  BodyLockRead lock(ww->sys.GetBodyLockInterface(), hit.mBodyID);
+  if (lock.Succeeded()) n = lock.GetBody().GetWorldSpaceSurfaceNormal(hit.mSubShapeID2, p);
+  float vals[7] = { hit.mFraction, float(p.GetX()), float(p.GetY()), float(p.GetZ()), n.GetX(), n.GetY(), n.GetZ() };
+  for (int i = 0; i < 7; i++) out[i] = vals[i];
+  return int(hit.mBodyID.GetIndexAndSequenceNumber());
+}
+
+// a ball of radius r let down from (x, from_y, z) at most `depth`: the height of what it lands
+// on, or -FLT_MAX when it lands on nothing in mask
+JPH_SHIM float jph_top_at(void *w, float x, float z, float r, float from_y, float depth, int mask) {
+  World *ww = static_cast<World *>(w);
+  SphereShape ball(r);
+  ball.SetEmbedded();
+  RShapeCast cast = RShapeCast::sFromWorldTransform(&ball, Vec3::sOne(), RMat44::sTranslation(RVec3(x, from_y, z)),
+                                                    Vec3(0.0f, -depth, 0.0f));
+  ClosestHitCollisionCollector<CastShapeCollector> hit;
+  MaskFilter f(mask);
+  ww->sys.GetNarrowPhaseQuery().CastShape(cast, ShapeCastSettings(), RVec3::sZero(), hit, {}, f);
+  if (!hit.HadHit()) return -FLT_MAX;
+  return float(hit.mHit.mContactPointOn2.GetY());
+}
+
+// the bodies (in mask) a ball at (x, y, z) touches, up to cap into ids; returns how many
+JPH_SHIM int jph_overlap(void *w, float x, float y, float z, float r, int mask, int *ids, int cap) {
+  World *ww = static_cast<World *>(w);
+  SphereShape ball(r);
+  ball.SetEmbedded();
+  AllHitCollisionCollector<CollideShapeCollector> hits;
+  MaskFilter f(mask);
+  ww->sys.GetNarrowPhaseQuery().CollideShape(&ball, Vec3::sOne(), RMat44::sTranslation(RVec3(x, y, z)),
+                                             CollideShapeSettings(), RVec3::sZero(), hits, {}, f);
+  hits.Sort();
+  int n = 0;
+  for (const CollideShapeResult &h : hits.mHits) {
+    int id = int(h.mBodyID2.GetIndexAndSequenceNumber());
+    bool seen = false;
+    for (int i = 0; i < n; i++) seen = seen || ids[i] == id;
+    if (!seen && n < cap) ids[n++] = id;
+  }
+  return n;
+}
+
+// an upright body r wide, from y0 to y1, at (x, z): how far to move it in x and z to stand clear
+// of what it overlaps in mask, into out[2]; returns how many things it overlapped
+JPH_SHIM int jph_push(void *w, float x, float z, float r, float y0, float y1, int mask, float *out) {
+  World *ww = static_cast<World *>(w);
+  float half = (y1 - y0) * 0.5f;
+  CylinderShape body(half, r, 0.0f);
+  body.SetEmbedded();
+  AllHitCollisionCollector<CollideShapeCollector> hits;
+  MaskFilter f(mask);
+  CollideShapeSettings cs;
+  cs.mBackFaceMode = EBackFaceMode::CollideWithBackFaces;
+  ww->sys.GetNarrowPhaseQuery().CollideShape(&body, Vec3::sOne(), RMat44::sTranslation(RVec3(x, y0 + half, z)),
+                                             cs, RVec3::sZero(), hits, {}, f);
+  hits.Sort();
+  float px = 0.0f, pz = 0.0f;
+  for (const CollideShapeResult &h : hits.mHits) {
+    Vec3 a = h.mPenetrationAxis;
+    Vec3 flat(a.GetX(), 0.0f, a.GetZ());
+    if (flat.LengthSq() < 1.0e-8f) continue;
+    flat = flat.Normalized() * h.mPenetrationDepth;
+    px -= flat.GetX();
+    pz -= flat.GetZ();
+  }
+  out[0] = px;
+  out[1] = pz;
+  return int(hits.mHits.size());
+}
diff --git a/packages/ludic.physics/native/shim/jph_shapes.inl b/packages/ludic.physics/native/shim/jph_shapes.inl
new file mode 100644
index 00000000..0e842a8b
--- /dev/null
+++ b/packages/ludic.physics/native/shim/jph_shapes.inl
@@ -0,0 +1,42 @@
+// jph_shapes.inl - shapes, each a handle the package keeps and frees; null when Jolt refused it
+
+JPH_SHIM void *jph_shape_box(float hx, float hy, float hz) {
+  return keep(BoxShapeSettings(Vec3(hx, hy, hz), 0.02f).Create());
+}
+JPH_SHIM void *jph_shape_sphere(float r) { return keep(SphereShapeSettings(r).Create()); }
+JPH_SHIM void *jph_shape_capsule(float half_h, float r) { return keep(CapsuleShapeSettings(half_h, r).Create()); }
+JPH_SHIM void *jph_shape_cylinder(float half_h, float r) {
+  float cr = r < 0.05f ? r * 0.5f : 0.02f;
+  return keep(CylinderShapeSettings(half_h, r, cr).Create());
+}
+
+// a dome r across and h tall standing on y = 0: the cap of a sphere whose centre is below
+JPH_SHIM void *jph_shape_dome(float r, float h) {
+  float big = (r * r + h * h) / (2.0f * h);
+  RefConst<Shape> s = SphereShapeSettings(big).Create().Get();
+  return keep(RotatedTranslatedShapeSettings(Vec3(0.0f, h - big, 0.0f), Quat::sIdentity(), s).Create());
+}
+
+// a shape moved and turned about y, inside its body
+JPH_SHIM void *jph_shape_offset(void *s, float x, float y, float z, float yaw) {
+  return keep(RotatedTranslatedShapeSettings(Vec3(x, y, z), yaw_q(yaw), static_cast<Shape *>(s)).Create());
+}
+
+// n x n heights (row z, column x) at ox, oz, cell metres apart; FLT_MAX is a hole
+JPH_SHIM void *jph_shape_heightfield(const float *h, int n, float ox, float oz, float cell) {
+  HeightFieldShapeSettings st(h, Vec3(ox, 0.0f, oz), Vec3(cell, 1.0f, cell), uint32(n));
+  return keep(st.Create());
+}
+
+// a triangle mesh for a dock or a cabin: nv vertices (x, y, z) and nt triangles (i, j, k)
+JPH_SHIM void *jph_shape_mesh(const float *v, int nv, const int *tri, int nt) {
+  VertexList vs;
+  vs.reserve(nv);
+  for (int i = 0; i < nv; i++) vs.push_back(Float3(v[i * 3], v[i * 3 + 1], v[i * 3 + 2]));
+  IndexedTriangleList ts;
+  ts.reserve(nt);
+  for (int i = 0; i < nt; i++) ts.push_back(IndexedTriangle(tri[i * 3], tri[i * 3 + 1], tri[i * 3 + 2]));
+  return keep(MeshShapeSettings(vs, ts).Create());
+}
+
+JPH_SHIM void jph_shape_free(void *s) { if (s) static_cast<Shape *>(s)->Release(); }
diff --git a/packages/ludic.physics/native/shim/jph_shim.cpp b/packages/ludic.physics/native/shim/jph_shim.cpp
new file mode 100644
index 00000000..09d1857a
--- /dev/null
+++ b/packages/ludic.physics/native/shim/jph_shim.cpp
@@ -0,0 +1,126 @@
+// jph_shim.cpp - ludic.physics' door into Jolt Physics (phase 16). It keeps the shim rules of
+// packages/README.md: int, long, float and opaque handles cross; no struct by value, no callback
+// into Ludic - contacts are recorded during the step and drained after it; errors are return
+// codes. Jolt's job threads stay in here.
+#include <Jolt/Jolt.h>
+#include <Jolt/RegisterTypes.h>
+#include <Jolt/Core/Factory.h>
+#include <Jolt/Core/TempAllocator.h>
+#include <Jolt/Core/JobSystemThreadPool.h>
+#include <Jolt/Physics/PhysicsSettings.h>
+#include <Jolt/Physics/PhysicsSystem.h>
+#include <Jolt/Physics/Collision/Shape/BoxShape.h>
+#include <Jolt/Physics/Collision/Shape/SphereShape.h>
+#include <Jolt/Physics/Collision/Shape/CapsuleShape.h>
+#include <Jolt/Physics/Collision/Shape/CylinderShape.h>
+#include <Jolt/Physics/Collision/Shape/HeightFieldShape.h>
+#include <Jolt/Physics/Collision/Shape/MeshShape.h>
+#include <Jolt/Physics/Collision/Shape/RotatedTranslatedShape.h>
+#include <Jolt/Physics/Collision/RayCast.h>
+#include <Jolt/Physics/Collision/CastResult.h>
+#include <Jolt/Physics/Collision/ShapeCast.h>
+#include <Jolt/Physics/Collision/CollideShape.h>
+#include <Jolt/Physics/Collision/CollisionCollectorImpl.h>
+#include <Jolt/Physics/Body/BodyCreationSettings.h>
+#include <mutex>
+#include <thread>
+#include <cfloat>
+
+#if defined(_WIN32)
+#define JPH_SHIM extern "C" __declspec(dllexport)
+#else
+#define JPH_SHIM extern "C" __attribute__((visibility("default")))
+#endif
+
+using namespace JPH;
+
+// ---- layers: the ground, the still things on it, and what moves --------------------------
+namespace {
+const ObjectLayer L_GROUND = 0, L_STATIC = 1, L_MOVING = 2, L_COUNT = 3;
+const BroadPhaseLayer BP_STILL(0), BP_MOVING(1);
+
+class BPLayers final : public BroadPhaseLayerInterface {
+public:
+  uint GetNumBroadPhaseLayers() const override { return 2; }
+  BroadPhaseLayer GetBroadPhaseLayer(ObjectLayer l) const override { return l == L_MOVING ? BP_MOVING : BP_STILL; }
+#if defined(JPH_EXTERNAL_PROFILE) || defined(JPH_PROFILE_ENABLED)
+  const char *GetBroadPhaseLayerName(BroadPhaseLayer l) const override { return l == BP_MOVING ? "moving" : "still"; }
+#endif
+};
+class ObjVsBP final : public ObjectVsBroadPhaseLayerFilter {
+public:
+  bool ShouldCollide(ObjectLayer l, BroadPhaseLayer bp) const override { return l == L_MOVING || bp == BP_MOVING; }
+};
+class ObjPairs final : public ObjectLayerPairFilter {
+public:
+  bool ShouldCollide(ObjectLayer a, ObjectLayer b) const override { return a == L_MOVING || b == L_MOVING; }
+};
+// a query's layers as a bit mask (1 ground, 2 static, 4 moving)
+class MaskFilter final : public ObjectLayerFilter {
+public:
+  explicit MaskFilter(int m) : mask(m) {}
+  bool ShouldCollide(ObjectLayer l) const override { return (mask >> l) & 1; }
+  int mask;
+};
+
+// ---- contacts: recorded on Jolt's threads, drained on the game's -------------------------
+const int C_CAP = 512;
+struct Contact { uint32 a, b; float x, y, z, speed; };
+
+class Contacts final : public ContactListener {
+public:
+  void OnContactAdded(const Body &b1, const Body &b2, const ContactManifold &m, ContactSettings &) override {
+    Vec3 n = m.mWorldSpaceNormal;
+    float speed = (b2.GetLinearVelocity() - b1.GetLinearVelocity()).Dot(n);
+    RVec3 p = m.GetWorldSpaceContactPointOn1(0);
+    std::lock_guard<std::mutex> g(lock);
+    if (count >= C_CAP) { lost++; return; }
+    items[count++] = { b1.GetID().GetIndexAndSequenceNumber(), b2.GetID().GetIndexAndSequenceNumber(),
+                       float(p.GetX()), float(p.GetY()), float(p.GetZ()), std::abs(speed) };
+  }
+  std::mutex lock;
+  Contact items[C_CAP];
+  int count = 0;
+  int lost = 0;
+};
+
+struct World {
+  World(int max_bodies, int threads)
+      : temp(16 * 1024 * 1024), jobs(cMaxPhysicsJobs, cMaxPhysicsBarriers, threads) {
+    sys.Init(max_bodies, 0, max_bodies, max_bodies * 4, bpl, ovb, olp);
+    sys.SetContactListener(&contacts);
+  }
+  TempAllocatorImpl temp;
+  JobSystemThreadPool jobs;
+  BPLayers bpl;
+  ObjVsBP ovb;
+  ObjPairs olp;
+  PhysicsSystem sys;
+  Contacts contacts;
+};
+
+bool g_ready = false;
+BodyID bid(int id) { return BodyID(uint32(id)); }
+Quat yaw_q(float yaw) { return Quat::sRotation(Vec3::sAxisY(), yaw); }
+void *keep(const ShapeSettings::ShapeResult &r) {
+  if (r.HasError()) return nullptr;
+  Shape *s = const_cast<Shape *>(r.Get().GetPtr());
+  s->AddRef();
+  return s;
+}
+} // namespace
+
+// Jolt's allocator, factory and types, once per process; 0 when ready
+JPH_SHIM int jph_init(void) {
+  if (g_ready) return 0;
+  RegisterDefaultAllocator();
+  Factory::sInstance = new Factory();
+  RegisterTypes();
+  g_ready = true;
+  return 0;
+}
+
+#include "jph_world.inl"
+#include "jph_shapes.inl"
+#include "jph_bodies.inl"
+#include "jph_queries.inl"
diff --git a/packages/ludic.physics/native/shim/jph_world.inl b/packages/ludic.physics/native/shim/jph_world.inl
new file mode 100644
index 00000000..41e2253b
--- /dev/null
+++ b/packages/ludic.physics/native/shim/jph_world.inl
@@ -0,0 +1,49 @@
+// jph_world.inl - a physics world: made, stepped, emptied of its contacts, freed
+
+// a world for up to max_bodies; threads 0 means the machine's cores less one, at most four
+JPH_SHIM void *jph_world_new(int max_bodies, int threads) {
+  if (jph_init() != 0) return nullptr;
+  if (threads <= 0) {
+    int hw = int(std::thread::hardware_concurrency()) - 1;
+    threads = hw < 1 ? 1 : (hw > 4 ? 4 : hw);
+  }
+  World *w = new World(max_bodies, threads);
+  w->sys.SetGravity(Vec3(0.0f, -9.81f, 0.0f));
+  return w;
+}
+
+JPH_SHIM void jph_world_free(void *w) { delete static_cast<World *>(w); }
+
+JPH_SHIM void jph_world_gravity(void *w, float y) { static_cast<World *>(w)->sys.SetGravity(Vec3(0.0f, y, 0.0f)); }
+
+// one step of dt seconds in `sub` collision steps; 0 when Jolt reports nothing wrong
+JPH_SHIM int jph_world_step(void *w, float dt, int sub) {
+  World *ww = static_cast<World *>(w);
+  EPhysicsUpdateError e = ww->sys.Update(dt, sub < 1 ? 1 : sub, &ww->temp, &ww->jobs);
+  return e == EPhysicsUpdateError::None ? 0 : -int(e);
+}
+
+// after adding many still bodies at once, before the first query
+JPH_SHIM void jph_world_optimize(void *w) { static_cast<World *>(w)->sys.OptimizeBroadPhase(); }
+
+JPH_SHIM int jph_world_bodies(void *w) { return int(static_cast<World *>(w)->sys.GetNumBodies()); }
+
+// up to cap contacts into ids (a, b pairs) and vals (x, y, z, speed); returns how many
+JPH_SHIM int jph_contacts_drain(void *w, int *ids, float *vals, int cap) {
+  Contacts &c = static_cast<World *>(w)->contacts;
+  std::lock_guard<std::mutex> g(c.lock);
+  int n = c.count < cap ? c.count : cap;
+  for (int i = 0; i < n; i++) {
+    ids[i * 2] = int(c.items[i].a);
+    ids[i * 2 + 1] = int(c.items[i].b);
+    vals[i * 4] = c.items[i].x;
+    vals[i * 4 + 1] = c.items[i].y;
+    vals[i * 4 + 2] = c.items[i].z;
+    vals[i * 4 + 3] = c.items[i].speed;
+  }
+  for (int i = n; i < c.count; i++) c.items[i - n] = c.items[i];
+  c.count -= n;
+  return n;
+}
+
+JPH_SHIM int jph_contacts_lost(void *w) { return static_cast<World *>(w)->contacts.lost; }
diff --git a/packages/ludic.physics/package.ludic b/packages/ludic.physics/package.ludic
new file mode 100644
index 00000000..26eaf1b3
--- /dev/null
+++ b/packages/ludic.physics/package.ludic
@@ -0,0 +1,7 @@
+# ludic.physics - rigid bodies, still shapes, queries and buoyancy over Jolt Physics (phase 16),
+# built here from a pinned tag (native/build.sh). Uses ludic.base and nothing else.
+package "ludic.physics"
+version "0.1.0"
+kind source
+native "macos-arm64" "lib/macos-arm64/libludicjolt.dylib"
+native "windows-x64" "lib/windows-x64/ludicjolt.dll"
diff --git a/packages/ludic.physics/queries.ludic b/packages/ludic.physics/queries.ludic
new file mode 100644
index 00000000..5febb98f
--- /dev/null
+++ b/packages/ludic.physics/queries.ludic
@@ -0,0 +1,64 @@
+# queries.ludic - what is there, asked of the world as it stood after the last step. Each takes a
+# mask of what it sees (PHYS_M_*), so a body's own questions can leave out the ground or itself.
+
+export property PhysHit {
+  body: int = -1        # -1: nothing
+  fraction: float = 1.0 # how far along the ray, 0 to 1
+  x: float = 0.0
+  y: float = 0.0
+  z: float = 0.0
+  nx: float = 0.0
+  ny: float = 1.0
+  nz: float = 0.0
+}
+
+# a ray from (ox, oy, oz) along (dx, dy, dz), whose length is the reach
+export function phys_ray(physics_st: mut PhysicsState, ox: float, oy: float, oz: float, dx: float, dy: float, dz: float, mask: int) -> PhysHit {
+  let h = new PhysHit
+  if physics_st.ph_world == null { return h }
+  ph_buffers(physics_st)
+  h.body = jph_ray(physics_st.ph_world, ox, oy, oz, dx, dy, dz, mask, physics_st.ph_vals)
+  if h.body == -1 { return h }
+  let v = physics_st.ph_vals
+  h.fraction = v[0]
+  h.x = v[1];  h.y = v[2];  h.z = v[3]
+  h.nx = v[4];  h.ny = v[5];  h.nz = v[6]
+  return h
+}
+
+# a ball r across let down from from_y at most depth metres: the height of what it lands on, or
+# PHYS_NONE - what CharacterGround.top_at asks, with the ground left out of the mask
+export function phys_top_at(physics_st: PhysicsState, x: float, z: float, r: float, from_y: float, depth: float, mask: int) -> float {
+  if physics_st.ph_world == null { return PHYS_NONE }
+  let y = jph_top_at(physics_st.ph_world, x, z, r, from_y, depth, mask)
+  if y < PHYS_NONE { return PHYS_NONE }
+  return y
+}
+
+# the bodies a ball at (x, y, z) touches, in the order Jolt sorts them
+export function phys_overlap(physics_st: mut PhysicsState, x: float, y: float, z: float, r: float, mask: int) -> []int {
+  let out = new []int
+  if physics_st.ph_world == null { return out }
+  ph_buffers(physics_st)
+  let n = jph_overlap(physics_st.ph_world, x, y, z, r, mask, physics_st.ph_ids, PH_CAP)
+  for i in 0 .. n { push(out, physics_st.ph_ids[i]) }
+  return out
+}
+
+export property PhysPush {
+  x: float = 0.0     # how far to move in x and z to stand clear
+  z: float = 0.0
+  n: int = 0         # how many things were in the way
+}
+
+# an upright body r wide from y0 to y1 at (x, z): the move that stands it clear - what
+# CharacterGround.push asks
+export function phys_push(physics_st: mut PhysicsState, x: float, z: float, r: float, y0: float, y1: float, mask: int) -> PhysPush {
+  let p = new PhysPush
+  if physics_st.ph_world == null { return p }
+  ph_buffers(physics_st)
+  p.n = jph_push(physics_st.ph_world, x, z, r, y0, y1, mask, physics_st.ph_vals)
+  p.x = physics_st.ph_vals[0]
+  p.z = physics_st.ph_vals[1]
+  return p
+}
diff --git a/packages/ludic.physics/shapes.ludic b/packages/ludic.physics/shapes.ludic
new file mode 100644
index 00000000..3d4cf6b4
--- /dev/null
+++ b/packages/ludic.physics/shapes.ludic
@@ -0,0 +1,71 @@
+# shapes.ludic - the world opened and closed, and shapes: made in C, kept here, named by an int.
+# A shape may be used by any number of bodies; phys_close frees them all.
+
+# a world for up to max_bodies (threads 0: the cores less one, at most four); false if C refused
+export function phys_open(physics_st: mut PhysicsState, max_bodies: int, threads: int) -> bool {
+  phys_close(physics_st)
+  ph_buffers(physics_st)
+  physics_st.ph_world = jph_world_new(max_bodies, threads)
+  physics_st.ph_acc = 0.0
+  physics_st.ph_steps = 0
+  return physics_st.ph_world != null
+}
+
+export function phys_is_open(physics_st: PhysicsState) -> bool { return physics_st.ph_world != null }
+
+# everything let go: the world, every body in it and every shape (a world swap, a quit)
+export function phys_close(physics_st: mut PhysicsState) -> void {
+  if physics_st.ph_world != null { jph_world_free(physics_st.ph_world) }
+  physics_st.ph_world = null
+  if physics_st.ph_shapes != null {
+    for i in 0 .. len(physics_st.ph_shapes) {
+      if physics_st.ph_shapes[i] != null { jph_shape_free(physics_st.ph_shapes[i]) }
+    }
+  }
+  physics_st.ph_shapes = new []pointer
+  physics_st.ph_floats = new []int
+  physics_st.ph_float_k = new []float
+  if physics_st.ph_facts != null { q_clear(physics_st.ph_facts) }
+}
+
+function ph_keep(physics_st: mut PhysicsState, s: pointer) -> int {
+  if s == null { return -1 }
+  ph_buffers(physics_st)
+  push(physics_st.ph_shapes, s)
+  return len(physics_st.ph_shapes) - 1
+}
+
+function ph_shape(physics_st: PhysicsState, id: int) -> pointer {
+  if physics_st.ph_shapes == null or id < 0 or id >= len(physics_st.ph_shapes) { return null }
+  return physics_st.ph_shapes[id]
+}
+
+# a crate: half extents in metres
+export function phys_box(physics_st: mut PhysicsState, hx: float, hy: float, hz: float) -> int { return ph_keep(physics_st, jph_shape_box(hx, hy, hz)) }
+export function phys_sphere(physics_st: mut PhysicsState, r: float) -> int { return ph_keep(physics_st, jph_shape_sphere(r)) }
+# a trunk or a body: upright, its straight part 2 * half_h tall, centred on its origin
+export function phys_capsule(physics_st: mut PhysicsState, half_h: float, r: float) -> int { return ph_keep(physics_st, jph_shape_capsule(half_h, r)) }
+# a post: upright, 2 * half_h tall, centred on its origin
+export function phys_cylinder(physics_st: mut PhysicsState, half_h: float, r: float) -> int { return ph_keep(physics_st, jph_shape_cylinder(half_h, r)) }
+# a boulder: a dome r across and h tall, standing on its origin
+export function phys_dome(physics_st: mut PhysicsState, r: float, h: float) -> int { return ph_keep(physics_st, jph_shape_dome(r, h)) }
+
+# another shape moved and turned about y inside its body (a post whose foot is its origin)
+export function phys_offset(physics_st: mut PhysicsState, shape: int, x: float, y: float, z: float, yaw: float) -> int {
+  let s = ph_shape(physics_st, shape)
+  if s == null { return -1 }
+  return ph_keep(physics_st, jph_shape_offset(s, x, y, z, yaw))
+}
+
+# n x n heights, row by row along z, the first at (ox, oz), cell metres apart; PHYS_HOLE is none
+export const PHYS_HOLE: float = 340282346638528859811704183484516925440.0
+export function phys_heightfield(physics_st: mut PhysicsState, h: []float, n: int, ox: float, oz: float, cell: float) -> int {
+  if len(h) < n * n { return -1 }
+  return ph_keep(physics_st, jph_shape_heightfield(h, n, ox, oz, cell))
+}
+
+# a dock or a cabin: v holds nv points (x, y, z), tri holds nt triangles (i, j, k)
+export function phys_mesh(physics_st: mut PhysicsState, v: []float, nv: int, tri: []int, nt: int) -> int {
+  if len(v) < nv * 3 or len(tri) < nt * 3 { return -1 }
+  return ph_keep(physics_st, jph_shape_mesh(v, nv, tri, nt))
+}
diff --git a/packages/ludic.physics/state.ludic b/packages/ludic.physics/state.ludic
new file mode 100644
index 00000000..b53b6a65
--- /dev/null
+++ b/packages/ludic.physics/state.ludic
@@ -0,0 +1,72 @@
+# state.ludic - the port, the kinds and layers, the facts, and what the world holds
+# what the package asks the world: the water's surface under a point and the current there.
+# Unbound there is no water.
+export port PhysWater {
+  height: fn(float, float) -> float = fn ph_dry
+  current_x: fn(float, float) -> float = fn ph_still
+  current_z: fn(float, float) -> float = fn ph_still
+}
+function ph_dry(x: float, z: float) -> float { return PHYS_NONE }
+function ph_still(x: float, z: float) -> float { return 0.0 }
+
+export const PHYS_NONE: float = -1000000000.0   # no height: no water, nothing under a query
+
+# what a query sees, as a mask: the ground, the still things on it, and what moves
+export const PHYS_M_GROUND: int = 1
+export const PHYS_M_STATIC: int = 2
+export const PHYS_M_MOVING: int = 4
+export const PHYS_M_SOLID: int = 3
+export const PHYS_M_ALL: int = 7
+
+# what happened
+export const PHYS_HIT: int = 0        # two bodies met at least ph_hard metres a second apart
+export const PHYS_SPLASH: int = 1     # a floating body reached the water
+
+export property PhysFact {
+  what: int = 0
+  a: int = -1
+  b: int = -1
+  x: float = 0.0
+  y: float = 0.0
+  z: float = 0.0
+  speed: float = 0.0
+}
+
+export const PH_CAP: int = 512        # contacts drained a step, queries' buffers
+
+export state PhysicsState {
+  ph_world: pointer = null
+  ph_shapes: []pointer = null     # a shape id is its place here; null once freed
+  ph_step: float = 0.016666668    # seconds a step
+  ph_sub: int = 1                 # collision steps in each
+  ph_max_steps: int = 4           # a slow frame catches up this far and no further
+  ph_acc: float = 0.0
+  ph_steps: int = 0
+  ph_hard: float = 2.0            # metres a second: a contact under this is not a fact
+  ph_floats: []int = null         # buoyant bodies' ids
+  ph_float_k: []float = null      # and theirs: buoyancy, linear drag, angular drag, wet
+  ph_ids: []int = null
+  ph_vals: []float = null
+  ph_facts: Queue<PhysFact> = null
+}
+
+# the facts, oldest first
+export function phys_facts(physics_st: mut PhysicsState) -> Queue<PhysFact> {
+  if physics_st.ph_facts == null { physics_st.ph_facts = queue_new("physics.facts") }
+  return physics_st.ph_facts
+}
+
+function ph_buffers(physics_st: mut PhysicsState) -> void {
+  if physics_st.ph_ids == null { physics_st.ph_ids = words(PH_CAP * 2) }
+  if physics_st.ph_vals == null { physics_st.ph_vals = floats(PH_CAP * 4) }
+  if physics_st.ph_shapes == null { physics_st.ph_shapes = new []pointer }
+  if physics_st.ph_floats == null { physics_st.ph_floats = new []int }
+  if physics_st.ph_float_k == null { physics_st.ph_float_k = new []float }
+}
+
+# a step's length, its collision steps, and how hard a contact must be to be a fact
+export function phys_config(physics_st: mut PhysicsState, step: float, sub: int, hard: float) -> void {
+  physics_st.ph_step = step
+  physics_st.ph_sub = sub
+  physics_st.ph_hard = hard
+}
diff --git a/packages/ludic.physics/system.ludic b/packages/ludic.physics/system.ludic
new file mode 100644
index 00000000..129db8f2
--- /dev/null
+++ b/packages/ludic.physics/system.ludic
@@ -0,0 +1,67 @@
+# system.ludic - the world stepped at a fixed rate in PH_SIMULATE, so frame rate never changes
+# the rules. It saves nothing: a thing at rest is saved by its owner and the world rebuilt on load.
+
+# as many fixed steps as the frame's time holds, at most ph_max_steps (a stall is not replayed)
+export function phys_tick(physics_st: mut PhysicsState, t: Tick) -> void { phys_advance(physics_st, t.dt) }
+
+export function phys_advance(physics_st: mut PhysicsState, dt: float) -> int {
+  if physics_st.ph_world == null { return 0 }
+  physics_st.ph_acc += dt
+  var n = 0
+  # a thousandth of a step of slack: 0.05 s is three steps even when float says 2.9999
+  while physics_st.ph_acc + physics_st.ph_step * 0.001 >= physics_st.ph_step and n < physics_st.ph_max_steps {
+    phys_step(physics_st)
+    physics_st.ph_acc -= physics_st.ph_step
+    n += 1
+  }
+  if physics_st.ph_acc > physics_st.ph_step or physics_st.ph_acc < 0.0 { physics_st.ph_acc = 0.0 }
+  return n
+}
+
+# exactly one step: the water, Jolt, then the contacts that matter as facts
+export function phys_step(physics_st: mut PhysicsState) -> void {
+  if physics_st.ph_world == null { return }
+  ph_buffers(physics_st)
+  ph_buoy_all(physics_st)
+  jph_world_step(physics_st.ph_world, physics_st.ph_step, physics_st.ph_sub)
+  physics_st.ph_steps += 1
+  ph_contacts(physics_st)
+}
+
+export function phys_steps(physics_st: PhysicsState) -> int { return physics_st.ph_steps }
+
+# m/s2 down is negative
+export function phys_gravity(physics_st: mut PhysicsState, y: float) -> void { if physics_st.ph_world != null { jph_world_gravity(physics_st.ph_world, y) } }
+
+function ph_contacts(physics_st: mut PhysicsState) -> void {
+  let n = jph_contacts_drain(physics_st.ph_world, physics_st.ph_ids, physics_st.ph_vals, PH_CAP)
+  for i in 0 .. n {
+    let speed = physics_st.ph_vals[i * 4 + 3]
+    if speed >= physics_st.ph_hard {
+      let f = new PhysFact
+      f.what = PHYS_HIT
+      f.a = physics_st.ph_ids[i * 2]
+      f.b = physics_st.ph_ids[i * 2 + 1]
+      f.x = physics_st.ph_vals[i * 4]
+      f.y = physics_st.ph_vals[i * 4 + 1]
+      f.z = physics_st.ph_vals[i * 4 + 2]
+      f.speed = speed
+      q_push(phys_facts(physics_st), f)
+    }
+  }
+}
+
+# contacts C had to drop because the buffer was full between two drains
+export function phys_contacts_lost(physics_st: PhysicsState) -> int {
+  if physics_st.ph_world == null { return 0 }
+  return jph_contacts_lost(physics_st.ph_world)
+}
+
+export function phys_reset(physics_st: mut PhysicsState) -> void { physics_st.ph_acc = 0.0 }
+
+export function phys_system() -> System {
+  let s = system_new("physics", PH_SIMULATE)
+  s.reset = fn phys_reset
+  s.tick = fn phys_tick
+  return s
+}
diff --git a/packages/ludic.physics/tests/motion_test.ludic b/packages/ludic.physics/tests/motion_test.ludic
new file mode 100644
index 00000000..1d7f2d12
--- /dev/null
+++ b/packages/ludic.physics/tests/motion_test.ludic
@@ -0,0 +1,130 @@
+# motion_test.ludic - what moves: a crate that floats in the fake lake and drifts with its current,
+# a hard landing as a fact and a soft one not, a kinematic body pushing a crate, fixed steps
+# whatever the frame, a world closed and opened again, and the same run twice to the bit
+import "ludic.physics"
+import "ludic.base"
+program MotionTest {
+  numbers float
+
+  function lake_h(x: float, z: float) -> float {
+    if x > 20.0 { return 0.0 }
+    return PHYS_NONE
+  }
+  function lake_cz(x: float, z: float) -> float { return 0.5 }
+  bind PhysWater {
+    height: fn lake_h
+    current_z: fn lake_cz
+  }
+
+  # a lake bed 4 m down east of x = 20, a bank at 0 west of it; 64 x 64, a metre apart
+  function world(physics_st: mut PhysicsState) -> void {
+    expect(phys_open(physics_st, 1024, 1))
+    let h = floats(64 * 64)
+    for j in 0 .. 64 {
+      for i in 0 .. 64 {
+        var y = 0.0
+        if i > 20 { y = -4.0 }
+        h[j * 64 + i] = y
+      }
+    }
+    phys_ground_add(physics_st, phys_heightfield(physics_st, h, 64, 0.0, -32.0, 1.0))
+  }
+
+  function sim(physics_st: mut PhysicsState, seconds: float) -> void {
+    for i in 0 .. int(seconds * 60.0) { phys_step(physics_st) }
+  }
+
+  function count(fs: []PhysFact, what: int) -> int {
+    var n = 0
+    for i in 0 .. len(fs) { if fs[i].what == what { n += 1 } }
+    return n
+  }
+
+  test "a crate floats level at the surface, drifts with the current, and splashes once" (physics_st: mut PhysicsState) {
+    world(physics_st)
+    let crate = phys_body_add(physics_st, phys_box(physics_st, 0.5, 0.25, 0.5), 30.0, 2.0, 0.0, 0.0, 60.0)
+    phys_float(physics_st, crate, 1.2, 2.0, 2.0)
+    sim(physics_st, 8.0)
+    let p = phys_pose(physics_st, crate)
+    expect(p.y > -0.4)
+    expect(p.y < 0.4)
+    expect(Math.abs(p.vy) < 0.2)
+    expect(Math.abs(p.qw) > 0.99)
+    expect(p.z > 0.5)
+    expect_eq(count(q_drain(phys_facts(physics_st)), PHYS_SPLASH), 1)
+    phys_sink(physics_st, crate)
+    expect(not phys_floating(physics_st, crate))
+    sim(physics_st, 4.0)
+    expect(phys_pose(physics_st, crate).y < -2.5)
+    phys_close(physics_st)
+  }
+
+  test "a hard landing is a fact, a gentle one is not" (physics_st: mut PhysicsState) {
+    world(physics_st)
+    let b = phys_box(physics_st, 0.3, 0.3, 0.3)
+    let hard = phys_body_add(physics_st, b, 5.0, 6.0, 0.0, 0.0, 10.0)
+    let soft = phys_body_add(physics_st, b, 10.0, 0.35, 0.0, 0.0, 10.0)
+    sim(physics_st, 2.0)
+    let fs = q_drain(phys_facts(physics_st))
+    var hard_hits = 0
+    var soft_hits = 0
+    for i in 0 .. len(fs) {
+      if fs[i].what == PHYS_HIT and (fs[i].a == hard or fs[i].b == hard) { hard_hits += 1;  expect(fs[i].speed >= 2.0) }
+      if fs[i].what == PHYS_HIT and (fs[i].a == soft or fs[i].b == soft) { soft_hits += 1 }
+    }
+    expect(hard_hits >= 1)
+    expect_eq(soft_hits, 0)
+    expect_eq(phys_contacts_lost(physics_st), 0)
+    phys_close(physics_st)
+  }
+
+  test "a kinematic body walking into a crate shoves it along" (physics_st: mut PhysicsState) {
+    world(physics_st)
+    let me = phys_kinematic_add(physics_st, phys_capsule(physics_st, 0.6, 0.3), 2.0, 0.95, 0.0, 0.0)
+    let crate = phys_body_add(physics_st, phys_box(physics_st, 0.3, 0.3, 0.3), 4.0, 0.3, 0.0, 0.0, 5.0)
+    sim(physics_st, 0.5)
+    for i in 0 .. 120 {
+      phys_set_kinematic_target(physics_st, me, 2.0 + float(i + 1) * 0.03, 0.95, 0.0, 0.0)
+      phys_step(physics_st)
+    }
+    expect(phys_pose(physics_st, crate).x > 5.0)
+    phys_close(physics_st)
+  }
+
+  test "steps are fixed whatever the frame, and a stall is not replayed" (physics_st: mut PhysicsState) {
+    world(physics_st)
+    expect_eq(phys_advance(physics_st, 0.05), 3)
+    expect_eq(phys_advance(physics_st, 0.01), 0)
+    expect_eq(phys_advance(physics_st, 0.01), 1)
+    expect_eq(phys_advance(physics_st, 2.0), 4)
+    expect_eq(phys_advance(physics_st, 0.0), 0)
+    expect_eq(phys_steps(physics_st), 8)
+    phys_close(physics_st)
+    expect(not phys_is_open(physics_st))
+    expect_eq(phys_advance(physics_st, 1.0), 0)
+    world(physics_st)
+    expect_eq(phys_steps(physics_st), 0)
+    phys_close(physics_st)
+  }
+
+  function tumble(physics_st: mut PhysicsState) -> PhysPose {
+    world(physics_st)
+    let b = phys_box(physics_st, 0.4, 0.2, 0.3)
+    var last = -1
+    for i in 0 .. 12 { last = phys_body_add(physics_st, b, 3.0 + float(i % 3) * 0.3, 2.0 + float(i) * 0.7, float(i % 2) * 0.2, float(i) * 0.4, 8.0) }
+    phys_impulse(physics_st, last, 20.0, 0.0, 5.0)
+    sim(physics_st, 4.0)
+    let p = phys_pose(physics_st, last)
+    phys_close(physics_st)
+    return p
+  }
+
+  test "the same pile, run twice, lands the same to the bit" (physics_st: mut PhysicsState) {
+    let a = tumble(physics_st)
+    let b = tumble(physics_st)
+    expect_eq(float_bits(a.x), float_bits(b.x))
+    expect_eq(float_bits(a.y), float_bits(b.y))
+    expect_eq(float_bits(a.z), float_bits(b.z))
+    expect_eq(float_bits(a.qw), float_bits(b.qw))
+  }
+}
diff --git a/packages/ludic.physics/tests/physics_test.ludic b/packages/ludic.physics/tests/physics_test.ludic
new file mode 100644
index 00000000..5accf4fe
--- /dev/null
+++ b/packages/ludic.physics/tests/physics_test.ludic
@@ -0,0 +1,97 @@
+# physics_test.ludic - the real Jolt, one world a test: a ground, a box that falls onto it and a
+# ray that finds it; still things added, asked about and removed; a boulder's top; a push out of a
+# post; a crate that floats in the fake lake; hard contacts as facts; fixed steps; the same run twice
+import "ludic.physics"
+import "ludic.base"
+program PhysicsTest {
+  numbers float
+
+  # the lake: everything east of x = 20, its surface at y = 0, drifting north at 0.5 m/s
+  function lake_h(x: float, z: float) -> float {
+    if x > 20.0 { return 0.0 }
+    return PHYS_NONE
+  }
+  function lake_cz(x: float, z: float) -> float { return 0.5 }
+  bind PhysWater {
+    height: fn lake_h
+    current_z: fn lake_cz
+  }
+
+  # 16 x 16 samples, a metre apart, flat at y = 0 from (-8, -8)
+  function ground(physics_st: mut PhysicsState) -> void {
+    expect(phys_open(physics_st, 1024, 1))
+    let h = floats(256)
+    for i in 0 .. 256 { h[i] = 0.0 }
+    let g = phys_heightfield(physics_st, h, 16, -8.0, -8.0, 1.0)
+    expect(g >= 0)
+    expect(phys_ground_add(physics_st, g) >= 0)
+  }
+
+  function sim(physics_st: mut PhysicsState, seconds: float) -> void {
+    let n = int(seconds * 60.0)
+    for i in 0 .. n { phys_step(physics_st) }
+  }
+
+  test "a box falls onto the ground, rests there, and a ray finds its top" (physics_st: mut PhysicsState) {
+    ground(physics_st)
+    let b = phys_box(physics_st, 0.5, 0.5, 0.5)
+    let id = phys_body_add(physics_st, b, 0.0, 4.0, 0.0, 0.0, 20.0)
+    expect(id != -1)
+    sim(physics_st, 3.0)
+    let p = phys_pose(physics_st, id)
+    expect(p != null)
+    expect(Math.abs(p.y - 0.5) < 0.05)
+    expect(Math.abs(p.vy) < 0.05)
+    let h = phys_ray(physics_st, 0.0, 10.0, 0.0, 0.0, -20.0, 0.0, PHYS_M_ALL)
+    expect_eq(h.body, id)
+    expect(Math.abs(h.y - 1.0) < 0.05)
+    expect(h.ny > 0.9)
+    let g = phys_ray(physics_st, 5.0, 10.0, 5.0, 0.0, -20.0, 0.0, PHYS_M_ALL)
+    expect(g.body != -1)
+    expect(g.body != id)
+    expect(Math.abs(g.y) < 0.01)
+    phys_close(physics_st)
+  }
+
+  test "a still thing is found, then removed for good" (physics_st: mut PhysicsState) {
+    ground(physics_st)
+    let post = phys_offset(physics_st, phys_cylinder(physics_st, 1.0, 0.3), 0.0, 1.0, 0.0, 0.0)
+    let id = phys_static_add(physics_st, post, 2.0, 0.0, 2.0, 0.0)
+    phys_settle(physics_st)
+    let n = phys_count(physics_st)
+    let near = phys_overlap(physics_st, 2.0, 1.0, 2.0, 0.5, PHYS_M_STATIC)
+    expect_eq(len(near), 1)
+    expect_eq(near[0], id)
+    expect_eq(len(phys_overlap(physics_st, 2.0, 1.0, 2.0, 0.5, PHYS_M_GROUND)), 0)
+    phys_remove(physics_st, id)
+    expect_eq(phys_count(physics_st), n - 1)
+    expect_eq(len(phys_overlap(physics_st, 2.0, 1.0, 2.0, 0.5, PHYS_M_STATIC)), 0)
+    expect(phys_pose(physics_st, id) == null)
+    phys_close(physics_st)
+  }
+
+  test "a boulder's top is where a foot would land, and the ground can be left out" (physics_st: mut PhysicsState) {
+    ground(physics_st)
+    phys_static_add(physics_st, phys_dome(physics_st, 1.0, 0.5), 3.0, 0.0, 3.0, 0.0)
+    phys_settle(physics_st)
+    let top = phys_top_at(physics_st, 3.0, 3.0, 0.05, 5.0, 10.0, PHYS_M_STATIC)
+    expect(Math.abs(top - 0.5) < 0.06)
+    expect(phys_top_at(physics_st, -3.0, -3.0, 0.05, 5.0, 10.0, PHYS_M_STATIC) == PHYS_NONE)
+    expect(Math.abs(phys_top_at(physics_st, -3.0, -3.0, 0.05, 5.0, 10.0, PHYS_M_SOLID)) < 0.06)
+    phys_close(physics_st)
+  }
+
+  test "a body standing in a post is pushed out of it, sideways" (physics_st: mut PhysicsState) {
+    ground(physics_st)
+    let post = phys_offset(physics_st, phys_cylinder(physics_st, 1.0, 0.3), 0.0, 1.0, 0.0, 0.0)
+    phys_static_add(physics_st, post, 0.0, 0.0, 0.0, 0.0)
+    phys_settle(physics_st)
+    let p = phys_push(physics_st, 0.5, 0.0, 0.4, 0.1, 1.8, PHYS_M_STATIC)
+    expect_eq(p.n, 1)
+    expect(p.x > 0.1)
+    expect(Math.abs(p.z) < 0.05)
+    let clear = phys_push(physics_st, 3.0, 0.0, 0.4, 0.1, 1.8, PHYS_M_STATIC)
+    expect_eq(clear.n, 0)
+    phys_close(physics_st)
+  }
+}
diff --git a/packages/ludic.physics/water.ludic b/packages/ludic.physics/water.ludic
new file mode 100644
index 00000000..6679ee00
--- /dev/null
+++ b/packages/ludic.physics/water.ludic
@@ -0,0 +1,65 @@
+# water.ludic - a body that floats: each step it is buoyed against the surface PhysWater gives,
+# carried by its current, and slowed by its drags. Reaching the water is a PHYS_SPLASH fact.
+const PH_K: int = 4
+
+# buoyancy 1 floats level with the surface, more rides higher; drags in 1/s
+export function phys_float(physics_st: mut PhysicsState, id: int, buoyancy: float, linear_drag: float, angular_drag: float) -> void {
+  ph_buffers(physics_st)
+  phys_sink(physics_st, id)
+  push(physics_st.ph_floats, id)
+  push(physics_st.ph_float_k, buoyancy)
+  push(physics_st.ph_float_k, linear_drag)
+  push(physics_st.ph_float_k, angular_drag)
+  push(physics_st.ph_float_k, 0.0)
+}
+
+# no longer floats (it is still a body)
+export function phys_sink(physics_st: mut PhysicsState, id: int) -> void {
+  if physics_st.ph_floats == null { return }
+  let ids = new []int
+  let ks = new []float
+  for i in 0 .. len(physics_st.ph_floats) {
+    if physics_st.ph_floats[i] != id {
+      push(ids, physics_st.ph_floats[i])
+      for k in 0 .. PH_K { push(ks, physics_st.ph_float_k[i * PH_K + k]) }
+    }
+  }
+  physics_st.ph_floats = ids
+  physics_st.ph_float_k = ks
+}
+
+export function phys_floating(physics_st: PhysicsState, id: int) -> bool {
+  if physics_st.ph_floats == null { return false }
+  for i in 0 .. len(physics_st.ph_floats) { if physics_st.ph_floats[i] == id { return true } }
+  return false
+}
+
+function ph_buoy_all(physics_st: mut PhysicsState) -> void {
+  let dt = physics_st.ph_step
+  for i in 0 .. len(physics_st.ph_floats) {
+    let id = physics_st.ph_floats[i]
+    if jph_body_read(physics_st.ph_world, id, physics_st.ph_vals) == 0 {
+      ph_buoy_one(physics_st, i, id, physics_st.ph_vals[0], physics_st.ph_vals[2], dt)
+    }
+  }
+}
+
+function ph_buoy_one(physics_st: mut PhysicsState, i: int, id: int, x: float, z: float, dt: float) -> void {
+  let sy = PhysWater.height(x, z)
+  if sy <= PHYS_NONE { return }
+  let k = i * PH_K
+  let cx = PhysWater.current_x(x, z)
+  let cz = PhysWater.current_z(x, z)
+  let kk = physics_st.ph_float_k
+  let wet = jph_body_buoyancy(physics_st.ph_world, id, sy, 0.0, 1.0, 0.0, kk[k], kk[k + 1], kk[k + 2], cx, cz, dt)
+  if wet == 1 and kk[k + 3] == 0.0 { ph_splash(physics_st, id, x, sy, z) }
+  physics_st.ph_float_k[k + 3] = float(wet)
+}
+
+function ph_splash(physics_st: mut PhysicsState, id: int, x: float, y: float, z: float) -> void {
+  let f = new PhysFact
+  f.what = PHYS_SPLASH
+  f.a = id
+  f.x = x;  f.y = y;  f.z = z
+  q_push(phys_facts(physics_st), f)
+}