fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 1// Time integration of the bead system.
2//
3// The conservative part (springs, external force pulses) is integrated with
4// velocity Verlet, optionally composed into Yoshida's fourth-order scheme.
5// Damping (friction and dashpots) is linear in the velocities, dv/dt = -A v,
6// and is applied as a separate exact step v <- exp(-tau A) v, once per audio
7// sample (Strang splitting over the sample interval Delta = 1/fs). The
8// splitting error is second order in Delta and scales with the damping rate,
9// which is small compared with the oscillation frequencies; the conservative
10// part is second order (Verlet) or fourth order (Yoshida) in the step h.
12import type { BeadSystem, Method, Modes, SimResult } from "./types";
14const YOSHIDA_W1 = 1 / (2 - Math.cbrt(2));
15const YOSHIDA_W0 = -Math.cbrt(2) * YOSHIDA_W1;
17export function stepCoefficients(method: Method): number[] {
18 return method === "yoshida4" ? [YOSHIDA_W1, YOSHIDA_W0, YOSHIDA_W1] : [1];
19}
21// Largest h*omega for which the method is stable on x'' = -omega^2 x.
22export function stabilityLimit(method: Method): number {
23 if (method === "verlet") return 2;
24 const cs = stepCoefficients(method);
25 const stable = (h: number) => {
26 // transfer matrix of one full step for omega = 1
27 let m = [1, 0, 0, 1]; // [x<-x, x<-v, v<-x, v<-v]
28 const mul = (a: number[], b: number[]) => [
29 a[0] * b[0] + a[1] * b[2],
30 a[0] * b[1] + a[1] * b[3],
31 a[2] * b[0] + a[3] * b[2],
32 a[2] * b[1] + a[3] * b[3],
33 ];
34 for (const c of cs) {
35 const kick = [1, 0, (-c * h) / 2, 1];
36 const drift = [1, c * h, 0, 1];
37 m = mul(kick, mul(drift, mul(kick, m)));
38 }
39 return Math.abs(m[0] + m[3]) < 2;
40 };
41 let h = 1e-3;
42 while (stable(h + 1e-3)) h += 1e-3;
43 return h;
44}
46// Error of the numerical oscillation frequency for x'' = -omega^2 x, in cents
47// (positive: the simulated tone is sharp).
48export function frequencyErrorCents(method: Method, omega: number, h: number): number {
49 let m = [1, 0, 0, 1];
50 const mul = (a: number[], b: number[]) => [
51 a[0] * b[0] + a[1] * b[2],
52 a[0] * b[1] + a[1] * b[3],
53 a[2] * b[0] + a[3] * b[2],
54 a[2] * b[1] + a[3] * b[3],
55 ];
56 for (const c of stepCoefficients(method)) {
57 const kick = [1, 0, (-c * h * omega * omega) / 2, 1];
58 const drift = [1, c * h, 0, 1];
59 m = mul(kick, mul(drift, mul(kick, m)));
60 }
61 const half = (m[0] + m[3]) / 2;
62 if (Math.abs(half) >= 1) return NaN;
63 const omegaNum = Math.acos(half) / h;
64 return 1200 * Math.log2(omegaNum / omega);
65}
67export interface SimOptions {
68 modes?: Modes | null;
69 onProgress?: (t: number) => void;
95f2216Minimal chaos example and perturbed-twin chaos test with separation plot and growth-rate fitJeremy Magland 70 // called at every frame record with the full-precision state (do not keep the arrays)
71 onFrame?: (frame: number, t: number, x: Float64Array, v: Float64Array) => void;
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 72 maxAnimBytes?: number;
73}
75export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
76 const t0 = performance.now();
77 const warnings: string[] = [];
78 const { n, mass, pinned } = sys;
79 const N3 = 3 * n;
80 const x = new Float64Array(sys.pos0);
81 const v = new Float64Array(N3);
82 const F = new Float64Array(N3);
83 const invm = new Float64Array(n);
84 for (let i = 0; i < n; i++) invm[i] = pinned[i] ? 0 : 1 / mass[i];
86 // ---- forces ------------------------------------------------------------
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 87 const { si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb } = sys;
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 88 const ns = si.length;
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 89 const nb = bi.length;
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 90 const internalForces = (): number => {
91 F.fill(0);
92 let pe = 0;
93 for (let s = 0; s < ns; s++) {
94 const i3 = 3 * si[s];
95 const j3 = 3 * sj[s];
96 const dx = x[j3] - x[i3];
97 const dy = x[j3 + 1] - x[i3 + 1];
98 const dz = x[j3 + 2] - x[i3 + 2];
99 const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
100 const e = L - sL0[s];
101 const k = sk[s];
102 const k3 = sk3[s];
103 const e2 = e * e;
104 pe += 0.5 * k * e2 + 0.25 * k3 * e2 * e2;
105 const g = (k * e + k3 * e2 * e) / L;
106 F[i3] += g * dx;
107 F[i3 + 1] += g * dy;
108 F[i3 + 2] += g * dz;
109 F[j3] -= g * dx;
110 F[j3 + 1] -= g * dy;
111 F[j3 + 2] -= g * dz;
112 }
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 113 for (let b = 0; b < nb; b++) {
114 const i3 = 3 * bi[b];
115 const j3 = 3 * bj[b];
116 const k3 = 3 * bk[b];
117 const kb = bkb[b];
118 const d0 = x[i3] - 2 * x[j3] + x[k3];
119 const d1 = x[i3 + 1] - 2 * x[j3 + 1] + x[k3 + 1];
120 const d2 = x[i3 + 2] - 2 * x[j3 + 2] + x[k3 + 2];
121 pe += 0.5 * kb * (d0 * d0 + d1 * d1 + d2 * d2);
122 F[i3] -= kb * d0;
123 F[i3 + 1] -= kb * d1;
124 F[i3 + 2] -= kb * d2;
125 F[j3] += 2 * kb * d0;
126 F[j3 + 1] += 2 * kb * d1;
127 F[j3 + 2] += 2 * kb * d2;
128 F[k3] -= kb * d0;
129 F[k3 + 1] -= kb * d1;
130 F[k3 + 2] -= kb * d2;
131 }
133 };
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 134 const pulses = sys.hits.filter(h => h.duration > 0);
135 const externalForces = (t: number) => {
136 for (const p of pulses) {
137 const u = (t - p.t) / p.duration;
138 if (u <= 0 || u >= 1) continue;
139 const w = 0.5 * (1 - Math.cos(2 * Math.PI * u));
140 for (const b of p.beads) {
141 F[3 * b] += w * p.force[0];
142 F[3 * b + 1] += w * p.force[1];
143 F[3 * b + 2] += w * p.force[2];
144 }
145 }
146 };
148 // ---- rest state checks ---------------------------------------------------
149 let pe = internalForces();
150 const peRest = pe; // energy stored by pre-tension; energies are reported relative to rest
151 let residualForce = 0;
152 for (let i = 0; i < n; i++) {
153 if (pinned[i]) continue;
154 const f = Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]);
155 if (f > residualForce) residualForce = f;
156 }
157 let maxTension = 0;
158 for (let s = 0; s < ns; s++) {
159 const i3 = 3 * si[s];
160 const j3 = 3 * sj[s];
161 const L = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
162 const e = L - sL0[s];
163 maxTension = Math.max(maxTension, Math.abs(sk[s] * e + sk3[s] * e * e * e));
164 }
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 165 // force scale: the larger of the rest tension and the force of a 1e-3 strain
166 // on any spring or bending triplet
167 let forceScale = maxTension;
168 for (let s = 0; s < ns; s++) forceScale = Math.max(forceScale, 1e-3 * sk[s] * sL0[s]);
169 for (let b = 0; b < nb; b++) {
170 const i3 = 3 * bi[b];
171 const j3 = 3 * bj[b];
172 const len = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
173 forceScale = Math.max(forceScale, 1e-3 * bkb[b] * len);
174 }
175 if (residualForce > 1e-6 * Math.max(forceScale, 1e-30)) {
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 176 warnings.push(
177 `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 178 `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 179 );
180 }
182 // ---- time step -----------------------------------------------------------
183 const omegaMax = estimateOmegaMax(sys);
184 const fs = sys.sampleRate;
185 const hLimit = stabilityLimit(sys.method) / omegaMax;
186 let substeps: number;
187 if (sys.dt > 0) {
188 substeps = Math.max(1, Math.round(1 / (fs * sys.dt)));
189 if (1 / (fs * substeps) > hLimit) {
190 warnings.push(
191 `dt = ${(1 / (fs * substeps)).toExponential(3)} s exceeds the stability limit ${hLimit.toExponential(3)} s.`
192 );
193 }
194 } else {
195 substeps = Math.max(1, Math.ceil(1 / (fs * sys.safety * hLimit)));
196 }
197 const h = 1 / (fs * substeps);
198 const nSamples = Math.max(1, Math.round(sys.duration * fs));
199 const nSteps = nSamples * substeps;
200 const cs = stepCoefficients(sys.method);
202 // ---- damping: v <- exp(-tau A) v ------------------------------------------
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 203 // dashpots with c = 0 (e.g. from bp_decay on untensioned springs) are dropped
204 const keepD = Array.from(sys.dc.keys()).filter(d => sys.dc[d] !== 0);
205 const di = Int32Array.from(keepD, d => sys.di[d]);
206 const dj = Int32Array.from(keepD, d => sys.dj[d]);
207 const dc = Float64Array.from(keepD, d => sys.dc[d]);
208 const { friction } = sys;
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 209 const nd = di.length;
210 let normA = 0; // infinity-norm bound of A
211 {
212 const row = new Float64Array(n);
213 for (let i = 0; i < n; i++) row[i] = friction[i];
214 for (let d = 0; d < nd; d++) {
215 row[di[d]] += 2 * dc[d];
216 row[dj[d]] += 2 * dc[d];
217 }
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 218 for (let b = 0; b < nb; b++) {
219 row[bi[b]] += 4 * bcb[b];
220 row[bj[b]] += 8 * bcb[b];
221 row[bk[b]] += 4 * bcb[b];
222 }
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 223 for (let i = 0; i < n; i++) normA = Math.max(normA, row[i] * invm[i]);
224 }
225 const hasDamping = normA > 0;
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 226 const frictionOnly = dc.every(c => c === 0) && bcb.every(c => c === 0); // then exp(-tau A) is diagonal
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 227 const term = new Float64Array(N3);
228 const tmp = new Float64Array(N3);
229 const fa = new Float64Array(n);
230 for (let i = 0; i < n; i++) fa[i] = friction[i] * invm[i];
231 const dai = new Float64Array(nd);
232 const daj = new Float64Array(nd);
233 for (let d = 0; d < nd; d++) {
234 dai[d] = dc[d] * invm[di[d]];
235 daj[d] = dc[d] * invm[dj[d]];
236 }
237 const applyA = (w: Float64Array, out: Float64Array) => {
238 for (let i = 0; i < n; i++) {
239 const a = fa[i];
240 out[3 * i] = a * w[3 * i];
241 out[3 * i + 1] = a * w[3 * i + 1];
242 out[3 * i + 2] = a * w[3 * i + 2];
243 }
244 for (let d = 0; d < nd; d++) {
245 const i3 = 3 * di[d];
246 const j3 = 3 * dj[d];
247 const ai = dai[d];
248 const aj = daj[d];
249 const r0 = w[i3] - w[j3];
250 const r1 = w[i3 + 1] - w[j3 + 1];
251 const r2 = w[i3 + 2] - w[j3 + 2];
252 out[i3] += ai * r0;
253 out[i3 + 1] += ai * r1;
254 out[i3 + 2] += ai * r2;
255 out[j3] -= aj * r0;
256 out[j3 + 1] -= aj * r1;
257 out[j3 + 2] -= aj * r2;
258 }
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 259 for (let b = 0; b < nb; b++) {
260 const cb = bcb[b];
261 if (cb === 0) continue;
262 const i3 = 3 * bi[b];
263 const j3 = 3 * bj[b];
264 const k3 = 3 * bk[b];
265 const ci = cb * invm[bi[b]];
266 const cj = -2 * cb * invm[bj[b]];
267 const ck = cb * invm[bk[b]];
268 const r0 = w[i3] - 2 * w[j3] + w[k3];
269 const r1 = w[i3 + 1] - 2 * w[j3 + 1] + w[k3 + 1];
270 const r2 = w[i3 + 2] - 2 * w[j3 + 2] + w[k3 + 2];
271 out[i3] += ci * r0;
272 out[i3 + 1] += ci * r1;
273 out[i3 + 2] += ci * r2;
274 out[j3] += cj * r0;
275 out[j3 + 1] += cj * r1;
276 out[j3 + 2] += cj * r2;
277 out[k3] += ck * r0;
278 out[k3 + 1] += ck * r1;
279 out[k3 + 2] += ck * r2;
280 }
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 281 };
282 const kinetic = () => {
283 let ke = 0;
284 for (let i = 0; i < n; i++) {
285 ke += 0.5 * mass[i] * (v[3 * i] ** 2 + v[3 * i + 1] ** 2 + v[3 * i + 2] ** 2);
286 }
287 return ke;
288 };
289 let dissipated = 0;
290 const applyDamping = (tau: number) => {
291 if (!hasDamping || tau === 0) return;
292 const keBefore = kinetic();
293 if (frictionOnly) {
294 for (let i = 0; i < n; i++) {
295 const f = Math.exp(-tau * friction[i] * invm[i]);
296 v[3 * i] *= f;
297 v[3 * i + 1] *= f;
298 v[3 * i + 2] *= f;
299 }
300 } else {
301 // Taylor series with scaling; |tau_s| * normA <= 0.5 so terms shrink fast.
302 const nsub = Math.max(1, Math.ceil((Math.abs(tau) * normA) / 0.5));
303 const ts = tau / nsub;
304 for (let r = 0; r < nsub; r++) {
305 let vmax = 0;
306 for (let q = 0; q < N3; q++) {
307 term[q] = v[q];
308 vmax = Math.max(vmax, Math.abs(v[q]));
309 }
310 if (vmax === 0) break;
311 for (let k = 1; k < 40; k++) {
312 applyA(term, tmp);
313 const c = -ts / k;
314 let tmax = 0;
315 for (let q = 0; q < N3; q++) {
316 const val = c * tmp[q];
317 term[q] = val;
318 v[q] += val;
319 const a = Math.abs(val);
320 if (a > tmax) tmax = a;
321 }
322 if (tmax <= 1e-17 * vmax) break;
323 }
324 }
325 }
326 dissipated += keBefore - kinetic();
327 };
329 // ---- recording -------------------------------------------------------------
330 const audio = new Float32Array(nSamples);
331 const lb = sys.listenBeads;
332 const [ldx, ldy, ldz] = sys.listenDir;
333 const pos0 = sys.pos0;
452405aListener quantity: acceleration (default), velocity or displacement; acceleration from forces with boundary-corrected damping; zero-mode threshold 1e-6 of top frequencyJeremy Magland 334 // Acceleration is (F - damping force) / m on each listener bead, with F the
335 // spring, bending and pulse forces at the current positions and time (the
336 // force array is current at every sample boundary). The damping force uses
337 // the dashpots and bending triplets that touch each listener bead.
338 const touchD = Array.from(lb, b => Array.from(di.keys()).filter(d => di[d] === b || dj[d] === b));
339 const touchB = Array.from(lb, b => Array.from(bi.keys()).filter(t => (bcb[t] !== 0 && (bi[t] === b || bj[t] === b || bk[t] === b))));
340 const dampForce = (v: Float64Array, q: number, b: number, c: number) => {
341 let f = -friction[b] * v[3 * b + c];
342 for (const d of touchD[q]) {
343 const other = di[d] === b ? dj[d] : di[d];
344 f -= dc[d] * (v[3 * b + c] - v[3 * other + c]);
345 }
346 for (const t of touchB[q]) {
347 const r = v[3 * bi[t] + c] - 2 * v[3 * bj[t] + c] + v[3 * bk[t] + c];
348 const w = bj[t] === b ? -2 : 1;
349 f -= bcb[t] * w * r;
350 }
351 return f;
352 };
353 // At a sample boundary the damping for the time pendingD has not been
354 // applied yet; vb = v - pendingD A v brings the velocities to the boundary
355 // (to second order) for the listener's velocity and damping force.
356 const vb = new Float64Array(N3);
357 const Av = new Float64Array(N3);
358 const listen = (pendingD: number) => {
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 359 let s = 0;
452405aListener quantity: acceleration (default), velocity or displacement; acceleration from forces with boundary-corrected damping; zero-mode threshold 1e-6 of top frequencyJeremy Magland 360 if (sys.listenQuantity !== "displacement") {
361 if (pendingD !== 0 && hasDamping) {
362 applyA(v, Av);
363 for (let q = 0; q < N3; q++) vb[q] = v[q] - pendingD * Av[q];
364 } else vb.set(v);
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 365 }
452405aListener quantity: acceleration (default), velocity or displacement; acceleration from forces with boundary-corrected damping; zero-mode threshold 1e-6 of top frequencyJeremy Magland 366 lb.forEach((b, q) => {
367 const dir = [ldx, ldy, ldz];
368 for (let c = 0; c < 3; c++) {
369 if (dir[c] === 0) continue;
370 let val: number;
371 if (sys.listenQuantity === "displacement") val = x[3 * b + c] - pos0[3 * b + c];
372 else {
373 if (sys.listenQuantity === "velocity") val = vb[3 * b + c];
374 else val = (F[3 * b + c] + dampForce(vb, q, b, c)) * invm[b];
375 }
376 s += val * dir[c];
377 }
378 });
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 379 return s / lb.length;
380 };
381 const frameStride = Math.max(1, Math.round(fs / sys.frameRate));
382 const nFrames = Math.floor((nSamples - 1) / frameStride) + 1;
383 const maxAnimBytes = opts.maxAnimBytes ?? 150e6;
384 const animStride = frameStride * Math.max(1, Math.ceil((nFrames * N3 * 4) / maxAnimBytes));
385 const frameTimes = new Float64Array(nFrames);
386 const energy = new Float64Array(nFrames);
387 const dissipatedRec = new Float64Array(nFrames);
388 const animTimes: number[] = [];
389 const animFrames: Float32Array[] = [];
391 const modes = opts.modes ?? null;
392 const nModes = modes ? modes.freqs.length : 0;
393 const modeEnergy = modes ? new Float32Array(nFrames * nModes) : null;
394 let bz: Float64Array | null = null;
395 let bv: Float64Array | null = null;
396 let om2: Float64Array | null = null;
397 if (modes) {
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 398 bz = new Float64Array(modes.free.length);
399 bv = new Float64Array(modes.free.length);
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 400 om2 = new Float64Array(nModes);
401 for (let k = 0; k < nModes; k++) om2[k] = (2 * Math.PI * modes.freqs[k]) ** 2;
402 }
403 const recordModes = (f: number) => {
404 if (!modes || !modeEnergy || !bz || !bv || !om2) return;
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 405 const free = modes.free;
406 const nf = free.length;
407 for (let r = 0; r < nf; r++) {
408 const b = free[r];
409 const sm = Math.sqrt(mass[b]);
410 bz[r] = sm * (x[3 * b + 2] - pos0[3 * b + 2]);
411 bv[r] = sm * v[3 * b + 2];
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 412 }
413 const psi = modes.psi;
414 for (let k = 0; k < nModes; k++) {
415 let q = 0;
416 let qd = 0;
418 for (let r = 0; r < nf; r++) {
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 419 q += psi[off + r] * bz[r];
420 qd += psi[off + r] * bv[r];
421 }
422 modeEnergy[f * nModes + k] = 0.5 * (qd * qd + om2[k] * q * q);
423 }
424 };
426 // ---- impulses, by audio sample ----------------------------------------------
427 const impulsesAt = new Map<number, typeof sys.hits>();
428 for (const hit of sys.hits) {
429 if (hit.duration > 0) continue;
430 const s = Math.round(hit.t * fs);
431 if (s >= nSamples) continue;
432 if (!impulsesAt.has(s)) impulsesAt.set(s, []);
433 impulsesAt.get(s)!.push(hit);
434 }
435 const tExcitationEnd = Math.max(0, ...sys.hits.map(hh => hh.t + hh.duration));
437 // ---- main loop ---------------------------------------------------------------
438 // Each audio sample interval Delta = substeps * h is one macro step:
439 // D(Delta/2) [substeps conservative steps] D(Delta/2). The damping halves of
440 // neighboring macro steps are merged unless something needs the velocities
441 // at the boundary (a frame record or an impulse).
442 const Delta = 1 / fs;
443 pe = internalForces();
444 externalForces(0);
445 let pendingD = 0;
446 let aborted: string | null = null;
447 let frame = 0;
448 let lastProgress = performance.now();
449 for (let s = 0; s <= nSteps; s++) {
450 const t = s * h;
451 const sample = s % substeps === 0 ? s / substeps : -1;
452 const isFrame = sample >= 0 && sample % frameStride === 0 && frame < nFrames;
453 const imps = sample >= 0 ? impulsesAt.get(sample) : undefined;
454 if ((isFrame || imps) && pendingD !== 0) {
455 applyDamping(pendingD);
456 pendingD = 0;
457 }
458 if (imps) {
459 for (const hit of imps) {
460 for (const b of hit.beads) {
461 if (pinned[b]) continue;
462 v[3 * b] += hit.velocity[0];
463 v[3 * b + 1] += hit.velocity[1];
464 v[3 * b + 2] += hit.velocity[2];
465 }
466 }
467 }
452405aListener quantity: acceleration (default), velocity or displacement; acceleration from forces with boundary-corrected damping; zero-mode threshold 1e-6 of top frequencyJeremy Magland 468 if (sample >= 0 && sample < nSamples) audio[sample] = listen(pendingD);
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 469 if (isFrame) {
470 frameTimes[frame] = t;
471 energy[frame] = kinetic() + pe - peRest;
472 dissipatedRec[frame] = dissipated;
473 recordModes(frame);
95f2216Minimal chaos example and perturbed-twin chaos test with separation plot and growth-rate fitJeremy Magland 474 opts.onFrame?.(frame, t, x, v);
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 475 if (sample % animStride === 0) {
476 const d = new Float32Array(N3);
477 for (let q = 0; q < N3; q++) d[q] = x[q] - pos0[q];
478 animTimes.push(t);
479 animFrames.push(d);
480 }
481 frame++;
482 }
483 if (s === nSteps) break;
484 if (sample >= 0) {
485 pendingD += Delta / 2;
486 applyDamping(pendingD);
487 pendingD = 0;
488 }
490 let tj = t;
491 for (const c of cs) {
492 const kh = (c * h) / 2;
493 for (let i = 0; i < n; i++) {
494 const im = invm[i] * kh;
495 v[3 * i] += im * F[3 * i];
496 v[3 * i + 1] += im * F[3 * i + 1];
497 v[3 * i + 2] += im * F[3 * i + 2];
498 }
499 const dh = c * h;
500 for (let q = 0; q < N3; q++) x[q] += dh * v[q];
501 tj += c * h;
502 pe = internalForces();
503 externalForces(tj);
504 for (let i = 0; i < n; i++) {
505 const im = invm[i] * kh;
506 v[3 * i] += im * F[3 * i];
507 v[3 * i + 1] += im * F[3 * i + 1];
508 v[3 * i + 2] += im * F[3 * i + 2];
509 }
510 }
511 if ((s + 1) % substeps === 0) pendingD = Delta / 2;
513 if ((s & 1023) === 0) {
514 if (!Number.isFinite(pe)) {
515 aborted = `The solution blew up at t = ${t.toFixed(4)} s. Reduce dt (or safety) or the strike strength.`;
516 break;
517 }
518 const now = performance.now();
519 if (opts.onProgress && now - lastProgress > 200) {
520 opts.onProgress(t);
521 lastProgress = now;
522 }
523 }
524 }
526 // energy drift after the last excitation
527 let energyDrift = 0;
528 let ref = -1;
529 for (let f = 0; f < frame; f++) {
530 if (frameTimes[f] < tExcitationEnd) continue;
531 const total = energy[f] + dissipatedRec[f];
532 if (ref < 0) ref = total;
533 else if (ref > 0) energyDrift = Math.max(energyDrift, Math.abs(total - ref) / ref);
534 }
536 return {
537 audio,
538 sampleRate: fs,
539 dt: h,
540 substeps,
541 omegaMax,
542 frameTimes: frameTimes.slice(0, frame),
543 energy: energy.slice(0, frame),
544 dissipated: dissipatedRec.slice(0, frame),
545 animTimes: Float64Array.from(animTimes),
546 animFrames,
547 modeEnergy: modeEnergy ? modeEnergy.slice(0, frame * nModes) : null,
548 residualForce,
549 maxTension,
550 energyDrift,
551 freqErrors: [1000, 2000, 5000, 10000, 20000]
552 .filter(hz => hz < fs / 2)
553 .map(hz => ({ hz, cents: frequencyErrorCents(sys.method, 2 * Math.PI * hz, h) })),
554 elapsedMs: performance.now() - t0,
555 aborted,
556 warnings,
557 };
558}
560// Highest frequency (rad/s) of the system linearized about the rest positions,
561// by power iteration on M^{-1/2} H M^{-1/2}, with a small margin.
562export function estimateOmegaMax(sys: BeadSystem): number {
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 563 const { n, si, sj, sk, sL0, sk3, pos0, mass, pinned } = sys;
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 564 const ns = si.length;
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 565 // per-spring 3x3 blocks B = keff d d^T + (f/L)(I - d d^T)
566 const B = new Float64Array(ns * 6);
567 for (let s = 0; s < ns; s++) {
568 const i3 = 3 * si[s];
569 const j3 = 3 * sj[s];
570 const dx = pos0[j3] - pos0[i3];
571 const dy = pos0[j3 + 1] - pos0[i3 + 1];
572 const dz = pos0[j3 + 2] - pos0[i3 + 2];
573 const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
574 const e = L - sL0[s];
575 const keff = sk[s] + 3 * sk3[s] * e * e;
576 const g = (sk[s] * e + sk3[s] * e * e * e) / L;
577 const u = [dx / L, dy / L, dz / L];
578 const a = keff - g;
579 B[6 * s] = a * u[0] * u[0] + g;
580 B[6 * s + 1] = a * u[1] * u[1] + g;
581 B[6 * s + 2] = a * u[2] * u[2] + g;
582 B[6 * s + 3] = a * u[0] * u[1];
583 B[6 * s + 4] = a * u[0] * u[2];
584 B[6 * s + 5] = a * u[1] * u[2];
585 }
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 586 const ism = new Float64Array(n);
587 for (let i = 0; i < n; i++) ism[i] = pinned[i] ? 0 : 1 / Math.sqrt(mass[i]);
588 let w = new Float64Array(3 * n);
589 let seed = 12345;
590 for (let q = 0; q < 3 * n; q++) {
591 seed = (seed * 1103515245 + 12345) & 0x7fffffff;
592 w[q] = ism[Math.floor(q / 3)] * (seed / 0x7fffffff - 0.5);
593 }
594 const u = new Float64Array(3 * n);
595 const y = new Float64Array(3 * n);
596 let lambda = 0;
597 for (let it = 0; it < 500; it++) {
598 let norm = 0;
599 for (let q = 0; q < 3 * n; q++) norm += w[q] * w[q];
600 norm = Math.sqrt(norm);
601 if (norm === 0) return 1;
602 for (let q = 0; q < 3 * n; q++) u[q] = (w[q] / norm) * ism[Math.floor(q / 3)];
603 y.fill(0);
604 for (let s = 0; s < ns; s++) {
605 const i3 = 3 * si[s];
606 const j3 = 3 * sj[s];
607 const r0 = u[i3] - u[j3];
608 const r1 = u[i3 + 1] - u[j3 + 1];
609 const r2 = u[i3 + 2] - u[j3 + 2];
610 const b = 6 * s;
611 const f0 = B[b] * r0 + B[b + 3] * r1 + B[b + 4] * r2;
612 const f1 = B[b + 3] * r0 + B[b + 1] * r1 + B[b + 5] * r2;
613 const f2 = B[b + 4] * r0 + B[b + 5] * r1 + B[b + 2] * r2;
614 y[i3] += f0;
615 y[i3 + 1] += f1;
616 y[i3 + 2] += f2;
617 y[j3] -= f0;
618 y[j3 + 1] -= f1;
619 y[j3 + 2] -= f2;
620 }
74b1dcdFree edges: bending triplets (sys.bends) with damping, lattice triplets from bp_hexdisk/bp_chain, bp_decay for bending, bp_anchor soft supports, zero-mode warning; plate and bar examplesJeremy Magland 621 for (let b = 0; b < sys.bi.length; b++) {
622 const i3 = 3 * sys.bi[b];
623 const j3 = 3 * sys.bj[b];
624 const k3 = 3 * sys.bk[b];
625 const kb = sys.bkb[b];
626 for (let q = 0; q < 3; q++) {
627 const r = kb * (u[i3 + q] - 2 * u[j3 + q] + u[k3 + q]);
628 y[i3 + q] += r;
629 y[j3 + q] -= 2 * r;
630 y[k3 + q] += r;
631 }
632 }
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 633 let rq = 0;
634 const wNext = new Float64Array(3 * n);
635 for (let q = 0; q < 3 * n; q++) {
636 wNext[q] = y[q] * ism[Math.floor(q / 3)];
637 rq += (w[q] / norm) * wNext[q];
638 }
639 w = wNext;
640 if (it > 20 && Math.abs(rq - lambda) < 1e-8 * Math.abs(rq)) {
641 lambda = rq;
642 break;
643 }
644 lambda = rq;
645 }
646 return Math.sqrt(Math.max(lambda, 1e-300)) * 1.02;
647}