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;
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;
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 ------------------------------------------------------------
87 const { si, sj, sk, sL0, sk3 } = sys;
88 const ns = si.length;
89 const internalForces = (): number => {
90 F.fill(0);
91 let pe = 0;
92 for (let s = 0; s < ns; s++) {
93 const i3 = 3 * si[s];
94 const j3 = 3 * sj[s];
95 const dx = x[j3] - x[i3];
96 const dy = x[j3 + 1] - x[i3 + 1];
97 const dz = x[j3 + 2] - x[i3 + 2];
98 const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
99 const e = L - sL0[s];
100 const k = sk[s];
101 const k3 = sk3[s];
102 const e2 = e * e;
103 pe += 0.5 * k * e2 + 0.25 * k3 * e2 * e2;
104 const g = (k * e + k3 * e2 * e) / L;
105 F[i3] += g * dx;
106 F[i3 + 1] += g * dy;
107 F[i3 + 2] += g * dz;
108 F[j3] -= g * dx;
109 F[j3 + 1] -= g * dy;
110 F[j3 + 2] -= g * dz;
111 }
112 return pe;
113 };
114 const pulses = sys.hits.filter(h => h.duration > 0);
115 const externalForces = (t: number) => {
116 for (const p of pulses) {
117 const u = (t - p.t) / p.duration;
118 if (u <= 0 || u >= 1) continue;
119 const w = 0.5 * (1 - Math.cos(2 * Math.PI * u));
120 for (const b of p.beads) {
121 F[3 * b] += w * p.force[0];
122 F[3 * b + 1] += w * p.force[1];
123 F[3 * b + 2] += w * p.force[2];
124 }
125 }
126 };
128 // ---- rest state checks ---------------------------------------------------
129 let pe = internalForces();
130 const peRest = pe; // energy stored by pre-tension; energies are reported relative to rest
131 let residualForce = 0;
132 for (let i = 0; i < n; i++) {
133 if (pinned[i]) continue;
134 const f = Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]);
135 if (f > residualForce) residualForce = f;
136 }
137 let maxTension = 0;
138 for (let s = 0; s < ns; s++) {
139 const i3 = 3 * si[s];
140 const j3 = 3 * sj[s];
141 const L = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
142 const e = L - sL0[s];
143 maxTension = Math.max(maxTension, Math.abs(sk[s] * e + sk3[s] * e * e * e));
144 }
145 if (residualForce > 1e-6 * Math.max(maxTension, 1e-30)) {
146 warnings.push(
147 `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
148 `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
149 );
150 }
152 // ---- time step -----------------------------------------------------------
153 const omegaMax = estimateOmegaMax(sys);
154 const fs = sys.sampleRate;
155 const hLimit = stabilityLimit(sys.method) / omegaMax;
156 let substeps: number;
157 if (sys.dt > 0) {
158 substeps = Math.max(1, Math.round(1 / (fs * sys.dt)));
159 if (1 / (fs * substeps) > hLimit) {
160 warnings.push(
161 `dt = ${(1 / (fs * substeps)).toExponential(3)} s exceeds the stability limit ${hLimit.toExponential(3)} s.`
162 );
163 }
164 } else {
165 substeps = Math.max(1, Math.ceil(1 / (fs * sys.safety * hLimit)));
166 }
167 const h = 1 / (fs * substeps);
168 const nSamples = Math.max(1, Math.round(sys.duration * fs));
169 const nSteps = nSamples * substeps;
170 const cs = stepCoefficients(sys.method);
172 // ---- damping: v <- exp(-tau A) v ------------------------------------------
173 const { di, dj, dc, friction } = sys;
174 const nd = di.length;
175 let normA = 0; // infinity-norm bound of A
176 {
177 const row = new Float64Array(n);
178 for (let i = 0; i < n; i++) row[i] = friction[i];
179 for (let d = 0; d < nd; d++) {
180 row[di[d]] += 2 * dc[d];
181 row[dj[d]] += 2 * dc[d];
182 }
183 for (let i = 0; i < n; i++) normA = Math.max(normA, row[i] * invm[i]);
184 }
185 const hasDamping = normA > 0;
186 const frictionOnly = dc.every(c => c === 0); // then exp(-tau A) is diagonal
187 const term = new Float64Array(N3);
188 const tmp = new Float64Array(N3);
189 const fa = new Float64Array(n);
190 for (let i = 0; i < n; i++) fa[i] = friction[i] * invm[i];
191 const dai = new Float64Array(nd);
192 const daj = new Float64Array(nd);
193 for (let d = 0; d < nd; d++) {
194 dai[d] = dc[d] * invm[di[d]];
195 daj[d] = dc[d] * invm[dj[d]];
196 }
197 const applyA = (w: Float64Array, out: Float64Array) => {
198 for (let i = 0; i < n; i++) {
199 const a = fa[i];
200 out[3 * i] = a * w[3 * i];
201 out[3 * i + 1] = a * w[3 * i + 1];
202 out[3 * i + 2] = a * w[3 * i + 2];
203 }
204 for (let d = 0; d < nd; d++) {
205 const i3 = 3 * di[d];
206 const j3 = 3 * dj[d];
207 const ai = dai[d];
208 const aj = daj[d];
209 const r0 = w[i3] - w[j3];
210 const r1 = w[i3 + 1] - w[j3 + 1];
211 const r2 = w[i3 + 2] - w[j3 + 2];
212 out[i3] += ai * r0;
213 out[i3 + 1] += ai * r1;
214 out[i3 + 2] += ai * r2;
215 out[j3] -= aj * r0;
216 out[j3 + 1] -= aj * r1;
217 out[j3 + 2] -= aj * r2;
218 }
219 };
220 const kinetic = () => {
221 let ke = 0;
222 for (let i = 0; i < n; i++) {
223 ke += 0.5 * mass[i] * (v[3 * i] ** 2 + v[3 * i + 1] ** 2 + v[3 * i + 2] ** 2);
224 }
225 return ke;
226 };
227 let dissipated = 0;
228 const applyDamping = (tau: number) => {
229 if (!hasDamping || tau === 0) return;
230 const keBefore = kinetic();
231 if (frictionOnly) {
232 for (let i = 0; i < n; i++) {
233 const f = Math.exp(-tau * friction[i] * invm[i]);
234 v[3 * i] *= f;
235 v[3 * i + 1] *= f;
236 v[3 * i + 2] *= f;
237 }
238 } else {
239 // Taylor series with scaling; |tau_s| * normA <= 0.5 so terms shrink fast.
240 const nsub = Math.max(1, Math.ceil((Math.abs(tau) * normA) / 0.5));
241 const ts = tau / nsub;
242 for (let r = 0; r < nsub; r++) {
243 let vmax = 0;
244 for (let q = 0; q < N3; q++) {
245 term[q] = v[q];
246 vmax = Math.max(vmax, Math.abs(v[q]));
247 }
248 if (vmax === 0) break;
249 for (let k = 1; k < 40; k++) {
250 applyA(term, tmp);
251 const c = -ts / k;
252 let tmax = 0;
253 for (let q = 0; q < N3; q++) {
254 const val = c * tmp[q];
255 term[q] = val;
256 v[q] += val;
257 const a = Math.abs(val);
258 if (a > tmax) tmax = a;
259 }
260 if (tmax <= 1e-17 * vmax) break;
261 }
262 }
263 }
264 dissipated += keBefore - kinetic();
265 };
267 // ---- recording -------------------------------------------------------------
268 const audio = new Float32Array(nSamples);
269 const lb = sys.listenBeads;
270 const [ldx, ldy, ldz] = sys.listenDir;
271 const pos0 = sys.pos0;
272 const listen = () => {
273 let s = 0;
274 for (const b of lb) {
275 s +=
276 (x[3 * b] - pos0[3 * b]) * ldx +
277 (x[3 * b + 1] - pos0[3 * b + 1]) * ldy +
278 (x[3 * b + 2] - pos0[3 * b + 2]) * ldz;
279 }
280 return s / lb.length;
281 };
282 const frameStride = Math.max(1, Math.round(fs / sys.frameRate));
283 const nFrames = Math.floor((nSamples - 1) / frameStride) + 1;
284 const maxAnimBytes = opts.maxAnimBytes ?? 150e6;
285 const animStride = frameStride * Math.max(1, Math.ceil((nFrames * N3 * 4) / maxAnimBytes));
286 const frameTimes = new Float64Array(nFrames);
287 const energy = new Float64Array(nFrames);
288 const dissipatedRec = new Float64Array(nFrames);
289 const animTimes: number[] = [];
290 const animFrames: Float32Array[] = [];
292 const modes = opts.modes ?? null;
293 const nModes = modes ? modes.freqs.length : 0;
294 const modeEnergy = modes ? new Float32Array(nFrames * nModes) : null;
295 let bz: Float64Array | null = null;
296 let bv: Float64Array | null = null;
297 let om2: Float64Array | null = null;
298 if (modes) {
299 bz = new Float64Array(modes.free.length);
300 bv = new Float64Array(modes.free.length);
301 om2 = new Float64Array(nModes);
302 for (let k = 0; k < nModes; k++) om2[k] = (2 * Math.PI * modes.freqs[k]) ** 2;
303 }
304 const recordModes = (f: number) => {
305 if (!modes || !modeEnergy || !bz || !bv || !om2) return;
306 const free = modes.free;
307 const nf = free.length;
308 for (let r = 0; r < nf; r++) {
309 const b = free[r];
310 const sm = Math.sqrt(mass[b]);
311 bz[r] = sm * (x[3 * b + 2] - pos0[3 * b + 2]);
312 bv[r] = sm * v[3 * b + 2];
313 }
314 const psi = modes.psi;
315 for (let k = 0; k < nModes; k++) {
316 let q = 0;
317 let qd = 0;
318 const off = k * nf;
319 for (let r = 0; r < nf; r++) {
320 q += psi[off + r] * bz[r];
321 qd += psi[off + r] * bv[r];
322 }
323 modeEnergy[f * nModes + k] = 0.5 * (qd * qd + om2[k] * q * q);
324 }
325 };
327 // ---- impulses, by audio sample ----------------------------------------------
328 const impulsesAt = new Map<number, typeof sys.hits>();
329 for (const hit of sys.hits) {
330 if (hit.duration > 0) continue;
331 const s = Math.round(hit.t * fs);
332 if (s >= nSamples) continue;
333 if (!impulsesAt.has(s)) impulsesAt.set(s, []);
334 impulsesAt.get(s)!.push(hit);
335 }
336 const tExcitationEnd = Math.max(0, ...sys.hits.map(hh => hh.t + hh.duration));
338 // ---- main loop ---------------------------------------------------------------
339 // Each audio sample interval Delta = substeps * h is one macro step:
340 // D(Delta/2) [substeps conservative steps] D(Delta/2). The damping halves of
341 // neighboring macro steps are merged unless something needs the velocities
342 // at the boundary (a frame record or an impulse).
343 const Delta = 1 / fs;
344 pe = internalForces();
345 externalForces(0);
346 let pendingD = 0;
347 let aborted: string | null = null;
348 let frame = 0;
349 let lastProgress = performance.now();
350 for (let s = 0; s <= nSteps; s++) {
351 const t = s * h;
352 const sample = s % substeps === 0 ? s / substeps : -1;
353 const isFrame = sample >= 0 && sample % frameStride === 0 && frame < nFrames;
354 const imps = sample >= 0 ? impulsesAt.get(sample) : undefined;
355 if ((isFrame || imps) && pendingD !== 0) {
356 applyDamping(pendingD);
357 pendingD = 0;
358 }
359 if (imps) {
360 for (const hit of imps) {
361 for (const b of hit.beads) {
362 if (pinned[b]) continue;
363 v[3 * b] += hit.velocity[0];
364 v[3 * b + 1] += hit.velocity[1];
365 v[3 * b + 2] += hit.velocity[2];
366 }
367 }
368 }
369 if (sample >= 0 && sample < nSamples) audio[sample] = listen();
370 if (isFrame) {
371 frameTimes[frame] = t;
372 energy[frame] = kinetic() + pe - peRest;
373 dissipatedRec[frame] = dissipated;
374 recordModes(frame);
375 opts.onFrame?.(frame, t, x, v);
376 if (sample % animStride === 0) {
377 const d = new Float32Array(N3);
378 for (let q = 0; q < N3; q++) d[q] = x[q] - pos0[q];
379 animTimes.push(t);
380 animFrames.push(d);
381 }
382 frame++;
383 }
384 if (s === nSteps) break;
385 if (sample >= 0) {
386 pendingD += Delta / 2;
387 applyDamping(pendingD);
388 pendingD = 0;
389 }
391 let tj = t;
392 for (const c of cs) {
393 const kh = (c * h) / 2;
394 for (let i = 0; i < n; i++) {
395 const im = invm[i] * kh;
396 v[3 * i] += im * F[3 * i];
397 v[3 * i + 1] += im * F[3 * i + 1];
398 v[3 * i + 2] += im * F[3 * i + 2];
399 }
400 const dh = c * h;
401 for (let q = 0; q < N3; q++) x[q] += dh * v[q];
402 tj += c * h;
403 pe = internalForces();
404 externalForces(tj);
405 for (let i = 0; i < n; i++) {
406 const im = invm[i] * kh;
407 v[3 * i] += im * F[3 * i];
408 v[3 * i + 1] += im * F[3 * i + 1];
409 v[3 * i + 2] += im * F[3 * i + 2];
410 }
411 }
412 if ((s + 1) % substeps === 0) pendingD = Delta / 2;
414 if ((s & 1023) === 0) {
415 if (!Number.isFinite(pe)) {
416 aborted = `The solution blew up at t = ${t.toFixed(4)} s. Reduce dt (or safety) or the strike strength.`;
417 break;
418 }
419 const now = performance.now();
420 if (opts.onProgress && now - lastProgress > 200) {
421 opts.onProgress(t);
422 lastProgress = now;
423 }
424 }
425 }
427 // energy drift after the last excitation
428 let energyDrift = 0;
429 let ref = -1;
430 for (let f = 0; f < frame; f++) {
431 if (frameTimes[f] < tExcitationEnd) continue;
432 const total = energy[f] + dissipatedRec[f];
433 if (ref < 0) ref = total;
434 else if (ref > 0) energyDrift = Math.max(energyDrift, Math.abs(total - ref) / ref);
435 }
437 return {
438 audio,
439 sampleRate: fs,
440 dt: h,
441 substeps,
442 omegaMax,
443 frameTimes: frameTimes.slice(0, frame),
444 energy: energy.slice(0, frame),
445 dissipated: dissipatedRec.slice(0, frame),
446 animTimes: Float64Array.from(animTimes),
447 animFrames,
448 modeEnergy: modeEnergy ? modeEnergy.slice(0, frame * nModes) : null,
449 residualForce,
450 maxTension,
451 energyDrift,
452 freqErrors: [1000, 2000, 5000, 10000, 20000]
453 .filter(hz => hz < fs / 2)
454 .map(hz => ({ hz, cents: frequencyErrorCents(sys.method, 2 * Math.PI * hz, h) })),
455 elapsedMs: performance.now() - t0,
456 aborted,
457 warnings,
458 };
459}
461// Highest frequency (rad/s) of the system linearized about the rest positions,
462// by power iteration on M^{-1/2} H M^{-1/2}, with a small margin.
463export function estimateOmegaMax(sys: BeadSystem): number {
464 const { n, si, sj, sk, sL0, sk3, pos0, mass, pinned } = sys;
465 const ns = si.length;
466 // per-spring 3x3 blocks B = keff d d^T + (f/L)(I - d d^T)
467 const B = new Float64Array(ns * 6);
468 for (let s = 0; s < ns; s++) {
469 const i3 = 3 * si[s];
470 const j3 = 3 * sj[s];
471 const dx = pos0[j3] - pos0[i3];
472 const dy = pos0[j3 + 1] - pos0[i3 + 1];
473 const dz = pos0[j3 + 2] - pos0[i3 + 2];
474 const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
475 const e = L - sL0[s];
476 const keff = sk[s] + 3 * sk3[s] * e * e;
477 const g = (sk[s] * e + sk3[s] * e * e * e) / L;
478 const u = [dx / L, dy / L, dz / L];
479 const a = keff - g;
480 B[6 * s] = a * u[0] * u[0] + g;
481 B[6 * s + 1] = a * u[1] * u[1] + g;
482 B[6 * s + 2] = a * u[2] * u[2] + g;
483 B[6 * s + 3] = a * u[0] * u[1];
484 B[6 * s + 4] = a * u[0] * u[2];
485 B[6 * s + 5] = a * u[1] * u[2];
486 }
487 const ism = new Float64Array(n);
488 for (let i = 0; i < n; i++) ism[i] = pinned[i] ? 0 : 1 / Math.sqrt(mass[i]);
489 let w = new Float64Array(3 * n);
490 let seed = 12345;
491 for (let q = 0; q < 3 * n; q++) {
492 seed = (seed * 1103515245 + 12345) & 0x7fffffff;
493 w[q] = ism[Math.floor(q / 3)] * (seed / 0x7fffffff - 0.5);
494 }
495 const u = new Float64Array(3 * n);
496 const y = new Float64Array(3 * n);
497 let lambda = 0;
498 for (let it = 0; it < 500; it++) {
499 let norm = 0;
500 for (let q = 0; q < 3 * n; q++) norm += w[q] * w[q];
501 norm = Math.sqrt(norm);
502 if (norm === 0) return 1;
503 for (let q = 0; q < 3 * n; q++) u[q] = (w[q] / norm) * ism[Math.floor(q / 3)];
504 y.fill(0);
505 for (let s = 0; s < ns; s++) {
506 const i3 = 3 * si[s];
507 const j3 = 3 * sj[s];
508 const r0 = u[i3] - u[j3];
509 const r1 = u[i3 + 1] - u[j3 + 1];
510 const r2 = u[i3 + 2] - u[j3 + 2];
511 const b = 6 * s;
512 const f0 = B[b] * r0 + B[b + 3] * r1 + B[b + 4] * r2;
513 const f1 = B[b + 3] * r0 + B[b + 1] * r1 + B[b + 5] * r2;
514 const f2 = B[b + 4] * r0 + B[b + 5] * r1 + B[b + 2] * r2;
515 y[i3] += f0;
516 y[i3 + 1] += f1;
517 y[i3 + 2] += f2;
518 y[j3] -= f0;
519 y[j3 + 1] -= f1;
520 y[j3 + 2] -= f2;
521 }
522 let rq = 0;
523 const wNext = new Float64Array(3 * n);
524 for (let q = 0; q < 3 * n; q++) {
525 wNext[q] = y[q] * ism[Math.floor(q / 3)];
526 rq += (w[q] / norm) * wNext[q];
527 }
528 w = wNext;
529 if (it > 20 && Math.abs(rq - lambda) < 1e-8 * Math.abs(rq)) {
530 lambda = rq;
531 break;
532 }
533 lambda = rq;
534 }
535 return Math.sqrt(Math.max(lambda, 1e-300)) * 1.02;
536}