Woodblock strike simulator: Python reference model and interactive web app
Physical model (woodblock.py, woodblock.md): Hunt-Crossley mallet contact,
free-free bar modes two-way coupled to a Helmholtz cavity, retarded-monopole
radiation, anti-aliased decimation.
Web app (webapp/): TypeScript port of the model running in a Web Worker with
tiered caching (solve / radiate / output), three.js geometry view, parameter
panel with localStorage persistence, waveform/spectrum/spectrogram/contact
force/mode shape visualizations, loop playback, and WAV export. Physics
validated against the Python reference diagnostics and listener targets.
32 changed files+5989−0
webapp/.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+node_modules/
2+dist/
webapp/README.mdadded+27−0View file
@@ -0,0 +1,27 @@
1+# Woodblock Simulator — web app
2+
3+Interactive browser port of the physical model in `../woodblock.py` (see `../woodblock.md` for the physics). The full chain runs client-side in a Web Worker: Hunt-Crossley mallet contact, free-free bar modes two-way coupled to a Helmholtz cavity, retarded-monopole radiation to the listener point, and Kaiser-windowed anti-aliased decimation to 48 kHz.
4+
5+## Usage
6+
7+```
8+npm install
9+npm run dev # development server
10+npm run build # production build in dist/
11+npm run selftest # validate the physics port against the Python reference values
12+```
13+
14+The self-test reproduces the reference diagnostics (partials 1003/2766/5423/8964/13390/18702 Hz, Helmholtz 484 Hz, contact 0.301 ms at 794 N, 109.7 dB SPL) and the four validation listener positions from `woodblock.md` §7.
15+
16+## How rendering is tiered
17+
18+Changing any parameter triggers a debounced fast preview (oversample 8, 24 radiating elements); the Render button runs the full-quality pass (oversample 16, 48 elements). While the current result is a preview, the Play and Download buttons and the status badge say so. The worker keys its caches by tier: geometry, material, damping, mallet, and cavity parameters invalidate the modal solve; listener and slit positions and block height only re-run radiation; pre-delay, normalization, and bit depth only re-run the output stage. Small LRU caches mean that moving the microphone after a full render reuses the cached solve and returns a full-quality result in a fraction of a second.
19+
20+The Loop button repeats the strike continuously, swapping in each new render as parameters change, so edits are audible without pressing Play. Parameters persist in localStorage; "Reset to defaults" restores the values from `woodblock.py`.
21+
22+## Layout
23+
24+- `src/params.ts` — parameter schema, defaults, tier assignment, persistence
25+- `src/physics/` — basis, RK4 solver, radiation, decimation (straight ports of the Python)
26+- `src/worker/renderWorker.ts` — off-thread rendering with tiered caching
27+- `src/components/` — 3D scene (three.js), parameter panel, plots (waveform, spectrum, spectrogram, contact force, animated mode shapes)
webapp/index.htmladded+13−0View file
@@ -0,0 +1,13 @@
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+ <link rel="icon" href="data:image/svg+xml,<svg xmlns=%22http://www.w3.org/2000/svg%22 viewBox=%220 0 100 100%22><text y=%22.9em%22 font-size=%2290%22>🪵</text></svg>" />
7+ <title>Woodblock Simulator</title>
8+ </head>
9+ <body>
10+ <div id="root"></div>
11+ <script type="module" src="/src/main.tsx"></script>
12+ </body>
13+</html>
webapp/package-lock.jsonadded+2361−0View file
This diff is 2,366 lines long and is not shown.
webapp/package.jsonadded+26−0View file
@@ -0,0 +1,26 @@
1+{
2+ "name": "woodblock-webapp",
3+ "private": true,
4+ "version": "0.1.0",
5+ "type": "module",
6+ "scripts": {
7+ "dev": "vite",
8+ "build": "tsc --noEmit && vite build",
9+ "preview": "vite preview",
10+ "selftest": "tsx src/physics/selftest.ts"
11+ },
12+ "dependencies": {
13+ "react": "^18.3.1",
14+ "react-dom": "^18.3.1",
15+ "three": "^0.168.0"
16+ },
17+ "devDependencies": {
18+ "@types/react": "^18.3.5",
19+ "@types/react-dom": "^18.3.0",
20+ "@types/three": "^0.168.0",
21+ "@vitejs/plugin-react": "^4.3.1",
22+ "tsx": "^4.19.0",
23+ "typescript": "^5.5.4",
24+ "vite": "^5.4.3"
25+ }
26+}
webapp/src/App.tsxadded+193−0View file
@@ -0,0 +1,193 @@
1+import { useCallback, useEffect, useRef, useState } from 'react';
2+import { DEFAULTS, loadParams, saveParams, type WoodblockParams } from './params';
3+import { useRenderer } from './hooks/useRenderer';
4+import { encodeWav } from './lib/wav';
5+import ParamPanel from './components/ParamPanel';
6+import Scene3D from './components/Scene3D';
7+import DiagnosticsPanel from './components/DiagnosticsPanel';
8+import WaveformPlot from './components/plots/WaveformPlot';
9+import SpectrumPlot from './components/plots/SpectrumPlot';
10+import SpectrogramPlot from './components/plots/SpectrogramPlot';
11+import ContactForcePlot from './components/plots/ContactForcePlot';
12+import ModeShapes from './components/plots/ModeShapes';
13+
14+export default function App() {
15+ const [params, setParams] = useState<WoodblockParams>(() => loadParams());
16+ const { result, busy, stage, progress, renderFull } = useRenderer(params);
17+ const audioCtxRef = useRef<AudioContext | null>(null);
18+
19+ useEffect(() => saveParams(params), [params]);
20+
21+ const onChange = useCallback((key: keyof WoodblockParams, value: number | boolean) => {
22+ setParams((p) => ({ ...p, [key]: value }));
23+ }, []);
24+
25+ const reset = useCallback(() => setParams({ ...DEFAULTS }), []);
26+
27+ const [looping, setLooping] = useState(false);
28+ const loopSrcRef = useRef<AudioBufferSourceNode | null>(null);
29+
30+ const getAudioCtx = useCallback((sampleRate: number) => {
31+ let ctx = audioCtxRef.current;
32+ if (!ctx || ctx.state === 'closed') {
33+ ctx = new AudioContext({ sampleRate });
34+ audioCtxRef.current = ctx;
35+ }
36+ void ctx.resume();
37+ return ctx;
38+ }, []);
39+
40+ const makeBuffer = useCallback((ctx: AudioContext, r: NonNullable<typeof result>) => {
41+ const buf = ctx.createBuffer(1, r.y.length, r.sampleRate);
42+ buf.copyToChannel(new Float32Array(r.y), 0);
43+ return buf;
44+ }, []);
45+
46+ const play = useCallback(() => {
47+ if (!result) return;
48+ const ctx = getAudioCtx(result.sampleRate);
49+ const src = ctx.createBufferSource();
50+ src.buffer = makeBuffer(ctx, result);
51+ src.connect(ctx.destination);
52+ src.start();
53+ }, [result, getAudioCtx, makeBuffer]);
54+
55+ // Loop mode: repeat the strike continuously, swapping in each new render as
56+ // parameters change so edits are audible without pressing Play.
57+ useEffect(() => {
58+ if (!looping || !result) return;
59+ const ctx = getAudioCtx(result.sampleRate);
60+ const src = ctx.createBufferSource();
61+ src.buffer = makeBuffer(ctx, result);
62+ src.loop = true;
63+ src.connect(ctx.destination);
64+ src.start();
65+ loopSrcRef.current = src;
66+ return () => {
67+ loopSrcRef.current = null;
68+ try {
69+ src.stop();
70+ } catch {
71+ // already stopped
72+ }
73+ src.disconnect();
74+ };
75+ }, [looping, result, getAudioCtx, makeBuffer]);
76+
77+ const download = useCallback(() => {
78+ if (!result) return;
79+ const blob = encodeWav(result.y, result.sampleRate, params.bit_depth === 16 ? 16 : 24);
80+ const url = URL.createObjectURL(blob);
81+ const a = document.createElement('a');
82+ a.href = url;
83+ a.download = result.quality === 'full' ? 'woodblock.wav' : 'woodblock_preview.wav';
84+ a.click();
85+ setTimeout(() => URL.revokeObjectURL(url), 5000);
86+ }, [result, params.bit_depth]);
87+
88+ return (
89+ <div className="app">
90+ <div className="topbar">
91+ <h1>Woodblock Simulator</h1>
92+ <button className="primary" onClick={renderFull} disabled={busy}>
93+ Render
94+ </button>
95+ <button
96+ className={`play ${result?.quality === 'preview' ? 'preview' : ''}`}
97+ onClick={play}
98+ disabled={!result}
99+ title={
100+ result?.quality === 'preview'
101+ ? 'Playing the fast preview render; press Render for full quality'
102+ : 'Play the rendered strike'
103+ }
104+ >
105+ ▶ Play{result?.quality === 'preview' ? ' (preview)' : ''}
106+ </button>
107+ <button
108+ className={`loop ${looping ? 'active' : ''}`}
109+ onClick={() => setLooping((v) => !v)}
110+ disabled={!result && !looping}
111+ title="Repeat the strike continuously so parameter changes are audible immediately"
112+ >
113+ ⟳ Loop{looping ? ' on' : ''}
114+ </button>
115+ <button onClick={download} disabled={!result}>
116+ Download .wav{result?.quality === 'preview' ? ' (preview)' : ''}
117+ </button>
118+ <button onClick={reset}>Reset to defaults</button>
119+ <div className="status">
120+ {result && (
121+ <span className={`badge ${result.quality}`}>
122+ {result.quality === 'full' ? 'full quality' : 'preview quality — press Render'}
123+ </span>
124+ )}
125+ {busy ? (
126+ <>
127+ <span>{stage}…</span>
128+ <div className="progressbar">
129+ <div style={{ width: `${Math.round(progress * 100)}%` }} />
130+ </div>
131+ </>
132+ ) : (
133+ <span>ready</span>
134+ )}
135+ </div>
136+ </div>
137+ <div className="main">
138+ <div className="sidebar">
139+ <ParamPanel params={params} onChange={onChange} />
140+ </div>
141+ <div className="content">
142+ <div className="card wide">
143+ <h2>Geometry — block, strike, slit, microphone</h2>
144+ <Scene3D params={params} />
145+ </div>
146+ {result && (
147+ <>
148+ <div className="card">
149+ <h2>Diagnostics</h2>
150+ <DiagnosticsPanel result={result} />
151+ </div>
152+ <div className="card">
153+ <h2>Contact force</h2>
154+ <ContactForcePlot
155+ force={result.force}
156+ forceFsHz={result.forceFsHz}
157+ peakForceN={result.diag.peakForceN}
158+ contactMs={result.diag.contactMs}
159+ />
160+ </div>
161+ <div className="card">
162+ <h2>Waveform</h2>
163+ <WaveformPlot y={result.y} sampleRate={result.sampleRate} />
164+ </div>
165+ <div className="card">
166+ <h2>Spectrum</h2>
167+ <SpectrumPlot
168+ y={result.y}
169+ sampleRate={result.sampleRate}
170+ partialsHz={result.diag.partialsHz}
171+ helmholtzHz={result.diag.helmholtzHz}
172+ />
173+ </div>
174+ <div className="card">
175+ <h2>Spectrogram</h2>
176+ <SpectrogramPlot y={result.y} sampleRate={result.sampleRate} />
177+ </div>
178+ <div className="card">
179+ <h2>Mode shapes</h2>
180+ <ModeShapes
181+ modeShapes={result.modeShapes}
182+ shapePts={result.shapePts}
183+ partialsHz={result.diag.partialsHz}
184+ strikePosition={params.strike_position}
185+ />
186+ </div>
187+ </>
188+ )}
189+ </div>
190+ </div>
191+ </div>
192+ );
193+}
webapp/src/components/DiagnosticsPanel.tsxadded+56−0View file
@@ -0,0 +1,56 @@
1+import type { RenderResponse } from '../worker/renderWorker';
2+
3+export default function DiagnosticsPanel({ result }: { result: RenderResponse }) {
4+ const d = result.diag;
5+ return (
6+ <div className="diag">
7+ <table>
8+ <tbody>
9+ <tr>
10+ <th>Partials</th>
11+ <td className="partials">
12+ {d.partialsHz.map((f) => f.toFixed(0)).join(', ')} Hz
13+ </td>
14+ </tr>
15+ <tr>
16+ <th>Ratios</th>
17+ <td className="partials">{d.ratios.map((r) => r.toFixed(2)).join(', ')}</td>
18+ </tr>
19+ <tr>
20+ <th>Helmholtz</th>
21+ <td>{d.helmholtzHz.toFixed(0)} Hz</td>
22+ </tr>
23+ <tr>
24+ <th>Contact</th>
25+ <td>
26+ {d.contactMs.toFixed(3)} ms, peak {d.peakForceN.toFixed(0)} N
27+ </td>
28+ </tr>
29+ <tr>
30+ <th>Listener</th>
31+ <td>{(d.listenerDistM * 100).toFixed(0)} cm from block centre</td>
32+ </tr>
33+ <tr>
34+ <th>Peak pressure</th>
35+ <td>
36+ {d.peakPa.toFixed(3)} Pa ({d.peakSplDb.toFixed(1)} dB SPL)
37+ </td>
38+ </tr>
39+ <tr>
40+ <th>Quality</th>
41+ <td>
42+ {result.quality} · oversample {d.oversample} · {d.elements} elements
43+ </td>
44+ </tr>
45+ </tbody>
46+ </table>
47+ <div className="cachenote">
48+ solve: {d.solveCached ? 'cached' : `${(d.solveMs / 1000).toFixed(2)} s`} · radiation:{' '}
49+ {d.radiateCached ? 'cached' : `${(d.radiateMs / 1000).toFixed(2)} s`}
50+ </div>
51+ {d.warnings.map((w, i) => (
52+ <div key={i} className="warning">⚠ {w}</div>
53+ ))}
54+ </div>
55+ );
56+}
webapp/src/components/ParamPanel.tsxadded+123−0View file
@@ -0,0 +1,123 @@
1+import { PARAM_GROUPS, type ParamSpec, type WoodblockParams } from '../params';
2+
3+interface Props {
4+ params: WoodblockParams;
5+ onChange: (key: keyof WoodblockParams, value: number | boolean) => void;
6+}
7+
8+const SLIDER_STEPS = 1000;
9+
10+function sliderFromValue(spec: ParamSpec, v: number): number {
11+ if (spec.log) {
12+ const t = (Math.log(v) - Math.log(spec.min)) / (Math.log(spec.max) - Math.log(spec.min));
13+ return Math.round(t * SLIDER_STEPS);
14+ }
15+ return Math.round(((v - spec.min) / (spec.max - spec.min)) * SLIDER_STEPS);
16+}
17+
18+function valueFromSlider(spec: ParamSpec, s: number): number {
19+ const t = s / SLIDER_STEPS;
20+ let v: number;
21+ if (spec.log) {
22+ v = Math.exp(Math.log(spec.min) + t * (Math.log(spec.max) - Math.log(spec.min)));
23+ // keep three significant digits so the displayed number stays readable
24+ const mag = 10 ** Math.floor(Math.log10(Math.abs(v)));
25+ v = Math.round((v / mag) * 100) / 100 * mag;
26+ } else {
27+ v = spec.min + t * (spec.max - spec.min);
28+ if (spec.step) v = Math.round(v / spec.step) * spec.step;
29+ }
30+ if (spec.integer) v = Math.round(v);
31+ return Math.min(spec.max, Math.max(spec.min, v));
32+}
33+
34+function formatValue(v: number): string {
35+ if (v === 0) return '0';
36+ const a = Math.abs(v);
37+ if (a >= 1e5 || a < 1e-3) return v.toExponential(2);
38+ return String(Math.round(v * 1e6) / 1e6);
39+}
40+
41+function ParamRow({ spec, params, onChange }: Props & { spec: ParamSpec }) {
42+ const v = params[spec.key] as number;
43+ return (
44+ <div className="param-row" title={spec.help ?? ''}>
45+ <div className="param-label">
46+ <span className="name">
47+ {spec.label}
48+ {spec.unit ? <span style={{ color: 'var(--muted)' }}> ({spec.unit})</span> : null}
49+ </span>
50+ </div>
51+ <div className="param-controls">
52+ <input
53+ type="range"
54+ min={0}
55+ max={SLIDER_STEPS}
56+ value={sliderFromValue(spec, v)}
57+ onChange={(e) => onChange(spec.key, valueFromSlider(spec, Number(e.target.value)))}
58+ />
59+ <input
60+ type="number"
61+ value={formatValue(v)}
62+ step={spec.step ?? 'any'}
63+ onChange={(e) => {
64+ const nv = Number(e.target.value);
65+ if (Number.isFinite(nv)) onChange(spec.key, spec.integer ? Math.round(nv) : nv);
66+ }}
67+ />
68+ </div>
69+ </div>
70+ );
71+}
72+
73+export default function ParamPanel({ params, onChange }: Props) {
74+ return (
75+ <div>
76+ {PARAM_GROUPS.map((g, gi) => (
77+ <details key={g.name} className="group" open={gi < 2}>
78+ <summary>{g.name}</summary>
79+ <div className="group-body">
80+ {g.specs.map((spec) => (
81+ <ParamRow key={spec.key} spec={spec} params={params} onChange={onChange} />
82+ ))}
83+ {g.name === 'Output' && (
84+ <>
85+ <div className="param-row">
86+ <div className="param-label">
87+ <span className="name">Normalize peak</span>
88+ </div>
89+ <div className="param-controls">
90+ <input
91+ type="checkbox"
92+ checked={params.normalize}
93+ onChange={(e) => onChange('normalize', e.target.checked)}
94+ />
95+ <span className="param-help">
96+ {params.normalize
97+ ? `scaled to ${params.peak_dbfs} dBFS`
98+ : `absolute: ${params.full_scale_pa} Pa = 0 dBFS`}
99+ </span>
100+ </div>
101+ </div>
102+ <div className="param-row">
103+ <div className="param-label">
104+ <span className="name">WAV bit depth</span>
105+ </div>
106+ <div className="param-controls">
107+ <select
108+ value={params.bit_depth}
109+ onChange={(e) => onChange('bit_depth', Number(e.target.value))}
110+ >
111+ <option value={16}>16-bit</option>
112+ <option value={24}>24-bit</option>
113+ </select>
114+ </div>
115+ </div>
116+ </>
117+ )}
118+ </div>
119+ </details>
120+ ))}
121+ </div>
122+ );
123+}
webapp/src/components/Scene3D.tsxadded+269−0View file
@@ -0,0 +1,269 @@
1+// 3D view of the geometry: block, cavity hint, slit, strike point + mallet,
2+// listener, radiating surface elements. Model coordinates are block-centred
3+// with +z up out of the struck face; three.js is y-up, so we map
4+// (x, y, z)_model -> (x, z, -y)_three, a proper rotation.
5+
6+import { useEffect, useRef } from 'react';
7+import * as THREE from 'three';
8+import { OrbitControls } from 'three/examples/jsm/controls/OrbitControls.js';
9+import type { WoodblockParams } from '../params';
10+
11+function m2t(x: number, y: number, z: number): [number, number, number] {
12+ return [x, z, -y];
13+}
14+
15+function makeLabel(text: string, color: string): THREE.Sprite {
16+ const canvas = document.createElement('canvas');
17+ canvas.width = 320;
18+ canvas.height = 64;
19+ const ctx = canvas.getContext('2d')!;
20+ ctx.font = '600 36px system-ui, sans-serif';
21+ ctx.fillStyle = color;
22+ ctx.textAlign = 'center';
23+ ctx.textBaseline = 'middle';
24+ ctx.fillText(text, 160, 32);
25+ const tex = new THREE.CanvasTexture(canvas);
26+ const mat = new THREE.SpriteMaterial({ map: tex, depthTest: false, transparent: true });
27+ const sprite = new THREE.Sprite(mat);
28+ sprite.scale.set(0.15, 0.03, 1);
29+ return sprite;
30+}
31+
32+interface Props {
33+ params: WoodblockParams;
34+}
35+
36+export default function Scene3D({ params }: Props) {
37+ const mountRef = useRef<HTMLDivElement | null>(null);
38+ const sceneRef = useRef<{
39+ scene: THREE.Scene;
40+ group: THREE.Group | null;
41+ renderer: THREE.WebGLRenderer;
42+ camera: THREE.PerspectiveCamera;
43+ } | null>(null);
44+
45+ // one-time scene setup
46+ useEffect(() => {
47+ const mount = mountRef.current!;
48+ const renderer = new THREE.WebGLRenderer({ antialias: true });
49+ renderer.setPixelRatio(window.devicePixelRatio);
50+ renderer.setClearColor(0x14161a);
51+ mount.appendChild(renderer.domElement);
52+
53+ const scene = new THREE.Scene();
54+ const camera = new THREE.PerspectiveCamera(45, 1, 0.01, 20);
55+ camera.position.set(0.45, 0.35, 0.6);
56+
57+ const controls = new OrbitControls(camera, renderer.domElement);
58+ controls.target.set(0, 0.06, 0);
59+ controls.enableDamping = true;
60+
61+ scene.add(new THREE.AmbientLight(0xffffff, 0.55));
62+ const dir = new THREE.DirectionalLight(0xffffff, 1.1);
63+ dir.position.set(0.5, 1.0, 0.6);
64+ scene.add(dir);
65+ const dir2 = new THREE.DirectionalLight(0x8899ff, 0.3);
66+ dir2.position.set(-0.6, 0.3, -0.5);
67+ scene.add(dir2);
68+
69+ const grid = new THREE.GridHelper(2, 40, 0x3a4048, 0x262b33);
70+ grid.position.y = -0.101;
71+ scene.add(grid);
72+
73+ // model axes: x red (along bar), y green (across width), z blue (up)
74+ const axisLen = 0.14;
75+ const axes: [THREE.Vector3, number, string][] = [
76+ [new THREE.Vector3(...m2t(axisLen, 0, 0)), 0xe57373, '+x'],
77+ [new THREE.Vector3(...m2t(0, axisLen, 0)), 0x81c784, '+y'],
78+ [new THREE.Vector3(...m2t(0, 0, axisLen)), 0x64b5f6, '+z'],
79+ ];
80+ for (const [v, color, label] of axes) {
81+ const geo = new THREE.BufferGeometry().setFromPoints([new THREE.Vector3(0, 0, 0), v]);
82+ scene.add(new THREE.Line(geo, new THREE.LineBasicMaterial({ color })));
83+ const s = makeLabel(label, '#' + color.toString(16).padStart(6, '0'));
84+ s.position.copy(v).multiplyScalar(1.15);
85+ scene.add(s);
86+ }
87+
88+ sceneRef.current = { scene, group: null, renderer, camera };
89+
90+ let raf = 0;
91+ const animate = () => {
92+ raf = requestAnimationFrame(animate);
93+ controls.update();
94+ renderer.render(scene, camera);
95+ };
96+ animate();
97+
98+ const resize = () => {
99+ const w = mount.clientWidth;
100+ const h = mount.clientHeight;
101+ renderer.setSize(w, h, false);
102+ camera.aspect = w / h;
103+ camera.updateProjectionMatrix();
104+ };
105+ resize();
106+ const ro = new ResizeObserver(resize);
107+ ro.observe(mount);
108+
109+ return () => {
110+ cancelAnimationFrame(raf);
111+ ro.disconnect();
112+ controls.dispose();
113+ renderer.dispose();
114+ mount.removeChild(renderer.domElement);
115+ sceneRef.current = null;
116+ };
117+ }, []);
118+
119+ // rebuild the geometry group whenever relevant parameters change
120+ useEffect(() => {
121+ const ctx = sceneRef.current;
122+ if (!ctx) return;
123+ if (ctx.group) {
124+ ctx.scene.remove(ctx.group);
125+ ctx.group.traverse((o) => {
126+ const mesh = o as THREE.Mesh;
127+ if (mesh.geometry) mesh.geometry.dispose();
128+ });
129+ }
130+
131+ const g = new THREE.Group();
132+ const L = params.length_m;
133+ const W = params.width_m;
134+ const Hb = params.block_height_m;
135+
136+ // block body
137+ const wood = new THREE.MeshStandardMaterial({
138+ color: 0xa06f45,
139+ roughness: 0.65,
140+ metalness: 0.0,
141+ transparent: true,
142+ opacity: 0.92,
143+ });
144+ const block = new THREE.Mesh(new THREE.BoxGeometry(L, Hb, W), wood);
145+ g.add(block);
146+ const edges = new THREE.LineSegments(
147+ new THREE.EdgesGeometry(new THREE.BoxGeometry(L, Hb, W)),
148+ new THREE.LineBasicMaterial({ color: 0x5c4026 }),
149+ );
150+ g.add(edges);
151+
152+ // cavity hint: darker inner box over the cavity wall span
153+ const lo = params.cavity_span_lo;
154+ const hi = params.cavity_span_hi;
155+ const cavLen = Math.max(0.001, (hi - lo) * L);
156+ const cavCentre = ((lo + hi) / 2 - 0.5) * L;
157+ const cavity = new THREE.Mesh(
158+ new THREE.BoxGeometry(cavLen, Hb * 0.55, W * 0.72),
159+ new THREE.MeshStandardMaterial({ color: 0x241a10, roughness: 0.9 }),
160+ );
161+ cavity.position.set(...m2t(cavCentre, 0, 0));
162+ g.add(cavity);
163+
164+ // slit opening marker: dark patch + small sphere at the monopole point
165+ const slitPos = new THREE.Vector3(...m2t(params.slit_x_m, params.slit_y_m, params.slit_z_m));
166+ const slitArea = params.slit_area_m2;
167+ const slitLen = Math.min(L * 0.6, Math.sqrt(slitArea * 10)); // display proportions only
168+ const slitWid = slitArea / slitLen;
169+ const slit = new THREE.Mesh(
170+ new THREE.BoxGeometry(slitLen, slitWid, 0.002),
171+ new THREE.MeshStandardMaterial({ color: 0x0d0a06, roughness: 1.0 }),
172+ );
173+ slit.position.copy(slitPos);
174+ // orient the patch to lie in the wall nearest the slit point
175+ if (Math.abs(params.slit_z_m) >= Math.abs(params.slit_y_m)) slit.rotation.x = Math.PI / 2;
176+ g.add(slit);
177+ const slitDot = new THREE.Mesh(
178+ new THREE.SphereGeometry(0.005, 16, 12),
179+ new THREE.MeshStandardMaterial({ color: 0xba68c8, emissive: 0x4a1f55 }),
180+ );
181+ slitDot.position.copy(slitPos);
182+ g.add(slitDot);
183+ const slitLabel = makeLabel('slit', '#ba68c8');
184+ slitLabel.position.copy(slitPos).add(new THREE.Vector3(0, 0.03, 0));
185+ g.add(slitLabel);
186+
187+ // strike point + mallet
188+ const xc = (params.strike_position - 0.5) * L;
189+ const strikeM = m2t(xc, 0, Hb / 2);
190+ const strikeDot = new THREE.Mesh(
191+ new THREE.SphereGeometry(0.006, 16, 12),
192+ new THREE.MeshStandardMaterial({ color: 0xef5350, emissive: 0x5a1210 }),
193+ );
194+ strikeDot.position.set(...strikeM);
195+ g.add(strikeDot);
196+
197+ const malletBall = new THREE.Mesh(
198+ new THREE.SphereGeometry(0.014, 20, 16),
199+ new THREE.MeshStandardMaterial({ color: 0xd7ccc8, roughness: 0.5 }),
200+ );
201+ malletBall.position.set(strikeM[0], strikeM[1] + 0.05, strikeM[2]);
202+ g.add(malletBall);
203+ const shaft = new THREE.Mesh(
204+ new THREE.CylinderGeometry(0.004, 0.004, 0.16, 12),
205+ new THREE.MeshStandardMaterial({ color: 0x8d6e63, roughness: 0.7 }),
206+ );
207+ shaft.position.set(strikeM[0] + 0.055, strikeM[1] + 0.115, strikeM[2]);
208+ shaft.rotation.z = -0.75;
209+ g.add(shaft);
210+ const strikeLabel = makeLabel('strike', '#ef5350');
211+ strikeLabel.position.set(strikeM[0] - 0.01, strikeM[1] + 0.09, strikeM[2]);
212+ g.add(strikeLabel);
213+
214+ // radiating surface elements on both faces
215+ const M = params.radiating_elements;
216+ const pts: number[] = [];
217+ for (let j = 0; j < M; j++) {
218+ const x = ((j + 0.5) / M - 0.5) * L;
219+ pts.push(...m2t(x, 0, Hb / 2 + 0.001));
220+ pts.push(...m2t(x, 0, -Hb / 2 - 0.001));
221+ }
222+ const ptsGeo = new THREE.BufferGeometry();
223+ ptsGeo.setAttribute('position', new THREE.Float32BufferAttribute(pts, 3));
224+ g.add(new THREE.Points(ptsGeo, new THREE.PointsMaterial({ color: 0xffe082, size: 0.004 })));
225+
226+ // listener: microphone marker + line from block centre
227+ const lis = new THREE.Vector3(
228+ ...m2t(params.listener_x_m, params.listener_y_m, params.listener_z_m),
229+ );
230+ const micHead = new THREE.Mesh(
231+ new THREE.SphereGeometry(0.016, 20, 16),
232+ new THREE.MeshStandardMaterial({ color: 0x4dd0e1, emissive: 0x0c3b42, roughness: 0.35 }),
233+ );
234+ micHead.position.copy(lis);
235+ g.add(micHead);
236+ const micBody = new THREE.Mesh(
237+ new THREE.CylinderGeometry(0.008, 0.008, 0.05, 12),
238+ new THREE.MeshStandardMaterial({ color: 0x455a64, roughness: 0.6 }),
239+ );
240+ micBody.position.copy(lis).add(new THREE.Vector3(0, -0.033, 0));
241+ g.add(micBody);
242+ const lineGeo = new THREE.BufferGeometry().setFromPoints([new THREE.Vector3(0, 0, 0), lis]);
243+ const line = new THREE.Line(
244+ lineGeo,
245+ new THREE.LineDashedMaterial({ color: 0x4dd0e1, dashSize: 0.012, gapSize: 0.009 }),
246+ );
247+ line.computeLineDistances();
248+ g.add(line);
249+ const dist = Math.hypot(params.listener_x_m, params.listener_y_m, params.listener_z_m);
250+ const micLabel = makeLabel(`mic ${(dist * 100).toFixed(0)} cm`, '#4dd0e1');
251+ micLabel.position.copy(lis).add(new THREE.Vector3(0, 0.045, 0));
252+ g.add(micLabel);
253+
254+ ctx.scene.add(g);
255+ ctx.group = g;
256+ }, [
257+ params.length_m, params.width_m, params.block_height_m,
258+ params.cavity_span_lo, params.cavity_span_hi,
259+ params.slit_x_m, params.slit_y_m, params.slit_z_m, params.slit_area_m2,
260+ params.strike_position, params.radiating_elements,
261+ params.listener_x_m, params.listener_y_m, params.listener_z_m,
262+ ]);
263+
264+ return (
265+ <div className="scene3d" ref={mountRef}>
266+ <div className="scene-hint">drag to orbit · scroll to zoom · right-drag to pan</div>
267+ </div>
268+ );
269+}
webapp/src/components/plots/ContactForcePlot.tsxadded+73−0View file
@@ -0,0 +1,73 @@
1+// Contact-force pulse f_c(t). Its width sets brightness: harder mallets give
2+// a shorter pulse and therefore more high-frequency drive.
3+
4+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
5+
6+interface Props {
7+ force: Float32Array;
8+ forceFsHz: number;
9+ peakForceN: number;
10+ contactMs: number;
11+}
12+
13+export default function ContactForcePlot({ force, forceFsHz, peakForceN, contactMs }: Props) {
14+ const ref = useCanvas(
15+ (ctx, w, h) => {
16+ drawFrame(ctx, w, h);
17+ const padL = 44;
18+ const padB = 16;
19+ const pw = w - padL - 8;
20+ const ph = h - padB - 10;
21+
22+ // window: from 0 to a bit past the end of contact, at least 1 ms
23+ let lastIdx = 0;
24+ for (let i = 0; i < force.length; i++) if (force[i] > 0) lastIdx = i;
25+ const winMs = Math.max(1.0, ((lastIdx / forceFsHz) * 1e3) * 1.4);
26+ const nWin = Math.min(force.length, Math.round((winMs / 1e3) * forceFsHz));
27+
28+ let peak = 0;
29+ for (let i = 0; i < nWin; i++) peak = Math.max(peak, force[i]);
30+ if (peak <= 0) peak = 1;
31+
32+ ctx.font = '10px system-ui';
33+ ctx.fillStyle = PLOT_COLORS.axis;
34+ ctx.textAlign = 'center';
35+ const step = winMs > 4 ? 1.0 : winMs > 2 ? 0.5 : 0.2;
36+ for (let t = 0; t <= winMs + 1e-9; t += step) {
37+ const x = padL + (t / winMs) * pw;
38+ ctx.strokeStyle = PLOT_COLORS.grid;
39+ ctx.beginPath();
40+ ctx.moveTo(x, 10);
41+ ctx.lineTo(x, 10 + ph);
42+ ctx.stroke();
43+ ctx.fillText(`${t.toFixed(1)}ms`, x, h - 4);
44+ }
45+ ctx.textAlign = 'right';
46+ ctx.fillText(`${peak.toFixed(0)} N`, padL - 4, 16);
47+ ctx.fillText('0', padL - 4, 10 + ph);
48+
49+ ctx.strokeStyle = PLOT_COLORS.accent;
50+ ctx.lineWidth = 1.5;
51+ ctx.beginPath();
52+ for (let px = 0; px <= pw; px++) {
53+ const i = Math.min(nWin - 1, Math.round((px / pw) * (nWin - 1)));
54+ const yy = 10 + ph - (force[i] / peak) * ph;
55+ if (px === 0) ctx.moveTo(padL, yy);
56+ else ctx.lineTo(padL + px, yy);
57+ }
58+ ctx.stroke();
59+ ctx.lineWidth = 1;
60+
61+ ctx.fillStyle = PLOT_COLORS.axis;
62+ ctx.textAlign = 'left';
63+ ctx.fillText(
64+ `contact ${contactMs.toFixed(3)} ms · peak ${peakForceN.toFixed(0)} N`,
65+ padL + 6,
66+ 16,
67+ );
68+ },
69+ [force, forceFsHz, peakForceN, contactMs],
70+ 170,
71+ );
72+ return <canvas ref={ref} />;
73+}
webapp/src/components/plots/ModeShapes.tsxadded+130−0View file
@@ -0,0 +1,130 @@
1+// Animated free-free bar mode shapes. One mode is shown large and oscillating;
2+// the strike position is marked so nodes/antinodes can be read against it.
3+
4+import { useEffect, useRef, useState } from 'react';
5+import { PLOT_COLORS } from './canvasUtils';
6+
7+interface Props {
8+ modeShapes: Float32Array; // n x shapePts
9+ shapePts: number;
10+ partialsHz: number[];
11+ strikePosition: number;
12+}
13+
14+export default function ModeShapes({ modeShapes, shapePts, partialsHz, strikePosition }: Props) {
15+ const n = partialsHz.length;
16+ const [selected, setSelected] = useState(0);
17+ const canvasRef = useRef<HTMLCanvasElement | null>(null);
18+ const sel = Math.min(selected, Math.max(0, n - 1));
19+
20+ useEffect(() => {
21+ const canvas = canvasRef.current;
22+ if (!canvas || n === 0) return;
23+ let raf = 0;
24+ const t0 = performance.now();
25+
26+ const draw = () => {
27+ raf = requestAnimationFrame(draw);
28+ const dpr = window.devicePixelRatio || 1;
29+ const w = canvas.clientWidth;
30+ const h = 180;
31+ canvas.width = Math.round(w * dpr);
32+ canvas.height = Math.round(h * dpr);
33+ canvas.style.height = `${h}px`;
34+ const ctx = canvas.getContext('2d')!;
35+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
36+ ctx.fillStyle = PLOT_COLORS.bg;
37+ ctx.fillRect(0, 0, w, h);
38+
39+ const padL = 14;
40+ const pw = w - 2 * padL;
41+ const midY = h / 2;
42+ const amp = h * 0.32;
43+
44+ // oscillation rate: slow, slightly faster for higher modes
45+ const t = (performance.now() - t0) / 1000;
46+ const rate = 0.8 + sel * 0.25;
47+ const osc = Math.sin(2 * Math.PI * rate * t);
48+
49+ // static envelope
50+ ctx.strokeStyle = PLOT_COLORS.grid;
51+ ctx.setLineDash([3, 3]);
52+ for (const sgn of [1, -1]) {
53+ ctx.beginPath();
54+ for (let i = 0; i < shapePts; i++) {
55+ const x = padL + (i / (shapePts - 1)) * pw;
56+ const y = midY - sgn * modeShapes[sel * shapePts + i] * amp;
57+ if (i === 0) ctx.moveTo(x, y);
58+ else ctx.lineTo(x, y);
59+ }
60+ ctx.stroke();
61+ }
62+ ctx.setLineDash([]);
63+
64+ // rest axis
65+ ctx.strokeStyle = PLOT_COLORS.grid;
66+ ctx.beginPath();
67+ ctx.moveTo(padL, midY);
68+ ctx.lineTo(padL + pw, midY);
69+ ctx.stroke();
70+
71+ // strike position
72+ const sx = padL + strikePosition * pw;
73+ ctx.strokeStyle = PLOT_COLORS.warn;
74+ ctx.setLineDash([4, 3]);
75+ ctx.beginPath();
76+ ctx.moveTo(sx, 14);
77+ ctx.lineTo(sx, h - 14);
78+ ctx.stroke();
79+ ctx.setLineDash([]);
80+ ctx.fillStyle = PLOT_COLORS.warn;
81+ ctx.font = '10px system-ui';
82+ ctx.textAlign = 'center';
83+ ctx.fillText('strike', sx, 11);
84+
85+ // animated shape as a thick bar
86+ ctx.strokeStyle = PLOT_COLORS.trace;
87+ ctx.lineWidth = 3;
88+ ctx.beginPath();
89+ for (let i = 0; i < shapePts; i++) {
90+ const x = padL + (i / (shapePts - 1)) * pw;
91+ const y = midY - modeShapes[sel * shapePts + i] * amp * osc;
92+ if (i === 0) ctx.moveTo(x, y);
93+ else ctx.lineTo(x, y);
94+ }
95+ ctx.stroke();
96+ ctx.lineWidth = 1;
97+
98+ // value of the mode at the strike point (drive efficiency)
99+ const idx = Math.round(strikePosition * (shapePts - 1));
100+ const phiS = Math.abs(modeShapes[sel * shapePts + idx]);
101+ ctx.fillStyle = PLOT_COLORS.axis;
102+ ctx.textAlign = 'left';
103+ ctx.fillText(
104+ `mode ${sel + 1} · ${partialsHz[sel].toFixed(0)} Hz · |φ| at strike: ${(phiS * 100).toFixed(0)}%`,
105+ padL,
106+ h - 4,
107+ );
108+ };
109+ draw();
110+ return () => cancelAnimationFrame(raf);
111+ }, [modeShapes, shapePts, partialsHz, strikePosition, sel, n]);
112+
113+ if (n === 0) return null;
114+ return (
115+ <div>
116+ <div className="mode-list">
117+ {partialsHz.map((f, i) => (
118+ <button
119+ key={i}
120+ className={i === sel ? 'active' : ''}
121+ onClick={() => setSelected(i)}
122+ >
123+ {f.toFixed(0)} Hz
124+ </button>
125+ ))}
126+ </div>
127+ <canvas ref={canvasRef} />
128+ </div>
129+ );
130+}
webapp/src/components/plots/SpectrogramPlot.tsxadded+64−0View file
@@ -0,0 +1,64 @@
1+// Time-frequency view of the rendered strike; per-partial decay is visible as
2+// horizontal streaks of different lengths.
3+
4+import { stft } from '../../lib/fft';
5+import { useCanvas, drawFrame, PLOT_COLORS, colormap } from './canvasUtils';
6+
7+interface Props {
8+ y: Float32Array;
9+ sampleRate: number;
10+}
11+
12+export default function SpectrogramPlot({ y, sampleRate }: Props) {
13+ const ref = useCanvas(
14+ (ctx, w, h) => {
15+ drawFrame(ctx, w, h);
16+ const padL = 34;
17+ const padB = 16;
18+ const pw = w - padL - 6;
19+ const ph = h - padB - 6;
20+
21+ const { frames, bins, db, floorDb } = stft(y, 512, 128);
22+
23+ const img = ctx.createImageData(frames, bins);
24+ for (let f = 0; f < frames; f++) {
25+ for (let b = 0; b < bins; b++) {
26+ const v = db[f * bins + b];
27+ const t = Math.min(1, Math.max(0, (v - floorDb) / -floorDb));
28+ const [r, g, bl] = colormap(t);
29+ // flip vertically: low frequency at the bottom
30+ const row = bins - 1 - b;
31+ const idx = (row * frames + f) * 4;
32+ img.data[idx] = r;
33+ img.data[idx + 1] = g;
34+ img.data[idx + 2] = bl;
35+ img.data[idx + 3] = 255;
36+ }
37+ }
38+ const off = new OffscreenCanvas(frames, bins);
39+ const octx = off.getContext('2d')!;
40+ octx.putImageData(img, 0, 0);
41+ ctx.imageSmoothingEnabled = true;
42+ ctx.drawImage(off, padL, 6, pw, ph);
43+
44+ // axes
45+ ctx.fillStyle = PLOT_COLORS.axis;
46+ ctx.font = '10px system-ui';
47+ ctx.textAlign = 'right';
48+ const nyq = sampleRate / 2;
49+ for (let kf = 0; kf <= nyq; kf += 4000) {
50+ const yy = 6 + ph - (kf / nyq) * ph;
51+ ctx.fillText(kf >= 1000 ? `${kf / 1000}k` : `${kf}`, padL - 4, yy + 3);
52+ }
53+ ctx.textAlign = 'center';
54+ const durS = y.length / sampleRate;
55+ for (let t = 0; t <= durS + 1e-9; t += 0.1) {
56+ const x = padL + (t / durS) * pw;
57+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
58+ }
59+ },
60+ [y, sampleRate],
61+ 200,
62+ );
63+ return <canvas ref={ref} />;
64+}
webapp/src/components/plots/SpectrumPlot.tsxadded+111−0View file
@@ -0,0 +1,111 @@
1+// Magnitude spectrum of the rendered strike, log-frequency axis, with the
2+// modal partial frequencies and the Helmholtz resonance overlaid.
3+
4+import { magnitudeSpectrum } from '../../lib/fft';
5+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
6+
7+interface Props {
8+ y: Float32Array;
9+ sampleRate: number;
10+ partialsHz: number[];
11+ helmholtzHz: number;
12+}
13+
14+const F_LO = 100;
15+const F_HI = 24000;
16+const DB_RANGE = 80;
17+
18+export default function SpectrumPlot({ y, sampleRate, partialsHz, helmholtzHz }: Props) {
19+ const ref = useCanvas(
20+ (ctx, w, h) => {
21+ drawFrame(ctx, w, h);
22+ const padL = 34;
23+ const padB = 18;
24+ const pw = w - padL - 6;
25+ const ph = h - padB - 8;
26+ const fx = (f: number) =>
27+ padL + ((Math.log10(f) - Math.log10(F_LO)) / (Math.log10(F_HI) - Math.log10(F_LO))) * pw;
28+
29+ // frequency gridlines
30+ ctx.font = '10px system-ui';
31+ ctx.textAlign = 'center';
32+ for (const f of [100, 200, 500, 1000, 2000, 5000, 10000, 20000]) {
33+ const x = fx(f);
34+ ctx.strokeStyle = PLOT_COLORS.grid;
35+ ctx.beginPath();
36+ ctx.moveTo(x, 8);
37+ ctx.lineTo(x, 8 + ph);
38+ ctx.stroke();
39+ ctx.fillStyle = PLOT_COLORS.axis;
40+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, x, h - 5);
41+ }
42+ ctx.textAlign = 'right';
43+ for (let db = 0; db >= -DB_RANGE; db -= 20) {
44+ const yy = 8 + (-db / DB_RANGE) * ph;
45+ ctx.strokeStyle = PLOT_COLORS.grid;
46+ ctx.beginPath();
47+ ctx.moveTo(padL, yy);
48+ ctx.lineTo(padL + pw, yy);
49+ ctx.stroke();
50+ ctx.fillStyle = PLOT_COLORS.axis;
51+ ctx.fillText(`${db}`, padL - 4, yy + 3);
52+ }
53+
54+ const { mag, size } = magnitudeSpectrum(y, 16384);
55+ let peak = 0;
56+ for (let i = 1; i < mag.length; i++) peak = Math.max(peak, mag[i]);
57+ if (peak <= 0) return;
58+
59+ // overlays first so the trace draws on top
60+ ctx.strokeStyle = PLOT_COLORS.good;
61+ ctx.setLineDash([4, 3]);
62+ if (helmholtzHz >= F_LO && helmholtzHz <= F_HI) {
63+ ctx.beginPath();
64+ ctx.moveTo(fx(helmholtzHz), 8);
65+ ctx.lineTo(fx(helmholtzHz), 8 + ph);
66+ ctx.stroke();
67+ ctx.fillStyle = PLOT_COLORS.good;
68+ ctx.textAlign = 'left';
69+ ctx.fillText(`H ${helmholtzHz.toFixed(0)}`, fx(helmholtzHz) + 3, 18);
70+ }
71+ ctx.strokeStyle = PLOT_COLORS.accent;
72+ for (const f of partialsHz) {
73+ if (f < F_LO || f > F_HI) continue;
74+ ctx.beginPath();
75+ ctx.moveTo(fx(f), 8);
76+ ctx.lineTo(fx(f), 8 + ph);
77+ ctx.stroke();
78+ }
79+ ctx.setLineDash([]);
80+
81+ ctx.strokeStyle = PLOT_COLORS.trace;
82+ ctx.beginPath();
83+ let started = false;
84+ for (let px = 0; px <= pw; px++) {
85+ const f = 10 ** (Math.log10(F_LO) + (px / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
86+ const bin = Math.min(mag.length - 1, Math.round((f / sampleRate) * size));
87+ // take the max over the bins this pixel spans, so narrow partials survive
88+ const fNext = 10 ** (Math.log10(F_LO) + ((px + 1) / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
89+ const binNext = Math.min(mag.length - 1, Math.max(bin + 1, Math.round((fNext / sampleRate) * size)));
90+ let m = 0;
91+ for (let b = bin; b < binNext; b++) m = Math.max(m, mag[b]);
92+ const db = 20 * Math.log10(m / peak + 1e-12);
93+ const yy = 8 + Math.min(ph, (-db / DB_RANGE) * ph);
94+ if (!started) {
95+ ctx.moveTo(padL + px, yy);
96+ started = true;
97+ } else {
98+ ctx.lineTo(padL + px, yy);
99+ }
100+ }
101+ ctx.stroke();
102+
103+ ctx.fillStyle = PLOT_COLORS.accent;
104+ ctx.textAlign = 'left';
105+ ctx.fillText('partials', padL + 4, 18);
106+ },
107+ [y, sampleRate, partialsHz, helmholtzHz],
108+ 200,
109+ );
110+ return <canvas ref={ref} />;
111+}
webapp/src/components/plots/WaveformPlot.tsxadded+64−0View file
@@ -0,0 +1,64 @@
1+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
2+
3+interface Props {
4+ y: Float32Array;
5+ sampleRate: number;
6+}
7+
8+export default function WaveformPlot({ y, sampleRate }: Props) {
9+ const ref = useCanvas(
10+ (ctx, w, h) => {
11+ drawFrame(ctx, w, h);
12+ const padL = 34;
13+ const padB = 16;
14+ const pw = w - padL - 6;
15+ const ph = h - padB - 6;
16+ const midY = 6 + ph / 2;
17+
18+ ctx.strokeStyle = PLOT_COLORS.grid;
19+ ctx.beginPath();
20+ ctx.moveTo(padL, midY);
21+ ctx.lineTo(padL + pw, midY);
22+ ctx.stroke();
23+
24+ // time gridlines every 0.1 s
25+ const durS = y.length / sampleRate;
26+ ctx.fillStyle = PLOT_COLORS.axis;
27+ ctx.font = '10px system-ui';
28+ ctx.textAlign = 'center';
29+ for (let t = 0; t <= durS + 1e-9; t += 0.1) {
30+ const x = padL + (t / durS) * pw;
31+ ctx.strokeStyle = PLOT_COLORS.grid;
32+ ctx.beginPath();
33+ ctx.moveTo(x, 6);
34+ ctx.lineTo(x, 6 + ph);
35+ ctx.stroke();
36+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
37+ }
38+ ctx.textAlign = 'right';
39+ ctx.fillText('+1', padL - 4, 12);
40+ ctx.fillText('-1', padL - 4, 6 + ph);
41+
42+ // min/max envelope per pixel column
43+ ctx.strokeStyle = PLOT_COLORS.trace;
44+ ctx.beginPath();
45+ const perPx = y.length / pw;
46+ for (let px = 0; px < pw; px++) {
47+ const i0 = Math.floor(px * perPx);
48+ const i1 = Math.min(y.length, Math.max(i0 + 1, Math.floor((px + 1) * perPx)));
49+ let mn = Infinity, mx = -Infinity;
50+ for (let i = i0; i < i1; i++) {
51+ if (y[i] < mn) mn = y[i];
52+ if (y[i] > mx) mx = y[i];
53+ }
54+ if (!Number.isFinite(mn)) continue;
55+ ctx.moveTo(padL + px + 0.5, midY - mx * (ph / 2));
56+ ctx.lineTo(padL + px + 0.5, midY - mn * (ph / 2) + 0.5);
57+ }
58+ ctx.stroke();
59+ },
60+ [y, sampleRate],
61+ 170,
62+ );
63+ return <canvas ref={ref} />;
64+}
webapp/src/components/plots/canvasUtils.tsadded+66−0View file
@@ -0,0 +1,66 @@
1+import { useEffect, useRef } from 'react';
2+
3+// Canvas hook: resizes to element width at devicePixelRatio and calls draw.
4+export function useCanvas(
5+ draw: (ctx: CanvasRenderingContext2D, w: number, h: number) => void,
6+ deps: unknown[],
7+ height: number,
8+) {
9+ const ref = useRef<HTMLCanvasElement | null>(null);
10+
11+ useEffect(() => {
12+ const canvas = ref.current;
13+ if (!canvas) return;
14+ const render = () => {
15+ const dpr = window.devicePixelRatio || 1;
16+ const w = canvas.clientWidth;
17+ canvas.width = Math.max(1, Math.round(w * dpr));
18+ canvas.height = Math.max(1, Math.round(height * dpr));
19+ canvas.style.height = `${height}px`;
20+ const ctx = canvas.getContext('2d')!;
21+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
22+ ctx.clearRect(0, 0, w, height);
23+ draw(ctx, w, height);
24+ };
25+ render();
26+ const ro = new ResizeObserver(render);
27+ ro.observe(canvas);
28+ return () => ro.disconnect();
29+ // eslint-disable-next-line react-hooks/exhaustive-deps
30+ }, deps);
31+
32+ return ref;
33+}
34+
35+export const PLOT_COLORS = {
36+ bg: '#181b20',
37+ grid: '#2a2f38',
38+ axis: '#8b93a1',
39+ trace: '#4dd0e1',
40+ accent: '#e8a33d',
41+ good: '#6fbf73',
42+ warn: '#e57373',
43+};
44+
45+export function drawFrame(ctx: CanvasRenderingContext2D, w: number, h: number) {
46+ ctx.fillStyle = PLOT_COLORS.bg;
47+ ctx.fillRect(0, 0, w, h);
48+}
49+
50+// Viridis-like colormap for the spectrogram, t in [0,1].
51+export function colormap(t: number): [number, number, number] {
52+ const stops: [number, number, number][] = [
53+ [13, 8, 60], [84, 39, 143], [140, 41, 129], [196, 66, 90],
54+ [237, 105, 37], [251, 155, 6], [247, 209, 61], [252, 255, 164],
55+ ];
56+ const x = Math.min(0.9999, Math.max(0, t)) * (stops.length - 1);
57+ const i = Math.floor(x);
58+ const f = x - i;
59+ const a = stops[i];
60+ const b = stops[i + 1];
61+ return [
62+ Math.round(a[0] + f * (b[0] - a[0])),
63+ Math.round(a[1] + f * (b[1] - a[1])),
64+ Math.round(a[2] + f * (b[2] - a[2])),
65+ ];
66+}
webapp/src/hooks/useRenderer.tsadded+78−0View file
@@ -0,0 +1,78 @@
1+// Hook owning the render worker: debounced auto-preview on parameter change,
2+// explicit full-quality render, and a latest-wins request queue.
3+
4+import { useCallback, useEffect, useRef, useState } from 'react';
5+import type { WoodblockParams } from '../params';
6+import type { ProgressMessage, RenderResponse } from '../worker/renderWorker';
7+
8+export interface RenderState {
9+ result: RenderResponse | null;
10+ busy: boolean;
11+ stage: string;
12+ progress: number; // 0..1 of the current stage
13+}
14+
15+export function useRenderer(params: WoodblockParams) {
16+ const workerRef = useRef<Worker | null>(null);
17+ const nextId = useRef(1);
18+ const inFlight = useRef<number | null>(null);
19+ const queued = useRef<{ params: WoodblockParams; quality: 'preview' | 'full' } | null>(null);
20+ const [state, setState] = useState<RenderState>({
21+ result: null,
22+ busy: false,
23+ stage: '',
24+ progress: 0,
25+ });
26+
27+ useEffect(() => {
28+ const worker = new Worker(new URL('../worker/renderWorker.ts', import.meta.url), {
29+ type: 'module',
30+ });
31+ workerRef.current = worker;
32+ worker.onmessage = (ev: MessageEvent<RenderResponse | ProgressMessage>) => {
33+ const msg = ev.data;
34+ if (msg.type === 'progress') {
35+ setState((s) => ({ ...s, stage: msg.stage, progress: msg.frac }));
36+ return;
37+ }
38+ inFlight.current = null;
39+ const q = queued.current;
40+ if (q) {
41+ queued.current = null;
42+ // show the result we have, then immediately start the queued request
43+ setState({ result: msg, busy: true, stage: 'queued', progress: 0 });
44+ send(q.params, q.quality);
45+ } else {
46+ setState({ result: msg, busy: false, stage: '', progress: 0 });
47+ }
48+ };
49+ return () => {
50+ worker.terminate();
51+ workerRef.current = null;
52+ };
53+ // eslint-disable-next-line react-hooks/exhaustive-deps
54+ }, []);
55+
56+ const send = useCallback((p: WoodblockParams, quality: 'preview' | 'full') => {
57+ const worker = workerRef.current;
58+ if (!worker) return;
59+ if (inFlight.current !== null) {
60+ queued.current = { params: p, quality };
61+ return;
62+ }
63+ const id = nextId.current++;
64+ inFlight.current = id;
65+ setState((s) => ({ ...s, busy: true, stage: 'starting', progress: 0 }));
66+ worker.postMessage({ type: 'render', id, params: p, quality });
67+ }, []);
68+
69+ const renderFull = useCallback(() => send(params, 'full'), [params, send]);
70+
71+ // Debounced auto-preview whenever parameters change (including on mount).
72+ useEffect(() => {
73+ const t = setTimeout(() => send(params, 'preview'), 250);
74+ return () => clearTimeout(t);
75+ }, [params, send]);
76+
77+ return { ...state, renderFull };
78+}
webapp/src/lib/fft.tsadded+91−0View file
@@ -0,0 +1,91 @@
1+// Minimal iterative radix-2 complex FFT, in-place on separate re/im arrays.
2+
3+export function fft(re: Float64Array, im: Float64Array): void {
4+ const n = re.length;
5+ if ((n & (n - 1)) !== 0) throw new Error('FFT length must be a power of 2');
6+ // bit-reversal permutation
7+ for (let i = 1, j = 0; i < n; i++) {
8+ let bit = n >> 1;
9+ for (; j & bit; bit >>= 1) j ^= bit;
10+ j ^= bit;
11+ if (i < j) {
12+ const tr = re[i]; re[i] = re[j]; re[j] = tr;
13+ const ti = im[i]; im[i] = im[j]; im[j] = ti;
14+ }
15+ }
16+ for (let len = 2; len <= n; len <<= 1) {
17+ const ang = (-2 * Math.PI) / len;
18+ const wr = Math.cos(ang);
19+ const wi = Math.sin(ang);
20+ for (let i = 0; i < n; i += len) {
21+ let cwr = 1, cwi = 0;
22+ for (let k = 0; k < len / 2; k++) {
23+ const ur = re[i + k], ui = im[i + k];
24+ const vr = re[i + k + len / 2] * cwr - im[i + k + len / 2] * cwi;
25+ const vi = re[i + k + len / 2] * cwi + im[i + k + len / 2] * cwr;
26+ re[i + k] = ur + vr;
27+ im[i + k] = ui + vi;
28+ re[i + k + len / 2] = ur - vr;
29+ im[i + k + len / 2] = ui - vi;
30+ const nwr = cwr * wr - cwi * wi;
31+ cwi = cwr * wi + cwi * wr;
32+ cwr = nwr;
33+ }
34+ }
35+ }
36+}
37+
38+// Magnitude spectrum of a real signal with a Hann window, zero-padded to the
39+// next power of two (at least minSize). Returns bins 0..N/2.
40+export function magnitudeSpectrum(x: Float32Array | Float64Array, minSize = 16384): {
41+ mag: Float64Array;
42+ size: number;
43+} {
44+ let size = minSize;
45+ while (size < x.length) size <<= 1;
46+ const re = new Float64Array(size);
47+ const im = new Float64Array(size);
48+ const n = x.length;
49+ for (let i = 0; i < n; i++) {
50+ const w = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (n - 1)));
51+ re[i] = x[i] * w;
52+ }
53+ fft(re, im);
54+ const half = size / 2;
55+ const mag = new Float64Array(half + 1);
56+ for (let i = 0; i <= half; i++) mag[i] = Math.hypot(re[i], im[i]);
57+ return { mag, size };
58+}
59+
60+// Simple STFT magnitude in dB for the spectrogram, frames x bins.
61+export function stft(
62+ x: Float32Array,
63+ win = 512,
64+ hop = 128,
65+): { frames: number; bins: number; db: Float32Array; floorDb: number } {
66+ const nFrames = Math.max(1, Math.floor((x.length - win) / hop) + 1);
67+ const bins = win / 2 + 1;
68+ const db = new Float32Array(nFrames * bins);
69+ const re = new Float64Array(win);
70+ const im = new Float64Array(win);
71+ const hann = new Float64Array(win);
72+ for (let i = 0; i < win; i++) hann[i] = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (win - 1)));
73+ let maxDb = -Infinity;
74+ for (let f = 0; f < nFrames; f++) {
75+ const off = f * hop;
76+ for (let i = 0; i < win; i++) {
77+ re[i] = (off + i < x.length ? x[off + i] : 0) * hann[i];
78+ im[i] = 0;
79+ }
80+ fft(re, im);
81+ for (let b = 0; b < bins; b++) {
82+ const m = Math.hypot(re[b], im[b]);
83+ const v = 20 * Math.log10(m + 1e-12);
84+ db[f * bins + b] = v;
85+ if (v > maxDb) maxDb = v;
86+ }
87+ }
88+ // normalize so 0 dB is the frame-bin maximum
89+ for (let i = 0; i < db.length; i++) db[i] -= maxDb;
90+ return { frames: nFrames, bins, db, floorDb: -90 };
91+}
webapp/src/lib/wav.tsadded+40−0View file
@@ -0,0 +1,40 @@
1+// Mono PCM WAV encoding (16- or 24-bit), matching write_wav in woodblock.py.
2+
3+export function encodeWav(y: Float32Array, sampleRate: number, bitDepth: 16 | 24): Blob {
4+ const bytesPerSample = bitDepth / 8;
5+ const dataSize = y.length * bytesPerSample;
6+ const buf = new ArrayBuffer(44 + dataSize);
7+ const dv = new DataView(buf);
8+ const writeStr = (off: number, s: string) => {
9+ for (let i = 0; i < s.length; i++) dv.setUint8(off + i, s.charCodeAt(i));
10+ };
11+ writeStr(0, 'RIFF');
12+ dv.setUint32(4, 36 + dataSize, true);
13+ writeStr(8, 'WAVE');
14+ writeStr(12, 'fmt ');
15+ dv.setUint32(16, 16, true);
16+ dv.setUint16(20, 1, true); // PCM
17+ dv.setUint16(22, 1, true); // mono
18+ dv.setUint32(24, sampleRate, true);
19+ dv.setUint32(28, sampleRate * bytesPerSample, true);
20+ dv.setUint16(32, bytesPerSample, true);
21+ dv.setUint16(34, bitDepth, true);
22+ writeStr(36, 'data');
23+ dv.setUint32(40, dataSize, true);
24+
25+ let off = 44;
26+ for (let i = 0; i < y.length; i++) {
27+ const v = Math.max(-1, Math.min(1, y[i]));
28+ if (bitDepth === 16) {
29+ dv.setInt16(off, Math.round(v * 32767), true);
30+ off += 2;
31+ } else {
32+ const iv = Math.round(v * 8388607);
33+ dv.setUint8(off, iv & 0xff);
34+ dv.setUint8(off + 1, (iv >> 8) & 0xff);
35+ dv.setUint8(off + 2, (iv >> 16) & 0xff);
36+ off += 3;
37+ }
38+ }
39+ return new Blob([buf], { type: 'audio/wav' });
40+}
webapp/src/main.tsxadded+10−0View file
@@ -0,0 +1,10 @@
1+import React from 'react';
2+import { createRoot } from 'react-dom/client';
3+import App from './App';
4+import './style.css';
5+
6+createRoot(document.getElementById('root')!).render(
7+ <React.StrictMode>
8+ <App />
9+ </React.StrictMode>,
10+);
webapp/src/params.tsadded+259−0View file
@@ -0,0 +1,259 @@
1+// Parameter definitions for the woodblock model, mirroring WoodblockConfig in
2+// woodblock.py. Each parameter is assigned a tier that determines how much
3+// recomputation a change requires:
4+// solve -> full re-simulation (ODE solve) + radiation + output
5+// radiate -> radiation + output only (cached modal solve is reused)
6+// output -> output stage only (decimation product is reused)
7+
8+export interface WoodblockParams {
9+ // Body geometry
10+ length_m: number;
11+ width_m: number;
12+ wall_thickness_m: number;
13+ block_height_m: number;
14+ // Wood material
15+ density: number;
16+ youngs_modulus: number;
17+ hollow_softening: number;
18+ participating_mass_kg: number;
19+ inharmonicity: number;
20+ n_modes: number;
21+ max_mode_hz: number;
22+ // Damping
23+ loss_factor: number;
24+ loss_slope: number;
25+ // Mallet / strike
26+ strike_position: number;
27+ mallet_mass: number;
28+ mallet_velocity: number;
29+ mallet_stiffness: number;
30+ contact_exponent: number;
31+ contact_damping: number;
32+ // Air cavity
33+ cavity_volume_m3: number;
34+ slit_area_m2: number;
35+ slit_length_m: number;
36+ cavity_damping: number;
37+ cavity_coupling: number;
38+ cavity_span_lo: number;
39+ cavity_span_hi: number;
40+ // Acoustic field
41+ listener_x_m: number;
42+ listener_y_m: number;
43+ listener_z_m: number;
44+ slit_x_m: number;
45+ slit_y_m: number;
46+ slit_z_m: number;
47+ radiating_elements: number;
48+ speed_of_sound: number;
49+ air_density: number;
50+ // Rendering
51+ sample_rate: number;
52+ oversample: number;
53+ duration_s: number;
54+ pre_delay_ms: number;
55+ normalize: boolean;
56+ peak_dbfs: number;
57+ full_scale_pa: number;
58+ bit_depth: number;
59+}
60+
61+export const DEFAULTS: WoodblockParams = {
62+ length_m: 0.2,
63+ width_m: 0.055,
64+ wall_thickness_m: 0.01,
65+ block_height_m: 0.05,
66+ density: 850.0,
67+ youngs_modulus: 1.6e10,
68+ hollow_softening: 0.9,
69+ participating_mass_kg: 0.35,
70+ inharmonicity: 1.0,
71+ n_modes: 26,
72+ max_mode_hz: 20000.0,
73+ loss_factor: 0.01,
74+ loss_slope: 1.2e-7,
75+ strike_position: 0.42,
76+ mallet_mass: 0.025,
77+ mallet_velocity: 2.5,
78+ mallet_stiffness: 1.0e8,
79+ contact_exponent: 1.5,
80+ contact_damping: 0.8,
81+ cavity_volume_m3: 2.6e-4,
82+ slit_area_m2: 8.0e-4,
83+ slit_length_m: 0.012,
84+ cavity_damping: 0.06,
85+ cavity_coupling: 1.0,
86+ cavity_span_lo: 0.3,
87+ cavity_span_hi: 0.7,
88+ listener_x_m: 0.15,
89+ listener_y_m: 0.25,
90+ listener_z_m: 0.3,
91+ slit_x_m: 0.0,
92+ slit_y_m: 0.0275,
93+ slit_z_m: 0.0,
94+ radiating_elements: 48,
95+ speed_of_sound: 343.0,
96+ air_density: 1.204,
97+ sample_rate: 48000,
98+ oversample: 16,
99+ duration_s: 0.45,
100+ pre_delay_ms: 3.0,
101+ normalize: true,
102+ peak_dbfs: -3.0,
103+ full_scale_pa: 20.0,
104+ bit_depth: 24,
105+};
106+
107+export type Tier = 'solve' | 'radiate' | 'output';
108+
109+// Keys whose change invalidates the cached ODE solve. Note that
110+// block_height_m, the listener/slit positions, and radiating_elements only
111+// enter the radiation post-process, so they are deliberately absent here.
112+export const SOLVE_KEYS: (keyof WoodblockParams)[] = [
113+ 'length_m', 'width_m', 'wall_thickness_m',
114+ 'density', 'youngs_modulus', 'hollow_softening', 'participating_mass_kg',
115+ 'inharmonicity', 'n_modes', 'max_mode_hz',
116+ 'loss_factor', 'loss_slope',
117+ 'strike_position', 'mallet_mass', 'mallet_velocity', 'mallet_stiffness',
118+ 'contact_exponent', 'contact_damping',
119+ 'cavity_volume_m3', 'slit_area_m2', 'slit_length_m', 'cavity_damping',
120+ 'cavity_coupling', 'cavity_span_lo', 'cavity_span_hi',
121+ 'speed_of_sound', 'air_density', 'sample_rate', 'duration_s',
122+];
123+
124+export const RADIATE_KEYS: (keyof WoodblockParams)[] = [
125+ 'block_height_m',
126+ 'listener_x_m', 'listener_y_m', 'listener_z_m',
127+ 'slit_x_m', 'slit_y_m', 'slit_z_m',
128+];
129+
130+export const OUTPUT_KEYS: (keyof WoodblockParams)[] = [
131+ 'pre_delay_ms', 'normalize', 'peak_dbfs', 'full_scale_pa',
132+];
133+
134+export interface ParamSpec {
135+ key: keyof WoodblockParams;
136+ label: string;
137+ unit?: string;
138+ min: number;
139+ max: number;
140+ step?: number;
141+ log?: boolean;
142+ integer?: boolean;
143+ help?: string;
144+}
145+
146+export interface ParamGroup {
147+ name: string;
148+ specs: ParamSpec[];
149+}
150+
151+export const PARAM_GROUPS: ParamGroup[] = [
152+ {
153+ name: 'Body geometry',
154+ specs: [
155+ { key: 'length_m', label: 'Length', unit: 'm', min: 0.08, max: 0.5, step: 0.001, help: 'Long axis of the block; pitch scales as 1/L^2.' },
156+ { key: 'width_m', label: 'Width', unit: 'm', min: 0.02, max: 0.15, step: 0.001, help: 'Across the grain; sets source area and level.' },
157+ { key: 'wall_thickness_m', label: 'Wall thickness', unit: 'm', min: 0.003, max: 0.03, step: 0.0005, help: 'Bending depth of the vibrating top wall. Pitch scales linearly with this.' },
158+ { key: 'block_height_m', label: 'Block height', unit: 'm', min: 0.02, max: 0.12, step: 0.001, help: 'Outer height; sets the dipole arm between the antiphase faces. Radiation-only.' },
159+ ],
160+ },
161+ {
162+ name: 'Wood material',
163+ specs: [
164+ { key: 'density', label: 'Density', unit: 'kg/m³', min: 300, max: 1400, step: 10, help: 'Rosewood ~850, maple ~700. Pitch scales as 1/sqrt(rho).' },
165+ { key: 'youngs_modulus', label: "Young's modulus", unit: 'Pa', min: 2e9, max: 5e10, log: true, help: 'Along the grain. Pitch scales as sqrt(E).' },
166+ { key: 'hollow_softening', label: 'Hollow softening', min: 0.5, max: 1.0, step: 0.01, help: 'Stiffness removed by the cut slot; overall pitch trim.' },
167+ { key: 'participating_mass_kg', label: 'Participating mass', unit: 'kg', min: 0.05, max: 1.0, step: 0.01, help: 'Lumped modal mass; amplitude scales as 1/mass.' },
168+ { key: 'inharmonicity', label: 'Inharmonicity', min: 0.8, max: 1.3, step: 0.01, help: '>1 stretches the upper partials (more knock).' },
169+ ],
170+ },
171+ {
172+ name: 'Damping',
173+ specs: [
174+ { key: 'loss_factor', label: 'Loss factor η₀', min: 0.001, max: 0.08, step: 0.001, help: 'Frequency-independent wood loss; sets T60.' },
175+ { key: 'loss_slope', label: 'Loss slope η₁', unit: 's', min: 1e-8, max: 2e-6, log: true, help: 'Extra HF damping; the main brightness control.' },
176+ ],
177+ },
178+ {
179+ name: 'Mallet / strike',
180+ specs: [
181+ { key: 'strike_position', label: 'Strike position', min: 0.05, max: 0.95, step: 0.01, help: '0..1 along the bar. Avoid exactly 0.5 (node of every antisymmetric mode).' },
182+ { key: 'mallet_mass', label: 'Mallet mass', unit: 'kg', min: 0.005, max: 0.1, step: 0.001, help: 'Contact time and level.' },
183+ { key: 'mallet_velocity', label: 'Mallet velocity', unit: 'm/s', min: 0.2, max: 8.0, step: 0.1, help: 'Dynamics and timbre (the contact is nonlinear).' },
184+ { key: 'mallet_stiffness', label: 'Mallet stiffness', unit: 'N/m^α', min: 1e7, max: 3e10, log: true, help: 'Harder = shorter contact, brighter click.' },
185+ { key: 'contact_exponent', label: 'Contact exponent', min: 1.1, max: 3.0, step: 0.05, help: '1.5 Hertzian (wood/plastic), 2-3 for felt.' },
186+ { key: 'contact_damping', label: 'Contact damping λ', unit: 's/m', min: 0.0, max: 3.0, step: 0.05, help: 'Hunt-Crossley loss; damps rebound, softens attack.' },
187+ ],
188+ },
189+ {
190+ name: 'Air cavity',
191+ specs: [
192+ { key: 'cavity_volume_m3', label: 'Cavity volume', unit: 'm³', min: 2e-5, max: 2e-3, log: true, help: 'Helmholtz pitch.' },
193+ { key: 'slit_area_m2', label: 'Slit area', unit: 'm²', min: 5e-5, max: 5e-3, log: true, help: 'Helmholtz pitch and cavity output.' },
194+ { key: 'slit_length_m', label: 'Slit length', unit: 'm', min: 0.002, max: 0.05, step: 0.001, help: 'Wall thickness at the slot.' },
195+ { key: 'cavity_damping', label: 'Cavity damping ζ', min: 0.005, max: 0.5, step: 0.005, help: 'Hollow-resonance Q.' },
196+ { key: 'cavity_coupling', label: 'Cavity coupling', min: 0.0, max: 3.0, step: 0.05, help: '0 = dry bar, 1 = nominal, >1 = hollower.' },
197+ { key: 'cavity_span_lo', label: 'Cavity span low', min: 0.0, max: 0.5, step: 0.01, help: 'Fraction of the bar where the cavity begins.' },
198+ { key: 'cavity_span_hi', label: 'Cavity span high', min: 0.5, max: 1.0, step: 0.01, help: 'Fraction of the bar where the cavity ends.' },
199+ ],
200+ },
201+ {
202+ name: 'Listener / slit position',
203+ specs: [
204+ { key: 'listener_x_m', label: 'Listener x', unit: 'm', min: -1.0, max: 1.0, step: 0.005, help: 'Along the bar axis. Radiation-only: no re-solve needed.' },
205+ { key: 'listener_y_m', label: 'Listener y', unit: 'm', min: -1.0, max: 1.0, step: 0.005, help: 'Across the width.' },
206+ { key: 'listener_z_m', label: 'Listener z', unit: 'm', min: -1.0, max: 1.0, step: 0.005, help: 'Up, out of the struck face.' },
207+ { key: 'slit_x_m', label: 'Slit x', unit: 'm', min: -0.25, max: 0.25, step: 0.001, help: 'Slit opening position (monopole source point).' },
208+ { key: 'slit_y_m', label: 'Slit y', unit: 'm', min: -0.1, max: 0.1, step: 0.001 },
209+ { key: 'slit_z_m', label: 'Slit z', unit: 'm', min: -0.1, max: 0.1, step: 0.001 },
210+ ],
211+ },
212+ {
213+ name: 'Output',
214+ specs: [
215+ { key: 'duration_s', label: 'Duration', unit: 's', min: 0.1, max: 2.0, step: 0.05 },
216+ { key: 'pre_delay_ms', label: 'Pre-delay', unit: 'ms', min: 0, max: 50, step: 0.5 },
217+ { key: 'peak_dbfs', label: 'Peak level', unit: 'dBFS', min: -24, max: 0, step: 0.5, help: 'Used when normalize is on.' },
218+ { key: 'full_scale_pa', label: 'Full scale', unit: 'Pa', min: 1, max: 200, step: 1, help: 'Used when normalize is off: this many pascals maps to 0 dBFS.' },
219+ ],
220+ },
221+ {
222+ name: 'Advanced / quality',
223+ specs: [
224+ { key: 'n_modes', label: 'Max modes', min: 2, max: 40, step: 1, integer: true },
225+ { key: 'max_mode_hz', label: 'Max mode freq', unit: 'Hz', min: 2000, max: 24000, step: 500 },
226+ { key: 'radiating_elements', label: 'Radiating elements', min: 8, max: 128, step: 1, integer: true, help: 'Surface discretisation; accuracy, not timbre.' },
227+ { key: 'oversample', label: 'Oversample', min: 4, max: 32, step: 1, integer: true, help: 'Solver rate multiple; 16 for render, 8 for preview.' },
228+ { key: 'speed_of_sound', label: 'Speed of sound', unit: 'm/s', min: 300, max: 400, step: 1 },
229+ { key: 'air_density', label: 'Air density', unit: 'kg/m³', min: 0.9, max: 1.5, step: 0.01 },
230+ ],
231+ },
232+];
233+
234+const STORAGE_KEY = 'woodblock-webapp-params-v1';
235+
236+export function loadParams(): WoodblockParams {
237+ try {
238+ const raw = localStorage.getItem(STORAGE_KEY);
239+ if (!raw) return { ...DEFAULTS };
240+ const parsed = JSON.parse(raw);
241+ const out = { ...DEFAULTS };
242+ for (const k of Object.keys(DEFAULTS) as (keyof WoodblockParams)[]) {
243+ if (typeof parsed[k] === typeof DEFAULTS[k]) {
244+ (out as Record<string, unknown>)[k] = parsed[k];
245+ }
246+ }
247+ return out;
248+ } catch {
249+ return { ...DEFAULTS };
250+ }
251+}
252+
253+export function saveParams(p: WoodblockParams): void {
254+ try {
255+ localStorage.setItem(STORAGE_KEY, JSON.stringify(p));
256+ } catch {
257+ // storage unavailable; persistence is a convenience only
258+ }
259+}
webapp/src/physics/basis.tsadded+145−0View file
@@ -0,0 +1,145 @@
1+// Modal basis: free-free Euler-Bernoulli bar, mass-normalized bending modes.
2+// Port of ModalBasis in woodblock.py.
3+
4+import type { WoodblockParams } from '../params';
5+
6+// Roots of cos(bL)cosh(bL) = 1 for the free-free beam.
7+const BETA_L = [4.73004074, 7.85320462, 10.99560784, 14.13716549];
8+
9+export function betaL(n: number): number {
10+ if (n <= BETA_L.length) return BETA_L[n - 1];
11+ return ((2 * n + 1) * Math.PI) / 2.0;
12+}
13+
14+// Free-free mode shape on s in [0,1], unnormalized.
15+export function modeShape(bl: number, s: Float64Array): Float64Array {
16+ const sigma = (Math.cosh(bl) - Math.cos(bl)) / (Math.sinh(bl) - Math.sin(bl));
17+ const out = new Float64Array(s.length);
18+ for (let i = 0; i < s.length; i++) {
19+ const x = bl * s[i];
20+ out[i] = Math.cosh(x) + Math.cos(x) - sigma * (Math.sinh(x) + Math.sin(x));
21+ }
22+ return out;
23+}
24+
25+function trapz(y: Float64Array, ds: number): number {
26+ let sum = 0;
27+ for (let i = 0; i < y.length; i++) sum += y[i];
28+ sum -= 0.5 * (y[0] + y[y.length - 1]);
29+ return sum * ds;
30+}
31+
32+function interp1(x: number, xs: Float64Array, ys: Float64Array): number {
33+ if (x <= xs[0]) return ys[0];
34+ const n = xs.length;
35+ if (x >= xs[n - 1]) return ys[n - 1];
36+ // uniform grid
37+ const dx = xs[1] - xs[0];
38+ const j = Math.min(n - 2, Math.floor((x - xs[0]) / dx));
39+ const t = (x - xs[j]) / dx;
40+ return ys[j] * (1 - t) + ys[j + 1] * t;
41+}
42+
43+export interface ModalBasisData {
44+ n: number;
45+ f: Float64Array; // modal frequencies (Hz)
46+ omega: Float64Array;
47+ zeta: Float64Array;
48+ phiC: Float64Array; // mode shape at strike point (1/sqrt(kg))
49+ A: Float64Array; // cavity volume coupling (m^2/sqrt(kg))
50+ mass: number;
51+ // Mass-normalized shapes sampled on the dense s grid, for later
52+ // interpolation at arbitrary radiating-element positions.
53+ sGrid: Float64Array;
54+ shapes: Float64Array[]; // length n, each sampled on sGrid
55+}
56+
57+const NS = 4096;
58+
59+export function buildBasis(cfg: WoodblockParams): ModalBasisData {
60+ const L = cfg.length_m;
61+ const W = cfg.width_m;
62+ const H = cfg.wall_thickness_m;
63+ const area = W * H;
64+ const mass = cfg.participating_mass_kg || cfg.density * area * L;
65+ const cLong = Math.sqrt(cfg.youngs_modulus / cfg.density);
66+
67+ // f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)), I = W H^3/12
68+ const prefactor = (cLong * H) / (2.0 * Math.PI * Math.sqrt(12.0) * L * L);
69+
70+ const s = new Float64Array(NS);
71+ for (let i = 0; i < NS; i++) s[i] = i / (NS - 1);
72+ const ds = 1.0 / (NS - 1);
73+
74+ const freqs: number[] = [];
75+ const phiC: number[] = [];
76+ const A: number[] = [];
77+ const shapes: Float64Array[] = [];
78+ const lo = cfg.cavity_span_lo;
79+ const hi = cfg.cavity_span_hi;
80+
81+ for (let n = 1; n <= cfg.n_modes; n++) {
82+ const bl = betaL(n);
83+ let f = prefactor * bl * bl * cfg.hollow_softening;
84+ if (n > 1 && cfg.inharmonicity !== 1.0) {
85+ const f0 = prefactor * betaL(1) ** 2 * cfg.hollow_softening;
86+ f = f0 * (f / f0) ** cfg.inharmonicity;
87+ }
88+ if (f > cfg.max_mode_hz) break;
89+
90+ const raw = modeShape(bl, s);
91+ const sq = new Float64Array(NS);
92+ for (let i = 0; i < NS; i++) sq[i] = raw[i] * raw[i];
93+ const rms = Math.sqrt(trapz(sq, ds));
94+ const norm = 1.0 / (rms * Math.sqrt(mass));
95+ const shape = new Float64Array(NS);
96+ for (let i = 0; i < NS; i++) shape[i] = raw[i] * norm;
97+
98+ freqs.push(f);
99+ phiC.push(interp1(cfg.strike_position, s, shape));
100+ const masked = new Float64Array(NS);
101+ for (let i = 0; i < NS; i++) masked[i] = s[i] >= lo && s[i] <= hi ? shape[i] : 0.0;
102+ A.push(W * L * trapz(masked, ds));
103+ shapes.push(shape);
104+ }
105+
106+ const nKept = freqs.length;
107+ const fArr = Float64Array.from(freqs);
108+ const omega = new Float64Array(nKept);
109+ const zeta = new Float64Array(nKept);
110+ const AArr = new Float64Array(nKept);
111+ for (let i = 0; i < nKept; i++) {
112+ omega[i] = 2 * Math.PI * fArr[i];
113+ zeta[i] = 0.5 * (cfg.loss_factor + cfg.loss_slope * omega[i]);
114+ AArr[i] = A[i] * cfg.cavity_coupling;
115+ }
116+
117+ return {
118+ n: nKept,
119+ f: fArr,
120+ omega,
121+ zeta,
122+ phiC: Float64Array.from(phiC),
123+ A: AArr,
124+ mass,
125+ sGrid: s,
126+ shapes,
127+ };
128+}
129+
130+// Interpolate the stored mode shapes at M radiating-element centres.
131+// Returns phiElem as an array of n Float64Arrays of length M, plus sElem.
132+export function elementShapes(basis: ModalBasisData, M: number): {
133+ phiElem: Float64Array[];
134+ sElem: Float64Array;
135+} {
136+ const sElem = new Float64Array(M);
137+ for (let j = 0; j < M; j++) sElem[j] = (j + 0.5) / M;
138+ const phiElem: Float64Array[] = [];
139+ for (let n = 0; n < basis.n; n++) {
140+ const row = new Float64Array(M);
141+ for (let j = 0; j < M; j++) row[j] = interp1(sElem[j], basis.sGrid, basis.shapes[n]);
142+ phiElem.push(row);
143+ }
144+ return { phiElem, sElem };
145+}
webapp/src/physics/radiate.tsadded+105−0View file
@@ -0,0 +1,105 @@
1+// Radiation: retarded free-field pressure at the listener point.
2+// Port of radiate() and check_setup() in woodblock.py. The top and bottom
3+// faces move together, so each strip is a local dipole; the surface velocity
4+// sum is computed once per strip and reused for both faces with opposite sign.
5+
6+import type { WoodblockParams } from '../params';
7+import type { ModalBasisData } from './basis';
8+import { elementShapes } from './basis';
9+
10+export interface RadiateResult {
11+ p: Float64Array; // pascals at the listener, at fsHi
12+ warnings: string[];
13+}
14+
15+export function radiate(
16+ cfg: WoodblockParams,
17+ Qdd: Float64Array[],
18+ Xdd: Float64Array,
19+ basis: ModalBasisData,
20+ fsHi: number,
21+ M: number,
22+ onProgress?: (frac: number) => void,
23+): RadiateResult {
24+ const rho0 = cfg.air_density;
25+ const c = cfg.speed_of_sound;
26+ const N = basis.n;
27+ const T = Xdd.length;
28+ const L = cfg.length_m;
29+ const W = cfg.width_m;
30+ const halfH = 0.5 * cfg.block_height_m;
31+ const patch = (W * L) / M;
32+
33+ const { phiElem, sElem } = elementShapes(basis, M);
34+
35+ const lx = cfg.listener_x_m;
36+ const ly = cfg.listener_y_m;
37+ const lz = cfg.listener_z_m;
38+
39+ const distTop = new Float64Array(M);
40+ const distBot = new Float64Array(M);
41+ let maxDist = 0;
42+ for (let j = 0; j < M; j++) {
43+ const x = (sElem[j] - 0.5) * L;
44+ const dx = x - lx;
45+ const dy = 0.0 - ly;
46+ const dzT = halfH - lz;
47+ const dzB = -halfH - lz;
48+ distTop[j] = Math.sqrt(dx * dx + dy * dy + dzT * dzT);
49+ distBot[j] = Math.sqrt(dx * dx + dy * dy + dzB * dzB);
50+ maxDist = Math.max(maxDist, distTop[j], distBot[j]);
51+ }
52+ const sx = cfg.slit_x_m - lx;
53+ const sy = cfg.slit_y_m - ly;
54+ const sz = cfg.slit_z_m - lz;
55+ const slitDist = Math.sqrt(sx * sx + sy * sy + sz * sz);
56+ maxDist = Math.max(maxDist, slitDist);
57+
58+ const maxDelay = Math.ceil((maxDist / c) * fsHi) + 2;
59+ const out = new Float64Array(T + maxDelay);
60+ const wdd = new Float64Array(T);
61+
62+ for (let j = 0; j < M; j++) {
63+ // Surface acceleration of strip j, shared by the +z and -z sources
64+ wdd.fill(0);
65+ for (let n = 0; n < N; n++) {
66+ const phi = phiElem[n][j];
67+ if (phi === 0) continue;
68+ const q = Qdd[n];
69+ for (let t = 0; t < T; t++) wdd[t] += phi * q[t];
70+ }
71+ const ampTop = (patch * rho0) / (4.0 * Math.PI * distTop[j]);
72+ const ampBot = -(patch * rho0) / (4.0 * Math.PI * distBot[j]);
73+ const kTop = Math.round((distTop[j] / c) * fsHi);
74+ const kBot = Math.round((distBot[j] / c) * fsHi);
75+ for (let t = 0; t < T; t++) {
76+ out[kTop + t] += ampTop * wdd[t];
77+ out[kBot + t] += ampBot * wdd[t];
78+ }
79+ if (onProgress) onProgress((j + 1) / M);
80+ }
81+
82+ // Slit monopole
83+ const ampSlit = (cfg.slit_area_m2 * rho0) / (4.0 * Math.PI * slitDist);
84+ const kSlit = Math.round((slitDist / c) * fsHi);
85+ for (let t = 0; t < T; t++) out[kSlit + t] += ampSlit * Xdd[t];
86+
87+ // Setup warnings (check_setup in the reference)
88+ const warnings: string[] = [];
89+ const spacing = L / M;
90+ const fMax = basis.f[basis.n - 1];
91+ if (spacing > (0.25 * c) / fMax) {
92+ const need = Math.ceil(L / ((0.25 * c) / fMax));
93+ warnings.push(
94+ `element spacing is coarse at ${fMax.toFixed(0)} Hz; consider ${need} radiating elements`,
95+ );
96+ }
97+ const d = Math.sqrt(lx * lx + ly * ly + lz * lz);
98+ if (d < 0.5 * L) {
99+ warnings.push(
100+ `listener is ${(d * 100).toFixed(0)} cm from the block centre, inside the source region; point-source elements are crude here`,
101+ );
102+ }
103+
104+ return { p: out.subarray(0, T).slice(), warnings };
105+}
webapp/src/physics/render.tsadded+114−0View file
@@ -0,0 +1,114 @@
1+// Rendering orchestration: solve -> radiate -> output stage. Mirrors render()
2+// in woodblock.py, split into stages so the worker can cache intermediates.
3+
4+import type { WoodblockParams } from '../params';
5+import { buildBasis, type ModalBasisData } from './basis';
6+import { simulate, type SolveResult } from './solver';
7+import { radiate } from './radiate';
8+import { decimate } from './resample';
9+
10+export interface Diagnostics {
11+ modesKept: number;
12+ partialsHz: number[];
13+ ratios: number[];
14+ helmholtzHz: number;
15+ contactMs: number;
16+ peakForceN: number;
17+ peakPa: number;
18+ peakSplDb: number;
19+ listenerDistM: number;
20+ warnings: string[];
21+}
22+
23+export interface OutputStageResult {
24+ y: Float64Array;
25+ peakPa: number;
26+ peakSplDb: number;
27+ clipWarning: string | null;
28+}
29+
30+// Steps 3-5 of §8: decimate, pre-roll, fade, normalize or absolute scale.
31+export function outputStage(cfg: WoodblockParams, p: Float64Array, oversample: number): OutputStageResult {
32+ let peakPa = 0;
33+ for (let i = 0; i < p.length; i++) {
34+ const a = Math.abs(p[i]);
35+ if (a > peakPa) peakPa = a;
36+ }
37+ const peakSplDb = 20.0 * Math.log10(peakPa / 20e-6 + 1e-30);
38+
39+ let y = decimate(p, oversample);
40+
41+ const pad = Math.floor((cfg.pre_delay_ms * 1e-3) * cfg.sample_rate);
42+ const total = Math.floor(cfg.duration_s * cfg.sample_rate);
43+ if (pad > 0) {
44+ const padded = new Float64Array(Math.min(pad + y.length, total));
45+ padded.set(y.subarray(0, padded.length - pad), pad);
46+ y = padded;
47+ } else if (y.length > total) {
48+ y = y.subarray(0, total).slice();
49+ }
50+
51+ const fade = Math.min(256, Math.floor(y.length / 8));
52+ if (fade > 1) {
53+ for (let i = 0; i < fade; i++) {
54+ const g = 1.0 - i / (fade - 1);
55+ y[y.length - fade + i] *= g * g;
56+ }
57+ }
58+
59+ let clipWarning: string | null = null;
60+ if (cfg.normalize) {
61+ let peak = 0;
62+ for (let i = 0; i < y.length; i++) peak = Math.max(peak, Math.abs(y[i]));
63+ if (peak > 0) {
64+ const g = (10 ** (cfg.peak_dbfs / 20.0)) / peak;
65+ for (let i = 0; i < y.length; i++) y[i] *= g;
66+ }
67+ } else {
68+ let peak = 0;
69+ for (let i = 0; i < y.length; i++) {
70+ y[i] /= cfg.full_scale_pa;
71+ peak = Math.max(peak, Math.abs(y[i]));
72+ }
73+ if (peak > 1.0) {
74+ clipWarning = `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${cfg.full_scale_pa.toFixed(2)}`;
75+ }
76+ }
77+
78+ return { y, peakPa, peakSplDb, clipWarning };
79+}
80+
81+export interface FullRenderResult {
82+ y: Float64Array;
83+ basis: ModalBasisData;
84+ solve: SolveResult;
85+ diag: Diagnostics;
86+}
87+
88+// One-shot render, used by the self-test; the worker calls the stages
89+// individually so it can cache them.
90+export function renderAll(cfg: WoodblockParams, oversample?: number, elements?: number): FullRenderResult {
91+ const os = oversample ?? cfg.oversample;
92+ const M = elements ?? cfg.radiating_elements;
93+ const basis = buildBasis(cfg);
94+ const solve = simulate(cfg, basis, os);
95+ const rad = radiate(cfg, solve.Qdd, solve.Xdd, basis, solve.fsHi, M);
96+ const out = outputStage(cfg, rad.p, os);
97+
98+ const lx = cfg.listener_x_m, ly = cfg.listener_y_m, lz = cfg.listener_z_m;
99+ const warnings = [...rad.warnings];
100+ if (out.clipWarning) warnings.push(out.clipWarning);
101+ const diag: Diagnostics = {
102+ modesKept: basis.n,
103+ partialsHz: Array.from(basis.f),
104+ ratios: Array.from(basis.f).map((f) => f / basis.f[0]),
105+ helmholtzHz: solve.helmholtzHz,
106+ contactMs: solve.contactMs,
107+ peakForceN: solve.peakForceN,
108+ peakPa: out.peakPa,
109+ peakSplDb: out.peakSplDb,
110+ listenerDistM: Math.sqrt(lx * lx + ly * ly + lz * lz),
111+ warnings,
112+ };
113+ return { y: out.y, basis, solve, diag };
114+}
webapp/src/physics/resample.tsadded+56−0View file
@@ -0,0 +1,56 @@
1+// Anti-aliased decimation: Kaiser-windowed-sinc FIR lowpass at ~0.45 * fs_lo,
2+// evaluated only at the retained output samples (polyphase, so the cost is
3+// taps * T / R rather than taps * T). Replaces scipy.resample_poly(p, 1, R).
4+
5+function besselI0(x: number): number {
6+ // Series expansion; converges quickly for the argument range used here.
7+ let sum = 1.0;
8+ let term = 1.0;
9+ const half = x / 2.0;
10+ for (let k = 1; k < 64; k++) {
11+ term *= (half / k) * (half / k);
12+ sum += term;
13+ if (term < 1e-16 * sum) break;
14+ }
15+ return sum;
16+}
17+
18+export function designDecimator(R: number, beta = 8.0): Float64Array {
19+ // Cutoff at 0.45 of the output Nyquist-setting rate: fc/fs_hi = 0.45/R.
20+ const taps = 12 * R + 1; // odd length, linear phase
21+ const h = new Float64Array(taps);
22+ const M2 = (taps - 1) / 2;
23+ const fcNorm = 0.45 / R; // cycles per input sample
24+ const i0beta = besselI0(beta);
25+ let sum = 0;
26+ for (let i = 0; i < taps; i++) {
27+ const m = i - M2;
28+ const sinc = m === 0 ? 2 * fcNorm : Math.sin(2 * Math.PI * fcNorm * m) / (Math.PI * m);
29+ const r = m / M2;
30+ const w = besselI0(beta * Math.sqrt(Math.max(0, 1 - r * r))) / i0beta;
31+ h[i] = sinc * w;
32+ sum += h[i];
33+ }
34+ for (let i = 0; i < taps; i++) h[i] /= sum; // unity DC gain
35+ return h;
36+}
37+
38+export function decimate(p: Float64Array, R: number): Float64Array {
39+ if (R === 1) return p.slice();
40+ const h = designDecimator(R);
41+ const taps = h.length;
42+ const M2 = (taps - 1) / 2;
43+ const T = p.length;
44+ const outLen = Math.floor(T / R);
45+ const y = new Float64Array(outLen);
46+ for (let m = 0; m < outLen; m++) {
47+ const centre = m * R;
48+ let acc = 0;
49+ // y[m] = sum_i h[i] * p[centre + M2 - i], zero outside the buffer
50+ const iLo = Math.max(0, centre + M2 - (T - 1));
51+ const iHi = Math.min(taps - 1, centre + M2);
52+ for (let i = iLo; i <= iHi; i++) acc += h[i] * p[centre + M2 - i];
53+ y[m] = acc;
54+ }
55+ return y;
56+}
webapp/src/physics/selftest.tsadded+42−0View file
@@ -0,0 +1,42 @@
1+// Validation against the Python reference (woodblock.py at defaults):
2+// modes kept : 6
3+// partials Hz: 1003, 2766, 5423, 8964, 13390, 18702
4+// helmholtz : 484 Hz
5+// contact : 0.301 ms, peak 794 N
6+// peak press.: 6.132 Pa (109.7 dB SPL)
7+// Run with: npm run selftest
8+
9+import { DEFAULTS } from '../params';
10+import { renderAll } from './render';
11+
12+// The app's default mallet stiffness is softer than woodblock.py's; pin the
13+// reference value here so the comparison against the Python output holds.
14+const REFERENCE = { ...DEFAULTS, mallet_stiffness: 3.0e8 };
15+
16+const t0 = performance.now();
17+const { y, basis, diag } = renderAll(REFERENCE);
18+const t1 = performance.now();
19+
20+console.log(`modes kept : ${basis.n}`);
21+console.log('partials Hz: ' + Array.from(basis.f).map((f) => f.toFixed(0)).join(', '));
22+console.log('ratios : ' + Array.from(basis.f).map((f) => (f / basis.f[0]).toFixed(2)).join(', '));
23+console.log(`helmholtz : ${diag.helmholtzHz.toFixed(0)} Hz`);
24+console.log(`contact : ${diag.contactMs.toFixed(3)} ms, peak ${diag.peakForceN.toFixed(0)} N`);
25+console.log(`peak press.: ${diag.peakPa.toFixed(3)} Pa (${diag.peakSplDb.toFixed(1)} dB SPL)`);
26+console.log(`length : ${(y.length / REFERENCE.sample_rate).toFixed(3)} s @ ${REFERENCE.sample_rate} Hz`);
27+for (const w of diag.warnings) console.log(`WARNING : ${w}`);
28+console.log(`render time: ${((t1 - t0) / 1000).toFixed(2)} s`);
29+
30+// Validation listener positions from the spec (§7), peak SPL targets.
31+const targets: [number, number, number, number][] = [
32+ [0.15, 0.25, 0.30, 109.7],
33+ [0.00, 0.00, 0.40, 100.9],
34+ [0.00, 0.40, 0.00, 90.9],
35+ [0.00, 0.05, 0.02, 125.4],
36+];
37+for (const [x, yy, z, target] of targets) {
38+ const cfg = { ...REFERENCE, listener_x_m: x, listener_y_m: yy, listener_z_m: z };
39+ const r = renderAll(cfg);
40+ const ok = Math.abs(r.diag.peakSplDb - target) < 0.5 ? 'ok' : 'MISMATCH';
41+ console.log(`listener (${x}, ${yy}, ${z}): ${r.diag.peakSplDb.toFixed(1)} dB SPL (target ${target}) ${ok}`);
42+}
webapp/src/physics/solver.tsadded+157−0View file
@@ -0,0 +1,157 @@
1+// Time-domain RK4 integration of the coupled mallet / body / cavity system.
2+// Port of simulate() in woodblock.py. State: [q(N), qdot(N), xi, xidot, y, ydot].
3+
4+import type { WoodblockParams } from '../params';
5+import type { ModalBasisData } from './basis';
6+
7+export interface SolveResult {
8+ Qdd: Float64Array[]; // N arrays of length T: modal accelerations
9+ Xdd: Float64Array; // slit-plug acceleration
10+ force: Float64Array; // contact force f_c(t)
11+ fsHi: number;
12+ helmholtzHz: number;
13+ peakForceN: number;
14+ contactMs: number;
15+}
16+
17+export function simulate(
18+ cfg: WoodblockParams,
19+ basis: ModalBasisData,
20+ oversample: number,
21+ onProgress?: (frac: number) => void,
22+): SolveResult {
23+ const N = basis.n;
24+ if (N === 0) throw new Error('No modes below max_mode_hz -- check geometry/material.');
25+
26+ const fsHi = cfg.sample_rate * oversample;
27+ const dt = 1.0 / fsHi;
28+ const nSteps = Math.floor(cfg.duration_s * fsHi);
29+
30+ const om2 = new Float64Array(N);
31+ const twoZetaOm = new Float64Array(N);
32+ for (let i = 0; i < N; i++) {
33+ om2[i] = basis.omega[i] * basis.omega[i];
34+ twoZetaOm[i] = 2.0 * basis.zeta[i] * basis.omega[i];
35+ }
36+ const phiC = basis.phiC;
37+ const A = basis.A;
38+
39+ const rho0 = cfg.air_density;
40+ const c = cfg.speed_of_sound;
41+ const V = cfg.cavity_volume_m3;
42+ const S = cfg.slit_area_m2;
43+ // End correction on the slit: ~0.85 * equivalent radius per open end.
44+ const Leff = cfg.slit_length_m + 2.0 * 0.85 * Math.sqrt(S / Math.PI);
45+ const omegaH = c * Math.sqrt(S / (V * Leff));
46+ const twoZetaHOm = 2.0 * cfg.cavity_damping * omegaH;
47+ const omegaH2 = omegaH * omegaH;
48+ const pScale = (rho0 * c * c) / V;
49+ const cavDrive = (c * c) / (V * Leff);
50+
51+ const K = cfg.mallet_stiffness;
52+ const alpha = cfg.contact_exponent;
53+ const mu = cfg.contact_damping;
54+ const invMh = 1.0 / cfg.mallet_mass;
55+
56+ const dim = 2 * N + 4;
57+ const IX = 2 * N;
58+ const IXD = 2 * N + 1;
59+ const IY = 2 * N + 2;
60+ const IYD = 2 * N + 3;
61+
62+ const state = new Float64Array(dim);
63+ state[IYD] = cfg.mallet_velocity;
64+
65+ let lastFc = 0.0;
66+ let lastXidd = 0.0;
67+
68+ const qddScratch = new Float64Array(N);
69+
70+ function deriv(s: Float64Array, d: Float64Array): void {
71+ // Mallet contact
72+ let w = 0.0;
73+ let wd = 0.0;
74+ for (let i = 0; i < N; i++) {
75+ w += phiC[i] * s[i];
76+ wd += phiC[i] * s[N + i];
77+ }
78+ const y = s[IY];
79+ const yd = s[IYD];
80+ const delta = y - w;
81+ let fc = 0.0;
82+ if (delta > 0.0) {
83+ // Hunt-Crossley: f = K d^a (1 + lambda * ddot); lambda is in s/m so the
84+ // correction stays dimensionless and O(1) over a sub-ms collision.
85+ fc = K * Math.pow(delta, alpha) * (1.0 + mu * (yd - wd));
86+ if (fc < 0.0) fc = 0.0;
87+ }
88+ lastFc = fc;
89+
90+ // Cavity pressure from combined volume change of walls and slit plug
91+ let Aq = 0.0;
92+ for (let i = 0; i < N; i++) Aq += A[i] * s[i];
93+ const pc = pScale * (S * s[IX] + Aq);
94+
95+ for (let i = 0; i < N; i++) {
96+ const qdd = -twoZetaOm[i] * s[N + i] - om2[i] * s[i] + fc * phiC[i] - A[i] * pc;
97+ d[i] = s[N + i];
98+ d[N + i] = qdd;
99+ qddScratch[i] = qdd;
100+ }
101+ const xidd = -twoZetaHOm * s[IXD] - omegaH2 * s[IX] - cavDrive * Aq;
102+ d[IX] = s[IXD];
103+ d[IXD] = xidd;
104+ d[IY] = yd;
105+ d[IYD] = -fc * invMh;
106+ lastXidd = xidd;
107+ }
108+
109+ const Qdd: Float64Array[] = [];
110+ for (let i = 0; i < N; i++) Qdd.push(new Float64Array(nSteps));
111+ const Xdd = new Float64Array(nSteps);
112+ const force = new Float64Array(nSteps);
113+
114+ const k1 = new Float64Array(dim);
115+ const k2 = new Float64Array(dim);
116+ const k3 = new Float64Array(dim);
117+ const k4 = new Float64Array(dim);
118+ const tmp = new Float64Array(dim);
119+
120+ const progressEvery = Math.max(1, Math.floor(nSteps / 50));
121+
122+ for (let k = 0; k < nSteps; k++) {
123+ deriv(state, k1);
124+ for (let i = 0; i < N; i++) Qdd[i][k] = qddScratch[i];
125+ Xdd[k] = lastXidd;
126+ force[k] = lastFc;
127+
128+ for (let i = 0; i < dim; i++) tmp[i] = state[i] + 0.5 * dt * k1[i];
129+ deriv(tmp, k2);
130+ for (let i = 0; i < dim; i++) tmp[i] = state[i] + 0.5 * dt * k2[i];
131+ deriv(tmp, k3);
132+ for (let i = 0; i < dim; i++) tmp[i] = state[i] + dt * k3[i];
133+ deriv(tmp, k4);
134+ for (let i = 0; i < dim; i++) {
135+ state[i] += (dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
136+ }
137+
138+ if (onProgress && k % progressEvery === 0) onProgress(k / nSteps);
139+ }
140+
141+ let peakForce = 0.0;
142+ let contactSamples = 0;
143+ for (let k = 0; k < nSteps; k++) {
144+ if (force[k] > peakForce) peakForce = force[k];
145+ if (force[k] > 0.0) contactSamples++;
146+ }
147+
148+ return {
149+ Qdd,
150+ Xdd,
151+ force,
152+ fsHi,
153+ helmholtzHz: omegaH / (2 * Math.PI),
154+ peakForceN: peakForce,
155+ contactMs: (contactSamples * dt) * 1e3,
156+ };
157+}
webapp/src/style.cssadded+235−0View file
@@ -0,0 +1,235 @@
1+:root {
2+ --bg: #14161a;
3+ --panel: #1d2026;
4+ --panel2: #23272f;
5+ --border: #31363f;
6+ --text: #d6dae2;
7+ --muted: #8b93a1;
8+ --accent: #e8a33d;
9+ --accent2: #4dd0e1;
10+ --good: #6fbf73;
11+ --warn: #e57373;
12+ font-family: system-ui, -apple-system, 'Segoe UI', Roboto, sans-serif;
13+}
14+
15+* { box-sizing: border-box; }
16+
17+html, body, #root {
18+ margin: 0;
19+ height: 100%;
20+ background: var(--bg);
21+ color: var(--text);
22+}
23+
24+.app {
25+ display: flex;
26+ flex-direction: column;
27+ height: 100vh;
28+}
29+
30+.topbar {
31+ display: flex;
32+ align-items: center;
33+ gap: 10px;
34+ padding: 8px 14px;
35+ background: var(--panel);
36+ border-bottom: 1px solid var(--border);
37+ flex-wrap: wrap;
38+}
39+
40+.topbar h1 {
41+ font-size: 16px;
42+ margin: 0 12px 0 0;
43+ font-weight: 600;
44+ color: var(--accent);
45+ white-space: nowrap;
46+}
47+
48+button {
49+ background: var(--panel2);
50+ color: var(--text);
51+ border: 1px solid var(--border);
52+ border-radius: 6px;
53+ padding: 6px 14px;
54+ font-size: 13px;
55+ cursor: pointer;
56+}
57+button:hover { border-color: var(--muted); }
58+button:disabled { opacity: 0.45; cursor: default; }
59+button.primary {
60+ background: var(--accent);
61+ color: #1a1204;
62+ border-color: var(--accent);
63+ font-weight: 600;
64+}
65+button.play {
66+ background: var(--good);
67+ color: #06210a;
68+ border-color: var(--good);
69+ font-weight: 600;
70+}
71+button.play.preview {
72+ background: transparent;
73+ color: var(--accent);
74+ border-color: var(--accent);
75+}
76+button.loop.active {
77+ background: var(--accent2);
78+ color: #062a2f;
79+ border-color: var(--accent2);
80+ font-weight: 600;
81+}
82+
83+.status {
84+ margin-left: auto;
85+ font-size: 12px;
86+ color: var(--muted);
87+ display: flex;
88+ align-items: center;
89+ gap: 8px;
90+ min-width: 180px;
91+ justify-content: flex-end;
92+}
93+
94+.badge {
95+ font-size: 11px;
96+ border-radius: 10px;
97+ padding: 2px 8px;
98+ border: 1px solid var(--border);
99+}
100+.badge.preview { color: var(--accent); border-color: var(--accent); }
101+.badge.full { color: var(--good); border-color: var(--good); }
102+
103+.progressbar {
104+ width: 90px;
105+ height: 6px;
106+ background: var(--panel2);
107+ border-radius: 3px;
108+ overflow: hidden;
109+}
110+.progressbar > div {
111+ height: 100%;
112+ background: var(--accent2);
113+ transition: width 0.1s linear;
114+}
115+
116+.main {
117+ display: flex;
118+ flex: 1;
119+ min-height: 0;
120+}
121+
122+.sidebar {
123+ width: 340px;
124+ min-width: 340px;
125+ overflow-y: auto;
126+ background: var(--panel);
127+ border-right: 1px solid var(--border);
128+ padding: 8px;
129+}
130+
131+.content {
132+ flex: 1;
133+ overflow-y: auto;
134+ padding: 10px;
135+ display: grid;
136+ grid-template-columns: 1fr 1fr;
137+ gap: 10px;
138+ align-content: start;
139+}
140+
141+.card {
142+ background: var(--panel);
143+ border: 1px solid var(--border);
144+ border-radius: 8px;
145+ padding: 8px 10px 10px;
146+ min-width: 0;
147+}
148+.card.wide { grid-column: 1 / -1; }
149+.card h2 {
150+ font-size: 12px;
151+ text-transform: uppercase;
152+ letter-spacing: 0.06em;
153+ color: var(--muted);
154+ margin: 0 0 6px;
155+ font-weight: 600;
156+}
157+.card canvas { display: block; width: 100%; }
158+
159+.scene3d { height: 380px; position: relative; }
160+.scene3d canvas { height: 100%; border-radius: 6px; }
161+.scene-hint {
162+ position: absolute;
163+ bottom: 10px;
164+ right: 14px;
165+ font-size: 11px;
166+ color: var(--muted);
167+ pointer-events: none;
168+}
169+
170+/* Parameter panel */
171+details.group {
172+ border: 1px solid var(--border);
173+ border-radius: 6px;
174+ margin-bottom: 6px;
175+ background: var(--panel2);
176+}
177+details.group summary {
178+ cursor: pointer;
179+ padding: 6px 10px;
180+ font-size: 13px;
181+ font-weight: 600;
182+ user-select: none;
183+}
184+details.group[open] summary { border-bottom: 1px solid var(--border); }
185+.group-body { padding: 6px 10px 8px; }
186+
187+.param-row { margin-bottom: 8px; }
188+.param-label {
189+ display: flex;
190+ justify-content: space-between;
191+ font-size: 12px;
192+ margin-bottom: 2px;
193+ align-items: baseline;
194+}
195+.param-label .name { color: var(--text); }
196+.param-controls { display: flex; gap: 8px; align-items: center; }
197+.param-controls input[type='range'] { flex: 1; accent-color: var(--accent2); }
198+.param-controls input[type='number'] {
199+ width: 88px;
200+ background: var(--bg);
201+ color: var(--text);
202+ border: 1px solid var(--border);
203+ border-radius: 4px;
204+ padding: 3px 5px;
205+ font-size: 12px;
206+}
207+.param-controls select, .param-controls input[type='checkbox'] {
208+ background: var(--bg);
209+ color: var(--text);
210+ border: 1px solid var(--border);
211+ border-radius: 4px;
212+ font-size: 12px;
213+}
214+.param-help { font-size: 11px; color: var(--muted); margin-top: 1px; }
215+
216+/* Diagnostics */
217+.diag table { border-collapse: collapse; font-size: 12px; width: 100%; }
218+.diag td, .diag th {
219+ padding: 2px 8px 2px 0;
220+ text-align: left;
221+ color: var(--text);
222+ font-weight: 400;
223+}
224+.diag th { color: var(--muted); white-space: nowrap; }
225+.diag .partials { font-variant-numeric: tabular-nums; }
226+.warning { color: var(--warn); font-size: 12px; margin-top: 4px; }
227+.cachenote { color: var(--muted); font-size: 11px; margin-top: 6px; }
228+
229+.mode-list { display: flex; gap: 6px; flex-wrap: wrap; margin-bottom: 6px; }
230+.mode-list button { padding: 3px 10px; font-size: 12px; }
231+.mode-list button.active { border-color: var(--accent2); color: var(--accent2); }
232+
233+@media (max-width: 1100px) {
234+ .content { grid-template-columns: 1fr; }
235+}
webapp/src/worker/renderWorker.tsadded+201−0View file
@@ -0,0 +1,201 @@
1+// Web Worker that runs the simulation off the main thread and caches
2+// intermediates so cheap parameter changes skip the expensive stages:
3+// solve tier -> ODE solve (cached on solve-tier params + oversample)
4+// radiate tier -> retarded-source sum (cached on radiate-tier params + M)
5+// output tier -> decimation + normalization (always cheap, always re-run)
6+
7+import { SOLVE_KEYS, RADIATE_KEYS, type WoodblockParams } from '../params';
8+import { buildBasis, type ModalBasisData } from '../physics/basis';
9+import { simulate, type SolveResult } from '../physics/solver';
10+import { radiate } from '../physics/radiate';
11+import { outputStage } from '../physics/render';
12+
13+export interface RenderRequest {
14+ type: 'render';
15+ id: number;
16+ params: WoodblockParams;
17+ quality: 'preview' | 'full';
18+}
19+
20+export interface RenderResponse {
21+ type: 'result';
22+ id: number;
23+ quality: 'preview' | 'full';
24+ y: Float32Array; // final signal, [-1, 1]
25+ sampleRate: number;
26+ force: Float32Array; // contact force, first forceMs at fsHi
27+ forceFsHz: number;
28+ modeShapes: Float32Array; // n x SHAPE_PTS, concatenated
29+ shapePts: number;
30+ diag: {
31+ modesKept: number;
32+ partialsHz: number[];
33+ ratios: number[];
34+ helmholtzHz: number;
35+ contactMs: number;
36+ peakForceN: number;
37+ peakPa: number;
38+ peakSplDb: number;
39+ listenerDistM: number;
40+ warnings: string[];
41+ solveMs: number; // 0 when the cached solve was reused
42+ radiateMs: number; // 0 when the cached radiation was reused
43+ solveCached: boolean;
44+ radiateCached: boolean;
45+ oversample: number;
46+ elements: number;
47+ };
48+}
49+
50+export interface ProgressMessage {
51+ type: 'progress';
52+ id: number;
53+ stage: 'solve' | 'radiate' | 'output';
54+ frac: number;
55+}
56+
57+const SHAPE_PTS = 129;
58+const FORCE_MS = 8;
59+
60+// Small LRU caches so preview- and full-quality intermediates coexist. A
61+// radiate-tier change right after a full render then reuses the full solve,
62+// and the "preview" is silently upgraded to full quality.
63+const solveCache = new Map<string, { basis: ModalBasisData; solve: SolveResult }>();
64+const radiateCache = new Map<string, { p: Float64Array; warnings: string[] }>();
65+const SOLVE_CACHE_MAX = 2;
66+const RADIATE_CACHE_MAX = 3;
67+
68+function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
69+ const v = cache.get(key);
70+ if (v !== undefined) {
71+ cache.delete(key);
72+ cache.set(key, v); // refresh recency
73+ }
74+ return v;
75+}
76+
77+function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
78+ cache.delete(key);
79+ cache.set(key, value);
80+ while (cache.size > max) {
81+ const oldest = cache.keys().next().value as string;
82+ cache.delete(oldest);
83+ }
84+}
85+
86+function keyFor(params: WoodblockParams, keys: readonly (keyof WoodblockParams)[]): string {
87+ return keys.map((k) => String(params[k])).join('|');
88+}
89+
90+self.onmessage = (ev: MessageEvent<RenderRequest>) => {
91+ const msg = ev.data;
92+ if (msg.type !== 'render') return;
93+ const { id, params } = msg;
94+ let quality = msg.quality;
95+
96+ const solveParamsKey = keyFor(params, SOLVE_KEYS);
97+ const fullSolveKey = solveParamsKey + `|os=${params.oversample}`;
98+ if (quality === 'preview' && solveCache.has(fullSolveKey)) {
99+ quality = 'full'; // full solve already cached; only cheap stages remain
100+ }
101+
102+ const oversample = quality === 'preview' ? Math.min(8, params.oversample) : params.oversample;
103+ const elements = quality === 'preview'
104+ ? Math.min(24, params.radiating_elements)
105+ : params.radiating_elements;
106+
107+ const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
108+
109+ const solveKey = solveParamsKey + `|os=${oversample}`;
110+ let solveMs = 0;
111+ let solveCached = true;
112+ let solveEntry = lruGet(solveCache, solveKey);
113+ if (!solveEntry) {
114+ solveCached = false;
115+ const t0 = performance.now();
116+ const basis = buildBasis(params);
117+ const solve = simulate(params, basis, oversample, (frac) =>
118+ post({ type: 'progress', id, stage: 'solve', frac }),
119+ );
120+ solveMs = performance.now() - t0;
121+ solveEntry = { basis, solve };
122+ lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
123+ }
124+ const { basis, solve } = solveEntry;
125+
126+ const radiateKey = solveKey + '||' + keyFor(params, RADIATE_KEYS) + `|M=${elements}`;
127+ let radiateMs = 0;
128+ let radiateCached = true;
129+ let radEntry = lruGet(radiateCache, radiateKey);
130+ if (!radEntry) {
131+ radiateCached = false;
132+ const t0 = performance.now();
133+ const rad = radiate(params, solve.Qdd, solve.Xdd, basis, solve.fsHi, elements, (frac) =>
134+ post({ type: 'progress', id, stage: 'radiate', frac }),
135+ );
136+ radiateMs = performance.now() - t0;
137+ radEntry = { p: rad.p, warnings: rad.warnings };
138+ lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
139+ }
140+
141+ const out = outputStage(params, radEntry.p, oversample);
142+
143+ // Contact force excerpt for plotting
144+ const forceLen = Math.min(solve.force.length, Math.round((FORCE_MS / 1000) * solve.fsHi));
145+ const force = new Float32Array(forceLen);
146+ for (let i = 0; i < forceLen; i++) force[i] = solve.force[i];
147+
148+ // Mode shapes decimated for the mode viewer (unnormalized visual scale)
149+ const modeShapes = new Float32Array(basis.n * SHAPE_PTS);
150+ for (let n = 0; n < basis.n; n++) {
151+ const src = basis.shapes[n];
152+ let peak = 0;
153+ for (let i = 0; i < src.length; i++) peak = Math.max(peak, Math.abs(src[i]));
154+ for (let i = 0; i < SHAPE_PTS; i++) {
155+ const j = Math.round((i / (SHAPE_PTS - 1)) * (src.length - 1));
156+ modeShapes[n * SHAPE_PTS + i] = peak > 0 ? src[j] / peak : 0;
157+ }
158+ }
159+
160+ const y32 = new Float32Array(out.y.length);
161+ for (let i = 0; i < out.y.length; i++) y32[i] = out.y[i];
162+
163+ const warnings = [...radEntry.warnings];
164+ if (out.clipWarning) warnings.push(out.clipWarning);
165+ const lx = params.listener_x_m, ly = params.listener_y_m, lz = params.listener_z_m;
166+
167+ const resp: RenderResponse = {
168+ type: 'result',
169+ id,
170+ quality,
171+ y: y32,
172+ sampleRate: params.sample_rate,
173+ force,
174+ forceFsHz: solve.fsHi,
175+ modeShapes,
176+ shapePts: SHAPE_PTS,
177+ diag: {
178+ modesKept: basis.n,
179+ partialsHz: Array.from(basis.f),
180+ ratios: Array.from(basis.f).map((f) => f / basis.f[0]),
181+ helmholtzHz: solve.helmholtzHz,
182+ contactMs: solve.contactMs,
183+ peakForceN: solve.peakForceN,
184+ peakPa: out.peakPa,
185+ peakSplDb: out.peakSplDb,
186+ listenerDistM: Math.sqrt(lx * lx + ly * ly + lz * lz),
187+ warnings,
188+ solveMs,
189+ radiateMs,
190+ solveCached,
191+ radiateCached,
192+ oversample,
193+ elements,
194+ },
195+ };
196+ (self as unknown as Worker).postMessage(resp, [
197+ y32.buffer,
198+ force.buffer,
199+ modeShapes.buffer,
200+ ]);
201+};
webapp/tsconfig.jsonadded+19−0View file
@@ -0,0 +1,19 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "useDefineForClassFields": true,
5+ "lib": ["ES2022", "DOM", "DOM.Iterable", "WebWorker"],
6+ "module": "ESNext",
7+ "skipLibCheck": true,
8+ "moduleResolution": "bundler",
9+ "resolveJsonModule": true,
10+ "isolatedModules": true,
11+ "noEmit": true,
12+ "jsx": "react-jsx",
13+ "strict": true,
14+ "noUnusedLocals": true,
15+ "noUnusedParameters": true,
16+ "noFallthroughCasesInSwitch": true
17+ },
18+ "include": ["src"]
19+}
webapp/vite.config.tsadded+8−0View file
@@ -0,0 +1,8 @@
1+import { defineConfig } from 'vite';
2+import react from '@vitejs/plugin-react';
3+
4+export default defineConfig({
5+ plugins: [react()],
6+ base: './',
7+ worker: { format: 'es' },
8+});
woodblock.mdadded+349−0View file
@@ -0,0 +1,349 @@
1+# Woodblock physical model — specification
2+
3+Companion to `woodblock.py`. This describes the physics, the equations, the state
4+vector, and the traps. The Python script is the reference implementation; where
5+this document and the script disagree, the script is right.
6+
7+The model produces **free-field acoustic pressure in pascals at a point in
8+space**. There is no microphone, no directivity, no gain staging.
9+
10+---
11+
12+## 1. Signal chain
13+
14+```
15+mallet (Hunt–Crossley contact)
16+ → body bending modes (free–free bar, modal reduction) ⇄ two-way coupled
17+ → air cavity (Helmholtz resonator) ⇄
18+ → distributed surface dipole sources + slit monopole
19+ → retarded free-field pressure at listener point
20+ → anti-aliased decimation → PCM
21+```
22+
23+Three physical sub-systems (body, cavity, mallet) integrate together as one
24+coupled ODE. Radiation is a post-process: it is linear and non-reacting, so it
25+runs once after the solve, not inside the loop.
26+
27+---
28+
29+## 2. Coordinate system and geometry
30+
31+Block-centred, metres. Needed for the 3D view and for radiation.
32+
33+| Axis | Direction |
34+|---|---|
35+| +x | along the bar's long axis |
36+| +y | across the width |
37+| +z | up, out of the struck face |
38+
39+- Block occupies `x ∈ [−L/2, +L/2]`, `y ∈ [−W/2, +W/2]`, `z ∈ [−H/2, +H/2]`
40+ where `L = length_m`, `W = width_m`, `H = block_height_m`.
41+- **Radiating faces** sit at `z = ±H/2`.
42+- **Strike point** is at `x_c = (strike_position − ½)·L` on the top face.
43+- **Slit opening** is a point at `(slit_x_m, slit_y_m, slit_z_m)`; default is on
44+ the side wall at `(0, +W/2, 0)`.
45+- **Listener** is at `(listener_x_m, listener_y_m, listener_z_m)`.
46+- `wall_thickness_m` (`h`) is the *bending depth* of the vibrating top wall. It
47+ is NOT the block height. It sets pitch; `block_height_m` sets the dipole arm.
48+
49+For the 3D view: draw the box, the slot, a marker at the strike point, a marker
50+at the slit, a marker at the listener, and optionally the sampled surface
51+elements along the top face.
52+
53+---
54+
55+## 3. Body: free–free Euler–Bernoulli bar
56+
57+### 3.1 Modal frequencies
58+
59+Eigenvalues `β_n L` are the roots of `cos(βL)·cosh(βL) = 1`:
60+
61+```
62+4.73004074, 7.85320462, 10.99560784, 14.13716549, then (2n+1)π/2
63+```
64+
65+The asymptotic form is exact to ~1e-9 from n = 5 up. With `c_L = √(E/ρ)`,
66+`I = W h³/12`, `A = W h`:
67+
68+```
69+f_n = κ · (β_n L)² / (2π L²) · √(E I / (ρ A))
70+ = κ · (β_n L)² · c_L · h / (2π √12 · L²)
71+```
72+
73+`κ = hollow_softening` accounts for stiffness removed by the slot.
74+
75+Pitch scales **linearly with h** and as **√(E/ρ)**. Verified: h = 6/10/16 mm
76+gives 602/1003/1606 Hz; E = 6e9 vs 1.6e10 gives 614 vs 1003 Hz.
77+
78+Optional stretch (`γ = inharmonicity`), applied to partials above the first:
79+
80+```
81+f_n ← f_1 · (f_n / f_1)^γ
82+```
83+
84+Discard modes above `max_mode_hz` (20 kHz). Defaults keep 6.
85+
86+### 3.2 Mode shapes
87+
88+On `s = x/L ∈ [0,1]`, with `σ_n = (cosh βL − cos βL)/(sinh βL − sin βL)`:
89+
90+```
91+φ̃_n(s) = cosh(β_n L s) + cos(β_n L s) − σ_n [ sinh(β_n L s) + sin(β_n L s) ]
92+```
93+
94+Mass-normalise so the modal mass is 1:
95+
96+```
97+φ_n(s) = φ̃_n(s) / ( √(∫₀¹ φ̃_n² ds) · √M_p )
98+```
99+
100+`M_p = participating_mass_kg`. Units: `φ_n` is `kg^−1/2`, modal coordinate `q_n`
101+is `m·√kg`, physical displacement `w = Σ q_n φ_n` is metres.
102+
103+> `M_p` is a **lumped correction, the least principled parameter in the model.**
104+> The bar formula otherwise locks frequency and modal mass together, and the
105+> thin-wall mass (`ρWhL` ≈ 93 g) makes the block far too compliant — 395 µm of
106+> deflection and 119 dB SPL at 42 cm. A real woodblock is 300–500 g with the
107+> whole box participating. `M_p = 0.35` gives 124 µm and 109.7 dB. It sets
108+> absolute level, so treat SPL as plausible, not calibrated.
109+
110+### 3.3 Damping
111+
112+```
113+ζ_n = ½ ( η₀ + η₁ ω_n )
114+```
115+
116+`η₀ = loss_factor` (frequency-independent wood loss), `η₁ = loss_slope` (extra
117+HF damping). Bigger = duller and shorter. `η₀`: 0.010 → T60 185 ms;
118+0.030 → 66 ms.
119+
120+---
121+
122+## 4. Mallet contact (the only nonlinearity)
123+
124+Surface displacement and velocity at the strike point:
125+
126+```
127+w = Σ_n q_n φ_n(s_c)
128+ẇ = Σ_n q̇_n φ_n(s_c)
129+δ = y − w (compression; y = mallet position)
130+```
131+
132+**Hunt–Crossley**, with `λ = contact_damping` in **s/m**:
133+
134+```
135+f_c = K δ^α (1 + λ (ẏ − ẇ)) if δ > 0, clamped to ≥ 0
136+f_c = 0 otherwise
137+
138+m_h ÿ = −f_c
139+```
140+
141+Initial conditions: `y = 0`, `ẏ = v₀ = mallet_velocity`, all else zero.
142+
143+> **Trap.** Do not write the loss term as `μ · d(δ^α)/dt`. That makes `μ` a
144+> relaxation *time*; at μ = 0.35 s against a 0.1 ms collision, a 3 µm
145+> compression produced **1.9 × 10⁷ N**. Every strike became a delta function and
146+> mallet stiffness had literally zero effect on the spectrum. The form above
147+> keeps the correction dimensionless and O(1).
148+
149+Expected behaviour, worth asserting in tests: contact duration scales roughly as
150+`K^−2/5`. `K` = 3e7 / 3e8 / 3e9 / 3e10 → 0.92 / 0.59 / 0.28 / 0.16 ms and
151+287 / 615 / 1289 / 3012 N. Because the contact is nonlinear, timbre varies with
152+strike velocity on its own — that is physics, not a scripted velocity layer.
153+
154+---
155+
156+## 5. Cavity: Helmholtz resonator
157+
158+Effective plug length with end corrections on both ends:
159+
160+```
161+L_eff = slit_length_m + 2 · 0.85 · √(S/π) S = slit_area_m2
162+ω_H = c √( S / (V L_eff) ) V = cavity_volume_m3
163+```
164+
165+Coupling coefficients — net volume each mode injects into the cavity, over the
166+wall span `[lo, hi] = cavity_wall_span`:
167+
168+```
169+A_n = χ · W L ∫_lo^hi φ_n(s) ds χ = cavity_coupling
170+```
171+
172+Cavity pressure, from wall motion plus slit-plug displacement `ξ`:
173+
174+```
175+p_c = (ρ₀ c² / V) ( S ξ + Σ_n A_n q_n )
176+```
177+
178+Plug equation (mass per unit area `ρ₀ L_eff`, driven by `−p_c`):
179+
180+```
181+ξ̈ = −2 ζ_H ω_H ξ̇ − ω_H² ξ − (c² / (V L_eff)) Σ_n A_n q_n
182+```
183+
184+`A_n` appears both in the cavity drive and in the modal load below. **That
185+symmetry is the correctness check** — if the coupling only goes one way, the
186+energy will not be conserved.
187+
188+---
189+
190+## 6. Coupled system
191+
192+State vector, length `2N + 4`:
193+
194+```
195+[ q₁..q_N , q̇₁..q̇_N , ξ , ξ̇ , y , ẏ ]
196+```
197+
198+Derivatives:
199+
200+```
201+q̈_n = −2 ζ_n ω_n q̇_n − ω_n² q_n + f_c φ_n(s_c) − A_n p_c
202+ξ̈ = −2 ζ_H ω_H ξ̇ − ω_H² ξ − (c²/(V L_eff)) Σ_m A_m q_m
203+ÿ = −f_c / m_h
204+```
205+
206+Integrate with **RK4** at `f_hi = sample_rate × oversample`. Record `q̈_n(t)` for
207+all modes and `ξ̈(t)` — those are the radiation inputs.
208+
209+The system is linear apart from `f_c`. Eigenvalues of the coupled linear part are
210+all strictly decaying (slowest ≈ −34 s⁻¹ at defaults); a zero or positive real
211+part means the coupling was assembled wrongly.
212+
213+**Convergence.** Relative to `oversample = 32`: `oversample = 8` differs by 5.2%,
214+`16` by 1.8%. Use 16 for render, 8 is acceptable for interactive preview. The
215+residual sits in the top mode near 18.7 kHz.
216+
217+---
218+
219+## 7. Radiation: pressure at a point
220+
221+Each cross-section of a bending bar **translates**, so the top and bottom faces
222+inject equal and opposite volume. Every slice is a local dipole, and the net
223+volume displacement of the bar is exactly zero.
224+
225+> **Trap.** This is why you cannot use a single lumped monopole term with
226+> `B_n = ∫φ_n ds` for the body: elastic modes are orthogonal to rigid
227+> translation, so that integral is **identically zero** and the body contributes
228+> nothing. Near-field emphasis and HF directivity must come from path lengths and
229+> delays, not from a gain constant.
230+
231+Discretise the top face into `M = radiating_elements` strips at
232+`s_j = (j − ½)/M`, `x_j = (s_j − ½)L`, each of area `patch = W L / M`. Every
233+strip contributes a pair of point sources:
234+
235+```
236+r_j⁺ = (x_j, 0, +H/2) sign +1
237+r_j⁻ = (x_j, 0, −H/2) sign −1
238+```
239+
240+with volume acceleration
241+
242+```
243+V̈_j(t) = ± patch · Σ_n φ_n(s_j) q̈_n(t)
244+```
245+
246+Plus one monopole at the slit: `V̈_slit(t) = S ξ̈(t)`.
247+
248+Sum retarded monopoles over all `2M + 1` sources:
249+
250+```
251+p(r, t) = Σ_j ρ₀ / (4π d_j) · V̈_j( t − d_j / c ) d_j = |r − r_j|
252+```
253+
254+Integer-sample delays are fine: at `f_hi` = 768 kHz one sample is 0.45 mm.
255+
256+**Element count.** Spacing must satisfy `L/M < ¼ · c/f_max`; at defaults that
257+needs M ≥ 44, so M = 48. Warn the user if `M` is too coarse, and if the listener
258+is inside the source region (`|r| < L/2`) where point-source strips are crude.
259+
260+### Validation — the app should reproduce these
261+
262+| Listener (x, y, z) m | Peak | Note |
263+|---|---|---|
264+| 0.15, 0.25, 0.30 | 109.7 dB | default, oblique |
265+| 0.00, 0.00, 0.40 | 100.9 dB | above centre; 2766 Hz partial −27 dB (symmetry null) |
266+| 0.00, 0.40, 0.00 | 90.9 dB | **edge-on**: faces equidistant, body radiation cancels — 2766 Hz −61 dB, 8964 Hz −90 dB, leaving near-pure cavity monopole |
267+| 0.00, 0.05, 0.02 | 125.4 dB | close to slit, cavity-dominated |
268+
269+Falloff is −7.8 dB per doubling close in, relaxing to −6.5 dB by 1.6 m: the
270+dipole-to-inverse-square transition. If your port shows a flat −6 dB everywhere,
271+the two faces are not in antiphase.
272+
273+---
274+
275+## 8. Sampling and output
276+
277+1. Solve at `f_hi = 48000 × oversample`.
278+2. Radiate → `p(t)` in pascals at `f_hi`.
279+3. **Anti-alias and decimate** by `oversample`. Python uses `scipy.resample_poly`;
280+ in the browser implement a windowed-sinc (Kaiser, β ≈ 8) FIR lowpass at
281+ ~0.45·48 kHz with roughly `12 × oversample` taps, then take every
282+ `oversample`-th sample. Do not decimate without the filter — the double
283+ differentiation puts real energy near `f_hi/2`.
284+4. Prepend `pre_delay_ms` of silence, apply a short squared fade at the tail.
285+5. Either normalise to `peak_dbfs`, or divide by `full_scale_pa` for absolute
286+ pressure. Report peak Pa and dB SPL (`20 log₁₀(p/20 µPa)`), and warn on clip.
287+
288+> **Trap.** Do not "remove DC" by subtracting the mean of the first N samples.
289+> The strike is at t = 0, so the head of the buffer is the loudest part of the
290+> file; subtracting it injects a large offset across the whole signal. The model
291+> has no DC — `p.mean()` is ~1e-11 — so no DC removal is needed at all.
292+
293+---
294+
295+## 9. Parameters
296+
297+| Parameter | Default | Effect |
298+|---|---|---|
299+| `length_m` | 0.200 | pitch ∝ 1/L² |
300+| `width_m` | 0.055 | source area / level |
301+| `wall_thickness_m` | 0.010 | **pitch ∝ h** |
302+| `block_height_m` | 0.050 | dipole arm; bigger = more low end |
303+| `density` | 850 | pitch ∝ 1/√ρ |
304+| `youngs_modulus` | 1.6e10 | pitch ∝ √E |
305+| `hollow_softening` | 0.90 | overall pitch trim |
306+| `participating_mass_kg` | 0.35 | level ∝ 1/M (lumped; see §3.2) |
307+| `inharmonicity` | 1.00 | >1 stretches partials (2766 → 3061 Hz at 1.10) |
308+| `loss_factor` | 0.010 | T60 (185 ms → 66 ms at 0.030) |
309+| `loss_slope` | 1.2e-7 | HF damping; main brightness control |
310+| `strike_position` | 0.42 | partial balance. **Avoid 0.5** — a node of every antisymmetric mode |
311+| `mallet_mass` | 0.025 | contact time, level |
312+| `mallet_velocity` | 2.5 | dynamics *and* timbre (nonlinear contact) |
313+| `mallet_stiffness` | 3.0e8 | hardness; centroid 1953 → 5550 Hz at 3e10 |
314+| `contact_exponent` | 1.5 | 1.5 Hertzian, 2–3 felt |
315+| `contact_damping` | 0.8 | Hunt–Crossley λ, **s/m not seconds** |
316+| `cavity_volume_m3` | 2.6e-4 | Helmholtz pitch |
317+| `slit_area_m2` | 8.0e-4 | Helmholtz pitch, cavity output |
318+| `slit_length_m` | 0.012 | Helmholtz pitch |
319+| `cavity_damping` | 0.06 | hollow-resonance Q |
320+| `cavity_coupling` | 1.0 | 0 = dry bar, >1 = hollower |
321+| `listener_*_m` | .15/.25/.30 | see §7 |
322+| `radiating_elements` | 48 | accuracy, not timbre (see §7) |
323+| `oversample` | 16 | accuracy vs speed |
324+
325+Defaults yield: partials 1003 / 2766 / 5423 / 8964 / 13390 / 18702 Hz (ratios
326+1.00, 2.76, 5.40, 8.93 — bar inharmonicity, **not** octaves), Helmholtz 484 Hz,
327+contact 0.30 ms at 794 N, T60 ≈ 149 ms, 109.7 dB SPL.
328+
329+---
330+
331+## 10. Browser implementation notes
332+
333+- **Cost.** Solve is `4 × T` derivative evaluations with N ≈ 6 (T = 345,600 at
334+ defaults) — cheap. Radiation is the bottleneck: `2M+1` sources × N × T ≈ 200 M
335+ multiply-adds. Compute `Σ_n φ_n(s_j) q̈_n(t)` **once per strip** and reuse it
336+ for the `+z` and `−z` sources with opposite sign — that halves the work.
337+- Run in a **Web Worker**, transfer the result as a `Float32Array`
338+ `ArrayBuffer`, and copy into an `AudioBuffer` at 48 kHz. Keep the solver in
339+ `Float64Array`.
340+- For live parameter dragging, preview at `oversample = 8` and
341+ `radiating_elements = 24`, then re-render at full quality on release.
342+- Nothing needs real-time streaming — one strike is one offline render of well
343+ under a second.
344+- Cache the modal basis; it only depends on geometry, material,
345+ `strike_position`, and `radiating_elements`, not on mallet or listener
346+ parameters. Moving the listener requires only re-running §7.
347+- Good things to visualise: mode shapes animated as `φ̃_n(s)` along the bar;
348+ partial frequencies as a spectrum overlay; the contact-force pulse `f_c(t)`
349+ (its width sets brightness); source strips and listener in the 3D view.
woodblock.pyadded+502−0View file
@@ -0,0 +1,502 @@
1+#!/usr/bin/env python3
2+"""
3+woodblock.py -- physically-based simulation of a single woodblock strike,
4+sampled as acoustic pressure at a point in space.
5+
6+Model chain:
7+
8+ mallet (Hunt-Crossley contact)
9+ -> body modes (free-free bar, modal reduction)
10+ <-> air cavity (Helmholtz resonator, two-way coupled)
11+ -> distributed surface sources + slit monopole
12+ -> free-field pressure at the listener point
13+ -> anti-aliased decimation -> .wav
14+
15+There is no microphone model: no directivity, no capsule response, no
16+gain staging. The output is p(r, t) in pascals, sampled at one point,
17+optionally normalized for convenience when writing the file.
18+
19+Everything interesting is in WoodblockConfig. Any field can be overridden
20+from the command line, e.g.
21+
22+ python woodblock.py --listener_z_m 0.05 --strike_position 0.30
23+ python woodblock.py --normalize 0 --full_scale_pa 50 -o absolute.wav
24+
25+Requires: numpy, scipy
26+"""
27+
28+from __future__ import annotations
29+
30+import argparse
31+import dataclasses
32+import math
33+import wave
34+from dataclasses import dataclass
35+
36+import numpy as np
37+from scipy.signal import resample_poly
38+
39+
40+# ----------------------------------------------------------------------------
41+# Configuration
42+# ----------------------------------------------------------------------------
43+
44+@dataclass
45+class WoodblockConfig:
46+ # --- Body geometry -------------------------------------------------------
47+ length_m: float = 0.200 # long axis of the block
48+ width_m: float = 0.055 # across the grain
49+ wall_thickness_m: float = 0.010 # the hollowing means the *top wall* is the
50+ # vibrating element, not the whole block.
51+ # Pitch scales linearly with this.
52+ block_height_m: float = 0.050 # outer height. Sets the spacing between the
53+ # antiphase top/bottom faces, i.e. the
54+ # dipole arm -- bigger = more low end.
55+
56+ # --- Wood material -------------------------------------------------------
57+ density: float = 850.0 # kg/m^3 (rosewood ~850, maple ~700)
58+ youngs_modulus: float = 1.6e10 # Pa, along the grain
59+ # Frequencies scale as sqrt(E/rho), so stiffer/lighter = brighter + higher.
60+
61+ hollow_softening: float = 0.90 # 0.5-1.0; the cut slot removes stiffness
62+ participating_mass_kg: float = 0.35 # modal mass used for normalization.
63+ # The bar formula sets the *frequencies*,
64+ # but a real box vibrates with far more
65+ # mass than an isolated thin wall, and the
66+ # two are otherwise locked together. This
67+ # is a lumped stand-in for that; set 0 to
68+ # fall back to the thin-wall beam mass.
69+ # Scales amplitude as 1/mass.
70+ inharmonicity: float = 1.00 # exponent on partial ratios. >1 stretches
71+ # the upper partials (more "knock")
72+ n_modes: int = 26
73+ max_mode_hz: float = 20000.0 # discard modes above this
74+
75+ # --- Damping: zeta_n = 0.5 * (eta0 + eta1 * omega_n) ---------------------
76+ loss_factor: float = 0.010 # eta0, frequency-independent wood loss
77+ loss_slope: float = 1.2e-7 # eta1, extra damping for high partials
78+ # (bigger = duller, faster decay)
79+
80+ # --- Mallet / strike -----------------------------------------------------
81+ strike_position: float = 0.42 # 0..1 along the bar. Avoid exactly 0.5:
82+ # dead centre is a node of every
83+ # antisymmetric mode and sounds too pure
84+ mallet_mass: float = 0.025 # kg
85+ mallet_velocity: float = 2.5 # m/s at first contact -- dynamics
86+ mallet_stiffness: float = 3.0e8 # K in f = K*delta^alpha. Bigger = harder,
87+ # shorter contact, brighter click.
88+ # ~1e7 soft rubber ... ~1e10 wood on wood
89+ contact_exponent: float = 1.5 # 1.5 Hertzian (wood/plastic), 2-3 for felt
90+ contact_damping: float = 0.8 # Hunt-Crossley lambda, in s/m (NOT a
91+ # time). 0 = lossless contact; 0.5-2 is
92+ # typical. Damps rebound, softens attack
93+
94+ # --- Air cavity (Helmholtz) ---------------------------------------------
95+ cavity_volume_m3: float = 2.6e-4 # enclosed air volume
96+ slit_area_m2: float = 8.0e-4 # opening cross-section
97+ slit_length_m: float = 0.012 # wall thickness at the slot
98+ cavity_damping: float = 0.06 # zeta_H; viscous + radiation losses
99+ cavity_coupling: float = 1.0 # scales body <-> cavity coupling.
100+ # 0 = dry bar, 1 = nominal, >1 = hollower
101+ cavity_wall_span: tuple = (0.30, 0.70) # fraction of bar facing the cavity
102+
103+ # --- Acoustic field ------------------------------------------------------
104+ # Block-centred coordinates, in metres: +x along the bar's long axis,
105+ # +y across the width, +z up out of the struck face. The block spans
106+ # x in [-length/2, +length/2] and its faces sit at z = +-block_height/2.
107+ listener_x_m: float = 0.15
108+ listener_y_m: float = 0.25
109+ listener_z_m: float = 0.30
110+ # The slit opening, treated as a monopole at this point.
111+ slit_x_m: float = 0.0
112+ slit_y_m: float = 0.0275
113+ slit_z_m: float = 0.0
114+ radiating_elements: int = 48 # surface discretisation along the bar.
115+ # Must be fine enough that element
116+ # spacing << wavelength (see check_setup)
117+ speed_of_sound: float = 343.0
118+ air_density: float = 1.204
119+
120+ # --- Rendering -----------------------------------------------------------
121+ sample_rate: int = 48000
122+ oversample: int = 16 # solver runs at oversample * sample_rate
123+ duration_s: float = 0.45
124+ pre_delay_ms: float = 3.0 # silence before the attack
125+ normalize: bool = True # scale so the peak hits peak_dbfs
126+ peak_dbfs: float = -3.0 # used only when normalize is True
127+ full_scale_pa: float = 20.0 # when normalize is False, this many
128+ # pascals maps to digital full scale
129+ bit_depth: int = 24 # 16 or 24
130+
131+ def __post_init__(self):
132+ assert 0.0 < self.strike_position < 1.0
133+ assert self.bit_depth in (16, 24)
134+
135+ @property
136+ def listener(self) -> np.ndarray:
137+ return np.array([self.listener_x_m, self.listener_y_m, self.listener_z_m])
138+
139+ @property
140+ def slit_point(self) -> np.ndarray:
141+ return np.array([self.slit_x_m, self.slit_y_m, self.slit_z_m])
142+
143+
144+# ----------------------------------------------------------------------------
145+# Modal basis: free-free Euler-Bernoulli bar
146+# ----------------------------------------------------------------------------
147+
148+# Roots of cos(bL)cosh(bL) = 1 for the free-free beam.
149+_BETA_L = np.array([4.73004074, 7.85320462, 10.99560784, 14.13716549])
150+
151+
152+def _beta_l(n: int) -> float:
153+ """n is 1-based. Above the 4th root the asymptotic form is exact to 1e-9."""
154+ if n <= len(_BETA_L):
155+ return float(_BETA_L[n - 1])
156+ return (2 * n + 1) * math.pi / 2.0
157+
158+
159+def _mode_shape(bl: float, s: np.ndarray) -> np.ndarray:
160+ """Free-free mode shape on s in [0,1], unnormalized."""
161+ sigma = (math.cosh(bl) - math.cos(bl)) / (math.sinh(bl) - math.sin(bl))
162+ x = bl * s
163+ return (np.cosh(x) + np.cos(x)) - sigma * (np.sinh(x) + np.sin(x))
164+
165+
166+class ModalBasis:
167+ """Mass-normalized bending modes and their couplings.
168+
169+ phi_c[n] : mode shape at the strike point (1/sqrt(kg))
170+ A[n] : volume displacement into the cavity (m^2/sqrt(kg))
171+ phi_elem[n,j] : mode shape at surface element j (1/sqrt(kg))
172+ """
173+
174+ def __init__(self, cfg: WoodblockConfig):
175+ L, W, H = cfg.length_m, cfg.width_m, cfg.wall_thickness_m
176+ area = W * H
177+ mass = cfg.participating_mass_kg or (cfg.density * area * L)
178+ c_long = math.sqrt(cfg.youngs_modulus / cfg.density)
179+
180+ # f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)), I = W H^3/12
181+ prefactor = c_long * H / (2.0 * math.pi * math.sqrt(12.0) * L * L)
182+
183+ s = np.linspace(0.0, 1.0, 4096)
184+ M = cfg.radiating_elements
185+ s_elem = (np.arange(M) + 0.5) / M # element centres along the bar
186+
187+ freqs, phi_c, A, phi_elem = [], [], [], []
188+ lo, hi = cfg.cavity_wall_span
189+ cav_mask = (s >= lo) & (s <= hi)
190+
191+ for n in range(1, cfg.n_modes + 1):
192+ bl = _beta_l(n)
193+ f = prefactor * bl * bl * cfg.hollow_softening
194+ if n > 1 and cfg.inharmonicity != 1.0:
195+ f0 = prefactor * _beta_l(1) ** 2 * cfg.hollow_softening
196+ f = f0 * (f / f0) ** cfg.inharmonicity
197+ if f > cfg.max_mode_hz:
198+ break
199+
200+ shape = _mode_shape(bl, s)
201+ # Normalize so that (1/L) * integral(phi^2 dx) = 1, then divide by
202+ # sqrt(total mass) to get the mass-normalized mode.
203+ rms = math.sqrt(np.trapezoid(shape ** 2, s))
204+ shape = shape / (rms * math.sqrt(mass))
205+
206+ freqs.append(f)
207+ phi_c.append(float(np.interp(cfg.strike_position, s, shape)))
208+ A.append(W * L * float(np.trapezoid(np.where(cav_mask, shape, 0.0), s)))
209+ phi_elem.append(np.interp(s_elem, s, shape))
210+
211+ self.f = np.array(freqs)
212+ self.omega = 2.0 * math.pi * self.f
213+ self.phi_c = np.array(phi_c)
214+ self.A = np.array(A) * cfg.cavity_coupling
215+ self.phi_elem = np.array(phi_elem) # (N, M)
216+ self.s_elem = s_elem
217+ self.zeta = 0.5 * (cfg.loss_factor + cfg.loss_slope * self.omega)
218+ self.n = len(freqs)
219+ self.mass = mass
220+
221+
222+# ----------------------------------------------------------------------------
223+# Time-domain solver
224+# ----------------------------------------------------------------------------
225+
226+_LAST_FC = [0.0] # scratch slot for contact-force diagnostics
227+
228+
229+def simulate(cfg: WoodblockConfig):
230+ """Integrate the coupled system.
231+
232+ Returns (modal accelerations (N,T), slit-plug acceleration (T,), basis, diag).
233+ """
234+ basis = ModalBasis(cfg)
235+ N = basis.n
236+ if N == 0:
237+ raise ValueError("No modes below max_mode_hz -- check geometry/material.")
238+
239+ fs_hi = cfg.sample_rate * cfg.oversample
240+ dt = 1.0 / fs_hi
241+ n_steps = int(cfg.duration_s * fs_hi)
242+
243+ om2 = basis.omega ** 2
244+ two_zeta_om = 2.0 * basis.zeta * basis.omega
245+ phi_c, A = basis.phi_c, basis.A
246+
247+ rho0, c = cfg.air_density, cfg.speed_of_sound
248+ V, S, Lslit = cfg.cavity_volume_m3, cfg.slit_area_m2, cfg.slit_length_m
249+ # End correction on the slit: ~0.85 * equivalent radius per open end.
250+ Leff = Lslit + 2.0 * 0.85 * math.sqrt(S / math.pi)
251+ omega_H = c * math.sqrt(S / (V * Leff))
252+ two_zeta_H_om = 2.0 * cfg.cavity_damping * omega_H
253+ p_scale = rho0 * c ** 2 / V
254+ cav_drive = c ** 2 / (V * Leff)
255+
256+ K = cfg.mallet_stiffness
257+ alpha = cfg.contact_exponent
258+ mu = cfg.contact_damping
259+ inv_mh = 1.0 / cfg.mallet_mass
260+
261+ # State: [q(N), qdot(N), xi, xidot, y, ydot]
262+ IX, IXD, IY, IYD = 2 * N, 2 * N + 1, 2 * N + 2, 2 * N + 3
263+ state = np.zeros(2 * N + 4)
264+ state[IYD] = cfg.mallet_velocity
265+
266+ def deriv(s):
267+ q, qd = s[:N], s[N:2 * N]
268+ xi, xid, y, yd = s[IX], s[IXD], s[IY], s[IYD]
269+
270+ # Mallet contact
271+ w = float(phi_c @ q)
272+ wd = float(phi_c @ qd)
273+ delta = y - w
274+ if delta > 0.0:
275+ # Hunt-Crossley: f = K d^a (1 + lambda * ddot). lambda has units of
276+ # s/m, so the correction is dimensionless and stays O(1) over a
277+ # sub-millisecond collision. (Writing the loss as mu*d(d^a)/dt
278+ # instead makes mu a relaxation *time*; any plausible-looking value
279+ # then dwarfs the elastic term and fires one huge impulse.)
280+ fc = K * delta ** alpha * (1.0 + mu * (yd - wd))
281+ if fc < 0.0:
282+ fc = 0.0
283+ else:
284+ fc = 0.0
285+ _LAST_FC[0] = fc
286+
287+ # Cavity pressure from combined volume change of walls and slit plug
288+ vol = S * xi + float(A @ q)
289+ pc = p_scale * vol
290+
291+ qdd = -two_zeta_om * qd - om2 * q + fc * phi_c - A * pc
292+ xidd = -two_zeta_H_om * xid - omega_H ** 2 * xi - cav_drive * float(A @ q)
293+ ydd = -fc * inv_mh
294+
295+ d = np.empty_like(s)
296+ d[:N] = qd
297+ d[N:2 * N] = qdd
298+ d[IX] = xid
299+ d[IXD] = xidd
300+ d[IY] = yd
301+ d[IYD] = ydd
302+ return d, qdd, xidd
303+
304+ Qdd = np.empty((N, n_steps))
305+ Xdd = np.empty(n_steps)
306+ force = np.zeros(n_steps)
307+
308+ for k in range(n_steps):
309+ k1, qdd, xidd = deriv(state)
310+ Qdd[:, k] = qdd
311+ Xdd[k] = xidd
312+ force[k] = _LAST_FC[0]
313+
314+ k2, _, _ = deriv(state + 0.5 * dt * k1)
315+ k3, _, _ = deriv(state + 0.5 * dt * k2)
316+ k4, _, _ = deriv(state + dt * k3)
317+ state += (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
318+
319+ contact = force > 0.0
320+ diag = {
321+ "contact_force": force,
322+ "peak_force_N": float(force.max()),
323+ "contact_ms": float(contact.sum() * dt * 1e3),
324+ "helmholtz_hz": omega_H / (2.0 * math.pi),
325+ }
326+ return Qdd, Xdd, basis, diag
327+
328+
329+# ----------------------------------------------------------------------------
330+# Radiation: free-field pressure at a point
331+# ----------------------------------------------------------------------------
332+
333+def radiate(cfg: WoodblockConfig, Qdd: np.ndarray, Xdd: np.ndarray,
334+ basis: ModalBasis, fs_hi: float):
335+ """Sum retarded monopole contributions from every source to the listener.
336+
337+ p(r,t) = sum_j rho0 / (4 pi d_j) * Qddot_j(t - d_j/c)
338+
339+ The struck face and the opposite face move together (each cross-section
340+ translates), so they inject equal and opposite volume -- each slice of the
341+ bar is a local dipole. Summed over the bar the net volume displacement is
342+ exactly zero, which is the correct statement that a bending bar has no
343+ monopole component. The audible near-field asymmetry comes entirely from
344+ the different path lengths and delays, not from any fudge factor.
345+ """
346+ rho0, c = cfg.air_density, cfg.speed_of_sound
347+ T = Qdd.shape[1]
348+ L, W = cfg.length_m, cfg.width_m
349+ half_h = 0.5 * cfg.block_height_m
350+ M = cfg.radiating_elements
351+ patch = W * L / M # area represented by one element
352+
353+ # Surface element positions along the bar, on both faces
354+ x_elem = (basis.s_elem - 0.5) * L
355+ src_pos, src_sign, src_mode = [], [], []
356+ for j in range(M):
357+ src_pos.append([x_elem[j], 0.0, +half_h]); src_sign.append(+1.0); src_mode.append(j)
358+ src_pos.append([x_elem[j], 0.0, -half_h]); src_sign.append(-1.0); src_mode.append(j)
359+ src_pos = np.array(src_pos)
360+
361+ r = cfg.listener
362+ dists = np.linalg.norm(src_pos - r, axis=1)
363+ slit_dist = float(np.linalg.norm(cfg.slit_point - r))
364+
365+ max_delay = int(math.ceil(max(dists.max(), slit_dist) / c * fs_hi)) + 2
366+ out = np.zeros(T + max_delay)
367+
368+ # Surface (body) sources -- volume acceleration = patch area * wddot
369+ for idx in range(len(dists)):
370+ j = src_mode[idx]
371+ wdd = basis.phi_elem[:, j] @ Qdd # (T,) surface acceleration
372+ amp = src_sign[idx] * patch * rho0 / (4.0 * math.pi * dists[idx])
373+ k = int(round(dists[idx] / c * fs_hi))
374+ out[k:k + T] += amp * wdd
375+
376+ # Slit monopole -- volume acceleration = slit area * plug acceleration
377+ amp = cfg.slit_area_m2 * rho0 / (4.0 * math.pi * slit_dist)
378+ k = int(round(slit_dist / c * fs_hi))
379+ out[k:k + T] += amp * Xdd
380+
381+ return out[:T]
382+
383+
384+def check_setup(cfg: WoodblockConfig, basis: ModalBasis) -> list:
385+ """Warn about configurations where the radiation model stops being valid."""
386+ msgs = []
387+ c = cfg.speed_of_sound
388+ spacing = cfg.length_m / cfg.radiating_elements
389+ f_max = basis.f[-1]
390+ if spacing > 0.25 * c / f_max:
391+ need = int(math.ceil(cfg.length_m / (0.25 * c / f_max)))
392+ msgs.append(f"element spacing is coarse at {f_max:.0f} Hz; "
393+ f"consider --radiating_elements {need}")
394+ d = float(np.linalg.norm(cfg.listener))
395+ if d < 0.5 * cfg.length_m:
396+ msgs.append(f"listener is {d*100:.0f} cm from the block centre, inside "
397+ f"the source region; point-source elements are crude here")
398+ return msgs
399+
400+
401+# ----------------------------------------------------------------------------
402+# Rendering
403+# ----------------------------------------------------------------------------
404+
405+def render(cfg: WoodblockConfig):
406+ Qdd, Xdd, basis, diag = simulate(cfg)
407+ fs_hi = cfg.sample_rate * cfg.oversample
408+
409+ p = radiate(cfg, Qdd, Xdd, basis, fs_hi) # pascals at the listener
410+ diag["peak_pa"] = float(np.abs(p).max())
411+ diag["peak_spl_db"] = 20.0 * math.log10(diag["peak_pa"] / 20e-6 + 1e-30)
412+ diag["warnings"] = check_setup(cfg, basis)
413+
414+ # Anti-aliased decimation to the target rate
415+ y = resample_poly(p, 1, cfg.oversample)
416+
417+ # Pre-roll so the attack transient does not sit against sample 0. (Do NOT
418+ # try to de-DC by subtracting the head of the buffer -- the strike is at
419+ # t=0, so the head is the loudest part of the file.)
420+ pad = int(cfg.pre_delay_ms * 1e-3 * cfg.sample_rate)
421+ if pad > 0:
422+ y = np.concatenate([np.zeros(pad), y])[:int(cfg.duration_s * cfg.sample_rate)]
423+
424+ fade = min(256, y.size // 8)
425+ if fade > 1:
426+ y[-fade:] *= np.linspace(1.0, 0.0, fade) ** 2
427+
428+ if cfg.normalize:
429+ peak = np.max(np.abs(y))
430+ if peak > 0:
431+ y = y / peak * (10.0 ** (cfg.peak_dbfs / 20.0))
432+ else:
433+ y = y / cfg.full_scale_pa
434+ if np.max(np.abs(y)) > 1.0:
435+ diag["warnings"].append(
436+ f"clipping: peak is {diag['peak_pa']:.2f} Pa but full_scale_pa "
437+ f"is {cfg.full_scale_pa:.2f}")
438+ return y, basis, diag
439+
440+
441+def write_wav(path: str, y: np.ndarray, cfg: WoodblockConfig) -> None:
442+ y = np.clip(y, -1.0, 1.0)
443+ with wave.open(path, "wb") as w:
444+ w.setnchannels(1)
445+ w.setsampwidth(cfg.bit_depth // 8)
446+ w.setframerate(cfg.sample_rate)
447+ if cfg.bit_depth == 16:
448+ w.writeframes((y * 32767.0).astype("<i2").tobytes())
449+ else:
450+ i32 = (y * 8388607.0).astype("<i4")
451+ raw = np.frombuffer(i32.tobytes(), dtype=np.uint8).reshape(-1, 4)
452+ w.writeframes(raw[:, :3].tobytes())
453+
454+
455+# ----------------------------------------------------------------------------
456+# CLI
457+# ----------------------------------------------------------------------------
458+
459+def build_parser() -> argparse.ArgumentParser:
460+ p = argparse.ArgumentParser(
461+ description="Simulate a woodblock strike and sample the pressure "
462+ "at a point in space.")
463+ p.add_argument("-o", "--output", default="woodblock.wav")
464+ defaults = WoodblockConfig()
465+ for f in dataclasses.fields(defaults):
466+ val = getattr(defaults, f.name)
467+ if isinstance(val, tuple):
468+ continue
469+ cast = (lambda s: bool(int(s))) if isinstance(val, bool) else type(val)
470+ p.add_argument(f"--{f.name}", type=cast, default=None,
471+ help=f"(default: {val})")
472+ return p
473+
474+
475+def main() -> None:
476+ args = build_parser().parse_args()
477+ overrides = {k: v for k, v in vars(args).items()
478+ if k != "output" and v is not None}
479+ cfg = WoodblockConfig(**overrides)
480+
481+ y, basis, diag = render(cfg)
482+ write_wav(args.output, y, cfg)
483+
484+ d = float(np.linalg.norm(cfg.listener))
485+ print(f"modes kept : {basis.n}")
486+ print("partials Hz: " +
487+ ", ".join(f"{f:.0f}" for f in basis.f[:6]) +
488+ ("..." if basis.n > 6 else ""))
489+ print("ratios : " +
490+ ", ".join(f"{f / basis.f[0]:.2f}" for f in basis.f[:6]))
491+ print(f"helmholtz : {diag['helmholtz_hz']:.0f} Hz")
492+ print(f"contact : {diag['contact_ms']:.3f} ms, peak {diag['peak_force_N']:.0f} N")
493+ print(f"listener : {d * 100:.0f} cm from block centre")
494+ print(f"peak press.: {diag['peak_pa']:.3f} Pa ({diag['peak_spl_db']:.1f} dB SPL)")
495+ print(f"length : {y.size / cfg.sample_rate:.3f} s @ {cfg.sample_rate} Hz")
496+ for w in diag["warnings"]:
497+ print(f"WARNING : {w}")
498+ print(f"wrote : {args.output}")
499+
500+
501+if __name__ == "__main__":
502+ main()