Beam step 2: Lyapunov exponent (real-time chart, note chart with energy); render loop shared in src/common/note-render.ts, used by fput-chain (output unchanged)
12 changed files+317−117
README.mdmodified+1−1View file
@@ -5,7 +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/common/`: code shared by the apps (help dialogs, chart cursor, WAV export, note player, 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.
8+- `src/common/`: code shared by the apps (help dialogs, chart cursor, WAV export, note player, note render loop, 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.
99 - `src/index.html`: landing page linking the apps.
1010
1111 Live site: https://chaos-apps--listening-for-chaos.vault1-sites.magland.org/ (built and deployed by `.github/workflows/pages.yml` on every push to main).
src/beam/README.mdmodified+5−4View file
@@ -1,10 +1,10 @@
11 # beam
22
3-Struck beam with bending stiffness, step 3 of "Listening for chaos" (after double-pendulum and fput-chain). Work in progress: step 1 (physics, animation, sped-up note with spectrogram).
3+Struck beam with bending stiffness, step 3 of "Listening for chaos" (after double-pendulum and fput-chain). Work in progress: step 1 (physics, animation, sped-up note with spectrogram), step 2 (Lyapunov exponent).
44
5-Current state: N masses (2 to 64, default 8) on vertical rails, moving sideways only, with tension k, bending κ and a stretching nonlinearity k3 (below). Ends (a select next to Masses; changing it returns the beam to rest): clamped–clamped (default), pinned–pinned, or cantilever (clamped at the left, free at the right). Strike gives energy E (log slider, 10 µJ to 100 J, default 2.51 J) as velocity, on one mass alone (default mass 2) or in the shape of one normal mode. The card draws the beam from the side, displacements scaled to fit (scale given), with rails as faint lines, a hatched wall for a clamp and a triangle for a pin, plus the time, energy and relative energy error. Speed 0.1x to 20x real time; Pause or space. Parameters (m, k, κ, k3, a, τ, stretching form) fold out with derived values (mode frequencies, ratios to mode 1, the nonlinear scale, integrator step); Apply returns the beam to rest.
5+Current state: N masses (2 to 64, default 8) on vertical rails, moving sideways only, with tension k, bending κ and a stretching nonlinearity k3 (below). Ends (a select next to Masses; changing it returns the beam to rest): clamped–clamped (default), pinned–pinned, or cantilever (clamped at the left, free at the right). Strike gives energy E (log slider, 10 µJ to 100 J, default 2.51 J) as velocity, on one mass alone (default mass 2) or in the shape of one normal mode. The card draws the beam from the side, displacements scaled to fit (scale given), with rails as faint lines, a hatched wall for a clamp and a triangle for a pin, plus the time, energy and relative energy error. A "Real time" card on the right has chart chips (`../common/traces.ts`): Lyapunov (on; the largest Lyapunov exponent over the last 10 s and since the strike, last 60 s shown, the first 5 s after a strike left out of the y range) and Energy error (off). The exponent is the two-trajectory method of `../common/lyapunov.ts` with the beam's shadow in `lyapunov.ts`: a hidden second beam 1e-8 away, displacements in units of A / sqrt(N) and momenta in units of m omega_1 A / sqrt(N), A the mode-1 amplitude at the nonlinear scale. Speed 0.1x to 20x real time; Pause or space. Parameters (m, k, κ, k3, a, τ, stretching form) fold out with derived values (mode frequencies, ratios to mode 1, the nonlinear scale, integrator step); Apply returns the beam to rest.
66
7-"Sped up · listen" (as in fput-chain): each Strike renders the same strike in a worker (`render-worker.ts`), sped up so the lowest mode sounds at the Pitch (5 to 1000 Hz, default 100), and plays it; the page starts muted. A damped note is the whole decay to -60 dB. A red warning shows when the top mode is above half the sample rate (it would alias). Signals (`signals.ts`): velocity (default), displacement or acceleration of mass j (default 2), or the sum of accelerations (the net force of the supports). Replay, Save WAV (named e.g. beam-clamped-local-N8-mass2-2.51J-v2-100Hz.wav), spectrogram over the whole note with dashed lines at the mode pitches. The status line gives the render time.
7+"Sped up · listen" (as in fput-chain): each Strike renders the same strike in a worker (`render-worker.ts`), sped up so the lowest mode sounds at the Pitch (5 to 1000 Hz, default 100), and plays it; the page starts muted. A damped note is the whole decay to -60 dB. A red warning shows when the top mode is above half the sample rate (it would alias). Signals (`signals.ts`): velocity (default), displacement or acceleration of mass j (default 2), or the sum of accelerations (the net force of the supports). Replay, Save WAV (named e.g. beam-clamped-local-N8-mass2-2.51J-v2-100Hz.wav), spectrogram over the whole note with dashed lines at the mode pitches. Chips: Spectrogram (on), Lyapunov (on; the note's exponent over the last 10 s of physical time, from a shadow run in the worker; the first 5 s of physical time left out of the y range; the shadow runs only while this chart is shown, as it doubles the cost) and Energy (log). The three span the whole note and share one cursor (playhead while playing; hover, click to pin). The render loop is shared with fput-chain (`../common/note-render.ts`). The status line gives the render time (default note 1.0 to 1.4 s with the shadow, 0.6 s without).
88
99 Live site: https://chaos-apps--listening-for-chaos.vault1-sites.magland.org/beam/ (part of the chaos-apps site).
1010
@@ -18,6 +18,7 @@ Live site: https://chaos-apps--listening-for-chaos.vault1-sites.magland.org/beam
1818 - Nonlinear scale: the energy of mode 1 alone at the amplitude where its stretching energy equals its linear energy (N = 8, local: clamped 176 mJ, pinned 8.7 mJ, cantilever 1.1 mJ).
1919 - Damping: drag proportional to mass, every mode's amplitude decays as exp(-t / τ); applied exactly in half steps around each Verlet step, dissipated energy in W.
2020 - Integrator: velocity Verlet, step 0.03 / omega_eff, omega_eff^2 = (Kmax + 12 k3 d^2) / m, Kmax the top eigenvalue of the stiffness matrix and d = (4 E / k3)^(1/4) (Kirchhoff: the k3 term over sqrt(nS)). Relative energy error <= 8e-5 over 200 mode-1 periods, all ends and both forms, 10 mJ to 100 J.
21+- Default note (clamped, local, 2.51 J on mass 2, τ 540 s, 100 Hz): λ ~0.27 /s at the strike (λ T1 ~1.3), falling into the noise of the 10 s window (±0.05 /s) by ~1.6 s of sound (E ~0.14 J), as the spectrogram resolves into lines. Kirchhoff at the same strike: ~0.04 /s, then noise. A 50 mJ strike: noise throughout. Cantilever: ~0.2 /s for the first ~0.3 s of its 1.3 s note.
2122 - First look (offline, N = 8, hammer on mass 2, no damping, two-trajectory Lyapunov exponent times the mode-1 period T1): local stretching is chaotic from ~1 J (clamped 0.57 at 1 J, 3.3 at 10 J; pinned 1.8 at 1 J; cantilever 1.2 at 0.1 J). Kirchhoff is much less chaotic at the same energy (clamped 0.04 at 1 J, 0.24 at 10 J; pinned ~0 up to 1 J, 0.09 at 10 J).
2223
23-Code: `beam.ts` (physics, modes, strikes), `main.ts` (page), `signals.ts`, `render-worker.ts`; shared: `../common/sound.ts`, `../common/spectrogram.ts`, `../common/wav.ts`, `../common/mute.ts`, `../common/format.ts`, `../common/canvas.ts`, `../common/help.ts`.
24+Code: `beam.ts` (physics, modes, strikes), `main.ts` (page), `signals.ts`, `render-worker.ts`, `lyapunov.ts` (shadow); shared: `../common/note-render.ts`, `../common/lyapunov.ts`, `../common/traces.ts`, `../common/cursor.ts`, `../common/sound.ts`, `../common/spectrogram.ts`, `../common/wav.ts`, `../common/mute.ts`, `../common/format.ts`, `../common/canvas.ts`, `../common/help.ts`.
src/beam/beam.tsmodified+18−7View file
@@ -134,6 +134,10 @@ export class Verlet {
134134 this.F = new Float64Array(P.N);
135135 forces(P, s.u, this.F);
136136 }
137+ // Call after changing u from outside.
138+ reset(s: State): void {
139+ forces(this.P, s.u, this.F);
140+ }
137141 // Force on each mass at the current displacements.
138142 get force(): Float64Array {
139143 return this.F;
@@ -254,16 +258,23 @@ export class Modes {
254258 }
255259 }
256260
257-// Energy scale of the nonlinearity: the energy of mode 1 alone at the
258-// amplitude where its stretching (quartic) energy equals its linear energy.
259-// Below it the beam rings as independent modes; above it the modes trade
260-// energy.
261-export function nonlinearScale(P: Params, modes: Modes): number {
261+// Scale of the nonlinearity: the amplitude A of mode 1 alone (its shape has
262+// unit norm) at which its stretching (quartic) energy equals its linear
263+// energy, and the energy of mode 1 at that amplitude. Below that energy the
264+// beam rings as independent modes; above it the modes trade energy.
265+function unitEnergies(P: Params, modes: Modes): { V2: number; V4: number } {
262266 const s = newState(P.N);
263267 s.u.set(modes.S[0]); // unit amplitude
264268 const V2 = potential({ ...P, k3: 0 }, s);
265- const V4 = potential(P, s) - V2;
266- return V4 > 0 ? (V2 * V2) / V4 : Infinity; // A^2 = V2 / V4 makes them equal
269+ return { V2, V4: potential(P, s) - V2 };
270+}
271+export function nonlinearAmplitude(P: Params, modes: Modes): number {
272+ const { V2, V4 } = unitEnergies(P, modes);
273+ return V4 > 0 ? Math.sqrt(V2 / V4) : Infinity; // V2 A^2 = V4 A^4
274+}
275+export function nonlinearScale(P: Params, modes: Modes): number {
276+ const { V2, V4 } = unitEnergies(P, modes);
277+ return V4 > 0 ? (V2 * V2) / V4 : Infinity; // V2 A^2
267278 }
268279
269280 // Integrator step (s) for motion with energy E: 0.03 / omega, with omega from
src/beam/index.htmlmodified+20−4View file
@@ -8,6 +8,7 @@
88 <body>
99 <main>
1010 <header><h1>Beam chaos</h1></header>
11+ <div class="grid">
1112 <section class="card">
1213 <canvas id="beam"></canvas>
1314 <div id="controls">
@@ -31,9 +32,9 @@
3132 <input id="n" type="number" min="2" max="64" step="1" value="8" />
3233 <span class="muted" style="margin-left: 12px">Ends</span>
3334 <select id="ends">
34- <option value="clamped" selected>clamped–clamped</option>
35- <option value="pinned">pinned–pinned</option>
36- <option value="cantilever">cantilever (clamped–free)</option>
35+ <option value="clamped" selected title="clamped at both ends">clamped</option>
36+ <option value="pinned" title="pinned at both ends">pinned</option>
37+ <option value="cantilever" title="clamped at the left, free at the right">cantilever</option>
3738 </select>
3839 <span class="muted" style="margin-left: 12px">Speed</span>
3940 <select id="speed">
@@ -69,6 +70,12 @@
6970 <div id="params"></div>
7071 </details>
7172 </section>
73+ <section class="card">
74+ <div class="card-head">Real time <button type="button" class="help" data-help="realtime" title="About the real-time charts" aria-label="About the real-time charts">?</button></div>
75+ <div id="toggles" class="chips"></div>
76+ <div id="traces"></div>
77+ </section>
78+ </div>
7279 <details class="card panel" id="sound" open>
7380 <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>
7481 <div class="row">
@@ -78,6 +85,9 @@
7885 <label><span class="label">Volume</span> <input id="volume" type="range" min="0" max="1" step="0.01" value="0.5" /></label>
7986 <button id="replay" type="button" disabled title="Play the last strike again with the current settings">Replay</button>
8087 <button id="save-wav" type="button" disabled title="Save the last note as a WAV file">Save WAV</button>
88+ <span id="sound-chips" class="chips">
89+ <label class="chip"><input id="spec-toggle" type="checkbox" checked /><span>Spectrogram</span></label>
90+ </span>
8191 </div>
8292 <div class="row">
8393 <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>
@@ -87,6 +97,7 @@
8797 <div class="row"><span id="sound-status" class="muted"></span></div>
8898 <p id="sound-hint" class="muted">Strike to hear the whole decay.</p>
8999 <figure id="spec-fig"><div id="spectrogram"></div></figure>
100+ <div id="sound-traces"></div>
90101 </details>
91102 </main>
92103 <dialog id="help">
@@ -104,10 +115,15 @@
104115 <p>Mass, tension k, bending κ (either may be 0, not both), stretching k3, the rest spacing a (drawing only), and the decay time τ of the amplitude (drag on the masses; 0 for none). Apply restarts the beam at rest.</p>
105116 <p>Stretching "local" stiffens each spring by its own stretch, as in the FPUT chain. "Kirchhoff" gives all springs one shared tension set by the total stretch, as for a beam whose ends are held apart. The two agree when every spring is stretched alike. The nonlinear scale below is the energy of mode 1 at which its stretching energy equals its linear energy.</p>
106117 </template>
118+ <template id="help-realtime">
119+ <h2>Real time</h2>
120+ <p>Charts of the beam on the left, on the same clock. Click a chip to show or hide one; hover to read values, click to pin the cursor.</p>
121+ <p>Lyapunov is the rate (1/s) at which nearby motions separate, measured with a hidden second beam kept a tiny distance away: positive for chaos, near 0 for regular motion. Shown over the last 10 s and since the strike; the first 5 s are left out of the y range. Energy error should stay tiny.</p>
122+ </template>
107123 <template id="help-sound">
108124 <h2>Sped up · listen</h2>
109125 <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>
110- <p>Save WAV saves the last note as heard (16-bit mono, normalized). Signal chooses what is heard. The spectrogram spans the whole note, with dashed lines at the normal-mode pitches; a warning shows when the top mode is above half the sample rate and would alias.</p>
126+ <p>Save WAV saves the last note as heard (16-bit mono, normalized). Signal chooses what is heard. The charts span the whole note: the spectrogram, with dashed lines at the normal-mode pitches, the Lyapunov exponent over the last 10 s of the beam's own time, and the energy. A warning shows when the top mode is above half the sample rate and would alias.</p>
111127 </template>
112128 <script type="module" src="./main.ts"></script>
113129 </body>
src/beam/lyapunov.tsadded+42−0View file
@@ -0,0 +1,42 @@
1+// The beam's Lyapunov shadow (see ../common/lyapunov.ts): a second beam kept
2+// D0 from the beam.
3+//
4+// Distance: displacements in units of ur = A / sqrt(N), with A the mode-1
5+// amplitude at which stretching matches the linear energy (the size of a
6+// nonlinear motion, per mass), and momenta in units of m omega_1 ur (the
7+// momentum of mode 1 at that size), so both count comparably.
8+
9+import { D0, type Shadow } from "../common/lyapunov";
10+import { Verlet, copyState, nonlinearAmplitude, type Modes, type Params, type State } from "./beam";
11+
12+// A shadow of the beam's state s, for Lyapunov.start. s is read at every
13+// step, so it must stay the beam's state object until the next start.
14+export function shadowOf(P: Params, modes: Modes, s: State): Shadow {
15+ const ur = nonlinearAmplitude(P, modes) / Math.sqrt(P.N), pr = P.m * modes.omega[0] * ur;
16+ const sh = copyState(s);
17+ // Initial separation: equal parts in all 2N coordinates.
18+ const e = D0 / Math.sqrt(2 * P.N);
19+ for (let j = 0; j < P.N; j++) {
20+ sh.u[j] += e * ur;
21+ sh.p[j] += e * pr;
22+ }
23+ const verlet = new Verlet(P, sh);
24+
25+ return {
26+ step: (h) => verlet.step(sh, h),
27+ separation: () => {
28+ let d2 = 0;
29+ for (let j = 0; j < s.u.length; j++) {
30+ d2 += ((sh.u[j] - s.u[j]) / ur) ** 2 + ((sh.p[j] - s.p[j]) / pr) ** 2;
31+ }
32+ return Math.sqrt(d2);
33+ },
34+ pullBack: (r) => {
35+ for (let j = 0; j < s.u.length; j++) {
36+ sh.u[j] = s.u[j] + (sh.u[j] - s.u[j]) * r;
37+ sh.p[j] = s.p[j] + (sh.p[j] - s.p[j]) * r;
38+ }
39+ verlet.reset(sh); // forces at the moved displacements
40+ },
41+ };
42+}
src/beam/main.tsmodified+91−7View file
@@ -5,13 +5,18 @@ import {
55 type Ends, type Params, type State, type Stretch,
66 } from "./beam";
77 import { fit } from "../common/canvas";
8+import { newCursor } from "../common/cursor";
89 import { fmtDuration, fmtEnergy, fmtLength } from "../common/format";
910 import { setupHelp } from "../common/help";
11+import { Lyapunov } from "../common/lyapunov";
1012 import { setupMute } from "../common/mute";
13+import { DIAG_EVERY, END_RATIO, LYAP_WINDOW, UNDAMPED_LENGTH } from "../common/note-render";
1114 import { Player } from "../common/sound";
1215 import { NoteSpectrogram } from "../common/spectrogram";
16+import { TracePanel } from "../common/traces";
1317 import { saveFile, wavBlob } from "../common/wav";
14-import { END_RATIO, UNDAMPED_LENGTH, type RenderRequest, type RenderResult } from "./render-worker";
18+import { shadowOf } from "./lyapunov";
19+import type { RenderRequest, RenderResult } from "./render-worker";
1520 import { SIGNALS, signalByKey } from "./signals";
1621
1722 setupHelp();
@@ -31,6 +36,7 @@ const pauseButton = document.getElementById("pause") as HTMLButtonElement;
3136 // tau: amplitude decay time from drag; 540 s makes the default strike ring
3237 // down in about 8 s of sound at the default pitch.
3338 let P: Params = { N: 8, m: 1, k: 0, kappa: 20, k3: 4000, stretch: "local", ends: "clamped", a: 0.5, tau: 540 };
39+const SAMPLE_DT = 0.01; // s between plotted samples
3440
3541 // The beam. It sits at rest until the first Strike; changing parameters puts
3642 // it back at rest.
@@ -85,14 +91,43 @@ document.getElementById("strike")!.addEventListener("click", () => {
8591 showDerived();
8692 t = 0;
8793 ampl = 0;
94+ lyap.start(shadowOf(P, modes, s), t);
95+ traces.clear();
96+ traces.sample(t);
8897 replayButton.disabled = false;
8998 renderNote();
9099 pauseButton.disabled = false;
91100 setRunning(true);
92101 });
93102
103+// Energy balance: energy plus what drag has dissipated, against the strike's.
94104 const relErr = () => (E0 > 0 ? (energy(P, s) + s.W - E0) / E0 : 0);
95105
106+// Largest Lyapunov exponent of the animated beam, measured from each strike
107+// (two-trajectory method, shadow in lyapunov.ts): over the last 10 s and
108+// since the strike.
109+const lyap = new Lyapunov(LYAP_WINDOW);
110+// Time after a strike left out of the Lyapunov charts' y range: the early
111+// values swing widely before they settle.
112+const LYAP_SKIP = 5;
113+
114+// Real-time charts of the beam, sampled every SAMPLE_DT.
115+const traces = new TracePanel(
116+ document.getElementById("toggles")!,
117+ document.getElementById("traces")!,
118+ [
119+ {
120+ label: "Lyapunov", unit: "1/s", window: 60, skip: () => LYAP_SKIP, shown: true,
121+ series: [
122+ { label: "λ 10 s", value: () => lyap.windowed() },
123+ { label: "λ since strike", value: () => lyap.mean() },
124+ ],
125+ },
126+ { label: "Energy error", series: [{ label: "energy balance error", value: relErr }] },
127+ ],
128+ 10,
129+);
130+
96131 // Parameters: numbers in the form, stretching form in a select; Apply
97132 // restarts at rest. N and the ends are set above the form.
98133 const paramInputs = {
@@ -184,13 +219,41 @@ soundPanel.addEventListener("toggle", () => {
184219 if (!soundPanel.open) player.stop();
185220 });
186221
222+// Displays of the note: spectrogram, Lyapunov exponent and energy, sharing
223+// one time cursor (the playhead while playing; hover / click to pin).
224+const soundCursor = newCursor();
225+const specFig = document.getElementById("spec-fig")!;
226+const specToggle = document.getElementById("spec-toggle") as HTMLInputElement;
227+specToggle.addEventListener("change", () => (specFig.hidden = !specToggle.checked));
228+let noteEnergy = NaN, noteLambda = NaN;
229+const soundTraces = new TracePanel(
230+ document.getElementById("sound-chips")!,
231+ document.getElementById("sound-traces")!,
232+ [
233+ {
234+ label: "Lyapunov", unit: "1/s", shown: true,
235+ skip: () => LYAP_SKIP / speedup(), // the same physical time as in the animation, in s of sound
236+ series: [{ label: "λ 10 s", value: () => noteLambda }],
237+ },
238+ // Log scale: the decay is exponential, a straight line.
239+ { label: "Energy", unit: "J", log: true, series: [{ label: "E", value: () => noteEnergy }] },
240+ ],
241+ 10, 960, soundCursor,
242+);
187243 let spec: NoteSpectrogram | null = null;
188-let note: { x: Float32Array; rate: number } | null = null;
189-let fed = 0; // samples of the note given to the spectrogram
244+let note: { x: Float32Array; energy: Float32Array; lambda: Float32Array; rate: number } | null = null;
245+let fed = 0; // samples of the note given to the displays
190246 function feedDisplays(upTo: number) {
191247 if (!note || !spec) return;
192248 const end = Math.min(note.x.length, upTo);
193- for (; fed < end; fed++) spec.push(note.x[fed]);
249+ for (; fed < end; fed++) {
250+ spec.push(note.x[fed]);
251+ if (fed % DIAG_EVERY === 0) {
252+ noteEnergy = note.energy[fed / DIAG_EVERY];
253+ noteLambda = note.lambda[fed / DIAG_EVERY];
254+ soundTraces.sample(fed / note.rate);
255+ }
256+ }
194257 }
195258
196259 const worker = new Worker(new URL("./render-worker.ts", import.meta.url), { type: "module" });
@@ -213,6 +276,9 @@ function renderNote() {
213276 const req: RenderRequest = {
214277 id: ++requestId, P, state: lastStrike, speedup: speedup(), sampleRate: ctx.sampleRate,
215278 maxLength: 2 * predictedNoteLength(), signal: sig.key, j: j - 1,
279+ // The note's Lyapunov shadow doubles the computation, so it runs only
280+ // while its chart is shown.
281+ lyapunov: soundTraces.shownLabels().includes("Lyapunov"),
216282 };
217283 noteInfo = `${sig.label}${sig.perMass ? ` ${j}` : ""}`;
218284 noteFile =
@@ -229,15 +295,20 @@ worker.onmessage = (e: MessageEvent<RenderResult>) => {
229295 const rate = player.sampleRate;
230296 const x = e.data.samples;
231297 Player.normalize(x, rate);
232- note = { x, rate };
298+ note = { x, energy: e.data.energy, lambda: e.data.lambda, rate };
233299 fed = 0;
234300 const span = x.length / rate;
235301 // Spectrogram window: 2048 samples down to a lowest-mode pitch of 30 Hz;
236302 // below that longer, in proportion (as in the FPUT chain).
237303 const n = Math.min(16384, 2 ** Math.ceil(Math.log2((2048 * 30) / Math.min(30, pitch()))));
238- spec = new NoteSpectrogram({ dt: 1 / rate, n, nfft: 2 * n, hop: n / 4, refFreqs: modePitches(), width: 960 });
304+ spec = new NoteSpectrogram(
305+ { dt: 1 / rate, n, nfft: 2 * n, hop: n / 4, refFreqs: modePitches(), width: 960 },
306+ soundCursor,
307+ );
239308 spec.start(span);
240309 document.getElementById("spectrogram")!.replaceChildren(spec.canvas);
310+ soundTraces.clear();
311+ soundTraces.setFixed({ t0: 0, span });
241312 const secs = (performance.now() - renderStart) / 1000;
242313 soundStatus.textContent = `${noteInfo} · ${fmtDuration(span)} · computed in ${secs.toFixed(1)} s`;
243314 player.play(x);
@@ -256,6 +327,8 @@ function setParams(next: Params) {
256327 E0 = 0;
257328 t = 0;
258329 ampl = 0;
330+ lyap.stop();
331+ traces.clear();
259332 nInput.value = String(P.N);
260333 indexInput.max = String(P.N);
261334 indexInput.value = String(strikeIndex());
@@ -389,6 +462,7 @@ function drawBeam() {
389462 // Animation: the simulation advances `speed` seconds per second of real time.
390463 let last = performance.now();
391464 let carry = 0; // simulated time owed but not yet stepped
465+let sinceSample = 0;
392466 let playing = false; // a note was playing at the last frame
393467 function frame(now: number) {
394468 const dt = Math.min(0.1, (now - last) / 1000);
@@ -397,20 +471,30 @@ function frame(now: number) {
397471 carry += Number(speedSelect.value) * dt;
398472 while (carry >= H) {
399473 verlet.step(s, H);
474+ lyap.step(H, t + H);
400475 carry -= H;
401476 t += H;
477+ sinceSample += H;
478+ if (sinceSample >= SAMPLE_DT) {
479+ sinceSample -= SAMPLE_DT;
480+ traces.sample(t);
481+ }
402482 }
403483 }
404484 let umax = 0;
405485 for (const u of s.u) umax = Math.max(umax, Math.abs(u));
406486 ampl = Math.max(umax, 0.998 * ampl);
407487 drawBeam();
488+ traces.draw();
408489 if (note && soundPanel.open) {
490+ // Feed the displays up to the audio clock; playhead unless pinned.
409491 const at = player.position();
410492 if (at !== null) feedDisplays(Math.round(at * note.rate));
411493 else if (playing) feedDisplays(note.x.length); // ended
494+ if (!soundCursor.pinned && (at !== null || playing)) soundCursor.t = at;
412495 playing = at !== null;
413- spec?.draw();
496+ if (!specFig.hidden) spec?.draw();
497+ soundTraces.draw();
414498 }
415499 requestAnimationFrame(frame);
416500 }
src/beam/render-worker.tsmodified+21−32View file
@@ -1,17 +1,11 @@
1-// Renders a note: the beam from a given state, sped up, sampled as audio.
2-// Each audio sample advances the physics by speedup / sampleRate seconds, in
3-// equal substeps no longer than stepFor the current energy (re-chosen every
4-// DIAG_EVERY samples, so the steps lengthen as a damped note decays). A
5-// damped note ends once its energy has fallen to END_RATIO of the start
6-// (-60 dB); an undamped one after UNDAMPED_LENGTH seconds.
1+// Renders a note: the beam from a given state, sped up, sampled as audio
2+// (the loop, energy and Lyapunov exponent: ../common/note-render.ts).
73
84 import { Modes, Verlet, copyState, energy, stepFor, type Params, type State } from "./beam";
5+import { renderNote, type Note } from "../common/note-render";
6+import { shadowOf } from "./lyapunov";
97 import { signalByKey } from "./signals";
108
11-export const END_RATIO = 1e-6;
12-export const UNDAMPED_LENGTH = 10; // s of sound
13-export const DIAG_EVERY = 256; // samples between energy checks
14-
159 export interface RenderRequest {
1610 id: number;
1711 P: Params;
@@ -21,35 +15,30 @@ export interface RenderRequest {
2115 maxLength: number; // s of sound, at most
2216 signal: string; // signal key
2317 j: number; // mass index (0-based) for per-mass signals
18+ lyapunov: boolean; // run the Lyapunov shadow
2419 }
2520
26-export interface RenderResult {
21+export interface RenderResult extends Note {
2722 id: number;
28- samples: Float32Array; // raw signal, in the signal's units
2923 }
3024
3125 self.onmessage = (e: MessageEvent<RenderRequest>) => {
3226 const r = e.data;
27+ const { P } = r;
3328 const s = copyState(r.state);
34- const modes = new Modes(r.P);
35- const verlet = new Verlet(r.P, s);
29+ const modes = new Modes(P);
30+ const verlet = new Verlet(P, s);
3631 const signal = signalByKey(r.signal);
37- const dtSample = r.speedup / r.sampleRate;
38- let sub = 1, h = dtSample;
39- const E0 = energy(r.P, s);
40- const nMax = Math.round((r.P.tau > 0 ? r.maxLength : UNDAMPED_LENGTH) * r.sampleRate);
41- const out = new Float32Array(nMax);
42- let n = 0;
43- for (; n < nMax; n++) {
44- if (n % DIAG_EVERY === 0) {
45- const E = energy(r.P, s);
46- if (r.P.tau > 0 && E <= END_RATIO * E0) break;
47- sub = Math.ceil(dtSample / stepFor(r.P, modes, E));
48- h = dtSample / sub;
49- }
50- out[n] = signal.value(r.P, s, verlet.force, r.j);
51- for (let k = 0; k < sub; k++) verlet.step(s, h);
52- }
53- const result: RenderResult = { id: r.id, samples: out.slice(0, n) };
54- (self as unknown as Worker).postMessage(result, [result.samples.buffer]);
32+ const note = renderNote(
33+ {
34+ energy: () => energy(P, s),
35+ stepFor: (E) => stepFor(P, modes, E),
36+ step: (h) => verlet.step(s, h),
37+ signal: () => signal.value(P, s, verlet.force, r.j),
38+ shadow: () => shadowOf(P, modes, s),
39+ },
40+ { speedup: r.speedup, sampleRate: r.sampleRate, damped: P.tau > 0, maxLength: r.maxLength, lyapunov: r.lyapunov },
41+ );
42+ const result: RenderResult = { id: r.id, ...note };
43+ (self as unknown as Worker).postMessage(result, [result.samples.buffer, result.energy.buffer, result.lambda.buffer]);
5544 };
src/beam/style.cssmodified+4−1View file
@@ -4,7 +4,10 @@
44 .row + .row { margin-top: 10px; }
55 canvas { display: block; width: 100%; }
66 [hidden] { display: none !important; }
7-#beam { aspect-ratio: 2.4; }
7+#beam { aspect-ratio: 1.4; }
8+.card-head { font-size: 12px; margin-bottom: 12px; display: flex; align-items: center; }
9+#traces canvas, #sound-traces canvas { height: auto; margin-top: 12px; }
10+#sound-chips { margin-left: auto; }
811 input[type="number"] { width: 4.5em; font-family: var(--mono); }
912 #energy-value { min-width: 4.5em; text-align: right; }
1013 details.panel { margin-top: 16px; border: 1px solid var(--line); padding: 16px; }
src/common/note-render.tsadded+77−0View file
@@ -0,0 +1,77 @@
1+// Renders a note in a worker: a struck system run from a given state, sped
2+// up, sampled as audio. Each audio sample advances the physics by
3+// speedup / sampleRate seconds, in equal substeps no longer than the
4+// system's step for the current energy (re-chosen every DIAG_EVERY samples,
5+// so the steps lengthen as a damped note decays). On request a Lyapunov
6+// shadow runs along (LYAP_WINDOW of physical time; it doubles the cost). The
7+// energy and the windowed exponent (NaN without the shadow) are recorded
8+// once per DIAG_EVERY samples. A damped note ends once its energy has fallen
9+// to END_RATIO of the start (-60 dB); an undamped one after UNDAMPED_LENGTH
10+// seconds.
11+//
12+// Each app supplies its system as a NoteSystem; anything else it records
13+// (mode energies, section crossings) goes in its step and diag hooks.
14+
15+import { Lyapunov, type Shadow } from "./lyapunov";
16+
17+export const END_RATIO = 1e-6;
18+export const UNDAMPED_LENGTH = 10; // s of sound
19+export const DIAG_EVERY = 256; // samples per diagnostic value (energy, exponent)
20+export const LYAP_WINDOW = 10; // s of physical time
21+
22+export interface NoteSystem {
23+ energy(): number;
24+ stepFor(E: number): number; // integrator step (s) at energy E
25+ step(h: number, n: number): void; // advance by h during sample n
26+ signal(): number; // value heard now
27+ shadow(): Shadow; // a fresh Lyapunov shadow of the current state
28+ diag?(k: number): void; // called at diagnostic sample k (sample k DIAG_EVERY)
29+}
30+
31+export interface NoteOptions {
32+ speedup: number;
33+ sampleRate: number;
34+ damped: boolean;
35+ maxLength: number; // s of sound, at most (damped notes)
36+ lyapunov: boolean; // run the Lyapunov shadow
37+}
38+
39+export interface Note {
40+ samples: Float32Array; // raw signal, in the signal's units
41+ energy: Float32Array; // energy (J) at samples 0, DIAG_EVERY, 2 DIAG_EVERY, ...
42+ lambda: Float32Array; // largest Lyapunov exponent over the last LYAP_WINDOW s (1/s), same samples
43+}
44+
45+export function renderNote(sys: NoteSystem, o: NoteOptions): Note {
46+ const dtSample = o.speedup / o.sampleRate;
47+ let sub = 1, h = dtSample;
48+ const E0 = sys.energy();
49+ const nMax = Math.round((o.damped ? o.maxLength : UNDAMPED_LENGTH) * o.sampleRate);
50+ const out = new Float32Array(nMax);
51+ const en = new Float32Array(Math.ceil(nMax / DIAG_EVERY));
52+ const lam = new Float32Array(en.length);
53+ const lyap = new Lyapunov(LYAP_WINDOW);
54+ let t = 0; // physical time (s)
55+ if (o.lyapunov) lyap.start(sys.shadow(), t);
56+ let n = 0;
57+ for (; n < nMax; n++) {
58+ if (n % DIAG_EVERY === 0) {
59+ const k = n / DIAG_EVERY;
60+ const E = sys.energy();
61+ en[k] = E;
62+ lam[k] = o.lyapunov ? lyap.windowed() : NaN;
63+ sys.diag?.(k);
64+ if (o.damped && E <= END_RATIO * E0) break;
65+ sub = Math.ceil(dtSample / sys.stepFor(E));
66+ h = dtSample / sub;
67+ }
68+ out[n] = sys.signal();
69+ for (let i = 0; i < sub; i++) {
70+ sys.step(h, n);
71+ t += h;
72+ if (o.lyapunov) lyap.step(h, t);
73+ }
74+ }
75+ const nDiag = Math.ceil(n / DIAG_EVERY);
76+ return { samples: out.slice(0, n), energy: en.slice(0, nDiag), lambda: lam.slice(0, nDiag) };
77+}
src/fput-chain/README.mdmodified+1−1View file
@@ -4,7 +4,7 @@ Struck FPUT-beta chain, step 2 of "Listening for chaos" (after double-pendulum-c
44
55 Current state: N masses (2 to 64, default 2) between two fixed walls, joined by springs with tension k d + k3 d^3. The chain sits at rest until Strike, which gives it energy E (log slider, 10 µJ to 100 J, default 2.51 J) as velocity, either on one mass j alone (default mass 1, a hammer blow at one end) or in the shape of one normal mode n. The left card animates the chain (top: to scale, rest spacing a; below: the same displacements drawn sideways as a profile, scaled to fit, with the scale given) with the time and the relative energy error. Speed: 0.1x to 20x real time (default 1x). The time and current energy show at the top left. The Ghost chip (as in the pendulum app) adds a second chain started from the same state with every mass moved by the same amount, 10^k / sqrt(N) m (total offset 10^k m, k = -12 to -1, default -10; moving only the struck mass would put the offset along the motion itself, which does not grow), run in lockstep without renormalization and drawn in translucent orange underneath; restarted on Strike, on "Restart ghost" (from the chain's current state), when switched on or when the offset changes; a readout gives the separation (root sum of squared displacement differences, m) and the time since the ghost started, and a "Separation" chart (log, last 60 s; offered only with the ghost on) plots it. A "Lyapunov" chart (last 60 s; the first 5 s after a strike left out of its y range) shows the largest Lyapunov exponent over the last 10 s and since the strike (`../common/lyapunov.ts`, two-trajectory Benettin method shared with the pendulum app, with the chain's shadow in `lyapunov.ts`: a hidden shadow chain 1e-8 away, displacements in units of d0 and momenta in units of sqrt(m k) d0, renormalized every step). Checked (N = 4, mass-1 strike, 300 s): without damping 0.1 J ~0.05 to 0.1 /s (weak chaos; N = 4 onset is near 3 k d0^2 = 75 mJ), 0.5 J ~0.5 /s, 2.51 J ~0.9 /s, 25 J ~1.7 /s; with tau 540 s, 0.1 J falls to ~0.01 /s and 2.51 J falls from ~1.0 to ~0.3 /s by 300 s. The slider is labelled "Strike energy". Parameters (m, k, k3, a, tau) fold out underneath with derived values (mode frequencies, d0, nonlinear energy scale, integrator step); Apply or changing N returns the chain to rest. The "Real time" card on the right has toggle chips for its charts. "Mode history" (off by default) is a scrolling heat map of each mode's share of the harmonic energy, log10(E_n / sum E) over 6 decades (inferno colors, mode 1 at the bottom), over the last 300 s (`LiveModeHistory` in `modehistory.ts`); the cascade of energy up the modes shows here (e.g. N = 32, mode-1 strike at 10 mJ: the top mode's share goes from ~1e-28 to ~1e-12 in 10 periods of mode 1 and ~3e-5 in 50 without damping; with tau 540 s the damping wins and the spread stalls near mode 17). "Mode bars" shows the current harmonic energy of each normal mode relative to E on a log scale (12 decades), with the number of modes sharing the energy (inverse participation, 1 / sum (E_n / sum E)^2). The other chips are charts over time (from `../common/traces.ts` and `../common/cursor.ts`, with per-group colors and a legend capped at 8 series): displacements, velocities, spring stretches, mode energies (log, last 60 s; values below 1e-12 of E left out), energies (E, T, V, W dissipated) and the energy balance error ((E + W - E0) / E0); others last 10 s; sampled every 10 ms; hover for values, click to pin. Mode bars, displacements and mode energies are on by default. Masses, springs and modes share one color per index (four fixed colors for up to 4, a blue-to-red sweep beyond).
66
7-"Sped up · listen" (open by default, below the cards; otherwise as in the pendulum app): the header says how much faster it runs ("the whole decay, 321× faster: 62 min of motion in 12 s"; without damping "321× faster than the chain above"). The page starts muted (the mute button is tinted red and reads "Muted"; muted notes are still computed and shown); the audio context is created at the first Strike; closing the panel stops the note. While it is open each Strike plays a note, the same strike rendered in a Web Worker (`render-worker.ts`) sped up so the lowest linear mode sounds at the Pitch (5 to 1000 Hz, default 100; below ~20 Hz the lowest modes are below hearing but show in the spectrogram; the readout gives the lowest and highest mode pitches). A damped note is the whole decay, ending when its energy has fallen to 1e-6 of the strike's (-60 dB; predicted length (tau/2) ln 1e6 of physical time, 11.6 s of sound by default, measured 11.3 s; capped at twice that); undamped notes last 10 s. "Strike to hear the whole decay." shows until the first note; the status line gives the signal and note length ("computing…" while rendering). Replay plays the last strike again with the current settings (to compare signals). Signal (`signals.ts`): velocity (default), displacement or acceleration of mass j (default 1, the end mass, which is never at a node), force on the left wall (the first spring's tension), or sum of accelerations (the net wall force, in which even modes cancel). In one mass's motion mode n is weighted by sqrt(E_n) / omega_n, sqrt(E_n) or sqrt(E_n) omega_n. Each note is made zero-mean, scaled to a 0.9 peak (all signals equally loud; Volume sets the level) and faded out over its last 20 ms (`../common/sound.ts`). Save WAV saves the last note as heard (16-bit mono PCM, peak at 0.9; `../common/wav.ts`; named e.g. fput-N4-mass1-2.51J-v1-100Hz.wav). Chips (as in the pendulum app, plus Modes): Spectrogram (on; `NoteSpectrogram` in `../common/spectrogram.ts`: 2048-sample Hann, zero-pad 4096, hop 512, for pitches down to 30 Hz; below, the window lengthens in proportion (up to 16384 samples at 5 Hz, zero-pad 2x, hop a quarter window) so neighbouring low modes stay resolved; dashed lines at the linear mode pitches, drawn only while at least 6 px apart, since crowded together their dashes merged into a solid on-off bar at the bottom), Modes (off by default; the note's mode-energy history, each mode's share of the energy over the whole note, from the render worker, one column per 256 samples; `ModeHistory` in `modehistory.ts`; a cascade that takes tens of mode-1 periods occupies only the first few percent of a whole-decay note, so Mode history in the Real time panel shows it better), Lyapunov (the note's largest Lyapunov exponent over the last 10 s of physical time, from a shadow run in the worker alongside the note; the first 5 s of physical time left out of the y range), Energy (log scale, so the decay is a straight line) and λ vs E (`../common/scatter.ts`: one point per 256 samples, accumulated over notes, log energy axis, colored by time into the note on a log scale, earlier notes faded, hover for values, Clear button). The first three span the whole note (strike at left, end at right) and fill in as the audio clock reaches each sample; one shared cursor (playhead while playing; hover for values, click to pin). Pause (or space) pauses the animation and the sound (audio context suspended). The note's Lyapunov shadow runs only while the Lyapunov chart or λ vs E is shown (it doubles the computation). Render time in the browser pane (with the Lyapunov shadow): ~0.6 s for N = 4, ~2 s for N = 32 (default strike). First look (N = 4, 2.51 J): λ ~1.3 /s at the strike, falling to ~0 by ~4 s of sound as the spectrogram resolves into lines; on λ vs E the note's points sit near 0 below ~0.1 J and rise to ~1.4 /s at the top.
7+"Sped up · listen" (open by default, below the cards; otherwise as in the pendulum app): the header says how much faster it runs ("the whole decay, 321× faster: 62 min of motion in 12 s"; without damping "321× faster than the chain above"). The page starts muted (the mute button is tinted red and reads "Muted"; muted notes are still computed and shown); the audio context is created at the first Strike; closing the panel stops the note. While it is open each Strike plays a note, the same strike rendered in a Web Worker (`render-worker.ts`, on the loop shared with the beam app in `../common/note-render.ts`) sped up so the lowest linear mode sounds at the Pitch (5 to 1000 Hz, default 100; below ~20 Hz the lowest modes are below hearing but show in the spectrogram; the readout gives the lowest and highest mode pitches). A damped note is the whole decay, ending when its energy has fallen to 1e-6 of the strike's (-60 dB; predicted length (tau/2) ln 1e6 of physical time, 11.6 s of sound by default, measured 11.3 s; capped at twice that); undamped notes last 10 s. "Strike to hear the whole decay." shows until the first note; the status line gives the signal and note length ("computing…" while rendering). Replay plays the last strike again with the current settings (to compare signals). Signal (`signals.ts`): velocity (default), displacement or acceleration of mass j (default 1, the end mass, which is never at a node), force on the left wall (the first spring's tension), or sum of accelerations (the net wall force, in which even modes cancel). In one mass's motion mode n is weighted by sqrt(E_n) / omega_n, sqrt(E_n) or sqrt(E_n) omega_n. Each note is made zero-mean, scaled to a 0.9 peak (all signals equally loud; Volume sets the level) and faded out over its last 20 ms (`../common/sound.ts`). Save WAV saves the last note as heard (16-bit mono PCM, peak at 0.9; `../common/wav.ts`; named e.g. fput-N4-mass1-2.51J-v1-100Hz.wav). Chips (as in the pendulum app, plus Modes): Spectrogram (on; `NoteSpectrogram` in `../common/spectrogram.ts`: 2048-sample Hann, zero-pad 4096, hop 512, for pitches down to 30 Hz; below, the window lengthens in proportion (up to 16384 samples at 5 Hz, zero-pad 2x, hop a quarter window) so neighbouring low modes stay resolved; dashed lines at the linear mode pitches, drawn only while at least 6 px apart, since crowded together their dashes merged into a solid on-off bar at the bottom), Modes (off by default; the note's mode-energy history, each mode's share of the energy over the whole note, from the render worker, one column per 256 samples; `ModeHistory` in `modehistory.ts`; a cascade that takes tens of mode-1 periods occupies only the first few percent of a whole-decay note, so Mode history in the Real time panel shows it better), Lyapunov (the note's largest Lyapunov exponent over the last 10 s of physical time, from a shadow run in the worker alongside the note; the first 5 s of physical time left out of the y range), Energy (log scale, so the decay is a straight line) and λ vs E (`../common/scatter.ts`: one point per 256 samples, accumulated over notes, log energy axis, colored by time into the note on a log scale, earlier notes faded, hover for values, Clear button). The first three span the whole note (strike at left, end at right) and fill in as the audio clock reaches each sample; one shared cursor (playhead while playing; hover for values, click to pin). Pause (or space) pauses the animation and the sound (audio context suspended). The note's Lyapunov shadow runs only while the Lyapunov chart or λ vs E is shown (it doubles the computation). Render time in the browser pane (with the Lyapunov shadow): ~0.6 s for N = 4, ~2 s for N = 32 (default strike). First look (N = 4, 2.51 J): λ ~1.3 /s at the strike, falling to ~0 by ~4 s of sound as the spectrogram resolves into lines; on λ vs E the note's points sit near 0 below ~0.1 J and rise to ~1.4 /s at the top.
88
99 Poincaré section (shown only for N = 2 and N = 3, below the sped-up panel, open; adapted from the pendulum app): its own energy slider (log, 10 µJ to 100 J, default 2.51 J), same parameters, no damping. N = 2: each time mass 1 passes its rest position moving right, mass 2's (u2, p2) is plotted. N = 3: only the symmetric family (u1 = u3, p1 = p3; only modes 1 and 3 move, since the cubic force couples odd modes only to odd modes), which is two-dimensional: each time mass 2 passes its rest position moving right, mass 1's (u1, p1) is plotted. On either section the springs are stretched by 0 or +-x, so V = vFactor (k x^2 / 2 + k3 x^4 / 4) (vFactor 2 or 4), the plotted momentum counts once (N = 2) or twice (N = 3) in the kinetic energy, the outline is where the crossing mass's momentum vanishes, and that momentum follows from E. 42 orbits seeded on a grid, 300 crossings each, in a worker (`poincare-core.ts`, the chain's section, on `../common/section.ts` and `../common/section-worker.ts`; `poincare.ts`, the panel's configuration for `../common/section-panel.ts`; integrator step stepFor(E), energy error < 5e-5; ~1 s per section offline, a few s in the pane), drawn as they arrive; per-orbit colors or Color by λ (each orbit's Lyapunov exponent over its run, log scale); hover to pick out an orbit with its λ; click inside the outline to add an orbit there, and with Click to strike (on) also strike the chain from that exact state (and play it, if the sound panel is open). Follow (as in the pendulum app): "real time" makes the section's energy track the chain's as it decays (recomputed after a 3% change, at most every 2 s, and at each strike); "sped up" (which opens the sped-up panel) shows the playing note's crossings (found in the render worker, for N = 3 only while exactly symmetric, and drawn as the audio clock reaches them) as a fading trail of the last 150, with the axes fitted to the note's current energy (refit every 2%) and each dot fitted from its own energy; nothing is computed live. "Precompute sections" (shown only with Follow "sped up" until complete) builds an atlas (`../common/atlas.ts`, configured in `atlas.ts`): sections at 32 energies log-spaced from 1 mJ to 100 J (a factor ~1.45 apart), 42 orbits x 150 crossings, highest first, in its own worker (~22 s in the pane for N = 3), stored in IndexedDB per N, m, k, k3 and deleted after 7 days; the sound view shows the atlas section nearest the current energy, fitted to its allowed region (exact in the linear regime, below ~25 mJ, where sections are scaled copies; an approximation across the gaps above). Without an atlas only the outline is drawn behind the trail. Setting the section's energy by hand pauses everything and clears the rings and the trail; Follow takes over again on Resume. The chain's own crossings show as dark rings when within 5% of the section's energy and, for N = 3, only while the chain is exactly symmetric (else the legend says "chain not symmetric, not shown"; strikes on mass 2 or on mode 1 or 3 are symmetric, mode shapes are made exactly mirror-symmetric for this). Findings: the two-mass chain appears regular at every energy tried (10 mJ to 100 J): nested closed curves and island chains, no chaotic sea; the grid orbits' λ falls as the run lengthens (at 10 J, median 0.045 /s over 300 crossings, 0.007 /s over 3000), as for regular orbits. A likely reason: in mode coordinates (u = S Q) the potential is (k/2) rho^2 + (k3/8) rho^4 with rho^2 = Q1^2 + 3 Q2^2, a function of one combination of the two modes. The N = 3 symmetric family is chaotic: grid orbits with λ > 0.1 /s: 0 of 42 at 25 mJ, 10 at 0.25 J (islands in a chaotic sea), 37 at 2.51 J. It is unstable to breaking the symmetry: an exactly symmetric start stays symmetric (the arithmetic is exactly mirror-symmetric), but an asymmetry of 1e-12 grows to full size within ~100 s at 0.25 to 25 J. The ghost and the Lyapunov shadow offset all masses equally, so after a symmetric strike they too stay symmetric and measure the exponent within the family.
1010
src/fput-chain/main.tsmodified+2−1View file
@@ -10,7 +10,8 @@ import { shadowOf } from "./lyapunov";
1010 import { LiveModeHistory, ModeHistory } from "./modehistory";
1111 import { newPoincarePanel, type PoincarePanel } from "./poincare";
1212 import { sectionDef } from "./poincare-core";
13-import { END_RATIO, DIAG_EVERY, UNDAMPED_LENGTH, type RenderRequest, type RenderResult } from "./render-worker";
13+import { END_RATIO, DIAG_EVERY, UNDAMPED_LENGTH } from "../common/note-render";
14+import type { RenderRequest, RenderResult } from "./render-worker";
1415 import { ScatterPlot } from "../common/scatter";
1516 import { SIGNALS, signalByKey } from "./signals";
1617 import { Player } from "../common/sound";
src/fput-chain/render-worker.tsmodified+35−59View file
@@ -1,25 +1,14 @@
1-// Renders a note: the chain from a given state, sped up, sampled as audio.
2-// Each audio sample advances the physics by speedup / sampleRate seconds, in
3-// equal substeps no longer than stepFor the current energy (re-chosen every
4-// DIAG_EVERY samples, so the steps lengthen as a damped note decays). On
5-// request a Lyapunov shadow runs along (as in the animation: 10 s window of
6-// physical time; it doubles the cost). The energy, the windowed exponent
7-// (NaN without the shadow) and the mode energies are sent with the samples,
8-// one value (or set) per DIAG_EVERY samples. A damped note ends
9-// once its energy has fallen to END_RATIO of the start (-60 dB); an undamped
10-// one after UNDAMPED_LENGTH seconds.
1+// Renders a note: the chain from a given state, sped up, sampled as audio
2+// (the loop is ../common/note-render.ts). Besides the samples, energy and
3+// Lyapunov exponent, the mode energies are sent, N values per diagnostic
4+// sample, and for N = 2 and 3 the Poincare section crossings.
115
126 import { Modes, Verlet, copyState, energy, stepFor, type Params, type State } from "./chain";
13-import { Lyapunov } from "../common/lyapunov";
7+import { renderNote } from "../common/note-render";
148 import { shadowOf } from "./lyapunov";
159 import { sectionDef } from "./poincare-core";
1610 import { signalByKey } from "./signals";
1711
18-export const END_RATIO = 1e-6;
19-export const UNDAMPED_LENGTH = 10; // s of sound
20-export const DIAG_EVERY = 256; // samples per diagnostic value (energy, exponent)
21-export const LYAP_WINDOW = 10; // s of physical time
22-
2312 export interface RenderRequest {
2413 id: number;
2514 P: Params;
@@ -45,53 +34,40 @@ export interface RenderResult {
4534
4635 self.onmessage = (e: MessageEvent<RenderRequest>) => {
4736 const r = e.data;
37+ const { P } = r;
4838 const s = copyState(r.state);
49- const verlet = new Verlet(r.P, s);
39+ const verlet = new Verlet(P, s);
5040 const signal = signalByKey(r.signal);
51- const dtSample = r.speedup / r.sampleRate;
52- let sub = 1, h = dtSample;
53- const E0 = energy(r.P, s);
54- const nMax = Math.round((r.P.tau > 0 ? r.maxLength : UNDAMPED_LENGTH) * r.sampleRate);
55- let out = new Float32Array(nMax);
56- const en = new Float32Array(Math.ceil(nMax / DIAG_EVERY));
57- const lam = new Float32Array(en.length);
58- const modes = new Modes(r.P);
59- const Ek = new Float64Array(r.P.N);
60- const modeE = new Float32Array(en.length * r.P.N);
61- const lyap = new Lyapunov(LYAP_WINDOW);
62- let t = 0; // physical time (s)
63- if (r.lyapunov) lyap.start(shadowOf(r.P, s), t);
64- const sec = sectionDef(r.P.N);
41+ const modes = new Modes(P);
42+ const Ek = new Float64Array(P.N);
43+ const modeE: number[] = [];
44+ const sec = sectionDef(P.N);
6545 const crossings: number[] = [];
66- let n = 0;
67- for (; n < nMax; n++) {
68- if (n % DIAG_EVERY === 0) {
69- const E = energy(r.P, s);
70- en[n / DIAG_EVERY] = E;
71- lam[n / DIAG_EVERY] = r.lyapunov ? lyap.windowed() : NaN;
72- modes.energies(s, Ek);
73- modeE.set(Ek, (n / DIAG_EVERY) * r.P.N);
74- if (r.P.tau > 0 && E <= END_RATIO * E0) break;
75- sub = Math.ceil(dtSample / stepFor(r.P, E));
76- h = dtSample / sub;
77- }
78- out[n] = signal.value(r.P, s, verlet.force, r.j);
79- for (let k = 0; k < sub; k++) {
80- const c = sec ? s.u[sec.cross] : 0, x = sec ? s.u[sec.plot] : 0, q = sec ? s.p[sec.plot] : 0;
81- verlet.step(s, h);
82- t += h;
83- if (r.lyapunov) lyap.step(h, t);
84- if (sec && c < 0 && s.u[sec.cross] >= 0 && (r.P.N !== 3 || (s.u[0] === s.u[2] && s.p[0] === s.p[2]))) {
85- const f = -c / (s.u[sec.cross] - c);
86- crossings.push(n, x + f * (s.u[sec.plot] - x), q + f * (s.p[sec.plot] - q));
87- }
88- }
89- }
90- out = out.slice(0, n);
91- const nDiag = Math.ceil(n / DIAG_EVERY);
46+ const note = renderNote(
47+ {
48+ energy: () => energy(P, s),
49+ stepFor: (E) => stepFor(P, E),
50+ step: (h, n) => {
51+ const c = sec ? s.u[sec.cross] : 0, x = sec ? s.u[sec.plot] : 0, q = sec ? s.p[sec.plot] : 0;
52+ verlet.step(s, h);
53+ if (sec && c < 0 && s.u[sec.cross] >= 0 && (P.N !== 3 || (s.u[0] === s.u[2] && s.p[0] === s.p[2]))) {
54+ const f = -c / (s.u[sec.cross] - c);
55+ crossings.push(n, x + f * (s.u[sec.plot] - x), q + f * (s.p[sec.plot] - q));
56+ }
57+ },
58+ signal: () => signal.value(P, s, verlet.force, r.j),
59+ shadow: () => shadowOf(P, s),
60+ diag: () => {
61+ modes.energies(s, Ek);
62+ for (const x of Ek) modeE.push(x);
63+ },
64+ },
65+ { speedup: r.speedup, sampleRate: r.sampleRate, damped: P.tau > 0, maxLength: r.maxLength, lyapunov: r.lyapunov },
66+ );
9267 const result: RenderResult = {
93- id: r.id, samples: out, energy: en.slice(0, nDiag), lambda: lam.slice(0, nDiag),
94- modeEnergy: modeE.slice(0, nDiag * r.P.N), crossings: Float32Array.from(crossings),
68+ id: r.id, ...note,
69+ modeEnergy: Float32Array.from(modeE.slice(0, note.energy.length * P.N)),
70+ crossings: Float32Array.from(crossings),
9571 };
9672 (self as unknown as Worker).postMessage(result, [
9773 result.samples.buffer, result.energy.buffer, result.lambda.buffer, result.modeEnergy.buffer, result.crossings.buffer,