Bead-and-spring percussion simulator: numbl-scripted setup, Verlet/Yoshida solver with exact damping, linear modes, browser app
34 changed files+6863−0
.claude/launch.jsonadded+6−0View file
@@ -0,0 +1,6 @@
1+{
2+ "version": "0.0.1",
3+ "configurations": [
4+ { "name": "bead-percussion", "runtimeExecutable": "npx", "runtimeArgs": ["vite", "--port", "5191", "--strictPort"], "port": 5191 }
5+ ]
6+}
.gitignoreadded+3−0View file
@@ -0,0 +1,3 @@
1+node_modules
2+dist
3+*.log
README.mdadded+63−0View file
@@ -0,0 +1,63 @@
1+# bead-percussion
2+
3+A browser app for idealized percussion: a network of beads (point masses) joined by springs, set up by a MATLAB-syntax script that runs in the browser with [numbl](https://numbl.org), integrated accurately in time, and heard through the displacement of a listener bead.
4+
5+The aim is to reproduce the phenomena of drums and cymbals (resonant modes, nonlinear coupling between them, transfer of energy to higher frequencies, amplitude-dependent detuning) in a system that is physical but much simpler than the continuum equations of a membrane or a shell. Nothing here discretizes a PDE. The bead network is the system.
6+
7+## Running
8+
9+```bash
10+npm install
11+npm run dev # the app, at the printed localhost URL
12+npm test # solver tests (Node)
13+```
14+
15+## The model
16+
17+Each bead has a mass and a position in 3D. The script defines:
18+
19+- *Springs* between pairs of beads, with stiffness k, rest length L0 and optional cubic stiffness k3: the force along the spring is k e + k3 e^3, with e the current length minus L0. A spring whose rest length is shorter than the distance between its beads pulls with a force at rest, and this pre-tension gives a flat net its transverse stiffness, as tension does for a drum head. Because the springs are real springs in 3D, the geometric nonlinearity (a transverse displacement stretches the springs) is exact rather than a truncated expansion.
20+- *Dashpots* between pairs of beads, with force c (v_j - v_i) on bead i. Short-wavelength motion has larger relative velocities between neighbours, so this damping grows with frequency.
21+- *Friction* on each bead, force -c0 v, which damps all frequencies alike.
22+- *Pinned* beads, held fixed.
23+- *Hits*: an instant velocity given to a set of beads, or a raised-cosine force pulse of given peak and duration. Any number, at any times.
24+- A *listener*: a bead (or the average of several) and a direction; the displacement along that direction is the audio signal.
25+
26+Bending stiffness (for plates and cymbals) is not in yet.
27+
28+The `matlab/helpers/bp_*.m` functions build common geometries (`bp_hexdisk`: a triangular lattice clipped to a disk, with its rim beads; `bp_chain`: a line of beads), list edges (`bp_edges`), make springs with a given tension at rest (`bp_springs`) and find beads near a point (`bp_near`). They are ordinary scripts; the app lists them, and a script can ignore them and build `sys` directly. The full list of `sys` fields is in the app under "sys fields".
29+
30+## Time integration
31+
32+The conservative part (springs and force pulses) is integrated with velocity Verlet, or with Yoshida's fourth-order composition of it (`sys.method = 'yoshida4'`). Damping is linear in the velocities, dv/dt = -A v, and is applied as a separate step v <- exp(-tau A) v, computed to round-off with a scaled Taylor series, once per audio sample (Strang splitting). The splitting error scales with the damping rate times the sample interval, which is small for the damping of interest here.
33+
34+The time step is chosen automatically. We estimate the highest frequency of the system linearized about rest by power iteration, take `sys.safety` (default 0.5) of the method's stability limit, and round so that a whole number of steps fits in each audio sample (48 kHz by default). The in-plane (stretching) modes of a net are much stiffer than its transverse modes, so they, not the audible frequencies, usually set the step; a smaller pre-stretch makes this worse.
35+
36+Each run reports:
37+
38+- the integrator's frequency error at 1 to 20 kHz, in cents (exact for a linear oscillator; with the default step for the drum example, Verlet is about 3.4 cents sharp at 5 kHz and 0.14 cents at 1 kHz);
39+- the energy drift after the last excitation, with the energy removed by damping added back (it is counted exactly);
40+- optionally, the relative difference of the listener signal from a run at dt/2 over [0, t] for several t. In a chaotic run this difference grows with t no matter how small dt is, so the short windows are the meaningful ones.
41+
42+## Linear modes
43+
44+For a net that rests flat (all beads at one z), small transverse motion decouples from in-plane motion, and its stiffness matrix is a graph Laplacian weighted by spring tension over length. The app computes its eigenpairs with numbl's `eig` (up to 1500 free beads; about 1 s for the 721 free beads of the drum example) and records the energy of every mode over time. The "energy by linear mode" plot shows the share of the quadratic energy in log-spaced frequency bands: in a linear system each band only decays, so any band that gains energy is fed by the nonlinearity.
45+
46+The "small-strike reference" option runs the same system with every hit scaled down by 1000 and the output scaled back up, which is what the system would do if it were linear. The app plays it for comparison and reports the difference from the actual run.
47+
48+## Examples
49+
50+| Example | What it shows |
51+| --- | --- |
52+| `drum` | A 15 cm drum head (823 beads, fundamental 237 Hz) struck by a soft mallet. Nearly linear at first (2% from the small-strike reference over 10 ms); the slight pitch change at large amplitude accumulates to 25% by 0.3 s. About 7x real time. |
53+| `cascade` | The same net with a 0.5% pre-stretch, weaker dashpots and a hard strike. The share of energy in modes above 1 kHz grows from 4% (small strike) to about 25% within 30 ms. About 10x real time. |
54+| `string` | 120 beads on a string, struck slightly off-axis; fast. |
55+
56+Timings are for a laptop (Chrome, single thread).
57+
58+## Limitations and next steps
59+
60+- No bending stiffness yet, so no plates or cymbals; the plan is hinge (dihedral angle) elements on pairs of triangles, with a domed rest shape.
61+- The solver is plain TypeScript in one Web Worker. Cost scales with the number of springs times the number of steps; nets of several thousand beads will need WASM or WebGPU.
62+- A bead net has an upper frequency set by its spacing (about 4 kHz for the 1 cm drum example), so a dense high-frequency spectrum needs many beads.
63+- Linear modes are computed only for flat nets.
index.htmladded+12−0View file
@@ -0,0 +1,12 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
6+ <title>bead-percussion</title>
7+ </head>
8+ <body>
9+ <div id="root"></div>
10+ <script type="module" src="/src/main.tsx"></script>
11+ </body>
12+</html>
matlab/examples/cascade.madded+37−0View file
@@ -0,0 +1,37 @@
1+% A loose, lightly damped net, struck hard: energy cascades from the struck
2+% low modes to high ones over tens of milliseconds.
3+%
4+% Compared with drum.m, the springs are stretched only 0.5% at rest, so the
5+% stretching caused by the motion itself is large relative to the rest
6+% tension, and the dashpots are 3x weaker, so high modes keep their energy
7+% long enough to pass it on. Turn on "small-strike reference" to hear and
8+% measure the difference from the same strike at small amplitude.
9+
10+R = 0.15;
11+a = 0.01;
12+sigma = 0.25;
13+T = 2000;
14+strain = 0.005;
15+
16+[pos, tri, rim] = bp_hexdisk(R, a);
17+E = bp_edges(tri);
18+f0 = T * a / sqrt(3);
19+k = f0 / (strain * a);
20+
21+sys.pos = pos;
22+sys.mass = sigma * (sqrt(3)/2) * a^2;
23+sys.pinned = rim;
24+sys.springs = bp_springs(pos, E, k, f0);
25+sys.dashpots = [E, 3e-5 * ones(size(E, 1), 1)];
26+sys.friction = 1e-4;
27+
28+hit.t = 0;
29+hit.beads = bp_near(pos, [0.05 0 0], 0.03);
30+hit.force = [0 0 -12]; % peak force per bead (N); try 1 to 50
31+hit.duration = 1e-3;
32+sys.hits = hit;
33+
34+sys.listen = struct('beads', bp_near(pos, [-0.04 0.03 0], 0), 'dir', [0 0 1]);
35+sys.duration = 1.5;
36+
37+fprintf('%d beads, %d springs, %d struck\n', size(pos, 1), size(E, 1), numel(hit.beads));
matlab/examples/drum.madded+42−0View file
@@ -0,0 +1,42 @@
1+% Drum head: a triangular net of beads under tension, pinned at the rim,
2+% struck by a soft mallet. Units are SI. The net stands in for a membrane
3+% with areal density sigma and tension T; its fundamental is near
4+% 2.405/(2*pi*R)*sqrt(T/sigma).
5+
6+R = 0.15; % radius (m)
7+a = 0.01; % bead spacing (m)
8+sigma = 0.25; % areal density (kg/m^2)
9+T = 2000; % membrane tension (N/m)
10+strain = 0.02; % spring stretch at rest. Smaller is closer to a real
11+ % membrane and more nonlinear, but stiffer in-plane, so
12+ % the time step shrinks.
13+
14+[pos, tri, rim] = bp_hexdisk(R, a);
15+E = bp_edges(tri);
16+f0 = T * a / sqrt(3); % spring tension that gives membrane tension T
17+k = f0 / (strain * a); % spring stiffness
18+
19+sys.pos = pos;
20+sys.mass = sigma * (sqrt(3)/2) * a^2; % lattice area per bead
21+sys.pinned = rim;
22+sys.springs = bp_springs(pos, E, k, f0);
23+sys.dashpots = [E, 1e-4 * ones(size(E, 1), 1)]; % damping that grows with frequency
24+sys.friction = 1e-4; % damping that does not
25+
26+% A soft mallet: a 1 ms raised-cosine push on the beads within 3 cm.
27+hit.t = 0;
28+hit.beads = bp_near(pos, [0.05 0 0], 0.03);
29+hit.force = [0 0 -0.4]; % peak force per bead (N)
30+hit.duration = 1e-3;
31+sys.hits = hit;
32+% For more strikes, make sys.hits a struct array, e.g.
33+% sys.hits(2) = struct('t', 0.5, 'beads', bp_near(pos, [0 0.06 0], 0.02), ...
34+% 'force', [0 0 -0.6], 'duration', 1e-3);
35+% (every hit in the array needs the same fields). A hit may instead give an
36+% instant velocity: struct('t', 0, 'beads', idx, 'velocity', [0 0 -3]).
37+
38+sys.listen = struct('beads', bp_near(pos, [-0.04 0.03 0], 0), 'dir', [0 0 1]);
39+sys.duration = 2;
40+
41+fprintf('%d beads (%d pinned), %d springs\n', size(pos, 1), numel(rim), size(E, 1));
42+fprintf('expected fundamental: %.1f Hz\n', 2.405 / (2*pi*R) * sqrt(T / sigma));
matlab/examples/string.madded+34−0View file
@@ -0,0 +1,34 @@
1+% A beaded string, pinned at both ends and struck near one end. Small and
2+% fast to run. Its transverse modes are
3+% f_k = (1/pi) sqrt(f0/(m a)) sin(k pi / (2(N+1))),
4+% nearly harmonic for k << N. A hard strike stretches the string, raising
5+% the tension and with it the pitch while the amplitude is large; the
6+% string can also swing in y as well as z.
7+
8+N = 120; % free beads
9+L = 0.65; % length (m)
10+mu = 2e-3; % mass per length (kg/m)
11+T = 60; % tension (N)
12+strain = 0.01; % spring stretch at rest
13+
14+[pos, E] = bp_chain(N + 2, L);
15+a = L / (N + 1);
16+k = T / (strain * a);
17+
18+sys.pos = pos;
19+sys.mass = mu * a;
20+sys.pinned = [1, N + 2];
21+sys.springs = bp_springs(pos, E, k, T);
22+sys.dashpots = [E, 2e-4 * ones(size(E, 1), 1)];
23+sys.friction = 2e-5;
24+
25+hit.t = 0;
26+hit.beads = bp_near(pos, [0.12 0 0], 0.01);
27+hit.force = [0 0.2 1] * 0.8; % N per bead; slightly off-axis
28+hit.duration = 1.5e-3;
29+sys.hits = hit;
30+
31+sys.listen = struct('beads', N, 'dir', [0 0 1]);
32+sys.duration = 2.5;
33+
34+fprintf('fundamental (small amplitude): %.1f Hz\n', sqrt(T / mu) / (2 * L));
matlab/helpers/bp_chain.madded+8−0View file
@@ -0,0 +1,8 @@
1+function [pos, E] = bp_chain(n, len)
2+% BP_CHAIN A straight line of n beads along x, from 0 to len.
3+% [pos, E] = bp_chain(n, len) returns positions (n-by-3) and the n-1
4+% neighbor pairs E. The end beads are 1 and n.
5+x = linspace(0, len, n)';
6+pos = [x, zeros(n, 2)];
7+E = [(1:n-1)', (2:n)'];
8+end
matlab/helpers/bp_edges.madded+6−0View file
@@ -0,0 +1,6 @@
1+function E = bp_edges(tri)
2+% BP_EDGES Unique bead pairs (edges) of a triangle list.
3+% E = bp_edges(tri) returns an m-by-2 matrix with E(:,1) < E(:,2).
4+E = sort([tri(:, [1 2]); tri(:, [2 3]); tri(:, [3 1])], 2);
5+E = unique(E, 'rows');
6+end
matlab/helpers/bp_hexdisk.madded+26−0View file
@@ -0,0 +1,26 @@
1+function [pos, tri, rim] = bp_hexdisk(R, a)
2+% BP_HEXDISK Triangular lattice of beads with spacing a, clipped to a disk of radius R.
3+% [pos, tri, rim] = bp_hexdisk(R, a)
4+% pos n-by-3 bead positions in the plane z = 0, centered at the origin
5+% tri triangles of the lattice (rows of three bead indices)
6+% rim beads with fewer than six lattice neighbors (the edge); pin these
7+% for a drum head. Every other bead then sits at an equilibrium
8+% under uniform spring tension.
9+m = ceil(R / a) + 2;
10+[I, J] = ndgrid(-2*m:2*m, -m:m);
11+X = a * (I + J/2);
12+Y = a * (sqrt(3)/2) * J;
13+inside = X.^2 + Y.^2 <= R^2 * (1 + 1e-12);
14+id = zeros(size(I));
15+id(inside) = 1:nnz(inside);
16+pos = [X(inside), Y(inside), zeros(nnz(inside), 1)];
17+A = id(1:end-1, 1:end-1);
18+B = id(2:end, 1:end-1);
19+C = id(1:end-1, 2:end);
20+D = id(2:end, 2:end);
21+tri = [A(:), B(:), C(:); B(:), D(:), C(:)];
22+tri = tri(all(tri > 0, 2), :);
23+E = bp_edges(tri);
24+deg = accumarray(E(:), 1, [size(pos, 1), 1]);
25+rim = find(deg < 6);
26+end
matlab/helpers/bp_near.madded+13−0View file
@@ -0,0 +1,13 @@
1+function idx = bp_near(pos, p, r)
2+% BP_NEAR Beads near a point.
3+% idx = bp_near(pos, p, r) returns the beads within distance r of point p
4+% (1-by-3). With r = 0, or if no bead is that close, returns the nearest bead.
5+if numel(p) == 2
6+ p = [p(:)', 0];
7+end
8+d = vecnorm(pos - p(:)', 2, 2);
9+idx = find(d <= r);
10+if isempty(idx)
11+ [~, idx] = min(d);
12+end
13+end
matlab/helpers/bp_springs.madded+19−0View file
@@ -0,0 +1,19 @@
1+function S = bp_springs(pos, E, k, f0, k3)
2+% BP_SPRINGS Springs along edges E with stiffness k and tension f0 at rest.
3+% S = bp_springs(pos, E, k, f0) returns rows [i j k L0] for sys.springs,
4+% choosing each rest length L0 so the spring pulls with force f0 at the
5+% bead positions pos: L0 = L - f0/k. k and f0 may be scalars or one per edge.
6+% S = bp_springs(pos, E, k, f0, k3) appends a cubic stiffness column.
7+L = vecnorm(pos(E(:,1), :) - pos(E(:,2), :), 2, 2);
8+m = size(E, 1);
9+k = k(:) .* ones(m, 1);
10+f0 = f0(:) .* ones(m, 1);
11+L0 = L - f0 ./ k;
12+if any(L0 <= 0)
13+ error('bp_springs: tension f0 is too large for stiffness k (rest length would be <= 0)');
14+end
15+S = [E, k, L0];
16+if nargin > 4
17+ S = [S, k3(:) .* ones(m, 1)];
18+end
19+end
package-lock.jsonadded+4014−0View file
This diff is 4,019 lines long and is not shown.
package.jsonadded+33−0View file
@@ -0,0 +1,33 @@
1+{
2+ "name": "bead-percussion",
3+ "version": "0.1.0",
4+ "private": true,
5+ "type": "module",
6+ "scripts": {
7+ "dev": "vite",
8+ "build": "tsc --noEmit && vite build",
9+ "preview": "vite preview",
10+ "test": "vitest run",
11+ "bench": "tsx tests/bench.ts"
12+ },
13+ "dependencies": {
14+ "@codemirror/language": "^6.12.4",
15+ "@codemirror/legacy-modes": "^6.5.4",
16+ "@codemirror/state": "^6.7.6",
17+ "@codemirror/view": "^6.43.13",
18+ "codemirror": "^6.0.2",
19+ "numbl": "^0.4.21",
20+ "react": "^19.3.0",
21+ "react-dom": "^19.3.0"
22+ },
23+ "devDependencies": {
24+ "@types/node": "^26.6.4",
25+ "@types/react": "^19.3.0",
26+ "@types/react-dom": "^19.3.0",
27+ "@vitejs/plugin-react": "^6.1.2",
28+ "tsx": "^4.23.15",
29+ "typescript": "^7.0.2",
30+ "vite": "^8.3.3",
31+ "vitest": "^5.0.3"
32+ }
33+}
src/App.tsxadded+367−0View file
@@ -0,0 +1,367 @@
1+import { useCallback, useEffect, useRef, useState } from "react";
2+import { wavBlob } from "./sim/analysis";
3+import { BeadView } from "./ui/BeadView";
4+import { Editor } from "./ui/Editor";
5+import { Energy, ModeEnergy, Spectrogram, Spectrum, Waveform } from "./ui/Plots";
6+import type { FromWorker, RunPayload, RunRequest } from "./worker";
7+
8+const exampleSources = import.meta.glob("../matlab/examples/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
9+const helperSources = import.meta.glob("../matlab/helpers/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
10+const examples = Object.entries(exampleSources)
11+ .map(([path, source]) => ({ name: path.split("/").pop()!.replace(/\.m$/, ""), source }))
12+ .sort((a, b) => (a.name === "drum" ? -1 : b.name === "drum" ? 1 : a.name.localeCompare(b.name)));
13+const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source }));
14+
15+const STORE_KEY = "bead-percussion:script";
16+const load = () => {
17+ try {
18+ return localStorage.getItem(STORE_KEY);
19+ } catch {
20+ return null;
21+ }
22+};
23+const save = (s: string) => {
24+ try {
25+ localStorage.setItem(STORE_KEY, s);
26+ } catch {
27+ /* storage unavailable */
28+ }
29+};
30+
31+export default function App() {
32+ const [script, setScript] = useState(() => load() ?? examples[0].source);
33+ const [running, setRunning] = useState(false);
34+ const [stage, setStage] = useState("");
35+ const [progress, setProgress] = useState<number | null>(null);
36+ const [log, setLog] = useState("");
37+ const [error, setError] = useState<string | null>(null);
38+ const [payload, setPayload] = useState<RunPayload | null>(null);
39+ const [fmax, setFmax] = useState(10000);
40+ const [cursor, setCursor] = useState<number | null>(null);
41+ const [viewTime, setViewTime] = useState<number | null>(null);
42+ const [dtCheck, setDtCheck] = useState(false);
43+ const [linearRef, setLinearRef] = useState(false);
44+ const [playing, setPlaying] = useState<"main" | "linear" | null>(null);
45+ const [runId, setRunId] = useState(0); // remounts panels with per-run state
46+ const worker = useRef<Worker | null>(null);
47+ const audioCtx = useRef<AudioContext | null>(null);
48+ const source = useRef<AudioBufferSourceNode | null>(null);
49+ const startedAt = useRef(0);
50+
51+ useEffect(() => save(script), [script]);
52+
53+ const stopAudio = useCallback(() => {
54+ source.current?.stop();
55+ source.current = null;
56+ setCursor(null);
57+ setPlaying(null);
58+ }, []);
59+
60+ const run = useCallback(
61+ () => {
62+ worker.current?.terminate();
63+ stopAudio();
64+ const w = new Worker(new URL("./worker.ts", import.meta.url), { type: "module" });
65+ worker.current = w;
66+ setRunning(true);
67+ setError(null);
68+ setLog("");
69+ setProgress(null);
70+ setStage("Starting");
71+ w.onmessage = (ev: MessageEvent<FromWorker>) => {
72+ const m = ev.data;
73+ if (m.type === "log") setLog(l => l + m.text);
74+ else if (m.type === "stage") {
75+ setStage(m.text);
76+ setProgress(null);
77+ } else if (m.type === "progress") setProgress(m.t / m.duration);
78+ else if (m.type === "error") {
79+ setError(m.message);
80+ setRunning(false);
81+ w.terminate();
82+ } else if (m.type === "result") {
83+ setPayload(m.payload);
84+ setRunId(i => i + 1);
85+ const fsHalf = m.payload.result.sampleRate / 2;
86+ const top = m.payload.modeFreqs ? m.payload.modeFreqs[m.payload.modeFreqs.length - 1] * 1.15 : 10000;
87+ setFmax([1000, 2000, 5000, 10000, 20000].find(f => f >= top && f <= fsHalf) ?? Math.min(20000, fsHalf));
88+ setRunning(false);
89+ w.terminate();
90+ }
91+ };
92+ w.onerror = e => {
93+ setError(e.message || "Worker error");
94+ setRunning(false);
95+ };
96+ const req: RunRequest = { type: "run", script, dtCheck, linear: linearRef };
97+ w.postMessage(req);
98+ },
99+ [script, stopAudio, dtCheck, linearRef]
100+ );
101+
102+ const stop = () => {
103+ worker.current?.terminate();
104+ worker.current = null;
105+ setRunning(false);
106+ setStage("Stopped");
107+ };
108+
109+ const play = (which: "main" | "linear") => {
110+ if (!payload) return;
111+ stopAudio();
112+ const { sampleRate } = payload.result;
113+ const audio = which === "linear" && payload.linear ? payload.linear.audio : payload.result.audio;
114+ const ctx = audioCtx.current ?? new AudioContext();
115+ audioCtx.current = ctx;
116+ // both signals share the main run's gain, so the A/B comparison keeps levels
117+ let peak = 0;
118+ for (const v of payload.result.audio) peak = Math.max(peak, Math.abs(v));
119+ const buf = ctx.createBuffer(1, audio.length, sampleRate);
120+ const ch = buf.getChannelData(0);
121+ const g = peak > 0 ? 0.9 / peak : 1;
122+ for (let i = 0; i < audio.length; i++) ch[i] = Math.max(-1, Math.min(1, g * audio[i]));
123+ const src = ctx.createBufferSource();
124+ src.buffer = buf;
125+ src.connect(ctx.destination);
126+ src.onended = () => {
127+ if (source.current === src) {
128+ source.current = null;
129+ setCursor(null);
130+ setPlaying(null);
131+ }
132+ };
133+ src.start();
134+ source.current = src;
135+ setPlaying(which);
136+ startedAt.current = ctx.currentTime;
137+ const T = audio.length / sampleRate;
138+ const tick = () => {
139+ if (source.current !== src) return;
140+ const t = ctx.currentTime - startedAt.current;
141+ setCursor(Math.min(t, T));
142+ if (t < T) requestAnimationFrame(tick);
143+ };
144+ requestAnimationFrame(tick);
145+ };
146+
147+ const download = () => {
148+ if (!payload) return;
149+ const url = URL.createObjectURL(wavBlob(payload.result.audio, payload.result.sampleRate));
150+ const a = document.createElement("a");
151+ a.href = url;
152+ a.download = "bead-percussion.wav";
153+ a.click();
154+ setTimeout(() => URL.revokeObjectURL(url), 1000);
155+ };
156+
157+ const r = payload?.result;
158+ return (
159+ <div className="app">
160+ <header>
161+ <h1>bead-percussion</h1>
162+ <span className="sub">beads and springs, set up in numbl, heard through one bead</span>
163+ </header>
164+ <div className="columns">
165+ <section className="left">
166+ <div className="controls">
167+ <select
168+ value=""
169+ onChange={e => {
170+ const ex = examples.find(x => x.name === e.target.value);
171+ if (ex && (script === ex.source || confirm(`Replace the script with the "${ex.name}" example?`))) setScript(ex.source);
172+ }}
173+ >
174+ <option value="" disabled>
175+ Load example…
176+ </option>
177+ {examples.map(ex => (
178+ <option key={ex.name} value={ex.name}>
179+ {ex.name}
180+ </option>
181+ ))}
182+ </select>
183+ <span className="grow" />
184+ {running ? (
185+ <button onClick={stop}>Stop</button>
186+ ) : (
187+ <button className="primary" onClick={run} title="Ctrl/Cmd+Enter">
188+ Run
189+ </button>
190+ )}
191+ </div>
192+ <div className="controls options">
193+ <label title="Run again at dt/2 and compare the listener signals">
194+ <input type="checkbox" checked={dtCheck} onChange={e => setDtCheck(e.target.checked)} /> dt/2 check
195+ </label>
196+ <label title="Run again with every hit scaled by 1e-3 (then scaled back up): what the system would do if it were linear">
197+ <input type="checkbox" checked={linearRef} onChange={e => setLinearRef(e.target.checked)} /> small-strike reference
198+ </label>
199+ </div>
200+ <Editor value={script} onChange={setScript} onRun={() => !running && run()} />
201+ {running && (
202+ <div className="status">
203+ {stage}
204+ {progress !== null && ` ${(100 * progress).toFixed(0)}%`}
205+ <div className="bar">
206+ <div style={{ width: `${100 * (progress ?? 0)}%` }} />
207+ </div>
208+ </div>
209+ )}
210+ {error && <pre className="error">{error}</pre>}
211+ {log && <pre className="log">{log}</pre>}
212+ {payload && <Report p={payload} />}
213+ <details className="helpers">
214+ <summary>Helper functions (bp_*.m)</summary>
215+ {helpers.map(h => (
216+ <div key={h.name}>
217+ <h4>{h.name}</h4>
218+ <pre>{h.source}</pre>
219+ </div>
220+ ))}
221+ </details>
222+ <details className="helpers">
223+ <summary>sys fields</summary>
224+ <pre>{SYS_DOC}</pre>
225+ </details>
226+ </section>
227+ <section className="right">
228+ {!payload && <div className="empty">Run a script to hear and see the result.</div>}
229+ {payload && r && (
230+ <>
231+ <div className="controls">
232+ <button className="primary" onClick={() => (playing === "main" ? stopAudio() : play("main"))}>
233+ {playing === "main" ? "Stop" : "Play"}
234+ </button>
235+ {payload.linear && (
236+ <button onClick={() => (playing === "linear" ? stopAudio() : play("linear"))} title="The same strike scaled down 1000x, then scaled back up">
237+ {playing === "linear" ? "Stop" : "Play small-strike reference"}
238+ </button>
239+ )}
240+ <button onClick={download}>Download WAV</button>
241+ <span className="grow" />
242+ <label>
243+ freq range{" "}
244+ <select value={fmax} onChange={e => setFmax(+e.target.value)}>
245+ {[1000, 2000, 5000, 10000, 20000]
246+ .filter(f => f <= r.sampleRate / 2)
247+ .map(f => (
248+ <option key={f} value={f}>
249+ {f / 1000} kHz
250+ </option>
251+ ))}
252+ </select>
253+ </label>
254+ </div>
255+ <h3>Listener signal</h3>
256+ <Waveform p={payload} cursor={cursor ?? viewTime} />
257+ <h3>Spectrogram</h3>
258+ <Spectrogram p={payload} fmax={fmax} />
259+ <h3>Spectrum, two windows (detuning shows as peaks shifting)</h3>
260+ <Spectrum key={runId} p={payload} fmax={fmax} />
261+ <h3>Energy by linear mode (the cascade)</h3>
262+ <ModeEnergy key={runId} p={payload} />
263+ <h3>Energy</h3>
264+ <Energy p={payload} />
265+ <h3>Beads</h3>
266+ <BeadView key={runId} p={payload} onTime={setViewTime} />
267+ </>
268+ )}
269+ </section>
270+ </div>
271+ </div>
272+ );
273+}
274+
275+function Report({ p }: { p: RunPayload }) {
276+ const r = p.result;
277+ const nPinned = p.sys.pinned.reduce((a, b) => a + b, 0);
278+ const rows: [string, string][] = [
279+ ["beads", `${p.sys.n} (${nPinned} pinned)`],
280+ ["springs, dashpots", `${p.sys.si.length}, ${p.sys.nDashpots}`],
281+ ["method", p.sys.method],
282+ ["dt", `${r.dt.toExponential(3)} s (${r.substeps} per audio sample)`],
283+ ["highest linear frequency", `${(r.omegaMax / (2 * Math.PI)).toFixed(0)} Hz`],
284+ ["frequency error", r.freqErrors.map(e => `${e.cents.toFixed(e.cents < 0.1 ? 3 : 2)}¢ at ${e.hz / 1000} kHz`).join(", ")],
285+ ["energy drift after excitation", r.energyDrift.toExponential(2)],
286+ ["run time", `${(r.elapsedMs / 1000).toFixed(1)} s (${(r.elapsedMs / 1000 / p.sys.duration).toFixed(1)}x real time)`],
287+ ];
288+ if (p.modeFreqs) rows.push(["lowest modes", Array.from(p.modeFreqs.slice(0, 5)).map(f => f.toFixed(1)).join(", ") + " Hz"]);
289+ return (
290+ <div className="report">
291+ {r.aborted && <div className="error">{r.aborted}</div>}
292+ {[...p.notes, ...r.warnings].map((w, i) => (
293+ <div key={i} className="warn">
294+ {w}
295+ </div>
296+ ))}
297+ <table>
298+ <tbody>
299+ {rows.map(([k, v]) => (
300+ <tr key={k}>
301+ <td>{k}</td>
302+ <td>{v}</td>
303+ </tr>
304+ ))}
305+ </tbody>
306+ </table>
307+ {p.linear && (
308+ <>
309+ <div className="caption">
310+ Listener signal vs the small-strike reference (how nonlinear the run is): relative RMS difference over [0, t]
311+ </div>
312+ <table>
313+ <tbody>
314+ {p.linear.windows.map(w => (
315+ <tr key={w.t}>
316+ <td>t = {w.t} s</td>
317+ <td>{w.rel.toExponential(2)}</td>
318+ </tr>
319+ ))}
320+ </tbody>
321+ </table>
322+ </>
323+ )}
324+ {p.check && (
325+ <>
326+ <div className="caption">
327+ Listener signal, dt vs dt/2 ({p.check.dt.toExponential(2)} s): relative RMS difference over [0, t]
328+ </div>
329+ <table>
330+ <tbody>
331+ {p.check.windows.map(w => (
332+ <tr key={w.t}>
333+ <td>t = {w.t} s</td>
334+ <td>{w.rel.toExponential(2)}</td>
335+ </tr>
336+ ))}
337+ </tbody>
338+ </table>
339+ </>
340+ )}
341+ </div>
342+ );
343+}
344+
345+const SYS_DOC = `Required
346+ sys.pos n-by-3 bead positions (m)
347+ sys.mass scalar or n-by-1 (kg)
348+ sys.listen struct: beads (index or indices, averaged), dir (default [0 0 1])
349+Elements
350+ sys.springs rows [i j k L0] or [i j k L0 k3]; force k*e + k3*e^3, e = length - L0
351+ sys.dashpots rows [i j c]; force c*(v_j - v_i) on bead i
352+ sys.friction scalar or n-by-1; force -c0*v
353+ sys.pinned bead indices (or logical mask) held fixed
354+Excitation
355+ sys.hits struct or struct array, each with
356+ beads, t (default 0), and either
357+ velocity [vx vy vz]: added instantly at t, or
358+ force [fx fy fz] + duration: raised-cosine pulse per bead
359+Solver (optional)
360+ sys.duration seconds (default 2)
361+ sys.method 'verlet' (default) or 'yoshida4'
362+ sys.safety fraction of the stability limit for automatic dt (default 0.5)
363+ sys.dt fixed step (rounded to divide the audio sample interval)
364+ sys.sampleRate audio rate (default 48000)
365+ sys.frameRate energy/mode/animation records per second (default 500)
366+ sys.modes compute linear transverse modes (default true; flat systems)
367+Bead indices are 1-based.`;
src/main.tsxadded+5−0View file
@@ -0,0 +1,5 @@
1+import { createRoot } from "react-dom/client";
2+import App from "./App";
3+import "./style.css";
4+
5+createRoot(document.getElementById("root")!).render(<App />);
src/sim/analysis.tsadded+141−0View file
@@ -0,0 +1,141 @@
1+// Signal analysis for display: FFT, spectrogram, spectra, mode-energy bands.
2+
3+export function fft(re: Float64Array, im: Float64Array) {
4+ const n = re.length;
5+ for (let i = 1, j = 0; i < n; i++) {
6+ let bit = n >> 1;
7+ for (; j & bit; bit >>= 1) j ^= bit;
8+ j ^= bit;
9+ if (i < j) {
10+ [re[i], re[j]] = [re[j], re[i]];
11+ [im[i], im[j]] = [im[j], im[i]];
12+ }
13+ }
14+ for (let len = 2; len <= n; len <<= 1) {
15+ const ang = (-2 * Math.PI) / len;
16+ const wr = Math.cos(ang);
17+ const wi = Math.sin(ang);
18+ for (let i = 0; i < n; i += len) {
19+ let cr = 1;
20+ let ci = 0;
21+ for (let k = 0; k < len / 2; k++) {
22+ const a = i + k;
23+ const b = a + len / 2;
24+ const tr = re[b] * cr - im[b] * ci;
25+ const ti = re[b] * ci + im[b] * cr;
26+ re[b] = re[a] - tr;
27+ im[b] = im[a] - ti;
28+ re[a] += tr;
29+ im[a] += ti;
30+ const nr = cr * wr - ci * wi;
31+ ci = cr * wi + ci * wr;
32+ cr = nr;
33+ }
34+ }
35+ }
36+}
37+
38+export interface Spectrogram {
39+ nFrames: number;
40+ nBins: number;
41+ hop: number;
42+ win: number;
43+ db: Float32Array; // frame-major, nFrames x nBins
44+ maxDb: number;
45+}
46+
47+export function spectrogram(x: Float32Array, win = 2048, hop = 256): Spectrogram {
48+ const nFrames = Math.max(1, Math.floor((x.length - win) / hop) + 1);
49+ const nBins = win / 2;
50+ const db = new Float32Array(nFrames * nBins);
51+ const re = new Float64Array(win);
52+ const im = new Float64Array(win);
53+ const w = new Float64Array(win);
54+ for (let i = 0; i < win; i++) w[i] = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / win);
55+ let maxDb = -Infinity;
56+ for (let f = 0; f < nFrames; f++) {
57+ for (let i = 0; i < win; i++) {
58+ re[i] = (x[f * hop + i] ?? 0) * w[i];
59+ im[i] = 0;
60+ }
61+ fft(re, im);
62+ for (let b = 0; b < nBins; b++) {
63+ const v = 10 * Math.log10(re[b] * re[b] + im[b] * im[b] + 1e-300);
64+ db[f * nBins + b] = v;
65+ if (v > maxDb) maxDb = v;
66+ }
67+ }
68+ return { nFrames, nBins, hop, win, db, maxDb };
69+}
70+
71+// Power spectrum (dB) of x[start, start+len) with a Hann window, zero padded to a power of two.
72+export function spectrum(x: Float32Array, start: number, len: number): { db: Float64Array; df: number; nfft: number } {
73+ len = Math.max(2, Math.min(len, x.length - start));
74+ let nfft = 1;
75+ while (nfft < 2 * len) nfft <<= 1;
76+ const re = new Float64Array(nfft);
77+ const im = new Float64Array(nfft);
78+ for (let i = 0; i < len; i++) re[i] = x[start + i] * (0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (len - 1)));
79+ fft(re, im);
80+ const db = new Float64Array(nfft / 2);
81+ for (let b = 0; b < nfft / 2; b++) db[b] = 10 * Math.log10(re[b] * re[b] + im[b] * im[b] + 1e-300);
82+ return { db, df: 1 / nfft, nfft };
83+}
84+
85+// Sum mode energies into log-spaced frequency bands.
86+export function modeBands(
87+ freqs: Float64Array,
88+ modeEnergy: Float32Array,
89+ nFrames: number,
90+ nBands = 48
91+): { edges: Float64Array; e: Float32Array; nBands: number } {
92+ const nModes = freqs.length;
93+ const fmin = Math.max(freqs[0], 1e-3) * 0.97;
94+ const fmax = freqs[nModes - 1] * 1.03;
95+ const edges = new Float64Array(nBands + 1);
96+ for (let b = 0; b <= nBands; b++) edges[b] = fmin * Math.pow(fmax / fmin, b / nBands);
97+ const band = new Int32Array(nModes);
98+ for (let k = 0; k < nModes; k++) {
99+ band[k] = Math.min(nBands - 1, Math.max(0, Math.floor((Math.log(freqs[k] / fmin) / Math.log(fmax / fmin)) * nBands)));
100+ }
101+ const e = new Float32Array(nFrames * nBands);
102+ for (let f = 0; f < nFrames; f++) {
103+ for (let k = 0; k < nModes; k++) e[f * nBands + band[k]] += modeEnergy[f * nModes + k];
104+ }
105+ return { edges, e, nBands };
106+}
107+
108+// Relative RMS difference between two signals over [0, end).
109+export function relDiff(a: Float32Array, b: Float32Array, end: number): number {
110+ let num = 0;
111+ let den = 0;
112+ for (let i = 0; i < Math.min(end, a.length, b.length); i++) {
113+ num += (a[i] - b[i]) ** 2;
114+ den += a[i] ** 2;
115+ }
116+ return den > 0 ? Math.sqrt(num / den) : 0;
117+}
118+
119+export function wavBlob(x: Float32Array, fs: number): Blob {
120+ let peak = 0;
121+ for (const v of x) peak = Math.max(peak, Math.abs(v));
122+ const g = peak > 0 ? 0.9 / peak : 1;
123+ const buf = new ArrayBuffer(44 + 2 * x.length);
124+ const dv = new DataView(buf);
125+ const w = (o: number, s: string) => [...s].forEach((c, i) => dv.setUint8(o + i, c.charCodeAt(0)));
126+ w(0, "RIFF");
127+ dv.setUint32(4, 36 + 2 * x.length, true);
128+ w(8, "WAVE");
129+ w(12, "fmt ");
130+ dv.setUint32(16, 16, true);
131+ dv.setUint16(20, 1, true);
132+ dv.setUint16(22, 1, true);
133+ dv.setUint32(24, fs, true);
134+ dv.setUint32(28, 2 * fs, true);
135+ dv.setUint16(32, 2, true);
136+ dv.setUint16(34, 16, true);
137+ w(36, "data");
138+ dv.setUint32(40, 2 * x.length, true);
139+ for (let i = 0; i < x.length; i++) dv.setInt16(44 + 2 * i, Math.round(32767 * Math.max(-1, Math.min(1, g * x[i]))), true);
140+ return new Blob([buf], { type: "audio/wav" });
141+}
src/sim/modes.tsadded+66−0View file
@@ -0,0 +1,66 @@
1+// Linear transverse modes of a flat bead system.
2+//
3+// When every bead rests in a plane z = const and every spring lies in that
4+// plane, small transverse (z) motion decouples from in-plane motion, and its
5+// stiffness matrix is a graph Laplacian weighted by spring tension over length:
6+// K = sum_s (f_s / L_s) (e_i - e_j)(e_i - e_j)^T. We return the eigenpairs of
7+// M^{-1/2} K M^{-1/2} over the free beads.
8+
9+import type { BeadSystem, Modes } from "./types";
10+
11+export type SymEig = (A: Float64Array, n: number) => { values: Float64Array; vectors: Float64Array };
12+
13+export function isPlanar(sys: BeadSystem): boolean {
14+ const { n, pos0 } = sys;
15+ let zmin = Infinity;
16+ let zmax = -Infinity;
17+ let span = 0;
18+ for (let i = 0; i < n; i++) {
19+ zmin = Math.min(zmin, pos0[3 * i + 2]);
20+ zmax = Math.max(zmax, pos0[3 * i + 2]);
21+ span = Math.max(span, Math.abs(pos0[3 * i]), Math.abs(pos0[3 * i + 1]));
22+ }
23+ return zmax - zmin <= 1e-12 * Math.max(span, 1e-30);
24+}
25+
26+export function transverseStiffness(sys: BeadSystem): { K: Float64Array; free: Int32Array } {
27+ const { n, pinned, si, sj, sk, sL0, sk3, pos0, mass } = sys;
28+ const freeIdx = new Int32Array(n).fill(-1);
29+ const free: number[] = [];
30+ for (let i = 0; i < n; i++) if (!pinned[i]) freeIdx[i] = free.push(i) - 1;
31+ const nf = free.length;
32+ const K = new Float64Array(nf * nf);
33+ for (let s = 0; s < si.length; s++) {
34+ const i3 = 3 * si[s];
35+ const j3 = 3 * sj[s];
36+ const L = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
37+ const e = L - sL0[s];
38+ const g = (sk[s] * e + sk3[s] * e * e * e) / L;
39+ const a = freeIdx[si[s]];
40+ const b = freeIdx[sj[s]];
41+ if (a >= 0) K[a * nf + a] += g;
42+ if (b >= 0) K[b * nf + b] += g;
43+ if (a >= 0 && b >= 0) {
44+ K[a * nf + b] -= g;
45+ K[b * nf + a] -= g;
46+ }
47+ }
48+ for (let a = 0; a < nf; a++) {
49+ for (let b = 0; b < nf; b++) K[a * nf + b] /= Math.sqrt(mass[free[a]] * mass[free[b]]);
50+ }
51+ return { K, free: Int32Array.from(free) };
52+}
53+
54+export function computeModes(sys: BeadSystem, eig: SymEig): Modes {
55+ const { K, free } = transverseStiffness(sys);
56+ const nf = free.length;
57+ const { values, vectors } = eig(K, nf);
58+ const order = Array.from(values.keys()).sort((p, q) => values[p] - values[q]);
59+ const freqs = new Float64Array(nf);
60+ const psi = new Float64Array(nf * nf);
61+ order.forEach((src, k) => {
62+ freqs[k] = Math.sqrt(Math.max(values[src], 0)) / (2 * Math.PI);
63+ psi.set(vectors.subarray(src * nf, (src + 1) * nf), k * nf);
64+ });
65+ return { freqs, psi, free };
66+}
src/sim/parse.tsadded+210−0View file
@@ -0,0 +1,210 @@
1+// Convert the `sys` struct produced by a numbl script into a BeadSystem.
2+// The script uses 1-based bead indices; they become 0-based here.
3+
4+import type { RuntimeValue } from "numbl";
5+import type { BeadSystem, Hit, Method } from "./types";
6+
7+type Mat = { rows: number; cols: number; data: Float64Array };
8+
9+class SysError extends Error {}
10+
11+function asMat(v: RuntimeValue | undefined, name: string): Mat | null {
12+ if (v === undefined) return null;
13+ if (typeof v === "number") return { rows: 1, cols: 1, data: Float64Array.of(v) };
14+ if (typeof v === "boolean") return { rows: 1, cols: 1, data: Float64Array.of(v ? 1 : 0) };
15+ if (typeof v === "object" && v !== null && (v as { kind?: string }).kind === "tensor") {
16+ const t = v as unknown as { data: Float64Array; shape: number[]; imag?: Float64Array };
17+ if (t.imag) throw new SysError(`${name} must be real`);
18+ const rows = t.shape[0] ?? 1;
19+ const cols = t.shape.slice(1).reduce((a, b) => a * b, 1);
20+ return { rows, cols, data: t.data };
21+ }
22+ throw new SysError(`${name} must be numeric`);
23+}
24+
25+function scalar(v: RuntimeValue | undefined, name: string, def: number): number {
26+ const m = asMat(v, name);
27+ if (!m || m.data.length === 0) return def;
28+ if (m.data.length !== 1) throw new SysError(`${name} must be a scalar`);
29+ return m.data[0];
30+}
31+
32+function vec3(v: RuntimeValue | undefined, name: string, def: [number, number, number]): [number, number, number] {
33+ const m = asMat(v, name);
34+ if (!m || m.data.length === 0) return def;
35+ if (m.data.length !== 3) throw new SysError(`${name} must have 3 elements`);
36+ return [m.data[0], m.data[1], m.data[2]];
37+}
38+
39+function str(v: RuntimeValue | undefined, name: string, def: string): string {
40+ if (v === undefined) return def;
41+ if (typeof v === "string") return v;
42+ if (typeof v === "object" && v !== null && (v as { kind?: string }).kind === "char") {
43+ return (v as unknown as { value: string }).value;
44+ }
45+ throw new SysError(`${name} must be a string`);
46+}
47+
48+function perBead(v: RuntimeValue | undefined, name: string, n: number, def: number): Float64Array {
49+ const out = new Float64Array(n);
50+ const m = asMat(v, name);
51+ if (!m || m.data.length === 0) return out.fill(def);
52+ if (m.data.length === 1) return out.fill(m.data[0]);
53+ if (m.data.length !== n) throw new SysError(`${name} must be a scalar or have one entry per bead (${n})`);
54+ out.set(m.data);
55+ return out;
56+}
57+
58+function indices(v: RuntimeValue | undefined, name: string, n: number): Int32Array {
59+ const m = asMat(v, name);
60+ if (!m) return new Int32Array(0);
61+ const isLogical = (v as { _isLogical?: boolean })._isLogical || typeof v === "boolean";
62+ if (isLogical) {
63+ if (m.data.length !== n) throw new SysError(`${name}: a logical mask must have one entry per bead (${n})`);
64+ const out: number[] = [];
65+ m.data.forEach((b, i) => b && out.push(i));
66+ return Int32Array.from(out);
67+ }
68+ const out = new Int32Array(m.data.length);
69+ m.data.forEach((b, q) => {
70+ if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`${name}: bead index ${b} is not in 1..${n}`);
71+ out[q] = b - 1;
72+ });
73+ return out;
74+}
75+
76+function fields(v: RuntimeValue | undefined, name: string): Map<string, RuntimeValue>[] {
77+ if (v === undefined) return [];
78+ const k = (v as { kind?: string }).kind;
79+ if (k === "struct") return [(v as unknown as { fields: Map<string, RuntimeValue> }).fields];
80+ if (k === "struct_array") {
81+ return (v as unknown as { elements: { fields: Map<string, RuntimeValue> }[] }).elements.map(e => e.fields);
82+ }
83+ throw new SysError(`${name} must be a struct or struct array`);
84+}
85+
86+const KNOWN = new Set([
87+ "pos", "mass", "pinned", "springs", "dashpots", "friction", "hits", "listen",
88+ "duration", "sampleRate", "dt", "safety", "method", "frameRate", "modes",
89+]);
90+
91+export function parseSystem(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warnings: string[] } {
92+ try {
93+ return parseInner(sysValue);
94+ } catch (e) {
95+ if (e instanceof SysError) throw new Error(e.message);
96+ throw e;
97+ }
98+}
99+
100+function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warnings: string[] } {
101+ if (sysValue === undefined) throw new SysError("The script must define a struct named sys.");
102+ const [f] = fields(sysValue, "sys");
103+ if (!f) throw new SysError("sys must be a single struct");
104+ const warnings: string[] = [];
105+ for (const key of f.keys()) if (!KNOWN.has(key)) warnings.push(`sys.${key} is not used by the solver.`);
106+
107+ const P = asMat(f.get("pos"), "sys.pos");
108+ if (!P || P.cols !== 3 || P.rows < 1) throw new SysError("sys.pos must be an n-by-3 matrix of bead positions");
109+ const n = P.rows;
110+ const pos0 = new Float64Array(3 * n);
111+ for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) pos0[3 * i + c] = P.data[i + c * n];
112+
113+ if (!f.has("mass")) throw new SysError("sys.mass is required (scalar or one per bead)");
114+ const mass = perBead(f.get("mass"), "sys.mass", n, 0);
115+ for (let i = 0; i < n; i++) if (!(mass[i] > 0)) throw new SysError(`sys.mass must be positive (bead ${i + 1})`);
116+
117+ const pinned = new Uint8Array(n);
118+ for (const b of indices(f.get("pinned"), "sys.pinned", n)) pinned[b] = 1;
119+
120+ const S = asMat(f.get("springs"), "sys.springs");
121+ const ns = S && S.data.length ? S.rows : 0;
122+ if (S && ns && S.cols !== 4 && S.cols !== 5) {
123+ throw new SysError("sys.springs must have 4 or 5 columns: [i j k L0] or [i j k L0 k3]");
124+ }
125+ const si = new Int32Array(ns);
126+ const sj = new Int32Array(ns);
127+ const sk = new Float64Array(ns);
128+ const sL0 = new Float64Array(ns);
129+ const sk3 = new Float64Array(ns);
130+ for (let s = 0; s < ns; s++) {
131+ const col = (c: number) => S!.data[s + c * ns];
132+ const i = col(0);
133+ const j = col(1);
134+ for (const b of [i, j]) {
135+ if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.springs row ${s + 1}: bead index ${b} is not in 1..${n}`);
136+ }
137+ if (i === j) throw new SysError(`sys.springs row ${s + 1} joins bead ${i} to itself`);
138+ si[s] = i - 1;
139+ sj[s] = j - 1;
140+ sk[s] = col(2);
141+ sL0[s] = col(3);
142+ sk3[s] = S!.cols === 5 ? col(4) : 0;
143+ if (sk[s] < 0 || sk3[s] < 0 || sL0[s] < 0) throw new SysError(`sys.springs row ${s + 1}: k, L0 and k3 must be nonnegative`);
144+ }
145+
146+ const D = asMat(f.get("dashpots"), "sys.dashpots");
147+ const nd = D && D.data.length ? D.rows : 0;
148+ if (D && nd && D.cols !== 3) throw new SysError("sys.dashpots must have 3 columns: [i j c]");
149+ const di = new Int32Array(nd);
150+ const dj = new Int32Array(nd);
151+ const dc = new Float64Array(nd);
152+ for (let d = 0; d < nd; d++) {
153+ const col = (c: number) => D!.data[d + c * nd];
154+ for (const b of [col(0), col(1)]) {
155+ if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.dashpots row ${d + 1}: bead index ${b} is not in 1..${n}`);
156+ }
157+ di[d] = col(0) - 1;
158+ dj[d] = col(1) - 1;
159+ dc[d] = col(2);
160+ if (dc[d] < 0) throw new SysError(`sys.dashpots row ${d + 1}: c must be nonnegative`);
161+ }
162+ const friction = perBead(f.get("friction"), "sys.friction", n, 0);
163+
164+ const hits: Hit[] = fields(f.get("hits"), "sys.hits").map((h, q) => {
165+ const nm = `sys.hits(${q + 1})`;
166+ const beads = indices(h.get("beads"), `${nm}.beads`, n);
167+ if (beads.length === 0) throw new SysError(`${nm}.beads is required`);
168+ const duration = scalar(h.get("duration"), `${nm}.duration`, 0);
169+ const hit: Hit = {
170+ t: scalar(h.get("t"), `${nm}.t`, 0),
171+ beads,
172+ velocity: vec3(h.get("velocity"), `${nm}.velocity`, [0, 0, 0]),
173+ force: vec3(h.get("force"), `${nm}.force`, [0, 0, 0]),
174+ duration,
175+ };
176+ // In a struct array every hit has every field, so treat empty [] as absent.
177+ const given = (k: string) => (asMat(h.get(k), `${nm}.${k}`)?.data.length ?? 0) > 0;
178+ if (duration > 0 && given("velocity")) throw new SysError(`${nm}: give velocity (impulse) or force with duration (pulse), not both`);
179+ if (duration === 0 && given("force")) throw new SysError(`${nm}.force needs a positive duration`);
180+ return hit;
181+ });
182+ if (hits.length === 0) warnings.push("sys.hits is empty: nothing excites the system.");
183+
184+ const [L] = fields(f.get("listen"), "sys.listen");
185+ if (!L) throw new SysError("sys.listen is required, e.g. struct('beads', 10, 'dir', [0 0 1])");
186+ const listenBeads = indices(L.get("beads") ?? L.get("bead"), "sys.listen.beads", n);
187+ if (listenBeads.length === 0) throw new SysError("sys.listen.beads is required");
188+ const dir = vec3(L.get("dir"), "sys.listen.dir", [0, 0, 1]);
189+ const dn = Math.hypot(...dir);
190+ if (dn === 0) throw new SysError("sys.listen.dir must be nonzero");
191+
192+ const method = str(f.get("method"), "sys.method", "verlet") as Method;
193+ if (method !== "verlet" && method !== "yoshida4") throw new SysError("sys.method must be 'verlet' or 'yoshida4'");
194+
195+ const sys: BeadSystem = {
196+ n, pos0, mass, pinned, si, sj, sk, sL0, sk3, di, dj, dc, friction, hits, listenBeads,
197+ listenDir: [dir[0] / dn, dir[1] / dn, dir[2] / dn],
198+ duration: scalar(f.get("duration"), "sys.duration", 2),
199+ sampleRate: scalar(f.get("sampleRate"), "sys.sampleRate", 48000),
200+ dt: scalar(f.get("dt"), "sys.dt", 0),
201+ safety: scalar(f.get("safety"), "sys.safety", 0.5),
202+ method,
203+ frameRate: scalar(f.get("frameRate"), "sys.frameRate", 500),
204+ modes: scalar(f.get("modes"), "sys.modes", 1) !== 0,
205+ };
206+ if (!(sys.duration > 0)) throw new SysError("sys.duration must be positive");
207+ if (!(sys.safety > 0 && sys.safety <= 1)) throw new SysError("sys.safety must be in (0, 1]");
208+ if (ns === 0) warnings.push("sys.springs is empty: the beads do not interact.");
209+ return { sys, warnings };
210+}
src/sim/runScript.tsadded+32−0View file
@@ -0,0 +1,32 @@
1+// Run a numbl script with the bp_* helpers on its path, and turn its `sys`
2+// struct into a BeadSystem. Shared by the browser worker and the tests.
3+
4+import { executeCode, RTV } from "numbl";
5+import { parseSystem } from "./parse";
6+import type { SymEig } from "./modes";
7+import type { BeadSystem } from "./types";
8+
9+export interface HelperFile {
10+ name: string;
11+ source: string;
12+}
13+
14+export function runScript(
15+ source: string,
16+ helpers: HelperFile[],
17+ onOutput: (text: string) => void
18+): { sys: BeadSystem; warnings: string[] } {
19+ const result = executeCode(source, { onOutput, displayResults: true }, helpers, "script.m");
20+ return parseSystem(result.variableValues.sys);
21+}
22+
23+// Symmetric eigendecomposition through numbl's eig.
24+export const numblEig: SymEig = (A, n) => {
25+ const r = executeCode("[V, D] = eig(A); d = diag(D);", {
26+ initialVariableValues: { A: RTV.tensor(A, [n, n]) },
27+ });
28+ const V = r.variableValues.V as unknown as { data: Float64Array };
29+ const d = r.variableValues.d as unknown as { data: Float64Array } | number;
30+ const values = typeof d === "number" ? Float64Array.of(d) : d.data;
31+ return { values, vectors: V.data ?? Float64Array.of(1) };
32+};
src/sim/solver.tsadded+533−0View file
@@ -0,0 +1,533 @@
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.
11+
12+import type { BeadSystem, Method, Modes, SimResult } from "./types";
13+
14+const YOSHIDA_W1 = 1 / (2 - Math.cbrt(2));
15+const YOSHIDA_W0 = -Math.cbrt(2) * YOSHIDA_W1;
16+
17+export function stepCoefficients(method: Method): number[] {
18+ return method === "yoshida4" ? [YOSHIDA_W1, YOSHIDA_W0, YOSHIDA_W1] : [1];
19+}
20+
21+// Largest h*omega for which the method is stable on x'' = -omega^2 x.
22+export 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+}
45+
46+// Error of the numerical oscillation frequency for x'' = -omega^2 x, in cents
47+// (positive: the simulated tone is sharp).
48+export 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+}
66+
67+export interface SimOptions {
68+ modes?: Modes | null;
69+ onProgress?: (t: number) => void;
70+ maxAnimBytes?: number;
71+}
72+
73+export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
74+ const t0 = performance.now();
75+ const warnings: string[] = [];
76+ const { n, mass, pinned } = sys;
77+ const N3 = 3 * n;
78+ const x = new Float64Array(sys.pos0);
79+ const v = new Float64Array(N3);
80+ const F = new Float64Array(N3);
81+ const invm = new Float64Array(n);
82+ for (let i = 0; i < n; i++) invm[i] = pinned[i] ? 0 : 1 / mass[i];
83+
84+ // ---- forces ------------------------------------------------------------
85+ const { si, sj, sk, sL0, sk3 } = sys;
86+ const ns = si.length;
87+ const internalForces = (): number => {
88+ F.fill(0);
89+ let pe = 0;
90+ for (let s = 0; s < ns; s++) {
91+ const i3 = 3 * si[s];
92+ const j3 = 3 * sj[s];
93+ const dx = x[j3] - x[i3];
94+ const dy = x[j3 + 1] - x[i3 + 1];
95+ const dz = x[j3 + 2] - x[i3 + 2];
96+ const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
97+ const e = L - sL0[s];
98+ const k = sk[s];
99+ const k3 = sk3[s];
100+ const e2 = e * e;
101+ pe += 0.5 * k * e2 + 0.25 * k3 * e2 * e2;
102+ const g = (k * e + k3 * e2 * e) / L;
103+ F[i3] += g * dx;
104+ F[i3 + 1] += g * dy;
105+ F[i3 + 2] += g * dz;
106+ F[j3] -= g * dx;
107+ F[j3 + 1] -= g * dy;
108+ F[j3 + 2] -= g * dz;
109+ }
110+ return pe;
111+ };
112+ const pulses = sys.hits.filter(h => h.duration > 0);
113+ const externalForces = (t: number) => {
114+ for (const p of pulses) {
115+ const u = (t - p.t) / p.duration;
116+ if (u <= 0 || u >= 1) continue;
117+ const w = 0.5 * (1 - Math.cos(2 * Math.PI * u));
118+ for (const b of p.beads) {
119+ F[3 * b] += w * p.force[0];
120+ F[3 * b + 1] += w * p.force[1];
121+ F[3 * b + 2] += w * p.force[2];
122+ }
123+ }
124+ };
125+
126+ // ---- rest state checks ---------------------------------------------------
127+ let pe = internalForces();
128+ const peRest = pe; // energy stored by pre-tension; energies are reported relative to rest
129+ let residualForce = 0;
130+ for (let i = 0; i < n; i++) {
131+ if (pinned[i]) continue;
132+ const f = Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]);
133+ if (f > residualForce) residualForce = f;
134+ }
135+ let maxTension = 0;
136+ for (let s = 0; s < ns; s++) {
137+ const i3 = 3 * si[s];
138+ const j3 = 3 * sj[s];
139+ const L = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
140+ const e = L - sL0[s];
141+ maxTension = Math.max(maxTension, Math.abs(sk[s] * e + sk3[s] * e * e * e));
142+ }
143+ if (residualForce > 1e-6 * Math.max(maxTension, 1e-30)) {
144+ warnings.push(
145+ `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
146+ `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
147+ );
148+ }
149+
150+ // ---- time step -----------------------------------------------------------
151+ const omegaMax = estimateOmegaMax(sys);
152+ const fs = sys.sampleRate;
153+ const hLimit = stabilityLimit(sys.method) / omegaMax;
154+ let substeps: number;
155+ if (sys.dt > 0) {
156+ substeps = Math.max(1, Math.round(1 / (fs * sys.dt)));
157+ if (1 / (fs * substeps) > hLimit) {
158+ warnings.push(
159+ `dt = ${(1 / (fs * substeps)).toExponential(3)} s exceeds the stability limit ${hLimit.toExponential(3)} s.`
160+ );
161+ }
162+ } else {
163+ substeps = Math.max(1, Math.ceil(1 / (fs * sys.safety * hLimit)));
164+ }
165+ const h = 1 / (fs * substeps);
166+ const nSamples = Math.max(1, Math.round(sys.duration * fs));
167+ const nSteps = nSamples * substeps;
168+ const cs = stepCoefficients(sys.method);
169+
170+ // ---- damping: v <- exp(-tau A) v ------------------------------------------
171+ const { di, dj, dc, friction } = sys;
172+ const nd = di.length;
173+ let normA = 0; // infinity-norm bound of A
174+ {
175+ const row = new Float64Array(n);
176+ for (let i = 0; i < n; i++) row[i] = friction[i];
177+ for (let d = 0; d < nd; d++) {
178+ row[di[d]] += 2 * dc[d];
179+ row[dj[d]] += 2 * dc[d];
180+ }
181+ for (let i = 0; i < n; i++) normA = Math.max(normA, row[i] * invm[i]);
182+ }
183+ const hasDamping = normA > 0;
184+ const frictionOnly = nd === 0;
185+ const term = new Float64Array(N3);
186+ const tmp = new Float64Array(N3);
187+ const fa = new Float64Array(n);
188+ for (let i = 0; i < n; i++) fa[i] = friction[i] * invm[i];
189+ const dai = new Float64Array(nd);
190+ const daj = new Float64Array(nd);
191+ for (let d = 0; d < nd; d++) {
192+ dai[d] = dc[d] * invm[di[d]];
193+ daj[d] = dc[d] * invm[dj[d]];
194+ }
195+ const applyA = (w: Float64Array, out: Float64Array) => {
196+ for (let i = 0; i < n; i++) {
197+ const a = fa[i];
198+ out[3 * i] = a * w[3 * i];
199+ out[3 * i + 1] = a * w[3 * i + 1];
200+ out[3 * i + 2] = a * w[3 * i + 2];
201+ }
202+ for (let d = 0; d < nd; d++) {
203+ const i3 = 3 * di[d];
204+ const j3 = 3 * dj[d];
205+ const ai = dai[d];
206+ const aj = daj[d];
207+ const r0 = w[i3] - w[j3];
208+ const r1 = w[i3 + 1] - w[j3 + 1];
209+ const r2 = w[i3 + 2] - w[j3 + 2];
210+ out[i3] += ai * r0;
211+ out[i3 + 1] += ai * r1;
212+ out[i3 + 2] += ai * r2;
213+ out[j3] -= aj * r0;
214+ out[j3 + 1] -= aj * r1;
215+ out[j3 + 2] -= aj * r2;
216+ }
217+ };
218+ const kinetic = () => {
219+ let ke = 0;
220+ for (let i = 0; i < n; i++) {
221+ ke += 0.5 * mass[i] * (v[3 * i] ** 2 + v[3 * i + 1] ** 2 + v[3 * i + 2] ** 2);
222+ }
223+ return ke;
224+ };
225+ let dissipated = 0;
226+ const applyDamping = (tau: number) => {
227+ if (!hasDamping || tau === 0) return;
228+ const keBefore = kinetic();
229+ if (frictionOnly) {
230+ for (let i = 0; i < n; i++) {
231+ const f = Math.exp(-tau * friction[i] * invm[i]);
232+ v[3 * i] *= f;
233+ v[3 * i + 1] *= f;
234+ v[3 * i + 2] *= f;
235+ }
236+ } else {
237+ // Taylor series with scaling; |tau_s| * normA <= 0.5 so terms shrink fast.
238+ const nsub = Math.max(1, Math.ceil((Math.abs(tau) * normA) / 0.5));
239+ const ts = tau / nsub;
240+ for (let r = 0; r < nsub; r++) {
241+ let vmax = 0;
242+ for (let q = 0; q < N3; q++) {
243+ term[q] = v[q];
244+ vmax = Math.max(vmax, Math.abs(v[q]));
245+ }
246+ if (vmax === 0) break;
247+ for (let k = 1; k < 40; k++) {
248+ applyA(term, tmp);
249+ const c = -ts / k;
250+ let tmax = 0;
251+ for (let q = 0; q < N3; q++) {
252+ const val = c * tmp[q];
253+ term[q] = val;
254+ v[q] += val;
255+ const a = Math.abs(val);
256+ if (a > tmax) tmax = a;
257+ }
258+ if (tmax <= 1e-17 * vmax) break;
259+ }
260+ }
261+ }
262+ dissipated += keBefore - kinetic();
263+ };
264+
265+ // ---- recording -------------------------------------------------------------
266+ const audio = new Float32Array(nSamples);
267+ const lb = sys.listenBeads;
268+ const [ldx, ldy, ldz] = sys.listenDir;
269+ const pos0 = sys.pos0;
270+ const listen = () => {
271+ let s = 0;
272+ for (const b of lb) {
273+ s +=
274+ (x[3 * b] - pos0[3 * b]) * ldx +
275+ (x[3 * b + 1] - pos0[3 * b + 1]) * ldy +
276+ (x[3 * b + 2] - pos0[3 * b + 2]) * ldz;
277+ }
278+ return s / lb.length;
279+ };
280+ const frameStride = Math.max(1, Math.round(fs / sys.frameRate));
281+ const nFrames = Math.floor((nSamples - 1) / frameStride) + 1;
282+ const maxAnimBytes = opts.maxAnimBytes ?? 150e6;
283+ const animStride = frameStride * Math.max(1, Math.ceil((nFrames * N3 * 4) / maxAnimBytes));
284+ const frameTimes = new Float64Array(nFrames);
285+ const energy = new Float64Array(nFrames);
286+ const dissipatedRec = new Float64Array(nFrames);
287+ const animTimes: number[] = [];
288+ const animFrames: Float32Array[] = [];
289+
290+ const modes = opts.modes ?? null;
291+ const nModes = modes ? modes.freqs.length : 0;
292+ const modeEnergy = modes ? new Float32Array(nFrames * nModes) : null;
293+ let bz: Float64Array | null = null;
294+ let bv: Float64Array | null = null;
295+ let om2: Float64Array | null = null;
296+ if (modes) {
297+ bz = new Float64Array(modes.free.length);
298+ bv = new Float64Array(modes.free.length);
299+ om2 = new Float64Array(nModes);
300+ for (let k = 0; k < nModes; k++) om2[k] = (2 * Math.PI * modes.freqs[k]) ** 2;
301+ }
302+ const recordModes = (f: number) => {
303+ if (!modes || !modeEnergy || !bz || !bv || !om2) return;
304+ const free = modes.free;
305+ const nf = free.length;
306+ for (let r = 0; r < nf; r++) {
307+ const b = free[r];
308+ const sm = Math.sqrt(mass[b]);
309+ bz[r] = sm * (x[3 * b + 2] - pos0[3 * b + 2]);
310+ bv[r] = sm * v[3 * b + 2];
311+ }
312+ const psi = modes.psi;
313+ for (let k = 0; k < nModes; k++) {
314+ let q = 0;
315+ let qd = 0;
316+ const off = k * nf;
317+ for (let r = 0; r < nf; r++) {
318+ q += psi[off + r] * bz[r];
319+ qd += psi[off + r] * bv[r];
320+ }
321+ modeEnergy[f * nModes + k] = 0.5 * (qd * qd + om2[k] * q * q);
322+ }
323+ };
324+
325+ // ---- impulses, by audio sample ----------------------------------------------
326+ const impulsesAt = new Map<number, typeof sys.hits>();
327+ for (const hit of sys.hits) {
328+ if (hit.duration > 0) continue;
329+ const s = Math.round(hit.t * fs);
330+ if (s >= nSamples) continue;
331+ if (!impulsesAt.has(s)) impulsesAt.set(s, []);
332+ impulsesAt.get(s)!.push(hit);
333+ }
334+ const tExcitationEnd = Math.max(0, ...sys.hits.map(hh => hh.t + hh.duration));
335+
336+ // ---- main loop ---------------------------------------------------------------
337+ // Each audio sample interval Delta = substeps * h is one macro step:
338+ // D(Delta/2) [substeps conservative steps] D(Delta/2). The damping halves of
339+ // neighboring macro steps are merged unless something needs the velocities
340+ // at the boundary (a frame record or an impulse).
341+ const Delta = 1 / fs;
342+ pe = internalForces();
343+ externalForces(0);
344+ let pendingD = 0;
345+ let aborted: string | null = null;
346+ let frame = 0;
347+ let lastProgress = performance.now();
348+ for (let s = 0; s <= nSteps; s++) {
349+ const t = s * h;
350+ const sample = s % substeps === 0 ? s / substeps : -1;
351+ const isFrame = sample >= 0 && sample % frameStride === 0 && frame < nFrames;
352+ const imps = sample >= 0 ? impulsesAt.get(sample) : undefined;
353+ if ((isFrame || imps) && pendingD !== 0) {
354+ applyDamping(pendingD);
355+ pendingD = 0;
356+ }
357+ if (imps) {
358+ for (const hit of imps) {
359+ for (const b of hit.beads) {
360+ if (pinned[b]) continue;
361+ v[3 * b] += hit.velocity[0];
362+ v[3 * b + 1] += hit.velocity[1];
363+ v[3 * b + 2] += hit.velocity[2];
364+ }
365+ }
366+ }
367+ if (sample >= 0 && sample < nSamples) audio[sample] = listen();
368+ if (isFrame) {
369+ frameTimes[frame] = t;
370+ energy[frame] = kinetic() + pe - peRest;
371+ dissipatedRec[frame] = dissipated;
372+ recordModes(frame);
373+ if (sample % animStride === 0) {
374+ const d = new Float32Array(N3);
375+ for (let q = 0; q < N3; q++) d[q] = x[q] - pos0[q];
376+ animTimes.push(t);
377+ animFrames.push(d);
378+ }
379+ frame++;
380+ }
381+ if (s === nSteps) break;
382+ if (sample >= 0) {
383+ pendingD += Delta / 2;
384+ applyDamping(pendingD);
385+ pendingD = 0;
386+ }
387+
388+ let tj = t;
389+ for (const c of cs) {
390+ const kh = (c * h) / 2;
391+ for (let i = 0; i < n; i++) {
392+ const im = invm[i] * kh;
393+ v[3 * i] += im * F[3 * i];
394+ v[3 * i + 1] += im * F[3 * i + 1];
395+ v[3 * i + 2] += im * F[3 * i + 2];
396+ }
397+ const dh = c * h;
398+ for (let q = 0; q < N3; q++) x[q] += dh * v[q];
399+ tj += c * h;
400+ pe = internalForces();
401+ externalForces(tj);
402+ for (let i = 0; i < n; i++) {
403+ const im = invm[i] * kh;
404+ v[3 * i] += im * F[3 * i];
405+ v[3 * i + 1] += im * F[3 * i + 1];
406+ v[3 * i + 2] += im * F[3 * i + 2];
407+ }
408+ }
409+ if ((s + 1) % substeps === 0) pendingD = Delta / 2;
410+
411+ if ((s & 1023) === 0) {
412+ if (!Number.isFinite(pe)) {
413+ aborted = `The solution blew up at t = ${t.toFixed(4)} s. Reduce dt (or safety) or the strike strength.`;
414+ break;
415+ }
416+ const now = performance.now();
417+ if (opts.onProgress && now - lastProgress > 200) {
418+ opts.onProgress(t);
419+ lastProgress = now;
420+ }
421+ }
422+ }
423+
424+ // energy drift after the last excitation
425+ let energyDrift = 0;
426+ let ref = -1;
427+ for (let f = 0; f < frame; f++) {
428+ if (frameTimes[f] < tExcitationEnd) continue;
429+ const total = energy[f] + dissipatedRec[f];
430+ if (ref < 0) ref = total;
431+ else if (ref > 0) energyDrift = Math.max(energyDrift, Math.abs(total - ref) / ref);
432+ }
433+
434+ return {
435+ audio,
436+ sampleRate: fs,
437+ dt: h,
438+ substeps,
439+ omegaMax,
440+ frameTimes: frameTimes.slice(0, frame),
441+ energy: energy.slice(0, frame),
442+ dissipated: dissipatedRec.slice(0, frame),
443+ animTimes: Float64Array.from(animTimes),
444+ animFrames,
445+ modeEnergy: modeEnergy ? modeEnergy.slice(0, frame * nModes) : null,
446+ residualForce,
447+ maxTension,
448+ energyDrift,
449+ freqErrors: [1000, 2000, 5000, 10000, 20000]
450+ .filter(hz => hz < fs / 2)
451+ .map(hz => ({ hz, cents: frequencyErrorCents(sys.method, 2 * Math.PI * hz, h) })),
452+ elapsedMs: performance.now() - t0,
453+ aborted,
454+ warnings,
455+ };
456+}
457+
458+// Highest frequency (rad/s) of the system linearized about the rest positions,
459+// by power iteration on M^{-1/2} H M^{-1/2}, with a small margin.
460+export function estimateOmegaMax(sys: BeadSystem): number {
461+ const { n, si, sj, sk, sL0, sk3, pos0, mass, pinned } = sys;
462+ const ns = si.length;
463+ // per-spring 3x3 blocks B = keff d d^T + (f/L)(I - d d^T)
464+ const B = new Float64Array(ns * 6);
465+ for (let s = 0; s < ns; s++) {
466+ const i3 = 3 * si[s];
467+ const j3 = 3 * sj[s];
468+ const dx = pos0[j3] - pos0[i3];
469+ const dy = pos0[j3 + 1] - pos0[i3 + 1];
470+ const dz = pos0[j3 + 2] - pos0[i3 + 2];
471+ const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
472+ const e = L - sL0[s];
473+ const keff = sk[s] + 3 * sk3[s] * e * e;
474+ const g = (sk[s] * e + sk3[s] * e * e * e) / L;
475+ const u = [dx / L, dy / L, dz / L];
476+ const a = keff - g;
477+ B[6 * s] = a * u[0] * u[0] + g;
478+ B[6 * s + 1] = a * u[1] * u[1] + g;
479+ B[6 * s + 2] = a * u[2] * u[2] + g;
480+ B[6 * s + 3] = a * u[0] * u[1];
481+ B[6 * s + 4] = a * u[0] * u[2];
482+ B[6 * s + 5] = a * u[1] * u[2];
483+ }
484+ const ism = new Float64Array(n);
485+ for (let i = 0; i < n; i++) ism[i] = pinned[i] ? 0 : 1 / Math.sqrt(mass[i]);
486+ let w = new Float64Array(3 * n);
487+ let seed = 12345;
488+ for (let q = 0; q < 3 * n; q++) {
489+ seed = (seed * 1103515245 + 12345) & 0x7fffffff;
490+ w[q] = ism[Math.floor(q / 3)] * (seed / 0x7fffffff - 0.5);
491+ }
492+ const u = new Float64Array(3 * n);
493+ const y = new Float64Array(3 * n);
494+ let lambda = 0;
495+ for (let it = 0; it < 500; it++) {
496+ let norm = 0;
497+ for (let q = 0; q < 3 * n; q++) norm += w[q] * w[q];
498+ norm = Math.sqrt(norm);
499+ if (norm === 0) return 1;
500+ for (let q = 0; q < 3 * n; q++) u[q] = (w[q] / norm) * ism[Math.floor(q / 3)];
501+ y.fill(0);
502+ for (let s = 0; s < ns; s++) {
503+ const i3 = 3 * si[s];
504+ const j3 = 3 * sj[s];
505+ const r0 = u[i3] - u[j3];
506+ const r1 = u[i3 + 1] - u[j3 + 1];
507+ const r2 = u[i3 + 2] - u[j3 + 2];
508+ const b = 6 * s;
509+ const f0 = B[b] * r0 + B[b + 3] * r1 + B[b + 4] * r2;
510+ const f1 = B[b + 3] * r0 + B[b + 1] * r1 + B[b + 5] * r2;
511+ const f2 = B[b + 4] * r0 + B[b + 5] * r1 + B[b + 2] * r2;
512+ y[i3] += f0;
513+ y[i3 + 1] += f1;
514+ y[i3 + 2] += f2;
515+ y[j3] -= f0;
516+ y[j3 + 1] -= f1;
517+ y[j3 + 2] -= f2;
518+ }
519+ let rq = 0;
520+ const wNext = new Float64Array(3 * n);
521+ for (let q = 0; q < 3 * n; q++) {
522+ wNext[q] = y[q] * ism[Math.floor(q / 3)];
523+ rq += (w[q] / norm) * wNext[q];
524+ }
525+ w = wNext;
526+ if (it > 20 && Math.abs(rq - lambda) < 1e-8 * Math.abs(rq)) {
527+ lambda = rq;
528+ break;
529+ }
530+ lambda = rq;
531+ }
532+ return Math.sqrt(Math.max(lambda, 1e-300)) * 1.02;
533+}
src/sim/types.tsadded+77−0View file
@@ -0,0 +1,77 @@
1+// The bead system as the solver sees it. Indices are 0-based here; the
2+// script uses 1-based indices and parse.ts converts them.
3+
4+export type Method = "verlet" | "yoshida4";
5+
6+export interface Hit {
7+ t: number; // start time (s)
8+ beads: Int32Array;
9+ // Impulse: velocity added to each bead at time t (m/s). Used when duration == 0.
10+ velocity: [number, number, number];
11+ // Force pulse: raised-cosine force on each bead, peak `force` (N), lasting `duration` (s).
12+ force: [number, number, number];
13+ duration: number;
14+}
15+
16+export interface BeadSystem {
17+ n: number;
18+ pos0: Float64Array; // 3n, interleaved x,y,z (m)
19+ mass: Float64Array; // n (kg)
20+ pinned: Uint8Array; // n, 1 = held fixed
21+ // Springs: f = k e + k3 e^3, e = L - L0
22+ si: Int32Array;
23+ sj: Int32Array;
24+ sk: Float64Array;
25+ sL0: Float64Array;
26+ sk3: Float64Array;
27+ // Dashpots on relative velocity: F_i = c (v_j - v_i)
28+ di: Int32Array;
29+ dj: Int32Array;
30+ dc: Float64Array;
31+ friction: Float64Array; // n, F_i = -c0_i v_i (N s/m)
32+ hits: Hit[];
33+ listenBeads: Int32Array;
34+ listenDir: [number, number, number]; // unit vector
35+ duration: number; // s
36+ sampleRate: number; // Hz, audio output
37+ dt: number; // s; 0 = automatic
38+ safety: number; // fraction of the stability limit used when dt is automatic
39+ method: Method;
40+ frameRate: number; // Hz, snapshots for animation and mode energies
41+ modes: boolean; // compute linear transverse modes and track their energies
42+}
43+
44+export interface SimProgress {
45+ t: number;
46+ duration: number;
47+}
48+
49+export interface Modes {
50+ freqs: Float64Array; // Hz, ascending
51+ // psi: orthonormal eigenvectors of M^{-1/2} K M^{-1/2} over free beads,
52+ // column-major (nFree x nModes)
53+ psi: Float64Array;
54+ free: Int32Array; // bead index of each row of psi
55+}
56+
57+export interface SimResult {
58+ audio: Float32Array; // listener displacement at sampleRate (m)
59+ sampleRate: number;
60+ dt: number;
61+ substeps: number; // solver steps per audio sample
62+ omegaMax: number; // rad/s, highest linearized frequency
63+ frameTimes: Float64Array; // times of the energy records (s)
64+ energy: Float64Array; // mechanical energy above the rest state at each frame (J)
65+ animTimes: Float64Array; // times of the stored displacement snapshots (s)
66+ animFrames: Float32Array[]; // displacement (3n) at each snapshot (m)
67+ dissipated: Float64Array; // cumulative energy removed by damping (J)
68+ // mode energies (nModes per frame, frame-major), present when modes were computed
69+ modeEnergy: Float32Array | null;
70+ residualForce: number; // max |F| on free beads at rest (N)
71+ maxTension: number; // max spring force at rest (N)
72+ energyDrift: number; // max relative change in E + dissipated after the last excitation
73+ freqErrors: { hz: number; cents: number }[]; // integrator frequency error at a few frequencies
74+ elapsedMs: number;
75+ aborted: string | null; // reason if the run stopped early
76+ warnings: string[];
77+}
src/style.cssadded+101−0View file
@@ -0,0 +1,101 @@
1+:root {
2+ --bg: #fafaf9;
3+ --panel: #ffffff;
4+ --text: #1c1917;
5+ --muted: #78716c;
6+ --border: #e7e5e4;
7+ --grid: #f0eeec;
8+ --grid-strong: #c8c3be;
9+ --axis: #a8a29e;
10+ --series: #2563eb;
11+ --series2: #d97706;
12+ --accent: #16a34a;
13+ --hit: #dc2626;
14+ --warn-bg: #fef3c7;
15+ --warn-text: #78350f;
16+ --err-bg: #fee2e2;
17+ --err-text: #7f1d1d;
18+ color-scheme: light;
19+}
20+@media (prefers-color-scheme: dark) {
21+ :root:not([data-theme="light"]) {
22+ --bg: #171615;
23+ --panel: #1f1e1c;
24+ --text: #e7e5e4;
25+ --muted: #a8a29e;
26+ --border: #353230;
27+ --grid: #2a2826;
28+ --grid-strong: #57534e;
29+ --axis: #6b6560;
30+ --series: #60a5fa;
31+ --series2: #fbbf24;
32+ --accent: #4ade80;
33+ --hit: #f87171;
34+ --warn-bg: #3b2f12;
35+ --warn-text: #fcd34d;
36+ --err-bg: #3f1717;
37+ --err-text: #fca5a5;
38+ color-scheme: dark;
39+ }
40+}
41+* { box-sizing: border-box; }
42+html, body { margin: 0; }
43+body {
44+ background: var(--bg);
45+ color: var(--text);
46+ font: 14px/1.45 system-ui, -apple-system, "Segoe UI", sans-serif;
47+}
48+header { padding: 10px 16px; border-bottom: 1px solid var(--border); display: flex; gap: 12px; align-items: baseline; flex-wrap: wrap; }
49+h1 { font-size: 17px; margin: 0; }
50+.sub { color: var(--muted); }
51+h3 { font-size: 13px; font-weight: 600; margin: 18px 0 6px; color: var(--muted); }
52+h4 { margin: 10px 0 4px; font-size: 13px; }
53+.columns { display: grid; grid-template-columns: minmax(360px, 5fr) 7fr; gap: 16px; padding: 12px 16px 40px; }
54+@media (max-width: 900px) { .columns { grid-template-columns: 1fr; } }
55+.left, .right { min-width: 0; }
56+.controls { display: flex; gap: 8px; align-items: center; flex-wrap: wrap; margin-bottom: 8px; }
57+.grow { flex: 1; }
58+button, select {
59+ font: inherit; color: var(--text); background: var(--panel);
60+ border: 1px solid var(--border); border-radius: 6px; padding: 4px 10px; cursor: pointer;
61+}
62+button:hover:not(:disabled) { border-color: var(--axis); }
63+button:disabled { opacity: 0.5; cursor: default; }
64+button.primary { background: var(--series); color: white; border-color: var(--series); }
65+.editor { border: 1px solid var(--border); border-radius: 6px; overflow: hidden; background: var(--panel); }
66+.editor .cm-editor { height: 520px; font-size: 13px; }
67+.editor .cm-editor.cm-focused { outline: none; }
68+.editor .cm-gutters { background: var(--panel); color: var(--muted); border-right: 1px solid var(--border); }
69+.editor .cm-content { caret-color: var(--text); }
70+.editor .cm-activeLine, .editor .cm-activeLineGutter { background: transparent; }
71+pre { white-space: pre-wrap; word-break: break-word; font: 12px/1.4 ui-monospace, SFMono-Regular, Menlo, monospace; margin: 8px 0; }
72+.log { background: var(--panel); border: 1px solid var(--border); border-radius: 6px; padding: 8px; max-height: 160px; overflow: auto; }
73+.error { background: var(--err-bg); color: var(--err-text); border-radius: 6px; padding: 8px; margin: 8px 0; }
74+.warn { background: var(--warn-bg); color: var(--warn-text); border-radius: 6px; padding: 6px 8px; margin: 6px 0; font-size: 13px; }
75+.status { margin: 8px 0; color: var(--muted); }
76+.bar { height: 4px; background: var(--grid); border-radius: 2px; margin-top: 4px; overflow: hidden; }
77+.bar > div { height: 100%; background: var(--series); }
78+.report { margin-top: 10px; }
79+.report table { border-collapse: collapse; width: 100%; font-size: 13px; }
80+.report td { padding: 3px 6px; border-bottom: 1px solid var(--border); vertical-align: top; }
81+.report td:first-child { color: var(--muted); white-space: nowrap; width: 1%; }
82+.caption { margin: 10px 0 4px; font-size: 13px; color: var(--muted); }
83+.helpers { margin-top: 10px; }
84+.helpers summary { cursor: pointer; color: var(--muted); }
85+.plot { width: 100%; height: 220px; display: block; background: var(--panel); border: 1px solid var(--border); border-radius: 6px; }
86+.plot.short { height: 150px; }
87+.plot.tall { height: 380px; }
88+.grab { cursor: grab; touch-action: none; }
89+.legend { font-size: 12px; color: var(--muted); margin: 4px 0; display: flex; gap: 6px; align-items: center; flex-wrap: wrap; }
90+.legend input[type="range"] { width: 120px; }
91+.swatch { display: inline-block; width: 14px; height: 3px; border-radius: 2px; }
92+.swatch.s1 { background: var(--series); }
93+.swatch.s2 { background: var(--series2); }
94+.ring { display: inline-block; width: 10px; height: 10px; border-radius: 50%; border: 2px solid; }
95+.ring.hit { border-color: var(--hit); }
96+.ring.listen { border-color: var(--accent); }
97+.mono { font-family: ui-monospace, monospace; font-size: 12px; min-width: 64px; }
98+.note, .empty { color: var(--muted); padding: 12px; border: 1px dashed var(--border); border-radius: 6px; }
99+.empty { padding: 40px; text-align: center; }
100+.options { font-size: 13px; color: var(--muted); gap: 14px; }
101+.options label { display: flex; gap: 4px; align-items: center; cursor: pointer; }
src/ui/BeadView.tsxadded+180−0View file
@@ -0,0 +1,180 @@
1+import { useEffect, useMemo, useRef, useState } from "react";
2+import type { RunPayload } from "../worker";
3+
4+// Orthographic 3D view of the beads, displaced by a chosen snapshot.
5+export function BeadView({ p, onTime }: { p: RunPayload; onTime: (t: number | null) => void }) {
6+ const { pos0, n, si, sj, pinned, listenBeads, hitBeads } = p.sys;
7+ const { animFrames, animTimes } = p.result;
8+ const canvasRef = useRef<HTMLCanvasElement>(null);
9+ const [frame, setFrame] = useState(0);
10+ const [playing, setPlaying] = useState(false);
11+ const [speed, setSpeed] = useState(1 / 64);
12+ const [gainExp, setGainExp] = useState(0); // log2 of extra exaggeration
13+ const view = useRef({ az: -0.6, el: 0.55 });
14+ const [, redraw] = useState(0);
15+
16+ const geom = useMemo(() => {
17+ let cx = 0, cy = 0, cz = 0;
18+ for (let i = 0; i < n; i++) {
19+ cx += pos0[3 * i];
20+ cy += pos0[3 * i + 1];
21+ cz += pos0[3 * i + 2];
22+ }
23+ cx /= n; cy /= n; cz /= n;
24+ let ext = 0;
25+ for (let i = 0; i < n; i++) ext = Math.max(ext, Math.hypot(pos0[3 * i] - cx, pos0[3 * i + 1] - cy, pos0[3 * i + 2] - cz));
26+ let dmax = 0;
27+ for (const fr of animFrames) for (let q = 0; q < fr.length; q++) dmax = Math.max(dmax, Math.abs(fr[q]));
28+ return { c: [cx, cy, cz], ext: ext || 1, dmax: dmax || 1 };
29+ }, [pos0, n, animFrames]);
30+ const listen = useMemo(() => new Set(listenBeads), [listenBeads]);
31+ const hits = useMemo(() => new Set(hitBeads), [hitBeads]);
32+
33+ useEffect(() => {
34+ if (!playing || animTimes.length < 2) return;
35+ const dtFrame = animTimes[1] - animTimes[0];
36+ let last = performance.now();
37+ let acc = frame;
38+ let raf = 0;
39+ const tick = (now: number) => {
40+ acc += ((now - last) / 1000) * speed / dtFrame;
41+ last = now;
42+ if (acc >= animTimes.length) acc = 0;
43+ setFrame(Math.floor(acc));
44+ raf = requestAnimationFrame(tick);
45+ };
46+ raf = requestAnimationFrame(tick);
47+ return () => cancelAnimationFrame(raf);
48+ // eslint-disable-next-line react-hooks/exhaustive-deps
49+ }, [playing, speed, animTimes]);
50+
51+ useEffect(() => {
52+ onTime(animTimes.length ? animTimes[frame] : null);
53+ }, [frame, animTimes, onTime]);
54+
55+ useEffect(() => {
56+ const c = canvasRef.current;
57+ if (!c) return;
58+ const draw = () => {
59+ const dpr = window.devicePixelRatio || 1;
60+ const w = c.clientWidth;
61+ const h = c.clientHeight;
62+ c.width = Math.round(w * dpr);
63+ c.height = Math.round(h * dpr);
64+ const ctx = c.getContext("2d")!;
65+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
66+ ctx.clearRect(0, 0, w, h);
67+ const style = getComputedStyle(c);
68+ const css = (v: string) => style.getPropertyValue(v).trim();
69+ const { az, el } = view.current;
70+ const [ca, sa, ce, se] = [Math.cos(az), Math.sin(az), Math.cos(el), Math.sin(el)];
71+ const scale = (0.42 * Math.min(w, h)) / geom.ext;
72+ const exag = ((0.12 * geom.ext) / geom.dmax) * Math.pow(2, gainExp);
73+ const d = animFrames[Math.min(frame, animFrames.length - 1)];
74+ const sx = new Float32Array(n);
75+ const sy = new Float32Array(n);
76+ for (let i = 0; i < n; i++) {
77+ const X = pos0[3 * i] - geom.c[0] + (d ? exag * d[3 * i] : 0);
78+ const Y = pos0[3 * i + 1] - geom.c[1] + (d ? exag * d[3 * i + 1] : 0);
79+ const Z = pos0[3 * i + 2] - geom.c[2] + (d ? exag * d[3 * i + 2] : 0);
80+ const x1 = ca * X - sa * Y;
81+ const y1 = sa * X + ca * Y;
82+ sx[i] = w / 2 + scale * x1;
83+ sy[i] = h / 2 - scale * (ce * Z - se * y1);
84+ }
85+ if (si.length <= 30000) {
86+ ctx.strokeStyle = css("--grid-strong");
87+ ctx.lineWidth = 0.6;
88+ ctx.beginPath();
89+ for (let s = 0; s < si.length; s++) {
90+ ctx.moveTo(sx[si[s]], sy[si[s]]);
91+ ctx.lineTo(sx[sj[s]], sy[sj[s]]);
92+ }
93+ ctx.stroke();
94+ }
95+ const r = Math.max(1.2, Math.min(3, (scale * geom.ext) / Math.sqrt(n) / 3));
96+ for (let i = 0; i < n; i++) {
97+ let col: string;
98+ if (pinned[i]) col = css("--muted");
99+ else if (d) {
100+ const z = Math.max(-1, Math.min(1, d[3 * i + 2] / geom.dmax));
101+ col = z >= 0 ? `rgba(214, 96, 77, ${0.25 + 0.75 * z})` : `rgba(67, 147, 195, ${0.25 - 0.75 * z})`;
102+ } else col = css("--series");
103+ ctx.fillStyle = col;
104+ ctx.beginPath();
105+ ctx.arc(sx[i], sy[i], r, 0, 2 * Math.PI);
106+ ctx.fill();
107+ }
108+ const ring = (set: Set<number>, color: string) => {
109+ ctx.strokeStyle = color;
110+ ctx.lineWidth = 2;
111+ for (const i of set) {
112+ ctx.beginPath();
113+ ctx.arc(sx[i], sy[i], r + 3, 0, 2 * Math.PI);
114+ ctx.stroke();
115+ }
116+ };
117+ ring(hits, css("--hit"));
118+ ring(listen, css("--accent"));
119+ };
120+ draw();
121+ const ro = new ResizeObserver(draw);
122+ ro.observe(c);
123+ return () => ro.disconnect();
124+ });
125+
126+ const drag = useRef<{ x: number; y: number } | null>(null);
127+ const t = animTimes.length ? animTimes[Math.min(frame, animTimes.length - 1)] : 0;
128+ return (
129+ <div>
130+ <canvas
131+ ref={canvasRef}
132+ className="plot tall grab"
133+ onPointerDown={e => {
134+ drag.current = { x: e.clientX, y: e.clientY };
135+ (e.target as HTMLElement).setPointerCapture(e.pointerId);
136+ }}
137+ onPointerMove={e => {
138+ if (!drag.current) return;
139+ view.current.az += (e.clientX - drag.current.x) * 0.01;
140+ view.current.el = Math.max(-1.5, Math.min(1.5, view.current.el + (e.clientY - drag.current.y) * 0.01));
141+ drag.current = { x: e.clientX, y: e.clientY };
142+ redraw(v => v + 1);
143+ }}
144+ onPointerUp={() => (drag.current = null)}
145+ />
146+ <div className="controls">
147+ <button onClick={() => setPlaying(!playing)} disabled={animTimes.length < 2}>
148+ {playing ? "Pause" : "Play"}
149+ </button>
150+ <input
151+ type="range"
152+ className="grow"
153+ min={0}
154+ max={Math.max(0, animTimes.length - 1)}
155+ value={frame}
156+ onChange={e => setFrame(+e.target.value)}
157+ />
158+ <span className="mono">{(t * 1000).toFixed(1)} ms</span>
159+ <label>
160+ speed{" "}
161+ <select value={speed} onChange={e => setSpeed(+e.target.value)}>
162+ {[1 / 256, 1 / 64, 1 / 16, 1 / 4].map(s => (
163+ <option key={s} value={s}>
164+ 1/{1 / s}
165+ </option>
166+ ))}
167+ </select>
168+ </label>
169+ <label>
170+ exaggerate{" "}
171+ <input type="range" min={-3} max={3} step={0.5} value={gainExp} onChange={e => setGainExp(+e.target.value)} />
172+ </label>
173+ </div>
174+ <div className="legend">
175+ Drag to rotate. Displacement exaggerated about {(((0.12 * geom.ext) / geom.dmax) * Math.pow(2, gainExp)).toPrecision(2)}x.
176+ Red up, blue down; <span className="ring hit" /> struck beads, <span className="ring listen" /> listener.
177+ </div>
178+ </div>
179+ );
180+}
src/ui/Editor.tsxadded+42−0View file
@@ -0,0 +1,42 @@
1+import { StreamLanguage } from "@codemirror/language";
2+import { octave } from "@codemirror/legacy-modes/mode/octave";
3+import { EditorState } from "@codemirror/state";
4+import { keymap } from "@codemirror/view";
5+import { basicSetup, EditorView } from "codemirror";
6+import { useEffect, useRef } from "react";
7+
8+export function Editor({ value, onChange, onRun }: { value: string; onChange: (s: string) => void; onRun: () => void }) {
9+ const host = useRef<HTMLDivElement>(null);
10+ const view = useRef<EditorView | null>(null);
11+ const cb = useRef({ onChange, onRun });
12+ cb.current = { onChange, onRun };
13+
14+ useEffect(() => {
15+ const v = new EditorView({
16+ parent: host.current!,
17+ state: EditorState.create({
18+ doc: value,
19+ extensions: [
20+ keymap.of([{ key: "Mod-Enter", run: () => (cb.current.onRun(), true) }]),
21+ basicSetup,
22+ StreamLanguage.define(octave),
23+ EditorView.updateListener.of(u => {
24+ if (u.docChanged) cb.current.onChange(u.state.doc.toString());
25+ }),
26+ ],
27+ }),
28+ });
29+ view.current = v;
30+ return () => v.destroy();
31+ // eslint-disable-next-line react-hooks/exhaustive-deps
32+ }, []);
33+
34+ useEffect(() => {
35+ const v = view.current;
36+ if (v && v.state.doc.toString() !== value) {
37+ v.dispatch({ changes: { from: 0, to: v.state.doc.length, insert: value } });
38+ }
39+ }, [value]);
40+
41+ return <div ref={host} className="editor" />;
42+}
src/ui/Plots.tsxadded+266−0View file
@@ -0,0 +1,266 @@
1+import { useEffect, useMemo, useRef, useState } from "react";
2+import { modeBands, spectrogram, spectrum } from "../sim/analysis";
3+import type { RunPayload } from "../worker";
4+import { drawAxes, heat, niceTicks, setupCanvas } from "./plot";
5+
6+function useCanvas(draw: (c: HTMLCanvasElement) => void, deps: unknown[]) {
7+ const ref = useRef<HTMLCanvasElement>(null);
8+ useEffect(() => {
9+ const c = ref.current;
10+ if (!c) return;
11+ const run = () => draw(c);
12+ run();
13+ const ro = new ResizeObserver(run);
14+ ro.observe(c);
15+ const mq = window.matchMedia("(prefers-color-scheme: dark)");
16+ mq.addEventListener("change", run);
17+ return () => {
18+ ro.disconnect();
19+ mq.removeEventListener("change", run);
20+ };
21+ // eslint-disable-next-line react-hooks/exhaustive-deps
22+ }, deps);
23+ return ref;
24+}
25+
26+export function Waveform({ p, cursor }: { p: RunPayload; cursor: number | null }) {
27+ const { audio, sampleRate } = p.result;
28+ const ref = useCanvas(
29+ c => {
30+ const f = setupCanvas(c);
31+ let peak = 0;
32+ for (const v of audio) peak = Math.max(peak, Math.abs(v));
33+ peak = peak || 1;
34+ const T = audio.length / sampleRate;
35+ drawAxes(f, [0, T], [-peak * 1e3, peak * 1e3], { xlabel: "time (s)", ylabel: "disp (mm)" });
36+ const { ctx, x0, y0, pw, ph } = f;
37+ ctx.fillStyle = f.css("--series");
38+ const per = audio.length / pw;
39+ for (let px = 0; px < pw; px++) {
40+ let lo = Infinity;
41+ let hi = -Infinity;
42+ const a = Math.floor(px * per);
43+ const b = Math.max(a + 1, Math.floor((px + 1) * per));
44+ for (let i = a; i < b && i < audio.length; i++) {
45+ lo = Math.min(lo, audio[i]);
46+ hi = Math.max(hi, audio[i]);
47+ }
48+ const yhi = y0 + ph / 2 - (hi / peak) * (ph / 2);
49+ const ylo = y0 + ph / 2 - (lo / peak) * (ph / 2);
50+ ctx.fillRect(x0 + px, yhi, 1, Math.max(1, ylo - yhi));
51+ }
52+ if (cursor !== null) {
53+ ctx.fillStyle = f.css("--accent");
54+ ctx.fillRect(x0 + (cursor / T) * pw, y0, 1.5, ph);
55+ }
56+ },
57+ [audio, cursor]
58+ );
59+ return <canvas ref={ref} className="plot short" />;
60+}
61+
62+export function Spectrogram({ p, fmax }: { p: RunPayload; fmax: number }) {
63+ const { audio, sampleRate } = p.result;
64+ const sg = useMemo(() => spectrogram(audio), [audio]);
65+ const ref = useCanvas(
66+ c => {
67+ const f = setupCanvas(c);
68+ const T = audio.length / sampleRate;
69+ const { ctx, x0, y0, pw, ph } = f;
70+ const binHz = sampleRate / sg.win;
71+ const nb = Math.min(sg.nBins, Math.ceil(fmax / binHz));
72+ const img = ctx.createImageData(sg.nFrames, nb);
73+ const range = 80;
74+ for (let fr = 0; fr < sg.nFrames; fr++) {
75+ for (let b = 0; b < nb; b++) {
76+ const t = (sg.db[fr * sg.nBins + b] - (sg.maxDb - range)) / range;
77+ const [r, g, bl] = heat(t);
78+ const o = 4 * ((nb - 1 - b) * sg.nFrames + fr);
79+ img.data[o] = r;
80+ img.data[o + 1] = g;
81+ img.data[o + 2] = bl;
82+ img.data[o + 3] = 255;
83+ }
84+ }
85+ const off = new OffscreenCanvas(sg.nFrames, nb);
86+ off.getContext("2d")!.putImageData(img, 0, 0);
87+ ctx.imageSmoothingEnabled = true;
88+ const tEnd = ((sg.nFrames - 1) * sg.hop + sg.win) / sampleRate;
89+ ctx.drawImage(off, x0, y0, (pw * tEnd) / T, ph);
90+ drawAxes(f, [0, T], [0, (nb * binHz) / 1000], { grid: false, xlabel: "time (s)", ylabel: "kHz" });
91+ },
92+ [sg, fmax]
93+ );
94+ return <canvas ref={ref} className="plot" />;
95+}
96+
97+export function Energy({ p }: { p: RunPayload }) {
98+ const { frameTimes, energy, dissipated } = p.result;
99+ const ref = useCanvas(
100+ c => {
101+ const f = setupCanvas(c);
102+ const T = p.sys.duration;
103+ let emax = 0;
104+ for (let i = 0; i < energy.length; i++) emax = Math.max(emax, energy[i] + dissipated[i]);
105+ emax = emax || 1;
106+ const unit = emax < 1e-3 ? 1e6 : emax < 1 ? 1e3 : 1;
107+ const uname = unit === 1e6 ? "µJ" : unit === 1e3 ? "mJ" : "J";
108+ drawAxes(f, [0, T], [0, emax * unit * 1.05], { xlabel: "time (s)", ylabel: `energy (${uname})` });
109+ const { ctx, x0, y0, pw, ph } = f;
110+ const line = (ys: (i: number) => number, color: string) => {
111+ ctx.strokeStyle = color;
112+ ctx.lineWidth = 1.5;
113+ ctx.beginPath();
114+ for (let i = 0; i < frameTimes.length; i++) {
115+ const px = x0 + (frameTimes[i] / T) * pw;
116+ const py = y0 + ph - (ys(i) / (emax * 1.05)) * ph;
117+ if (i) ctx.lineTo(px, py);
118+ else ctx.moveTo(px, py);
119+ }
120+ ctx.stroke();
121+ };
122+ line(i => energy[i] + dissipated[i], f.css("--series2"));
123+ line(i => energy[i], f.css("--series"));
124+ },
125+ [p]
126+ );
127+ return (
128+ <>
129+ <canvas ref={ref} className="plot short" />
130+ <div className="legend">
131+ <span className="swatch s1" /> in the beads <span className="swatch s2" /> plus dissipated by damping
132+ </div>
133+ </>
134+ );
135+}
136+
137+export function ModeEnergy({ p }: { p: RunPayload }) {
138+ const { modeEnergy, frameTimes } = p.result;
139+ const [tmax, setTmax] = useState(p.sys.duration);
140+ const freqs = p.modeFreqs;
141+ const bands = useMemo(
142+ () => (freqs && modeEnergy ? modeBands(freqs, modeEnergy, frameTimes.length) : null),
143+ [freqs, modeEnergy, frameTimes]
144+ );
145+ const ref = useCanvas(
146+ c => {
147+ if (!bands) return;
148+ const f = setupCanvas(c);
149+ const { ctx, x0, y0, pw, ph } = f;
150+ let nF = frameTimes.length;
151+ while (nF > 2 && frameTimes[nF - 1] > tmax) nF--;
152+ const { nBands, e, edges } = bands;
153+ // log energy per band, relative to the total linear energy at each time
154+ const img = ctx.createImageData(nF, nBands);
155+ const range = 8; // decades
156+ for (let fr = 0; fr < nF; fr++) {
157+ let tot = 0;
158+ for (let b = 0; b < nBands; b++) tot += e[fr * nBands + b];
159+ for (let b = 0; b < nBands; b++) {
160+ const v = tot > 0 ? Math.log10(e[fr * nBands + b] / tot + 1e-300) : -range;
161+ const [r, g, bl] = heat(1 + v / range);
162+ const o = 4 * ((nBands - 1 - b) * nF + fr);
163+ img.data[o] = r;
164+ img.data[o + 1] = g;
165+ img.data[o + 2] = bl;
166+ img.data[o + 3] = e[fr * nBands + b] > 0 ? 255 : 0; // bands with no modes stay blank
167+ }
168+ }
169+ const off = new OffscreenCanvas(nF, nBands);
170+ off.getContext("2d")!.putImageData(img, 0, 0);
171+ ctx.imageSmoothingEnabled = false;
172+ ctx.drawImage(off, x0, y0, pw, ph);
173+ const lf0 = Math.log10(edges[0]);
174+ const lf1 = Math.log10(edges[nBands]);
175+ const decs: number[] = [];
176+ for (let d = Math.ceil(lf0); d <= lf1; d++) for (const m of [1, 2, 5]) if (Math.log10(m) + d >= lf0 && Math.log10(m) + d <= lf1) decs.push(Math.log10(m) + d);
177+ drawAxes(f, [0, frameTimes[nF - 1]], [lf0, lf1], {
178+ grid: false,
179+ xlabel: "time (s)",
180+ ylabel: "mode freq (Hz)",
181+ yticks: decs.length ? decs : niceTicks(lf0, lf1, 4),
182+ yfmt: v => String(Math.round(Math.pow(10, v))),
183+ });
184+ },
185+ [bands, tmax]
186+ );
187+ if (!bands) return <div className="note">Mode energies are shown when linear modes are computed (flat systems only).</div>;
188+ return (
189+ <>
190+ <canvas ref={ref} className="plot" />
191+ <div className="legend">
192+ Share of the linear (quadratic) energy in each band of transverse modes, 8 decades of color. Show the first{" "}
193+ <select value={tmax} onChange={e => setTmax(+e.target.value)}>
194+ {[0.03, 0.1, 0.3, 1]
195+ .filter(t => t < p.sys.duration)
196+ .concat([p.sys.duration])
197+ .map(t => (
198+ <option key={t} value={t}>
199+ {t} s
200+ </option>
201+ ))}
202+ </select>
203+ </div>
204+ </>
205+ );
206+}
207+
208+export function Spectrum({ p, fmax }: { p: RunPayload; fmax: number }) {
209+ const { audio, sampleRate } = p.result;
210+ const T = audio.length / sampleRate;
211+ const [early, setEarly] = useState(0);
212+ const [late, setLate] = useState(Math.max(0, T - 0.25));
213+ const len = Math.round(Math.min(0.25, T) * sampleRate);
214+ const sE = useMemo(() => spectrum(audio, Math.round(early * sampleRate), len), [audio, early, len, sampleRate]);
215+ const sL = useMemo(() => spectrum(audio, Math.round(late * sampleRate), len), [audio, late, len, sampleRate]);
216+ const ref = useCanvas(
217+ c => {
218+ const f = setupCanvas(c);
219+ const { ctx, x0, y0, pw, ph } = f;
220+ let top = -Infinity;
221+ for (const s of [sE, sL]) for (const v of s.db) top = Math.max(top, v);
222+ const yr: [number, number] = [top - 90, top + 3];
223+ drawAxes(f, [0, fmax / 1000], [yr[0] - top, yr[1] - top], { xlabel: "kHz", ylabel: "dB" });
224+ if (p.modeFreqs) {
225+ ctx.strokeStyle = f.css("--grid-strong");
226+ for (const fm of p.modeFreqs) {
227+ if (fm > fmax) break;
228+ const px = x0 + (fm / fmax) * pw;
229+ ctx.beginPath();
230+ ctx.moveTo(px, y0);
231+ ctx.lineTo(px, y0 + 6);
232+ ctx.stroke();
233+ }
234+ }
235+ const line = (s: typeof sE, color: string) => {
236+ ctx.strokeStyle = color;
237+ ctx.lineWidth = 1.2;
238+ ctx.beginPath();
239+ const nb = Math.min(s.db.length, Math.ceil(fmax / (s.df * sampleRate)));
240+ for (let b = 0; b < nb; b++) {
241+ const px = x0 + ((b * s.df * sampleRate) / fmax) * pw;
242+ const py = y0 + ph - ((s.db[b] - yr[0]) / (yr[1] - yr[0])) * ph;
243+ if (b) ctx.lineTo(px, Math.min(py, y0 + ph));
244+ else ctx.moveTo(px, Math.min(py, y0 + ph));
245+ }
246+ ctx.stroke();
247+ };
248+ line(sL, f.css("--series2"));
249+ line(sE, f.css("--series"));
250+ },
251+ [sE, sL, fmax, p]
252+ );
253+ const maxStart = Math.max(0, T - len / sampleRate);
254+ return (
255+ <>
256+ <canvas ref={ref} className="plot" />
257+ <div className="legend">
258+ <span className="swatch s1" /> window at{" "}
259+ <input type="range" min={0} max={maxStart} step={0.005} value={early} onChange={e => setEarly(+e.target.value)} /> {early.toFixed(3)} s
260+ <span className="swatch s2" /> window at{" "}
261+ <input type="range" min={0} max={maxStart} step={0.005} value={late} onChange={e => setLate(+e.target.value)} /> {late.toFixed(3)} s
262+ {p.modeFreqs && <span> · ticks at top: linear mode frequencies</span>}
263+ </div>
264+ </>
265+ );
266+}
src/ui/plot.tsadded+126−0View file
@@ -0,0 +1,126 @@
1+// Minimal canvas plotting helpers.
2+
3+export interface Frame {
4+ ctx: CanvasRenderingContext2D;
5+ w: number;
6+ h: number;
7+ // plot area
8+ x0: number;
9+ y0: number;
10+ pw: number;
11+ ph: number;
12+ css: (name: string) => string;
13+}
14+
15+export function setupCanvas(canvas: HTMLCanvasElement, margin = { l: 56, r: 12, t: 10, b: 32 }): Frame {
16+ const dpr = window.devicePixelRatio || 1;
17+ const w = canvas.clientWidth;
18+ const h = canvas.clientHeight;
19+ canvas.width = Math.round(w * dpr);
20+ canvas.height = Math.round(h * dpr);
21+ const ctx = canvas.getContext("2d")!;
22+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
23+ ctx.clearRect(0, 0, w, h);
24+ const style = getComputedStyle(canvas);
25+ return {
26+ ctx,
27+ w,
28+ h,
29+ x0: margin.l,
30+ y0: margin.t,
31+ pw: Math.max(10, w - margin.l - margin.r),
32+ ph: Math.max(10, h - margin.t - margin.b),
33+ css: name => style.getPropertyValue(name).trim(),
34+ };
35+}
36+
37+export function niceTicks(lo: number, hi: number, count = 5): number[] {
38+ if (!(hi > lo)) return [lo];
39+ const step0 = (hi - lo) / count;
40+ const mag = Math.pow(10, Math.floor(Math.log10(step0)));
41+ const step = [1, 2, 5, 10].map(m => m * mag).find(s => s >= step0) ?? 10 * mag;
42+ const out: number[] = [];
43+ for (let v = Math.ceil(lo / step) * step; v <= hi + 1e-9 * step; v += step) out.push(Math.abs(v) < 1e-12 * step ? 0 : v);
44+ return out;
45+}
46+
47+export function fmt(v: number): string {
48+ if (v === 0) return "0";
49+ const a = Math.abs(v);
50+ if (a >= 1e4 || a < 1e-3) return v.toExponential(0);
51+ return String(+v.toPrecision(3));
52+}
53+
54+export function drawAxes(
55+ f: Frame,
56+ xr: [number, number],
57+ yr: [number, number],
58+ opts: { grid?: boolean; xlabel?: string; ylabel?: string; xticks?: number[]; yticks?: number[]; yfmt?: (v: number) => string; xfmt?: (v: number) => string; ylog?: boolean } = {}
59+) {
60+ const { ctx, x0, y0, pw, ph } = f;
61+ ctx.font = "11px system-ui, sans-serif";
62+ ctx.strokeStyle = f.css("--grid");
63+ ctx.fillStyle = f.css("--muted");
64+ ctx.lineWidth = 1;
65+ const xt = opts.xticks ?? niceTicks(xr[0], xr[1]);
66+ const yt = opts.yticks ?? niceTicks(yr[0], yr[1], 4);
67+ ctx.textAlign = "center";
68+ ctx.textBaseline = "top";
69+ for (const v of xt) {
70+ const px = x0 + ((v - xr[0]) / (xr[1] - xr[0])) * pw;
71+ if (px < x0 - 0.5 || px > x0 + pw + 0.5) continue;
72+ if (opts.grid !== false) {
73+ ctx.beginPath();
74+ ctx.moveTo(px, y0);
75+ ctx.lineTo(px, y0 + ph);
76+ ctx.stroke();
77+ }
78+ ctx.fillText((opts.xfmt ?? fmt)(v), px, y0 + ph + 4);
79+ }
80+ ctx.textAlign = "right";
81+ ctx.textBaseline = "middle";
82+ for (const v of yt) {
83+ const py = y0 + ph - ((v - yr[0]) / (yr[1] - yr[0])) * ph;
84+ if (py < y0 - 0.5 || py > y0 + ph + 0.5) continue;
85+ if (opts.grid !== false) {
86+ ctx.beginPath();
87+ ctx.moveTo(x0, py);
88+ ctx.lineTo(x0 + pw, py);
89+ ctx.stroke();
90+ }
91+ ctx.fillText((opts.yfmt ?? fmt)(v), x0 - 5, py);
92+ }
93+ ctx.strokeStyle = f.css("--axis");
94+ ctx.strokeRect(x0 + 0.5, y0 + 0.5, pw - 1, ph - 1);
95+ if (opts.xlabel) {
96+ ctx.textAlign = "right";
97+ ctx.textBaseline = "bottom";
98+ ctx.fillText(opts.xlabel, x0 + pw, f.h - 1);
99+ }
100+ if (opts.ylabel) {
101+ ctx.save();
102+ ctx.translate(11, y0 + ph / 2);
103+ ctx.rotate(-Math.PI / 2);
104+ ctx.textAlign = "center";
105+ ctx.textBaseline = "middle";
106+ ctx.fillText(opts.ylabel, 0, 0);
107+ ctx.restore();
108+ }
109+}
110+
111+// Sequential color map (dark to bright), t in [0, 1].
112+export function heat(t: number): [number, number, number] {
113+ t = Math.max(0, Math.min(1, t));
114+ const stops: [number, number, number][] = [
115+ [12, 7, 35],
116+ [72, 18, 110],
117+ [150, 40, 110],
118+ [218, 85, 60],
119+ [250, 170, 50],
120+ [252, 250, 190],
121+ ];
122+ const s = t * (stops.length - 1);
123+ const i = Math.min(stops.length - 2, Math.floor(s));
124+ const u = s - i;
125+ return [0, 1, 2].map(c => Math.round(stops[i][c] + u * (stops[i + 1][c] - stops[i][c]))) as [number, number, number];
126+}
src/worker.tsadded+143−0View file
@@ -0,0 +1,143 @@
1+// Web Worker: runs the numbl script, computes linear modes, and simulates.
2+
3+import { relDiff } from "./sim/analysis";
4+import { computeModes, isPlanar } from "./sim/modes";
5+import { numblEig, runScript } from "./sim/runScript";
6+import { simulate } from "./sim/solver";
7+import type { BeadSystem, SimResult } from "./sim/types";
8+
9+const helperSources = import.meta.glob("../matlab/helpers/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
10+const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source }));
11+
12+const MAX_MODE_BEADS = 1500;
13+
14+export interface RunPayload {
15+ sys: {
16+ n: number;
17+ pos0: Float64Array;
18+ pinned: Uint8Array;
19+ si: Int32Array;
20+ sj: Int32Array;
21+ nDashpots: number;
22+ listenBeads: Int32Array;
23+ hitBeads: Int32Array;
24+ method: string;
25+ duration: number;
26+ };
27+ result: SimResult;
28+ modeFreqs: Float64Array | null;
29+ notes: string[];
30+ check: { dt: number; windows: { t: number; rel: number }[] } | null;
31+ // the same run with every hit scaled by 1e-3, its output scaled back up
32+ linear: { audio: Float32Array; windows: { t: number; rel: number }[] } | null;
33+}
34+
35+export interface RunRequest {
36+ type: "run";
37+ script: string;
38+ dtCheck: boolean;
39+ linear: boolean;
40+}
41+
42+export type FromWorker =
43+ | { type: "log"; text: string }
44+ | { type: "stage"; text: string }
45+ | { type: "progress"; t: number; duration: number; label: string }
46+ | { type: "result"; payload: RunPayload }
47+ | { type: "error"; message: string };
48+
49+const post = (m: FromWorker, transfer: Transferable[] = []) => (self as unknown as Worker).postMessage(m, transfer);
50+
51+const windowsOf = (duration: number) =>
52+ [0.01, 0.03, 0.1, 0.3, 1, duration].filter((t, i, a) => t <= duration && a.indexOf(t) === i);
53+
54+self.onmessage = (ev: MessageEvent<RunRequest>) => {
55+ const { script, dtCheck, linear } = ev.data;
56+ try {
57+ post({ type: "stage", text: "Running script" });
58+ const { sys, warnings } = runScript(script, helpers, text => post({ type: "log", text }));
59+ const notes = [...warnings];
60+
61+ let modes = null;
62+ const nFree = sys.pinned.reduce((a, p) => a + (p ? 0 : 1), 0);
63+ if (sys.modes) {
64+ if (!isPlanar(sys)) notes.push("Linear modes skipped: the beads do not rest in a plane z = const.");
65+ else if (nFree > MAX_MODE_BEADS) notes.push(`Linear modes skipped: ${nFree} free beads (limit ${MAX_MODE_BEADS}).`);
66+ else {
67+ post({ type: "stage", text: `Computing linear modes (${nFree} free beads)` });
68+ modes = computeModes(sys, numblEig);
69+ }
70+ }
71+
72+ post({ type: "stage", text: "Simulating" });
73+ const result = simulate(sys, {
74+ modes,
75+ onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Simulating" }),
76+ });
77+
78+ let checkOut: RunPayload["check"] = null;
79+ if (dtCheck && !result.aborted) {
80+ post({ type: "stage", text: "Simulating again with dt/2" });
81+ const half: BeadSystem = { ...sys, dt: result.dt / 2 };
82+ const r2 = simulate(half, {
83+ maxAnimBytes: 0,
84+ onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Simulating with dt/2" }),
85+ });
86+ const fs = result.sampleRate;
87+ const windows = windowsOf(sys.duration).map(t => ({ t, rel: relDiff(r2.audio, result.audio, Math.round(t * fs)) }));
88+ checkOut = { dt: r2.dt, windows };
89+ }
90+
91+ let linearOut: RunPayload["linear"] = null;
92+ if (linear && !result.aborted) {
93+ post({ type: "stage", text: "Simulating the small-strike reference" });
94+ const s = 1e-3;
95+ const small: BeadSystem = {
96+ ...sys,
97+ dt: result.dt,
98+ hits: sys.hits.map(h => ({
99+ ...h,
100+ velocity: [h.velocity[0] * s, h.velocity[1] * s, h.velocity[2] * s],
101+ force: [h.force[0] * s, h.force[1] * s, h.force[2] * s],
102+ })),
103+ };
104+ const r3 = simulate(small, {
105+ maxAnimBytes: 0,
106+ onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Small-strike reference" }),
107+ });
108+ const audio = r3.audio.map(v => v / s);
109+ const fs = result.sampleRate;
110+ const windows = windowsOf(sys.duration).map(t => ({ t, rel: relDiff(audio, result.audio, Math.round(t * fs)) }));
111+ linearOut = { audio, windows };
112+ }
113+
114+ const hitBeads = Int32Array.from(new Set(sys.hits.flatMap(h => Array.from(h.beads))));
115+ const payload: RunPayload = {
116+ sys: {
117+ n: sys.n,
118+ pos0: sys.pos0,
119+ pinned: sys.pinned,
120+ si: sys.si,
121+ sj: sys.sj,
122+ nDashpots: sys.di.length,
123+ listenBeads: sys.listenBeads,
124+ hitBeads,
125+ method: sys.method,
126+ duration: sys.duration,
127+ },
128+ result,
129+ modeFreqs: modes ? modes.freqs : null,
130+ notes,
131+ check: checkOut,
132+ linear: linearOut,
133+ };
134+ const transfer: Transferable[] = [result.audio.buffer, ...result.animFrames.map(f => f.buffer)];
135+ if (result.modeEnergy) transfer.push(result.modeEnergy.buffer);
136+ if (linearOut) transfer.push(linearOut.audio.buffer);
137+ post({ type: "result", payload }, transfer);
138+ } catch (e) {
139+ const err = e as { message?: string; line?: number; file?: string };
140+ const where = err.line ? ` (line ${err.line})` : "";
141+ post({ type: "error", message: `${err.message ?? String(e)}${where}` });
142+ }
143+};
tests/bench.tsadded+17−0View file
@@ -0,0 +1,17 @@
1+import { readFileSync, readdirSync } from "fs";
2+import { runScript, numblEig } from "../src/sim/runScript";
3+import { simulate } from "../src/sim/solver";
4+import { computeModes, isPlanar } from "../src/sim/modes";
5+const helpers = readdirSync("matlab/helpers").map(name => ({ name, source: readFileSync("matlab/helpers/" + name, "utf8") }));
6+const src = readFileSync(process.argv[2], "utf8");
7+let t = performance.now();
8+const { sys, warnings } = runScript(src, helpers, s => process.stdout.write(s));
9+console.log("script ms", (performance.now() - t).toFixed(0), warnings);
10+sys.duration = Number(process.argv[3] ?? 0.2);
11+t = performance.now();
12+const modes = isPlanar(sys) ? computeModes(sys, numblEig) : null;
13+console.log("modes ms", (performance.now() - t).toFixed(0), modes?.freqs.slice(0, 6));
14+const r = simulate(sys, { modes });
15+console.log({ dt: r.dt, substeps: r.substeps, omegaMax: r.omegaMax, ms: r.elapsedMs.toFixed(0), drift: r.energyDrift, resid: r.residualForce, maxT: r.maxTension, warnings: r.warnings, aborted: r.aborted });
16+const msPerSec = r.elapsedMs / sys.duration; console.log("x real time:", (msPerSec / 1000).toFixed(1));
17+let peak = 0; for (const a of r.audio) peak = Math.max(peak, Math.abs(a)); console.log("peak listener disp", peak);
tests/cascade.tsadded+23−0View file
@@ -0,0 +1,23 @@
1+// Share of linear (transverse) energy above fc over time: nonlinear strike vs the same strike scaled down.
2+// Usage: tsx tests/cascade.ts script.m duration fc scale...
3+import { readFileSync, readdirSync } from "fs";
4+import { numblEig, runScript } from "../src/sim/runScript";
5+import { simulate } from "../src/sim/solver";
6+import { computeModes } from "../src/sim/modes";
7+const helpers = readdirSync("matlab/helpers").map(name => ({ name, source: readFileSync("matlab/helpers/" + name, "utf8") }));
8+const { sys } = runScript(readFileSync(process.argv[2], "utf8"), helpers, () => {});
9+sys.duration = Number(process.argv[3]);
10+const fc = Number(process.argv[4]);
11+const modes = computeModes(sys, numblEig);
12+const nm = modes.freqs.length;
13+const scaleHits = (s: number) => sys.hits.map(h => ({ ...h, velocity: h.velocity.map(v => v * s) as [number, number, number], force: h.force.map(v => v * s) as [number, number, number] }));
14+const share = (r: ReturnType<typeof simulate>, t: number) => {
15+ const f = r.frameTimes.findIndex(x => x >= t); if (f < 0) return NaN;
16+ let hi = 0, tot = 0; for (let k = 0; k < nm; k++) { const e = r.modeEnergy![f * nm + k]; tot += e; if (modes.freqs[k] > fc) hi += e; } return hi / tot;
17+};
18+const ts = [0.002, 0.01, 0.03, 0.1, 0.3].filter(t => t < sys.duration);
19+for (const sc of process.argv.slice(5).map(Number)) {
20+ const big = simulate({ ...sys, hits: scaleHits(sc) }, { modes });
21+ const small = simulate({ ...sys, hits: scaleHits(sc * 1e-3) }, { modes });
22+ console.log(`scale ${sc}: share of energy above ${fc} Hz (linear -> actual): ` + ts.map(t => `${t}s ${(100*share(small, t)).toFixed(1)}%->${(100*share(big, t)).toFixed(1)}%`).join(" ") + ` ${big.aborted ?? ""}`);
23+}
tests/nonlin.tsadded+21−0View file
@@ -0,0 +1,21 @@
1+// How nonlinear is a run? Compare the listener signal with the same strike
2+// scaled down by 1000 (then scaled back up). Usage: tsx tests/nonlin.ts script.m [duration] [v1 v2 ...]
3+import { readFileSync, readdirSync } from "fs";
4+import { runScript } from "../src/sim/runScript";
5+import { simulate } from "../src/sim/solver";
6+import { relDiff, spectrum } from "../src/sim/analysis";
7+const helpers = readdirSync("matlab/helpers").map(name => ({ name, source: readFileSync("matlab/helpers/" + name, "utf8") }));
8+const { sys } = runScript(readFileSync(process.argv[2], "utf8"), helpers, () => {});
9+sys.duration = Number(process.argv[3] ?? 0.3);
10+const scales = process.argv.slice(4).map(Number);
11+for (const sc of scales.length ? scales : [1]) {
12+ const hits = sys.hits.map(h => ({ ...h, velocity: h.velocity.map(v => v * sc) as [number, number, number], force: h.force.map(v => v * sc) as [number, number, number] }));
13+ const big = simulate({ ...sys, hits });
14+ const small = simulate({ ...sys, hits: hits.map(h => ({ ...h, velocity: h.velocity.map(v => v * 1e-3) as [number, number, number], force: h.force.map(v => v * 1e-3) as [number, number, number] })) });
15+ const lin = small.audio.map(v => v * 1e3);
16+ const fs = big.sampleRate;
17+ const win = [0.01, 0.03, 0.1, 0.3].filter(t => t <= sys.duration).map(t => `${t}s:${relDiff(lin, big.audio, Math.round(t * fs)).toFixed(3)}`);
18+ // spectral centroid, first 50 ms
19+ const cen = (x: Float32Array) => { const s = spectrum(x, 0, Math.round(0.05 * fs)); let a = 0, b = 0; s.db.forEach((d, i) => { const p = 10 ** (d / 10); a += p * i * s.df * fs; b += p; }); return a / b; };
20+ console.log(`scale ${sc}: dt ${big.dt.toExponential(2)} drift ${big.energyDrift.toExponential(1)} ${big.aborted ?? ""} | nonlin diff ${win.join(" ")} | centroid lin ${cen(lin).toFixed(0)} Hz, actual ${cen(big.audio).toFixed(0)} Hz | ${(big.elapsedMs/1000).toFixed(1)}s`);
21+}
tests/solver.test.tsadded+171−0View file
@@ -0,0 +1,171 @@
1+import { readFileSync, readdirSync } from "fs";
2+import { describe, expect, it } from "vitest";
3+import { computeModes } from "../src/sim/modes";
4+import { numblEig, runScript } from "../src/sim/runScript";
5+import { frequencyErrorCents, simulate, stabilityLimit } from "../src/sim/solver";
6+import type { BeadSystem } from "../src/sim/types";
7+
8+const helpers = readdirSync("matlab/helpers").map(name => ({
9+ name,
10+ source: readFileSync(`matlab/helpers/${name}`, "utf8"),
11+}));
12+const build = (src: string) => runScript(src, helpers, () => {});
13+
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))).
16+const STRING = (extra: string) => `
17+N = 15; a = 0.01; m = 1e-3; f0 = 10; k = 1e4;
18+[pos, E] = bp_chain(N + 2, (N + 1) * a);
19+sys.pos = pos;
20+sys.mass = m;
21+sys.pinned = [1, N + 2];
22+sys.springs = bp_springs(pos, E, k, f0);
23+hit.beads = 5; hit.velocity = [0 0 1e-4];
24+sys.hits = hit;
25+sys.listen = struct('beads', 4);
26+sys.duration = 0.2;
27+${extra}
28+`;
29+
30+function 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+}
54+
55+describe("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);
58+
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+ });
66+
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+ });
74+
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+});
81+
82+describe("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;
86+
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+ });
94+
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+ });
100+
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+});
114+
115+describe("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+ });
132+
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+});
149+
150+describe("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+ });
158+
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+});
165+
166+describe("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+});
tsconfig.jsonadded+17−0View file
@@ -0,0 +1,17 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "lib": ["ES2022", "DOM", "DOM.Iterable", "WebWorker"],
5+ "module": "ESNext",
6+ "moduleResolution": "bundler",
7+ "jsx": "react-jsx",
8+ "strict": true,
9+ "noUnusedLocals": true,
10+ "noUnusedParameters": true,
11+ "skipLibCheck": true,
12+ "types": ["vite/client", "node"],
13+ "noEmit": true,
14+ "allowImportingTsExtensions": true
15+ },
16+ "include": ["src", "tests", "vite.config.ts"]
17+}
vite.config.tsadded+9−0View file
@@ -0,0 +1,9 @@
1+import react from "@vitejs/plugin-react";
2+import { defineConfig } from "vitest/config";
3+
4+export default defineConfig({
5+ base: "./",
6+ plugins: [react()],
7+ worker: { format: "es" },
8+ test: { testTimeout: 30000 },
9+});