/ proof-of-concept / bead-percussion
Sign in
proof-of-concept / bead-percussion
main
Branches
main
bead-percussion / tests / solver.test.ts
486 lines · 20.0 KBCodeBlameHistory
0330c41Damping by decay times: bp_decay helper, per-mode decay-time plot, examples convertedJeremy Magland 4import { computeModes, modalDamping } from "../src/sim/modes";
185968eRevert "Dome example: chaos that keeps its pitch..."Jeremy Magland 5import { numblEig, runScript } from "../src/sim/runScript";
95f2216Minimal chaos example and perturbed-twin chaos test with separation plot and growth-rate fitJeremy Magland 6import { growthRate } from "../src/sim/analysis";
fc06d37Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser appJeremy Magland 7import { frequencyErrorCents, simulate, stabilityLimit } from "../src/sim/solver";
95f2216Minimal chaos example and perturbed-twin chaos test with separation plot and growth-rate fitJeremy Magland 8import { runTwin } from "../src/sim/twin";
11const helpers = readdirSync("matlab/helpers").map(name => ({
12 name,
13 source: readFileSync(`matlab/helpers/${name}`, "utf8"),
14}));
15const build = (src: string) => runScript(src, helpers, () => {});
17// A string: N free beads between two pinned ends, tension f0, spacing a.
18// Transverse modes: omega_k = 2 sqrt(f0/(m a)) sin(k pi / (2(N+1))).
19const STRING = (extra: string) => `
20N = 15; a = 0.01; m = 1e-3; f0 = 10; k = 1e4;
21[pos, E] = bp_chain(N + 2, (N + 1) * a);
22sys.pos = pos;
23sys.mass = m;
24sys.pinned = [1, N + 2];
25sys.springs = bp_springs(pos, E, k, f0);
26hit.beads = 5; hit.velocity = [0 0 1e-4];
27sys.hits = hit;
28sys.listen = struct('beads', 4);
29sys.duration = 0.2;
30${extra}
31`;
33function peakFrequency(x: Float32Array, fs: number, fLo: number, fHi: number): number {
34 // direct DFT magnitude on a fine grid, refined by a parabola
35 let best = 0;
36 let bestF = 0;
37 const mag = (f: number) => {
38 let re = 0;
39 let im = 0;
40 for (let i = 0; i < x.length; i++) {
41 const w = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (x.length - 1));
42 re += w * x[i] * Math.cos((2 * Math.PI * f * i) / fs);
43 im -= w * x[i] * Math.sin((2 * Math.PI * f * i) / fs);
44 }
45 return Math.hypot(re, im);
46 };
47 for (let f = fLo; f <= fHi; f += 0.5) {
48 const m = mag(f);
49 if (m > best) {
50 best = m;
51 bestF = f;
52 }
53 }
54 const [a, b, c] = [mag(bestF - 0.5), best, mag(bestF + 0.5)];
55 return bestF + (0.5 * 0.5 * (a - c)) / (a - 2 * b + c);
58describe("string", () => {
59 const { sys } = build(STRING(""));
60 const f1 = (2 * Math.sqrt(10 / (1e-3 * 0.01)) * Math.sin(Math.PI / 32)) / (2 * Math.PI);
62 it("has the analytic transverse modes", () => {
63 const modes = computeModes(sys, numblEig);
64 for (let k = 1; k <= 15; k++) {
65 const fk = (2 * Math.sqrt(10 / (1e-3 * 0.01)) * Math.sin((k * Math.PI) / 32)) / (2 * Math.PI);
66 expect(modes.freqs[k - 1]).toBeCloseTo(fk, 6);
67 }
68 });
70 it("rings at the fundamental for a small strike", () => {
71 const s: BeadSystem = { ...sys, hits: [{ ...sys.hits[0], beads: Int32Array.of(8) }], listenBeads: Int32Array.of(8) };
72 const r = simulate(s);
73 expect(r.aborted).toBeNull();
74 const f = peakFrequency(r.audio, r.sampleRate, f1 * 0.8, f1 * 1.2);
75 expect(Math.abs(f - f1) / f1).toBeLessThan(2e-3);
76 });
78 it("is at equilibrium at rest", () => {
79 const r = simulate({ ...sys, duration: 0.001 });
80 expect(r.residualForce).toBeLessThan(1e-9);
81 expect(r.warnings).toEqual([]);
82 });
83});
85describe("energy", () => {
ac09b0fDefault damping: every mode falls 60 dB in 2.5 s; examples at 2.5 s; twin separation relative to recent displacement; stronger chaos strikeJeremy Magland 86 const { sys } = build(STRING("sys.hits.velocity = [0 0.3 2]; sys.friction = 0;"));
88 simulate({ ...sys, method, dt, duration: 0.05 }).energyDrift;
90 it("is conserved to second order with Verlet", () => {
91 const d1 = drift("verlet", 1 / 48000 / 8);
92 const d2 = drift("verlet", 1 / 48000 / 16);
93 expect(d1).toBeLessThan(1e-3);
94 expect(d1 / d2).toBeGreaterThan(3);
95 expect(d1 / d2).toBeLessThan(5);
96 });
98 it("is conserved to fourth order with Yoshida", () => {
99 const d1 = drift("yoshida4", 1 / 48000 / 8);
100 const d2 = drift("yoshida4", 1 / 48000 / 16);
101 expect(d1 / d2).toBeGreaterThan(12);
102 });
104 it("balances with the dissipated energy when damped", () => {
105 const { sys: damped } = build(
106 STRING("sys.hits.velocity = [0 0.3 2]; sys.friction = 2e-3; sys.dashpots = [E, 1e-3*ones(size(E,1),1)];")
107 );
108 // The damping step's energy loss is counted exactly, so E + dissipated
109 // drifts only by the conservative integrator's O(h^2) error.
110 const r1 = simulate({ ...damped, duration: 0.1, dt: 1 / 48000 / 8 });
111 const r2 = simulate({ ...damped, duration: 0.1, dt: 1 / 48000 / 16 });
112 expect(r1.dissipated[r1.dissipated.length - 1]).toBeGreaterThan(0.3 * r1.energy[0]);
113 expect(r1.energyDrift).toBeLessThan(1e-4);
114 expect(r1.energyDrift / r2.energyDrift).toBeGreaterThan(3);
115 });
116});
118describe("damping", () => {
119 it("matches friction when a dashpot ties a bead to a pinned bead", () => {
120 // one free bead between two pinned beads: friction c0 on it is the same
121 // force as two dashpots of c0/2 to its pinned neighbours
122 const base = STRING("").replace("N = 15", "N = 1");
123 const { sys: a } = build(base + "sys.friction = [0; 0.02; 0]; sys.hits.beads = 2; sys.listen.beads = 2;");
124 const { sys: b } = build(base + "sys.dashpots = [1 2 0.01; 2 3 0.01]; sys.hits.beads = 2; sys.listen.beads = 2;");
125 const ra = simulate({ ...a, duration: 0.05 });
126 const rb = simulate({ ...b, duration: 0.05 });
127 let err = 0;
128 let peak = 0;
129 for (let i = 0; i < ra.audio.length; i++) {
130 err = Math.max(err, Math.abs(ra.audio[i] - rb.audio[i]));
131 peak = Math.max(peak, Math.abs(ra.audio[i]));
132 }
133 expect(err / peak).toBeLessThan(1e-6);
134 });
136 it("decays a mode at the predicted rate", () => {
137 // one free bead, friction c0: amplitude decays as exp(-c0/(2m) t)
138 const base = STRING("").replace("N = 15", "N = 1");
139 const { sys } = build(base + "sys.friction = 0.02; sys.hits.beads = 2; sys.listen.beads = 2;");
140 const r = simulate({ ...sys, duration: 0.2 });
141 const rms = (from: number, to: number) => {
142 let s = 0;
143 for (let i = from; i < to; i++) s += r.audio[i] ** 2;
144 return Math.sqrt(s / (to - from));
145 };
146 const fs = r.sampleRate;
147 const ratio = rms(0.15 * fs, 0.2 * fs) / rms(0.05 * fs, 0.1 * fs);
148 const expected = Math.exp((-0.02 / (2 * 1e-3)) * 0.1);
149 expect(Math.abs(ratio - expected) / expected).toBeLessThan(0.02);
150 });
151});
153describe("integrator", () => {
154 it("reports Verlet's known frequency error", () => {
155 // Verlet: omega_num h = 2 asin(omega h / 2)
156 const w = 2 * Math.PI * 5000;
157 const h = 1e-5;
158 const exact = 1200 * Math.log2((2 * Math.asin((w * h) / 2)) / (w * h));
159 expect(frequencyErrorCents("verlet", w, h)).toBeCloseTo(exact, 9);
160 });
162 it("has the expected stability limits", () => {
163 expect(stabilityLimit("verlet")).toBe(2);
164 expect(stabilityLimit("yoshida4")).toBeGreaterThan(1.5);
165 expect(stabilityLimit("yoshida4")).toBeLessThan(1.6);
166 });
167});
169describe("script errors", () => {
170 it("names the field and bead", () => {
171 expect(() => build(STRING("sys.pinned = [1 99];"))).toThrow(/sys.pinned: bead index 99/);
172 expect(() => build(STRING("sys.springs = [1 2 3];"))).toThrow(/4 or 5 columns/);
173 });
174});
176describe("chaos example", () => {
177 const src = readFileSync("matlab/examples/chaos.m", "utf8");
178 const rate = (strength: number) => {
179 const { sys } = build(src.replace(/^s = 1;/m, `s = ${strength};`));
180 const s1 = { ...sys, duration: 1, frameRate: 1000 };
181 const saved: Float64Array[] = [];
182 const r = simulate(s1, { onFrame: (_f, _t, x) => saved.push(x.slice()) });
183 const { sep } = runTwin(s1, saved, r.dt, 1e-9);
184 return growthRate(r.frameTimes, sep, 1e-9);
185 };
187 it("is chaotic when struck hard", () => {
188 const g = rate(1);
189 expect(g).not.toBeNull();
190 expect(g!).toBeGreaterThan(50); // e-folding time under 20 ms
191 });
193 it("is regular when struck softly", () => {
196});
198describe("bp_decay", () => {
199 // String struck in the shape of mode p; its amplitude should decay at
200 // sigma0 + sigma1 w_p^2, the law through T60 = 0.5 s at 100 Hz and 0.1 s at 1 kHz.
201 const shaped = (p: number) =>
202 build(
203 STRING(`
204sys.hits = [];
205for i = 1:N
206 sys.hits(i).t = 0;
207 sys.hits(i).beads = i + 1;
208 sys.hits(i).velocity = [0 0 1e-5 * sin(${p} * pi * i / (N + 1))];
209end
212sys.listen.dir = [0 0 1];
213`).replace("sys.listen = struct('beads', 4);", "")
214 );
216 it("gives each mode the decay rate of the chosen law", () => {
217 const s = [Math.log(1000) / 0.5, Math.log(1000) / 0.1];
218 const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 1000) ** 2];
219 const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
220 const sigma0 = s[0] - sigma1 * w2[0];
221 const { sys } = shaped(1);
222 const modes = computeModes(sys, numblEig);
223 const md = modalDamping(sys, modes);
224 for (let k = 0; k < modes.freqs.length; k++) {
225 const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[k]) ** 2;
226 expect(Math.abs(md[k] - expected) / expected).toBeLessThan(1e-9);
227 }
228 for (const p of [1, 3]) {
229 const { sys: sp } = shaped(p);
230 const r = simulate({ ...sp, duration: 0.3, frameRate: 1000 });
231 // energy decays at 2 sigma; compare energy at 0.05 s and 0.25 s (averaged over a period)
232 const avg = (t: number) => {
233 let e = 0;
234 let c = 0;
235 for (let i = 0; i < r.frameTimes.length; i++) if (Math.abs(r.frameTimes[i] - t) < 0.01) (e += r.energy[i]), c++;
236 return e / c;
237 };
238 const rate = Math.log(avg(0.05) / avg(0.25)) / 0.2 / 2;
239 const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[p - 1]) ** 2;
240 expect(Math.abs(rate - expected) / expected).toBeLessThan(0.02);
241 }
242 });
244 it("rejects a law the dashpots cannot make", () => {
247});
249describe("default damping", () => {
250 it("makes every mode fall 60 dB in 2.5 s when the script sets none", () => {
251 const { sys, warnings } = build(STRING(""));
252 expect(warnings.join(" ")).toMatch(/default/);
253 const modes = computeModes(sys, numblEig);
254 for (const s of modalDamping(sys, modes)) expect(Math.log(1000) / s).toBeCloseTo(2.5, 9);
255 });
257 it("is off when the script sets friction = 0", () => {
258 const { sys } = build(STRING("sys.friction = 0;"));
259 expect(Array.from(sys.friction).every(c => c === 0)).toBe(true);
260 });
261});
263describe("bending", () => {
264 // Stiff string: tension plus bending triplets. With theta = k pi/(N+1),
265 // omega_k^2 = [g (2 - 2 cos theta) + kb (2 - 2 cos theta)^2] / m exactly.
266 const STIFF = (extra: string) =>
267 STRING(`
268[pos, E, trip] = bp_chain(N + 2, (N + 1) * a);
269sys.bends = [trip, 50 * ones(size(trip, 1), 1)];
270${extra}`);
272 it("gives the stiff string its exact modes", () => {
273 const { sys } = build(STIFF(""));
274 const modes = computeModes(sys, numblEig);
275 const g = 10 / 0.01;
276 for (let k = 1; k <= 15; k++) {
277 const c = 2 - 2 * Math.cos((k * Math.PI) / 16);
278 const fk = Math.sqrt((g * c + 50 * c * c) / 1e-3) / (2 * Math.PI);
279 expect(modes.freqs[k - 1] / fk).toBeCloseTo(1, 10);
280 }
281 });
283 it("conserves energy to second order", () => {
284 const { sys } = build(STIFF("sys.hits.velocity = [0 0.3 2]; sys.friction = 0;"));
285 const d1 = simulate({ ...sys, dt: 1 / 48000 / 8, duration: 0.05 }).energyDrift;
286 const d2 = simulate({ ...sys, dt: 1 / 48000 / 16, duration: 0.05 }).energyDrift;
287 expect(d1).toBeLessThan(1e-3);
288 expect(d1 / d2).toBeGreaterThan(3);
289 expect(d1 / d2).toBeLessThan(5);
290 });
292 it("follows the bp_decay law with bending damping", () => {
295 const md = modalDamping(sys, modes);
296 const s = [Math.log(1000) / 0.5, Math.log(1000) / 0.1];
297 const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 1000) ** 2];
298 const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
299 const sigma0 = s[0] - sigma1 * w2[0];
300 for (let k = 0; k < modes.freqs.length; k++) {
301 const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[k]) ** 2;
302 expect(md[k] / expected).toBeCloseTo(1, 9);
303 }
304 // and the solver's damping step agrees: energy balance still exact
305 const r1 = simulate({ ...sys, hits: [{ ...sys.hits[0], velocity: [0, 0.3, 2] }], dt: 1 / 48000 / 8, duration: 0.05 });
306 expect(r1.energyDrift).toBeLessThan(1e-3);
307 });
309 it("leaves a straight chain in equilibrium after a rigid rotation", () => {
310 const { sys } = build(
311 STIFF("R3 = [cos(0.7) -sin(0.7) 0; sin(0.7) cos(0.7) 0; 0 0 1] * [1 0 0; 0 cos(0.4) -sin(0.4); 0 sin(0.4) cos(0.4)]; sys.pos = pos * R3';")
312 );
313 const r = simulate({ ...sys, duration: 0.001 });
314 expect(r.residualForce).toBeLessThan(1e-9);
315 });
316});
318describe("free objects", () => {
319 it("bp_anchor holds a bead with the chosen out-of-plane stiffness, at equilibrium", () => {
320 const { sys } = build(`
321[pos, springs, anchors] = bp_anchor([0 0 0], 1, 400, 0.02);
322sys.pos = pos; sys.mass = 1e-3; sys.pinned = anchors; sys.springs = springs;
323hit.beads = 1; hit.velocity = [0 0 1e-4]; sys.hits = hit; sys.listen = struct('beads', 1);
324sys.friction = 0; sys.duration = 0.01;`);
325 const modes = computeModes(sys, numblEig);
326 expect(modes.freqs[0]).toBeCloseTo(Math.sqrt(400 / 1e-3) / (2 * Math.PI), 8);
327 expect(simulate(sys).residualForce).toBeLessThan(1e-9);
328 });
330 it("bar example has the overtones of a free bar", () => {
331 const { sys } = build(readFileSync("matlab/examples/bar.m", "utf8"));
332 const f = computeModes(sys, numblEig).freqs;
333 // two low modes are the bar bouncing on its supports
334 expect(f[1]).toBeLessThan(100);
335 const f1 = f[2];
336 [2.756, 5.404, 8.933].forEach((ratio, q) => expect(Math.abs(f[3 + q] / f1 / ratio - 1)).toBeLessThan(0.01));
337 });
338});
340describe("listener quantity", () => {
341 // Stiff, heavily damped string (dashpots and bending damping), listened to
342 // as displacement, velocity and acceleration at the same bead.
343 const run = (q: string) => {
344 const { sys } = build(
345 STRING(`
346[pos, E, trip] = bp_chain(N + 2, (N + 1) * a);
347sys.bends = [trip, 50 * ones(size(trip, 1), 1)];
350sys.listen.quantity = '${q}';`)
351 );
352 return simulate({ ...sys, duration: 0.05 });
353 };
355 it("defaults to acceleration", () => {
356 expect(build(STRING("")).sys.listenQuantity).toBe("acceleration");
357 });
359 it("gives acceleration and velocity consistent with each other and the displacement", () => {
360 const x = run("displacement").audio;
361 const v = run("velocity").audio;
362 const acc = run("acceleration").audio;
363 const fs = 48000;
364 const rel = (p: (i: number) => number, ref: Float32Array) => {
365 let num = 0;
366 let den = 0;
367 for (let i = 1; i < ref.length - 1; i++) {
368 num += (p(i) - ref[i]) ** 2;
369 den += ref[i] ** 2;
370 }
371 return Math.sqrt(num / den);
372 };
373 // central differences of the recorded signals; their own error is second
374 // order in 1/fs (about 5e-4 for v -> a here, and 1e-5 for x -> v, which is
375 // limited by the single-precision audio)
376 expect(rel(i => ((x[i + 1] - x[i - 1]) * fs) / 2, v)).toBeLessThan(1e-4);
377 expect(rel(i => ((v[i + 1] - v[i - 1]) * fs) / 2, acc)).toBeLessThan(1e-3);
378 });
380 it("rejects an unknown quantity", () => {
381 expect(() => build(STRING("sys.listen.quantity = 'jerk';"))).toThrow(/quantity/);
382 });
383});
385describe("hinges", () => {
386 // A small dome: hexagonal lattice lifted to a paraboloid, clamped at the top,
387 // with stretching springs (no tension) and hinges across interior edges.
388 const DOME = (extra: string, H = 0.01) => `
389R = 0.06; a = 0.012; H = ${H}; m = 1e-3; kh = 0.5; k = 2e4;
390[pos, tri] = bp_hexdisk(R, a);
391pos(:, 3) = H * (1 - sum(pos(:, 1:2).^2, 2) / R^2);
392E = bp_edges(tri);
393sys.pos = pos; sys.mass = m;
394sys.pinned = bp_near(pos, [0 0 H], 1.05 * a);
395sys.springs = bp_springs(pos, E, k, 0);
396hh = bp_hinges(tri);
397sys.hinges = [hh, kh * ones(size(hh, 1), 1)];
398hit.beads = bp_near(pos, [0.04 0 0], 0); hit.velocity = [0 0 1e-3]; sys.hits = hit;
399sys.listen = struct('beads', bp_near(pos, [-0.03 0.02 0], 0));
400sys.friction = 0; sys.duration = 0.05;
401${extra}`;
403 it("bp_hinges lists every interior edge once, with the two opposite beads", () => {
404 const r = executeCode(
405 `[pos, tri] = bp_hexdisk(0.06, 0.012); E = bp_edges(tri); H = bp_hinges(tri);
406 nint = size(E, 1) - sum(accumarray(sort([tri(:,[1 2]); tri(:,[2 3]); tri(:,[3 1])], 2) * [1e4; 1], 1) == 1);`,
407 {},
408 helpers,
409 "t.m"
410 );
411 const H = r.variableValues.H as unknown as { data: Float64Array; shape: number[] };
412 const tri = r.variableValues.tri as unknown as { data: Float64Array; shape: number[] };
413 const nh = H.shape[0];
414 const nt = tri.shape[0];
415 const tris = new Set<string>();
416 for (let t = 0; t < nt; t++) tris.add([0, 1, 2].map(c => tri.data[t + c * nt]).sort((p, q) => p - q).join(","));
417 for (let h = 0; h < nh; h++) {
418 const [i, j, k, l] = [0, 1, 2, 3].map(c => H.data[h + c * nh]);
419 expect(tris.has([i, j, k].sort((p, q) => p - q).join(","))).toBe(true);
420 expect(tris.has([l, j, k].sort((p, q) => p - q).join(","))).toBe(true);
421 }
422 expect(nh).toBe(r.variableValues.nint);
423 });
425 it("leaves a dome, and a rotated dome, unstressed at rest", () => {
426 for (const rot of ["", "Q = [cos(0.7) 0 sin(0.7); 0 1 0; -sin(0.7) 0 cos(0.7)]; sys.pos = sys.pos * Q';"]) {
427 const { sys } = build(DOME(rot));
428 const r = simulate({ ...sys, duration: 0.001 });
429 expect(r.residualForce).toBeLessThan(1e-9);
430 expect(r.warnings).toEqual([]);
431 }
432 });
434 it("conserves energy to second order", () => {
435 const { sys } = build(DOME("sys.hits.velocity = [1 0.5 3];"));
436 const d1 = simulate({ ...sys, dt: 1 / 48000 / 8 }).energyDrift;
437 const d2 = simulate({ ...sys, dt: 1 / 48000 / 16 }).energyDrift;
438 expect(d1).toBeLessThan(1e-3);
439 expect(d1 / d2).toBeGreaterThan(3);
440 expect(d1 / d2).toBeLessThan(5);
441 });
443 it("gives a flat net the same transverse modes in the fast and the full analysis", () => {
444 const { sys } = build(DOME("", 0));
445 const fz = computeModes(sys, numblEig, "transverse").freqs;
446 const full = Array.from(computeModes(sys, numblEig, "full").freqs);
447 for (const f of fz) expect(Math.min(...full.map(g => Math.abs(g / f - 1)))).toBeLessThan(1e-8);
448 });
450 it("gives every 3D mode of a dome the bp_decay law", () => {
451 const { sys } = build(DOME("sys = bp_decay(sys, [100 0.4; 3000 0.1]);"));
452 const modes = computeModes(sys, numblEig, "full");
453 const sigma = modalDamping(sys, modes);
454 const s = [Math.log(1000) / 0.4, Math.log(1000) / 0.1];
455 const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 3000) ** 2];
456 const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
457 const sigma0 = s[0] - sigma1 * w2[0];
458 modes.freqs.forEach((f, k) => expect(sigma[k] / (sigma0 + sigma1 * (2 * Math.PI * f) ** 2)).toBeCloseTo(1, 8));
459 });
461 it.each([6, -3])("keeps a dome's energy in one 3D mode (%i), decaying at the predicted rate", kSel => {
462 const { sys } = build(DOME("sys = bp_decay(sys, [100 0.4; 3000 0.1]);"));
463 const modes = computeModes(sys, numblEig, "full");
464 const sigma = modalDamping(sys, modes);
465 const k = kSel >= 0 ? kSel : modes.freqs.length + kSel; // -3: one of the highest, stretching modes
466 const nd = modes.dof.length;
467 // initial velocity in the shape of mode k: v = eps * psi_k / sqrt(m)
468 const v = new Float64Array(3 * sys.n);
469 for (let r = 0; r < nd; r++) v[modes.dof[r]] = (1e-4 * modes.psi[k * nd + r]) / Math.sqrt(sys.mass[(modes.dof[r] / 3) | 0]);
470 const hits = Array.from({ length: sys.n }, (_, b) => ({
471 t: 0,
472 beads: Int32Array.of(b),
473 velocity: [v[3 * b], v[3 * b + 1], v[3 * b + 2]] as [number, number, number],
474 force: [0, 0, 0] as [number, number, number],
475 duration: 0,
476 })).filter(h => !sys.pinned[h.beads[0]]);
477 const r = simulate({ ...sys, hits, duration: 0.1, frameRate: 1000, dt: 1 / 48000 / 8 }, { modes });
478 const nm = modes.freqs.length;
479 const at = (t: number) => r.frameTimes.findIndex(x => x >= t);
480 const share = (f: number) => r.modeEnergy![f * nm + k] / r.modeEnergy!.slice(f * nm, f * nm + nm).reduce((p, q) => p + q, 0);
481 expect(share(at(0.05))).toBeGreaterThan(0.999);
482 const e = (f: number) => r.modeEnergy![f * nm + k];
483 const rate = Math.log(e(at(0.02)) / e(at(0.08))) / (r.frameTimes[at(0.08)] - r.frameTimes[at(0.02)]) / 2;
484 expect(Math.abs(rate / sigma[k] - 1)).toBeLessThan(0.02);
485 });
486});
↑↓ move↵openescclose