Bead drum: hex-lattice round head, felt mallet contact, top view and cross-section; shared instrument page
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_012zsyegMfQnzJEzgMiFLzew
14 changed files+1244−261
README.mdmodified+12−1View file
@@ -22,10 +22,21 @@ npm run build # static site in dist/
2222 - 3D model view (three.js): Setup shows the beads at rest with supports (wall = clamped, translucent plane + ring = sliding), pinned beads (cubes), excitation (pluck shape, strike velocity arrows, hammer force arrows) and pickup (bead + axis arrow). Replay animates the saved frames ("Replay frames" per second, default 2000) with adjustable exaggeration, on a timeline shared with the 2D side/top views and space-time image.
2323 - Checks shown per run: energy balance (energy + dissipated - work in - initial) relative to peak energy, largest strain. Linear transverse and longitudinal modes (Sturm bisection on the tridiagonal stiffness) are marked on the spectrum.
2424
25-Tests (`test/string.test.ts`): modes vs the discrete-string formula; small mode-1 motion vs `cos(omega_1 t)`; energy balance with and without damping/forcing; decay at the requested T60.
25+String tests (`test/string.test.ts`): modes vs the discrete-string formula; small mode-1 motion vs `cos(omega_1 t)`; energy balance with and without damping/forcing; decay at the requested T60.
26+
27+## Bead drum
28+
29+- Round head: beads on a triangular (hex) lattice in the x-z plane, cut to a circle; every missing neighbour of an edge bead is a fixed rim point. Vibration along y.
30+- Inputs are continuum values: diameter, rings (spacing a = R / (rings + 1)), surface density, tension tau (N/m), in-plane stiffness E t (N/m). Lattice: spring tension f = tau a / sqrt(3), stiffness k = sqrt(3) E t / 2 (Poisson ratio 1/3), bead mass = density * sqrt(3)/2 a^2.
31+- Excitation: felt mallet (point mass, contact force K delta^p on a Gaussian patch, integrated with the beads; contact time and rebounds come out of the simulation), strike, force pulse, press (initial dent).
32+- Damping as for the string (dashpots calibrated to the springs' transverse stiffness f / a).
33+- Extra views: top-down displacement map and cross-section along the diameter through the strike point; mallet contact force.
34+- Presets: a fast tom with in-plane stiffness softened ~70x (about 8x slower than real time for 571 beads in Chromium) and a realistic Mylar one (~6x more steps).
35+- Tests (`test/drum.test.ts`): every bead has six springs; a small centre press rings within 1.3% of the ideal membrane fundamental 2.405 c / (2 pi R); the mallet bounces off and energy balances.
2636
2737 ## Layout
2838
2939 - `src/sim/` generic 3D bead-spring system, solver, worker.
3040 - `src/instruments/<id>/` model (parameters to system), presets, page. Register new instruments in `src/instruments/registry.ts`.
41+- `src/ui/instrumentPage.ts` the shared page (controls, 3D model, setup, run, audio, plots); an instrument supplies setup rows and extra cards.
3142 - `src/ui/` reusable controls form, line plot, spectrogram, chain displacement viewer, 3D scene, shared playback timeline, FFT, WAV.
src/instruments/bead-drum/model.tsadded+255−0View file
@@ -0,0 +1,255 @@
1+import type { Anchor, BeadSystem, Excitation } from '../../sim/types';
2+import { dampingFromT60 } from '../bead-string/model';
3+
4+// A round drumhead: beads on a triangular (hex) lattice in the x-z plane,
5+// cut to a circle, held by a rim of fixed lattice points just outside it.
6+// Transverse vibration is along y, as on the string.
7+
8+export type DrumExcitation = 'mallet' | 'strike' | 'force' | 'press';
9+
10+export interface DrumParams {
11+ // geometry
12+ diameter: number; // m, rim to rim
13+ rings: number; // beads from centre to edge (lattice spacing a = R / (rings + 1))
14+ // membrane properties (continuum values; the lattice is chosen to match them)
15+ density: number; // kg/m^2, surface density
16+ tension: number; // N/m, isotropic membrane tension
17+ inPlane: number; // N/m, in-plane stiffness E t
18+ // damping
19+ t60Low: number;
20+ fLow: number;
21+ t60High: number;
22+ fHigh: number;
23+ // excitation
24+ excitation: DrumExcitation;
25+ excR: number; // fraction of radius, 0 = centre
26+ excAngle: number; // degrees, in the x-z plane from +x
27+ excWidth: number; // m, Gaussian sigma of the contact patch (0 = one bead)
28+ malletMass: number; // kg
29+ malletSpeed: number; // m/s
30+ feltK: number; // N/m^p
31+ feltP: number;
32+ strikeVel: number; // m/s
33+ forcePeak: number; // N
34+ forceMs: number; // ms
35+ pressAmp: number; // m
36+ // pickup
37+ pickupMode: 'bead' | 'sum';
38+ pickupR: number;
39+ pickupAngle: number;
40+ pickupAxis: 'y' | 'x' | 'z';
41+ // simulation
42+ duration: number;
43+ fs: number;
44+ cfl: number;
45+ frameRate: number;
46+}
47+
48+export interface Lattice {
49+ a: number; // spacing (m)
50+ pos: [number, number][]; // bead positions (x, z)
51+ rim: [number, number][]; // anchor positions (x, z)
52+ edges: [number, number][]; // node pairs; node >= pos.length is rim[node - n]
53+}
54+
55+const key = (i: number, j: number) => `${i},${j}`;
56+
57+export function hexLattice(R: number, rings: number): Lattice {
58+ const a = R / (rings + 1);
59+ const at = (i: number, j: number): [number, number] => [a * (i + 0.5 * j), a * (Math.sqrt(3) / 2) * j];
60+ const inside = (i: number, j: number) => Math.hypot(...at(i, j)) <= R - 0.5 * a + 1e-9;
61+ const index = new Map<string, number>();
62+ const pos: [number, number][] = [];
63+ const M = rings + 3;
64+ for (let j = -M; j <= M; j++)
65+ for (let i = -M - Math.abs(j); i <= M + Math.abs(j); i++)
66+ if (inside(i, j)) {
67+ index.set(key(i, j), pos.length);
68+ pos.push(at(i, j));
69+ }
70+ const n = pos.length;
71+ const rimIndex = new Map<string, number>();
72+ const rim: [number, number][] = [];
73+ const edges: [number, number][] = [];
74+ const nb = [[1, 0], [0, 1], [-1, 1]];
75+ const all = [[1, 0], [0, 1], [-1, 1], [-1, 0], [0, -1], [1, -1]];
76+ for (const [k, b] of index) {
77+ const [i, j] = k.split(',').map(Number);
78+ // bead-bead springs: each pair once
79+ for (const [di, dj] of nb) {
80+ const o = index.get(key(i + di, j + dj));
81+ if (o !== undefined) edges.push([b, o]);
82+ }
83+ // bead-rim springs: every missing neighbour is a rim point
84+ for (const [di, dj] of all) {
85+ const kk = key(i + di, j + dj);
86+ if (index.has(kk)) continue;
87+ let r = rimIndex.get(kk);
88+ if (r === undefined) {
89+ r = rim.length;
90+ rimIndex.set(kk, r);
91+ rim.push(at(i + di, j + dj));
92+ }
93+ edges.push([b, n + r]);
94+ }
95+ }
96+ return { a, pos, rim, edges };
97+}
98+
99+export interface DrumDerived {
100+ lattice: Lattice;
101+ n: number;
102+ beadMass: number;
103+ springK: number;
104+ springTension: number;
105+ restLength: number;
106+ prestretch: number;
107+ waveSpeedT: number;
108+ waveSpeedP: number; // in-plane longitudinal (P) wave speed
109+ idealF01: number; // ideal clamped membrane fundamental
110+ sigma0: number;
111+ sigma1: number;
112+ dashpot: number;
113+ warnings: string[];
114+}
115+
116+// Triangular spring lattice, spacing a, spring stiffness k, spring tension f:
117+// membrane tension tau = sqrt(3) f / a, 2D Young's modulus E t = 2 k / sqrt(3)
118+// (Poisson ratio 1/3), area per bead sqrt(3)/2 a^2.
119+export function deriveDrum(p: DrumParams): DrumDerived {
120+ const R = p.diameter / 2;
121+ const lattice = hexLattice(R, Math.max(1, Math.round(p.rings)));
122+ const a = lattice.a;
123+ const beadMass = p.density * (Math.sqrt(3) / 2) * a * a;
124+ const springK = (Math.sqrt(3) / 2) * p.inPlane;
125+ const springTension = (p.tension * a) / Math.sqrt(3);
126+ const restLength = a - springTension / springK;
127+ const warnings: string[] = [];
128+ if (!(restLength > 0)) warnings.push('Tension too high for this in-plane stiffness: the springs would need a negative rest length.');
129+ const { sigma0, sigma1, note } = dampingFromT60(p);
130+ if (note) warnings.push(note);
131+ const waveSpeedT = Math.sqrt(Math.max(p.tension, 0) / p.density);
132+ const lambdaPlus2Mu = (3 * Math.sqrt(3) * springK) / 4;
133+ return {
134+ lattice,
135+ n: lattice.pos.length,
136+ beadMass,
137+ springK,
138+ springTension,
139+ restLength,
140+ prestretch: a / restLength - 1,
141+ waveSpeedT,
142+ waveSpeedP: Math.sqrt(lambdaPlus2Mu / p.density),
143+ idealF01: (2.404826 * waveSpeedT) / (2 * Math.PI * R),
144+ sigma0,
145+ sigma1,
146+ // each spring's transverse stiffness is f / a, so this damps transverse modes by sigma1 omega^2
147+ dashpot: 2 * sigma1 * (springTension > 0 ? springTension / a : springK),
148+ warnings,
149+ };
150+}
151+
152+export function buildDrum(d: DrumDerived): BeadSystem {
153+ const { pos, rim, edges } = d.lattice;
154+ const n = pos.length;
155+ const rest = new Float64Array(3 * n);
156+ pos.forEach(([x, z], i) => {
157+ rest[3 * i] = x;
158+ rest[3 * i + 2] = z;
159+ });
160+ const anchors: Anchor[] = rim.map(([x, z]) => ({ pos: [x, 0, z], slideBead: -1 }));
161+ const ns = edges.length;
162+ return {
163+ n,
164+ rest,
165+ mass: new Float64Array(n).fill(d.beadMass),
166+ pinned: new Uint8Array(n),
167+ anchors,
168+ si: Int32Array.from(edges.map((e) => e[0])),
169+ sj: Int32Array.from(edges.map((e) => e[1])),
170+ sk: new Float64Array(ns).fill(d.springK),
171+ sL0: new Float64Array(ns).fill(d.restLength),
172+ sc: new Float64Array(ns).fill(d.dashpot),
173+ sigma0: d.sigma0,
174+ };
175+}
176+
177+const polar = (R: number, frac: number, deg: number): [number, number] => {
178+ const th = (deg * Math.PI) / 180;
179+ return [R * frac * Math.cos(th), R * frac * Math.sin(th)];
180+};
181+
182+export function nearestBead(d: DrumDerived, x: number, z: number): number {
183+ let best = 0;
184+ let bd = Infinity;
185+ d.lattice.pos.forEach(([px, pz], i) => {
186+ const q = (px - x) ** 2 + (pz - z) ** 2;
187+ if (q < bd) {
188+ bd = q;
189+ best = i;
190+ }
191+ });
192+ return best;
193+}
194+
195+// Gaussian patch around the excitation point; width 0 = the nearest bead
196+function patch(p: DrumParams, d: DrumDerived): Float64Array {
197+ const [cx, cz] = polar(p.diameter / 2, p.excR, p.excAngle);
198+ const n = d.n;
199+ const w = new Float64Array(n);
200+ if (!(p.excWidth > 0)) {
201+ w[nearestBead(d, cx, cz)] = 1;
202+ return w;
203+ }
204+ d.lattice.pos.forEach(([x, z], i) => {
205+ const r2 = ((x - cx) ** 2 + (z - cz) ** 2) / (p.excWidth * p.excWidth);
206+ w[i] = r2 < 9 ? Math.exp(-0.5 * r2) : 0;
207+ });
208+ if (!w.some((v) => v > 0)) w[nearestBead(d, cx, cz)] = 1;
209+ return w;
210+}
211+
212+export function excitationPoint(p: DrumParams): [number, number] {
213+ return polar(p.diameter / 2, p.excR, p.excAngle);
214+}
215+
216+export function buildDrumExcitation(p: DrumParams, d: DrumDerived): Excitation {
217+ const n = d.n;
218+ const u0 = new Float64Array(3 * n);
219+ const v0 = new Float64Array(3 * n);
220+ const w = patch(p, d);
221+ let max = 0;
222+ let sum = 0;
223+ for (let i = 0; i < n; i++) {
224+ max = Math.max(max, w[i]);
225+ sum += w[i];
226+ }
227+ // everything acts downward (-y), onto the head from above
228+ if (p.excitation === 'press') {
229+ for (let i = 0; i < n; i++) u0[3 * i + 1] = (-p.pressAmp * w[i]) / max;
230+ return { u0, v0 };
231+ }
232+ if (p.excitation === 'strike') {
233+ for (let i = 0; i < n; i++) v0[3 * i + 1] = (-p.strikeVel * w[i]) / max;
234+ return { u0, v0 };
235+ }
236+ const weights = w.map((v) => v / sum);
237+ if (p.excitation === 'force') {
238+ const shape = new Float64Array(3 * n);
239+ for (let i = 0; i < n; i++) shape[3 * i + 1] = -weights[i];
240+ return { u0, v0, force: { shape, F0: p.forcePeak, tau: p.forceMs / 1000 } };
241+ }
242+ return {
243+ u0,
244+ v0,
245+ mallet: { weights, dir: [0, -1, 0], mass: p.malletMass, K: p.feltK, p: p.feltP, v0: p.malletSpeed },
246+ };
247+}
248+
249+export function buildDrumPickup(p: DrumParams, d: DrumDerived): Float64Array {
250+ const w = new Float64Array(3 * d.n);
251+ const axis = p.pickupAxis === 'x' ? 0 : p.pickupAxis === 'y' ? 1 : 2;
252+ if (p.pickupMode === 'sum') for (let i = 0; i < d.n; i++) w[3 * i + axis] = 1;
253+ else w[3 * nearestBead(d, ...polar(p.diameter / 2, p.pickupR, p.pickupAngle)) + axis] = 1;
254+ return w;
255+}
src/instruments/bead-drum/page.tsadded+206−0View file
@@ -0,0 +1,206 @@
1+import type { Section } from '../../ui/controls';
2+import { LinePlot, css, fmtNum } from '../../ui/plot';
3+import { MembraneView } from '../../ui/membraneview';
4+import { mountInstrumentPage, type Built } from '../../ui/instrumentPage';
5+import {
6+ buildDrum,
7+ buildDrumExcitation,
8+ buildDrumPickup,
9+ deriveDrum,
10+ excitationPoint,
11+ nearestBead,
12+ type DrumParams,
13+} from './model';
14+import { DRUM_PRESETS } from './presets';
15+
16+const isExc = (t: DrumParams['excitation']) => (p: DrumParams) => p.excitation === t;
17+
18+const SECTIONS: Section<DrumParams>[] = [
19+ {
20+ title: 'Geometry',
21+ params: [
22+ { key: 'diameter', label: 'Diameter', kind: 'number', unit: 'cm', scale: 100 },
23+ { key: 'rings', label: 'Rings of beads', kind: 'int', min: 1, help: 'Beads from the centre to the rim. The lattice spacing is radius / (rings + 1); bead count grows as about 3 × rings².' },
24+ ],
25+ },
26+ {
27+ title: 'Membrane',
28+ params: [
29+ { key: 'density', label: 'Surface density', kind: 'number', unit: 'kg/m²', help: 'Mass per area. Mylar drumheads: about 0.26 kg/m².' },
30+ { key: 'tension', label: 'Tension', kind: 'number', unit: 'N/m', help: 'Force per unit length in the head (same in every direction). Sets the pitch.' },
31+ { key: 'inPlane', label: 'In-plane stiffness E·t', kind: 'number', unit: 'N/m', help: "Young's modulus times thickness. Mylar 0.19 mm: about 8.6e5 N/m. Lower values run faster but exaggerate nonlinearity." },
32+ ],
33+ },
34+ {
35+ title: 'Damping (0 = off)',
36+ params: [
37+ { key: 't60Low', label: 'T60 low', kind: 'number', unit: 's' },
38+ { key: 'fLow', label: 'at', kind: 'number', unit: 'Hz' },
39+ { key: 't60High', label: 'T60 high', kind: 'number', unit: 's' },
40+ { key: 'fHigh', label: 'at', kind: 'number', unit: 'Hz' },
41+ ],
42+ },
43+ {
44+ title: 'Excitation',
45+ params: [
46+ {
47+ key: 'excitation',
48+ label: 'Type',
49+ kind: 'select',
50+ options: [
51+ ['mallet', 'Mallet (felt contact)'],
52+ ['strike', 'Strike (initial velocity)'],
53+ ['force', 'Force pulse'],
54+ ['press', 'Press (initial dent, released)'],
55+ ],
56+ },
57+ { key: 'excR', label: 'Position', kind: 'number', unit: '× radius', help: '0 = centre, 1 = rim.' },
58+ { key: 'excAngle', label: 'Angle', kind: 'number', unit: '°' },
59+ { key: 'excWidth', label: 'Patch width', kind: 'number', unit: 'mm', scale: 1000, help: 'Gaussian width of the contact patch; 0 = a single bead.' },
60+ { key: 'malletMass', label: 'Mallet mass', kind: 'number', unit: 'g', scale: 1000, show: isExc('mallet') },
61+ { key: 'malletSpeed', label: 'Mallet speed', kind: 'number', unit: 'm/s', show: isExc('mallet') },
62+ { key: 'feltK', label: 'Felt stiffness K', kind: 'number', unit: 'N/mᵖ', show: isExc('mallet'), help: 'Contact force K·δ^p, δ = felt compression. Harder felt = larger K: shorter contact, brighter sound.' },
63+ { key: 'feltP', label: 'Felt exponent p', kind: 'number', show: isExc('mallet'), help: 'Typical felt: 2 to 3.5.' },
64+ { key: 'strikeVel', label: 'Velocity', kind: 'number', unit: 'm/s', show: isExc('strike') },
65+ { key: 'forcePeak', label: 'Peak force', kind: 'number', unit: 'N', show: isExc('force') },
66+ { key: 'forceMs', label: 'Duration', kind: 'number', unit: 'ms', show: isExc('force') },
67+ { key: 'pressAmp', label: 'Depth', kind: 'number', unit: 'mm', scale: 1000, show: isExc('press') },
68+ ],
69+ },
70+ {
71+ title: 'Pickup (acceleration)',
72+ params: [
73+ { key: 'pickupMode', label: 'Pickup', kind: 'select', options: [['sum', 'Sum over all beads'], ['bead', 'One bead']], help: 'The sum of accelerations is the monopole term of the radiated sound.' },
74+ { key: 'pickupR', label: 'Position', kind: 'number', unit: '× radius', show: (p) => p.pickupMode === 'bead' },
75+ { key: 'pickupAngle', label: 'Angle', kind: 'number', unit: '°', show: (p) => p.pickupMode === 'bead' },
76+ { key: 'pickupAxis', label: 'Component', kind: 'select', options: [['y', 'y (transverse)'], ['x', 'x (in-plane)'], ['z', 'z (in-plane)']] },
77+ ],
78+ },
79+ {
80+ title: 'Simulation',
81+ params: [
82+ { key: 'duration', label: 'Duration', kind: 'number', unit: 's' },
83+ { key: 'fs', label: 'Sample rate', kind: 'select', options: [['44100', '44.1 kHz'], ['48000', '48 kHz'], ['96000', '96 kHz']] },
84+ { key: 'cfl', label: 'Step safety', kind: 'number', help: 'Fraction of the leapfrog stability limit dt = 2/ω_max.' },
85+ { key: 'frameRate', label: 'Replay frames', kind: 'number', unit: 'per s', help: 'How often displacements of all beads are saved for replay.' },
86+ ],
87+ },
88+];
89+
90+// separate contact intervals in the contact-force record
91+function contacts(cf: Float32Array, fs: number): { start: number; ms: number; peak: number }[] {
92+ const out: { start: number; ms: number; peak: number }[] = [];
93+ let i = 0;
94+ while (i < cf.length) {
95+ if (cf[i] > 0) {
96+ const a = i;
97+ let peak = 0;
98+ while (i < cf.length && cf[i] > 0) peak = Math.max(peak, cf[i++]);
99+ out.push({ start: a / fs, ms: ((i - a) / fs) * 1000, peak });
100+ } else i++;
101+ }
102+ return out;
103+}
104+
105+export function mountBeadDrum(root: HTMLElement) {
106+ const build = (p: DrumParams): Built | null => {
107+ const d = deriveDrum(p);
108+ if (!(d.restLength > 0) || d.n < 1) return null;
109+ const sys = buildDrum(d);
110+ const exc = buildDrumExcitation(p, d);
111+ const pickup = buildDrumPickup(p, d);
112+ return {
113+ sys,
114+ exc,
115+ pickup,
116+ scene: {
117+ sys,
118+ exc,
119+ excKind: p.excitation,
120+ pickup,
121+ pickupLabel: `${p.pickupMode === 'sum' ? 'sum over all beads' : 'one bead'}, ${p.pickupAxis} acceleration`,
122+ hoop: { radius: p.diameter / 2 },
123+ },
124+ };
125+ };
126+
127+ mountInstrumentPage<DrumParams>(root, {
128+ title: 'Bead drum',
129+ wavName: 'bead-drum.wav',
130+ presets: DRUM_PRESETS,
131+ sections: SECTIONS,
132+ setup(p) {
133+ const d = deriveDrum(p);
134+ const rows: [string, string][] = [
135+ ['Beads', `${d.n} (rim points ${d.lattice.rim.length})`],
136+ ['Lattice spacing', `${fmtNum(d.lattice.a * 1000)} mm`],
137+ ['Bead mass', `${fmtNum(d.beadMass * 1000)} g`],
138+ ['Spring stiffness, tension', `${fmtNum(d.springK)} N/m, ${fmtNum(d.springTension)} N`],
139+ ['Rest length L0', `${fmtNum(d.restLength * 1000)} mm (pre-stretch ${fmtNum(d.prestretch * 100)}%)`],
140+ ['Wave speed (transverse, in-plane)', `${fmtNum(d.waveSpeedT)}, ${fmtNum(d.waveSpeedP)} m/s`],
141+ ['Ideal membrane fundamental', `${fmtNum(d.idealF01)} Hz`],
142+ ];
143+ const warnings = [...d.warnings];
144+ if (p.excR >= 0.97) warnings.push('The excitation point is at the rim, where the head barely moves.');
145+ if (d.n > 3000) warnings.push(`${d.n} beads: expect a slow run.`);
146+ return { rows, warnings, built: build(p) };
147+ },
148+ summaryRows(p, _b, r) {
149+ const rows: [string, string][] = [['Largest spring strain', `${fmtNum(r.maxStretch * 100)}% (pre-stretch ${fmtNum(deriveDrum(p).prestretch * 100)}%)`]];
150+ if (r.contactForce) {
151+ const c = contacts(r.contactForce, r.fs);
152+ rows.push([
153+ 'Mallet contact',
154+ c.length
155+ ? `${fmtNum(c[0].ms)} ms, peak ${fmtNum(c[0].peak)} N${c.length > 1 ? ` (+${c.length - 1} rebound${c.length > 2 ? 's' : ''})` : ''}`
156+ : 'no contact',
157+ ]);
158+ }
159+ return rows;
160+ },
161+ extraCards: `
162+ <section class="card" data-r="out" hidden><h2>Drumhead (top view and cross-section)</h2><div data-r="membrane"></div></section>
163+ <section class="card" data-r="out-contact" hidden>
164+ <h2>Mallet contact force</h2>
165+ <div class="canvas-wrap short"><canvas data-r="contact"></canvas></div>
166+ <p class="note" data-r="contact-note"></p>
167+ </section>`,
168+ mountExtras(q, timeline) {
169+ const view = new MembraneView(q('membrane'), timeline);
170+ const contactPlot = new LinePlot(q('contact'));
171+ return {
172+ show(p, _built, r) {
173+ const d = deriveDrum(p);
174+ view.set(
175+ {
176+ pos: d.lattice.pos,
177+ radius: p.diameter / 2,
178+ spacing: d.lattice.a,
179+ strike: excitationPoint(p),
180+ pickup: p.pickupMode === 'bead' ? nearestBead(d, (p.diameter / 2) * p.pickupR * Math.cos((p.pickupAngle * Math.PI) / 180), (p.diameter / 2) * p.pickupR * Math.sin((p.pickupAngle * Math.PI) / 180)) : null,
181+ sectionAngle: p.excR > 0.01 ? p.excAngle : 0,
182+ },
183+ r.frames,
184+ );
185+ const card = q('out-contact');
186+ card.hidden = !r.contactForce;
187+ if (!r.contactForce) return;
188+ const cf = r.contactForce;
189+ const c = contacts(cf, r.fs);
190+ const end = c.length ? c[c.length - 1].start + c[c.length - 1].ms / 1000 : 0.01;
191+ const nShow = Math.min(cf.length, Math.ceil((end * 1.3 + 0.002) * r.fs));
192+ const t = new Float64Array(nShow);
193+ for (let i = 0; i < nShow; i++) t[i] = (i / r.fs) * 1000;
194+ contactPlot.set([{ x: t, y: cf.subarray(0, nShow), color: css('--series-2'), label: 'force (N)' }], {
195+ xLabel: 'time (ms)',
196+ yLabel: 'contact force (N)',
197+ yRange: [0, Math.max(...cf.subarray(0, nShow)) * 1.08 || 1],
198+ });
199+ q('contact-note').textContent = c.length
200+ ? `${c.length} contact${c.length > 1 ? 's' : ''}: ` + c.map((k) => `${fmtNum(k.ms)} ms from ${fmtNum(k.start * 1000)} ms, peak ${fmtNum(k.peak)} N`).join('; ') + `. During contact the head takes energy from the mallet and gives much of it back as the mallet rebounds, so the energy plot spikes first; the mallet then leaves at constant speed (no gravity).`
201+ : 'The mallet never touched the head.';
202+ },
203+ };
204+ },
205+ });
206+}
src/instruments/bead-drum/presets.tsadded+70−0View file
@@ -0,0 +1,70 @@
1+import type { DrumParams } from './model';
2+
3+// Mylar-like head: 0.19 mm film, surface density 0.26 kg/m^2, E t ~ 8.6e5 N/m.
4+const fast: DrumParams = {
5+ diameter: 0.3,
6+ rings: 12,
7+ density: 0.26,
8+ tension: 1500, // c = 76 m/s, ideal f01 ~ 194 Hz
9+ inPlane: 12000, // softened: in-plane waves only ~3x faster than transverse
10+ t60Low: 1.2,
11+ fLow: 200,
12+ t60High: 0.25,
13+ fHigh: 3000,
14+ excitation: 'mallet',
15+ excR: 0.35,
16+ excAngle: 0,
17+ excWidth: 0.008,
18+ malletMass: 0.02,
19+ malletSpeed: 3,
20+ feltK: 1e8,
21+ feltP: 2.5,
22+ strikeVel: 1,
23+ forcePeak: 50,
24+ forceMs: 1,
25+ pressAmp: 0.003,
26+ pickupMode: 'sum',
27+ pickupR: 0.5,
28+ pickupAngle: 120,
29+ pickupAxis: 'y',
30+ duration: 1.5,
31+ fs: 48000,
32+ cfl: 0.5,
33+ frameRate: 2000,
34+};
35+
36+export interface DrumPreset {
37+ name: string;
38+ description: string;
39+ params: DrumParams;
40+}
41+
42+export const DRUM_PRESETS: DrumPreset[] = [
43+ {
44+ name: 'Tom, fast (softened in-plane)',
45+ description:
46+ '30 cm head, ~470 beads, 1500 N/m, ideal fundamental ~194 Hz. In-plane springs softened about 70x for speed: small-amplitude pitch is unchanged, but large hits glide and couple more than a real head.',
47+ params: fast,
48+ },
49+ {
50+ name: 'Tom, realistic Mylar (slow)',
51+ description:
52+ 'Same head with Mylar-like in-plane stiffness (E·t = 8.6e5 N/m). Accurate in-plane behaviour, about 6x more time steps.',
53+ params: { ...fast, inPlane: 8.6e5 },
54+ },
55+ {
56+ name: 'Hard hit near the rim',
57+ description: 'Fast head, light hard mallet at 4 m/s close to the rim: short contact, bright attack, many modes.',
58+ params: { ...fast, excR: 0.8, malletMass: 0.01, malletSpeed: 4, feltK: 1e9, feltP: 2.8, excWidth: 0.005 },
59+ },
60+ {
61+ name: 'Soft mallet in the centre',
62+ description: 'Heavy soft mallet at the centre: long contact, mostly the axisymmetric modes, dark sound.',
63+ params: { ...fast, excR: 0, malletMass: 0.05, malletSpeed: 1.5, feltK: 2e6, feltP: 2.3, excWidth: 0.015 },
64+ },
65+ {
66+ name: 'Small coarse drum',
67+ description: '6 rings (~130 beads): runs quickly; the coarse lattice is audibly more inharmonic than a real head.',
68+ params: { ...fast, diameter: 0.2, rings: 6, excWidth: 0.01 },
69+ },
70+];
src/instruments/bead-string/model.tsmodified+3−1View file
@@ -65,7 +65,9 @@ export interface Derived {
6565
6666 // Linear damping sigma(omega) = sigma0 + sigma1 omega^2 for transverse modes,
6767 // fitted to two decay times T60 = 3 ln(10) / sigma.
68-export function dampingFromT60(p: StringParams): { sigma0: number; sigma1: number; note?: string } {
68+export function dampingFromT60(
69+ p: Pick<StringParams, 't60Low' | 'fLow' | 't60High' | 'fHigh'>,
70+): { sigma0: number; sigma1: number; note?: string } {
6971 if (!(p.t60Low > 0) && !(p.t60High > 0)) return { sigma0: 0, sigma1: 0 };
7072 const ln1000 = 3 * Math.LN10;
7173 if (!(p.t60High > 0) || !(p.fHigh > p.fLow)) return { sigma0: ln1000 / p.t60Low, sigma1: 0 };
src/instruments/bead-string/page.tsmodified+80−244View file
@@ -1,14 +1,7 @@
1-import { renderControls, type Section } from '../../ui/controls';
2-import { LinePlot, css, fmtNum } from '../../ui/plot';
3-import { Spectrogram } from '../../ui/spectrogram';
1+import type { Section } from '../../ui/controls';
2+import { css, fmtNum } from '../../ui/plot';
43 import { ChainView } from '../../ui/chainview';
5-import { BeadScene, type SceneSpec } from '../../ui/scene3d';
6-import { Timeline } from '../../ui/timeline';
7-import { spectrum } from '../../ui/fft';
8-import { encodeWav } from '../../ui/wav';
9-import { omegaMaxBound } from '../../sim/solver';
10-import { runInWorker, type Running } from '../../sim/run';
11-import type { BeadSystem, Excitation, SimResult } from '../../sim/types';
4+import { mountInstrumentPage, type Built } from '../../ui/instrumentPage';
125 import { buildExcitation, buildPickup, buildSystem, derive, parsePins, type StringParams } from './model';
136 import { PRESETS } from './presets';
147
@@ -94,242 +87,85 @@ const SECTIONS: Section<StringParams>[] = [
9487 ];
9588
9689 export function mountBeadString(root: HTMLElement) {
97- root.innerHTML = `
98- <div class="instrument">
99- <aside class="panel">
100- <a href="#/" class="back">← All instruments</a>
101- <h1>String of beads</h1>
102- <label class="ctl preset"><span class="ctl-name">Preset</span><select data-r="preset"></select></label>
103- <p class="sub" data-r="preset-desc"></p>
104- <div data-r="controls"></div>
105- <div class="run-row">
106- <button class="btn" data-r="run">Run simulation</button>
107- <button class="btn-small" data-r="cancel" hidden>Cancel</button>
108- </div>
109- <progress data-r="progress" max="1" value="0" hidden></progress>
110- </aside>
111- <main class="results">
112- <section class="card"><h2>Model</h2><div data-r="scene"></div></section>
113- <section class="card"><h2>Setup</h2><div data-r="derived" class="stats"></div><ul data-r="warnings" class="warnings"></ul></section>
114- <section class="card" data-r="out" hidden>
115- <h2>Run</h2><div data-r="summary" class="stats"></div>
116- <div class="audio-row"><audio controls data-r="audio"></audio><a class="btn-small" data-r="download" download="bead-string.wav">Download .wav</a></div>
117- <p class="note" data-r="audio-note"></p>
118- </section>
119- <section class="card" data-r="out" hidden>
120- <div class="card-head"><h2>Pickup acceleration</h2>
121- <div class="seg" data-r="wave-range"><button data-v="0" class="on">Full</button><button data-v="0.05">First 50 ms</button><button data-v="0.005">First 5 ms</button></div></div>
122- <div class="canvas-wrap"><canvas data-r="wave"></canvas></div>
123- </section>
124- <section class="card" data-r="out" hidden>
125- <div class="card-head"><h2>Spectrogram</h2>
126- <label class="inline">up to <select data-r="fmax"><option>1000</option><option>2500</option><option selected>5000</option><option>10000</option><option>24000</option></select> Hz</label></div>
127- <div class="canvas-wrap tall"><canvas data-r="spec"></canvas></div>
128- </section>
129- <section class="card" data-r="out" hidden>
130- <div class="card-head"><h2>Spectrum</h2>
131- <div class="legend"><span><i class="sw" style="background:var(--series-1)"></i>pickup signal</span><span><i class="dash" style="border-color:var(--series-1)"></i>linear transverse modes</span><span><i class="dash" style="border-color:var(--series-2)"></i>linear longitudinal modes</span></div></div>
132- <div class="canvas-wrap"><canvas data-r="spectrum"></canvas></div>
133- </section>
134- <section class="card" data-r="out" hidden><h2>Displacements (2D)</h2><div data-r="viewer"></div></section>
135- <section class="card" data-r="out" hidden>
136- <h2>Energy</h2>
137- <div class="canvas-wrap short"><canvas data-r="energy"></canvas></div>
138- <p class="note">Kinetic plus spring energy, minus the energy stored by pre-tension at equilibrium. Without damping or forcing it should stay flat; the run summary gives the error of the energy balance (energy + dissipated − work in − initial).</p>
139- </section>
140- </main>
141- </div>`;
142- const q = <T extends HTMLElement>(r: string) => root.querySelector(`[data-r="${r}"]`) as T;
143- const presetSel = q<HTMLSelectElement>('preset');
144- PRESETS.forEach((p, i) => presetSel.add(new Option(p.name, String(i))));
145- const presetDesc = q<HTMLParagraphElement>('preset-desc');
146- const runBtn = q<HTMLButtonElement>('run');
147- const cancelBtn = q<HTMLButtonElement>('cancel');
148- const progress = q<HTMLProgressElement>('progress');
149-
150- const controls = renderControls(q('controls'), SECTIONS, PRESETS[0].params, () => updateDerived());
151- presetDesc.textContent = PRESETS[0].description;
152- presetSel.addEventListener('change', () => {
153- const pr = PRESETS[Number(presetSel.value)];
154- presetDesc.textContent = pr.description;
155- controls.set(pr.params);
156- });
157-
158- const wave = new LinePlot(q('wave'));
159- const spec = new Spectrogram(q('spec'));
160- const specPlot = new LinePlot(q('spectrum'));
161- const energyPlot = new LinePlot(q('energy'));
162- const timeline = new Timeline();
163- const viewer = new ChainView(q('viewer'), timeline);
164- const scene = new BeadScene(q('scene'), timeline);
165- let waveRange = 0;
166- let last: { p: StringParams; r: SimResult } | null = null;
167- let running: Running | null = null;
168- let audioUrl: string | null = null;
169-
170- const sceneSpec = (p: StringParams, sys: BeadSystem, exc: Excitation): SceneSpec => ({
171- sys,
172- exc,
173- excKind: p.excitation,
174- pickup: buildPickup(p),
175- pickupLabel: `${p.pickupMode === 'sum' ? 'sum over all beads' : `bead ${Math.round(p.pickupBead)}`}, ${p.pickupAxis} acceleration`,
176- });
177-
178- function updateDerived() {
179- const p = controls.get();
180- const d = derive(p);
181- const tr = d.transverseHz;
182- const rows: [string, string][] = [
183- ['Length (wall to wall)', `${fmtNum(d.length)} m`],
184- ['Rest length L0', `${fmtNum(d.restLength * 1000)} mm (pre-stretch ${fmtNum(d.prestretch * 100)}%)`],
185- ['Linear density', `${fmtNum(d.linearDensity * 1000)} g/m`],
186- ['Wave speed (transverse, longitudinal)', `${fmtNum(d.waveSpeedT)}, ${fmtNum(d.waveSpeedL)} m/s`],
187- ['Lowest transverse mode', tr.length ? `${fmtNum(tr[0])} Hz` : 'none'],
188- ['Overtone ratios', tr.slice(1, 7).map((f) => fmtNum(f / tr[0])).join(', ') || '–'],
189- ['Lowest longitudinal mode', d.longitudinalHz.length ? `${fmtNum(d.longitudinalHz[0])} Hz` : 'none'],
190- ];
191- const warnings = [...d.warnings];
192- if (d.restLength > 0) {
193- const sys = buildSystem(p, d);
194- scene.showSetup(sceneSpec(p, sys, buildExcitation(p)));
195- const wmax = omegaMaxBound(sys);
196- const spp = Math.max(1, Math.ceil(1 / (p.fs * ((p.cfl * 2) / wmax))));
197- rows.push(['Highest frequency bound', `${fmtNum(wmax / (2 * Math.PI))} Hz`]);
198- rows.push(['Time step', `${fmtNum(1 / (p.fs * spp))} s (${spp} per sample)`]);
199- const work = spp * p.fs * p.duration * sys.si.length;
200- rows.push(['Spring updates', fmtNum(work)]);
201- const frameMB = (Math.round(p.duration * p.frameRate) * 3 * p.n * 4) / 1e6;
202- if (frameMB > 400) warnings.push(`Displacement frames would take ${fmtNum(frameMB)} MB; lower the frame rate or duration.`);
203- if (work > 3e9) warnings.push('This run is large (more than 3e9 spring updates) and may take minutes.');
204- }
205- const pins = parsePins(p.pins, p.n);
206- if (pins.includes(Math.round(p.excBead) - 1)) warnings.push('The excited bead is pinned.');
207- if (p.pickupMode === 'bead' && pins.includes(Math.round(p.pickupBead) - 1)) warnings.push('The pickup bead is pinned: its acceleration is zero.');
208- if (tr.length && tr[0] < 20) warnings.push(`The lowest mode (${fmtNum(tr[0])} Hz) is below hearing; you will mostly hear overtones.`);
209- q('derived').innerHTML = rows.map(([k, v]) => `<div><span>${k}</span><b>${v}</b></div>`).join('');
210- q('warnings').innerHTML = warnings.map((w) => `<li>${w}</li>`).join('');
211- runBtn.disabled = !(d.restLength > 0) || !(p.n >= 1);
212- }
213- updateDerived();
214-
215- q('wave-range').addEventListener('click', (e) => {
216- const b = (e.target as HTMLElement).closest('button');
217- if (!b) return;
218- waveRange = Number(b.dataset.v);
219- q('wave-range').querySelectorAll('button').forEach((x) => x.classList.toggle('on', x === b));
220- drawWave();
221- });
222- q<HTMLSelectElement>('fmax').addEventListener('change', (e) => {
223- spec.fmax = Number((e.target as HTMLSelectElement).value);
224- spec.draw();
225- drawSpectrum();
226- });
227-
228- function drawWave() {
229- if (!last) return;
230- const { r } = last;
231- const t = new Float64Array(r.audio.length);
232- for (let i = 0; i < t.length; i++) t[i] = i / r.fs;
233- wave.set([{ x: t, y: r.audio, color: css('--series-1'), label: 'a' }], {
234- xLabel: 'time (s)',
235- yLabel: 'acceleration (m/s²)',
236- xRange: waveRange ? [0, Math.min(waveRange, t[t.length - 1])] : undefined,
237- });
238- }
239-
240- function drawSpectrum() {
241- if (!last) return;
242- const { p, r } = last;
243- const d = derive(p);
244- const s = spectrum(r.audio, r.fs);
245- let top = -Infinity;
246- for (const v of s.db) top = Math.max(top, v);
247- const db = s.db.map((v) => Math.max(v - top, -100));
248- const fmax = Math.min(spec.fmax, r.fs / 2);
249- const markers = [
250- ...d.transverseHz.map((f, i) => ({ x: f, color: css('--series-1'), label: `transverse mode ${i + 1}: ${fmtNum(f)} Hz` })),
251- ...d.longitudinalHz.map((f, i) => ({ x: f, color: css('--series-2'), label: `longitudinal mode ${i + 1}: ${fmtNum(f)} Hz` })),
252- ];
253- specPlot.set([{ x: s.freq, y: db, color: css('--series-1'), label: 'level (dB)' }], {
254- xLabel: 'frequency (Hz)',
255- yLabel: 'level (dB re max)',
256- xRange: [0, fmax],
257- yRange: [-100, 5],
258- markers,
259- });
260- }
261-
262- function show(p: StringParams, sys: BeadSystem, exc: Excitation, r: SimResult) {
263- last = { p, r };
264- root.querySelectorAll<HTMLElement>('[data-r="out"]').forEach((e) => (e.hidden = false));
265- const d = derive(p);
266- let peak = 0;
267- for (const v of r.audio) peak = Math.max(peak, Math.abs(v));
268- const emax = Math.max(...r.energy);
269- let bal = 0;
270- for (const b of r.balance) bal = Math.max(bal, Math.abs(b));
271- const rt = p.duration / r.wallSeconds;
272- q('summary').innerHTML = [
273- ['Time step', `${fmtNum(r.dt)} s, ${r.stepsPerSample} per sample`],
274- ['Wall time', `${fmtNum(r.wallSeconds)} s (${rt >= 1 ? fmtNum(rt) + '× faster than' : fmtNum(1 / rt) + '× slower than'} real time)`],
275- ['Energy balance error', emax > 0 ? `${fmtNum((bal / emax) * 100)}% of peak energy` : '–'],
276- ['Largest spring strain', `${fmtNum(r.maxStretch * 100)}% (pre-stretch ${fmtNum(d.prestretch * 100)}%)`],
277- ['Peak pickup acceleration', `${fmtNum(peak)} m/s²`],
278- ]
279- .map(([k, v]) => `<div><span>${k}</span><b>${v}</b></div>`)
280- .join('');
281- if (audioUrl) URL.revokeObjectURL(audioUrl);
282- audioUrl = URL.createObjectURL(encodeWav(r.audio, r.fs));
283- q<HTMLAudioElement>('audio').src = audioUrl;
284- q<HTMLAnchorElement>('download').href = audioUrl;
285- q('audio-note').textContent = `16-bit mono, ${r.fs} Hz, normalized so the peak is −1 dBFS.`;
286- drawWave();
287- spec.setSignal(r.audio, r.fs);
288- drawSpectrum();
289- timeline.setTimes(r.frameTimes);
290- scene.showReplay({ ...sceneSpec(p, sys, exc), frames: r.frames });
291- const rest = new Float64Array(p.n);
292- for (let i = 0; i < p.n; i++) rest[i] = sys.rest[3 * i];
293- viewer.set(
294- {
295- n: p.n,
296- restX: rest,
297- leftWall: p.leftEnd === 'free' ? null : 0,
298- rightWall: p.rightEnd === 'free' ? null : (p.n + 1) * p.spacing,
299- slideLeft: p.leftEnd === 'sliding',
300- slideRight: p.rightEnd === 'sliding',
301- pinned: sys.pinned,
302- },
303- { frames: r.frames, times: r.frameTimes },
304- );
305- energyPlot.set([{ x: r.frameTimes, y: r.energy, color: css('--series-1'), label: 'energy (J)' }], {
306- xLabel: 'time (s)',
307- yLabel: 'energy (J)',
308- yRange: [0, emax * 1.05 || 1],
309- });
310- }
311-
312- runBtn.addEventListener('click', async () => {
313- const p = controls.get();
90+ const build = (p: StringParams): Built | null => {
31491 const d = derive(p);
92+ if (!(d.restLength > 0) || !(p.n >= 1)) return null;
31593 const sys = buildSystem(p, d);
31694 const exc = buildExcitation(p);
317- const opt = { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: buildPickup(p) };
318- runBtn.disabled = true;
319- cancelBtn.hidden = false;
320- progress.hidden = false;
321- progress.value = 0;
322- running = runInWorker(sys, exc, opt, (f) => (progress.value = f));
323- try {
324- show(p, sys, exc, await running.promise);
325- } catch (e) {
326- if ((e as Error).message !== 'cancelled') alert(`Simulation failed: ${(e as Error).message}`);
327- } finally {
328- running = null;
329- runBtn.disabled = false;
330- cancelBtn.hidden = true;
331- progress.hidden = true;
332- }
95+ const pickup = buildPickup(p);
96+ return {
97+ sys,
98+ exc,
99+ pickup,
100+ scene: {
101+ sys,
102+ exc,
103+ excKind: p.excitation,
104+ pickup,
105+ pickupLabel: `${p.pickupMode === 'sum' ? 'sum over all beads' : `bead ${Math.round(p.pickupBead)}`}, ${p.pickupAxis} acceleration`,
106+ },
107+ };
108+ };
109+
110+ mountInstrumentPage<StringParams>(root, {
111+ title: 'String of beads',
112+ wavName: 'bead-string.wav',
113+ presets: PRESETS,
114+ sections: SECTIONS,
115+ setup(p) {
116+ const d = derive(p);
117+ const tr = d.transverseHz;
118+ const rows: [string, string][] = [
119+ ['Length (wall to wall)', `${fmtNum(d.length)} m`],
120+ ['Rest length L0', `${fmtNum(d.restLength * 1000)} mm (pre-stretch ${fmtNum(d.prestretch * 100)}%)`],
121+ ['Linear density', `${fmtNum(d.linearDensity * 1000)} g/m`],
122+ ['Wave speed (transverse, longitudinal)', `${fmtNum(d.waveSpeedT)}, ${fmtNum(d.waveSpeedL)} m/s`],
123+ ['Lowest transverse mode', tr.length ? `${fmtNum(tr[0])} Hz` : 'none'],
124+ ['Overtone ratios', tr.slice(1, 7).map((f) => fmtNum(f / tr[0])).join(', ') || '–'],
125+ ['Lowest longitudinal mode', d.longitudinalHz.length ? `${fmtNum(d.longitudinalHz[0])} Hz` : 'none'],
126+ ];
127+ const warnings = [...d.warnings];
128+ const pins = parsePins(p.pins, p.n);
129+ if (pins.includes(Math.round(p.excBead) - 1)) warnings.push('The excited bead is pinned.');
130+ if (p.pickupMode === 'bead' && pins.includes(Math.round(p.pickupBead) - 1)) warnings.push('The pickup bead is pinned: its acceleration is zero.');
131+ if (tr.length && tr[0] < 20) warnings.push(`The lowest mode (${fmtNum(tr[0])} Hz) is below hearing; you will mostly hear overtones.`);
132+ return { rows, warnings, built: build(p) };
133+ },
134+ summaryRows: (p, _b, r) => [['Largest spring strain', `${fmtNum(r.maxStretch * 100)}% (pre-stretch ${fmtNum(derive(p).prestretch * 100)}%)`]],
135+ spectrumMarkers(p) {
136+ const d = derive(p);
137+ return {
138+ markers: [
139+ ...d.transverseHz.map((f, i) => ({ x: f, color: css('--series-1'), label: `transverse mode ${i + 1}: ${fmtNum(f)} Hz` })),
140+ ...d.longitudinalHz.map((f, i) => ({ x: f, color: css('--series-2'), label: `longitudinal mode ${i + 1}: ${fmtNum(f)} Hz` })),
141+ ],
142+ legend:
143+ '<span><i class="dash" style="border-color:var(--series-1)"></i>linear transverse modes</span>' +
144+ '<span><i class="dash" style="border-color:var(--series-2)"></i>linear longitudinal modes</span>',
145+ };
146+ },
147+ extraCards: '<section class="card" data-r="out" hidden><h2>Displacements (2D)</h2><div data-r="viewer"></div></section>',
148+ mountExtras(q, timeline) {
149+ const viewer = new ChainView(q('viewer'), timeline);
150+ return {
151+ show(p, built, r) {
152+ const { sys } = built;
153+ const rest = new Float64Array(p.n);
154+ for (let i = 0; i < p.n; i++) rest[i] = sys.rest[3 * i];
155+ viewer.set(
156+ {
157+ n: p.n,
158+ restX: rest,
159+ leftWall: p.leftEnd === 'free' ? null : 0,
160+ rightWall: p.rightEnd === 'free' ? null : (p.n + 1) * p.spacing,
161+ slideLeft: p.leftEnd === 'sliding',
162+ slideRight: p.rightEnd === 'sliding',
163+ pinned: sys.pinned,
164+ },
165+ { frames: r.frames, times: r.frameTimes },
166+ );
167+ },
168+ };
169+ },
333170 });
334- cancelBtn.addEventListener('click', () => running?.cancel());
335171 }
src/instruments/registry.tsmodified+7−1View file
@@ -1,4 +1,5 @@
11 import { mountBeadString } from './bead-string/page';
2+import { mountBeadDrum } from './bead-drum/page';
23
34 export interface Instrument {
45 id: string;
@@ -14,7 +15,12 @@ export const INSTRUMENTS: Instrument[] = [
1415 blurb: 'Point masses joined by pre-tensioned springs, moving in 3D. Pluck, strike or hammer it; clamp, slide or free the ends; pin beads.',
1516 mount: mountBeadString,
1617 },
17- { id: 'bead-drum', name: 'Bead drum', blurb: 'A membrane: a tensioned 2D net of beads with a clamped rim.' },
18+ {
19+ id: 'bead-drum',
20+ name: 'Bead drum',
21+ blurb: 'A round drumhead: a tensioned hex net of beads in a clamped rim. Hit it with a felt mallet, strike, push or press it.',
22+ mount: mountBeadDrum,
23+ },
1824 { id: 'bead-bar', name: 'Bar', blurb: 'A stiff chain with bending springs, free ends: xylophone-like overtones.' },
1925 { id: 'bead-plate', name: 'Plate', blurb: 'A 2D net with bending, free or clamped edge.' },
2026 ];
src/sim/solver.tsmodified+71−7View file
@@ -1,9 +1,9 @@
1-import type { BeadSystem, Excitation, SimOptions, SimResult } from './types';
1+import type { BeadSystem, Excitation, Mallet, SimOptions, SimResult } from './types';
22
33 // Upper bound on the highest angular frequency of the linearized system
44 // (Gershgorin on M^-1 K at rest). A spring under tension T has stiffness k along
55 // its axis and T/l across it, so each spring contributes 2 max(k, |T|/l) / m.
6-export function omegaMaxBound(sys: BeadSystem): number {
6+export function omegaMaxBound(sys: BeadSystem, mallet?: Mallet): number {
77 const { n, si, sj, sk, sL0 } = sys;
88 const rowSum = new Float64Array(n);
99 const p = (node: number, c: number) =>
@@ -19,8 +19,20 @@ export function omegaMaxBound(sys: BeadSystem): number {
1919 if (a < n) rowSum[a] += 2 * kmax;
2020 if (b < n) rowSum[b] += 2 * kmax;
2121 }
22- let w2 = 0;
23- for (let i = 0; i < n; i++) if (!sys.pinned[i]) w2 = Math.max(w2, rowSum[i] / sys.mass[i]);
22+ // a mallet adds a rank-one contact stiffness k w w^T between it and the beads,
23+ // linearized at the largest compression a rigid surface would see
24+ let kc = 0;
25+ if (mallet) {
26+ const { K, p, mass, v0 } = mallet;
27+ const dmax = (((p + 1) * mass * v0 * v0) / (2 * K)) ** (1 / (p + 1));
28+ kc = p * K * dmax ** (p - 1);
29+ }
30+ let w2 = mallet ? kc / mallet.mass : 0;
31+ for (let i = 0; i < n; i++) {
32+ if (sys.pinned[i]) continue;
33+ const extra = mallet ? kc * mallet.weights[i] : 0;
34+ w2 = Math.max(w2, (rowSum[i] + extra) / sys.mass[i]);
35+ }
2436 return Math.sqrt(w2);
2537 }
2638
@@ -40,7 +52,8 @@ export function simulate(
4052 const nn = n + anchors.length;
4153 const ns = si.length;
4254
43- const omegaMax = omegaMaxBound(sys);
55+ const mallet = exc.mallet;
56+ const omegaMax = omegaMaxBound(sys, mallet);
4457 const dtMax = (opt.cfl * 2) / omegaMax;
4558 const stepsPerSample = Math.max(1, Math.ceil(1 / (opt.fs * dtMax)));
4659 const dt = 1 / (opt.fs * stepsPerSample);
@@ -83,12 +96,48 @@ export function simulate(
8396 ? force.F0 * 0.5 * (1 - Math.cos((2 * Math.PI * t) / force.tau))
8497 : 0;
8598
99+ // mallet state: position s and velocity along dir; contact force and its potential
100+ let ms = 0;
101+ let mv = mallet ? mallet.v0 : 0;
102+ let ma = 0;
103+ let fContact = 0;
104+ let eContact = 0;
105+ const mw = mallet?.weights;
106+ const [mdx, mdy, mdz] = mallet?.dir ?? [0, 0, 0];
107+ const contact = () => {
108+ fContact = 0;
109+ eContact = 0;
110+ if (!mallet) return;
111+ let h = 0;
112+ for (let i = 0; i < n; i++) {
113+ const wi = mw![i];
114+ if (wi === 0) continue;
115+ const k = 3 * i;
116+ h += wi * ((X[k] - sys.rest[k]) * mdx + (X[k + 1] - sys.rest[k + 1]) * mdy + (X[k + 2] - sys.rest[k + 2]) * mdz);
117+ }
118+ const d = ms - h;
119+ if (d > 0) {
120+ const dp = mallet.K * d ** mallet.p;
121+ fContact = dp;
122+ eContact = (dp * d) / (mallet.p + 1);
123+ for (let i = 0; i < n; i++) {
124+ const wi = mw![i];
125+ if (wi === 0) continue;
126+ F[3 * i] += fContact * wi * mdx;
127+ F[3 * i + 1] += fContact * wi * mdy;
128+ F[3 * i + 2] += fContact * wi * mdz;
129+ }
130+ }
131+ ma = -fContact / mallet.mass;
132+ };
133+
86134 // Fills F with total force on beads; returns damping power (W, >= 0)
87135 // and external power (W), both at velocity V.
88136 let pDamp = 0;
89137 let pExt = 0;
90138 const computeForces = (t: number) => {
91139 F.fill(0);
140+ contact();
92141 pDamp = 0;
93142 for (let s = 0; s < ns; s++) {
94143 const a = 3 * si[s];
@@ -170,6 +219,7 @@ export function simulate(
170219 e += 0.5 * sk[s] * d * d;
171220 maxStretch = Math.max(maxStretch, Math.abs(d) / sL0[s]);
172221 }
222+ if (mallet) e += 0.5 * mallet.mass * mv * mv + eContact;
173223 return e;
174224 };
175225
@@ -178,6 +228,8 @@ export function simulate(
178228 const frameTimes = new Float64Array(nFrames);
179229 const energy = new Float64Array(nFrames);
180230 const balance = new Float64Array(nFrames);
231+ const contactForce = mallet ? new Float32Array(nSamples) : undefined;
232+ const malletPos = mallet ? new Float32Array(nFrames) : undefined;
181233 const w = opt.pickup;
182234
183235 // energy is reported relative to the equilibrium state (pre-tension energy subtracted)
@@ -187,6 +239,7 @@ export function simulate(
187239 const l = Math.hypot(pa(sj[s], 0) - pa(si[s], 0), pa(sj[s], 1) - pa(si[s], 1), pa(sj[s], 2) - pa(si[s], 2));
188240 Eeq += 0.5 * sk[s] * (l - sL0[s]) ** 2;
189241 }
242+ computeForces(0); // contact energy at t = 0
190243 const E0 = energyNow();
191244 let dissipated = 0;
192245 let workIn = 0;
@@ -200,12 +253,16 @@ export function simulate(
200253 let acc = 0;
201254 for (let k = 0; k < 3 * n; k++) acc += w[k] * A[k];
202255 audio[smp] = acc;
256+ if (contactForce) contactForce[smp] = fContact;
203257 if (smp % samplesPerFrame === 0 && frame < nFrames) {
258+ if (malletPos) malletPos[frame] = ms;
204259 const off = frame * 3 * n;
205260 for (let k = 0; k < 3 * n; k++) frames[off + k] = X[k] - sys.rest[k];
206261 frameTimes[frame] = t;
207- energy[frame] = energyNow() - Eeq;
208- balance[frame] = energy[frame] + Eeq + dissipated - workIn - E0;
262+ const eTotal = energyNow();
263+ // the plot shows the instrument's own energy; the balance includes the mallet
264+ energy[frame] = eTotal - Eeq - (mallet ? 0.5 * mallet.mass * mv * mv + eContact : 0);
265+ balance[frame] = eTotal + dissipated - workIn - E0;
209266 frame++;
210267 }
211268 if (smp === nSamples - 1) break;
@@ -216,6 +273,10 @@ export function simulate(
216273 V[k] += 0.5 * dt * A[k];
217274 X[k] += dt * V[k];
218275 }
276+ if (mallet) {
277+ mv += 0.5 * dt * ma;
278+ ms += dt * mv;
279+ }
219280 syncAnchors();
220281 t += dt;
221282 // forces at new position with half-step velocity; powers use it too
@@ -223,6 +284,7 @@ export function simulate(
223284 dissipated += dt * pDamp;
224285 workIn += dt * pExt;
225286 for (let k = 0; k < 3 * n; k++) V[k] += 0.5 * dt * A[k];
287+ if (mallet) mv += 0.5 * dt * ma;
226288 syncAnchors();
227289 }
228290 // at a sample time: recompute with the full-step velocity
@@ -242,5 +304,7 @@ export function simulate(
242304 omegaMax,
243305 wallSeconds: (performance.now() - t0) / 1000,
244306 maxStretch,
307+ contactForce,
308+ malletPos,
245309 };
246310 }
src/sim/types.tsmodified+15−0View file
@@ -28,6 +28,19 @@ export interface Excitation {
2828 v0: Float64Array; // 3n initial velocity
2929 // External force F(t) = shape * F0 * 0.5 (1 - cos(2 pi t / tau)), 0 <= t <= tau
3030 force?: { shape: Float64Array; F0: number; tau: number };
31+ // A mallet: point mass moving along dir, touching the beads through a felt
32+ // spring with force K * delta^p, delta = s - sum_i w_i (x_i - rest_i) . dir.
33+ // Starts at s = 0 (just touching the undeformed surface) with speed v0 along dir.
34+ mallet?: Mallet;
35+}
36+
37+export interface Mallet {
38+ weights: Float64Array; // n, sums to 1 (contact patch)
39+ dir: [number, number, number]; // unit
40+ mass: number; // kg
41+ K: number; // N / m^p
42+ p: number;
43+ v0: number; // m/s, > 0 toward the surface
3144 }
3245
3346 export interface SimOptions {
@@ -50,4 +63,6 @@ export interface SimResult {
5063 omegaMax: number;
5164 wallSeconds: number;
5265 maxStretch: number; // largest |l - L0| / L0 seen, at frame times
66+ contactForce?: Float32Array; // mallet contact force at audio rate (N)
67+ malletPos?: Float32Array; // mallet position s along dir at frame times (m)
5368 }
src/style.cssmodified+5−0View file
@@ -164,3 +164,8 @@ audio { max-width: 100%; }
164164 .seg button:disabled { opacity: 0.45; cursor: default; }
165165 .btn-small.icon { padding: 4px 8px; }
166166 .viewer-bar[hidden] { display: none; }
167+.membrane-row { display: grid; grid-template-columns: minmax(0, 1fr) minmax(0, 1.4fr); gap: 12px; align-items: center; }
168+@media (max-width: 1100px) { .membrane-row { grid-template-columns: 1fr; } }
169+.map-wrap { position: relative; aspect-ratio: 1; max-width: 100%; max-height: 420px; }
170+.map-wrap canvas { width: 100%; height: 100%; display: block; }
171+.profile-wrap { height: 260px; }
src/ui/instrumentPage.tsadded+255−0View file
@@ -0,0 +1,255 @@
1+import { renderControls, type Section } from './controls';
2+import { LinePlot, css, fmtNum, type Marker } from './plot';
3+import { Spectrogram } from './spectrogram';
4+import { spectrum } from './fft';
5+import { encodeWav } from './wav';
6+import { BeadScene, type SceneSpec } from './scene3d';
7+import { Timeline } from './timeline';
8+import { omegaMaxBound } from '../sim/solver';
9+import { runInWorker, type Running } from '../sim/run';
10+import type { BeadSystem, Excitation, SimResult } from '../sim/types';
11+
12+// The page every instrument shares: presets and control panel, 3D model,
13+// setup summary, run summary, audio, waveform, spectrogram, spectrum, energy.
14+// An instrument adds its own setup rows, extra result cards and views.
15+
16+export interface SimParams {
17+ duration: number;
18+ fs: number;
19+ cfl: number;
20+ frameRate: number;
21+}
22+
23+export interface Built {
24+ sys: BeadSystem;
25+ exc: Excitation;
26+ pickup: Float64Array;
27+ scene: SceneSpec;
28+}
29+
30+export interface Preset<P> {
31+ name: string;
32+ description: string;
33+ params: P;
34+}
35+
36+export interface Extras<P> {
37+ show(p: P, built: Built, r: SimResult): void;
38+}
39+
40+export interface InstrumentConfig<P extends SimParams> {
41+ title: string;
42+ wavName: string;
43+ presets: Preset<P>[];
44+ sections: Section<P>[];
45+ // instrument-specific summary of the setup; built = null when the parameters are invalid
46+ setup(p: P): { rows: [string, string][]; warnings: string[]; built: Built | null };
47+ summaryRows?(p: P, built: Built, r: SimResult): [string, string][];
48+ spectrumMarkers?(p: P): { markers: Marker[]; legend: string };
49+ extraCards?: string; // cards with data-r="out" hidden, shown after a run
50+ mountExtras?(q: <T extends HTMLElement>(r: string) => T, timeline: Timeline): Extras<P>;
51+}
52+
53+const statRows = (rows: [string, string][]) => rows.map(([k, v]) => `<div><span>${k}</span><b>${v}</b></div>`).join('');
54+
55+export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg: InstrumentConfig<P>) {
56+ root.innerHTML = `
57+ <div class="instrument">
58+ <aside class="panel">
59+ <a href="#/" class="back">← All instruments</a>
60+ <h1>${cfg.title}</h1>
61+ <label class="ctl preset"><span class="ctl-name">Preset</span><select data-r="preset"></select></label>
62+ <p class="sub" data-r="preset-desc"></p>
63+ <div data-r="controls"></div>
64+ <div class="run-row">
65+ <button class="btn" data-r="run">Run simulation</button>
66+ <button class="btn-small" data-r="cancel" hidden>Cancel</button>
67+ </div>
68+ <progress data-r="progress" max="1" value="0" hidden></progress>
69+ <p class="error" data-r="error" hidden></p>
70+ </aside>
71+ <main class="results">
72+ <section class="card"><h2>Model</h2><div data-r="scene"></div></section>
73+ <section class="card"><h2>Setup</h2><div data-r="derived" class="stats"></div><ul data-r="warnings" class="warnings"></ul></section>
74+ <section class="card" data-r="out" hidden>
75+ <h2>Run</h2><div data-r="summary" class="stats"></div>
76+ <div class="audio-row"><audio controls data-r="audio"></audio><a class="btn-small" data-r="download" download="${cfg.wavName}">Download .wav</a></div>
77+ <p class="note" data-r="audio-note"></p>
78+ </section>
79+ <section class="card" data-r="out" hidden>
80+ <div class="card-head"><h2>Pickup acceleration</h2>
81+ <div class="seg" data-r="wave-range"><button data-v="0" class="on">Full</button><button data-v="0.05">First 50 ms</button><button data-v="0.005">First 5 ms</button></div></div>
82+ <div class="canvas-wrap"><canvas data-r="wave"></canvas></div>
83+ </section>
84+ <section class="card" data-r="out" hidden>
85+ <div class="card-head"><h2>Spectrogram</h2>
86+ <label class="inline">up to <select data-r="fmax"><option>1000</option><option>2500</option><option selected>5000</option><option>10000</option><option>24000</option></select> Hz</label></div>
87+ <div class="canvas-wrap tall"><canvas data-r="spec"></canvas></div>
88+ </section>
89+ <section class="card" data-r="out" hidden>
90+ <div class="card-head"><h2>Spectrum</h2><div class="legend" data-r="spec-legend"></div></div>
91+ <div class="canvas-wrap"><canvas data-r="spectrum"></canvas></div>
92+ </section>
93+ ${cfg.extraCards ?? ''}
94+ <section class="card" data-r="out" hidden>
95+ <h2>Energy</h2>
96+ <div class="canvas-wrap short"><canvas data-r="energy"></canvas></div>
97+ <p class="note">Kinetic plus spring energy of the instrument, minus the energy stored by pre-tension at equilibrium. Without damping or forcing it stays flat; the run summary gives the error of the energy balance (energy + dissipated − work in − initial, mallet included).</p>
98+ </section>
99+ </main>
100+ </div>`;
101+ const q = <T extends HTMLElement>(r: string) => root.querySelector(`[data-r="${r}"]`) as T;
102+ const presetSel = q<HTMLSelectElement>('preset');
103+ cfg.presets.forEach((p, i) => presetSel.add(new Option(p.name, String(i))));
104+ const presetDesc = q<HTMLParagraphElement>('preset-desc');
105+ const runBtn = q<HTMLButtonElement>('run');
106+ const cancelBtn = q<HTMLButtonElement>('cancel');
107+ const progress = q<HTMLProgressElement>('progress');
108+ const errorBox = q<HTMLParagraphElement>('error');
109+
110+ const wave = new LinePlot(q('wave'));
111+ const spec = new Spectrogram(q('spec'));
112+ const specPlot = new LinePlot(q('spectrum'));
113+ const energyPlot = new LinePlot(q('energy'));
114+ const timeline = new Timeline();
115+ const scene = new BeadScene(q('scene'), timeline);
116+ const extras = cfg.mountExtras?.(q, timeline);
117+ let waveRange = 0;
118+ let last: { p: P; r: SimResult } | null = null;
119+ let running: Running | null = null;
120+ let audioUrl: string | null = null;
121+
122+ const controls = renderControls(q('controls'), cfg.sections, cfg.presets[0].params, () => updateSetup());
123+ presetDesc.textContent = cfg.presets[0].description;
124+ presetSel.addEventListener('change', () => {
125+ const pr = cfg.presets[Number(presetSel.value)];
126+ presetDesc.textContent = pr.description;
127+ controls.set(pr.params);
128+ });
129+
130+ function updateSetup() {
131+ const p = controls.get();
132+ const { rows, warnings, built } = cfg.setup(p);
133+ if (built) {
134+ scene.showSetup(built.scene);
135+ const wmax = omegaMaxBound(built.sys, built.exc.mallet);
136+ const spp = Math.max(1, Math.ceil(1 / (p.fs * ((p.cfl * 2) / wmax))));
137+ rows.push(['Highest frequency bound', `${fmtNum(wmax / (2 * Math.PI))} Hz`]);
138+ rows.push(['Time step', `${fmtNum(1 / (p.fs * spp))} s (${spp} per sample)`]);
139+ const work = spp * p.fs * p.duration * built.sys.si.length;
140+ rows.push(['Spring updates', fmtNum(work)]);
141+ const frameMB = (Math.round(p.duration * p.frameRate) * 3 * built.sys.n * 4) / 1e6;
142+ if (frameMB > 400) warnings.push(`Replay frames would take ${fmtNum(frameMB)} MB; lower the replay frame rate or the duration.`);
143+ if (work > 3e9) warnings.push('This run is large (more than 3e9 spring updates) and may take minutes.');
144+ }
145+ q('derived').innerHTML = statRows(rows);
146+ q('warnings').innerHTML = warnings.map((w) => `<li>${w}</li>`).join('');
147+ runBtn.disabled = !built || !!running;
148+ }
149+ updateSetup();
150+
151+ q('wave-range').addEventListener('click', (e) => {
152+ const b = (e.target as HTMLElement).closest('button');
153+ if (!b) return;
154+ waveRange = Number(b.dataset.v);
155+ q('wave-range').querySelectorAll('button').forEach((x) => x.classList.toggle('on', x === b));
156+ drawWave();
157+ });
158+ q<HTMLSelectElement>('fmax').addEventListener('change', (e) => {
159+ spec.fmax = Number((e.target as HTMLSelectElement).value);
160+ spec.draw();
161+ drawSpectrum();
162+ });
163+
164+ function drawWave() {
165+ if (!last) return;
166+ const { r } = last;
167+ const t = new Float64Array(r.audio.length);
168+ for (let i = 0; i < t.length; i++) t[i] = i / r.fs;
169+ wave.set([{ x: t, y: r.audio, color: css('--series-1'), label: 'a' }], {
170+ xLabel: 'time (s)',
171+ yLabel: 'acceleration (m/s²)',
172+ xRange: waveRange ? [0, Math.min(waveRange, t[t.length - 1])] : undefined,
173+ });
174+ }
175+
176+ function drawSpectrum() {
177+ if (!last) return;
178+ const { p, r } = last;
179+ const s = spectrum(r.audio, r.fs);
180+ let top = -Infinity;
181+ for (const v of s.db) top = Math.max(top, v);
182+ const db = s.db.map((v) => Math.max(v - top, -100));
183+ const m = cfg.spectrumMarkers?.(p);
184+ q('spec-legend').innerHTML =
185+ `<span><i class="sw" style="background:var(--series-1)"></i>pickup signal</span>` + (m?.legend ?? '');
186+ specPlot.set([{ x: s.freq, y: db, color: css('--series-1'), label: 'level (dB)' }], {
187+ xLabel: 'frequency (Hz)',
188+ yLabel: 'level (dB re max)',
189+ xRange: [0, Math.min(spec.fmax, r.fs / 2)],
190+ yRange: [-100, 5],
191+ markers: m?.markers,
192+ });
193+ }
194+
195+ function show(p: P, built: Built, r: SimResult) {
196+ last = { p, r };
197+ root.querySelectorAll<HTMLElement>('[data-r="out"]').forEach((e) => (e.hidden = false));
198+ let peak = 0;
199+ for (const v of r.audio) peak = Math.max(peak, Math.abs(v));
200+ const emax = Math.max(...r.energy);
201+ let bal = 0;
202+ for (const b of r.balance) bal = Math.max(bal, Math.abs(b));
203+ const rt = p.duration / r.wallSeconds;
204+ q('summary').innerHTML = statRows([
205+ ['Time step', `${fmtNum(r.dt)} s, ${r.stepsPerSample} per sample`],
206+ ['Wall time', `${fmtNum(r.wallSeconds)} s (${rt >= 1 ? fmtNum(rt) + '× faster than' : fmtNum(1 / rt) + '× slower than'} real time)`],
207+ ['Energy balance error', emax > 0 ? `${fmtNum((bal / emax) * 100)}% of peak energy` : '–'],
208+ ...(cfg.summaryRows?.(p, built, r) ?? []),
209+ ['Peak pickup acceleration', `${fmtNum(peak)} m/s²`],
210+ ]);
211+ if (audioUrl) URL.revokeObjectURL(audioUrl);
212+ audioUrl = URL.createObjectURL(encodeWav(r.audio, r.fs));
213+ q<HTMLAudioElement>('audio').src = audioUrl;
214+ q<HTMLAnchorElement>('download').href = audioUrl;
215+ q('audio-note').textContent = `16-bit mono, ${r.fs} Hz, normalized so the peak is −1 dBFS.`;
216+ drawWave();
217+ spec.setSignal(r.audio, r.fs);
218+ drawSpectrum();
219+ timeline.setTimes(r.frameTimes);
220+ scene.showReplay({ ...built.scene, frames: r.frames, malletPos: r.malletPos });
221+ extras?.show(p, built, r);
222+ energyPlot.set([{ x: r.frameTimes, y: r.energy, color: css('--series-1'), label: 'energy (J)' }], {
223+ xLabel: 'time (s)',
224+ yLabel: 'energy (J)',
225+ yRange: [0, emax * 1.05 || 1],
226+ });
227+ }
228+
229+ runBtn.addEventListener('click', async () => {
230+ const p = controls.get();
231+ const { built } = cfg.setup(p);
232+ if (!built) return;
233+ const opt = { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: built.pickup };
234+ runBtn.disabled = true;
235+ errorBox.hidden = true;
236+ cancelBtn.hidden = false;
237+ progress.hidden = false;
238+ progress.value = 0;
239+ running = runInWorker(built.sys, built.exc, opt, (f) => (progress.value = f));
240+ try {
241+ show(p, built, await running.promise);
242+ } catch (e) {
243+ if ((e as Error).message !== 'cancelled') {
244+ errorBox.textContent = `The simulation failed: ${(e as Error).message}`;
245+ errorBox.hidden = false;
246+ }
247+ } finally {
248+ running = null;
249+ runBtn.disabled = false;
250+ cancelBtn.hidden = true;
251+ progress.hidden = true;
252+ }
253+ });
254+ cancelBtn.addEventListener('click', () => running?.cancel());
255+}
src/ui/membraneview.tsadded+149−0View file
@@ -0,0 +1,149 @@
1+import { LinePlot, css, fmtNum, setupCanvas } from './plot';
2+import { rampLUT } from './spectrogram';
3+import { maxAbsPerComponent } from './chainview';
4+import { timelineBar, type Timeline } from './timeline';
5+
6+// 2D views of a flat membrane in the x-z plane vibrating along y: a top-down
7+// map of y displacement, and a cross-section along a chosen diameter.
8+
9+export interface MembraneGeometry {
10+ pos: [number, number][]; // bead (x, z)
11+ radius: number; // rim radius
12+ spacing: number; // lattice spacing
13+ strike: [number, number] | null; // excitation point
14+ pickup: number | null; // pickup bead, null = sum
15+ sectionAngle: number; // degrees, direction of the cross-section diameter
16+}
17+
18+export class MembraneView {
19+ private g: MembraneGeometry | null = null;
20+ private frames: Float32Array | null = null;
21+ private maxY = 1;
22+ private section: number[] = []; // bead indices along the diameter, sorted
23+ private sectionS: Float64Array = new Float64Array(0);
24+ private map: HTMLCanvasElement;
25+ private profile: LinePlot;
26+ private note: HTMLParagraphElement;
27+
28+ constructor(root: HTMLElement, private tl: Timeline) {
29+ root.innerHTML = `
30+ <div data-r="bar"></div>
31+ <div class="membrane-row">
32+ <div class="map-wrap"><canvas data-r="map"></canvas></div>
33+ <div class="canvas-wrap profile-wrap"><canvas data-r="profile"></canvas></div>
34+ </div>
35+ <p class="note" data-r="note"></p>`;
36+ const q = <T extends HTMLElement>(r: string) => root.querySelector(`[data-r="${r}"]`) as T;
37+ q('bar').replaceWith(timelineBar(tl));
38+ this.map = q('map');
39+ this.profile = new LinePlot(q('profile'));
40+ this.note = q('note');
41+ tl.onChange(() => this.draw());
42+ new ResizeObserver(() => this.draw()).observe(this.map);
43+ }
44+
45+ set(g: MembraneGeometry, frames: Float32Array) {
46+ this.g = g;
47+ this.frames = frames;
48+ this.maxY = Math.max(maxAbsPerComponent(frames)[1], 1e-30);
49+ const th = (g.sectionAngle * Math.PI) / 180;
50+ const ux = Math.cos(th);
51+ const uz = Math.sin(th);
52+ const near = g.pos
53+ .map(([x, z], i) => ({ i, s: x * ux + z * uz, d: Math.abs(-x * uz + z * ux) }))
54+ .filter((b) => b.d < 0.51 * g.spacing)
55+ .sort((a, b) => a.s - b.s);
56+ this.section = near.map((b) => b.i);
57+ this.sectionS = Float64Array.from(near.map((b) => b.s));
58+ this.draw();
59+ }
60+
61+ private draw() {
62+ const g = this.g;
63+ const f = this.frames;
64+ if (!g || !f) return;
65+ const n = g.pos.length;
66+ const off = this.tl.frame * 3 * n;
67+
68+ // top-down map
69+ const { ctx, w, h } = setupCanvas(this.map);
70+ ctx.clearRect(0, 0, w, h);
71+ const size = Math.min(w, h);
72+ const s = (0.46 * size) / g.radius;
73+ const cx = w / 2;
74+ const cy = h / 2;
75+ const lut = rampLUT('--diverging');
76+ const cell = 0.58 * g.spacing * s; // hexagon circumradius that tiles the lattice
77+ for (let i = 0; i < n; i++) {
78+ const v = f[off + 3 * i + 1] / this.maxY;
79+ const li = Math.round((0.5 + 0.5 * Math.max(-1, Math.min(1, v))) * 255) * 3;
80+ ctx.fillStyle = `rgb(${lut[li]},${lut[li + 1]},${lut[li + 2]})`;
81+ const [x, z] = g.pos[i];
82+ const px = cx + x * s;
83+ const pz = cy + z * s;
84+ ctx.beginPath();
85+ for (let k = 0; k < 6; k++) {
86+ const a = (Math.PI / 3) * k + Math.PI / 6;
87+ const X = px + cell * Math.cos(a);
88+ const Z = pz + cell * Math.sin(a);
89+ if (k) ctx.lineTo(X, Z);
90+ else ctx.moveTo(X, Z);
91+ }
92+ ctx.closePath();
93+ ctx.fill();
94+ }
95+ ctx.strokeStyle = css('--text-muted');
96+ ctx.lineWidth = 2;
97+ ctx.beginPath();
98+ ctx.arc(cx, cy, g.radius * s, 0, 2 * Math.PI);
99+ ctx.stroke();
100+ // section line
101+ const th = (g.sectionAngle * Math.PI) / 180;
102+ ctx.strokeStyle = css('--text-secondary');
103+ ctx.setLineDash([4, 4]);
104+ ctx.lineWidth = 1;
105+ ctx.beginPath();
106+ ctx.moveTo(cx - g.radius * s * Math.cos(th), cy - g.radius * s * Math.sin(th));
107+ ctx.lineTo(cx + g.radius * s * Math.cos(th), cy + g.radius * s * Math.sin(th));
108+ ctx.stroke();
109+ ctx.setLineDash([]);
110+ if (g.strike) {
111+ const [x, z] = g.strike;
112+ ctx.strokeStyle = css('--series-2');
113+ ctx.lineWidth = 2.5;
114+ const r = 6;
115+ ctx.beginPath();
116+ ctx.moveTo(cx + x * s - r, cy + z * s - r);
117+ ctx.lineTo(cx + x * s + r, cy + z * s + r);
118+ ctx.moveTo(cx + x * s + r, cy + z * s - r);
119+ ctx.lineTo(cx + x * s - r, cy + z * s + r);
120+ ctx.stroke();
121+ }
122+ if (g.pickup !== null) {
123+ const [x, z] = g.pos[g.pickup];
124+ ctx.strokeStyle = css('--series-3');
125+ ctx.lineWidth = 2.5;
126+ ctx.beginPath();
127+ ctx.arc(cx + x * s, cy + z * s, 7, 0, 2 * Math.PI);
128+ ctx.stroke();
129+ }
130+ ctx.fillStyle = css('--text-muted');
131+ ctx.font = '11px system-ui, sans-serif';
132+ ctx.fillText('x →', w - 34, h - 8);
133+ ctx.fillText('z ↓', 6, 14);
134+
135+ // cross-section
136+ const y = new Float64Array(this.section.length);
137+ this.section.forEach((i, k) => (y[k] = f[off + 3 * i + 1] * 1000));
138+ const m = this.maxY * 1000 * 1.1;
139+ this.profile.set([{ x: this.sectionS, y, color: css('--series-1'), label: 'y (mm)' }], {
140+ xLabel: `position along the dashed diameter (m)`,
141+ yLabel: 'y displacement (mm)',
142+ xRange: [-g.radius, g.radius],
143+ yRange: [-m, m],
144+ });
145+ this.note.textContent =
146+ `Colour: vertical (y) displacement, blue down and red up, scaled to ±${fmtNum(this.maxY * 1000)} mm (the largest over the run). ` +
147+ `${g.strike ? 'Orange ×: excitation point. ' : ''}${g.pickup !== null ? 'Green ring: pickup bead.' : 'Pickup: sum over all beads.'}`;
148+ }
149+}
src/ui/scene3d.tsmodified+68−7View file
@@ -11,13 +11,15 @@ import { timelineBar, type Timeline } from './timeline';
1111 export interface SceneSpec {
1212 sys: BeadSystem;
1313 exc: Excitation;
14- excKind: 'pluck' | 'strike' | 'hammer';
14+ excKind: 'pluck' | 'strike' | 'hammer' | 'mallet' | 'force' | 'press';
1515 pickup: Float64Array; // 3n weights
1616 pickupLabel: string;
17+ hoop?: { radius: number }; // round rim in the x-z plane, centred on the origin (instead of walls)
1718 }
1819
1920 export interface ReplaySpec extends SceneSpec {
2021 frames: Float32Array; // nFrames x 3n, matches the timeline
22+ malletPos?: Float32Array; // mallet position along its direction, per frame
2123 }
2224
2325 type Mode = 'setup' | 'replay';
@@ -56,6 +58,7 @@ export class BeadScene {
5658 private slideRings: { mesh: THREE.Mesh; bead: number }[] = [];
5759 private hammerArrows: THREE.Group | null = null;
5860 private pickupArrow: { arrow: THREE.ArrowHelper; bead: number } | null = null;
61+ private mallet: { group: THREE.Group; contact: THREE.Vector3; dir: THREE.Vector3; restS: number } | null = null;
5962 private spec: SceneSpec | null = null;
6063 private frames: Float32Array | null = null;
6164
@@ -167,7 +170,7 @@ export class BeadScene {
167170 };
168171 // back off far enough that the model (about 1.4 wide with its supports) fits horizontally
169172 const half = Math.tan(THREE.MathUtils.degToRad(this.camera.fov / 2));
170- const fit = 0.75 / (half * Math.min(this.camera.aspect, 1)) + 0.3;
173+ const fit = 0.75 / (half * this.camera.aspect) + 0.3;
171174 const dist = Math.max(this.view === 'persp' ? 1.2 : 1.45, fit);
172175 this.camera.position.set(...dirs[v]).setLength(dist);
173176 this.camera.up.set(0, 1, 0);
@@ -211,6 +214,7 @@ export class BeadScene {
211214 this.beads = this.pins = this.springs = this.hammerArrows = null;
212215 this.slideRings = [];
213216 this.pickupArrow = null;
217+ this.mallet = null;
214218 this.spec = spec;
215219 this.frames = this.mode === 'replay' && this.replaySpec ? this.replaySpec.frames : null;
216220 if (!spec) return this.render();
@@ -257,7 +261,10 @@ export class BeadScene {
257261 const beadMat = new THREE.MeshStandardMaterial({ roughness: 0.55, metalness: 0.1 });
258262 this.beads = new THREE.InstancedMesh(new THREE.SphereGeometry(r, 20, 14), beadMat, Math.max(1, n));
259263 this.beads.count = n;
260- for (let i = 0; i < n; i++) this.beads.setColorAt(i, pickBead(i) ? cPick : cBead);
264+ // a pickup over every bead is named in the legend rather than painting the whole model
265+ let nPick = 0;
266+ for (let i = 0; i < n; i++) if (pickBead(i)) nPick++;
267+ for (let i = 0; i < n; i++) this.beads.setColorAt(i, pickBead(i) && nPick < n ? cPick : cBead);
261268 this.model.add(this.beads);
262269 if (pinned.length) {
263270 this.pins = new THREE.InstancedMesh(new THREE.BoxGeometry(2.4 * r, 2.4 * r, 2.4 * r), new THREE.MeshStandardMaterial({ color: cInk, roughness: 0.6 }), pinned.length);
@@ -273,7 +280,16 @@ export class BeadScene {
273280
274281 // supports
275282 const wallSize = 0.16;
276- sys.anchors.forEach((a) => {
283+ if (spec.hoop) {
284+ const hoop = new THREE.Mesh(
285+ new THREE.TorusGeometry(spec.hoop.radius * this.scale, 0.012, 10, 96),
286+ new THREE.MeshStandardMaterial({ color: cMuted, roughness: 0.7 }),
287+ );
288+ hoop.rotation.x = Math.PI / 2;
289+ hoop.position.set(-this.center.x * this.scale, -this.center.y * this.scale, -this.center.z * this.scale);
290+ this.model.add(hoop);
291+ }
292+ if (!spec.hoop) sys.anchors.forEach((a) => {
277293 const at = this.nodePos(sys, sys.n + sys.anchors.indexOf(a), new THREE.Vector3());
278294 const side = at.x < 0 ? -1 : 1;
279295 if (a.slideBead < 0) {
@@ -308,7 +324,8 @@ export class BeadScene {
308324 const sw = (c: THREE.Color, shape = 'dot') => `<i class="glyph ${shape}" style="--c:#${c.getHexString()}"></i>`;
309325 legend.push(`<span>${sw(cBead)}bead</span>`);
310326 if (pinned.length) legend.push(`<span>${sw(cInk, 'square')}pinned bead</span>`);
311- if (sys.anchors.some((a) => a.slideBead < 0)) legend.push(`<span>${sw(cMuted, 'square')}clamped end (wall)</span>`);
327+ if (spec.hoop) legend.push(`<span>${sw(cMuted, 'ring')}clamped rim</span>`);
328+ else if (sys.anchors.some((a) => a.slideBead < 0)) legend.push(`<span>${sw(cMuted, 'square')}clamped end (wall)</span>`);
312329 if (sys.anchors.some((a) => a.slideBead >= 0)) legend.push(`<span>${sw(cInk, 'ring')}sliding end (ring on a frictionless plane)</span>`);
313330 const size0 = modelSize(sys);
314331 const dirs = (field: Float64Array) => {
@@ -346,14 +363,53 @@ export class BeadScene {
346363 if (m > 0) this.model.add(this.arrows(sys, exc.v0, m, cExc, false));
347364 legend.push(`<span>${sw(cExc, 'arrow')}strike: initial velocity (${fmtNum(m)} m/s peak)</span>`);
348365 }
349- if (spec.excKind === 'hammer' && exc.force) {
366+ if (this.mode === 'setup' && spec.excKind === 'press') {
367+ const m = dirs(exc.u0);
368+ if (m > 0) {
369+ const E = (0.5 * GLYPH * size0) / m;
370+ const pts: THREE.Vector3[] = [];
371+ for (let i = 0; i < n; i++) {
372+ if (Math.hypot(exc.u0[3 * i], exc.u0[3 * i + 1], exc.u0[3 * i + 2]) < 0.02 * m) continue;
373+ const w = new THREE.Vector3(exc.u0[3 * i], exc.u0[3 * i + 1], exc.u0[3 * i + 2]).multiplyScalar(E * this.scale);
374+ pts.push(this.nodePos(sys, i, new THREE.Vector3()).add(w));
375+ }
376+ const dots = new THREE.Points(new THREE.BufferGeometry().setFromPoints(pts), new THREE.PointsMaterial({ color: cExc, size: 1.2 * r }));
377+ this.model.add(dots);
378+ legend.push(`<span>${sw(cExc)}press: initial shape, drawn ×${fmtNum(E)} (${fmtNum(m * 1000)} mm deep)</span>`);
379+ }
380+ }
381+ if (exc.mallet) {
382+ const ml = exc.mallet;
383+ const dir = new THREE.Vector3(...ml.dir).normalize();
384+ const contact = new THREE.Vector3();
385+ for (let i = 0; i < n; i++) if (ml.weights[i] > 0) contact.addScaledVector(this.nodePos(sys, i, new THREE.Vector3()), ml.weights[i]);
386+ const head = Math.max(1.6 * r, 0.025);
387+ const group = new THREE.Group();
388+ const felt = new THREE.Mesh(new THREE.SphereGeometry(head, 24, 16), new THREE.MeshStandardMaterial({ color: cExc, roughness: 0.9 }));
389+ felt.position.copy(dir).multiplyScalar(-head);
390+ const stick = new THREE.Mesh(
391+ new THREE.CylinderGeometry(0.15 * head, 0.15 * head, 0.3, 10),
392+ new THREE.MeshStandardMaterial({ color: cMuted, roughness: 0.6 }),
393+ );
394+ stick.position.copy(dir).multiplyScalar(-(2 * head + 0.15));
395+ stick.quaternion.setFromUnitVectors(new THREE.Vector3(0, 1, 0), dir);
396+ group.add(felt, stick);
397+ this.model.add(group);
398+ // in setup the mallet hovers above the head
399+ const restS = -0.6 * GLYPH;
400+ this.mallet = { group, contact, dir, restS };
401+ legend.push(
402+ `<span>${sw(cExc)}mallet: ${fmtNum(ml.mass * 1000)} g at ${fmtNum(ml.v0)} m/s, felt F = K δ^${fmtNum(ml.p)}${this.mode === 'replay' ? ' (position exaggerated with the head)' : ''}</span>`,
403+ );
404+ }
405+ if ((spec.excKind === 'hammer' || spec.excKind === 'force') && exc.force) {
350406 const m = dirs(exc.force.shape);
351407 if (m > 0) {
352408 this.hammerArrows = this.arrows(sys, exc.force.shape, m, cExc, true);
353409 this.model.add(this.hammerArrows);
354410 }
355411 legend.push(
356- `<span>${sw(cExc, 'arrow')}hammer: force pulse, ${fmtNum(exc.force.F0)} N peak over ${fmtNum(exc.force.tau * 1000)} ms${this.mode === 'replay' ? ' (shown while it acts)' : ''}</span>`,
412+ `<span>${sw(cExc, 'arrow')}${spec.excKind === 'hammer' ? 'hammer' : 'force'}: force pulse, ${fmtNum(exc.force.F0)} N peak over ${fmtNum(exc.force.tau * 1000)} ms${this.mode === 'replay' ? ' (shown while it acts)' : ''}</span>`,
357413 );
358414 }
359415
@@ -442,6 +498,11 @@ export class BeadScene {
442498 mesh.position.z = P[bead].z;
443499 }
444500 if (this.pickupArrow) this.pickupArrow.arrow.position.copy(P[this.pickupArrow.bead]);
501+ if (this.mallet) {
502+ const mp = this.mode === 'replay' ? this.replaySpec?.malletPos : undefined;
503+ const sNorm = mp ? mp[this.tl.frame] * E : this.mallet.restS;
504+ this.mallet.group.position.copy(this.mallet.contact).addScaledVector(this.mallet.dir, sNorm);
505+ }
445506 if (this.hammerArrows && spec.exc.force) {
446507 const f = spec.exc.force;
447508 const t = frames ? this.tl.times[this.tl.frame] : f.tau / 2;
test/drum.test.tsadded+48−0View file
@@ -0,0 +1,48 @@
1+import { describe, expect, it } from 'vitest';
2+import { buildDrum, buildDrumExcitation, buildDrumPickup, deriveDrum, hexLattice, type DrumParams } from '../src/instruments/bead-drum/model';
3+import { DRUM_PRESETS } from '../src/instruments/bead-drum/presets';
4+import { simulate } from '../src/sim/solver';
5+import { spectrum } from '../src/ui/fft';
6+
7+const base = (over: Partial<DrumParams> = {}): DrumParams => ({ ...DRUM_PRESETS[0].params, ...over });
8+
9+describe('hex lattice', () => {
10+ it('gives every bead six springs and is symmetric', () => {
11+ const L = hexLattice(0.15, 8);
12+ const deg = new Array(L.pos.length).fill(0);
13+ for (const [i, j] of L.edges) {
14+ deg[i]++;
15+ if (j < L.pos.length) deg[j]++;
16+ }
17+ expect(deg.every((d) => d === 6)).toBe(true);
18+ let cx = 0, cz = 0;
19+ for (const [x, z] of L.pos) { cx += x; cz += z; }
20+ expect(Math.abs(cx) + Math.abs(cz)).toBeLessThan(1e-9);
21+ });
22+});
23+
24+function run(p: DrumParams) {
25+ const d = deriveDrum(p);
26+ const r = simulate(buildDrum(d), buildDrumExcitation(p, d), { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: buildDrumPickup(p, d) });
27+ return { d, r };
28+}
29+
30+describe('drum', () => {
31+ it('small centre press rings near the ideal membrane fundamental', () => {
32+ const { d, r } = run(base({ excitation: 'press', pressAmp: 1e-5, excR: 0, excWidth: 0.03, t60Low: 0, t60High: 0, duration: 0.25, frameRate: 100 }));
33+ const s = spectrum(r.audio, r.fs);
34+ let best = 0;
35+ for (let k = 1; k < s.freq.length && s.freq[k] < 400; k++) if (s.db[k] > s.db[best] || best === 0) best = k;
36+ console.log('f01 measured', s.freq[best].toFixed(1), 'ideal', d.idealF01.toFixed(1));
37+ expect(Math.abs(s.freq[best] / d.idealF01 - 1)).toBeLessThan(0.05);
38+ });
39+ it('mallet bounces off and energy balances', () => {
40+ const { r } = run(base({ t60Low: 0, t60High: 0, duration: 0.05, frameRate: 1000 }));
41+ const cf = r.contactForce!;
42+ const nz = cf.reduce((c, v) => c + (v > 0 ? 1 : 0), 0);
43+ console.log('contact ms', ((nz / r.fs) * 1000).toFixed(2), 'peak N', Math.max(...cf).toFixed(1), 'balance', Math.max(...r.balance.map(Math.abs)) / Math.max(...r.energy));
44+ expect(cf[cf.length - 1]).toBe(0);
45+ expect(r.malletPos![r.malletPos!.length - 1]).toBeLessThan(0);
46+ expect(Math.max(...r.balance.map(Math.abs)) / Math.max(...r.energy)).toBeLessThan(1e-2);
47+ });
48+});