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 { angleDiff, hingeAngle } from "./hinge";
13import type { BeadSystem, Method, Modes, SimResult } from "./types";
15const YOSHIDA_W1 = 1 / (2 - Math.cbrt(2));
16const YOSHIDA_W0 = -Math.cbrt(2) * YOSHIDA_W1;
18export function stepCoefficients(method: Method): number[] {
19 return method === "yoshida4" ? [YOSHIDA_W1, YOSHIDA_W0, YOSHIDA_W1] : [1];
20}
22// Largest h*omega for which the method is stable on x'' = -omega^2 x.
23export function stabilityLimit(method: Method): number {
24 if (method === "verlet") return 2;
25 const cs = stepCoefficients(method);
26 const stable = (h: number) => {
27 // transfer matrix of one full step for omega = 1
28 let m = [1, 0, 0, 1]; // [x<-x, x<-v, v<-x, v<-v]
29 const mul = (a: number[], b: number[]) => [
30 a[0] * b[0] + a[1] * b[2],
31 a[0] * b[1] + a[1] * b[3],
32 a[2] * b[0] + a[3] * b[2],
33 a[2] * b[1] + a[3] * b[3],
34 ];
35 for (const c of cs) {
36 const kick = [1, 0, (-c * h) / 2, 1];
37 const drift = [1, c * h, 0, 1];
38 m = mul(kick, mul(drift, mul(kick, m)));
39 }
40 return Math.abs(m[0] + m[3]) < 2;
41 };
42 let h = 1e-3;
43 while (stable(h + 1e-3)) h += 1e-3;
44 return h;
45}
47// Error of the numerical oscillation frequency for x'' = -omega^2 x, in cents
48// (positive: the simulated tone is sharp).
49export function frequencyErrorCents(method: Method, omega: number, h: number): number {
50 let m = [1, 0, 0, 1];
51 const mul = (a: number[], b: number[]) => [
52 a[0] * b[0] + a[1] * b[2],
53 a[0] * b[1] + a[1] * b[3],
54 a[2] * b[0] + a[3] * b[2],
55 a[2] * b[1] + a[3] * b[3],
56 ];
57 for (const c of stepCoefficients(method)) {
58 const kick = [1, 0, (-c * h * omega * omega) / 2, 1];
59 const drift = [1, c * h, 0, 1];
60 m = mul(kick, mul(drift, mul(kick, m)));
61 }
62 const half = (m[0] + m[3]) / 2;
63 if (Math.abs(half) >= 1) return NaN;
64 const omegaNum = Math.acos(half) / h;
65 return 1200 * Math.log2(omegaNum / omega);
66}
68export interface SimOptions {
69 modes?: Modes | null;
70 onProgress?: (t: number) => void;
71 // called at every frame record with the full-precision state (do not keep the arrays)
72 onFrame?: (frame: number, t: number, x: Float64Array, v: Float64Array) => void;
73 maxAnimBytes?: number;
74}
76export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
77 const t0 = performance.now();
78 const warnings: string[] = [];
79 const { n, mass, pinned } = sys;
80 const N3 = 3 * n;
81 const x = new Float64Array(sys.pos0);
82 const v = new Float64Array(N3);
83 const F = new Float64Array(N3);
84 const invm = new Float64Array(n);
85 for (let i = 0; i < n; i++) invm[i] = pinned[i] ? 0 : 1 / mass[i];
87 // ---- forces ------------------------------------------------------------
88 const { si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb, hi, hj, hk, hl, hkh, hch, hth0 } = sys;
89 const nhg = hi.length;
90 const hg = new Float64Array(12 * nhg); // hinge gradients, refreshed where needed
91 const hb = (h: number) => [hi[h], hj[h], hk[h], hl[h]]; // the hinge's beads, in gradient order
92 const ns = si.length;
93 const nb = bi.length;
94 const internalForces = (): number => {
95 F.fill(0);
96 let pe = 0;
97 for (let s = 0; s < ns; s++) {
98 const i3 = 3 * si[s];
99 const j3 = 3 * sj[s];
100 const dx = x[j3] - x[i3];
101 const dy = x[j3 + 1] - x[i3 + 1];
102 const dz = x[j3 + 2] - x[i3 + 2];
103 const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
104 const e = L - sL0[s];
105 const k = sk[s];
106 const k3 = sk3[s];
107 const e2 = e * e;
108 pe += 0.5 * k * e2 + 0.25 * k3 * e2 * e2;
109 const g = (k * e + k3 * e2 * e) / L;
110 F[i3] += g * dx;
111 F[i3 + 1] += g * dy;
112 F[i3 + 2] += g * dz;
113 F[j3] -= g * dx;
114 F[j3 + 1] -= g * dy;
115 F[j3 + 2] -= g * dz;
116 }
117 for (let b = 0; b < nb; b++) {
118 const i3 = 3 * bi[b];
119 const j3 = 3 * bj[b];
120 const k3 = 3 * bk[b];
121 const kb = bkb[b];
122 const d0 = x[i3] - 2 * x[j3] + x[k3];
123 const d1 = x[i3 + 1] - 2 * x[j3 + 1] + x[k3 + 1];
124 const d2 = x[i3 + 2] - 2 * x[j3 + 2] + x[k3 + 2];
125 pe += 0.5 * kb * (d0 * d0 + d1 * d1 + d2 * d2);
126 F[i3] -= kb * d0;
127 F[i3 + 1] -= kb * d1;
128 F[i3 + 2] -= kb * d2;
129 F[j3] += 2 * kb * d0;
130 F[j3 + 1] += 2 * kb * d1;
131 F[j3 + 2] += 2 * kb * d2;
132 F[k3] -= kb * d0;
133 F[k3 + 1] -= kb * d1;
134 F[k3 + 2] -= kb * d2;
135 }
136 for (let h = 0; h < nhg; h++) {
137 const o = 12 * h;
138 const d = angleDiff(hingeAngle(x, hi[h], hj[h], hk[h], hl[h], hg, o), hth0[h]);
139 pe += 0.5 * hkh[h] * d * d;
140 const s = hkh[h] * d;
141 const i3 = 3 * hi[h], j3 = 3 * hj[h], k3 = 3 * hk[h], l3 = 3 * hl[h];
142 F[i3] -= s * hg[o];
143 F[i3 + 1] -= s * hg[o + 1];
144 F[i3 + 2] -= s * hg[o + 2];
145 F[j3] -= s * hg[o + 3];
146 F[j3 + 1] -= s * hg[o + 4];
147 F[j3 + 2] -= s * hg[o + 5];
148 F[k3] -= s * hg[o + 6];
149 F[k3 + 1] -= s * hg[o + 7];
150 F[k3 + 2] -= s * hg[o + 8];
151 F[l3] -= s * hg[o + 9];
152 F[l3 + 1] -= s * hg[o + 10];
153 F[l3 + 2] -= s * hg[o + 11];
154 }
155 return pe;
156 };
157 const pulses = sys.hits.filter(h => h.duration > 0);
158 const externalForces = (t: number) => {
159 for (const p of pulses) {
160 const u = (t - p.t) / p.duration;
161 if (u <= 0 || u >= 1) continue;
162 const w = 0.5 * (1 - Math.cos(2 * Math.PI * u));
163 for (const b of p.beads) {
164 F[3 * b] += w * p.force[0];
165 F[3 * b + 1] += w * p.force[1];
166 F[3 * b + 2] += w * p.force[2];
167 }
168 }
169 };
171 // ---- rest state checks ---------------------------------------------------
172 let pe = internalForces();
173 const peRest = pe; // energy stored by pre-tension; energies are reported relative to rest
174 let residualForce = 0;
175 for (let i = 0; i < n; i++) {
176 if (pinned[i]) continue;
177 const f = Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]);
178 if (f > residualForce) residualForce = f;
179 }
180 let maxTension = 0;
181 for (let s = 0; s < ns; s++) {
182 const i3 = 3 * si[s];
183 const j3 = 3 * sj[s];
184 const L = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
185 const e = L - sL0[s];
186 maxTension = Math.max(maxTension, Math.abs(sk[s] * e + sk3[s] * e * e * e));
187 }
188 // force scale: the larger of the rest tension and the force of a 1e-3 strain
189 // on any spring or bending triplet
190 let forceScale = maxTension;
191 for (let s = 0; s < ns; s++) forceScale = Math.max(forceScale, 1e-3 * sk[s] * sL0[s]);
192 for (let b = 0; b < nb; b++) {
193 const i3 = 3 * bi[b];
194 const j3 = 3 * bj[b];
195 const len = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
196 forceScale = Math.max(forceScale, 1e-3 * bkb[b] * len);
197 }
198 {
199 const g = new Float64Array(12);
200 for (let h = 0; h < nhg; h++) {
201 hingeAngle(x, hi[h], hj[h], hk[h], hl[h], g);
202 forceScale = Math.max(forceScale, 1e-3 * hkh[h] * Math.max(...g.map(Math.abs)));
203 }
204 }
205 if (residualForce > 1e-6 * Math.max(forceScale, 1e-30)) {
206 warnings.push(
207 `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
208 `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
209 );
210 }
212 // ---- time step -----------------------------------------------------------
213 const omegaMax = estimateOmegaMax(sys);
214 const fs = sys.sampleRate;
215 const hLimit = stabilityLimit(sys.method) / omegaMax;
216 let substeps: number;
217 if (sys.dt > 0) {
218 substeps = Math.max(1, Math.round(1 / (fs * sys.dt)));
219 if (1 / (fs * substeps) > hLimit) {
220 warnings.push(
221 `dt = ${(1 / (fs * substeps)).toExponential(3)} s exceeds the stability limit ${hLimit.toExponential(3)} s.`
222 );
223 }
224 } else {
225 substeps = Math.max(1, Math.ceil(1 / (fs * sys.safety * hLimit)));
226 }
227 const h = 1 / (fs * substeps);
228 const nSamples = Math.max(1, Math.round(sys.duration * fs));
229 const nSteps = nSamples * substeps;
230 const cs = stepCoefficients(sys.method);
232 // ---- damping: v <- exp(-tau A) v ------------------------------------------
233 // dashpots with c = ca = 0 are dropped
234 const keepD = Array.from(sys.dc.keys()).filter(d => sys.dc[d] !== 0 || sys.dca[d] !== 0);
235 const di = Int32Array.from(keepD, d => sys.di[d]);
236 const dj = Int32Array.from(keepD, d => sys.dj[d]);
237 const dc = Float64Array.from(keepD, d => sys.dc[d]);
238 const dca = Float64Array.from(keepD, d => sys.dca[d]);
239 const axialDamped = dca.some(c => c > 0);
240 const dd = new Float64Array(3 * dca.length); // dashpot directions, refreshed where needed
241 const { friction } = sys;
242 const nd = di.length;
243 let normA = 0; // infinity-norm bound of A
244 {
245 const row = new Float64Array(n);
246 for (let i = 0; i < n; i++) row[i] = friction[i];
247 for (let d = 0; d < nd; d++) {
248 row[di[d]] += 2 * dc[d];
249 row[dj[d]] += 2 * dc[d];
250 row[di[d]] += 2 * dca[d];
251 row[dj[d]] += 2 * dca[d];
252 }
253 for (let b = 0; b < nb; b++) {
254 row[bi[b]] += 4 * bcb[b];
255 row[bj[b]] += 8 * bcb[b];
256 row[bk[b]] += 4 * bcb[b];
257 }
258 for (let i = 0; i < n; i++) normA = Math.max(normA, row[i] * invm[i]);
259 }
260 const hingeDamped = hch.some(c => c > 0);
261 const hasDamping = normA > 0 || hingeDamped;
262 const frictionOnly = dc.every(c => c === 0) && !axialDamped && bcb.every(c => c === 0) && !hingeDamped; // then exp(-tau A) is diagonal
263 // Hinge damping depends on the configuration through grad theta, so its
264 // gradients (and its part of the norm bound) are refreshed at the current
265 // positions before each damping step; positions do not change during it.
266 let normAh = 0;
267 const hingeRow = new Float64Array(n);
268 const refreshHinges = () => {
269 if (axialDamped) {
270 for (let d = 0; d < dca.length; d++) {
271 const i3 = 3 * di[d];
272 const j3 = 3 * dj[d];
273 const e0 = x[j3] - x[i3];
274 const e1 = x[j3 + 1] - x[i3 + 1];
275 const e2 = x[j3 + 2] - x[i3 + 2];
276 const L = Math.sqrt(e0 * e0 + e1 * e1 + e2 * e2) || 1;
277 dd[3 * d] = e0 / L;
278 dd[3 * d + 1] = e1 / L;
279 dd[3 * d + 2] = e2 / L;
280 }
281 }
282 if (!hingeDamped) return;
283 hingeRow.fill(0);
284 const row = hingeRow;
285 for (let h = 0; h < nhg; h++) {
286 const o = 12 * h;
287 hingeAngle(x, hi[h], hj[h], hk[h], hl[h], hg, o);
288 let total = 0;
289 for (let q = 0; q < 12; q++) total += Math.abs(hg[o + q]);
290 const c = hch[h] * total;
291 row[hi[h]] += c * (Math.abs(hg[o]) + Math.abs(hg[o + 1]) + Math.abs(hg[o + 2]));
292 row[hj[h]] += c * (Math.abs(hg[o + 3]) + Math.abs(hg[o + 4]) + Math.abs(hg[o + 5]));
293 row[hk[h]] += c * (Math.abs(hg[o + 6]) + Math.abs(hg[o + 7]) + Math.abs(hg[o + 8]));
294 row[hl[h]] += c * (Math.abs(hg[o + 9]) + Math.abs(hg[o + 10]) + Math.abs(hg[o + 11]));
295 }
296 normAh = 0;
297 for (let i = 0; i < n; i++) normAh = Math.max(normAh, row[i] * invm[i]);
298 };
299 const term = new Float64Array(N3);
300 const tmp = new Float64Array(N3);
301 const fa = new Float64Array(n);
302 for (let i = 0; i < n; i++) fa[i] = friction[i] * invm[i];
303 const applyA = (w: Float64Array, out: Float64Array) => {
304 for (let i = 0; i < n; i++) {
305 const a = fa[i];
306 out[3 * i] = a * w[3 * i];
307 out[3 * i + 1] = a * w[3 * i + 1];
308 out[3 * i + 2] = a * w[3 * i + 2];
309 }
310 for (let d = 0; d < nd; d++) {
311 const i3 = 3 * di[d];
312 const j3 = 3 * dj[d];
313 const r0 = w[i3] - w[j3];
314 const r1 = w[i3 + 1] - w[j3 + 1];
315 const r2 = w[i3 + 2] - w[j3 + 2];
316 // force on i: -c r - ca (d . r) d, with r = w_i - w_j
317 let f0 = dc[d] * r0;
318 let f1 = dc[d] * r1;
319 let f2 = dc[d] * r2;
320 if (dca[d] !== 0) {
321 const u0 = dd[3 * d], u1 = dd[3 * d + 1], u2 = dd[3 * d + 2];
322 const s = dca[d] * (u0 * r0 + u1 * r1 + u2 * r2);
323 f0 += s * u0;
324 f1 += s * u1;
325 f2 += s * u2;
326 }
327 const mi = invm[di[d]];
328 const mj = invm[dj[d]];
329 out[i3] += mi * f0;
330 out[i3 + 1] += mi * f1;
331 out[i3 + 2] += mi * f2;
332 out[j3] -= mj * f0;
333 out[j3 + 1] -= mj * f1;
334 out[j3 + 2] -= mj * f2;
335 }
336 for (let b = 0; b < nb; b++) {
337 const cb = bcb[b];
338 if (cb === 0) continue;
339 const i3 = 3 * bi[b];
340 const j3 = 3 * bj[b];
341 const k3 = 3 * bk[b];
342 const ci = cb * invm[bi[b]];
343 const cj = -2 * cb * invm[bj[b]];
344 const ck = cb * invm[bk[b]];
345 const r0 = w[i3] - 2 * w[j3] + w[k3];
346 const r1 = w[i3 + 1] - 2 * w[j3 + 1] + w[k3 + 1];
347 const r2 = w[i3 + 2] - 2 * w[j3 + 2] + w[k3 + 2];
348 out[i3] += ci * r0;
349 out[i3 + 1] += ci * r1;
350 out[i3 + 2] += ci * r2;
351 out[j3] += cj * r0;
352 out[j3 + 1] += cj * r1;
353 out[j3 + 2] += cj * r2;
354 out[k3] += ck * r0;
355 out[k3 + 1] += ck * r1;
356 out[k3 + 2] += ck * r2;
357 }
358 if (hingeDamped) {
359 for (let h = 0; h < nhg; h++) {
360 const ch = hch[h];
361 if (ch === 0) continue;
362 const o = 12 * h;
363 const i3 = 3 * hi[h], j3 = 3 * hj[h], k3 = 3 * hk[h], l3 = 3 * hl[h];
364 let s =
365 hg[o] * w[i3] + hg[o + 1] * w[i3 + 1] + hg[o + 2] * w[i3 + 2] +
366 hg[o + 3] * w[j3] + hg[o + 4] * w[j3 + 1] + hg[o + 5] * w[j3 + 2] +
367 hg[o + 6] * w[k3] + hg[o + 7] * w[k3 + 1] + hg[o + 8] * w[k3 + 2] +
368 hg[o + 9] * w[l3] + hg[o + 10] * w[l3 + 1] + hg[o + 11] * w[l3 + 2];
369 s *= ch;
370 const ai = s * invm[hi[h]], aj = s * invm[hj[h]], ak = s * invm[hk[h]], al = s * invm[hl[h]];
371 out[i3] += ai * hg[o];
372 out[i3 + 1] += ai * hg[o + 1];
373 out[i3 + 2] += ai * hg[o + 2];
374 out[j3] += aj * hg[o + 3];
375 out[j3 + 1] += aj * hg[o + 4];
376 out[j3 + 2] += aj * hg[o + 5];
377 out[k3] += ak * hg[o + 6];
378 out[k3 + 1] += ak * hg[o + 7];
379 out[k3 + 2] += ak * hg[o + 8];
380 out[l3] += al * hg[o + 9];
381 out[l3 + 1] += al * hg[o + 10];
382 out[l3 + 2] += al * hg[o + 11];
383 }
384 }
385 };
386 const kinetic = () => {
387 let ke = 0;
388 for (let i = 0; i < n; i++) {
389 ke += 0.5 * mass[i] * (v[3 * i] ** 2 + v[3 * i + 1] ** 2 + v[3 * i + 2] ** 2);
390 }
391 return ke;
392 };
393 let dissipated = 0;
394 const applyDamping = (tau: number) => {
395 if (!hasDamping || tau === 0) return;
396 const keBefore = kinetic();
397 if (frictionOnly) {
398 for (let i = 0; i < n; i++) {
399 const f = Math.exp(-tau * friction[i] * invm[i]);
400 v[3 * i] *= f;
401 v[3 * i + 1] *= f;
402 v[3 * i + 2] *= f;
403 }
404 } else {
405 // Taylor series with scaling; |tau_s| * normA <= 0.5 so terms shrink fast.
406 refreshHinges();
407 const nsub = Math.max(1, Math.ceil((Math.abs(tau) * (normA + normAh)) / 0.5));
408 const ts = tau / nsub;
409 for (let r = 0; r < nsub; r++) {
410 let vmax = 0;
411 for (let q = 0; q < N3; q++) {
412 term[q] = v[q];
413 vmax = Math.max(vmax, Math.abs(v[q]));
414 }
415 if (vmax === 0) break;
416 for (let k = 1; k < 40; k++) {
417 applyA(term, tmp);
418 const c = -ts / k;
419 let tmax = 0;
420 for (let q = 0; q < N3; q++) {
421 const val = c * tmp[q];
422 term[q] = val;
423 v[q] += val;
424 const a = Math.abs(val);
425 if (a > tmax) tmax = a;
426 }
427 if (tmax <= 1e-17 * vmax) break;
428 }
429 }
430 }
431 dissipated += keBefore - kinetic();
432 };
434 // ---- recording -------------------------------------------------------------
435 const audio = new Float32Array(nSamples);
436 const lb = sys.listenBeads;
437 const [ldx, ldy, ldz] = sys.listenDir;
438 const pos0 = sys.pos0;
439 // Acceleration is (F - damping force) / m on each listener bead, with F the
440 // spring, bending and pulse forces at the current positions and time (the
441 // force array is current at every sample boundary). The damping force uses
442 // the dashpots and bending triplets that touch each listener bead.
443 const touchD = Array.from(lb, b => Array.from(di.keys()).filter(d => di[d] === b || dj[d] === b));
444 const touchB = Array.from(lb, b => Array.from(bi.keys()).filter(t => (bcb[t] !== 0 && (bi[t] === b || bj[t] === b || bk[t] === b))));
445 const touchH = Array.from(lb, b => Array.from(hi.keys()).filter(h => hch[h] !== 0 && hb(h).includes(b)));
446 const dampForce = (v: Float64Array, q: number, b: number, c: number) => {
447 let f = -friction[b] * v[3 * b + c];
448 for (const d of touchD[q]) {
449 const other = di[d] === b ? dj[d] : di[d];
450 f -= dc[d] * (v[3 * b + c] - v[3 * other + c]);
451 if (dca[d] !== 0) {
452 const rel = [0, 1, 2].reduce((acc, q) => acc + dd[3 * d + q] * (v[3 * b + q] - v[3 * other + q]), 0);
453 f -= dca[d] * rel * dd[3 * d + c];
454 }
455 }
456 for (const t of touchB[q]) {
457 const r = v[3 * bi[t] + c] - 2 * v[3 * bj[t] + c] + v[3 * bk[t] + c];
458 const w = bj[t] === b ? -2 : 1;
459 f -= bcb[t] * w * r;
460 }
461 for (const h of touchH[q]) {
462 const o = 12 * h;
463 const beads = hb(h);
464 let s = 0;
465 for (let p = 0; p < 4; p++) {
466 const b3 = 3 * beads[p];
467 s += hg[o + 3 * p] * v[b3] + hg[o + 3 * p + 1] * v[b3 + 1] + hg[o + 3 * p + 2] * v[b3 + 2];
468 }
469 f -= hch[h] * s * hg[o + 3 * beads.indexOf(b) + c];
470 }
471 return f;
472 };
473 // At a sample boundary the damping for the time pendingD has not been
474 // applied yet; vb = v - pendingD A v brings the velocities to the boundary
475 // (to second order) for the listener's velocity and damping force.
476 const vb = new Float64Array(N3);
477 const Av = new Float64Array(N3);
478 const listen = (pendingD: number) => {
479 let s = 0;
480 if (sys.listenQuantity !== "displacement") {
481 refreshHinges(); // hinge damping gradients at the current positions
482 if (pendingD !== 0 && hasDamping) {
483 applyA(v, Av);
484 for (let q = 0; q < N3; q++) vb[q] = v[q] - pendingD * Av[q];
485 } else vb.set(v);
486 }
487 lb.forEach((b, q) => {
488 const dir = [ldx, ldy, ldz];
489 for (let c = 0; c < 3; c++) {
490 if (dir[c] === 0) continue;
491 let val: number;
492 if (sys.listenQuantity === "displacement") val = x[3 * b + c] - pos0[3 * b + c];
493 else {
494 if (sys.listenQuantity === "velocity") val = vb[3 * b + c];
495 else val = (F[3 * b + c] + dampForce(vb, q, b, c)) * invm[b];
496 }
497 s += val * dir[c];
498 }
499 });
500 return s / lb.length;
501 };
502 const frameStride = Math.max(1, Math.round(fs / sys.frameRate));
503 const nFrames = Math.floor((nSamples - 1) / frameStride) + 1;
504 const maxAnimBytes = opts.maxAnimBytes ?? 150e6;
505 const animStride = frameStride * Math.max(1, Math.ceil((nFrames * N3 * 4) / maxAnimBytes));
506 const frameTimes = new Float64Array(nFrames);
507 const energy = new Float64Array(nFrames);
508 const dissipatedRec = new Float64Array(nFrames);
509 const animTimes: number[] = [];
510 const animFrames: Float32Array[] = [];
512 const modes = opts.modes ?? null;
513 const nModes = modes ? modes.freqs.length : 0;
514 const modeEnergy = modes ? new Float32Array(nFrames * nModes) : null;
515 let bz: Float64Array | null = null;
516 let bv: Float64Array | null = null;
517 let om2: Float64Array | null = null;
518 if (modes) {
519 bz = new Float64Array(modes.dof.length);
520 bv = new Float64Array(modes.dof.length);
521 om2 = new Float64Array(nModes);
522 for (let k = 0; k < nModes; k++) om2[k] = (2 * Math.PI * modes.freqs[k]) ** 2;
523 }
524 const recordModes = (f: number) => {
525 if (!modes || !modeEnergy || !bz || !bv || !om2) return;
526 const dof = modes.dof;
527 const nf = dof.length;
528 for (let r = 0; r < nf; r++) {
529 const q = dof[r];
530 const sm = Math.sqrt(mass[(q / 3) | 0]);
531 bz[r] = sm * (x[q] - pos0[q]);
532 bv[r] = sm * v[q];
533 }
534 const psi = modes.psi;
535 for (let k = 0; k < nModes; k++) {
536 let q = 0;
537 let qd = 0;
538 const off = k * nf;
539 for (let r = 0; r < nf; r++) {
540 q += psi[off + r] * bz[r];
541 qd += psi[off + r] * bv[r];
542 }
543 modeEnergy[f * nModes + k] = 0.5 * (qd * qd + om2[k] * q * q);
544 }
545 };
547 // ---- impulses, by audio sample ----------------------------------------------
548 const impulsesAt = new Map<number, typeof sys.hits>();
549 for (const hit of sys.hits) {
550 if (hit.duration > 0) continue;
551 const s = Math.round(hit.t * fs);
552 if (s >= nSamples) continue;
553 if (!impulsesAt.has(s)) impulsesAt.set(s, []);
554 impulsesAt.get(s)!.push(hit);
555 }
556 const tExcitationEnd = Math.max(0, ...sys.hits.map(hh => hh.t + hh.duration));
558 // ---- main loop ---------------------------------------------------------------
559 // Each audio sample interval Delta = substeps * h is one macro step:
560 // D(Delta/2) [substeps conservative steps] D(Delta/2). The damping halves of
561 // neighboring macro steps are merged unless something needs the velocities
562 // at the boundary (a frame record or an impulse).
563 const Delta = 1 / fs;
564 pe = internalForces();
565 externalForces(0);
566 let pendingD = 0;
567 let aborted: string | null = null;
568 let frame = 0;
569 let lastProgress = performance.now();
570 for (let s = 0; s <= nSteps; s++) {
571 const t = s * h;
572 const sample = s % substeps === 0 ? s / substeps : -1;
573 const isFrame = sample >= 0 && sample % frameStride === 0 && frame < nFrames;
574 const imps = sample >= 0 ? impulsesAt.get(sample) : undefined;
575 if ((isFrame || imps) && pendingD !== 0) {
576 applyDamping(pendingD);
577 pendingD = 0;
578 }
579 if (imps) {
580 for (const hit of imps) {
581 for (const b of hit.beads) {
582 if (pinned[b]) continue;
583 v[3 * b] += hit.velocity[0];
584 v[3 * b + 1] += hit.velocity[1];
585 v[3 * b + 2] += hit.velocity[2];
586 }
587 }
588 }
589 if (sample >= 0 && sample < nSamples) audio[sample] = listen(pendingD);
590 if (isFrame) {
591 frameTimes[frame] = t;
592 energy[frame] = kinetic() + pe - peRest;
593 dissipatedRec[frame] = dissipated;
594 recordModes(frame);
595 opts.onFrame?.(frame, t, x, v);
596 if (sample % animStride === 0) {
597 const d = new Float32Array(N3);
598 for (let q = 0; q < N3; q++) d[q] = x[q] - pos0[q];
599 animTimes.push(t);
600 animFrames.push(d);
601 }
602 frame++;
603 }
604 if (s === nSteps) break;
605 if (sample >= 0) {
606 pendingD += Delta / 2;
607 applyDamping(pendingD);
608 pendingD = 0;
609 }
611 let tj = t;
612 for (const c of cs) {
613 const kh = (c * h) / 2;
614 for (let i = 0; i < n; i++) {
615 const im = invm[i] * kh;
616 v[3 * i] += im * F[3 * i];
617 v[3 * i + 1] += im * F[3 * i + 1];
618 v[3 * i + 2] += im * F[3 * i + 2];
619 }
620 const dh = c * h;
621 for (let q = 0; q < N3; q++) x[q] += dh * v[q];
622 tj += c * h;
623 pe = internalForces();
624 externalForces(tj);
625 for (let i = 0; i < n; i++) {
626 const im = invm[i] * kh;
627 v[3 * i] += im * F[3 * i];
628 v[3 * i + 1] += im * F[3 * i + 1];
629 v[3 * i + 2] += im * F[3 * i + 2];
630 }
631 }
632 if ((s + 1) % substeps === 0) pendingD = Delta / 2;
634 if ((s & 1023) === 0) {
635 if (!Number.isFinite(pe)) {
636 aborted = `The solution blew up at t = ${t.toFixed(4)} s. Reduce dt (or safety) or the strike strength.`;
637 break;
638 }
639 const now = performance.now();
640 if (opts.onProgress && now - lastProgress > 200) {
641 opts.onProgress(t);
642 lastProgress = now;
643 }
644 }
645 }
647 // energy drift after the last excitation
648 let energyDrift = 0;
649 let ref = -1;
650 for (let f = 0; f < frame; f++) {
651 if (frameTimes[f] < tExcitationEnd) continue;
652 const total = energy[f] + dissipatedRec[f];
653 if (ref < 0) ref = total;
654 else if (ref > 0) energyDrift = Math.max(energyDrift, Math.abs(total - ref) / ref);
655 }
657 return {
658 audio,
659 sampleRate: fs,
660 dt: h,
661 substeps,
662 omegaMax,
663 frameTimes: frameTimes.slice(0, frame),
664 energy: energy.slice(0, frame),
665 dissipated: dissipatedRec.slice(0, frame),
666 animTimes: Float64Array.from(animTimes),
667 animFrames,
668 modeEnergy: modeEnergy ? modeEnergy.slice(0, frame * nModes) : null,
669 residualForce,
670 maxTension,
671 energyDrift,
672 freqErrors: [1000, 2000, 5000, 10000, 20000]
673 .filter(hz => hz < fs / 2)
674 .map(hz => ({ hz, cents: frequencyErrorCents(sys.method, 2 * Math.PI * hz, h) })),
675 elapsedMs: performance.now() - t0,
676 aborted,
677 warnings,
678 };
679}
681// Highest frequency (rad/s) of the system linearized about the rest positions,
682// by power iteration on M^{-1/2} H M^{-1/2}, with a small margin.
683export function estimateOmegaMax(sys: BeadSystem): number {
684 const { n, si, sj, sk, sL0, sk3, pos0, mass, pinned } = sys;
685 const ns = si.length;
686 // per-spring 3x3 blocks B = keff d d^T + (f/L)(I - d d^T)
687 const B = new Float64Array(ns * 6);
688 for (let s = 0; s < ns; s++) {
689 const i3 = 3 * si[s];
690 const j3 = 3 * sj[s];
691 const dx = pos0[j3] - pos0[i3];
692 const dy = pos0[j3 + 1] - pos0[i3 + 1];
693 const dz = pos0[j3 + 2] - pos0[i3 + 2];
694 const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
695 const e = L - sL0[s];
696 const keff = sk[s] + 3 * sk3[s] * e * e;
697 const g = (sk[s] * e + sk3[s] * e * e * e) / L;
698 const u = [dx / L, dy / L, dz / L];
699 const a = keff - g;
700 B[6 * s] = a * u[0] * u[0] + g;
701 B[6 * s + 1] = a * u[1] * u[1] + g;
702 B[6 * s + 2] = a * u[2] * u[2] + g;
703 B[6 * s + 3] = a * u[0] * u[1];
704 B[6 * s + 4] = a * u[0] * u[2];
705 B[6 * s + 5] = a * u[1] * u[2];
706 }
707 const ism = new Float64Array(n);
708 for (let i = 0; i < n; i++) ism[i] = pinned[i] ? 0 : 1 / Math.sqrt(mass[i]);
709 // hinges at rest: stiffness kh grad(theta) grad(theta)^T (theta = theta0 there)
710 const hg0 = new Float64Array(12 * sys.hi.length);
711 for (let h = 0; h < sys.hi.length; h++) hingeAngle(pos0, sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h], hg0.subarray(12 * h, 12 * h + 12));
712 let w = new Float64Array(3 * n);
713 let seed = 12345;
714 for (let q = 0; q < 3 * n; q++) {
715 seed = (seed * 1103515245 + 12345) & 0x7fffffff;
716 w[q] = ism[Math.floor(q / 3)] * (seed / 0x7fffffff - 0.5);
717 }
718 const u = new Float64Array(3 * n);
719 const y = new Float64Array(3 * n);
720 let lambda = 0;
721 for (let it = 0; it < 500; it++) {
722 let norm = 0;
723 for (let q = 0; q < 3 * n; q++) norm += w[q] * w[q];
724 norm = Math.sqrt(norm);
725 if (norm === 0) return 1;
726 for (let q = 0; q < 3 * n; q++) u[q] = (w[q] / norm) * ism[Math.floor(q / 3)];
727 y.fill(0);
728 for (let s = 0; s < ns; s++) {
729 const i3 = 3 * si[s];
730 const j3 = 3 * sj[s];
731 const r0 = u[i3] - u[j3];
732 const r1 = u[i3 + 1] - u[j3 + 1];
733 const r2 = u[i3 + 2] - u[j3 + 2];
734 const b = 6 * s;
735 const f0 = B[b] * r0 + B[b + 3] * r1 + B[b + 4] * r2;
736 const f1 = B[b + 3] * r0 + B[b + 1] * r1 + B[b + 5] * r2;
737 const f2 = B[b + 4] * r0 + B[b + 5] * r1 + B[b + 2] * r2;
738 y[i3] += f0;
739 y[i3 + 1] += f1;
740 y[i3 + 2] += f2;
741 y[j3] -= f0;
742 y[j3 + 1] -= f1;
743 y[j3 + 2] -= f2;
744 }
745 for (let h = 0; h < sys.hi.length; h++) {
746 const o = 12 * h;
747 const beads = [sys.hi[h], sys.hj[h], sys.hk[h], sys.hl[h]];
748 let s = 0;
749 for (let p = 0; p < 4; p++) for (let c = 0; c < 3; c++) s += hg0[o + 3 * p + c] * u[3 * beads[p] + c];
750 s *= sys.hkh[h];
751 for (let p = 0; p < 4; p++) for (let c = 0; c < 3; c++) y[3 * beads[p] + c] += s * hg0[o + 3 * p + c];
752 }
753 for (let b = 0; b < sys.bi.length; b++) {
754 const i3 = 3 * sys.bi[b];
755 const j3 = 3 * sys.bj[b];
756 const k3 = 3 * sys.bk[b];
757 const kb = sys.bkb[b];
758 for (let q = 0; q < 3; q++) {
759 const r = kb * (u[i3 + q] - 2 * u[j3 + q] + u[k3 + q]);
760 y[i3 + q] += r;
761 y[j3 + q] -= 2 * r;
762 y[k3 + q] += r;
763 }
764 }
765 let rq = 0;
766 const wNext = new Float64Array(3 * n);
767 for (let q = 0; q < 3 * n; q++) {
768 wNext[q] = y[q] * ism[Math.floor(q / 3)];
769 rq += (w[q] / norm) * wNext[q];
770 }
771 w = wNext;
772 if (it > 20 && Math.abs(rq - lambda) < 1e-8 * Math.abs(rq)) {
773 lambda = rq;
774 break;
775 }
776 lambda = rq;
777 }
778 return Math.sqrt(Math.max(lambda, 1e-300)) * 1.02;
779}