Plate: new app, square simply supported linear plate on a grid (exact sine modes, mallet or kick at a point, height map, sped-up note with spectrogram)
11 changed files+946−1
README.mdmodified+2−1View file
@@ -5,6 +5,7 @@ Interactive apps for "Listening for chaos": struck dynamical systems heard as so
55 - `src/double-pendulum/`: double pendulum ([README](src/double-pendulum/README.md)), served at `/double-pendulum/`.
66 - `src/fput-chain/`: FPUT-beta chain ([README](src/fput-chain/README.md)), served at `/fput-chain/`.
77 - `src/beam/`: beam with bending stiffness ([README](src/beam/README.md)), served at `/beam/`.
8+- `src/plate/`: square plate with bending stiffness, simply supported ([README](src/plate/README.md)), served at `/plate/`.
89 - `src/common/`: code shared by the apps (help dialogs, chart cursor, WAV export, note player, ghost controls, note render loop, mode bars and mode-energy heat maps, index colors, mute button, number formats, canvas sizing, λ vs E scatter, strip charts, whole-note spectrogram, Lyapunov exponent by the two-trajectory method with a per-app `Shadow`; Poincare section orbits, their worker, the precomputed atlas and the section panel, with a per-app `Section` and panel configuration) and `style.css`, the shared styles. Each app's `main.ts` imports `../common/style.css` and then its own `style.css`, so app rules come after the shared ones in dev and in the build.
910 - `src/index.html`: landing page linking the apps.
1011
@@ -17,4 +18,4 @@ npm install
1718 npm run dev
1819 ```
1920
20-The dev server serves every app (http://localhost:5173/double-pendulum/, http://localhost:5173/fput-chain/, http://localhost:5173/beam/); edits anywhere under `src/` reload the open pages.
21+The dev server serves every app (http://localhost:5173/double-pendulum/, http://localhost:5173/fput-chain/, http://localhost:5173/beam/, http://localhost:5173/plate/); edits anywhere under `src/` reload the open pages.
src/index.htmlmodified+1−0View file
@@ -34,6 +34,7 @@
3434 <li><a href="./double-pendulum/"><strong>Double pendulum</strong><div class="muted">Strike it, watch and hear the motion turn regular as the energy decays; Lyapunov exponent and Poincaré section.</div></a></li>
3535 <li><a href="./fput-chain/"><strong>FPUT chain</strong><div class="muted">A struck chain of N masses with cubic springs; energy spreading up the modes, Lyapunov exponent and Poincaré sections for N = 2 and 3.</div></a></li>
3636 <li><a href="./beam/"><strong>Beam</strong><div class="muted">A struck beam of N masses with bending stiffness and stretching: clamped, pinned or cantilever; inharmonic, cymbal-like overtones.</div></a></li>
37+ <li><a href="./plate/"><strong>Plate</strong><div class="muted">A struck square plate of N × N masses with bending stiffness, simply supported; linear for now, modes (m, n) at m² + n².</div></a></li>
3738 </ul>
3839 </main>
3940 </body>
src/plate/README.mdadded+9−0View file
@@ -0,0 +1,9 @@
1+# plate
2+
3+Struck square plate, step 4 of "Listening for chaos" (after double-pendulum, fput-chain and beam), the first two-dimensional system on the way to cymbal and gong. Work in progress: step 1 (2026-10-09), a linear, simply supported square plate with mallet strike, animation and sped-up note with spectrogram.
4+
5+Current state: N × N masses (N per side, 2 to 32, default 8) on vertical rails at the nodes of a square grid, spacing a, moving up and down only. Bending κ on the five-point Laplacian at each node, energy κ (L u)² / 2, force −κ L(L u). The edge nodes are supports at u = 0 with no bending term (mirror ghost), which makes the edge simply supported. Linear: no stretching yet. Damping as in the beam: drag (τ, all modes alike) and a frequency-dependent loss (dashpots alongside the grid springs, force −c L v, set by τ top at the top of the band). Parameters are given for the 8 × 8 plate (m 1 kg, κ 20 N/m, a 0.5 m, so the side is 4.5 m; τ 540 s, τ top 200 s); Masses per side refines the same plate (`refine`: m with the cell area, κ inversely, the loss coefficient c unchanged). Strike: a kick (pyramid velocity bump of half-width 0.5 m) or the shared Hertz mallet (`../common/mallet.ts`, M 1 kg, K 1000) over the same bump, at (x, y) m from the lower left corner (default (1.3, 1.1)); clicking the plate moves the strike point. The card shows the plate from above as a height map (red up, blue down, scaled to the recent largest |u|), with the strike point (cross) and the listening point (circle). The sound panel is the beam's, reduced: velocity, displacement or acceleration at a listening point (default (3.1, 1.7), bilinear between masses), pitch of the lowest mode, the whole decay as a note through the shared renderer and anti-alias filter, spectrogram with dashed lines at the distinct mode pitches, energy chart, Replay and Save WAV. Not yet: Lyapunov exponent, ghost, mode bars and histories, cross-section.
6+
7+Modes: L is diagonal in the discrete sine basis, so the modes are exact, S_mn ∝ sin(mπ(i + 1)/(N + 1)) sin(nπ(j + 1)/(N + 1)), with ω_mn = √(κ/m) (l_m + l_n), l_m = 4 sin²(mπ/(2(N + 1))) (`Modes`, sorted by frequency with (m, n)). For low modes ω_mn ∝ m² + n², the continuum plate; (m, n) and (n, m) are degenerate. Ratios to the lowest, against the continuum's 2.5 and 4.0: (1, 2) 2.309 and (2, 2) 3.618 at N = 4, 2.440 and 3.879 at N = 8, 2.483 and 3.966 at N = 16, 2.496 and 3.991 at N = 32; ω_11 stays at 1.054 to 1.089 rad/s across N. Check (2026-10-09, N = 8, kick at (1.3, 1.1), undamped, 4000 s, displacement at (3.1, 1.7)): 32 of the 33 distinct frequencies appear as spectral peaks within 0.5 cents of the Verlet-corrected prediction (2/h) asin(ωh/2) (FFT resolution 2.5 cents); the missing one, (7, 6), is barely excited there (peak 1e-4 of the largest). Energy wobble 1.2e-3 (bounded). Measured decay times match `decayTimes` (exact, L is diagonal in the modes) to the energy wobble: (1, 1) 498 s, (3, 5) 249 s, (8, 8) 146 s.
8+
9+Integrator: velocity Verlet at ω h = 0.2 at the top of the band (as the beam), linear so the same step at every energy. Drag and loss in half steps around it (Strang). The loss is an explicit step p → (I − x) p, x = (c dt / m) L, whose largest eigenvalue is 2 dt / τ top (about 3e-5 by default): it never adds energy, and each mode's decay rate is off by about x / 2 relative (substeps if 2 dt / τ top passes 0.01). Unlike the beam's Crank–Nicolson step (a tridiagonal solve), a Crank–Nicolson step on the grid needs an iterative solve; with Jacobi it cost 4 to 5 times the rest of the step, for no audible difference. Cost per step: 1.2 µs at N = 8, 4.7 µs at N = 16 (Node, this laptop); the default note (4.9 s of sound, about 3400 s of plate time) renders in 1.2 s at N = 8 and 10 s at N = 16. The step shrinks as 1 / N² (ω_max / ω_11 is 32 at N = 8, 117 at N = 16, 441 at N = 32), so the cost of a note grows as N⁴; the modal integrator planned later removes the step limit.
src/plate/index.htmladded+115−0View file
@@ -0,0 +1,115 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
6+ <title>Plate chaos</title>
7+ </head>
8+ <body>
9+ <main>
10+ <header><h1>Plate chaos</h1></header>
11+ <div class="grid">
12+ <section class="card">
13+ <canvas id="plate"></canvas>
14+ <div id="controls">
15+ <div class="row">
16+ <button id="strike" class="primary">Strike</button>
17+ <select id="strike-kind">
18+ <option value="at">kick at</option>
19+ <option value="mallet" selected>mallet at</option>
20+ </select>
21+ <span class="unit">x</span> <input id="strike-x" type="number" min="0" value="1.3" step="0.05" />
22+ <span class="unit">y</span> <input id="strike-y" type="number" min="0" value="1.1" step="0.05" /> <span class="unit">m</span>
23+ <button id="pause" disabled>Pause</button>
24+ <button type="button" class="help" data-help="plate" title="About the plate" aria-label="About the plate">?</button>
25+ </div>
26+ <div class="row">
27+ <span class="muted">Strike energy</span>
28+ <input id="energy" type="range" min="-5" max="2" step="0.01" value="0.4" />
29+ <span id="energy-value" class="num"></span>
30+ </div>
31+ <div class="row"><span id="strike-info" class="muted"></span></div>
32+ <div class="row">
33+ <span class="muted">Masses per side</span>
34+ <input id="n" type="number" min="2" max="32" step="1" value="8" />
35+ <span class="muted" style="margin-left: 12px">Speed</span>
36+ <select id="speed">
37+ <option value="0.25">0.25×</option>
38+ <option value="0.5">0.5×</option>
39+ <option value="1" selected>1× (real time)</option>
40+ <option value="2">2×</option>
41+ <option value="5">5×</option>
42+ <option value="10">10×</option>
43+ <option value="20">20×</option>
44+ </select>
45+ </div>
46+ </div>
47+ </section>
48+ <section class="card">
49+ <div class="card-head">Parameters <button type="button" class="help" data-help="params" title="About the parameters" aria-label="About the parameters">?</button></div>
50+ <div id="param-form">
51+ <label>m <input id="p-m" type="number" step="any" /> <span class="unit">kg</span></label>
52+ <label>κ <input id="p-kappa" type="number" step="any" /> <span class="unit">N/m (bending)</span></label>
53+ <label>a <input id="p-a" type="number" step="any" /> <span class="unit">m (spacing)</span></label>
54+ <label>τ <input id="p-tau" type="number" step="any" /> <span class="unit">s (decay, 0 = off)</span></label>
55+ <label>τ top <input id="p-tautop" type="number" step="any" /> <span class="unit">s (decay of the highest modes from the loss, 0 = off)</span></label>
56+ <label>mallet M <input id="p-mallet-m" type="number" step="any" /> <span class="unit">kg</span></label>
57+ <label>mallet K <input id="p-mallet-k" type="number" step="any" /> <span class="unit">N/m^1.5 (contact stiffness)</span></label>
58+ </div>
59+ <div id="param-actions" class="row">
60+ <button id="param-apply">Apply</button>
61+ <span id="param-status"></span>
62+ </div>
63+ <div id="params"></div>
64+ </section>
65+ </div>
66+ <details class="card panel" id="sound" open>
67+ <summary><span class="label">Sped up · listen</span> <span id="sound-clock" class="muted"></span> <button type="button" class="help" data-help="sound" title="About the sped-up run" aria-label="About the sped-up run">?</button></summary>
68+ <div class="row">
69+ <button id="mute" type="button" class="mute" title="Unmute" aria-label="Unmute"></button>
70+ <label><span class="label">Signal</span> <select id="signal"></select></label>
71+ <span><span class="unit">x</span> <input id="sound-x" type="number" min="0" step="0.05" value="3.1" aria-label="Listening point x (m)" />
72+ <span class="unit">y</span> <input id="sound-y" type="number" min="0" step="0.05" value="1.7" aria-label="Listening point y (m)" /> <span class="unit">m</span></span>
73+ <label><span class="label">Volume</span> <input id="volume" type="range" min="0" max="1" step="0.01" value="0.5" /></label>
74+ <button id="replay" type="button" disabled title="Play the last strike again with the current settings">Replay</button>
75+ <button id="save-wav" type="button" disabled title="Save the last note as a WAV file">Save WAV</button>
76+ <span id="sound-chips" class="chips">
77+ <label class="chip"><input id="spec-toggle" type="checkbox" checked /><span>Spectrogram</span></label>
78+ <label class="chip" title="Logarithmic frequency axis (off: linear from 0)"><input id="spec-log" type="checkbox" checked /><span>log f</span></label>
79+ </span>
80+ </div>
81+ <div class="row">
82+ <label class="grow"><span class="label">Pitch</span> <input id="pitch" type="range" min="0.699" max="3" step="0.001" value="2" aria-label="Pitch of the lowest mode" /></label>
83+ <span id="pitch-value" class="num"></span>
84+ </div>
85+ <div class="row"><span id="alias" class="muted" hidden></span></div>
86+ <div class="row"><span id="sound-status" class="muted"></span></div>
87+ <p id="sound-hint" class="muted">Strike to hear the whole decay.</p>
88+ <figure id="spec-fig"><div id="spectrogram"></div></figure>
89+ <div id="sound-traces"></div>
90+ </details>
91+ </main>
92+ <dialog id="help">
93+ <div id="help-body"></div>
94+ <div id="help-actions"><button id="help-close" type="button">Close</button></div>
95+ </dialog>
96+
97+ <template id="help-plate">
98+ <h2>The plate</h2>
99+ <p>A square plate as a grid of equal masses on vertical rails, moving up and down only, seen from above: red is up, blue down, white at rest (scaled to fit). Each mass resists bending (κ, on how far it sits off the mean of its four neighbours). The plate is simply supported: its edge is held at zero height but free to turn. It is linear for now, so a weak and a strong strike sound the same; the nonlinearity (von Kármán stretching) comes next.</p>
100+ <p>Strike gives the energy on the slider as a hammer blow (a short pyramid bump of velocity, one spacing of the 8 × 8 plate wide on each side) or with a mallet (mass M, contact stiffness K in Parameters) that presses into the plate over the same bump through a Hertz contact (force K δ^1.5) and lifts off. The point is in metres from the lower left corner; click the plate to move it there. The cross marks the strike point, the circle the listening point. The space bar pauses.</p>
101+ </template>
102+ <template id="help-params">
103+ <h2>Parameters</h2>
104+ <p>Mass, bending κ, the spacing a, the decay time τ of the amplitude (drag on the masses, all modes alike; 0 for none), and τ top, the decay time a frequency-dependent loss (dashpots alongside the grid springs) gives the highest modes (at 8 × 8 masses), less for lower ones (0 for none). Apply restarts the plate at rest.</p>
105+ <p>The values are for the plate with 8 masses per side. Masses per side refines the same plate: same size, strike energies and decay, with each mass scaled with the area of its cell and κ inversely, so the low modes stay put and more masses add higher ones.</p>
106+ <p>The normal modes are products of sines, numbered (m, n) by their half waves along x and y. Their frequencies go as m² + n² for the low ones (the continuum plate), so (1, 2) and (2, 1) share a pitch, and the ratios to the lowest, 2.5, 4, 5, 6.5, 8.5, ..., are not harmonic. The derived values list the lowest modes and how close the grid is to the continuum.</p>
107+ </template>
108+ <template id="help-sound">
109+ <h2>Sped up · listen</h2>
110+ <p>A copy of the same strike, run fast enough to be heard: Pitch sets where the lowest normal mode lands, and the whole decay plays in a few seconds. While this panel is open each Strike plays a note; Replay plays the last one again with the current settings.</p>
111+ <p>Signal chooses what is heard, at a listening point (x, y in metres, between masses by interpolation). Save WAV saves the last note as heard (16-bit mono, normalized). The spectrogram has dashed lines at the normal-mode pitches. Notes are anti-alias filtered (motion above 20 kHz at 48 kHz is removed rather than folded back); a note says when the top modes are above that and so are not heard.</p>
112+ </template>
113+ <script type="module" src="./main.ts"></script>
114+ </body>
115+</html>
src/plate/main.tsadded+476−0View file
@@ -0,0 +1,476 @@
1+import "../common/style.css";
2+import "./style.css";
3+import {
4+ Modes, REF_N, Verlet, copyState, decayTimes, energy, forces, hammerShape, isDamped, newState, plateSide, refine, stepFor, strikeAt,
5+ type Params, type State,
6+} from "./plate";
7+import { fit } from "../common/canvas";
8+import { newCursor } from "../common/cursor";
9+import { fmtDuration, fmtEnergy, fmtLength } from "../common/format";
10+import { setupHelp } from "../common/help";
11+import { setupMute } from "../common/mute";
12+import { PASS } from "../common/antialias";
13+import { malletStrike, type MalletParams } from "../common/mallet";
14+import { DIAG_EVERY, END_RATIO, UNDAMPED_LENGTH } from "../common/note-render";
15+import { Player } from "../common/sound";
16+import { NoteSpectrogram } from "../common/spectrogram";
17+import { TracePanel } from "../common/traces";
18+import { saveFile, wavBlob } from "../common/wav";
19+import type { RenderRequest, RenderResult } from "./render-worker";
20+import { SIGNALS, signalByKey } from "./signals";
21+
22+setupHelp();
23+
24+const plateCanvas = document.getElementById("plate") as HTMLCanvasElement;
25+const energyInput = document.getElementById("energy") as HTMLInputElement;
26+const energyValue = document.getElementById("energy-value")!;
27+const nInput = document.getElementById("n") as HTMLInputElement;
28+const kindSelect = document.getElementById("strike-kind") as HTMLSelectElement;
29+const strikeXInput = document.getElementById("strike-x") as HTMLInputElement;
30+const strikeYInput = document.getElementById("strike-y") as HTMLInputElement;
31+const strikeInfo = document.getElementById("strike-info")!;
32+const speedSelect = document.getElementById("speed") as HTMLSelectElement;
33+const pauseButton = document.getElementById("pause") as HTMLButtonElement;
34+
35+// Parameters, as the beam's defaults: bending kappa = 20 on unit masses,
36+// drag tau 540 s, the frequency-dependent loss tauTop 200 s. base: the page's
37+// parameters, for the plate at REF_N masses per side; P: the same plate
38+// refined to the chosen N (see refine).
39+let base: Params = { N: REF_N, m: 1, kappa: 20, a: 0.5, tau: 540, tauTop: 200 };
40+let P: Params = base;
41+// The mallet: a physical object, not refined with N.
42+let mallet: MalletParams = { M: 1, K: 1000 };
43+
44+// The plate. It sits at rest until the first Strike; changing parameters
45+// puts it back at rest.
46+let modes: Modes;
47+let s: State;
48+let verlet: Verlet;
49+let H: number; // integrator step (s)
50+let E0 = 0; // energy of the last strike
51+let t = 0;
52+let running = false;
53+let ampl = 0; // recent largest |u|, for scaling the colours
54+let lastStrike: State | null = null; // state right after the last strike (for the sound)
55+
56+const strikeEnergy = () => 10 ** Number(energyInput.value);
57+const showEnergy = () => (energyValue.textContent = fmtEnergy(strikeEnergy()));
58+showEnergy();
59+energyInput.addEventListener("input", showEnergy);
60+
61+// Pause stops everything: the animation and the sound.
62+function setRunning(r: boolean) {
63+ running = r;
64+ pauseButton.textContent = r ? "Pause" : "Resume";
65+ player.setPaused(!r);
66+}
67+pauseButton.addEventListener("click", () => setRunning(!running));
68+document.addEventListener("keydown", (e) => {
69+ const el = e.target as HTMLElement;
70+ if (e.key !== " " || pauseButton.disabled || el.closest("input, select, button, summary")) return;
71+ e.preventDefault();
72+ setRunning(!running);
73+});
74+
75+// A coordinate on the plate (m from the lower left corner), kept on it.
76+const coord = (el: HTMLInputElement) => Math.min(plateSide(P), Math.max(0, Number(el.value) || 0));
77+
78+// label names the strike in saved files.
79+let strikeLabel = "";
80+document.getElementById("strike")!.addEventListener("click", () => {
81+ const E = strikeEnergy();
82+ const x = coord(strikeXInput), y = coord(strikeYInput);
83+ strikeXInput.value = String(x);
84+ strikeYInput.value = String(y);
85+ const where = `x${x}y${y}m`;
86+ strikeInfo.textContent = "";
87+ // The hammer's half-width is one spacing of the plate at REF_N, so the
88+ // blow is the same at any N.
89+ if (kindSelect.value === "at") {
90+ strikeLabel = where;
91+ s = strikeAt(P, x, y, base.a, E);
92+ } else {
93+ strikeLabel = `mallet${mallet.M}kg${mallet.K}-${where}`;
94+ const phi = hammerShape(P, x, y, base.a), sum = phi.reduce((a, b) => a + b, 0);
95+ const r = malletStrike(
96+ { N: P.N * P.N, m: P.m, forces: (u, F) => forces(P, u, F), weights: phi.map((v) => v / sum), step: H },
97+ mallet, E,
98+ );
99+ s = { u: r.u, p: r.p, W: 0 };
100+ const share = energy(P, s) / E;
101+ strikeInfo.textContent =
102+ `mallet: contact ${(1000 * r.contactTime).toPrecision(3)} ms, ${(100 * share).toFixed(0)}% of the energy into the plate` +
103+ `${share < 0.995 ? `, ${(100 * (1 - share)).toFixed(0)}% leaves with the mallet` : ""}` +
104+ `; peak force ${r.peakForce.toPrecision(3)} N; shown from lift-off`;
105+ }
106+ lastStrike = copyState(s);
107+ verlet = new Verlet(P, s);
108+ E0 = energy(P, s);
109+ t = 0;
110+ ampl = 0;
111+ replayButton.disabled = false;
112+ renderNote();
113+ pauseButton.disabled = false;
114+ setRunning(true);
115+});
116+
117+// Clicking the plate moves the strike point there.
118+plateCanvas.addEventListener("click", (e) => {
119+ const g = geometry();
120+ const rect = plateCanvas.getBoundingClientRect();
121+ const L = plateSide(P);
122+ const fx = (e.clientX - rect.left - g.x0) / g.size, fy = 1 - (e.clientY - rect.top - g.y0) / g.size;
123+ if (fx < 0 || fx > 1 || fy < 0 || fy > 1) return;
124+ strikeXInput.value = (L * fx).toFixed(2);
125+ strikeYInput.value = (L * fy).toFixed(2);
126+});
127+
128+// Energy balance: energy plus what damping has dissipated, against the strike's.
129+const relErr = () => (E0 > 0 ? (energy(P, s) + s.W - E0) / E0 : 0);
130+
131+// Parameters: numbers in the form; Apply restarts at rest. N is set above.
132+const paramInputs = {
133+ m: document.getElementById("p-m") as HTMLInputElement,
134+ kappa: document.getElementById("p-kappa") as HTMLInputElement,
135+ a: document.getElementById("p-a") as HTMLInputElement,
136+ tau: document.getElementById("p-tau") as HTMLInputElement,
137+ tauTop: document.getElementById("p-tautop") as HTMLInputElement,
138+};
139+const paramStatus = document.getElementById("param-status")!;
140+const malletInputs = {
141+ M: document.getElementById("p-mallet-m") as HTMLInputElement,
142+ K: document.getElementById("p-mallet-k") as HTMLInputElement,
143+};
144+function showParams() {
145+ for (const [key, el] of Object.entries(paramInputs)) el.value = String(base[key as keyof typeof paramInputs]);
146+ malletInputs.M.value = String(mallet.M);
147+ malletInputs.K.value = String(mallet.K);
148+ showDerived();
149+}
150+// Values that follow from the parameters: the lowest modes with their (m, n)
151+// (degenerate pairs once), their ratios to the lowest against the
152+// continuum's (m^2 + n^2) / 2, and the decay times.
153+function showDerived() {
154+ const w1 = modes.omega[0], hz = (w: number) => w / (2 * Math.PI);
155+ const rows: string[] = [];
156+ for (let k = 0; k < modes.count && rows.length < 8; k++) {
157+ const m = modes.m[k], n = modes.n[k];
158+ if (m > n) continue; // (n, m) has the same frequency
159+ rows.push(
160+ ` (${m}, ${n})${m === n ? " " : " ×2 "}${hz(modes.omega[k]).toPrecision(3).padStart(7)} Hz ` +
161+ `ratio ${(modes.omega[k] / w1).toFixed(3)} continuum ${((m * m + n * n) / 2).toFixed(1)}`,
162+ );
163+ }
164+ const g = (x: number) => (x === 0 ? "0" : x.toPrecision(4));
165+ let decay = "none (no damping)";
166+ if (isDamped(P)) {
167+ const d = decayTimes(P, modes);
168+ decay = `${fmtDuration(d[0])} (mode 1, 1) … ${fmtDuration(d[modes.count - 1])} (mode ${P.N}, ${P.N}), amplitude`;
169+ }
170+ document.getElementById("params")!.textContent =
171+ (P.N === REF_N ? "" : `at ${P.N} × ${P.N}: m ${g(P.m)} kg, κ ${g(P.kappa)} N/m, a ${g(P.a)} m\n`) +
172+ `side ${fmtLength(plateSide(P))}, ${modes.count} modes, ${hz(w1).toPrecision(3)} … ${hz(modes.omegaMax).toPrecision(3)} Hz ` +
173+ `(top ${(modes.omegaMax / w1).toFixed(1)}× the lowest)\n` +
174+ `lowest modes (m, n):\n${rows.join("\n")}\n` +
175+ `decay time ${decay}\n` +
176+ `integrator step ${(1000 * H).toPrecision(3)} ms`;
177+}
178+
179+// Sound: the last strike again, sped up so the lowest mode sounds at the
180+// chosen pitch, rendered in a worker, as in the beam.
181+const soundPanel = document.getElementById("sound") as HTMLDetailsElement;
182+const replayButton = document.getElementById("replay") as HTMLButtonElement;
183+const saveButton = document.getElementById("save-wav") as HTMLButtonElement;
184+const signalSelect = document.getElementById("signal") as HTMLSelectElement;
185+const soundX = document.getElementById("sound-x") as HTMLInputElement;
186+const soundY = document.getElementById("sound-y") as HTMLInputElement;
187+const pitchInput = document.getElementById("pitch") as HTMLInputElement;
188+const soundStatus = document.getElementById("sound-status")!;
189+const soundHint = document.getElementById("sound-hint")!;
190+const soundClock = document.getElementById("sound-clock")!;
191+const aliasWarning = document.getElementById("alias")!;
192+for (const sig of SIGNALS) signalSelect.add(new Option(sig.label, sig.key));
193+signalSelect.value = "v";
194+
195+const player = new Player();
196+setupMute(document.getElementById("mute") as HTMLButtonElement, player);
197+document.getElementById("volume")!.addEventListener("input", (e) =>
198+ player.setVolume(Number((e.target as HTMLInputElement).value)),
199+);
200+const pitch = () => 10 ** Number(pitchInput.value); // Hz of the lowest mode
201+const speedup = () => (pitch() * 2 * Math.PI) / modes.omega[0];
202+// Distinct mode pitches (degenerate pairs once), for the spectrogram.
203+const modePitches = () => {
204+ const f: number[] = [];
205+ for (let k = 0; k < modes.count; k++) if (modes.m[k] <= modes.n[k]) f.push((modes.omega[k] / (2 * Math.PI)) * speedup());
206+ return f;
207+};
208+// The slowest mode's energy decays at 2 / tau1, so a drop by END_RATIO takes
209+// about (tau1 / 2) ln(1 / END_RATIO) of physical time.
210+const predictedNoteLength = () =>
211+ isDamped(P) ? ((Math.max(...decayTimes(P, modes)) / 2) * Math.log(1 / END_RATIO)) / speedup() : UNDAMPED_LENGTH;
212+function showPitch() {
213+ const S = speedup(), top = (modes.omegaMax / (2 * Math.PI)) * S;
214+ document.getElementById("pitch-value")!.textContent = `lowest mode ${pitch().toFixed(0)} Hz · highest ${top.toFixed(0)} Hz`;
215+ const edge = PASS * player.sampleRate;
216+ aliasWarning.hidden = top < edge;
217+ aliasWarning.textContent = `The highest modes are above ${(edge / 1000).toFixed(1)} kHz and are filtered out (inaudible, and they would alias); lower the pitch to hear them.`;
218+ const note = predictedNoteLength();
219+ soundClock.textContent = isDamped(P)
220+ ? `the whole decay, ${S.toFixed(0)}× faster: ${fmtDuration(note * S)} of motion in ${fmtDuration(note)}`
221+ : `${S.toFixed(0)}× faster than the plate above`;
222+}
223+pitchInput.addEventListener("input", showPitch);
224+soundPanel.addEventListener("toggle", () => {
225+ if (!soundPanel.open) player.stop();
226+});
227+
228+// Displays of the note: spectrogram and energy, sharing one time cursor.
229+const soundCursor = newCursor();
230+const specFig = document.getElementById("spec-fig")!;
231+const specToggle = document.getElementById("spec-toggle") as HTMLInputElement;
232+specToggle.addEventListener("change", () => (specFig.hidden = !specToggle.checked));
233+const specLog = document.getElementById("spec-log") as HTMLInputElement;
234+specLog.addEventListener("change", () => spec?.setLog(specLog.checked));
235+let noteEnergy = NaN;
236+const soundTraces = new TracePanel(
237+ document.getElementById("sound-chips")!,
238+ document.getElementById("sound-traces")!,
239+ [{ label: "Energy", unit: "J", log: true, series: [{ label: "E", value: () => noteEnergy }] }],
240+ 10, 960, soundCursor,
241+);
242+let spec: NoteSpectrogram | null = null;
243+let note: { x: Float32Array; energy: Float32Array; rate: number } | null = null;
244+let fed = 0; // samples of the note given to the displays
245+function feedDisplays(upTo: number) {
246+ if (!note || !spec) return;
247+ const end = Math.min(note.x.length, upTo);
248+ for (; fed < end; fed++) {
249+ spec.push(note.x[fed]);
250+ if (fed % DIAG_EVERY === 0) {
251+ noteEnergy = note.energy[fed / DIAG_EVERY];
252+ soundTraces.sample(fed / note.rate);
253+ }
254+ }
255+}
256+
257+const worker = new Worker(new URL("./render-worker.ts", import.meta.url), { type: "module" });
258+let requestId = 0;
259+let noteInfo = "";
260+let noteFile = "";
261+let renderStart = 0;
262+saveButton.addEventListener("click", () => {
263+ if (note) saveFile(wavBlob(note.x, note.rate), noteFile);
264+});
265+function renderNote() {
266+ if (!lastStrike || !soundPanel.open) return;
267+ const ctx = player.ensure();
268+ player.setPaused(!running);
269+ showPitch();
270+ const sig = signalByKey(signalSelect.value);
271+ const x = coord(soundX), y = coord(soundY);
272+ soundX.value = String(x);
273+ soundY.value = String(y);
274+ const req: RenderRequest = {
275+ id: ++requestId, P, state: lastStrike, speedup: speedup(), sampleRate: ctx.sampleRate,
276+ maxLength: 2 * predictedNoteLength(), signal: sig.key, x, y,
277+ };
278+ noteInfo = `${sig.label} (${x}, ${y}) m`;
279+ noteFile = `plate-ss-N${P.N}-${strikeLabel}-${energy(P, lastStrike).toPrecision(3)}J-${sig.key}x${x}y${y}m-${pitch().toFixed(0)}Hz.wav`;
280+ soundStatus.textContent = "computing…";
281+ soundHint.hidden = true;
282+ player.stop();
283+ renderStart = performance.now();
284+ worker.postMessage(req);
285+}
286+worker.onmessage = (e: MessageEvent<RenderResult>) => {
287+ if (e.data.id !== requestId) return; // superseded
288+ const rate = player.sampleRate;
289+ const x = e.data.samples;
290+ Player.normalize(x, rate);
291+ note = { x, energy: e.data.energy, rate };
292+ fed = 0;
293+ const span = x.length / rate;
294+ const n = Math.min(16384, 2 ** Math.ceil(Math.log2((2048 * 30) / Math.min(30, pitch()))));
295+ spec = new NoteSpectrogram(
296+ { dt: 1 / rate, n, nfft: 2 * n, hop: n / 4, refFreqs: modePitches(), width: 960, log: specLog.checked },
297+ soundCursor,
298+ );
299+ spec.start(span);
300+ document.getElementById("spectrogram")!.replaceChildren(spec.canvas);
301+ soundTraces.clear();
302+ soundTraces.setFixed({ t0: 0, span });
303+ const secs = (performance.now() - renderStart) / 1000;
304+ soundStatus.textContent = `${noteInfo} · ${fmtDuration(span)} · computed in ${secs.toFixed(1)} s`;
305+ player.play(x);
306+ saveButton.disabled = false;
307+ player.setPaused(!running);
308+};
309+replayButton.addEventListener("click", renderNote);
310+
311+// Set new parameters (for REF_N masses per side) and N: back to rest.
312+function setParams(next: Params, N: number) {
313+ base = next;
314+ P = refine(base, N);
315+ modes = new Modes(P);
316+ H = stepFor(P, modes);
317+ s = newState(P.N);
318+ verlet = new Verlet(P, s);
319+ E0 = 0;
320+ t = 0;
321+ ampl = 0;
322+ nInput.value = String(P.N);
323+ pauseButton.disabled = true;
324+ setRunning(false);
325+ pauseButton.textContent = "Pause";
326+ lastStrike = null;
327+ replayButton.disabled = true;
328+ saveButton.disabled = true;
329+ note = null;
330+ player.stop();
331+ for (const el of [strikeXInput, strikeYInput, soundX, soundY]) {
332+ el.max = String(plateSide(P));
333+ el.value = String(coord(el));
334+ }
335+ showPitch();
336+ showParams();
337+}
338+setParams(base, 8);
339+
340+nInput.addEventListener("change", () => {
341+ const N = Math.min(32, Math.max(2, Math.round(Number(nInput.value) || P.N)));
342+ setParams(base, N);
343+});
344+document.getElementById("param-apply")!.addEventListener("click", () => {
345+ const v = Object.fromEntries(Object.entries(paramInputs).map(([key, el]) => [key, Number(el.value)]));
346+ const mv = { M: Number(malletInputs.M.value), K: Number(malletInputs.K.value) };
347+ const ok = Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop"].includes(key) && x === 0))) &&
348+ Object.values(mv).every((x) => Number.isFinite(x) && x > 0);
349+ if (!ok) {
350+ paramStatus.textContent = "parameters must be positive (τ, τ top may be 0)";
351+ return;
352+ }
353+ paramStatus.textContent = "";
354+ mallet = mv;
355+ setParams({ ...base, ...v }, P.N);
356+});
357+
358+const TEXT = "#23262b", MUTED = "#7a7d82";
359+const MONO = "12px ui-monospace, SFMono-Regular, Menlo, monospace";
360+
361+// The plate seen from above, in the canvas: a square of side `size` px at
362+// (x0, y0), y up.
363+function geometry() {
364+ const W = plateCanvas.clientWidth, Ht = plateCanvas.clientHeight;
365+ const size = Math.min(W - 24, Ht - 48);
366+ return { W, Ht, size, x0: (W - size) / 2, y0: 26 };
367+}
368+// Height map: u interpolated bilinearly between the nodes (0 at the edge),
369+// red up and blue down, white at rest, scaled to the recent largest |u|;
370+// thin grid lines through the masses; the strike point (cross) and the
371+// listening point (circle).
372+const RES = 160; // pixels of the height map per side
373+const heat = document.createElement("canvas");
374+heat.width = heat.height = RES;
375+const heatCtx = heat.getContext("2d")!;
376+const heatImage = heatCtx.createImageData(RES, RES);
377+function drawPlate() {
378+ const ctx = fit(plateCanvas);
379+ const g = geometry();
380+ ctx.clearRect(0, 0, g.W, g.Ht);
381+ const { N } = P, L = plateSide(P);
382+ const at = (i: number, j: number) => (i < 0 || j < 0 || i >= N || j >= N ? 0 : s.u[i + N * j]);
383+ const scale = ampl > 0 ? 1 / ampl : 0;
384+ const d = heatImage.data;
385+ for (let py = 0; py < RES; py++) {
386+ const fy = ((RES - 1 - py + 0.5) / RES) * (N + 1) - 1; // node index, y up
387+ const j0 = Math.floor(fy), ay = fy - j0;
388+ for (let px = 0; px < RES; px++) {
389+ const fx = ((px + 0.5) / RES) * (N + 1) - 1;
390+ const i0 = Math.floor(fx), ax = fx - i0;
391+ const u = (1 - ax) * (1 - ay) * at(i0, j0) + ax * (1 - ay) * at(i0 + 1, j0) + (1 - ax) * ay * at(i0, j0 + 1) + ax * ay * at(i0 + 1, j0 + 1);
392+ const v = Math.max(-1, Math.min(1, u * scale)), q = 4 * (px + RES * py);
393+ // white to red (up) or blue (down)
394+ const c = v >= 0 ? [255, 255 - 200 * v, 255 - 210 * v] : [255 + 210 * v, 255 + 150 * v, 255 + 40 * v];
395+ d[q] = c[0];
396+ d[q + 1] = c[1];
397+ d[q + 2] = c[2];
398+ d[q + 3] = 255;
399+ }
400+ }
401+ heatCtx.putImageData(heatImage, 0, 0);
402+ ctx.imageSmoothingEnabled = true;
403+ ctx.drawImage(heat, g.x0, g.y0, g.size, g.size);
404+ // Grid lines through the masses, and the supported edge.
405+ const px = (x: number) => g.x0 + (x / L) * g.size, py = (y: number) => g.y0 + g.size - (y / L) * g.size;
406+ ctx.strokeStyle = "rgba(0, 0, 0, 0.06)";
407+ ctx.lineWidth = 1;
408+ ctx.beginPath();
409+ for (let i = 1; i <= N; i++) {
410+ const z = i * P.a;
411+ ctx.moveTo(px(z), py(0));
412+ ctx.lineTo(px(z), py(L));
413+ ctx.moveTo(px(0), py(z));
414+ ctx.lineTo(px(L), py(z));
415+ }
416+ ctx.stroke();
417+ ctx.strokeStyle = MUTED;
418+ ctx.lineWidth = 2;
419+ ctx.strokeRect(g.x0, g.y0, g.size, g.size);
420+ // Strike point and listening point.
421+ const sx = px(coord(strikeXInput)), sy = py(coord(strikeYInput));
422+ ctx.strokeStyle = TEXT;
423+ ctx.lineWidth = 1.5;
424+ ctx.beginPath();
425+ ctx.moveTo(sx - 6, sy - 6);
426+ ctx.lineTo(sx + 6, sy + 6);
427+ ctx.moveTo(sx - 6, sy + 6);
428+ ctx.lineTo(sx + 6, sy - 6);
429+ ctx.stroke();
430+ ctx.beginPath();
431+ ctx.arc(px(coord(soundX)), py(coord(soundY)), 6, 0, 2 * Math.PI);
432+ ctx.stroke();
433+ // Readouts.
434+ ctx.fillStyle = TEXT;
435+ ctx.font = MONO;
436+ ctx.textAlign = "left";
437+ ctx.fillText(`t ${t.toFixed(2)} s · E ${fmtEnergy(energy(P, s))}`, 6, 16);
438+ ctx.fillStyle = MUTED;
439+ ctx.textAlign = "right";
440+ ctx.fillText(`energy error ${Math.abs(relErr()).toExponential(1)}`, g.W - 6, 16);
441+ ctx.font = "11px system-ui, sans-serif";
442+ ctx.fillText(`side ${fmtLength(L)}${ampl > 0 ? ` · colour full at ±${fmtLength(ampl)}` : ""}`, g.W - 6, g.Ht - 6);
443+ ctx.textAlign = "left";
444+}
445+
446+// Animation: the simulation advances `speed` seconds per second of real time.
447+let last = performance.now();
448+let carry = 0;
449+let playing = false;
450+function frame(now: number) {
451+ const dt = Math.min(0.1, (now - last) / 1000);
452+ last = now;
453+ if (running) {
454+ carry += Number(speedSelect.value) * dt;
455+ while (carry >= H) {
456+ verlet.step(s, H);
457+ carry -= H;
458+ t += H;
459+ }
460+ }
461+ let umax = 0;
462+ for (const u of s.u) umax = Math.max(umax, Math.abs(u));
463+ ampl = Math.max(umax, 0.998 * ampl);
464+ drawPlate();
465+ if (note && soundPanel.open) {
466+ const at = player.position();
467+ if (at !== null) feedDisplays(Math.round(at * note.rate));
468+ else if (playing) feedDisplays(note.x.length);
469+ if (!soundCursor.pinned && (at !== null || playing)) soundCursor.t = at;
470+ playing = at !== null;
471+ if (!specFig.hidden) spec?.draw();
472+ soundTraces.draw();
473+ }
474+ requestAnimationFrame(frame);
475+}
476+requestAnimationFrame(frame);
src/plate/plate.tsadded+252−0View file
@@ -0,0 +1,252 @@
1+// A square plate as N x N equal masses on vertical rails at the nodes of a
2+// grid, spacing a, moving up and down only, by u (m), stored row by row:
3+// node (i, j) at index i + N j, at x = (i + 1) a, y = (j + 1) a. The edge
4+// nodes (i or j = -1 or N) are supports at u = 0: the plate is simply
5+// supported (held at its edge, free to turn there). SI units. Linear for now:
6+// - bending: kappa c^2 / 2 at each node, with c = (L u) the five-point
7+// Laplacian, 4 u(i, j) minus the four neighbours (kappa = D / a^2, D the
8+// flexural rigidity). A simply supported edge has no bending term: with
9+// the mirror ghost u(-2, j) = -u(0, j), its Laplacian is 0. The force is
10+// F = -kappa L (L u).
11+// Normal modes: L is diagonal in the discrete sine basis, so the modes are
12+// exact products of sines, S_mn(i, j) = (2 / (N + 1)) sin(m pi (i + 1) / (N +
13+// 1)) sin(n pi (j + 1) / (N + 1)), m, n = 1..N, with omega_mn = sqrt(kappa /
14+// m) (l_m + l_n), l_m = 4 sin^2(m pi / (2 (N + 1))). For low modes l_m ~ (m
15+// pi / (N + 1))^2, so omega_mn ~ m^2 + n^2: the continuum plate's spectrum.
16+// Optional damping, as in the beam (../beam/beam.ts): drag -gamma p on each
17+// mass (all modes decay as exp(-t / tau)), and a frequency-dependent loss,
18+// a dashpot alongside each grid spring, force -c (L v), which damps mode mn
19+// at rate c (l_m + l_n) / 2m, set by tauTop at the top of the band.
20+// State: displacements u, momenta p (kg m/s), and W, the energy (J)
21+// dissipated so far (so E + W is conserved).
22+
23+export interface Params {
24+ N: number; // masses per side (N^2 in all)
25+ m: number; // kg, each mass
26+ kappa: number; // N/m, bending
27+ a: number; // m, rest spacing (for positions and drawing; it does not enter the equations)
28+ tau: number; // s, amplitude decay time from drag (all modes); 0 = no damping
29+ tauTop: number; // s, amplitude decay time of the top of the band from the loss alone; 0 = off
30+}
31+
32+// Resolution: the page's parameters describe the plate at REF_N masses per
33+// side; refine() gives the same plate (same size, strike energies, decay
34+// times) with N. The spacing becomes a r, r = (REF_N + 1) / (N + 1), and
35+// each constant follows its continuum origin: m = rho h a^2 (x r^2), kappa =
36+// D / a^2 (/ r^2), the loss dashpots c = eta (unchanged, so a wave of given
37+// length decays at the same rate at any N). Drag does not change.
38+export const REF_N = 8;
39+export function refine(base: Params, N: number): Params {
40+ const r = (base.N + 1) / (N + 1);
41+ const P = { ...base, N, m: base.m * r * r, kappa: base.kappa / (r * r), a: base.a * r };
42+ // tauTop is set at the top of the band, which moves with N: keep c.
43+ if (base.tauTop > 0) P.tauTop = base.tauTop * (laplacianMax(base) / laplacianMax(P)) * r * r;
44+ return P;
45+}
46+
47+export const dragRate = (P: Params) => (P.tau > 0 ? 2 / P.tau : 0);
48+export const isDamped = (P: Params) => P.tau > 0 || P.tauTop > 0;
49+
50+// Eigenvalues of the one-dimensional second difference, l_m, m = 1..N.
51+const lam1 = (N: number, m: number) => 4 * Math.sin((m * Math.PI) / (2 * (N + 1))) ** 2;
52+// Largest eigenvalue of L: 2 l_N.
53+const laplacianMax = (P: Params) => 2 * lam1(P.N, P.N);
54+// Loss coefficient c (kg/s); 0 when off.
55+const lossCoefficient = (P: Params) => (P.tauTop > 0 ? (2 * P.m) / (P.tauTop * laplacianMax(P)) : 0);
56+
57+// (L x) into out, with x = 0 at the supports.
58+function laplacian(N: number, x: Float64Array, out: Float64Array): void {
59+ for (let j = 0, k = 0; j < N; j++) {
60+ for (let i = 0; i < N; i++, k++) {
61+ let s = 4 * x[k];
62+ if (i > 0) s -= x[k - 1];
63+ if (i < N - 1) s -= x[k + 1];
64+ if (j > 0) s -= x[k - N];
65+ if (j < N - 1) s -= x[k + N];
66+ out[k] = s;
67+ }
68+ }
69+}
70+
71+export interface State {
72+ u: Float64Array;
73+ p: Float64Array;
74+ W: number;
75+}
76+export const newState = (N: number): State => ({ u: new Float64Array(N * N), p: new Float64Array(N * N), W: 0 });
77+export const copyState = (s: State): State => ({ u: Float64Array.from(s.u), p: Float64Array.from(s.p), W: s.W });
78+
79+// Forces F = -kappa L (L u); c (if given) receives L u.
80+export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Array = new Float64Array(u.length)): void {
81+ laplacian(P.N, u, c);
82+ laplacian(P.N, c, F);
83+ for (let k = 0; k < F.length; k++) F[k] *= -P.kappa;
84+}
85+
86+export function kinetic(P: Params, s: State): number {
87+ let T = 0;
88+ for (const x of s.p) T += x * x;
89+ return T / (2 * P.m);
90+}
91+export function potential(P: Params, s: State): number {
92+ const c = new Float64Array(s.u.length);
93+ laplacian(P.N, s.u, c);
94+ let V = 0;
95+ for (const x of c) V += x * x;
96+ return (P.kappa * V) / 2;
97+}
98+export const energy = (P: Params, s: State) => kinetic(P, s) + potential(P, s);
99+
100+// Velocity Verlet, with drag and loss in half steps around it (Strang), as
101+// in the beam. Drag: p -> p exp(-gamma h / 2). Loss: explicit Euler steps of
102+// dp/dt = -(c / m) L p, p -> (I - x) p with x = (c dt / m) L. x is tiny here:
103+// its largest eigenvalue is 2 dt / tauTop (~3e-5 at the default step and
104+// tauTop), so each mode's decay rate is off by ~x / 2 relative, and since 1
105+// - x lies in [0, 1] the step never adds energy. The step is split into
106+// substeps if 2 dt / tauTop exceeds LOSS_X (only for a very short tauTop).
107+// (The beam solves Crank-Nicolson by the Thomas algorithm; on the grid that
108+// would need an iterative solve, several times the cost of the rest of the
109+// step, for no audible difference.) The kinetic energy each removes goes
110+// into W.
111+const LOSS_X = 0.01;
112+export class Verlet {
113+ private P: Params;
114+ private F: Float64Array;
115+ private c: number;
116+ private tmp: Float64Array;
117+ constructor(P: Params, s: State) {
118+ this.P = P;
119+ this.F = new Float64Array(P.N * P.N);
120+ this.tmp = new Float64Array(P.N * P.N);
121+ this.c = lossCoefficient(P);
122+ forces(P, s.u, this.F, this.tmp);
123+ }
124+ // Force on each mass at the current displacements.
125+ get force(): Float64Array {
126+ return this.F;
127+ }
128+ step(s: State, h: number): void {
129+ const { u, p } = s, F = this.F, n = u.length, m = this.P.m;
130+ const gamma = dragRate(this.P);
131+ if (gamma > 0) this.drag(s, Math.exp(-0.5 * gamma * h));
132+ if (this.c > 0) this.loss(s, h / 2);
133+ for (let k = 0; k < n; k++) p[k] += 0.5 * h * F[k];
134+ for (let k = 0; k < n; k++) u[k] += (h * p[k]) / m;
135+ forces(this.P, u, F, this.tmp);
136+ for (let k = 0; k < n; k++) p[k] += 0.5 * h * F[k];
137+ if (this.c > 0) this.loss(s, h / 2);
138+ if (gamma > 0) this.drag(s, Math.exp(-0.5 * gamma * h));
139+ }
140+ private loss(s: State, dt: number): void {
141+ const { N, m, tauTop } = this.P, p = s.p, Lp = this.tmp, n = p.length;
142+ const sub = Math.ceil((2 * dt) / tauTop / LOSS_X), b = (this.c * dt) / (m * sub);
143+ let T0 = 0;
144+ for (let k = 0; k < n; k++) T0 += p[k] * p[k];
145+ for (let q = 0; q < sub; q++) {
146+ laplacian(N, p, Lp);
147+ for (let k = 0; k < n; k++) p[k] -= b * Lp[k];
148+ }
149+ let T1 = 0;
150+ for (let k = 0; k < n; k++) T1 += p[k] * p[k];
151+ s.W += (T0 - T1) / (2 * m);
152+ }
153+ private drag(s: State, f: number): void {
154+ let p2 = 0;
155+ for (let k = 0; k < s.p.length; k++) {
156+ p2 += s.p[k] * s.p[k];
157+ s.p[k] *= f;
158+ }
159+ s.W += ((1 - f * f) * p2) / (2 * this.P.m);
160+ }
161+}
162+
163+// The normal modes, sorted by frequency (ties: by m, then n): mode numbers
164+// (m, n) and omega (rad/s).
165+export class Modes {
166+ N: number;
167+ count: number;
168+ m: Int32Array;
169+ n: Int32Array;
170+ omega: Float64Array;
171+ constructor(P: Params) {
172+ const N = P.N, w0 = Math.sqrt(P.kappa / P.m);
173+ const all: [number, number, number][] = [];
174+ for (let m = 1; m <= N; m++) for (let n = 1; n <= N; n++) all.push([w0 * (lam1(N, m) + lam1(N, n)), m, n]);
175+ all.sort((x, y) => x[0] - y[0] || x[1] - y[1] || x[2] - y[2]);
176+ this.N = N;
177+ this.count = all.length;
178+ this.omega = Float64Array.from(all, (x) => x[0]);
179+ this.m = Int32Array.from(all, (x) => x[1]);
180+ this.n = Int32Array.from(all, (x) => x[2]);
181+ }
182+ get omegaMax(): number {
183+ return this.omega[this.count - 1];
184+ }
185+}
186+
187+// Amplitude decay time (s) of each mode, in Modes order (Infinity if
188+// undamped): 1 / (1 / tau + (l_m + l_n) / (lambdaMax tauTop)). Exact: L is
189+// diagonal in the modes.
190+export function decayTimes(P: Params, modes: Modes): Float64Array {
191+ const lmax = laplacianMax(P);
192+ return Float64Array.from(modes.omega, (_, k) => {
193+ const q = lam1(P.N, modes.m[k]) + lam1(P.N, modes.n[k]);
194+ const rate = (P.tau > 0 ? 1 / P.tau : 0) + (P.tauTop > 0 ? q / (lmax * P.tauTop) : 0);
195+ return rate > 0 ? 1 / rate : Infinity;
196+ });
197+}
198+
199+// Integrator step (s): STEP / omega at the top of the band, as in the beam
200+// (linear, so the same at every energy).
201+const STEP = 0.2;
202+export const stepFor = (P: Params, modes: Modes) => STEP / modes.omegaMax;
203+
204+// Positions (m from the corner at the supports (-1, -1)): node (i, j) is at
205+// ((i + 1) a, (j + 1) a), the plate's side is (N + 1) a; nodeIndex gives the
206+// fractional index along one side at x.
207+export const plateSide = (P: Params) => (P.N + 1) * P.a;
208+export const nodeIndex = (P: Params, x: number) => x / P.a - 1;
209+
210+// The hammer's shape: a pyramid bump, (1 - |x_i - x| / w)(1 - |y_j - y| / w)
211+// on the masses it covers, (x, y) its centre and w (m) its half-width (at
212+// least one spacing); if it covers none, the nearest mass alone.
213+export function hammerShape(P: Params, x: number, y: number, w: number): Float64Array {
214+ const N = P.N, we = Math.max(w, P.a), phi = new Float64Array(N * N);
215+ let any = false;
216+ for (let j = 0; j < N; j++) {
217+ const fy = Math.max(0, 1 - Math.abs((j + 1) * P.a - y) / we);
218+ for (let i = 0; i < N; i++) {
219+ const v = fy * Math.max(0, 1 - Math.abs((i + 1) * P.a - x) / we);
220+ phi[i + N * j] = v;
221+ if (v > 0) any = true;
222+ }
223+ }
224+ if (!any) {
225+ const near = (z: number) => Math.min(N - 1, Math.max(0, Math.round(nodeIndex(P, z))));
226+ phi[near(x) + N * near(y)] = 1;
227+ }
228+ return phi;
229+}
230+
231+// A hammer blow at (x, y): velocity in the hammer's shape, energy E.
232+export function strikeAt(P: Params, x: number, y: number, w: number, E: number): State {
233+ const s = newState(P.N), phi = hammerShape(P, x, y, w);
234+ let S = 0;
235+ for (const v of phi) S += v * v;
236+ const scale = Math.sqrt((2 * P.m * E) / S);
237+ for (let k = 0; k < phi.length; k++) s.p[k] = phi[k] * scale;
238+ return s;
239+}
240+
241+// A per-node value v(k) at (x, y), bilinear between the four nearest nodes;
242+// the supports count as 0 (they do not move).
243+export function valueAt(P: Params, v: (k: number) => number, x: number, y: number): number {
244+ const N = P.N, L = plateSide(P);
245+ const fx = nodeIndex(P, Math.min(Math.max(x, 0), L)), fy = nodeIndex(P, Math.min(Math.max(y, 0), L));
246+ const i0 = Math.floor(fx), j0 = Math.floor(fy), ax = fx - i0, ay = fy - j0;
247+ const at = (i: number, j: number) => (i < 0 || j < 0 || i >= N || j >= N ? 0 : v(i + N * j));
248+ return (
249+ (1 - ax) * (1 - ay) * at(i0, j0) + ax * (1 - ay) * (ax ? at(i0 + 1, j0) : 0) +
250+ (1 - ax) * ay * (ay ? at(i0, j0 + 1) : 0) + ax * ay * (ax && ay ? at(i0 + 1, j0 + 1) : 0)
251+ );
252+}
src/plate/render-worker.tsadded+45−0View file
@@ -0,0 +1,45 @@
1+// Renders a note: the plate from a given state, sped up, sampled as audio
2+// (the loop and energy: ../common/note-render.ts). No Lyapunov shadow yet.
3+
4+import { Modes, Verlet, copyState, energy, isDamped, stepFor, type Params, type State } from "./plate";
5+import { renderNote, type Note } from "../common/note-render";
6+import { signalAt, signalByKey } from "./signals";
7+
8+export interface RenderRequest {
9+ id: number;
10+ P: Params;
11+ state: State;
12+ speedup: number;
13+ sampleRate: number;
14+ maxLength: number; // s of sound, at most
15+ signal: string; // signal key
16+ x: number; // m, listening point
17+ y: number;
18+}
19+
20+export interface RenderResult extends Note {
21+ id: number;
22+}
23+
24+self.onmessage = (e: MessageEvent<RenderRequest>) => {
25+ const r = e.data;
26+ const { P } = r;
27+ const s = copyState(r.state);
28+ const modes = new Modes(P);
29+ const verlet = new Verlet(P, s);
30+ const signal = signalByKey(r.signal);
31+ const note = renderNote(
32+ {
33+ energy: () => energy(P, s),
34+ stepFor: () => stepFor(P, modes),
35+ step: (h) => verlet.step(s, h),
36+ signal: () => signalAt(signal, P, s, verlet.force, r.x, r.y),
37+ shadow: () => {
38+ throw new Error("no Lyapunov shadow for the plate yet");
39+ },
40+ },
41+ { speedup: r.speedup, sampleRate: r.sampleRate, damped: isDamped(P), maxLength: r.maxLength, lyapunov: false },
42+ );
43+ const result: RenderResult = { id: r.id, ...note };
44+ (self as unknown as Worker).postMessage(result, [result.samples.buffer, result.energy.buffer, result.lambda.buffer]);
45+};
src/plate/signals.tsadded+23−0View file
@@ -0,0 +1,23 @@
1+// Quantities that can be heard, read at a point (x, y) of the plate (m from
2+// the corner), bilinear between the nearest masses, so the same listening
3+// point serves every N.
4+
5+import { valueAt, type Params, type State } from "./plate";
6+
7+export interface Signal {
8+ key: string;
9+ label: string;
10+ // Value of mass k now; F is the force on each mass at the current displacements.
11+ value: (P: Params, s: State, F: Float64Array, k: number) => number;
12+}
13+
14+export const SIGNALS: Signal[] = [
15+ { key: "v", label: "velocity at", value: (P, s, _F, k) => s.p[k] / P.m },
16+ { key: "u", label: "displacement at", value: (_P, s, _F, k) => s.u[k] },
17+ { key: "a", label: "acceleration at", value: (P, _s, F, k) => F[k] / P.m },
18+];
19+
20+export const signalByKey = (key: string) => SIGNALS.find((s) => s.key === key) ?? SIGNALS[0];
21+
22+export const signalAt = (sig: Signal, P: Params, s: State, F: Float64Array, x: number, y: number) =>
23+ valueAt(P, (k) => sig.value(P, s, F, k), x, y);
src/plate/style.cssadded+21−0View file
@@ -0,0 +1,21 @@
1+/* Styles for this app only, after ../common/style.css. */
2+
3+.row { display: flex; align-items: center; gap: 8px; flex-wrap: wrap; }
4+.row + .row { margin-top: 10px; }
5+canvas { display: block; width: 100%; }
6+[hidden] { display: none !important; }
7+#plate { aspect-ratio: 1; cursor: crosshair; }
8+.card-head { font-size: 12px; margin-bottom: 12px; display: flex; align-items: center; }
9+#sound-traces canvas { height: auto; margin-top: 12px; }
10+#sound-chips { margin-left: auto; }
11+input[type="number"] { width: 4.5em; font-family: var(--mono); }
12+#energy-value { min-width: 4.5em; text-align: right; }
13+details.panel { margin-top: 16px; border: 1px solid var(--line); padding: 16px; }
14+details.panel summary { display: flex; align-items: center; gap: 10px; }
15+#sound .row + .row { margin-top: 12px; }
16+#sound label { display: flex; align-items: center; gap: 8px; }
17+#sound .grow { flex: 1; min-width: 220px; }
18+#sound-hint { margin: 14px 0 0; }
19+#alias { color: #a33; }
20+#spectrogram canvas { width: 100%; height: auto; border-radius: 4px; }
21+#param-form { display: grid; grid-template-columns: repeat(auto-fill, minmax(150px, 1fr)); gap: 8px 16px; }
src/plate/vite-env.d.tsadded+1−0View file
@@ -0,0 +1 @@
1+/// <reference types="vite/client" />
vite.config.tsmodified+1−0View file
@@ -19,6 +19,7 @@ export default defineConfig({
1919 "double-pendulum": page("double-pendulum/index.html"),
2020 "fput-chain": page("fput-chain/index.html"),
2121 beam: page("beam/index.html"),
22+ plate: page("plate/index.html"),
2223 },
2324 // Name the chunk holding src/common/ "common" (otherwise it takes the
2425 // name of one of its modules).