1import { readFileSync, readdirSync } from "fs";
2import { describe, expect, it } from "vitest";
3import { computeModes } from "../src/sim/modes";
4import { numblEig, runScript } from "../src/sim/runScript";
5import { frequencyErrorCents, simulate, stabilityLimit } from "../src/sim/solver";
6import type { BeadSystem } from "../src/sim/types";
8const helpers = readdirSync("matlab/helpers").map(name => ({
9 name,
10 source: readFileSync(`matlab/helpers/${name}`, "utf8"),
11}));
12const build = (src: string) => runScript(src, helpers, () => {});
14// A string: N free beads between two pinned ends, tension f0, spacing a.
15// Transverse modes: omega_k = 2 sqrt(f0/(m a)) sin(k pi / (2(N+1))).
16const STRING = (extra: string) => `
17N = 15; a = 0.01; m = 1e-3; f0 = 10; k = 1e4;
18[pos, E] = bp_chain(N + 2, (N + 1) * a);
19sys.pos = pos;
20sys.mass = m;
21sys.pinned = [1, N + 2];
22sys.springs = bp_springs(pos, E, k, f0);
23hit.beads = 5; hit.velocity = [0 0 1e-4];
24sys.hits = hit;
25sys.listen = struct('beads', 4);
26sys.duration = 0.2;
27${extra}
28`;
30function peakFrequency(x: Float32Array, fs: number, fLo: number, fHi: number): number {
31 // direct DFT magnitude on a fine grid, refined by a parabola
32 let best = 0;
33 let bestF = 0;
34 const mag = (f: number) => {
35 let re = 0;
36 let im = 0;
37 for (let i = 0; i < x.length; i++) {
38 const w = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (x.length - 1));
39 re += w * x[i] * Math.cos((2 * Math.PI * f * i) / fs);
40 im -= w * x[i] * Math.sin((2 * Math.PI * f * i) / fs);
41 }
42 return Math.hypot(re, im);
43 };
44 for (let f = fLo; f <= fHi; f += 0.5) {
45 const m = mag(f);
46 if (m > best) {
47 best = m;
48 bestF = f;
49 }
50 }
51 const [a, b, c] = [mag(bestF - 0.5), best, mag(bestF + 0.5)];
52 return bestF + (0.5 * 0.5 * (a - c)) / (a - 2 * b + c);
53}
55describe("string", () => {
56 const { sys } = build(STRING(""));
57 const f1 = (2 * Math.sqrt(10 / (1e-3 * 0.01)) * Math.sin(Math.PI / 32)) / (2 * Math.PI);
59 it("has the analytic transverse modes", () => {
60 const modes = computeModes(sys, numblEig);
61 for (let k = 1; k <= 15; k++) {
62 const fk = (2 * Math.sqrt(10 / (1e-3 * 0.01)) * Math.sin((k * Math.PI) / 32)) / (2 * Math.PI);
63 expect(modes.freqs[k - 1]).toBeCloseTo(fk, 6);
64 }
65 });
67 it("rings at the fundamental for a small strike", () => {
68 const s: BeadSystem = { ...sys, hits: [{ ...sys.hits[0], beads: Int32Array.of(8) }], listenBeads: Int32Array.of(8) };
69 const r = simulate(s);
70 expect(r.aborted).toBeNull();
71 const f = peakFrequency(r.audio, r.sampleRate, f1 * 0.8, f1 * 1.2);
72 expect(Math.abs(f - f1) / f1).toBeLessThan(2e-3);
73 });
75 it("is at equilibrium at rest", () => {
76 const r = simulate({ ...sys, duration: 0.001 });
77 expect(r.residualForce).toBeLessThan(1e-9);
78 expect(r.warnings).toEqual([]);
79 });
80});
82describe("energy", () => {
83 const { sys } = build(STRING("sys.hits.velocity = [0 0.3 2];"));
84 const drift = (method: "verlet" | "yoshida4", dt: number) =>
85 simulate({ ...sys, method, dt, duration: 0.05 }).energyDrift;
87 it("is conserved to second order with Verlet", () => {
88 const d1 = drift("verlet", 1 / 48000 / 8);
89 const d2 = drift("verlet", 1 / 48000 / 16);
90 expect(d1).toBeLessThan(1e-3);
91 expect(d1 / d2).toBeGreaterThan(3);
92 expect(d1 / d2).toBeLessThan(5);
93 });
95 it("is conserved to fourth order with Yoshida", () => {
96 const d1 = drift("yoshida4", 1 / 48000 / 8);
97 const d2 = drift("yoshida4", 1 / 48000 / 16);
98 expect(d1 / d2).toBeGreaterThan(12);
99 });
101 it("balances with the dissipated energy when damped", () => {
102 const { sys: damped } = build(
103 STRING("sys.hits.velocity = [0 0.3 2]; sys.friction = 2e-3; sys.dashpots = [E, 1e-3*ones(size(E,1),1)];")
104 );
105 // The damping step's energy loss is counted exactly, so E + dissipated
106 // drifts only by the conservative integrator's O(h^2) error.
107 const r1 = simulate({ ...damped, duration: 0.1, dt: 1 / 48000 / 8 });
108 const r2 = simulate({ ...damped, duration: 0.1, dt: 1 / 48000 / 16 });
109 expect(r1.dissipated[r1.dissipated.length - 1]).toBeGreaterThan(0.3 * r1.energy[0]);
110 expect(r1.energyDrift).toBeLessThan(1e-4);
111 expect(r1.energyDrift / r2.energyDrift).toBeGreaterThan(3);
112 });
113});
115describe("damping", () => {
116 it("matches friction when a dashpot ties a bead to a pinned bead", () => {
117 // one free bead between two pinned beads: friction c0 on it is the same
118 // force as two dashpots of c0/2 to its pinned neighbours
119 const base = STRING("").replace("N = 15", "N = 1");
120 const { sys: a } = build(base + "sys.friction = [0; 0.02; 0]; sys.hits.beads = 2; sys.listen.beads = 2;");
121 const { sys: b } = build(base + "sys.dashpots = [1 2 0.01; 2 3 0.01]; sys.hits.beads = 2; sys.listen.beads = 2;");
122 const ra = simulate({ ...a, duration: 0.05 });
123 const rb = simulate({ ...b, duration: 0.05 });
124 let err = 0;
125 let peak = 0;
126 for (let i = 0; i < ra.audio.length; i++) {
127 err = Math.max(err, Math.abs(ra.audio[i] - rb.audio[i]));
128 peak = Math.max(peak, Math.abs(ra.audio[i]));
129 }
130 expect(err / peak).toBeLessThan(1e-6);
131 });
133 it("decays a mode at the predicted rate", () => {
134 // one free bead, friction c0: amplitude decays as exp(-c0/(2m) t)
135 const base = STRING("").replace("N = 15", "N = 1");
136 const { sys } = build(base + "sys.friction = 0.02; sys.hits.beads = 2; sys.listen.beads = 2;");
137 const r = simulate({ ...sys, duration: 0.2 });
138 const rms = (from: number, to: number) => {
139 let s = 0;
140 for (let i = from; i < to; i++) s += r.audio[i] ** 2;
141 return Math.sqrt(s / (to - from));
142 };
143 const fs = r.sampleRate;
144 const ratio = rms(0.15 * fs, 0.2 * fs) / rms(0.05 * fs, 0.1 * fs);
145 const expected = Math.exp((-0.02 / (2 * 1e-3)) * 0.1);
146 expect(Math.abs(ratio - expected) / expected).toBeLessThan(0.02);
147 });
148});
150describe("integrator", () => {
151 it("reports Verlet's known frequency error", () => {
152 // Verlet: omega_num h = 2 asin(omega h / 2)
153 const w = 2 * Math.PI * 5000;
154 const h = 1e-5;
155 const exact = 1200 * Math.log2((2 * Math.asin((w * h) / 2)) / (w * h));
156 expect(frequencyErrorCents("verlet", w, h)).toBeCloseTo(exact, 9);
157 });
159 it("has the expected stability limits", () => {
160 expect(stabilityLimit("verlet")).toBe(2);
161 expect(stabilityLimit("yoshida4")).toBeGreaterThan(1.5);
162 expect(stabilityLimit("yoshida4")).toBeLessThan(1.6);
163 });
164});
166describe("script errors", () => {
167 it("names the field and bead", () => {
168 expect(() => build(STRING("sys.pinned = [1 99];"))).toThrow(/sys.pinned: bead index 99/);
169 expect(() => build(STRING("sys.springs = [1 2 3];"))).toThrow(/4 or 5 columns/);
170 });
171});