1import { readFileSync, readdirSync } from "fs";
2import { describe, expect, it } from "vitest";
3import { computeModes, modalDamping } from "../src/sim/modes";
4import { numblEig, runScript } from "../src/sim/runScript";
5import { growthRate } from "../src/sim/analysis";
6import { frequencyErrorCents, simulate, stabilityLimit } from "../src/sim/solver";
7import { runTwin } from "../src/sim/twin";
8import type { BeadSystem } from "../src/sim/types";
10const helpers = readdirSync("matlab/helpers").map(name => ({
11 name,
12 source: readFileSync(`matlab/helpers/${name}`, "utf8"),
13}));
14const build = (src: string) => runScript(src, helpers, () => {});
16// A string: N free beads between two pinned ends, tension f0, spacing a.
17// Transverse modes: omega_k = 2 sqrt(f0/(m a)) sin(k pi / (2(N+1))).
18const STRING = (extra: string) => `
19N = 15; a = 0.01; m = 1e-3; f0 = 10; k = 1e4;
20[pos, E] = bp_chain(N + 2, (N + 1) * a);
21sys.pos = pos;
22sys.mass = m;
23sys.pinned = [1, N + 2];
24sys.springs = bp_springs(pos, E, k, f0);
25hit.beads = 5; hit.velocity = [0 0 1e-4];
26sys.hits = hit;
27sys.listen = struct('beads', 4);
28sys.duration = 0.2;
29${extra}
30`;
32function peakFrequency(x: Float32Array, fs: number, fLo: number, fHi: number): number {
33 // direct DFT magnitude on a fine grid, refined by a parabola
34 let best = 0;
35 let bestF = 0;
36 const mag = (f: number) => {
37 let re = 0;
38 let im = 0;
39 for (let i = 0; i < x.length; i++) {
40 const w = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (x.length - 1));
41 re += w * x[i] * Math.cos((2 * Math.PI * f * i) / fs);
42 im -= w * x[i] * Math.sin((2 * Math.PI * f * i) / fs);
43 }
44 return Math.hypot(re, im);
45 };
46 for (let f = fLo; f <= fHi; f += 0.5) {
47 const m = mag(f);
48 if (m > best) {
49 best = m;
50 bestF = f;
51 }
52 }
53 const [a, b, c] = [mag(bestF - 0.5), best, mag(bestF + 0.5)];
54 return bestF + (0.5 * 0.5 * (a - c)) / (a - 2 * b + c);
55}
57describe("string", () => {
58 const { sys } = build(STRING(""));
59 const f1 = (2 * Math.sqrt(10 / (1e-3 * 0.01)) * Math.sin(Math.PI / 32)) / (2 * Math.PI);
61 it("has the analytic transverse modes", () => {
62 const modes = computeModes(sys, numblEig);
63 for (let k = 1; k <= 15; k++) {
64 const fk = (2 * Math.sqrt(10 / (1e-3 * 0.01)) * Math.sin((k * Math.PI) / 32)) / (2 * Math.PI);
65 expect(modes.freqs[k - 1]).toBeCloseTo(fk, 6);
66 }
67 });
69 it("rings at the fundamental for a small strike", () => {
70 const s: BeadSystem = { ...sys, hits: [{ ...sys.hits[0], beads: Int32Array.of(8) }], listenBeads: Int32Array.of(8) };
71 const r = simulate(s);
72 expect(r.aborted).toBeNull();
73 const f = peakFrequency(r.audio, r.sampleRate, f1 * 0.8, f1 * 1.2);
74 expect(Math.abs(f - f1) / f1).toBeLessThan(2e-3);
75 });
77 it("is at equilibrium at rest", () => {
78 const r = simulate({ ...sys, duration: 0.001 });
79 expect(r.residualForce).toBeLessThan(1e-9);
80 expect(r.warnings).toEqual([]);
81 });
82});
84describe("energy", () => {
85 const { sys } = build(STRING("sys.hits.velocity = [0 0.3 2]; sys.friction = 0;"));
86 const drift = (method: "verlet" | "yoshida4", dt: number) =>
87 simulate({ ...sys, method, dt, duration: 0.05 }).energyDrift;
89 it("is conserved to second order with Verlet", () => {
90 const d1 = drift("verlet", 1 / 48000 / 8);
91 const d2 = drift("verlet", 1 / 48000 / 16);
92 expect(d1).toBeLessThan(1e-3);
93 expect(d1 / d2).toBeGreaterThan(3);
94 expect(d1 / d2).toBeLessThan(5);
95 });
97 it("is conserved to fourth order with Yoshida", () => {
98 const d1 = drift("yoshida4", 1 / 48000 / 8);
99 const d2 = drift("yoshida4", 1 / 48000 / 16);
100 expect(d1 / d2).toBeGreaterThan(12);
101 });
103 it("balances with the dissipated energy when damped", () => {
104 const { sys: damped } = build(
105 STRING("sys.hits.velocity = [0 0.3 2]; sys.friction = 2e-3; sys.dashpots = [E, 1e-3*ones(size(E,1),1)];")
106 );
107 // The damping step's energy loss is counted exactly, so E + dissipated
108 // drifts only by the conservative integrator's O(h^2) error.
109 const r1 = simulate({ ...damped, duration: 0.1, dt: 1 / 48000 / 8 });
110 const r2 = simulate({ ...damped, duration: 0.1, dt: 1 / 48000 / 16 });
111 expect(r1.dissipated[r1.dissipated.length - 1]).toBeGreaterThan(0.3 * r1.energy[0]);
112 expect(r1.energyDrift).toBeLessThan(1e-4);
113 expect(r1.energyDrift / r2.energyDrift).toBeGreaterThan(3);
114 });
115});
117describe("damping", () => {
118 it("matches friction when a dashpot ties a bead to a pinned bead", () => {
119 // one free bead between two pinned beads: friction c0 on it is the same
120 // force as two dashpots of c0/2 to its pinned neighbours
121 const base = STRING("").replace("N = 15", "N = 1");
122 const { sys: a } = build(base + "sys.friction = [0; 0.02; 0]; sys.hits.beads = 2; sys.listen.beads = 2;");
123 const { sys: b } = build(base + "sys.dashpots = [1 2 0.01; 2 3 0.01]; sys.hits.beads = 2; sys.listen.beads = 2;");
124 const ra = simulate({ ...a, duration: 0.05 });
125 const rb = simulate({ ...b, duration: 0.05 });
126 let err = 0;
127 let peak = 0;
128 for (let i = 0; i < ra.audio.length; i++) {
129 err = Math.max(err, Math.abs(ra.audio[i] - rb.audio[i]));
130 peak = Math.max(peak, Math.abs(ra.audio[i]));
131 }
132 expect(err / peak).toBeLessThan(1e-6);
133 });
135 it("decays a mode at the predicted rate", () => {
136 // one free bead, friction c0: amplitude decays as exp(-c0/(2m) t)
137 const base = STRING("").replace("N = 15", "N = 1");
138 const { sys } = build(base + "sys.friction = 0.02; sys.hits.beads = 2; sys.listen.beads = 2;");
139 const r = simulate({ ...sys, duration: 0.2 });
140 const rms = (from: number, to: number) => {
141 let s = 0;
142 for (let i = from; i < to; i++) s += r.audio[i] ** 2;
143 return Math.sqrt(s / (to - from));
144 };
145 const fs = r.sampleRate;
146 const ratio = rms(0.15 * fs, 0.2 * fs) / rms(0.05 * fs, 0.1 * fs);
147 const expected = Math.exp((-0.02 / (2 * 1e-3)) * 0.1);
148 expect(Math.abs(ratio - expected) / expected).toBeLessThan(0.02);
149 });
150});
152describe("integrator", () => {
153 it("reports Verlet's known frequency error", () => {
154 // Verlet: omega_num h = 2 asin(omega h / 2)
155 const w = 2 * Math.PI * 5000;
156 const h = 1e-5;
157 const exact = 1200 * Math.log2((2 * Math.asin((w * h) / 2)) / (w * h));
158 expect(frequencyErrorCents("verlet", w, h)).toBeCloseTo(exact, 9);
159 });
161 it("has the expected stability limits", () => {
162 expect(stabilityLimit("verlet")).toBe(2);
163 expect(stabilityLimit("yoshida4")).toBeGreaterThan(1.5);
164 expect(stabilityLimit("yoshida4")).toBeLessThan(1.6);
165 });
166});
168describe("script errors", () => {
169 it("names the field and bead", () => {
170 expect(() => build(STRING("sys.pinned = [1 99];"))).toThrow(/sys.pinned: bead index 99/);
171 expect(() => build(STRING("sys.springs = [1 2 3];"))).toThrow(/4 or 5 columns/);
172 });
173});
175describe("chaos example", () => {
176 const src = readFileSync("matlab/examples/chaos.m", "utf8");
177 const rate = (strength: number) => {
178 const { sys } = build(src.replace(/^s = 1;/m, `s = ${strength};`));
179 const s1 = { ...sys, duration: 1, frameRate: 1000 };
180 const saved: Float64Array[] = [];
181 const r = simulate(s1, { onFrame: (_f, _t, x) => saved.push(x.slice()) });
182 const { sep } = runTwin(s1, saved, r.dt, 1e-9);
183 return growthRate(r.frameTimes, sep, 1e-9);
184 };
186 it("is chaotic when struck hard", () => {
187 const g = rate(1);
188 expect(g).not.toBeNull();
189 expect(g!).toBeGreaterThan(50); // e-folding time under 20 ms
190 });
192 it("is regular when struck softly", () => {
193 expect(rate(0.025)).toBeNull();
194 });
195});
197describe("bp_decay", () => {
198 // String struck in the shape of mode p; its amplitude should decay at
199 // sigma0 + sigma1 w_p^2, the law through T60 = 0.5 s at 100 Hz and 0.1 s at 1 kHz.
200 const shaped = (p: number) =>
201 build(
202 STRING(`
203sys.hits = [];
204for i = 1:N
205 sys.hits(i).t = 0;
206 sys.hits(i).beads = i + 1;
207 sys.hits(i).velocity = [0 0 1e-5 * sin(${p} * pi * i / (N + 1))];
208end
209[sys.friction, sys.dashpots] = bp_decay(pos, m, sys.springs, [100 0.5; 1000 0.1]);
210sys.listen.beads = 2:N+1;
211sys.listen.dir = [0 0 1];
212`).replace("sys.listen = struct('beads', 4);", "")
213 );
215 it("gives each mode the decay rate of the chosen law", () => {
216 const s = [Math.log(1000) / 0.5, Math.log(1000) / 0.1];
217 const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 1000) ** 2];
218 const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
219 const sigma0 = s[0] - sigma1 * w2[0];
220 const { sys } = shaped(1);
221 const modes = computeModes(sys, numblEig);
222 const md = modalDamping(sys, modes);
223 for (let k = 0; k < modes.freqs.length; k++) {
224 const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[k]) ** 2;
225 expect(Math.abs(md[k] - expected) / expected).toBeLessThan(1e-9);
226 }
227 for (const p of [1, 3]) {
228 const { sys: sp } = shaped(p);
229 const r = simulate({ ...sp, duration: 0.3, frameRate: 1000 });
230 // energy decays at 2 sigma; compare energy at 0.05 s and 0.25 s (averaged over a period)
231 const avg = (t: number) => {
232 let e = 0;
233 let c = 0;
234 for (let i = 0; i < r.frameTimes.length; i++) if (Math.abs(r.frameTimes[i] - t) < 0.01) (e += r.energy[i]), c++;
235 return e / c;
236 };
237 const rate = Math.log(avg(0.05) / avg(0.25)) / 0.2 / 2;
238 const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[p - 1]) ** 2;
239 expect(Math.abs(rate - expected) / expected).toBeLessThan(0.02);
240 }
241 });
243 it("rejects a law the dashpots cannot make", () => {
244 expect(() => build(STRING("[f, d] = bp_decay(pos, m, sys.springs, [100 0.1; 1000 0.5]);"))).toThrow(/must not ring longer/);
245 });
246});
248describe("default damping", () => {
249 it("makes every mode fall 60 dB in 2.5 s when the script sets none", () => {
250 const { sys, warnings } = build(STRING(""));
251 expect(warnings.join(" ")).toMatch(/default/);
252 const modes = computeModes(sys, numblEig);
253 for (const s of modalDamping(sys, modes)) expect(Math.log(1000) / s).toBeCloseTo(2.5, 9);
254 });
256 it("is off when the script sets friction = 0", () => {
257 const { sys } = build(STRING("sys.friction = 0;"));
258 expect(Array.from(sys.friction).every(c => c === 0)).toBe(true);
259 });
260});
262describe("bending", () => {
263 // Stiff string: tension plus bending triplets. With theta = k pi/(N+1),
264 // omega_k^2 = [g (2 - 2 cos theta) + kb (2 - 2 cos theta)^2] / m exactly.
265 const STIFF = (extra: string) =>
266 STRING(`
267[pos, E, trip] = bp_chain(N + 2, (N + 1) * a);
268sys.bends = [trip, 50 * ones(size(trip, 1), 1)];
269${extra}`);
271 it("gives the stiff string its exact modes", () => {
272 const { sys } = build(STIFF(""));
273 const modes = computeModes(sys, numblEig);
274 const g = 10 / 0.01;
275 for (let k = 1; k <= 15; k++) {
276 const c = 2 - 2 * Math.cos((k * Math.PI) / 16);
277 const fk = Math.sqrt((g * c + 50 * c * c) / 1e-3) / (2 * Math.PI);
278 expect(modes.freqs[k - 1] / fk).toBeCloseTo(1, 10);
279 }
280 });
282 it("conserves energy to second order", () => {
283 const { sys } = build(STIFF("sys.hits.velocity = [0 0.3 2]; sys.friction = 0;"));
284 const d1 = simulate({ ...sys, dt: 1 / 48000 / 8, duration: 0.05 }).energyDrift;
285 const d2 = simulate({ ...sys, dt: 1 / 48000 / 16, duration: 0.05 }).energyDrift;
286 expect(d1).toBeLessThan(1e-3);
287 expect(d1 / d2).toBeGreaterThan(3);
288 expect(d1 / d2).toBeLessThan(5);
289 });
291 it("follows the bp_decay law with bending damping", () => {
292 const { sys } = build(STIFF("[sys.friction, sys.dashpots, sys.bends] = bp_decay(pos, m, sys.springs, [100 0.5; 1000 0.1], sys.bends);"));
293 const modes = computeModes(sys, numblEig);
294 const md = modalDamping(sys, modes);
295 const s = [Math.log(1000) / 0.5, Math.log(1000) / 0.1];
296 const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 1000) ** 2];
297 const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
298 const sigma0 = s[0] - sigma1 * w2[0];
299 for (let k = 0; k < modes.freqs.length; k++) {
300 const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[k]) ** 2;
301 expect(md[k] / expected).toBeCloseTo(1, 9);
302 }
303 // and the solver's damping step agrees: energy balance still exact
304 const r1 = simulate({ ...sys, hits: [{ ...sys.hits[0], velocity: [0, 0.3, 2] }], dt: 1 / 48000 / 8, duration: 0.05 });
305 expect(r1.energyDrift).toBeLessThan(1e-3);
306 });
308 it("leaves a straight chain in equilibrium after a rigid rotation", () => {
309 const { sys } = build(
310 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';")
311 );
312 const r = simulate({ ...sys, duration: 0.001 });
313 expect(r.residualForce).toBeLessThan(1e-9);
314 });
315});
317describe("free objects", () => {
318 it("bp_anchor holds a bead with the chosen out-of-plane stiffness, at equilibrium", () => {
319 const { sys } = build(`
320[pos, springs, anchors] = bp_anchor([0 0 0], 1, 400, 0.02);
321sys.pos = pos; sys.mass = 1e-3; sys.pinned = anchors; sys.springs = springs;
322hit.beads = 1; hit.velocity = [0 0 1e-4]; sys.hits = hit; sys.listen = struct('beads', 1);
323sys.friction = 0; sys.duration = 0.01;`);
324 const modes = computeModes(sys, numblEig);
325 expect(modes.freqs[0]).toBeCloseTo(Math.sqrt(400 / 1e-3) / (2 * Math.PI), 8);
326 expect(simulate(sys).residualForce).toBeLessThan(1e-9);
327 });
329 it("bar example has the overtones of a free bar", () => {
330 const { sys } = build(readFileSync("matlab/examples/bar.m", "utf8"));
331 const f = computeModes(sys, numblEig).freqs;
332 // two low modes are the bar bouncing on its supports
333 expect(f[1]).toBeLessThan(100);
334 const f1 = f[2];
335 [2.756, 5.404, 8.933].forEach((ratio, q) => expect(Math.abs(f[3 + q] / f1 / ratio - 1)).toBeLessThan(0.01));
336 });
337});