Drops, stacks and collides boxes and spheres under gravity and records where each one is at every frame, to be baked into keys.
How the precision is kept
Positions, quaternions and velocities are 64-bit numbers. The method is position-based dynamics with small fixed substeps (the file cites Müller and others, SCA 2020): the result depends on the frame rate, the number of substeps and the contact margins, not on the number type. The default gravity is 9.81 m/s².
Open item. Use 9.80665 as the default gravity, or say in the panel that 9.81 is a rounded value.
Calculations under this policy
12 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
The center, half sizes and turn of the box or sphere that stands for a placement.
scale = lengths of the matrix columns; half = max(1e-4, half box × scale); sphere radius = largest half
The scale is taken with Math.hypot per column and divided out of the rotation; a column of zero length is not guarded. Half sizes have an absolute floor of 1e-4.
anim/physics.js
lines 41 to 52
ModelCraft 0.8.1
/** A rigid body from a placement: its definition's bounds, its world matrix (whose scale goes into the size). */
export function bodyFromPlacement(id, W, bounds, { kind = 'active', shape = 'box', mass = 1, friction = 0.5, bounce = 0.1 } = {}) {
const sx = Math.hypot(W[0], W[1], W[2]), sy = Math.hypot(W[4], W[5], W[6]), sz = Math.hypot(W[8], W[9], W[10]);
const c = scl(add(bounds.min, bounds.max), 0.5), hl = scl(sub(bounds.max, bounds.min), 0.5);
const half = [Math.max(1e-4, hl[0] * sx), Math.max(1e-4, hl[1] * sy), Math.max(1e-4, hl[2] * sz)];
// the turn: W's columns without their scale (a mirrored placement is turned, not flipped)
const x = [W[0] / sx, W[1] / sx, W[2] / sx], y = [W[4] / sy, W[5] / sy, W[6] / sy], z = [W[8] / sz, W[9] / sz, W[10] / sz];
const q = matToQuat(x, y, z);
const centre = [W[0] * c[0] + W[4] * c[1] + W[8] * c[2] + W[12], W[1] * c[0] + W[5] * c[1] + W[9] * c[2] + W[13], W[2] * c[0] + W[6] * c[1] + W[10] * c[2] + W[14]];
return { id, kind, shape, half, r: Math.max(...half), x: centre, q, mass, friction, bounce, // a sphere: as big as the bounds' widest half
place: { scale: [sx, sy, sz], c } };
}
Rotation matrix to quaternion
64-bit
The unit quaternion of a rotation given by its three axes.
Picks one of four formulas by the largest diagonal term, so the square root it divides by stays well away from 0. The result is made unit afterward.
The placement's world matrix for a simulated position and turn.
axes = q applied to X, Y, Z, times the placement's scale; origin = body center − rotated, scaled center offset
The scale and the center offset measured when the body was made are put back, so a body that has not moved returns its matrix to rounding of the quaternion round trip.
anim/physics.js
lines 62 to 69
ModelCraft 0.8.1
/** The placement's world matrix with its body at state {x, q}. */
export function placementAt(body, st) {
const { scale: s, c } = body.place, q = st.q;
const X = rot(q, [1, 0, 0]), Y = rot(q, [0, 1, 0]), Z = rot(q, [0, 0, 1]);
const off = add(add(scl(X, s[0] * c[0]), scl(Y, s[1] * c[1])), scl(Z, s[2] * c[2]));
const t = sub(st.x, off);
return [X[0] * s[0], X[1] * s[0], X[2] * s[0], 0, Y[0] * s[1], Y[1] * s[1], Y[2] * s[1], 0, Z[0] * s[2], Z[1] * s[2], Z[2] * s[2], 0, t[0], t[1], t[2], 1];
}
Mass and inertia
64-bit
The inverse mass and inverse moments of inertia of a box or a sphere.
sphere: I = 0.4 m r²; box: I_x = m (b² + c²) ÷ 12, and likewise for y and z
Closed forms in 64-bit. The inverse mass is 1 ÷ mass with no guard for a mass of 0 in this function. 0.4 stands for 2 ÷ 5.
anim/physics.js
lines 71 to 86
ModelCraft 0.8.1
// ---- the solver ----
function prepare(b) {
const B = { ...b, x: b.x.slice(), q: qnorm(b.q.slice()), v: (b.v || [0, 0, 0]).slice(), w: (b.w || [0, 0, 0]).slice() };
B.dyn = b.kind === 'active';
B.invMass = B.dyn ? 1 / b.mass : 0;
if (b.shape === 'sphere') { B.r = b.r || Math.max(...b.half); const I = 0.4 * b.mass * B.r * B.r; B.invI = B.dyn ? [1 / I, 1 / I, 1 / I] : [0, 0, 0]; }
else {
const [a, c, d] = b.half.map((h) => 2 * h), m = b.mass;
B.invI = B.dyn ? [12 / (m * (c * c + d * d)), 12 / (m * (a * a + d * d)), 12 / (m * (a * a + c * c))] : [0, 0, 0];
// the points a box touches with: corners, edge middles, face middles (body coordinates)
B.pts = [];
for (const i of [-1, 0, 1]) for (const j of [-1, 0, 1]) for (const k of [-1, 0, 1]) if (i || j || k) B.pts.push([i * b.half[0], j * b.half[1], k * b.half[2]]);
}
B.bound = b.shape === 'sphere' ? B.r : len(b.half);
return B;
}
Depth of a point inside a box
64-bit
How far a point is inside a box and the outward normal of the nearest face.
local = q⁻¹ (p − center); depth = min over axes of (half + margin − |local|) − margin
The point is turned into the box's frame by the conjugate quaternion. The nearest face is chosen by comparing three depths exactly, so a point near an edge can switch faces from one substep to the next.
anim/physics.js
lines 89 to 97
ModelCraft 0.8.1
/** How far point p (world) is inside box B, and the normal out of B there (the nearest face). → {d, n} | null */
function inBox(B, p, margin = 0) {
const l = rotInv(B.q, sub(p, B.x)), h = B.half;
const ex = [h[0] + margin - Math.abs(l[0]), h[1] + margin - Math.abs(l[1]), h[2] + margin - Math.abs(l[2])];
if (ex[0] < 0 || ex[1] < 0 || ex[2] < 0) return null;
const a = ex[0] < ex[1] ? (ex[0] < ex[2] ? 0 : 2) : ex[1] < ex[2] ? 1 : 2;
const nl = [0, 0, 0]; nl[a] = l[a] >= 0 ? 1 : -1;
return { d: ex[a] - margin, n: rot(B.q, nl) };
}
Positional correction
64-bit
Pushes two bodies apart along a normal by a given distance, sharing it by their mass and inertia.
w = Σ (1 ÷ m + (r × n) · I⁻¹ (r × n)); λ = c ÷ w; x += n λ ÷ m; q += ½ (Δω, 0) q, made unit
The quaternion update is first order and is made unit again after each correction. Nothing is done when the combined inverse mass is 0 or less.
anim/physics.js
lines 137 to 154
ModelCraft 0.8.1
/** Move bodies apart along n by a positional correction of size c at points pa / pb (XPBD, compliance 0). → λ */
function applyPositional(A, B, pa, pb, n, c) {
const ra = sub(pa, A.x), rb = B ? sub(pb, B.x) : null;
const ca = cross(ra, n), cb = rb ? cross(rb, n) : null;
const wa = A.invMass + dot(ca, invIw(A, ca)), wb = B ? B.invMass + dot(cb, invIw(B, cb)) : 0;
const w = wa + wb;
if (w <= 0) return 0;
const lam = c / w, P = scl(n, lam);
const nudge = (X, r, Pp, s) => {
if (!X.dyn) return;
X.x = add(X.x, scl(Pp, s * X.invMass));
const dw = invIw(X, cross(r, scl(Pp, s)));
X.q = qnorm(add4(X.q, scl4(qmul([dw[0], dw[1], dw[2], 0], X.q), 0.5)));
};
nudge(A, ra, P, 1);
if (B) nudge(B, rb, P, -1);
return lam;
}
Step size and frame count
64-bit
How many frames are recorded and how long one substep is.
The frame count is rounded to a whole number. The step is two divisions, so for most frame rates it is not an exact binary fraction and time advances by a rounded step.
Pushes each touching point out of the surface and holds it still where friction can.
μ = √(μ_a μ_b); held when slide < μ |λ| × 1.5 + 1e-9
Each contact is measured again just before it is solved. The slide along the surface is a difference of the contact point's position before and after the substep. The factor 1.5 and the 1e-9 margin are fixed in the code.
anim/physics.js
lines 214 to 230
ModelCraft 0.8.1
// positions: push out, then hold still where friction can (static friction)
for (const c of cs) {
const A = c.a, O = c.b, g = c.at();
if (!g) continue;
c.p = g.p; c.n = g.n;
const pa = g.p, pb = g.p;
const lam = applyPositional(A, O, pa, pb, g.n, g.d);
c.lam = lam;
// how far the contact point slid along the surface during this substep
const ra = rotInv(A.pq || A.q, sub(c.p, A.px || A.x));
const pPrev = add(A.px || A.x, rot(A.pq || A.q, ra)), pNow = c.p;
let slide = sub(pNow, pPrev);
if (O && O.dyn) { const rb = rotInv(O.pq || O.q, sub(c.p, O.px || O.x)); slide = sub(slide, sub(add(O.x, rot(O.q, rb)), add(O.px || O.x, rot(O.pq || O.q, rb)))); }
const st = sub(slide, scl(c.n, dot(slide, c.n))), stl = len(st);
const mu = O ? Math.sqrt(A.friction * O.friction) : A.friction;
if (stl > 1e-12 && stl < mu * Math.abs(lam) * 1.5 + 1e-9) applyPositional(A, O, pa, pb, scl(st, -1 / stl), stl);
}
Velocities from positions
64-bit
Works each body's speed and spin back from how far it moved in the substep.
v = (x − x_prev) ÷ h; ω = 2 ÷ h × vector part of (q × q_prev⁻¹)
The speed is a difference of two positions divided by a small step, so it loses digits when the body barely moves; the settled test allows for that with fixed thresholds.
Slows sliding at each contact and sets the speed a body leaves a surface with.
friction impulse = min(μ |λ| ÷ h, sliding speed ÷ inverse mass); leaving speed = −e × min(incoming normal speed, 0); e = 0 when |incoming| < 2 |g · n| h
The normal speeds are corrected in 4 fixed sweeps over the contacts, each skipped when the difference is 1e-9 or less. The number of sweeps is fixed, so stacked bodies are not solved to a tolerance.
anim/physics.js
lines 238 to 254
ModelCraft 0.8.1
// dynamic friction, then the normal speeds (a few sweeps: one corner's correction changes the others')
const live = cs.filter((c) => c.lam && c.n);
for (const c of live) {
const A = c.a, O = c.b, vr = sub(velAt(A, c.p), velAt(O, c.p)), vn = dot(vr, c.n);
const vt = sub(vr, scl(c.n, vn)), vtl = len(vt);
// dynamic friction as an impulse: μ × this contact's push (λ / h), never more than stops the sliding
const mu = O ? Math.sqrt(A.friction * O.friction) : A.friction;
if (vtl > 1e-9) { const dir = scl(vt, -1 / vtl); applyImpulse(A, O, c.p, dir, Math.min((mu * Math.abs(c.lam)) / h, vtl / genMass(A, O, c.p, dir))); }
const vPrevN = dot(sub(add(A.vPrev || A.v, cross(A.wPrev || A.w, sub(c.p, A.x))), O ? add(O.vPrev || O.v, cross(O.wPrev || O.w, sub(c.p, O.x))) : [0, 0, 0]), c.n);
const e = Math.abs(vPrevN) < 2 * Math.abs(dot(gravity, c.n)) * h ? 0 : (O ? Math.max(A.bounce, O.bounce) : A.bounce);
// the normal speed after the contact: −e × the speed it came in with (the push out of the surface gives none of its own)
c.target = -e * Math.min(vPrevN, 0);
}
for (let sweep = 0; sweep < 4; sweep++) for (const c of live) {
const vn = dot(sub(velAt(c.a, c.p), velAt(c.b, c.p)), c.n);
if (Math.abs(c.target - vn) > 1e-9) applyVelocity(c.a, c.b, c.p, scl(c.n, c.target - vn));
}
Settled test
64-bit
Says whether every body has come to rest at the end of the run.
settled when every |v| < 0.02 and every |ω| < 0.05
Fixed absolute thresholds on the last substep's speeds; they do not scale with the size of the bodies.
Extended position-based dynamics for rigid bodies: explicit integration per substep, contacts solved as position corrections with static friction, then restitution and dynamic friction on the velocities in 4 sweeps.
Number types
64-bit numbers in plain arrays.
Tolerances
Contact search margin 1e-3 (absolute, meters); sizes floored at 1e-4; normal-speed correction skipped under 1e-9; a run counts as settled when every speed is under 0.02 and every spin under 0.05.
Checked by
bin/test-animation.mjs loads the development build's copy of this file and calls simulate, bodyFromPlacement and placementAt.
Limits
A box touches with 26 points (corners, edge middles, face middles), so box against box is approximate. Only boxes and spheres. The step is 1 ÷ frames per second ÷ substeps (20 by default). The default gravity 9.81 is a rounded value of standard gravity (9.80665).
Approach: 64-bit, stated limit
The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.