Poincare section panel; minimal text - listening-for-chaos/double-pendulum-chaos

listening-for-chaos / double-pendulum-chaos

Poincare section panel; minimal text

- Poincare section with its own energy slider (same parameters, no damping): allowed-region outline, 42 grid orbits computed in a worker, hover to pick out an orbit, click to add one, color by Lyapunov exponent - Title 'Double pendulum chaos'; explanatory captions and descriptions removed; status shown only for problems or while computing

Jeremy Magland <jeremy.magland@gmail.com> committed 10 days ago commit 5ca3fa433aff parent d3c8c3f Browse files

7 changed files+499−40

README.mdmodified+3−2View file

@@ -2,7 +2,7 @@

22

33 Interactive double pendulum, step 1 of "Listening for chaos". Work in progress.

44

5-Current state: the pendulum hangs at rest until Strike, which starts it from hanging rest with the energy on the slider (default 50 J), plus a random relative jitter of up to 1e-9 in each angle and momentum (so chaotic strikes differ from one another while regular ones sound alike); Pause freezes the animation. The Ghost toggle adds a second pendulum (translucent orange) started from the same state with the outer arm's angle offset by 10^-12 to 10^-1 rad (default 10^-10), restarted on Strike or Reset, with a readout of the angular separation and the time since the ghost started: it grows exponentially when the motion is chaotic and stays tiny when it is regular. Charts "Separation" (ghost separation, log scale, last 30 s) and "Lyapunov" (largest Lyapunov exponent over a 10 s window and since the strike, last 30 s). Parameters (masses, lengths, g, decay time, presets) fold out underneath. Toggle chips choose the plots: a spectrogram of the selected signal (frequency axis auto-scaled, dashed lines at the normal modes) and live traces of energies, energy balance error, angles, momenta and bob positions. The Sound panel at the bottom starts closed; opening it switches the sound engine on (closing it switches it off). While it is open, every Strike plays a note: a separate, sped-up copy of the motion from the same strike; a damped note ends once it has decayed by 60 dB. The mute button (speaker icon) starts muted: notes are still computed and shown, just silent. The panel also has "Sound of" (the signal, also used by the animation spectrogram), volume, pitch (where the small-amplitude lower mode lands, 30 to 1000 Hz, default 100 Hz), and optional whole-note displays that do not scroll (strike at the left edge, end of the note at the right): a spectrogram of the sound, the largest Lyapunov exponent of the motion being heard (10 s window, physical time), its energy (from 0), and a scatter of lambda against energy accumulated over notes (colored by time into the note, earlier notes faded, Clear button). Spectrogram on by default. Hovering over any chart (animation or sound) shows a time cursor and the values there; a click pins it. Charts in the same card share the cursor. The time axis spans the predicted note length ((tau/2) ln 1e6 of physical time, divided by the speedup), is extended if the note runs longer, and is fitted to the note when it ends.

5+Current state: the pendulum hangs at rest until Strike, which starts it from hanging rest with the energy on the slider (default 50 J), plus a random relative jitter of up to 1e-9 in each angle and momentum (so chaotic strikes differ from one another while regular ones sound alike); Pause freezes the animation. The Ghost toggle adds a second pendulum (translucent orange) started from the same state with the outer arm's angle offset by 10^-12 to 10^-1 rad (default 10^-10), restarted on Strike or Reset, with a readout of the angular separation and the time since the ghost started: it grows exponentially when the motion is chaotic and stays tiny when it is regular. Charts "Separation" (ghost separation, log scale, last 30 s) and "Lyapunov" (largest Lyapunov exponent over a 10 s window and since the strike, last 30 s). Parameters (masses, lengths, g, decay time, presets) fold out underneath. Toggle chips choose the plots: a spectrogram of the selected signal (frequency axis auto-scaled, dashed lines at the normal modes) and live traces of energies, energy balance error, angles, momenta and bob positions. The Poincare section panel (closed by default) has its own energy slider and uses the same masses and lengths, without damping: each time the inner arm swings down through vertical moving right, the outer arm's angle and momentum are plotted. 42 orbits start on a grid inside the allowed region (outlined) and leave 300 dots each, computed in a Web Worker and drawn as they arrive; regular orbits trace closed curves, chaotic ones scatter. Hover to pick out an orbit (with its Lyapunov exponent), click inside the outline to start a new orbit there, and optionally color orbits by Lyapunov exponent (log scale, blue regular to red chaotic). The Sound panel at the bottom starts closed; opening it switches the sound engine on (closing it switches it off). While it is open, every Strike plays a note: a separate, sped-up copy of the motion from the same strike; a damped note ends once it has decayed by 60 dB. The mute button (speaker icon) starts muted: notes are still computed and shown, just silent. The panel also has "Sound of" (the signal, also used by the animation spectrogram), volume, pitch (where the small-amplitude lower mode lands, 30 to 1000 Hz, default 100 Hz), and optional whole-note displays that do not scroll (strike at the left edge, end of the note at the right): a spectrogram of the sound, the largest Lyapunov exponent of the motion being heard (10 s window, physical time), its energy (from 0), and a scatter of lambda against energy accumulated over notes (colored by time into the note, earlier notes faded, Clear button). Spectrogram on by default. Hovering over any chart (animation or sound) shows a time cursor and the values there; a click pins it. Charts in the same card share the cursor. The time axis spans the predicted note length ((tau/2) ln 1e6 of physical time, divided by the speedup), is extended if the note runs longer, and is fitted to the note when it ends.

66

77 ```bash

88 npm install

@@ -17,8 +17,9 @@ npm run dev

1717 - Damping: linear drag on each bob, force -gamma m_i v_i (proportional to its mass; equal for equal bobs), so the generalized forces are -gamma p and every mode's amplitude decays with time constant tau = 2/gamma. Default tau = 540 s: about 10% amplitude per minute in the animation, about 2 s in the sound at x262. The state carries the dissipated energy W, so E + W is checked for conservation. With damping the energy slider sets the strike energy; with tau = 0 it sets the energy live.

1818 - Live energy change: both momenta are scaled by one factor at fixed angles (a push along the current motion). A target below the current potential energy can't be reached at once; until it can, kinetic energy is removed at a rate of 6 per second.

1919 - Lyapunov exponent: two-trajectory (Benettin) method in `src/lyapunov.ts`. A hidden shadow pendulum is kept 1e-8 away (angles in rad, momenta in units of sqrt((m1 + m2) l1^2 E_top)), moved back after every step, and the log growth averaged. The shadow keeps its own energy (forcing it onto the pendulum's energy by rescaling momenta inflates the separation near turning points). Checked (no damping, 300 s): 2 J and 10 J fall toward 0 (0.01/s), 30 J ~1.0/s, 50 J ~1.4/s; unchanged for shadow distances 1e-6 and 1e-10. With damping (tau 540 s) a 50 J strike's 10 s value falls from ~1.4 to ~0.1/s by 300 s.

20+- Poincare section (`src/poincare-core.ts`): on the section th1 = 0, a point (th2, p2) is reachable at energy E exactly when p2^2 <= 2 m2 l2^2 (E - V(th2)), V(th2) = m2 g l2 (1 - cos th2); p1 is the root with the inner arm moving right. Crossings are located by linear interpolation within a step. The step is 4x the animation's: over 300 crossings the energy error stays below ~2e-5, the regular/chaotic classification is unchanged, and regular orbits move by ~1e-4. Each orbit's exponent is the shadow method's mean over its run (regular orbits ~0.01/s or less).

2021 - Integrator: 2-stage Gauss-Legendre (order 4, symplectic), fixed step h chosen per parameter set by a short calibration run so the relative energy error stays below 1e-6 up to the top of the energy slider (twice the both-inverted energy); 3.13 ms for the defaults, fixed-point stage solve (started from the previous step's stage derivatives extrapolated in time; tolerance 1e-12, below which the energy error does not change). Angles are wrapped to (-pi, pi] after each step. The physics costs well under 1 ms of CPU per second; drawing dominates, so traces draw at most a min/max pair per pixel column and the spectrogram scrolls in new columns instead of redrawing.

2122

2223 Audio: a Web Worker (`src/audio-worker.ts`) integrates the audio copy with substeps of at most 3.2 ms of simulated time per sample, high-passes at 20 Hz, and streams 1024-sample chunks through a MessageChannel to an AudioWorklet (`src/audio-worklet.ts`) that plays them from a ring buffer, keeping about 0.25 s queued. The page shows how many times faster than real time the worker computes, and any dropouts. The worker also sends each chunk to the page, which feeds it to the audio spectrogram as the audio clock reaches it (2048-sample window, new column every 512 samples). The worker runs the same Lyapunov shadow on the note and sends one exponent and one energy value per 256 samples with each chunk; the page plots them in step with the spectrogram. The whole-note spectrogram keeps every column (pairs are averaged past 4096 columns) and draws only new pixel columns unless the axis or color scale changes.

2324

24-Code: `src/pendulum.ts` (physics), `src/sim.ts` (running simulation with energy target), `src/audio.ts` (audio setup), `src/signals.ts` (audible quantities), `src/lyapunov.ts` (Lyapunov exponent), `src/scatter.ts` (lambda vs energy), `src/cursor.ts` (shared hover/pin cursor), `src/audio-shared.ts` (constants shared with the worker), `src/traces.ts` (strip charts), `src/spectrogram.ts` (spectrogram), `src/main.ts` (page).

25+Code: `src/pendulum.ts` (physics), `src/sim.ts` (running simulation with energy target), `src/audio.ts` (audio setup), `src/signals.ts` (audible quantities), `src/lyapunov.ts` (Lyapunov exponent), `src/scatter.ts` (lambda vs energy), `src/cursor.ts` (shared hover/pin cursor), `src/poincare-core.ts`, `src/poincare-worker.ts`, `src/poincare.ts` (Poincare section), `src/audio-shared.ts` (constants shared with the worker), `src/traces.ts` (strip charts), `src/spectrogram.ts` (spectrogram), `src/main.ts` (page).

index.htmlmodified+26−21View file

@@ -3,7 +3,7 @@

33 <head>

44 <meta charset="UTF-8" />

55 <meta name="viewport" content="width=device-width, initial-scale=1.0" />

6- <title>Double pendulum</title>

6+ <title>Double pendulum chaos</title>

77 <style>

88 :root {

99 --bg: #f6f6f3;

@@ -25,9 +25,8 @@

2525 font: 14px/1.45 var(--sans);

2626 }

2727 main { max-width: 1040px; margin: 0 auto; padding: 28px 16px 48px; }

28- header { margin-bottom: 20px; }

28+ header { margin-bottom: 16px; }

2929 h1 { font: 400 28px/1.2 var(--serif); margin: 0 0 4px; letter-spacing: 0.01em; }

30- .sub { margin: 0; color: var(--muted); }

3130 .grid { display: grid; grid-template-columns: minmax(0, 1fr) minmax(0, 1fr); gap: 16px; align-items: start; }

3231 @media (max-width: 900px) { .grid { grid-template-columns: minmax(0, 1fr); } }

3332 .card { background: var(--card); border: 1px solid var(--line); border-radius: 10px; padding: 16px; }

@@ -74,10 +73,11 @@

7473 summary::-webkit-details-marker { display: none; }

7574 summary::before { content: "▸ "; }

7675 details[open] summary::before { content: "▾ "; }

77- details#audio { border-top: 1px solid var(--line); padding-top: 16px; }

78- details#audio summary { display: flex; align-items: baseline; gap: 10px; }

79- details#audio summary .label { font-size: 12px; color: var(--text); }

80- details#audio[open] summary { margin-bottom: 14px; }

76+ details.panel { border-top: 1px solid var(--line); padding-top: 16px; margin-top: 16px; }

77+ details.panel summary { display: flex; align-items: baseline; gap: 10px; }

78+ details.panel summary .label { font-size: 12px; color: var(--text); }

79+ details.panel[open] summary { margin-bottom: 14px; }

80+ #pc-energy { max-width: 360px; }

8181 #param-form { display: grid; grid-template-columns: repeat(auto-fill, minmax(140px, 1fr)); gap: 8px 16px; margin-top: 10px; }

8282 #param-form label { display: grid; grid-template-columns: 1.6em 5.5em auto; align-items: center; gap: 6px; font-size: 13px; }

8383 #param-form .unit { color: var(--muted); font-size: 12px; }

@@ -98,10 +98,8 @@

9898

9999 figure { margin: 14px 0 0; }

100100 figure canvas { width: 100%; height: auto; border-radius: 4px; }

101- #scatter-fig figcaption { justify-content: space-between; flex-wrap: nowrap; }

102- figcaption { margin-top: 4px; }

101+ .scatter-actions { justify-content: flex-end; margin-top: 6px; }

103102 #traces canvas, #audio-traces canvas { width: 100%; height: auto; margin-top: 12px; }

104- #plots-empty { margin: 14px 0 2px; }

105103

106104 #audio { margin-top: 16px; }

107105 #audio .row + .row { margin-top: 12px; }

@@ -113,8 +111,7 @@

113111 <body>

114112 <main>

115113 <header>

116- <h1>Double pendulum</h1>

117- <p class="sub">Raise the energy and listen for where the motion becomes chaotic.</p>

114+ <h1>Double pendulum chaos</h1>

118115 </header>

119116

120117 <div class="grid">

@@ -168,17 +165,15 @@

168165 <div class="chips" id="toggles">

169166 <label class="chip"><input id="spec-toggle" type="checkbox" /><span>Spectrogram</span></label>

170167 </div>

171- <p id="plots-empty" class="muted">Choose what to plot.</p>

172168 <figure id="spec-fig" hidden>

173169 <div id="spectrogram"></div>

174- <figcaption id="spec-label" class="muted"></figcaption>

175170 </figure>

176171 <div id="traces"></div>

177172 </section>

178173 </div>

179174

180- <details class="card" id="audio">

181- <summary><span class="label">Sound</span><span class="muted">Hear the motion, sped up. Each Strike plays a note.</span></summary>

175+ <details class="card panel" id="audio">

176+ <summary><span class="label">Sound</span></summary>

182177 <div class="row">

183178 <button id="mute" type="button" class="mute" title="Unmute" aria-label="Unmute"></button>

184179 <label><span class="label">Sound of</span> <select id="signal"></select></label>

@@ -194,15 +189,25 @@

194189 <div class="row"><span id="audio-status" class="muted"></span></div>

195190 <figure id="audio-spec-fig">

196191 <div id="audio-spectrogram"></div>

197- <figcaption id="audio-spec-label" class="muted"></figcaption>

198192 </figure>

199193 <div id="audio-traces"></div>

200194 <figure id="scatter-fig" hidden>

201195 <div id="scatter"></div>

202- <figcaption class="muted row">

203- <span>Each point is a moment of a note: its energy and Lyapunov exponent. Color is time into the note on a log scale: red just after the strike, then amber, blue late. Earlier notes are faded; strike several times to build up the picture.</span>

204- <button id="scatter-clear" type="button">Clear</button>

205- </figcaption>

196+ <div class="row scatter-actions"><button id="scatter-clear" type="button">Clear</button></div>

197+ </figure>

198+ </details>

199+

200+ <details class="card panel" id="poincare">

201+ <summary><span class="label">Poincaré section</span></summary>

202+ <div class="row">

203+ <span class="label">Energy</span>

204+ <input id="pc-energy" type="range" value="20" min="0" max="120" step="0.1" aria-label="Section energy" />

205+ <span id="pc-energy-value" class="num">20.0 J</span>

206+ <label class="chip"><input id="pc-lambda" type="checkbox" /><span>Color by λ</span></label>

207+ <span id="pc-status" class="muted"></span>

208+ </div>

209+ <figure>

210+ <div id="pc-plot"></div>

206211 </figure>

207212 </details>

208213 </main>

src/main.tsmodified+9−16View file

@@ -1,6 +1,7 @@

11 import { AudioOut } from "./audio";

22 import { END_RATIO, LAMBDA_EVERY } from "./audio-shared";

33 import { newCursor } from "./cursor";

4+import { PoincarePanel } from "./poincare";

45 import { Lyapunov } from "./lyapunov";

56 import { ScatterPlot } from "./scatter";

67 import {

@@ -72,10 +73,6 @@ const makeSpectrogram = () => {

7273 };

7374 let spectrogram = makeSpectrogram();

7475 specToggle.addEventListener("change", () => (specFig.hidden = !specToggle.checked));

75-const specLabel = document.getElementById("spec-label")!;

76-const showSpecLabel = () =>

77- (specLabel.textContent = `Spectrogram of ${signal.label}, simulated time. Dashed: normal modes.`);

78-showSpecLabel();

7976 let stepCount = 0;

8077

8178 // Energy balance: E plus the energy dissipated since the last strike or

@@ -158,9 +155,6 @@ const traces = new TracePanel(

158155 ],

159156 10, // visible time span (s)

160157 );

161-// Hint when no plot is selected.

162-const plotsEmpty = document.getElementById("plots-empty")!;

163-toggles.addEventListener("change", () => (plotsEmpty.hidden = specToggle.checked || traces.anyShown()));

164158

165159 // Sound: a separate, sped-up simulation. While the Sound panel is open, each

166160 // Strike starts a note from the same state as the animation (opening it while

@@ -178,15 +172,12 @@ const PITCH_MIN = 30, PITCH_MAX = 1000;

178172 const pitchHz = () => PITCH_MIN * (PITCH_MAX / PITCH_MIN) ** (Number(pitchInput.value) / 1000);

179173 pitchInput.value = String(Math.round((1000 * Math.log(100 / PITCH_MIN)) / Math.log(PITCH_MAX / PITCH_MIN)));

180174 const audioModeHz = () => modeHz.map((f) => (f * pitchHz()) / modeHz[0]);

181-const audioSpecLabel = document.getElementById("audio-spec-label")!;

182175 const showPitch = () => {

183176 const [lo, hi] = audioModeHz();

184177 const S = lo / modeHz[0];

185178 pitchValue.textContent =

186179 `${lo.toFixed(0)} Hz · upper ${hi.toFixed(0)} Hz · ×${S.toFixed(0)}` +

187180 (damped() ? ` · decay ${(P.tau / S).toFixed(1)} s` : "");

188- audioSpecLabel.textContent =

189- `Spectrogram of the whole note, strike (left) to end (right): ${signal.label}, sped up ×${(lo / modeHz[0]).toFixed(0)}. Dashed: normal modes.`;

190181 };

191182 showPitch();

192183

@@ -291,17 +282,16 @@ function feedAudioSpec() {

291282 const audio = new AudioOut(P, ALPHA, {

292283 onStats: (s) => {

293284 if (!audio.playing) return;

285+ // Only problems are reported.

294286 audioStatus.textContent =

295- `Computed ${s.speed.toFixed(1)}× faster than real time` +

296- (s.speed < 1.5 ? " (too slow for live audio)" : "") +

297- (s.underruns > 0 ? ` · ${s.underruns} dropouts` : "");

287+ (s.speed < 1.5 ? `Too slow for live audio (computed ${s.speed.toFixed(1)}× real time)` : "") +

288+ (s.underruns > 0 ? `${s.speed < 1.5 ? " · " : ""}${s.underruns} dropouts` : "");

298289 },

299290 onSamples: (data, lambda, energy) => {

300291 queue.push({ data, lambda, energy });

301292 received += data.length;

302293 },

303294 onEnded: () => {

304- audioStatus.textContent = "Note finished (decayed by 60 dB).";

305295 setNoteSpan(received / audio.sampleRate);

306296 scatter.setNoteLength(received / audio.sampleRate);

307297 },

@@ -330,7 +320,6 @@ audioPanel.addEventListener("toggle", async () => {

330320 await audio.enable();

331321 if (audio.sampleRate !== audioSpecRate) audioSpec = makeAudioSpec((audioSpecRate = audio.sampleRate));

332322 if (running && energy(P, sim.y) > 0) startNote();

333- else audioStatus.textContent = "Press Strike to play a note.";

334323 });

335324 // Speaker icon: crossed out when muted, with sound waves when audible.

336325 const SPEAKER = '<path d="M4 9h3l5-4v14l-5-4H4z" fill="currentColor" stroke="none"/>';

@@ -355,7 +344,6 @@ pitchInput.addEventListener("input", () => {

355344 });

356345 signalSelect.addEventListener("change", () => {

357346 signal = signalByKey(signalSelect.value);

358- showSpecLabel();

359347 showPitch();

360348 spectrogram.clear();

361349 audio.setSignal(signal.key);

@@ -460,6 +448,9 @@ document.getElementById("param-form")!.addEventListener("submit", (e) => {

460448 applyParams();

461449 });

462450

451+// Poincare section panel (own energy; same parameters, no damping).

452+const poincare = new PoincarePanel(P);

453+

463454 // Apply the parameters in the form: the pendulum returns to rest (waiting for

464455 // a Strike) and the strike energy keeps its fraction of the both-inverted

465456 // energy.

@@ -492,6 +483,7 @@ function applyParams() {

492483 audioSpec.setRefFreqs(audioModeHz());

493484 audio.setParams(P);

494485 audio.stopNote();

486+ poincare.setParams(P);

495487 }

496488

497489 // Keep the canvas backing store at device resolution.

@@ -595,6 +587,7 @@ function frame(now: number) {

595587 if (!audioSpecFig.hidden) audioSpec.draw();

596588 audioTraces.draw();

597589 if (!scatterFig.hidden) scatter.draw();

590+ poincare.draw();

598591 requestAnimationFrame(frame);

599592 }

600593 requestAnimationFrame(frame);

src/poincare-core.tsadded+104−0View file

@@ -0,0 +1,104 @@

1+// Poincare section of the undamped double pendulum at fixed energy E: the

2+// state is recorded each time the inner arm passes straight down moving right

3+// (th1 = 0, th1dot > 0), and plotted as the outer arm's (th2, p2).

4+//

5+// On the section, kinetic energy T = (b p1^2 - 2 c p1 p2 + a p2^2) / (2 det),

6+// with a = (m1 + m2) l1^2, b = m2 l2^2, c = m2 l1 l2 cos(th2), det = ab - c^2,

7+// must equal E - V(th2), V(th2) = m2 g l2 (1 - cos th2). Solving for p1 gives

8+// discriminant det (2 b (E - V) - p2^2), so a point is reachable exactly when

9+// p2^2 <= 2 b (E - V(th2)) (the outline of the allowed region), and the root

10+// with th1dot = sqrt(disc)/det > 0 is p1 = (c p2 + sqrt(disc)) / b.

11+

12+import { energy, type Params } from "./pendulum";

13+import { Lyapunov } from "./lyapunov";

14+import { Sim } from "./sim";

15+

16+export const CROSSINGS = 300; // crossings recorded per orbit

17+

18+const outerV = (P: Params, th2: number) => P.m2 * P.g * P.l2 * (1 - Math.cos(th2));

19+

20+// Largest |th2| on the section (pi if the outer arm can go over the top).

21+export function thetaMax(P: Params, E: number): number {

22+ const r = 1 - E / (P.m2 * P.g * P.l2);

23+ return r <= -1 ? Math.PI : Math.acos(Math.max(-1, Math.min(1, r)));

24+}

25+

26+// Largest |p2| on the section at angle th2 (0 outside the allowed range).

27+export function p2Max(P: Params, E: number, th2: number): number {

28+ return Math.sqrt(Math.max(0, 2 * P.m2 * P.l2 * P.l2 * (E - outerV(P, th2))));

29+}

30+

31+// Full state on the section, or null if (th2, p2) is not reachable at E.

32+export function sectionState(P: Params, E: number, th2: number, p2: number): Float64Array | null {

33+ const { m1, m2, l1, l2 } = P;

34+ const a = (m1 + m2) * l1 * l1, b = m2 * l2 * l2, c = m2 * l1 * l2 * Math.cos(th2);

35+ const det = a * b - c * c;

36+ const disc = det * (2 * b * (E - outerV(P, th2)) - p2 * p2);

37+ if (!(disc >= 0)) return null;

38+ return new Float64Array([0, th2, (c * p2 + Math.sqrt(disc)) / b, p2, 0]);

39+}

40+

41+// Starting points spread over the allowed region: an nx-by-ny grid in

42+// coordinates normalized to the region (u across th2, v across p2 at that th2).

43+export function seedGrid(P: Params, E: number, nx: number, ny: number): [number, number][] {

44+ const tm = thetaMax(P, E);

45+ const seeds: [number, number][] = [];

46+ for (let i = 0; i < nx; i++) {

47+ const th2 = tm * (-1 + (2 * (i + 0.5)) / nx);

48+ const pm = p2Max(P, E, th2);

49+ for (let j = 0; j < ny; j++) seeds.push([th2, pm * (-1 + (2 * (j + 0.5)) / ny)]);

50+ }

51+ return seeds;

52+}

53+

54+// One orbit followed from a point on the section, without damping, recording

55+// its crossings and measuring its Lyapunov exponent over the whole run.

56+export class SectionOrbit {

57+ readonly sim: Sim;

58+ readonly lyap = new Lyapunov(Infinity);

59+ points: number[]; // th2, p2 pairs (the start point first)

60+ crossings = 0;

61+ private h: number;

62+ private maxTime: number;

63+

64+ constructor(P: Params, E: number, h: number, th2: number, p2: number, maxTime: number) {

65+ const y = sectionState(P, E, th2, p2)!;

66+ this.sim = new Sim({ ...P, tau: 0 }, 0, 0.5);

67+ const E0 = energy(this.sim.P, y);

68+ this.sim.setState(y, E0, E0);

69+ this.lyap.start(this.sim);

70+ this.h = h;

71+ this.maxTime = maxTime;

72+ this.points = [th2, p2];

73+ }

74+

75+ get done(): boolean {

76+ return this.sim.t >= this.maxTime;

77+ }

78+

79+ // Integrate until the next crossing (or the time limit); returns whether a

80+ // crossing was recorded. The crossing is located by linear interpolation in

81+ // the step where th1 changes sign from negative to positive.

82+ advance(): boolean {

83+ const { sim, h } = this;

84+ const y = sim.y;

85+ while (sim.t < this.maxTime) {

86+ const th1 = y[0], th2 = y[1], p2 = y[3];

87+ sim.step(h);

88+ this.lyap.step(sim, h);

89+ if (th1 < 0 && y[0] >= 0 && y[0] - th1 < 1) {

90+ const f = -th1 / (y[0] - th1);

91+ const d2 = Math.atan2(Math.sin(y[1] - th2), Math.cos(y[1] - th2));

92+ const t2 = th2 + f * d2;

93+ this.points.push(Math.atan2(Math.sin(t2), Math.cos(t2)), p2 + f * (y[3] - p2));

94+ this.crossings++;

95+ return true;

96+ }

97+ }

98+ return false;

99+ }

100+

101+ lambda(): number {

102+ return this.lyap.mean();

103+ }

104+}

src/poincare-worker.tsadded+90−0View file

@@ -0,0 +1,90 @@

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 { stepSize, topEnergy, type Params } from "./pendulum";

8+import { CROSSINGS, SectionOrbit } from "./poincare-core";

9+

10+const BATCH_MS = 40;

11+// Step: 4x the animation's (which keeps the energy error below 1e-6). Over

12+// 300 crossings the energy error stays below ~2e-5 and the regular/chaotic

13+// classification is unchanged; regular orbits move by ~1e-4, far below a pixel.

14+const STEP_FACTOR = 4;

15+const MAX_TIME = 2000; // s; an orbit that stops crossing (e.g. inner arm rotating left) ends here

16+

17+export type ToPoincare =

18+ | { type: "run"; id: number; P: Params; E: number; seeds: [number, number][] }

19+ | { type: "add"; id: number; seed: [number, number] }

20+ | { type: "stop" };

21+

22+export type FromPoincare = {

23+ type: "points";

24+ id: number;

25+ orbits: { index: number; points: Float32Array; lambda: number; done: boolean }[];

26+ finished: boolean;

27+};

28+

29+let runId = -1;

30+let P: Params;

31+let E = 0;

32+let h = 0;

33+let hKey = "";

34+let orbits: SectionOrbit[] = [];

35+let sent: number[] = []; // points (pairs) already posted, per orbit

36+let doneSent: boolean[] = []; // whether the orbit's end has been posted

37+let running = false;

38+

39+onmessage = (e: MessageEvent<ToPoincare>) => {

40+ const msg = e.data;

41+ if (msg.type === "run") {

42+ runId = msg.id;

43+ P = { ...msg.P, tau: 0 };

44+ E = msg.E;

45+ const key = JSON.stringify([P.m1, P.m2, P.l1, P.l2, P.g]);

46+ if (key !== hKey) {

47+ h = STEP_FACTOR * stepSize(P, 2 * topEnergy(P));

48+ hKey = key;

49+ }

50+ orbits = msg.seeds.map(([a, b]) => new SectionOrbit(P, E, h, a, b, MAX_TIME));

51+ sent = orbits.map(() => 0);

52+ doneSent = orbits.map(() => false);

53+ } else if (msg.type === "add") {

54+ if (msg.id !== runId) return;

55+ orbits.push(new SectionOrbit(P, E, h, msg.seed[0], msg.seed[1], MAX_TIME));

56+ sent.push(0);

57+ doneSent.push(false);

58+ } else {

59+ orbits = [];

60+ sent = [];

61+ doneSent = [];

62+ }

63+ if (!running && orbits.length) {

64+ running = true;

65+ setTimeout(work, 0);

66+ }

67+};

68+

69+const active = (o: SectionOrbit) => o.crossings < CROSSINGS && !o.done;

70+

71+function work(): void {

72+ const t0 = performance.now();

73+ while (performance.now() - t0 < BATCH_MS && orbits.some(active)) {

74+ for (const o of orbits) if (active(o)) o.advance();

75+ }

76+ const out: FromPoincare["orbits"] = [];

77+ orbits.forEach((o, index) => {

78+ const n = o.points.length / 2;

79+ const done = !active(o);

80+ if (n > sent[index] || (done && !doneSent[index])) {

81+ out.push({ index, points: Float32Array.from(o.points.slice(2 * sent[index])), lambda: o.lambda(), done });

82+ sent[index] = n;

83+ doneSent[index] = done;

84+ }

85+ });

86+ const finished = !orbits.some(active);

87+ postMessage({ type: "points", id: runId, orbits: out, finished } satisfies FromPoincare);

88+ if (finished) running = false;

89+ else setTimeout(work, 0);

90+}

src/poincare.tsadded+266−0View file

@@ -0,0 +1,266 @@

1+// Poincare section panel: its own energy slider (same parameters as the

2+// pendulum above, no damping), the allowed region's outline, orbits seeded on

3+// a grid and computed in a worker, click to add an orbit, hover to highlight

4+// one, and coloring per orbit or by Lyapunov exponent.

5+

6+import { topEnergy, type Params } from "./pendulum";

7+import { p2Max, sectionState, seedGrid, thetaMax } from "./poincare-core";

8+import type { FromPoincare, ToPoincare } from "./poincare-worker";

9+

10+const WIDTH = 960, HEIGHT = 480; // CSS pixels

11+const TOP = 18; // legend row

12+const PAD = 10;

13+const FONT = "11px ui-monospace, SFMono-Regular, Menlo, monospace";

14+const GRID = [6, 7]; // seed grid: across th2, across p2

15+// Lambda coloring is logarithmic between these (1/s): regular orbits measure

16+// ~0.01 or less over a run (lambda -> 0 as the run grows), weakly chaotic ones

17+// ~0.2 to 0.3, strongly chaotic ones ~1 to 1.5.

18+const LAMBDA_LO = 0.01, LAMBDA_HI = 1.5;

19+

20+interface Orbit {

21+ points: number[]; // th2, p2 pairs

22+ lambda: number;

23+ done: boolean;

24+ clicked: boolean;

25+}

26+

27+export class PoincarePanel {

28+ private P: Params;

29+ private E: number;

30+ private canvas: HTMLCanvasElement;

31+ private dpr = window.devicePixelRatio || 1;

32+ private worker: Worker | null = null;

33+ private runId = 0;

34+ private orbits: Orbit[] = [];

35+ private finished = false;

36+ private changed = true;

37+ private hoverOrbit = -1;

38+ private hoverPoint: [number, number] | null = null;

39+ private byLambda = false;

40+ private open = false;

41+ private energyInput: HTMLInputElement;

42+ private energyValue: HTMLElement;

43+ private status: HTMLElement;

44+ private debounce = 0;

45+ // Axes (data ranges), set from E.

46+ private thMax = 1;

47+ private pMax = 1;

48+

49+ constructor(P: Params) {

50+ this.P = P;

51+ this.energyInput = document.getElementById("pc-energy") as HTMLInputElement;

52+ this.energyValue = document.getElementById("pc-energy-value")!;

53+ this.status = document.getElementById("pc-status")!;

54+ this.E = Number(this.energyInput.value);

55+ this.canvas = document.createElement("canvas");

56+ this.canvas.width = WIDTH * this.dpr;

57+ this.canvas.height = HEIGHT * this.dpr;

58+ this.canvas.style.cursor = "crosshair";

59+ document.getElementById("pc-plot")!.appendChild(this.canvas);

60+

61+ const panel = document.getElementById("poincare") as HTMLDetailsElement;

62+ panel.addEventListener("toggle", () => {

63+ this.open = panel.open;

64+ if (this.open) this.run();

65+ else this.send({ type: "stop" });

66+ });

67+ this.energyInput.addEventListener("input", () => {

68+ this.E = Number(this.energyInput.value);

69+ this.showEnergy();

70+ // Recompute shortly after the slider stops moving.

71+ clearTimeout(this.debounce);

72+ this.debounce = window.setTimeout(() => this.run(), 150);

73+ });

74+ const lambdaBox = document.getElementById("pc-lambda") as HTMLInputElement;

75+ lambdaBox.addEventListener("change", () => {

76+ this.byLambda = lambdaBox.checked;

77+ this.changed = true;

78+ });

79+ this.canvas.addEventListener("mousemove", (e) => this.onHover(e));

80+ this.canvas.addEventListener("mouseleave", () => {

81+ this.hoverOrbit = -1;

82+ this.hoverPoint = null;

83+ this.changed = true;

84+ });

85+ this.canvas.addEventListener("click", (e) => this.onClick(e));

86+ this.setParams(P);

87+ }

88+

89+ // New pendulum parameters: rescale the energy slider (keeping its fraction

90+ // of the both-inverted energy) and recompute if the panel is open.

91+ setParams(P: Params): void {

92+ const frac = this.E / topEnergy(this.P);

93+ this.P = P;

94+ const eMax = 2 * topEnergy(P);

95+ this.energyInput.max = String(eMax);

96+ this.energyInput.step = String(eMax / 1200);

97+ this.E = Math.min(eMax, frac * topEnergy(P));

98+ this.energyInput.value = String(this.E);

99+ this.showEnergy();

100+ if (this.open) this.run();

101+ }

102+

103+ private showEnergy(): void {

104+ this.energyValue.textContent = `${this.E.toFixed(1)} J`;

105+ }

106+

107+ private send(msg: ToPoincare): void {

108+ if (!this.worker) {

109+ this.worker = new Worker(new URL("./poincare-worker.ts", import.meta.url), { type: "module" });

110+ this.worker.onmessage = (e: MessageEvent<FromPoincare>) => this.receive(e.data);

111+ }

112+ this.worker.postMessage(msg);

113+ }

114+

115+ // Start computing the section at the current energy.

116+ private run(): void {

117+ this.runId++;

118+ this.orbits = [];

119+ this.finished = false;

120+ this.thMax = thetaMax(this.P, this.E);

121+ this.pMax = p2Max(this.P, this.E, 0);

122+ const seeds = this.E > 0 ? seedGrid(this.P, this.E, GRID[0], GRID[1]) : [];

123+ for (let i = 0; i < seeds.length; i++) this.orbits.push({ points: [], lambda: NaN, done: false, clicked: false });

124+ if (seeds.length) this.send({ type: "run", id: this.runId, P: this.P, E: this.E, seeds });

125+ this.changed = true;

126+ }

127+

128+ private receive(msg: FromPoincare): void {

129+ if (msg.id !== this.runId) return; // from an earlier energy

130+ for (const o of msg.orbits) {

131+ const orbit = this.orbits[o.index];

132+ if (!orbit) continue;

133+ for (const v of o.points) orbit.points.push(v);

134+ orbit.lambda = o.lambda;

135+ orbit.done = o.done;

136+ }

137+ this.finished = msg.finished;

138+ this.changed = true;

139+ }

140+

141+ // Data <-> CSS pixel coordinates.

142+ private px(th2: number): number {

143+ return PAD + ((th2 + this.thMax * 1.04) / (2 * this.thMax * 1.04)) * (WIDTH - 2 * PAD);

144+ }

145+ private py(p2: number): number {

146+ return TOP + PAD + ((this.pMax * 1.06 - p2) / (2 * this.pMax * 1.06)) * (HEIGHT - TOP - 2 * PAD);

147+ }

148+ private fromPx(x: number, y: number): [number, number] {

149+ const th2 = ((x - PAD) / (WIDTH - 2 * PAD)) * 2 * this.thMax * 1.04 - this.thMax * 1.04;

150+ const p2 = this.pMax * 1.06 - ((y - TOP - PAD) / (HEIGHT - TOP - 2 * PAD)) * 2 * this.pMax * 1.06;

151+ return [th2, p2];

152+ }

153+ private mouse(e: MouseEvent): [number, number] {

154+ const r = this.canvas.getBoundingClientRect();

155+ return [((e.clientX - r.left) / r.width) * WIDTH, ((e.clientY - r.top) / r.height) * HEIGHT];

156+ }

157+

158+ private onHover(e: MouseEvent): void {

159+ const [x, y] = this.mouse(e);

160+ this.hoverPoint = this.fromPx(x, y);

161+ let best = 8 * 8; // within 8 px

162+ this.hoverOrbit = -1;

163+ this.orbits.forEach((o, k) => {

164+ for (let i = 0; i < o.points.length; i += 2) {

165+ const dx = this.px(o.points[i]) - x, dy = this.py(o.points[i + 1]) - y;

166+ if (dx * dx + dy * dy < best) {

167+ best = dx * dx + dy * dy;

168+ this.hoverOrbit = k;

169+ }

170+ }

171+ });

172+ this.changed = true;

173+ }

174+

175+ // Add an orbit starting at the clicked point, if it is reachable.

176+ private onClick(e: MouseEvent): void {

177+ const [x, y] = this.mouse(e);

178+ const [th2, p2] = this.fromPx(x, y);

179+ if (!this.worker || !sectionState(this.P, this.E, th2, p2)) return;

180+ this.orbits.push({ points: [], lambda: NaN, done: false, clicked: true });

181+ this.finished = false;

182+ this.send({ type: "add", id: this.runId, seed: [th2, p2] });

183+ this.changed = true;

184+ }

185+

186+ private color(k: number, alpha: number): string {

187+ const o = this.orbits[k];

188+ if (this.byLambda) {

189+ if (!Number.isFinite(o.lambda)) return `rgba(150, 150, 150, ${alpha})`;

190+ // Blue (regular) -> amber -> red (chaotic).

191+ // (A regular orbit's lambda can come out slightly negative.)

192+ const u = o.lambda <= LAMBDA_LO ? 0 : Math.min(1, Math.log(o.lambda / LAMBDA_LO) / Math.log(LAMBDA_HI / LAMBDA_LO));

193+ const ramp = [[47, 93, 138], [227, 162, 26], [194, 65, 45]];

194+ const s = u * 2, i = Math.min(1, Math.floor(s)), w = s - i;

195+ const c = ramp[i].map((a, j) => Math.round(a + (ramp[i + 1][j] - a) * w));

196+ return `rgba(${c[0]}, ${c[1]}, ${c[2]}, ${alpha})`;

197+ }

198+ const hue = (k * 137.508) % 360; // golden angle: neighbouring orbits differ

199+ return `hsla(${hue}, 45%, 42%, ${alpha})`;

200+ }

201+

202+ draw(): void {

203+ if (!this.open || !this.changed) return;

204+ this.changed = false;

205+ const ctx = this.canvas.getContext("2d")!;

206+ ctx.setTransform(this.dpr, 0, 0, this.dpr, 0, 0);

207+ ctx.clearRect(0, 0, WIDTH, HEIGHT);

208+

209+ // Axes through 0 and the allowed region's outline.

210+ ctx.strokeStyle = "#ecece6";

211+ ctx.lineWidth = 1;

212+ ctx.beginPath();

213+ ctx.moveTo(this.px(0), TOP);

214+ ctx.lineTo(this.px(0), HEIGHT);

215+ ctx.moveTo(0, this.py(0));

216+ ctx.lineTo(WIDTH, this.py(0));

217+ ctx.stroke();

218+ ctx.strokeStyle = "#b9b8b0";

219+ ctx.beginPath();

220+ for (const sign of [1, -1]) {

221+ for (let i = 0; i <= 200; i++) {

222+ const th2 = -this.thMax + (2 * this.thMax * i) / 200;

223+ const x = this.px(th2), y = this.py(sign * p2Max(this.P, this.E, th2));

224+ if (i === 0) ctx.moveTo(x, y);

225+ else ctx.lineTo(x, y);

226+ }

227+ }

228+ ctx.stroke();

229+

230+ // Orbits; the hovered one on top and the others faded while hovering.

231+ const order = this.orbits.map((_, k) => k).filter((k) => k !== this.hoverOrbit);

232+ if (this.hoverOrbit >= 0) order.push(this.hoverOrbit);

233+ for (const k of order) {

234+ const o = this.orbits[k];

235+ const hovered = k === this.hoverOrbit;

236+ ctx.fillStyle = this.color(k, this.hoverOrbit >= 0 && !hovered ? 0.25 : 0.85);

237+ const size = hovered || o.clicked ? 2.5 : 1.8;

238+ for (let i = 0; i < o.points.length; i += 2) {

239+ ctx.fillRect(this.px(o.points[i]) - size / 2, this.py(o.points[i + 1]) - size / 2, size, size);

240+ }

241+ }

242+

243+ // Legend row.

244+ ctx.font = FONT;

245+ ctx.textBaseline = "top";

246+ ctx.fillStyle = "#2f5d8a";

247+ let legend: string;

248+ if (this.hoverOrbit >= 0) {

249+ const o = this.orbits[this.hoverOrbit];

250+ const kind = !Number.isFinite(o.lambda) ? "" : o.lambda > 0.1 ? "chaotic" : "regular";

251+ legend = `orbit ${this.hoverOrbit + 1} · λ ${Number.isFinite(o.lambda) ? o.lambda.toFixed(3) : "–"} 1/s ${kind ? `(${kind})` : ""}`;

252+ } else if (this.hoverPoint) {

253+ const [th2, p2] = this.hoverPoint;

254+ legend = `θ2 ${th2.toFixed(3)} rad · p2 ${p2.toFixed(3)} kg m²/s${sectionState(this.P, this.E, th2, p2) ? " · click to add an orbit here" : ""}`;

255+ } else {

256+ legend = `θ2 (rad) across · p2 (kg m²/s) up · ${this.orbits.length} orbits`;

257+ }

258+ ctx.fillText(legend, 4, 3);

259+ ctx.fillStyle = "#999";

260+ ctx.textAlign = "right";

261+ ctx.fillText(`θ2 ±${this.thMax.toFixed(2)} · p2 ±${this.pMax.toFixed(2)}`, WIDTH - 4, 3);

262+ ctx.textAlign = "left";

263+

264+ this.status.textContent = this.finished ? "" : "computing…";

265+ }

266+}

src/scatter.tsmodified+1−1View file

@@ -181,7 +181,7 @@ export class ScatterPlot {

181181 const legend =

182182 near >= 0

183183 ? `E ${E[near].toFixed(2)} J λ ${L[near].toFixed(3)} 1/s ${T[near].toFixed(2)} s into note`

184- : `λ (1/s, 10 s window) vs E (J) · ${this.noteLength.length} note${this.noteLength.length === 1 ? "" : "s"}`;

184+ : "λ (1/s) vs E (J)";

185185 ctx.fillText(legend, 4, 3);

186186 ctx.fillStyle = "#999";

187187 ctx.textAlign = "right";