Poincaré section for N = 2 and N = 3's symmetric family, with Follow (real time, sped up) and precomputed sections; default N = 2
9 changed files+990−11
src/fput-chain/README.mdmodified+4−2View file
@@ -2,10 +2,12 @@
22
33 Struck FPUT-beta chain, step 2 of "Listening for chaos" (after double-pendulum-chaos). Work in progress.
44
5-Current state: N masses (2 to 64, default 4) 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 5x 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 (`src/lyapunov.ts`, two-trajectory Benettin method as in the pendulum app: 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 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 `src/traces.ts` and `src/cursor.ts` of double-pendulum-chaos, 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).
5+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 5x 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 (`src/lyapunov.ts`, two-trajectory Benettin method as in the pendulum app: 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 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 `src/traces.ts` and `src/cursor.ts` of double-pendulum-chaos, 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
77 "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 (`src/render-worker.ts`) sped up so the lowest linear mode sounds at the Pitch (30 to 1000 Hz, default 100; 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 (`src/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 (`src/sound.ts`). Chips (as in the pendulum app): Spectrogram (on; pendulum's `NoteSpectrogram`: 2048-sample Hann, zero-pad 4096, hop 512; 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), 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 (`src/scatter.ts` from the pendulum app: 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). 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
9+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 (`src/poincare-core.ts`, `src/poincare-worker.ts`, `src/poincare.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 (`src/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.
10+
911 Live site: https://fput-chain-chaos--listening-for-chaos.vault1-sites.magland.org/ (built and deployed by `.github/workflows/pages.yml` on every push to main).
1012
1113 ```bash
@@ -22,4 +24,4 @@ npm run dev
2224 - Chaos for N = 4 (two-trajectory Lyapunov, 5000 / sqrt(k/m) s): none below E ~ 1 k d0^2; onset near 3 k d0^2; at the default 10 k d0^2, lambda ~0.1 sqrt(k/m) = 0.32 /s, energy shared by ~2.5 modes.
2325 - Symmetry: the quartic force couples odd modes only to odd modes, so a strike in an odd mode leaves the even modes at exactly zero (they stay off the bottom of the chart). A strike on one mass, other than the middle one, excites all modes.
2426
25-Code: `src/chain.ts` (physics, modes, strikes), `src/main.ts` (page), `src/traces.ts` and `src/cursor.ts` (charts), `src/signals.ts`, `src/render-worker.ts`, `src/sound.ts`, `src/spectrogram.ts`, `src/scatter.ts` (sound), `src/lyapunov.ts` (Lyapunov exponent, animation and sound).
27+Code: `src/chain.ts` (physics, modes, strikes), `src/main.ts` (page), `src/traces.ts` and `src/cursor.ts` (charts), `src/signals.ts`, `src/render-worker.ts`, `src/sound.ts`, `src/spectrogram.ts`, `src/scatter.ts` (sound), `src/lyapunov.ts` (Lyapunov exponent, animation and sound), `src/poincare*.ts` (section, N = 2 and 3).
src/fput-chain/atlas.tsadded+138−0View file
@@ -0,0 +1,138 @@
1+// Atlas of Poincare sections at a ladder of energies, so the section can
2+// follow a decaying note faster than sections can be computed (as in the
3+// double pendulum app). Computed on request in its own worker (highest
4+// energy first) and kept in the browser (IndexedDB) per parameter set,
5+// deleted after MAX_AGE. Falls back to memory only if storage is unavailable.
6+
7+import type { Params } from "./chain";
8+import { sectionDef, seedGrid } from "./poincare-core";
9+import type { FromPoincare, ToPoincare } from "./poincare-worker";
10+
11+export const ATLAS_SIZE = 32; // energies
12+const ATLAS_CROSSINGS = 150; // per orbit
13+const GRID = [6, 7];
14+const MAX_AGE = 7 * 24 * 3600 * 1000; // ms
15+const VERSION = "atlas-v1"; // change when the computation changes
16+const DB = "fput-chain-chaos", STORE = "atlas";
17+
18+export interface AtlasOrbit {
19+ points: Float32Array; // x, q pairs
20+ lambda: number;
21+}
22+export interface AtlasEntry {
23+ E: number;
24+ orbits: AtlasOrbit[];
25+}
26+interface Record {
27+ key: string;
28+ created: number;
29+ entries: AtlasEntry[];
30+}
31+
32+// Energies: log-spaced from 1 mJ to 100 J (the top of the energy sliders),
33+// a factor of ~1.45 apart. Below ~25 mJ (k d0^2 for the defaults) the motion
34+// is nearly linear and sections at different energies are scaled copies, so
35+// the lowest entry serves for any lower energy.
36+export function atlasEnergies(): number[] {
37+ const lo = 1e-3, hi = 100;
38+ return Array.from({ length: ATLAS_SIZE }, (_, i) => lo * (hi / lo) ** (i / (ATLAS_SIZE - 1)));
39+}
40+
41+const keyOf = (P: Params) => `${VERSION}|${P.N}|${P.m}|${P.k}|${P.k3}`;
42+
43+function openDb(): Promise<IDBDatabase> {
44+ return new Promise((resolve, reject) => {
45+ const req = indexedDB.open(DB, 1);
46+ req.onupgradeneeded = () => req.result.createObjectStore(STORE, { keyPath: "key" });
47+ req.onsuccess = () => resolve(req.result);
48+ req.onerror = () => reject(req.error);
49+ });
50+}
51+
52+function tx<T>(db: IDBDatabase, mode: IDBTransactionMode, fn: (s: IDBObjectStore) => IDBRequest<T>): Promise<T> {
53+ return new Promise((resolve, reject) => {
54+ const req = fn(db.transaction(STORE, mode).objectStore(STORE));
55+ req.onsuccess = () => resolve(req.result);
56+ req.onerror = () => reject(req.error);
57+ });
58+}
59+
60+export class Atlas {
61+ private db: Promise<IDBDatabase | null>;
62+ private memory = new Map<string, Record>();
63+ private worker: Worker | null = null;
64+ private job = 0; // increments to cancel a computation
65+
66+ constructor() {
67+ this.db = openDb()
68+ .then(async (db) => {
69+ // Delete expired atlases.
70+ const all = await tx<Record[]>(db, "readonly", (s) => s.getAll());
71+ for (const r of all) if (Date.now() - r.created > MAX_AGE) await tx(db, "readwrite", (s) => s.delete(r.key));
72+ return db;
73+ })
74+ .catch(() => null);
75+ }
76+
77+ // Cached entries for P (possibly incomplete), sorted by energy.
78+ async load(P: Params): Promise<AtlasEntry[]> {
79+ const key = keyOf(P);
80+ let r = this.memory.get(key);
81+ const db = await this.db;
82+ if (!r && db) r = await tx<Record | undefined>(db, "readonly", (s) => s.get(key)).catch(() => undefined);
83+ if (!r || Date.now() - r.created > MAX_AGE) return [];
84+ this.memory.set(key, r);
85+ return [...r.entries].sort((a, b) => a.E - b.E);
86+ }
87+
88+ cancel(): void {
89+ this.job++;
90+ }
91+
92+ // Compute the missing entries for P, highest energy first, calling onEntry
93+ // after each one is computed and stored.
94+ async compute(P: Params, onEntry: (entries: AtlasEntry[]) => void): Promise<void> {
95+ const job = ++this.job;
96+ const key = keyOf(P);
97+ const have = await this.load(P);
98+ const record: Record = this.memory.get(key) ?? { key, created: Date.now(), entries: have };
99+ this.memory.set(key, record);
100+ const todo = atlasEnergies()
101+ .filter((E) => !record.entries.some((e) => Math.abs(e.E - E) < 1e-9 * E))
102+ .sort((a, b) => b - a);
103+ for (const E of todo) {
104+ const entry = await this.section(P, E, job);
105+ if (!entry || job !== this.job) return;
106+ record.entries.push(entry);
107+ const db = await this.db;
108+ if (db) await tx(db, "readwrite", (s) => s.put(record)).catch(() => undefined);
109+ onEntry([...record.entries].sort((a, b) => a.E - b.E));
110+ }
111+ }
112+
113+ // One section at energy E, from the worker (null if cancelled).
114+ private section(P: Params, E: number, job: number): Promise<AtlasEntry | null> {
115+ if (!this.worker) this.worker = new Worker(new URL("./poincare-worker.ts", import.meta.url), { type: "module" });
116+ const worker = this.worker;
117+ const seeds = seedGrid(P, sectionDef(P.N)!, E, GRID[0], GRID[1]);
118+ const points: number[][] = seeds.map(() => []);
119+ const lambdas: number[] = seeds.map(() => NaN);
120+ return new Promise((resolve) => {
121+ worker.onmessage = (e: MessageEvent<FromPoincare>) => {
122+ if (job !== this.job) {
123+ worker.postMessage({ type: "stop" } satisfies ToPoincare);
124+ resolve(null);
125+ return;
126+ }
127+ for (const o of e.data.orbits) {
128+ for (const v of o.points) points[o.index].push(v);
129+ lambdas[o.index] = o.lambda;
130+ }
131+ if (e.data.finished) {
132+ resolve({ E, orbits: points.map((p, k) => ({ points: Float32Array.from(p), lambda: lambdas[k] })) });
133+ }
134+ };
135+ worker.postMessage({ type: "run", id: job, P, E, seeds, crossings: ATLAS_CROSSINGS } satisfies ToPoincare);
136+ });
137+ }
138+}
src/fput-chain/chain.tsmodified+5−0View file
@@ -150,6 +150,11 @@ export class Modes {
150150 this.omega[n - 1] = 2 * w0 * Math.sin((Math.PI * n) / (2 * (N + 1)));
151151 const row = new Float64Array(N);
152152 for (let j = 0; j < N; j++) row[j] = c * Math.sin((Math.PI * n * (j + 1)) / (N + 1));
153+ // Make the mirror symmetry exact (row[N-1-j] = (-1)^(n+1) row[j]; sin
154+ // rounds differently on either side), so odd-mode strikes keep odd
155+ // motions exactly odd.
156+ for (let j = 0; j < Math.floor(N / 2); j++) row[N - 1 - j] = n % 2 ? row[j] : -row[j];
157+ if (N % 2 && n % 2 === 0) row[(N - 1) / 2] = 0; // even modes have a node at the middle mass
153158 this.S.push(row);
154159 }
155160 }
src/fput-chain/index.htmlmodified+25−2View file
@@ -72,7 +72,10 @@
7272 #sound-hint { margin: 14px 0 0; }
7373 #spectrogram canvas { width: 100%; height: auto; border-radius: 4px; }
7474 #sound-traces canvas { height: auto; margin-top: 12px; }
75- #scatter canvas { width: 100%; height: auto; }
75+ #scatter canvas, #pc-plot canvas { width: 100%; height: auto; }
76+ #pc-energy { max-width: 360px; }
77+ #poincare .row + .row { margin-top: 10px; }
78+ #poincare label { display: flex; align-items: center; gap: 8px; }
7679 .scatter-actions { justify-content: flex-end; margin-top: 6px; }
7780 button.mute { min-width: 36px; height: 30px; padding: 0 9px; display: inline-flex; align-items: center; justify-content: center; gap: 5px; }
7881 button.mute svg { width: 18px; height: 18px; flex: none; }
@@ -126,7 +129,7 @@
126129 </div>
127130 <div class="row">
128131 <span class="muted">Masses</span>
129- <input id="n" type="number" min="2" max="64" step="1" value="4" />
132+ <input id="n" type="number" min="2" max="64" step="1" value="2" />
130133 <span class="muted" style="margin-left: 12px">Speed</span>
131134 <select id="speed">
132135 <option value="0.1">0.1×</option>
@@ -202,6 +205,26 @@
202205 <div class="row scatter-actions"><button id="scatter-clear" type="button">Clear</button></div>
203206 </figure>
204207 </details>
208+
209+ <details class="card panel" id="poincare" open hidden>
210+ <summary><span class="label">Poincaré section</span></summary>
211+ <div class="row">
212+ <span class="label">Energy</span>
213+ <input id="pc-energy" type="range" min="-5" max="2" step="0.01" value="0.4" aria-label="Section energy" />
214+ <span id="pc-energy-value" class="num"></span>
215+ <label title="Real time: the section's energy follows the chain's as it decays. Sped up: the sped-up run's crossings are shown as it plays, on the section at its current energy."><span class="label">Follow</span> <select id="pc-follow">
216+ <option value="off">off</option>
217+ <option value="realtime">real time</option>
218+ <option value="sound">sped up</option>
219+ </select></label>
220+ <button id="pc-atlas" type="button" title="Compute sections at 32 energies (1 mJ to 100 J), so Follow sped up can show the orbits behind the decaying note's trail. Kept in this browser for 7 days." hidden>Precompute sections</button>
221+ </div>
222+ <div class="row">
223+ <label class="chip"><input id="pc-lambda" type="checkbox" /><span>Color by λ</span></label>
224+ <label class="chip" title="Clicking a point also strikes the chain from that exact state, at this energy"><input id="pc-strike" type="checkbox" checked /><span>Click to strike</span></label>
225+ </div>
226+ <figure><div id="pc-plot"></div></figure>
227+ </details>
205228 </main>
206229 <script type="module" src="/src/main.ts"></script>
207230 </body>
src/fput-chain/main.tsmodified+55−5View file
@@ -3,6 +3,8 @@ import {
33 } from "./chain";
44 import { newCursor } from "./cursor";
55 import { Lyapunov } from "./lyapunov";
6+import { PoincarePanel } from "./poincare";
7+import { sectionDef } from "./poincare-core";
68 import { END_RATIO, DIAG_EVERY, UNDAMPED_LENGTH, type RenderRequest, type RenderResult } from "./render-worker";
79 import { ScatterPlot } from "./scatter";
810 import { SIGNALS, signalByKey } from "./signals";
@@ -26,7 +28,7 @@ const modeStats = document.getElementById("mode-stats")!;
2628 // tau: amplitude decay time from drag; 540 s barely shows in the animation
2729 // (~10% amplitude per minute) but makes a struck note ring down in about
2830 // 12 s of sound. 0 = no damping.
29-let P: Params = { N: 4, m: 1, k: 10, k3: 4000, a: 0.5, tau: 540 };
31+let P: Params = { N: 2, m: 1, k: 10, k3: 4000, a: 0.5, tau: 540 };
3032 const SAMPLE_DT = 0.01; // s between plotted samples
3133
3234 // The chain. It sits at rest until the first Strike; changing parameters
@@ -81,7 +83,13 @@ const strikeIndex = () => Math.min(P.N, Math.max(1, Math.round(Number(indexInput
8183 document.getElementById("strike")!.addEventListener("click", () => {
8284 const E = strikeEnergy();
8385 indexInput.value = String(strikeIndex());
84- s = kindSelect.value === "mode" ? strikeMode(P, modes, strikeIndex(), E) : strikeMass(P, strikeIndex(), E);
86+ strikeFrom(kindSelect.value === "mode" ? strikeMode(P, modes, strikeIndex(), E) : strikeMass(P, strikeIndex(), E));
87+});
88+
89+// Start the chain (and, with the sound panel open, a note) from a state: a
90+// strike, or a point clicked on the Poincare section.
91+function strikeFrom(state: State) {
92+ s = state;
8593 lastStrike = copyState(s);
8694 verlet = new Verlet(P, s);
8795 E0 = energy(P, s);
@@ -96,9 +104,11 @@ document.getElementById("strike")!.addEventListener("click", () => {
96104 traces.clear();
97105 modes.energies(s, Ek);
98106 traces.sample(t);
107+ poincare?.clearLive();
108+ poincare?.follow(E0, true);
99109 pauseButton.disabled = false;
100110 setRunning(true);
101-});
111+}
102112
103113 // The ghost: a second chain started from the same state with every mass
104114 // moved by the same tiny amount (total offset delta, so each by
@@ -370,7 +380,11 @@ const scatterFig = document.getElementById("scatter-fig")!;
370380 }
371381 document.getElementById("scatter-clear")!.addEventListener("click", () => scatter.clear());
372382 // The rendered note and how much of it the displays have been fed.
373-let note: { x: Float32Array; energy: Float32Array; lambda: Float32Array; rate: number } | null = null;
383+// Its section crossings go to the Poincare panel as the audio clock reaches
384+// them (next: index of the next crossing triple).
385+let note: {
386+ x: Float32Array; energy: Float32Array; lambda: Float32Array; crossings: Float32Array; next: number; rate: number;
387+} | null = null;
374388 let fed = 0;
375389 function feedDisplays(upTo: number) {
376390 if (!note || !spec) return;
@@ -383,6 +397,13 @@ function feedDisplays(upTo: number) {
383397 noteLambda = note.lambda[fed / DIAG_EVERY];
384398 soundTraces.sample(t);
385399 scatter.add(noteEnergy, noteLambda, t);
400+ poincare?.setNoteEnergy(noteEnergy);
401+ }
402+ const c = note.crossings;
403+ while (note.next < c.length && c[note.next] <= fed) {
404+ const atE = note.energy[Math.min(note.energy.length - 1, Math.floor(c[note.next] / DIAG_EVERY))];
405+ poincare?.addNoteDot(c[note.next + 1], c[note.next + 2], atE);
406+ note.next += 3;
386407 }
387408 }
388409 }
@@ -412,7 +433,8 @@ worker.onmessage = (e: MessageEvent<RenderResult>) => {
412433 const rate = player.sampleRate;
413434 const x = e.data.samples;
414435 Player.normalize(x, rate);
415- note = { x, energy: e.data.energy, lambda: e.data.lambda, rate };
436+ note = { x, energy: e.data.energy, lambda: e.data.lambda, crossings: e.data.crossings, next: 0, rate };
437+ poincare?.noteStarted(e.data.energy[0]);
416438 fed = 0;
417439 const span = x.length / rate;
418440 spec = new NoteSpectrogram(
@@ -457,9 +479,25 @@ function setParams(next: Params) {
457479 soundMass.value = String(Math.min(P.N, Number(soundMass.value) || 1));
458480 showPitch();
459481 buildTraces(shown);
482+ poincare?.setParams(P);
460483 showParams();
461484 }
485+// Poincare section (N = 2 only). Click to strike starts the chain from the
486+// clicked state at the section's energy.
487+let poincare: PoincarePanel | null = null;
462488 setParams(P);
489+poincare = new PoincarePanel(
490+ P,
491+ (state) => strikeFrom(state),
492+ () => {
493+ if (running) setRunning(false); // a hand-set section energy pauses everything
494+ },
495+ fmtEnergy,
496+);
497+// Follow "sped up" needs notes, so choosing it opens the sound panel.
498+document.getElementById("pc-follow")!.addEventListener("change", (e) => {
499+ if ((e.target as HTMLSelectElement).value === "sound") soundPanel.open = true;
500+});
463501
464502 nInput.addEventListener("change", () => {
465503 const N = Math.min(64, Math.max(2, Math.round(Number(nInput.value) || P.N)));
@@ -644,9 +682,19 @@ function frame(now: number) {
644682 const dt = Math.min(0.1, (now - last) / 1000);
645683 last = now;
646684 if (running) {
685+ const sec = sectionDef(P.N);
647686 carry += Number(speedSelect.value) * dt;
648687 while (carry >= H) {
688+ const c = sec ? s.u[sec.cross] : 0, x = sec ? s.u[sec.plot] : 0, q = sec ? s.p[sec.plot] : 0;
649689 verlet.step(s, H);
690+ // N = 2 and 3: the section's mass through its rest position moving
691+ // right is a point on the Poincare section (for N = 3, drawn only while
692+ // the chain is in the symmetric family).
693+ if (sec && c < 0 && s.u[sec.cross] >= 0) {
694+ const f = -c / (s.u[sec.cross] - c);
695+ const inside = P.N !== 3 || (s.u[0] === s.u[2] && s.p[0] === s.p[2]);
696+ poincare?.addLive(x + f * (s.u[sec.plot] - x), q + f * (s.p[sec.plot] - q), energy(P, s), inside);
697+ }
650698 ghost?.verlet.step(ghost.s, H);
651699 lyap.step(s, H, t + H);
652700 carry -= H;
@@ -663,6 +711,8 @@ function frame(now: number) {
663711 for (const u of s.u) umax = Math.max(umax, Math.abs(u));
664712 ampl = Math.max(umax, 0.998 * ampl);
665713 drawChain();
714+ if (running) poincare?.follow(energy(P, s));
715+ poincare?.draw();
666716 showSeparation();
667717 if (!barsFig.hidden) drawModes();
668718 traces.draw();
src/fput-chain/poincare-core.tsadded+136−0View file
@@ -0,0 +1,136 @@
1+// Poincare sections at fixed energy E, without damping, for chains whose
2+// motion (or part of it) has two degrees of freedom:
3+//
4+// N = 2: recorded each time mass 1 passes its rest position moving right
5+// (u1 = 0, p1 > 0), plotted as mass 2's (u2, p2).
6+//
7+// N = 3, symmetric family: motions with u1 = u3 and p1 = p3 stay that way
8+// (the cubic force couples odd modes only to odd modes; only modes 1 and 3
9+// move). Recorded each time mass 2 passes its rest position moving right
10+// (u2 = 0, p2 > 0), plotted as mass 1's (u1, p1) (mass 3's are the same).
11+// The family is unstable to breaking the symmetry: an exactly symmetric start
12+// stays symmetric, but a tiny asymmetry grows to full size (1e-12 within
13+// ~100 s at 0.25 to 25 J).
14+//
15+// On either section the springs are stretched by 0 or +-x, where x is the
16+// plotted displacement, so the potential energy is vFactor v(x) with
17+// v(d) = k d^2 / 2 + k3 d^4 / 4 (vFactor 2 for N = 2: springs 2 and 3;
18+// 4 for N = 3: all four springs). The plotted momentum q counts `mult` times
19+// in the kinetic energy (1 for N = 2; 2 for N = 3, masses 1 and 3), and the
20+// crossing mass's momentum c follows from the energy:
21+// c^2 / 2m = E - vFactor v(x) - mult q^2 / 2m. A point is reachable exactly
22+// when mult q^2 <= 2 m (E - vFactor v(x)) (the outline); c is the positive root.
23+
24+import { Verlet, energy, stepFor, type Params, type State } from "./chain";
25+import { Lyapunov } from "./lyapunov";
26+
27+export const CROSSINGS = 300; // crossings recorded per orbit
28+
29+export interface SectionDef {
30+ cross: number; // index of the mass whose rest crossing (moving right) defines the section
31+ plot: number; // index of the mass whose (u, p) is plotted
32+ vFactor: number;
33+ mult: number;
34+ uLabel: string;
35+ pLabel: string;
36+}
37+
38+// The section for N masses, or null if there is none.
39+export function sectionDef(N: number): SectionDef | null {
40+ if (N === 2) return { cross: 0, plot: 1, vFactor: 2, mult: 1, uLabel: "u2", pLabel: "p2" };
41+ if (N === 3) return { cross: 1, plot: 0, vFactor: 4, mult: 2, uLabel: "u1 = u3", pLabel: "p1 = p3" };
42+ return null;
43+}
44+
45+const spring = (P: Params, d: number) => 0.5 * P.k * d * d + 0.25 * P.k3 * d ** 4;
46+
47+// Largest |x| on the section: vFactor v(x) = E.
48+export function uMax(P: Params, D: SectionDef, E: number): number {
49+ const e = E / D.vFactor; // k x^2 / 2 + k3 x^4 / 4 = e
50+ return Math.sqrt((-P.k + Math.sqrt(P.k * P.k + 4 * P.k3 * e)) / P.k3);
51+}
52+
53+// Largest |q| on the section at x (0 outside the allowed range).
54+export function pMax(P: Params, D: SectionDef, E: number, x: number): number {
55+ return Math.sqrt(Math.max(0, (2 * P.m * (E - D.vFactor * spring(P, x))) / D.mult));
56+}
57+
58+// Full state on the section, or null if (x, q) is not reachable at E.
59+export function sectionState(P: Params, D: SectionDef, E: number, x: number, q: number): State | null {
60+ const c2 = 2 * P.m * (E - D.vFactor * spring(P, x)) - D.mult * q * q;
61+ if (!(c2 >= 0)) return null;
62+ const c = Math.sqrt(c2);
63+ if (P.N === 2) return { u: new Float64Array([0, x]), p: new Float64Array([c, q]), W: 0 };
64+ return { u: new Float64Array([x, 0, x]), p: new Float64Array([q, c, q]), W: 0 };
65+}
66+
67+// Starting points spread over the allowed region: an nx-by-ny grid in
68+// coordinates normalized to the region (across x, then across q at that x).
69+export function seedGrid(P: Params, D: SectionDef, E: number, nx: number, ny: number): [number, number][] {
70+ const um = uMax(P, D, E);
71+ const seeds: [number, number][] = [];
72+ for (let i = 0; i < nx; i++) {
73+ const x = um * (-1 + (2 * (i + 0.5)) / nx);
74+ const pm = pMax(P, D, E, x);
75+ for (let j = 0; j < ny; j++) seeds.push([x, pm * (-1 + (2 * (j + 0.5)) / ny)]);
76+ }
77+ return seeds;
78+}
79+
80+// One orbit followed from a point on the section, without damping, recording
81+// its crossings and measuring its Lyapunov exponent over the whole run. For
82+// N = 3 the orbit stays exactly symmetric (the shadow's offset is symmetric
83+// too), so its exponent is the one within the symmetric family.
84+export class SectionOrbit {
85+ private P: Params;
86+ private D: SectionDef;
87+ private s: State;
88+ private verlet: Verlet;
89+ readonly lyap = new Lyapunov(Infinity);
90+ points: number[]; // x, q pairs (the start point first)
91+ crossings = 0;
92+ private t = 0;
93+ private h: number;
94+ private maxTime: number;
95+
96+ constructor(P: Params, E: number, x: number, q: number, maxTime: number) {
97+ this.P = { ...P, tau: 0 };
98+ this.D = sectionDef(P.N)!;
99+ this.s = sectionState(this.P, this.D, E, x, q)!;
100+ this.verlet = new Verlet(this.P, this.s);
101+ this.h = stepFor(this.P, energy(this.P, this.s));
102+ this.lyap.start(this.P, this.s, 0);
103+ this.maxTime = maxTime;
104+ this.points = [x, q];
105+ }
106+
107+ get done(): boolean {
108+ return this.t >= this.maxTime;
109+ }
110+
111+ // Integrate until the next crossing (or the time limit); returns whether a
112+ // crossing was recorded. The crossing is located by linear interpolation in
113+ // the step where the crossing mass's displacement changes sign from
114+ // negative to positive.
115+ advance(): boolean {
116+ const { s, h } = this;
117+ const { cross, plot } = this.D;
118+ while (this.t < this.maxTime) {
119+ const c = s.u[cross], x = s.u[plot], q = s.p[plot];
120+ this.verlet.step(s, h);
121+ this.t += h;
122+ this.lyap.step(s, h, this.t);
123+ if (c < 0 && s.u[cross] >= 0) {
124+ const f = -c / (s.u[cross] - c);
125+ this.points.push(x + f * (s.u[plot] - x), q + f * (s.p[plot] - q));
126+ this.crossings++;
127+ return true;
128+ }
129+ }
130+ return false;
131+ }
132+
133+ lambda(): number {
134+ return this.lyap.mean();
135+ }
136+}
src/fput-chain/poincare-worker.tsadded+81−0View file
@@ -0,0 +1,81 @@
1+// Computes Poincare section orbits in a Web Worker (see poincare-core.ts).
2+// Orbits are advanced round-robin, one crossing at a time, and new points are
3+// posted every BATCH_MS so the whole picture fills in together. A new "run"
4+// (energy or parameters changed) replaces the current one; "add" appends an
5+// orbit started at a clicked point.
6+
7+import type { Params } from "./chain";
8+import { CROSSINGS, SectionOrbit } from "./poincare-core";
9+
10+const BATCH_MS = 40;
11+const MAX_TIME = 5000; // s; an orbit that stops crossing ends here
12+
13+export type ToPoincare =
14+ | { type: "run"; id: number; P: Params; E: number; seeds: [number, number][]; crossings?: number }
15+ | { type: "add"; id: number; seed: [number, number] }
16+ | { type: "stop" };
17+
18+export type FromPoincare = {
19+ type: "points";
20+ id: number;
21+ orbits: { index: number; points: Float32Array; lambda: number; done: boolean }[];
22+ finished: boolean;
23+};
24+
25+let runId = -1;
26+let P: Params;
27+let E = 0;
28+let orbits: SectionOrbit[] = [];
29+let sent: number[] = []; // points (pairs) already posted, per orbit
30+let doneSent: boolean[] = []; // whether the orbit's end has been posted
31+let running = false;
32+let limit = CROSSINGS; // crossings per orbit in this run
33+
34+onmessage = (e: MessageEvent<ToPoincare>) => {
35+ const msg = e.data;
36+ if (msg.type === "run") {
37+ runId = msg.id;
38+ limit = msg.crossings ?? CROSSINGS;
39+ P = msg.P;
40+ E = msg.E;
41+ orbits = msg.seeds.map(([a, b]) => new SectionOrbit(P, E, a, b, MAX_TIME));
42+ sent = orbits.map(() => 0);
43+ doneSent = orbits.map(() => false);
44+ } else if (msg.type === "add") {
45+ if (msg.id !== runId) return;
46+ orbits.push(new SectionOrbit(P, E, msg.seed[0], msg.seed[1], MAX_TIME));
47+ sent.push(0);
48+ doneSent.push(false);
49+ } else {
50+ orbits = [];
51+ sent = [];
52+ doneSent = [];
53+ }
54+ if (!running && orbits.length) {
55+ running = true;
56+ setTimeout(work, 0);
57+ }
58+};
59+
60+const active = (o: SectionOrbit) => o.crossings < limit && !o.done;
61+
62+function work(): void {
63+ const t0 = performance.now();
64+ while (performance.now() - t0 < BATCH_MS && orbits.some(active)) {
65+ for (const o of orbits) if (active(o)) o.advance();
66+ }
67+ const out: FromPoincare["orbits"] = [];
68+ orbits.forEach((o, index) => {
69+ const n = o.points.length / 2;
70+ const done = !active(o);
71+ if (n > sent[index] || (done && !doneSent[index])) {
72+ out.push({ index, points: Float32Array.from(o.points.slice(2 * sent[index])), lambda: o.lambda(), done });
73+ sent[index] = n;
74+ doneSent[index] = done;
75+ }
76+ });
77+ const finished = !orbits.some(active);
78+ postMessage({ type: "points", id: runId, orbits: out, finished } satisfies FromPoincare);
79+ if (finished) running = false;
80+ else setTimeout(work, 0);
81+}
src/fput-chain/poincare.tsadded+529−0View file
@@ -0,0 +1,529 @@
1+// Poincare section panel, for N = 2 and for N = 3's symmetric family only
2+// (see poincare-core.ts), adapted from the double pendulum app: its own
3+// energy slider (same parameters as the
4+// chain above, no damping), the allowed region's outline, orbits seeded on a
5+// grid and computed in a worker, click to add an orbit (or, with "Click to
6+// strike" on, also strike the chain above from that exact state), hover to
7+// highlight one, and coloring per orbit or by Lyapunov exponent. The animated
8+// chain's own crossings are drawn as dark rings ("you are here"), only those
9+// within LIVE_TOLERANCE of the section's energy (others would fall outside
10+// the allowed region), and for N = 3 only while the chain is in the
11+// symmetric family. Follow "real time" makes the section's energy track the
12+// chain's as it decays. Follow "sped up" shows the playing note's crossings
13+// instead, as a short fading trail, with the axes fitted to the note's
14+// current energy and, if an atlas has been computed ("Precompute sections",
15+// shown only with Follow "sped up"), the atlas section nearest that energy
16+// behind it (see showSound). Setting the section's energy by hand pauses
17+// everything (onUserEnergy), so Follow does not override it until Resume.
18+
19+import type { Params, State } from "./chain";
20+import { Atlas, ATLAS_SIZE, type AtlasEntry } from "./atlas";
21+import { pMax, sectionDef, sectionState, seedGrid, uMax, type SectionDef } from "./poincare-core";
22+import type { FromPoincare, ToPoincare } from "./poincare-worker";
23+
24+const WIDTH = 960, HEIGHT = 480; // CSS pixels
25+const TOP = 18; // legend row
26+const PAD = 10;
27+const FONT = "11px ui-monospace, SFMono-Regular, Menlo, monospace";
28+const GRID = [6, 7]; // seed grid: across x, across q
29+const LIVE_MAX = 400; // chain crossings kept
30+const LIVE_TOLERANCE = 0.05; // relative energy difference for drawing a chain crossing
31+const FOLLOW_CHANGE = 0.03; // relative energy change that triggers a recompute
32+const FOLLOW_MIN_MS = 2000; // and at most this often
33+const NOTE_TRAIL = 150; // note crossings shown (each fitted from its own energy)
34+const RESCALE_STEP = 0.02; // sound view: refit after this relative energy change
35+// Lambda coloring is logarithmic between these (1/s): regular orbits measure
36+// ~0.05 or less over a run (lambda -> 0 as the run grows), chaotic ones in the
37+// larger chains ~0.5 to 1.5.
38+const LAMBDA_LO = 0.01, LAMBDA_HI = 1.5;
39+
40+interface Orbit {
41+ points: number[]; // x, q pairs (the plotted mass's displacement and momentum)
42+ lambda: number;
43+ done: boolean;
44+ clicked: boolean;
45+}
46+
47+export class PoincarePanel {
48+ private P: Params;
49+ private D: SectionDef | null = null;
50+ private E: number;
51+ private panel: HTMLDetailsElement;
52+ private canvas: HTMLCanvasElement;
53+ private dpr = window.devicePixelRatio || 1;
54+ private worker: Worker | null = null;
55+ private runId = 0;
56+ private orbits: Orbit[] = [];
57+ private finished = false;
58+ private changed = true;
59+ private hoverOrbit = -1;
60+ private hoverPoint: [number, number] | null = null;
61+ private byLambda = false;
62+ private energyInput: HTMLInputElement;
63+ private energyValue: HTMLElement;
64+ private debounce = 0;
65+ private strikeBox: HTMLInputElement;
66+ private live: number[] = []; // the chain's crossings: x, q, E triples
67+ private outside = false; // N = 3: the chain has left the symmetric family
68+ private followSelect: HTMLSelectElement;
69+ private lastFollow = 0;
70+ private noteDots: number[] = []; // x, q, E triples of the playing note
71+ private atlas = new Atlas();
72+ private atlasEntries: AtlasEntry[] = [];
73+ private atlasComputing = false;
74+ private atlasButton: HTMLButtonElement;
75+ private soundView = false; // showing the sound view (Follow "sped up" during a note)
76+ private soundE = 0; // energy the sound view is fitted to
77+ private onStrike: (s: State, E: number) => void;
78+ private onUserEnergy: () => void;
79+ private fmtEnergy: (E: number) => string;
80+ // Axes (data ranges), set from E.
81+ private uM = 1;
82+ private pM = 1;
83+
84+ // onUserEnergy: the section's energy was set by hand (the page pauses).
85+ constructor(
86+ P: Params,
87+ onStrike: (s: State, E: number) => void,
88+ onUserEnergy: () => void,
89+ fmtEnergy: (E: number) => string,
90+ ) {
91+ this.P = P;
92+ this.onStrike = onStrike;
93+ this.onUserEnergy = onUserEnergy;
94+ this.fmtEnergy = fmtEnergy;
95+ this.followSelect = document.getElementById("pc-follow") as HTMLSelectElement;
96+ this.followSelect.addEventListener("change", () => {
97+ // Leaving the sound view: back to the section at the slider's energy.
98+ if (this.followMode() !== "sound" && this.soundView) this.run();
99+ this.showAtlasButton();
100+ this.changed = true;
101+ });
102+ this.atlasButton = document.getElementById("pc-atlas") as HTMLButtonElement;
103+ this.atlasButton.addEventListener("click", () => this.computeAtlas());
104+ this.panel = document.getElementById("poincare") as HTMLDetailsElement;
105+ this.strikeBox = document.getElementById("pc-strike") as HTMLInputElement;
106+ this.energyInput = document.getElementById("pc-energy") as HTMLInputElement;
107+ this.energyValue = document.getElementById("pc-energy-value")!;
108+ this.E = 10 ** Number(this.energyInput.value);
109+ this.showEnergy();
110+ this.canvas = document.createElement("canvas");
111+ this.canvas.width = WIDTH * this.dpr;
112+ this.canvas.height = HEIGHT * this.dpr;
113+ this.canvas.style.cursor = "crosshair";
114+ document.getElementById("pc-plot")!.appendChild(this.canvas);
115+
116+ this.panel.addEventListener("toggle", () => this.update());
117+ this.energyInput.addEventListener("input", () => {
118+ this.onUserEnergy(); // pauses, so Follow does not override this energy
119+ // The recent crossings (chain rings, note trail) belong to other energies.
120+ this.live = [];
121+ this.noteDots = [];
122+ this.E = 10 ** Number(this.energyInput.value);
123+ this.showEnergy();
124+ // Recompute shortly after the slider stops moving.
125+ clearTimeout(this.debounce);
126+ this.debounce = window.setTimeout(() => this.run(), 150);
127+ });
128+ const lambdaBox = document.getElementById("pc-lambda") as HTMLInputElement;
129+ lambdaBox.addEventListener("change", () => {
130+ this.byLambda = lambdaBox.checked;
131+ this.changed = true;
132+ });
133+ this.canvas.addEventListener("mousemove", (e) => this.onHover(e));
134+ this.canvas.addEventListener("mouseleave", () => {
135+ this.hoverOrbit = -1;
136+ this.hoverPoint = null;
137+ this.changed = true;
138+ });
139+ this.canvas.addEventListener("click", (e) => this.onClick(e));
140+ this.setParams(P);
141+ }
142+
143+ // Offered only where there is a section; computed only while shown and open.
144+ private get active(): boolean {
145+ return this.D !== null && this.panel.open;
146+ }
147+
148+ private update(): void {
149+ this.panel.hidden = this.D === null;
150+ if (this.active) this.run();
151+ else this.send({ type: "stop" });
152+ }
153+
154+ // New chain parameters (the panel shows only for N = 2 and 3).
155+ setParams(P: Params): void {
156+ this.P = P;
157+ this.D = sectionDef(P.N);
158+ this.live = [];
159+ this.noteDots = [];
160+ this.outside = false;
161+ this.atlas.cancel();
162+ this.atlasComputing = false;
163+ this.setAtlas([]);
164+ if (this.D) {
165+ this.atlas.load(P).then((entries) => {
166+ if (this.P === P) this.setAtlas(entries);
167+ });
168+ }
169+ this.update();
170+ }
171+
172+ followMode(): string {
173+ return this.followSelect.value;
174+ }
175+
176+ // A note started at energy E: with Follow "sped up", switch to the sound view.
177+ noteStarted(E: number): void {
178+ this.noteDots = [];
179+ this.changed = true;
180+ if (this.D && this.followMode() === "sound") this.showSound(E, true);
181+ }
182+
183+ // The playing note's current energy (Follow "sped up").
184+ setNoteEnergy(E: number): void {
185+ if (this.D && this.followMode() === "sound" && E > 0) this.showSound(E);
186+ }
187+
188+ // A crossing of the playing note, as the audio clock reaches it, with the
189+ // note's energy E then (the dot is fitted to the current energy's region
190+ // when drawn, like the section).
191+ addNoteDot(x: number, q: number, E: number): void {
192+ this.noteDots.push(x, q, E);
193+ if (this.noteDots.length > 3 * NOTE_TRAIL) this.noteDots.splice(0, this.noteDots.length - 3 * NOTE_TRAIL);
194+ if (this.followMode() === "sound") this.changed = true;
195+ }
196+
197+ // The chain's current energy: with Follow "real time", recompute the
198+ // section at it once it has changed by FOLLOW_CHANGE (at most every
199+ // FOLLOW_MIN_MS), or at once when `now` is set (a new strike).
200+ follow(E: number, now = false): void {
201+ if (!this.D || this.followMode() !== "realtime" || !(E > 0) || Math.abs(E - this.E) < 1e-9 * E) return;
202+ const t = performance.now();
203+ if (!now && (Math.abs(E - this.E) < FOLLOW_CHANGE * this.E || t - this.lastFollow < FOLLOW_MIN_MS)) return;
204+ this.lastFollow = t;
205+ this.E = E;
206+ this.energyInput.value = String(Math.log10(E));
207+ this.showEnergy();
208+ this.run();
209+ }
210+
211+ private setAtlas(entries: AtlasEntry[]): void {
212+ this.atlasEntries = entries;
213+ this.showAtlasButton();
214+ if (this.soundView) this.showSound(this.soundE, true);
215+ }
216+
217+ // The atlas only matters for Follow "sped up": the button is shown there
218+ // until the atlas is complete.
219+ private showAtlasButton(): void {
220+ const n = this.atlasEntries.length;
221+ this.atlasButton.hidden = this.followMode() !== "sound" || n === ATLAS_SIZE;
222+ this.atlasButton.disabled = this.atlasComputing;
223+ this.atlasButton.textContent = this.atlasComputing
224+ ? `Computing sections ${n}/${ATLAS_SIZE}`
225+ : n ? `Precompute sections (${n}/${ATLAS_SIZE})` : "Precompute sections";
226+ }
227+
228+ private computeAtlas(): void {
229+ this.atlasComputing = true;
230+ this.showAtlasButton();
231+ const P = this.P;
232+ this.atlas
233+ .compute(P, (entries) => {
234+ if (this.P === P) this.setAtlas(entries);
235+ })
236+ .finally(() => {
237+ if (this.P !== P) return;
238+ this.atlasComputing = false;
239+ this.showAtlasButton();
240+ });
241+ }
242+
243+ // Extent of the allowed region at energy E: largest |x| and |q|.
244+ private extent(E: number): [number, number] {
245+ return [uMax(this.P, this.D!, E), pMax(this.P, this.D!, E, 0)];
246+ }
247+
248+ // Sound view at the note's energy E (no computation): axes fitted to E's
249+ // allowed region, with the atlas section nearest E (if any) fitted to it,
250+ // each point (x, q) of the section at E_k drawn at
251+ // (x xmax(E)/xmax(E_k), q qmax(E)/qmax(E_k)). In the linear regime (low
252+ // energy) sections at different energies are scaled copies, so this is
253+ // exact; above it, it is an approximation over the gap between atlas
254+ // energies (a factor ~1.45). Without an atlas only the outline is drawn.
255+ // Refitted after a RESCALE_STEP change.
256+ private showSound(E: number, force = false): void {
257+ if (!(E > 0) || !this.panel.open) return;
258+ if (!force && this.soundView && Math.abs(Math.log(E / this.soundE)) < RESCALE_STEP) return;
259+ if (!this.soundView) {
260+ this.runId++; // ignore the interactive computation's remaining results
261+ this.send({ type: "stop" });
262+ this.soundView = true;
263+ }
264+ this.soundE = E;
265+ const [um, pm] = this.extent(E);
266+ this.orbits = [];
267+ if (this.atlasEntries.length) {
268+ let k = 0;
269+ for (let i = 1; i < this.atlasEntries.length; i++) {
270+ if (Math.abs(Math.log(this.atlasEntries[i].E / E)) < Math.abs(Math.log(this.atlasEntries[k].E / E))) k = i;
271+ }
272+ const entry = this.atlasEntries[k];
273+ const [uk, pk] = this.extent(entry.E);
274+ const fu = um / uk, fp = pm / pk;
275+ this.orbits = entry.orbits.map((o) => {
276+ const points = new Array<number>(o.points.length);
277+ for (let i = 0; i < o.points.length; i += 2) {
278+ points[i] = o.points[i] * fu;
279+ points[i + 1] = o.points[i + 1] * fp;
280+ }
281+ return { points, lambda: o.lambda, done: true, clicked: false };
282+ });
283+ }
284+ this.E = E;
285+ this.energyInput.value = String(Math.log10(E));
286+ this.showEnergy();
287+ this.uM = um;
288+ this.pM = pm;
289+ this.finished = true;
290+ this.changed = true;
291+ }
292+
293+ // A crossing of the animated chain at energy E (for N = 3, inside: whether
294+ // the chain is in the symmetric family; crossings outside it are not drawn).
295+ addLive(x: number, q: number, E: number, inside = true): void {
296+ this.outside = !inside;
297+ this.changed = true;
298+ if (!inside) return;
299+ this.live.push(x, q, E);
300+ if (this.live.length > 3 * LIVE_MAX) this.live.splice(0, this.live.length - 3 * LIVE_MAX);
301+ }
302+
303+ clearLive(): void {
304+ this.live = [];
305+ this.outside = false;
306+ this.changed = true;
307+ }
308+
309+ private showEnergy(): void {
310+ this.energyValue.textContent = this.fmtEnergy(this.E);
311+ }
312+
313+ private send(msg: ToPoincare): void {
314+ if (!this.worker) {
315+ if (msg.type === "stop") return;
316+ this.worker = new Worker(new URL("./poincare-worker.ts", import.meta.url), { type: "module" });
317+ this.worker.onmessage = (e: MessageEvent<FromPoincare>) => this.receive(e.data);
318+ }
319+ this.worker.postMessage(msg);
320+ }
321+
322+ // Start computing the section at the current energy.
323+ private run(): void {
324+ if (!this.active) return;
325+ this.soundView = false;
326+ this.soundE = 0;
327+ this.runId++;
328+ this.orbits = [];
329+ this.finished = false;
330+ const D = this.D!;
331+ this.uM = uMax(this.P, D, this.E);
332+ this.pM = pMax(this.P, D, this.E, 0);
333+ const seeds = seedGrid(this.P, D, this.E, GRID[0], GRID[1]);
334+ for (let i = 0; i < seeds.length; i++) this.orbits.push({ points: [], lambda: NaN, done: false, clicked: false });
335+ this.send({ type: "run", id: this.runId, P: this.P, E: this.E, seeds });
336+ this.changed = true;
337+ }
338+
339+ private receive(msg: FromPoincare): void {
340+ if (msg.id !== this.runId) return; // from an earlier energy
341+ for (const o of msg.orbits) {
342+ const orbit = this.orbits[o.index];
343+ if (!orbit) continue;
344+ for (const v of o.points) orbit.points.push(v);
345+ orbit.lambda = o.lambda;
346+ orbit.done = o.done;
347+ }
348+ this.finished = msg.finished;
349+ this.changed = true;
350+ }
351+
352+ // Data <-> CSS pixel coordinates.
353+ private px(u: number): number {
354+ return PAD + ((u + this.uM * 1.04) / (2 * this.uM * 1.04)) * (WIDTH - 2 * PAD);
355+ }
356+ private py(p: number): number {
357+ return TOP + PAD + ((this.pM * 1.06 - p) / (2 * this.pM * 1.06)) * (HEIGHT - TOP - 2 * PAD);
358+ }
359+ private fromPx(x: number, y: number): [number, number] {
360+ const u = ((x - PAD) / (WIDTH - 2 * PAD)) * 2 * this.uM * 1.04 - this.uM * 1.04;
361+ const p = this.pM * 1.06 - ((y - TOP - PAD) / (HEIGHT - TOP - 2 * PAD)) * 2 * this.pM * 1.06;
362+ return [u, p];
363+ }
364+ private mouse(e: MouseEvent): [number, number] {
365+ const r = this.canvas.getBoundingClientRect();
366+ return [((e.clientX - r.left) / r.width) * WIDTH, ((e.clientY - r.top) / r.height) * HEIGHT];
367+ }
368+
369+ private onHover(e: MouseEvent): void {
370+ const [x, y] = this.mouse(e);
371+ this.hoverPoint = this.fromPx(x, y);
372+ let best = 8 * 8; // within 8 px
373+ this.hoverOrbit = -1;
374+ this.orbits.forEach((o, k) => {
375+ for (let i = 0; i < o.points.length; i += 2) {
376+ const dx = this.px(o.points[i]) - x, dy = this.py(o.points[i + 1]) - y;
377+ if (dx * dx + dy * dy < best) {
378+ best = dx * dx + dy * dy;
379+ this.hoverOrbit = k;
380+ }
381+ }
382+ });
383+ this.changed = true;
384+ }
385+
386+ // Add an orbit starting at the clicked point, if it is reachable, and with
387+ // "Click to strike" on, strike the chain from that state.
388+ private onClick(e: MouseEvent): void {
389+ const [x, y] = this.mouse(e);
390+ const [u, p] = this.fromPx(x, y);
391+ const state = sectionState(this.P, this.D!, this.E, u, p);
392+ if (!state) return;
393+ if (this.strikeBox.checked) this.onStrike(state, this.E);
394+ if (this.soundView) return; // the sound view is not computed here
395+ this.orbits.push({ points: [], lambda: NaN, done: false, clicked: true });
396+ this.finished = false;
397+ this.send({ type: "add", id: this.runId, seed: [u, p] });
398+ this.changed = true;
399+ }
400+
401+ private color(k: number, alpha: number): string {
402+ const o = this.orbits[k];
403+ if (this.byLambda) {
404+ if (!Number.isFinite(o.lambda)) return `rgba(150, 150, 150, ${alpha})`;
405+ // Blue (regular) -> amber -> red (chaotic).
406+ // (A regular orbit's lambda can come out slightly negative.)
407+ const u = o.lambda <= LAMBDA_LO ? 0 : Math.min(1, Math.log(o.lambda / LAMBDA_LO) / Math.log(LAMBDA_HI / LAMBDA_LO));
408+ const ramp = [[47, 93, 138], [227, 162, 26], [194, 65, 45]];
409+ const s = u * 2, i = Math.min(1, Math.floor(s)), w = s - i;
410+ const c = ramp[i].map((a, j) => Math.round(a + (ramp[i + 1][j] - a) * w));
411+ return `rgba(${c[0]}, ${c[1]}, ${c[2]}, ${alpha})`;
412+ }
413+ const hue = (k * 137.508) % 360; // golden angle: neighbouring orbits differ
414+ return `hsla(${hue}, 45%, 42%, ${alpha})`;
415+ }
416+
417+ draw(): void {
418+ if (!this.active || !this.changed) return;
419+ this.changed = false;
420+ const D = this.D!;
421+ const ctx = this.canvas.getContext("2d")!;
422+ ctx.setTransform(this.dpr, 0, 0, this.dpr, 0, 0);
423+ ctx.clearRect(0, 0, WIDTH, HEIGHT);
424+
425+ // Axes through 0 and the allowed region's outline.
426+ ctx.strokeStyle = "#ecece6";
427+ ctx.lineWidth = 1;
428+ ctx.beginPath();
429+ ctx.moveTo(this.px(0), TOP);
430+ ctx.lineTo(this.px(0), HEIGHT);
431+ ctx.moveTo(0, this.py(0));
432+ ctx.lineTo(WIDTH, this.py(0));
433+ ctx.stroke();
434+ ctx.strokeStyle = "#b9b8b0";
435+ ctx.beginPath();
436+ for (const sign of [1, -1]) {
437+ for (let i = 0; i <= 200; i++) {
438+ const u = -this.uM + (2 * this.uM * i) / 200;
439+ const x = this.px(u), y = this.py(sign * pMax(this.P, D, this.E, u));
440+ if (i === 0) ctx.moveTo(x, y);
441+ else ctx.lineTo(x, y);
442+ }
443+ }
444+ ctx.stroke();
445+
446+ // Orbits; the hovered one on top and the others faded while hovering.
447+ const order = this.orbits.map((_, k) => k).filter((k) => k !== this.hoverOrbit);
448+ if (this.hoverOrbit >= 0) order.push(this.hoverOrbit);
449+ for (const k of order) {
450+ const o = this.orbits[k];
451+ const hovered = k === this.hoverOrbit;
452+ ctx.fillStyle = this.color(k, this.hoverOrbit >= 0 && !hovered ? 0.25 : 0.85);
453+ const size = hovered || o.clicked ? 2.5 : 1.8;
454+ for (let i = 0; i < o.points.length; i += 2) {
455+ ctx.fillRect(this.px(o.points[i]) - size / 2, this.py(o.points[i + 1]) - size / 2, size, size);
456+ }
457+ }
458+
459+ if (this.followMode() === "sound") {
460+ // The playing note's latest crossings: a fading trail, the newest
461+ // larger, each fitted from its own energy to the current one.
462+ const n = this.noteDots.length / 3;
463+ for (let k = 0; k < n; k++) {
464+ const latest = k === n - 1;
465+ const [uk, pk] = this.extent(this.noteDots[3 * k + 2]);
466+ ctx.globalAlpha = 0.2 + (0.8 * (k + 1)) / n;
467+ ctx.beginPath();
468+ ctx.arc(
469+ this.px((this.noteDots[3 * k] * this.uM) / uk),
470+ this.py((this.noteDots[3 * k + 1] * this.pM) / pk),
471+ latest ? 6 : 3.5, 0, 2 * Math.PI,
472+ );
473+ ctx.fillStyle = "#fff";
474+ ctx.fill();
475+ ctx.strokeStyle = "#c2412d";
476+ ctx.lineWidth = latest ? 2.5 : 1.5;
477+ ctx.stroke();
478+ }
479+ ctx.globalAlpha = 1;
480+ } else {
481+ // The chain's own crossings at about this energy: dark rings on white,
482+ // the latest larger.
483+ ctx.strokeStyle = "#23262b";
484+ ctx.fillStyle = "#fff";
485+ for (let i = 0; i < this.live.length; i += 3) {
486+ if (Math.abs(this.live[i + 2] - this.E) > LIVE_TOLERANCE * this.E) continue;
487+ const latest = i === this.live.length - 3;
488+ ctx.beginPath();
489+ ctx.arc(this.px(this.live[i]), this.py(this.live[i + 1]), latest ? 6 : 3.5, 0, 2 * Math.PI);
490+ ctx.fill();
491+ ctx.lineWidth = latest ? 2.5 : 1.5;
492+ ctx.stroke();
493+ }
494+ }
495+
496+ // Legend row.
497+ ctx.font = FONT;
498+ ctx.textBaseline = "top";
499+ ctx.fillStyle = "#2f5d8a";
500+ let legend: string;
501+ if (this.hoverOrbit >= 0) {
502+ const o = this.orbits[this.hoverOrbit];
503+ const kind = !Number.isFinite(o.lambda) ? "" : o.lambda > 0.1 ? "chaotic" : "regular";
504+ legend = `orbit ${this.hoverOrbit + 1} · λ ${Number.isFinite(o.lambda) ? o.lambda.toFixed(3) : "–"} 1/s ${kind ? `(${kind})` : ""}`;
505+ } else if (this.hoverPoint) {
506+ const [u, p] = this.hoverPoint;
507+ legend = `${D.uLabel} ${(100 * u).toFixed(2)} cm · ${D.pLabel} ${p.toFixed(3)} kg m/s${sectionState(this.P, D, this.E, u, p) ? " · click to add an orbit here" : ""}`;
508+ } else {
509+ legend = `${D.uLabel} (cm) across · ${D.pLabel} (kg m/s) up · ${this.orbits.length} orbits`;
510+ if (this.P.N === 3) legend += " · symmetric motions only";
511+ // Say why the chain's own crossings are missing.
512+ const lastE = this.live[this.live.length - 1];
513+ if (this.followMode() === "sound") {
514+ if (this.soundView && !this.atlasEntries.length) legend += " · no precomputed sections";
515+ } else if (this.outside) legend += " · chain not symmetric, not shown";
516+ else if (lastE > 0 && Math.abs(lastE - this.E) > LIVE_TOLERANCE * this.E) {
517+ legend += ` · chain at ${this.fmtEnergy(lastE)}, not shown`;
518+ }
519+ }
520+ ctx.fillText(legend, 4, 3);
521+ ctx.fillStyle = "#999";
522+ ctx.textAlign = "right";
523+ ctx.fillText(
524+ `${this.finished || this.soundView ? "" : "computing… · "}±${(100 * this.uM).toFixed(2)} cm · ±${this.pM.toFixed(3)} kg m/s`,
525+ WIDTH - 4, 3,
526+ );
527+ ctx.textAlign = "left";
528+ }
529+}
src/fput-chain/render-worker.tsmodified+17−2View file
@@ -10,6 +10,7 @@
1010
1111 import { Verlet, copyState, energy, stepFor, type Params, type State } from "./chain";
1212 import { Lyapunov } from "./lyapunov";
13+import { sectionDef } from "./poincare-core";
1314 import { signalByKey } from "./signals";
1415
1516 export const END_RATIO = 1e-6;
@@ -33,6 +34,9 @@ export interface RenderResult {
3334 samples: Float32Array; // raw signal, in the signal's units
3435 energy: Float32Array; // chain energy (J) at samples 0, DIAG_EVERY, 2 DIAG_EVERY, ...
3536 lambda: Float32Array; // largest Lyapunov exponent over the last LYAP_WINDOW s (1/s), same samples
37+ // Poincare section crossings (N = 2 and 3; for N = 3 only while exactly
38+ // symmetric): sample index, x, q; repeated.
39+ crossings: Float32Array;
3640 }
3741
3842 self.onmessage = (e: MessageEvent<RenderRequest>) => {
@@ -50,6 +54,8 @@ self.onmessage = (e: MessageEvent<RenderRequest>) => {
5054 const lyap = new Lyapunov(LYAP_WINDOW);
5155 let t = 0; // physical time (s)
5256 lyap.start(r.P, s, t);
57+ const sec = sectionDef(r.P.N);
58+ const crossings: number[] = [];
5359 let n = 0;
5460 for (; n < nMax; n++) {
5561 if (n % DIAG_EVERY === 0) {
@@ -62,13 +68,22 @@ self.onmessage = (e: MessageEvent<RenderRequest>) => {
6268 }
6369 out[n] = signal.value(r.P, s, verlet.force, r.j);
6470 for (let k = 0; k < sub; k++) {
71+ const c = sec ? s.u[sec.cross] : 0, x = sec ? s.u[sec.plot] : 0, q = sec ? s.p[sec.plot] : 0;
6572 verlet.step(s, h);
6673 t += h;
6774 lyap.step(s, h, t);
75+ 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]))) {
76+ const f = -c / (s.u[sec.cross] - c);
77+ crossings.push(n, x + f * (s.u[sec.plot] - x), q + f * (s.p[sec.plot] - q));
78+ }
6879 }
6980 }
7081 out = out.slice(0, n);
7182 const nDiag = Math.ceil(n / DIAG_EVERY);
72- const result: RenderResult = { id: r.id, samples: out, energy: en.slice(0, nDiag), lambda: lam.slice(0, nDiag) };
73- (self as unknown as Worker).postMessage(result, [result.samples.buffer, result.energy.buffer, result.lambda.buffer]);
83+ const result: RenderResult = {
84+ id: r.id, samples: out, energy: en.slice(0, nDiag), lambda: lam.slice(0, nDiag), crossings: Float32Array.from(crossings),
85+ };
86+ (self as unknown as Worker).postMessage(result, [
87+ result.samples.buffer, result.energy.buffer, result.lambda.buffer, result.crossings.buffer,
88+ ]);
7489 };