Two-oscillator chaos toy: model, experiments, figures, web app
10 changed files+1516−0
.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+__pycache__/
2+*.npz
README.mdadded+30−0View file
@@ -0,0 +1,30 @@
1+# chaos-toy-examples
2+
3+Simple systems for studying the kind of chaos found in the struck plate of `cymbal-elastodynamics/flat-plate-2`, where it is entangled with everything else a plate does.
4+
5+## The phenomenon
6+
7+In the plate, a weak strike gives regular, quasi-periodic vibration: two runs that differ by a tiny perturbation stay a tiny distance apart. Above a threshold in strike strength the motion becomes chaotic: such differences grow exponentially (e-folding in a few milliseconds at the strongest strike), energy cascades to high frequencies, and after a few tenths of a second two runs have unrelated waveforms. As damping drains the energy the plate returns to the regular regime, and the energy distribution among the modes that the chaotic phase left behind is frozen in: each strong strike rings out with its own pattern of partial levels, which is audible.
8+
9+Structurally the plate in modal form is a set of oscillators coupled through a quartic potential (the stretching energy), lightly damped and started by an impulse. The toys here keep that structure and nothing else.
10+
11+## Toys
12+
13+| Directory | System | Purpose |
14+| --- | --- | --- |
15+| `two-oscillators/` | two oscillators with a quartic coupling, $H = \tfrac12 (p_x^2 + p_y^2) + \tfrac12 (\omega_1^2 x^2 + \omega_2^2 y^2) + \beta x^2 y^2$ | the mechanism: how regular motion breaks up as the energy grows (Poincare sections, Lyapunov exponent vs energy), and why damping freezes a random outcome |
16+| `../fput-chain` (separate project) | the Fermi-Pasta-Ulam-Tsingou beta chain, struck, with frequency-dependent damping | the many-mode phenomena: energy cascade or not, threshold, the frozen spectrum, listening; web app at https://claude.ai/artifact/C6oaNiG2zCRX9VR7Wn4Wfx |
17+
18+With $\beta$ fixed, rescaling the variables by $\sqrt{E}$ shows that the two-oscillator dynamics depends on $\beta E$ only, the analogue of the plate's strike strength (or peak deflection over thickness).
19+
20+## Web application
21+
22+`web/index.html` is a single-page app (no build step) that simulates the two oscillators in the browser and plays them: strike energy, frequency ratio, strike direction, ring time, pitch and the difference between strikes A and B are adjustable; it shows the path in the (x, y) plane, the Poincare section of the undamped system at the strike energy, the energy share of 12 nearby strikes over time, the local Lyapunov exponent with the energy, and a spectrogram of any of the 12 strikes (a selector, or the arrow keys, flips between strikes on a shared color scale, draws the chosen strike's path in the bowl and highlights it in the energy-share plot, and plays it), with a playhead while a strike plays. Time is scaled so that the x oscillator sounds at the chosen pitch (220 Hz by default), and the sound is the pickup $p_x + p_y$. Published as a private artifact: https://claude.ai/artifact/BNxuYsQoFwNavTWnF7EHXg. Locally: `cd web && python -m http.server`.
23+
24+## Running
25+
26+Experiments run on the office workstation (`~/work/chaos-toy-examples`, Python 3.9 with numpy), plots locally:
27+
28+ cd two-oscillators
29+ WORKERS=24 python experiments.py all # results/*.npz
30+ python plots.py # results/*.png
two-oscillators/README.mdadded+30−0View file
@@ -0,0 +1,30 @@
1+# Two oscillators with a quartic coupling
2+
3+$$H = \tfrac12 (p_x^2 + p_y^2) + \tfrac12 (\omega_1^2 x^2 + \omega_2^2 y^2) + \beta x^2 y^2, \qquad \omega_1 = 1, \ \beta = 1,$$
4+
5+with $\omega_2 = 1$ (equal frequencies, the system of Pullen and Edmonds 1981) or $\omega_2 = 1.618$ (an inharmonic pair, closer to the plate). Since rescaling the variables by $\sqrt{E}$ leaves only $\beta E$, the energy $E$ is the single parameter, the analogue of the plate's strike strength.
6+
7+`model.py` integrates many orbits at once with a fourth-order symplectic scheme (Yoshida's composition of leapfrog steps; energy conserved to 1e-6 or better over 1000 time units), optionally with damping $p \to p\, e^{-\gamma \Delta t}$, and integrates the tangent dynamics for the largest Lyapunov exponent (Benettin's renormalization). An orbit counts as chaotic when its finite-time exponent over 2000 to 3000 time units exceeds 0.02; regular orbits give values near $\log(t)/t \approx 0.005$.
8+
9+## Results (2026-10-05)
10+
11+**Poincare sections** (`results/sections.png`; crossings of $y = 0$ with $p_y > 0$, in the $(x, p_x)$ plane). At low energy every orbit lies on a closed curve (a torus): the motion is quasi-periodic. As the energy grows, chaotic layers appear around the separatrices between families of tori, widen, and merge, until only small islands of regular motion remain.
12+
13+**Onset of chaos** (`results/lyapunov.png`; 200 orbits per energy, 25 energies). Fraction of orbits that are chaotic:
14+
15+| E | 0.42 | 0.56 | 0.75 | 1.0 | 1.3 | 1.8 | 2.4 | 3.2 | 5.6 | 10 |
16+| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
17+| $\omega_2 = 1$ | 0 | 0.03 | 0.14 | 0.48 | 0.69 | 0.87 | 0.91 | 0.94 | 1.00 | 1.00 |
18+| $\omega_2 = 1.618$ | 0 | 0 | 0 | 0 | 0 | 0.17 | 0.60 | 0.81 | 0.96 | 0.98 |
19+
20+The transition is gradual, not a sharp threshold, and the inharmonic pair needs about twice the energy before it starts: with equal frequencies the 1:1 resonance couples the two oscillators strongly from the start, while detuned oscillators exchange energy only through higher-order resonances. In the chaotic regime the Lyapunov exponent keeps growing with energy (median about 0.3 at E = 5 to 1.4 at E = 100 for $\omega_2 = 1.618$): stronger excitation, faster separation, as in the plate.
21+
22+**A damped strike and its frozen outcome** (`results/damped.png`; $\omega_2 = 1.618$, $\gamma = 0.005$; damping the momenta by $e^{-\gamma \Delta t}$ makes the energy decay as $e^{-\gamma t}$). Both oscillators start at rest with momentum $\sqrt{2 E_0}$ in a fixed direction; an ensemble of 16 strikes differs in that direction by $10^{-10}$ times the member index. At $E_0 = 0.5$ and $E_0 = 3$ the strike lands on a regular orbit: the local Lyapunov exponent decays to zero, and all 16 strikes end with the same split of the energy between the two oscillators (0.557 and 0.398 in the x oscillator). At $E_0 = 20$ the motion is chaotic for the first 300 to 600 time units (local exponent 0.4 to 0.8), the 16 strikes separate completely, and as the energy falls below about 1 each is captured by a regular orbit and keeps the split it happened to have. At the end the share of energy in the x oscillator ranges from 0.05 to 0.95 across strikes that differed by $10^{-10}$ (standard deviation 0.31).
23+
24+This is the plate's frozen fingerprint in its simplest form: a transient chaotic phase while the energy is high, capture onto a torus once damping has taken the energy below the onset of chaos, and a final state that is a random sample from a distribution, set by perturbations far below anything measurable. Note that $E_0 = 3$ lies where 81% of section orbits are chaotic, yet this particular strike is regular: whether a strike is chaotic depends on where in phase space it starts, not only on its energy.
25+
26+## Possible next steps
27+
28+- The distribution of the frozen split as a function of $E_0$ and of the damping rate (slower damping gives the chaotic phase longer to mix).
29+- Damping that grows with frequency (as in the plate), which makes the higher oscillator lose energy faster and biases the outcome.
30+- The many-mode version (FPUT chain) for the cascade.
two-oscillators/experiments.pyadded+104−0View file
@@ -0,0 +1,104 @@
1+"""Experiments with the two quartic-coupled oscillators. Each writes results/<name>.npz; plots.py draws them.
2+
3+ python experiments.py sections # Poincare sections at several energies (w2 = 1 and 1.618)
4+ python experiments.py lyapunov # finite-time largest Lyapunov exponent vs energy
5+ python experiments.py damped # a damped 'strike' and an ensemble of nearby strikes
6+ python experiments.py all # all three, in parallel processes
7+
8+beta = 1 throughout, so the energy E is the single parameter (beta E).
9+"""
10+
11+from __future__ import annotations
12+
13+import os
14+import sys
15+from multiprocessing import Pool
16+
17+os.environ.setdefault("OMP_NUM_THREADS", "1")
18+
19+import numpy as np # noqa: E402
20+
21+from model import System, run, section_initial_conditions # noqa: E402
22+
23+HERE = os.path.dirname(os.path.abspath(__file__))
24+RES = os.path.join(HERE, "results")
25+W2S = (1.0, 1.618)
26+
27+
28+def dt_for(E):
29+ """Step small enough for the stiffening at energy E (the effective frequency grows like sqrt(E))."""
30+ return min(0.01, 0.1 / np.sqrt(1 + 4 * E))
31+
32+
33+def _sections_one(args):
34+ w2, E = args
35+ sys_ = System(1.0, w2, 1.0)
36+ s0 = section_initial_conditions(sys_, E, 9)
37+ r = run(sys_, s0, 3000.0, dt=dt_for(E), lyap=True, section=True)
38+ pts = r["section"]
39+ return w2, E, r["lyap"], np.concatenate([np.full(len(p), i) for i, p in enumerate(pts)]), \
40+ np.concatenate([p for p in pts if len(p)] or [np.zeros((0, 2))])
41+
42+
43+def sections(pool):
44+ energies = (0.25, 0.5, 1.0, 1.5, 2.0, 4.0, 8.0)
45+ out = {}
46+ for w2, E, lam, owner, pts in pool.imap_unordered(_sections_one, [(w2, E) for w2 in W2S for E in energies]):
47+ key = f"w2_{w2:g}_E_{E:g}"
48+ out[key + "_lyap"], out[key + "_owner"], out[key + "_pts"] = lam, owner, pts
49+ print(f" sections w2 {w2:g} E {E:g}: {len(lam)} orbits", flush=True)
50+ np.savez(os.path.join(RES, "sections.npz"), energies=np.array(energies), w2s=np.array(W2S), **out)
51+
52+
53+def _lyap_one(args):
54+ w2, E = args
55+ sys_ = System(1.0, w2, 1.0)
56+ s0 = section_initial_conditions(sys_, E, 16)
57+ r = run(sys_, s0, 2000.0, dt=dt_for(E), lyap=True)
58+ return w2, E, r["lyap"]
59+
60+
61+def lyapunov(pool):
62+ energies = np.geomspace(0.1, 100, 25)
63+ lam = {}
64+ for w2, E, l in pool.imap_unordered(_lyap_one, [(w2, E) for w2 in W2S for E in energies]):
65+ lam[(w2, E)] = l
66+ print(f" lyapunov w2 {w2:g} E {E:.3g}: median {np.median(l):.3f}", flush=True)
67+ np.savez(os.path.join(RES, "lyapunov.npz"), energies=energies, w2s=np.array(W2S),
68+ **{f"w2_{w2:g}_E{i}": lam[(w2, E)] for w2 in W2S for i, E in enumerate(energies)})
69+
70+
71+def _damped_one(args):
72+ """A 'strike': both oscillators start at rest at the origin with momentum sqrt(2 E0) in direction a,
73+ then decay with p -> p exp(-gamma dt). An ensemble perturbs a by i * 1e-10. Records the mode energies
74+ and the running Lyapunov log-sum of member 0."""
75+ E0, gamma, n_members = args
76+ sys_ = System(1.0, 1.618, 1.0)
77+ a = 0.7 + 1e-10 * np.arange(n_members)
78+ p = np.sqrt(2 * E0)
79+ s0 = np.array([np.zeros(n_members), np.zeros(n_members), p * np.cos(a), p * np.sin(a)])
80+ r = run(sys_, s0, 1500.0, dt=dt_for(E0), gamma=gamma, lyap=True, renorm_every=10,
81+ record_every=int(round(0.5 / dt_for(E0))))
82+ Ex, Ey = sys_.mode_energies(np.moveaxis(r["states"], 1, 0))
83+ return E0, r["t"], Ex, Ey, sys_.energy(np.moveaxis(r["states"], 1, 0)), r["logsum"]
84+
85+
86+def damped(pool):
87+ gamma = 0.005
88+ out = {}
89+ for E0, t, Ex, Ey, Etot, logsum in pool.imap_unordered(_damped_one, [(E0, gamma, 16) for E0 in (0.5, 3.0, 20.0)]):
90+ k = f"E0_{E0:g}"
91+ out[k + "_t"], out[k + "_Ex"], out[k + "_Ey"], out[k + "_E"], out[k + "_logsum"] = t, Ex, Ey, Etot, logsum
92+ print(f" damped E0 {E0:g} done", flush=True)
93+ np.savez(os.path.join(RES, "damped.npz"), gamma=gamma, **out)
94+
95+
96+if __name__ == "__main__":
97+ os.makedirs(RES, exist_ok=True)
98+ which = sys.argv[1:] or ["all"]
99+ if "all" in which:
100+ which = ["sections", "lyapunov", "damped"]
101+ with Pool(int(os.environ.get("WORKERS", "16"))) as pool:
102+ for w in which:
103+ print(w, flush=True)
104+ dict(sections=sections, lyapunov=lyapunov, damped=damped)[w](pool)
two-oscillators/model.pyadded+158−0View file
@@ -0,0 +1,158 @@
1+"""Two oscillators coupled by a quartic potential, vectorized over many orbits.
2+
3+ H = (px^2 + py^2) / 2 + (w1^2 x^2 + w2^2 y^2) / 2 + beta x^2 y^2
4+
5+With w1 = w2 this is the system studied by Pullen and Edmonds (1981). Scaling x, y, px, py by sqrt(E)
6+shows that for fixed w1, w2 the dynamics depends only on beta E, the analogue of the strike strength
7+of the plate. The forces are
8+
9+ Fx = -w1^2 x - 2 beta x y^2, Fy = -w2^2 y - 2 beta x^2 y.
10+
11+Integration: 4th-order symplectic composition (Yoshida 1990) of the kick-drift-kick leapfrog, which
12+conserves energy to high accuracy over long times. Optional linear damping p -> p exp(-gamma dt) after
13+each step (dp/dt = -gamma p: the amplitude decays at rate gamma / 2 and the energy at rate gamma,
14+the analogue of the plate's sigma0). The tangent
15+dynamics (for Lyapunov exponents) is integrated with the same splitting, linearized.
16+"""
17+
18+from __future__ import annotations
19+
20+from dataclasses import dataclass
21+
22+import numpy as np
23+
24+_C = 2 ** (1 / 3)
25+YOSHIDA = (1 / (2 - _C), -_C / (2 - _C), 1 / (2 - _C))
26+
27+
28+@dataclass(frozen=True)
29+class System:
30+ w1: float = 1.0
31+ w2: float = 1.0
32+ beta: float = 1.0
33+
34+ def force(self, x, y):
35+ return -self.w1**2 * x - 2 * self.beta * x * y * y, -self.w2**2 * y - 2 * self.beta * x * x * y
36+
37+ def tangent_force(self, x, y, dx, dy):
38+ """Linearized force applied to a tangent vector (dx, dy) at (x, y)."""
39+ b = self.beta
40+ fxx = -self.w1**2 - 2 * b * y * y
41+ fyy = -self.w2**2 - 2 * b * x * x
42+ fxy = -4 * b * x * y
43+ return fxx * dx + fxy * dy, fxy * dx + fyy * dy
44+
45+ def energy(self, s):
46+ x, y, px, py = s
47+ return 0.5 * (px**2 + py**2) + 0.5 * (self.w1**2 * x**2 + self.w2**2 * y**2) + self.beta * x**2 * y**2
48+
49+ def mode_energies(self, s):
50+ """Linear (uncoupled) energies of the two oscillators."""
51+ x, y, px, py = s
52+ return 0.5 * (px**2 + self.w1**2 * x**2), 0.5 * (py**2 + self.w2**2 * y**2)
53+
54+
55+def _leapfrog(sys, s, d, h):
56+ """One kick-drift-kick step of size h on state s (4, M) and tangent d (4, M) or None."""
57+ x, y, px, py = s
58+ fx, fy = sys.force(x, y)
59+ px = px + 0.5 * h * fx
60+ py = py + 0.5 * h * fy
61+ if d is not None:
62+ dx, dy, dpx, dpy = d
63+ gx, gy = sys.tangent_force(x, y, dx, dy)
64+ dpx = dpx + 0.5 * h * gx
65+ dpy = dpy + 0.5 * h * gy
66+ dx = dx + h * dpx
67+ dy = dy + h * dpy
68+ x = x + h * px
69+ y = y + h * py
70+ fx, fy = sys.force(x, y)
71+ px = px + 0.5 * h * fx
72+ py = py + 0.5 * h * fy
73+ if d is not None:
74+ gx, gy = sys.tangent_force(x, y, dx, dy)
75+ dpx = dpx + 0.5 * h * gx
76+ dpy = dpy + 0.5 * h * gy
77+ d = np.array([dx, dy, dpx, dpy])
78+ return np.array([x, y, px, py]), d
79+
80+
81+def step(sys, s, d, dt, gamma=0.0):
82+ for c in YOSHIDA:
83+ s, d = _leapfrog(sys, s, d, c * dt)
84+ if gamma:
85+ f = np.exp(-gamma * dt)
86+ s[2:] *= f
87+ if d is not None:
88+ d[2:] *= f
89+ return s, d
90+
91+
92+def run(sys, s0, t_end, dt=0.01, gamma=0.0, lyap=False, renorm_every=10, record_every=None,
93+ section=False, d0=None):
94+ """Integrate orbits s0 (4, M) to t_end.
95+
96+ lyap: integrate a tangent vector and return the finite-time largest Lyapunov exponent of each orbit
97+ (Benettin: renormalize every renorm_every steps and average the log growth), plus the running
98+ log-growth sum at each record time if record_every is given.
99+ section: record crossings of y = 0 with py > 0 (Poincare section), as (x, px) per orbit.
100+ record_every: record the state every that many steps.
101+ """
102+ s = np.array(s0, dtype=float)
103+ M = s.shape[1]
104+ n = int(round(t_end / dt))
105+ d = None
106+ logsum = np.zeros(M)
107+ if lyap:
108+ if d0 is None:
109+ rng = np.random.default_rng(0)
110+ d = rng.standard_normal((4, M))
111+ else:
112+ d = np.array(d0, dtype=float)
113+ d /= np.linalg.norm(d, axis=0)
114+ rec_t, rec_s, rec_log = [], [], []
115+ sec = [[] for _ in range(M)] if section else None
116+ for i in range(n):
117+ y_old = s[1].copy()
118+ s_old = s.copy()
119+ s, d = step(sys, s, d, dt, gamma)
120+ if section:
121+ cross = np.nonzero((y_old < 0) & (s[1] >= 0) & (s[3] > 0))[0]
122+ for m in cross:
123+ a = -y_old[m] / (s[1, m] - y_old[m])
124+ sec[m].append(((1 - a) * s_old[0, m] + a * s[0, m], (1 - a) * s_old[2, m] + a * s[2, m]))
125+ if lyap and (i + 1) % renorm_every == 0:
126+ nrm = np.linalg.norm(d, axis=0)
127+ logsum += np.log(nrm)
128+ d /= nrm
129+ if record_every and (i + 1) % record_every == 0:
130+ rec_t.append((i + 1) * dt)
131+ rec_s.append(s.copy())
132+ if lyap:
133+ rec_log.append(logsum.copy())
134+ out = dict(s=s)
135+ if lyap:
136+ out["lyap"] = logsum / (n * dt)
137+ if record_every:
138+ out["t"] = np.array(rec_t)
139+ out["states"] = np.array(rec_s) # (T, 4, M)
140+ if lyap:
141+ out["logsum"] = np.array(rec_log) # (T, M)
142+ if section:
143+ out["section"] = [np.array(p) for p in sec]
144+ return out
145+
146+
147+def section_initial_conditions(sys, E, n_grid):
148+ """Initial conditions on the section y = 0, py > 0 at energy E: a grid over the allowed (x, px) ellipse."""
149+ xm = np.sqrt(2 * E) / sys.w1
150+ pm = np.sqrt(2 * E)
151+ xs = np.linspace(-xm, xm, n_grid + 2)[1:-1]
152+ ps = np.linspace(-pm, pm, n_grid + 2)[1:-1]
153+ X, P = np.meshgrid(xs, ps)
154+ X, P = X.ravel(), P.ravel()
155+ rest = E - 0.5 * P**2 - 0.5 * sys.w1**2 * X**2
156+ ok = rest > 1e-9 * E
157+ X, P, rest = X[ok], P[ok], rest[ok]
158+ return np.array([X, np.zeros_like(X), P, np.sqrt(2 * rest)])
two-oscillators/plots.pyadded+102−0View file
@@ -0,0 +1,102 @@
1+"""Figures from results/*.npz (written by experiments.py). Run: python plots.py"""
2+
3+from __future__ import annotations
4+
5+import os
6+
7+import matplotlib
8+
9+matplotlib.use("Agg")
10+import matplotlib.pyplot as plt # noqa: E402
11+import numpy as np # noqa: E402
12+
13+HERE = os.path.dirname(os.path.abspath(__file__))
14+RES = os.path.join(HERE, "results")
15+CHAOS = 0.02 # finite-time Lyapunov exponent above which an orbit counts as chaotic
16+
17+
18+def sections():
19+ d = np.load(os.path.join(RES, "sections.npz"))
20+ energies, w2s = d["energies"], d["w2s"]
21+ fig, ax = plt.subplots(len(w2s), len(energies), figsize=(3 * len(energies), 3.1 * len(w2s)))
22+ for i, w2 in enumerate(w2s):
23+ for j, E in enumerate(energies):
24+ k = f"w2_{w2:g}_E_{E:g}"
25+ lam, owner, pts = d[k + "_lyap"], d[k + "_owner"], d[k + "_pts"]
26+ a = ax[i, j]
27+ chaotic = lam[owner] > CHAOS
28+ reg = ~chaotic
29+ a.scatter(pts[reg, 0], pts[reg, 1], s=0.15, c=owner[reg] % 10, cmap="tab10", vmin=0, vmax=9, rasterized=True)
30+ a.scatter(pts[chaotic, 0], pts[chaotic, 1], s=0.15, c="k", alpha=0.5, rasterized=True)
31+ a.set_title(f"w2 = {w2:g}, E = {E:g}\n{np.mean(lam > CHAOS):.0%} of orbits chaotic", fontsize=9)
32+ a.set_xticks([])
33+ a.set_yticks([])
34+ if j == 0:
35+ a.set_ylabel("p_x (section y = 0, p_y > 0)")
36+ if i == len(w2s) - 1:
37+ a.set_xlabel("x")
38+ fig.suptitle("Poincare sections: colored curves are regular orbits (tori), black dots chaotic orbits", y=1.0)
39+ fig.tight_layout()
40+ fig.savefig(os.path.join(RES, "sections.png"), dpi=110)
41+
42+
43+def lyapunov():
44+ d = np.load(os.path.join(RES, "lyapunov.npz"))
45+ energies, w2s = d["energies"], d["w2s"]
46+ fig, ax = plt.subplots(1, 2, figsize=(12, 4.5))
47+ lines = []
48+ for w2 in w2s:
49+ lams = [d[f"w2_{w2:g}_E{i}"] for i in range(len(energies))]
50+ frac = [np.mean(l > CHAOS) for l in lams]
51+ med = [np.median(l[l > CHAOS]) if np.any(l > CHAOS) else np.nan for l in lams]
52+ ax[0].semilogx(energies, frac, "o-", label=f"w2 = {w2:g}")
53+ ax[1].loglog(energies, med, "o-", label=f"w2 = {w2:g}")
54+ lines.append((w2, energies, frac, med))
55+ ax[0].set(xlabel="energy E (beta = 1)", ylabel="fraction of orbits chaotic", title="Onset of chaos")
56+ ax[1].set(xlabel="energy E", ylabel="median Lyapunov exponent of chaotic orbits",
57+ title="Rate of separation in the chaotic regime")
58+ for a in ax:
59+ a.grid(True, which="both", alpha=0.3)
60+ a.legend()
61+ fig.tight_layout()
62+ fig.savefig(os.path.join(RES, "lyapunov.png"), dpi=110)
63+ return lines
64+
65+
66+def damped():
67+ d = np.load(os.path.join(RES, "damped.npz"))
68+ E0s = sorted({float(k.split("_")[1]) for k in d.files if k.startswith("E0_")})
69+ fig, ax = plt.subplots(3, len(E0s), figsize=(5.5 * len(E0s), 10), sharex=True)
70+ summary = []
71+ for j, E0 in enumerate(E0s):
72+ k = f"E0_{E0:g}"
73+ t, Ex, Ey, E, logsum = d[k + "_t"], d[k + "_Ex"], d[k + "_Ey"], d[k + "_E"], d[k + "_logsum"]
74+ frac = Ex / (Ex + Ey)
75+ ax[0, j].semilogy(t, E[:, 0], "k")
76+ ax[0, j].set(title=f"strike energy E0 = {E0:g}", ylabel="energy (member 0)")
77+ ax[1, j].plot(t, frac, lw=0.6)
78+ ax[1, j].set(ylabel="share of energy in the x oscillator\n(16 strikes differing by 1e-10)", ylim=(0, 1))
79+ # local Lyapunov exponent: slope of the log-growth sum over windows of 20 time units
80+ w = max(1, int(round(20 / (t[1] - t[0]))))
81+ loc = (logsum[w:, 0] - logsum[:-w, 0]) / (t[w:] - t[:-w])
82+ ax[2, j].plot(t[w:], loc, "k", lw=0.8)
83+ ax[2, j].axhline(0, color="0.6", lw=0.5)
84+ ax[2, j].set(xlabel="time", ylabel="local Lyapunov exponent (20-unit windows)")
85+ late = t > t[-1] * 0.8
86+ summary.append((E0, float(frac[late].mean(0).std()), float(frac[late].mean(0).min()), float(frac[late].mean(0).max())))
87+ fig.tight_layout()
88+ fig.savefig(os.path.join(RES, "damped.png"), dpi=110)
89+ return summary
90+
91+
92+if __name__ == "__main__":
93+ for name, f in (("sections", sections), ("lyapunov", lyapunov), ("damped", damped)):
94+ if os.path.exists(os.path.join(RES, name + ".npz")):
95+ out = f()
96+ if name == "lyapunov":
97+ for w2, energies, frac, med in out:
98+ print(f"w2 = {w2:g}:", " ".join(f"E {E:.2g}: {fr:.2f}" for E, fr in zip(energies, frac)))
99+ if name == "damped":
100+ for E0, sd, lo, hi in out:
101+ print(f"damped E0 {E0:g}: late share of energy in x across strikes: std {sd:.3f}, range {lo:.3f} to {hi:.3f}")
102+ print("wrote", os.path.join(RES, name + ".png"))
two-oscillators/results/damped.pngadded+0−0View file
Binary file not shown.
two-oscillators/results/lyapunov.pngadded+0−0View file
Binary file not shown.
two-oscillators/results/sections.pngadded+0−0View file
Binary file not shown.
web/index.htmladded+1090−0View file
@@ -0,0 +1,1090 @@
1+<meta charset="utf-8">
2+<title>Coupled Oscillator Chaos</title>
3+<link rel="preconnect" href="https://fonts.googleapis.com">
4+<link rel="preconnect" href="https://fonts.gstatic.com" crossorigin>
5+<link rel="stylesheet" href="https://fonts.googleapis.com/css2?family=Archivo:wdth,wght@75..100,400..800&family=IBM+Plex+Mono:wght@400;500&display=swap">
6+<style>
7+:root {
8+ --paper: #eef2f5;
9+ --panel: #f8fafb;
10+ --ink: #172334;
11+ --muted: #5b6b7d;
12+ --rule: #cbd5df;
13+ --grid: #dde4eb;
14+ --x: #2563a6;
15+ --y: #a23b72;
16+ --chaos: #e8b931;
17+ --chaos-wash: rgba(232, 185, 49, 0.18);
18+ --faint: rgba(23, 35, 52, 0.16);
19+ --focus: #2563a6;
20+ --sans: "Archivo", "Helvetica Neue", Arial, sans-serif;
21+ --mono: "IBM Plex Mono", ui-monospace, "SFMono-Regular", Menlo, monospace;
22+}
23+@media (prefers-color-scheme: dark) {
24+ :root:not([data-theme="light"]) {
25+ color-scheme: dark;
26+ --paper: #0f161f;
27+ --panel: #16202c;
28+ --ink: #e3e9f0;
29+ --muted: #93a3b5;
30+ --rule: #2a3747;
31+ --grid: #1f2a37;
32+ --x: #6fa8e8;
33+ --y: #d57bae;
34+ --chaos: #e8b931;
35+ --chaos-wash: rgba(232, 185, 49, 0.14);
36+ --faint: rgba(227, 233, 240, 0.16);
37+ --focus: #6fa8e8;
38+ }
39+}
40+:root[data-theme="dark"] {
41+ color-scheme: dark;
42+ --paper: #0f161f;
43+ --panel: #16202c;
44+ --ink: #e3e9f0;
45+ --muted: #93a3b5;
46+ --rule: #2a3747;
47+ --grid: #1f2a37;
48+ --x: #6fa8e8;
49+ --y: #d57bae;
50+ --chaos: #e8b931;
51+ --chaos-wash: rgba(232, 185, 49, 0.14);
52+ --faint: rgba(227, 233, 240, 0.16);
53+ --focus: #6fa8e8;
54+}
55+* { box-sizing: border-box; }
56+body {
57+ background: var(--paper);
58+ color: var(--ink);
59+ font-family: var(--sans);
60+ font-size: 15px;
61+ line-height: 1.5;
62+ padding-inline: 16px;
63+ padding-block: 20px 48px;
64+}
65+.wrap { max-width: 1280px; margin: 0 auto; display: grid; gap: 20px; }
66+header { display: grid; gap: 6px; }
67+h1 {
68+ font-size: clamp(26px, 4vw, 38px);
69+ font-weight: 750;
70+ font-stretch: 78%;
71+ letter-spacing: -0.01em;
72+ line-height: 1.05;
73+ margin: 0;
74+ text-wrap: balance;
75+}
76+h1 .eq { font-family: var(--mono); font-weight: 400; font-size: 0.42em; color: var(--muted); letter-spacing: 0; display: block; margin-top: 8px; font-stretch: 100%; }
77+.lede { max-width: 72ch; color: var(--muted); margin: 0; }
78+.layout { display: grid; grid-template-columns: 300px minmax(0, 1fr); gap: 20px; align-items: start; }
79+@media (max-width: 900px) { .layout { grid-template-columns: 1fr; } }
80+aside {
81+ background: var(--panel);
82+ border: 1px solid var(--rule);
83+ border-radius: 6px;
84+ padding: 16px;
85+ display: grid;
86+ gap: 16px;
87+ position: sticky;
88+ top: calc(env(safe-area-inset-top, 0px) + 12px);
89+}
90+@media (max-width: 900px) { aside { position: static; } }
91+.label {
92+ font-size: 11px;
93+ font-weight: 650;
94+ letter-spacing: 0.08em;
95+ text-transform: uppercase;
96+ color: var(--muted);
97+}
98+.presets { display: flex; flex-wrap: wrap; gap: 6px; }
99+.chip {
100+ font: inherit;
101+ font-size: 13px;
102+ padding: 4px 10px;
103+ border-radius: 999px;
104+ border: 1px solid var(--rule);
105+ background: var(--paper);
106+ color: var(--ink);
107+ cursor: pointer;
108+}
109+.chip:hover { border-color: var(--muted); }
110+.field { display: grid; gap: 4px; }
111+.field .row { display: flex; justify-content: space-between; align-items: baseline; gap: 8px; }
112+.field label { font-size: 14px; font-weight: 550; }
113+.val { font-family: var(--mono); font-size: 13px; font-variant-numeric: tabular-nums; color: var(--ink); white-space: nowrap; }
114+.hint { font-size: 12px; color: var(--muted); margin: 0; }
115+input[type="range"] { width: 100%; accent-color: var(--x); }
116+.player { display: grid; gap: 8px; }
117+.btns { display: grid; grid-template-columns: 1fr 1fr; gap: 8px; }
118+.btn {
119+ font: inherit;
120+ font-weight: 600;
121+ font-size: 14px;
122+ padding: 9px 10px;
123+ border-radius: 6px;
124+ border: 1px solid var(--rule);
125+ background: var(--paper);
126+ color: var(--ink);
127+ cursor: pointer;
128+ display: flex;
129+ align-items: center;
130+ justify-content: center;
131+ gap: 8px;
132+}
133+.btn:hover { border-color: var(--muted); }
134+.btn .sw { width: 10px; height: 10px; border-radius: 2px; flex: none; }
135+.btn.wide { grid-column: 1 / -1; }
136+.btn[aria-pressed="true"] { background: var(--ink); color: var(--paper); border-color: var(--ink); }
137+button:focus-visible, input:focus-visible { outline: 2px solid var(--focus); outline-offset: 2px; }
138+.stats { display: grid; grid-template-columns: auto 1fr; gap: 4px 12px; font-size: 13px; margin: 0; }
139+.stats dt { color: var(--muted); }
140+.stats dd { margin: 0; font-family: var(--mono); font-variant-numeric: tabular-nums; text-align: right; }
141+.status { font-size: 12px; color: var(--muted); min-height: 1.5em; font-family: var(--mono); }
142+main { display: grid; gap: 16px; min-width: 0; }
143+.pair { display: grid; grid-template-columns: 1fr 1fr; gap: 16px; }
144+@media (max-width: 640px) { .pair { grid-template-columns: 1fr; } }
145+figure { margin: 0; display: grid; gap: 6px; min-width: 0; }
146+figcaption { display: grid; gap: 2px; }
147+figcaption b { font-weight: 650; font-size: 15px; }
148+figcaption span { font-size: 13px; color: var(--muted); max-width: 80ch; }
149+.plot { position: relative; width: 100%; background: var(--panel); border: 1px solid var(--rule); border-radius: 6px; }
150+.plot canvas { position: absolute; inset: 0; width: 100%; height: 100%; display: block; }
151+.sq { aspect-ratio: 1 / 1; max-width: 100%; }
152+.wide { aspect-ratio: 3.2 / 1; min-height: 170px; }
153+.short { aspect-ratio: 4.6 / 1; min-height: 140px; }
154+.setup { display: grid; grid-template-columns: minmax(0, 1.25fr) minmax(0, 1fr); gap: 16px; align-items: start; }
155+@media (max-width: 760px) { .setup { grid-template-columns: 1fr; } }
156+.setupplot { aspect-ratio: 1.45 / 1; }
157+.setuptext { display: grid; gap: 10px; font-size: 14px; }
158+.setuptext p { margin: 0; max-width: 60ch; }
159+.setuptext .k { font-family: var(--mono); font-size: 13px; }
160+.setuptext .btn { justify-self: start; }
161+.strikes { display: flex; flex-wrap: wrap; gap: 6px; align-items: center; }
162+.sk { font: inherit; font-family: var(--mono); font-size: 13px; min-width: 38px; padding: 4px 8px; border-radius: 6px; border: 1px solid var(--rule); background: var(--paper); color: var(--ink); cursor: pointer; display: inline-flex; align-items: center; gap: 6px; justify-content: center; }
163+.sk .sw { width: 8px; height: 8px; border-radius: 2px; }
164+.sk[aria-pressed="true"] { background: var(--ink); color: var(--paper); border-color: var(--ink); }
165+.sk:disabled { opacity: 0.4; cursor: default; }
166+.strikes .btn { padding: 4px 12px; font-size: 13px; margin-left: auto; }
167+.legend { display: flex; flex-wrap: wrap; gap: 14px; font-size: 12px; color: var(--muted); }
168+.legend i { display: inline-block; width: 18px; height: 3px; vertical-align: middle; margin-right: 6px; border-radius: 2px; }
169+footer { color: var(--muted); font-size: 13px; max-width: 80ch; }
170+footer p { margin: 0 0 8px; }
171+@media (prefers-reduced-motion: reduce) { * { scroll-behavior: auto !important; } }
172+</style>
173+
174+<div class="wrap">
175+ <header>
176+ <h1>Coupled Oscillator Chaos
177+ <span class="eq">H = ½(p<sub>x</sub>² + p<sub>y</sub>²) + ½(x² + ω²y²) + x²y²</span>
178+ </h1>
179+ <p class="lede">One mass on springs that can swing two ways, x and y, which stiffen each other as they move: two oscillators with a quartic coupling, struck once and left to ring down. A gentle strike rings regularly; a strong one goes chaotic until damping takes the energy back below the threshold, and then it rings on with whatever split of energy the chaos left behind. Strikes A and B differ by a tiny change in direction. Listen to both.</p>
180+ </header>
181+
182+ <div class="layout">
183+ <aside aria-label="Controls">
184+ <div class="player">
185+ <div class="label">Listen</div>
186+ <div class="btns">
187+ <button class="btn" id="playA" type="button"><span class="sw" style="background: var(--x)"></span>Strike A</button>
188+ <button class="btn" id="playB" type="button"><span class="sw" style="background: var(--y)"></span>Strike B</button>
189+ <button class="btn" id="playAB" type="button">A, then B</button>
190+ <button class="btn" id="stop" type="button">Stop</button>
191+ </div>
192+ <div class="status" id="status">Computing…</div>
193+ </div>
194+
195+ <div class="field">
196+ <div class="label">Presets</div>
197+ <div class="presets">
198+ <button class="chip" type="button" data-preset="gentle">Gentle</button>
199+ <button class="chip" type="button" data-preset="edge">Edge of chaos</button>
200+ <button class="chip" type="button" data-preset="strong">Strong</button>
201+ <button class="chip" type="button" data-preset="equal">Equal pitches</button>
202+ </div>
203+ </div>
204+
205+ <div class="field">
206+ <div class="row"><label for="E0">Strike energy E</label><span class="val" id="E0v"></span></div>
207+ <input id="E0" type="range" min="0" max="1000" step="1">
208+ <p class="hint">The only parameter of the dynamics (β = 1). Chaos sets in around E ≈ 2 for ω = 1.618.</p>
209+ </div>
210+ <div class="field">
211+ <div class="row"><label for="w2">Frequency ratio ω</label><span class="val" id="w2v"></span></div>
212+ <input id="w2" type="range" min="1" max="2.5" step="0.001">
213+ </div>
214+ <div class="field">
215+ <div class="row"><label for="ang">Strike direction</label><span class="val" id="angv"></span></div>
216+ <input id="ang" type="range" min="5" max="85" step="0.5">
217+ <p class="hint">How the strike's momentum is split between x and y.</p>
218+ </div>
219+ <div class="field">
220+ <div class="row"><label for="dl">Difference between A and B</label><span class="val" id="dlv"></span></div>
221+ <input id="dl" type="range" min="-14" max="-2" step="0.5">
222+ </div>
223+ <div class="field">
224+ <div class="row"><label for="t60">Ring time (−60 dB)</label><span class="val" id="t60v"></span></div>
225+ <input id="t60" type="range" min="0.5" max="10" step="0.1">
226+ </div>
227+ <div class="field">
228+ <div class="row"><label for="f1">Pitch of x</label><span class="val" id="f1v"></span></div>
229+ <input id="f1" type="range" min="80" max="880" step="1">
230+ </div>
231+ <div class="field">
232+ <div class="row"><label for="dur">Length</label><span class="val" id="durv"></span></div>
233+ <input id="dur" type="range" min="1" max="8" step="0.5">
234+ </div>
235+
236+ <dl class="stats" aria-live="polite">
237+ <dt>Chaotic phase ends</dt><dd id="sChaos">–</dd>
238+ <dt>Peak separation rate</dt><dd id="sLam">–</dd>
239+ <dt>Final share in x, A</dt><dd id="sA">–</dd>
240+ <dt>Final share in x, B</dt><dd id="sB">–</dd>
241+ <dt>Spread over 12 strikes</dt><dd id="sSpread">–</dd>
242+ </dl>
243+ </aside>
244+
245+ <main>
246+ <figure>
247+ <figcaption><b>Physical setup</b><span>One mass, free to slide in a plane, held by springs whose anchors roll along the walls, so each spring pushes only along its own axis.</span></figcaption>
248+ <div class="setup">
249+ <div class="plot setupplot"><canvas id="cSetup"></canvas></div>
250+ <div class="setuptext">
251+ <p>The mass (m = 1) sits at rest in the middle until the strike: a sharp blow that gives it energy E in the chosen direction. It then oscillates along x at frequency 1 and along y at frequency ω.</p>
252+ <p>What couples the two motions is the stiffness. The springs along x get stiffer the farther the mass is from the x axis, and the springs along y get stiffer the farther it is from the y axis:</p>
253+ <p class="k">k<sub>x</sub> = 1 + 2y² k<sub>y</sub> = ω² + 2x²</p>
254+ <p>That is the x²y² term in the energy. A plate behaves the same way: bending it one way stretches its midplane, which stiffens it for every other way of bending. A light friction on the mass sets the ring time.</p>
255+ <button class="btn" id="slow" type="button" aria-pressed="false">Slow motion</button>
256+ </div>
257+ </div>
258+ </figure>
259+
260+ <div class="pair">
261+ <figure>
262+ <figcaption><b>Motion in the bowl</b><span>The same mass seen as a ball rolling in a bowl, V(x, y) = ½x² + ½ω²y² + x²y². Bands are heights of 10%, 25%, 50% and 100% of the strike energy; the ball cannot climb past the outer one. The line is the path of the selected strike (A unless you choose another below).</span></figcaption>
263+ <div class="plot sq"><canvas id="cPath"></canvas><canvas id="cPathO"></canvas></div>
264+ </figure>
265+ <figure>
266+ <figcaption><b>Phase space at the strike energy</b><span>Poincaré section of the undamped system (y = 0, p<sub>y</sub> > 0). Closed curves are regular motion, scattered dots chaos. The ring marks where the strike starts.</span></figcaption>
267+ <div class="plot sq"><canvas id="cSec"></canvas></div>
268+ </figure>
269+ </div>
270+
271+ <figure>
272+ <figcaption><b>Who gets the energy</b><span>Share of the energy in the x oscillator for 12 strikes, each differing from the next by the A/B difference. The shaded span is the chaotic phase of strike A.</span></figcaption>
273+ <div class="plot wide"><canvas id="cShare"></canvas><canvas id="cShareO"></canvas></div>
274+ <div class="legend"><span><i style="background: var(--x)"></i>Strike A</span><span><i style="background: var(--y)"></i>Strike B</span><span><i style="background: var(--faint)"></i>10 more strikes</span></div>
275+ </figure>
276+
277+ <figure>
278+ <figcaption><b>How fast nearby strikes separate</b><span>Local Lyapunov exponent of strike A (per second, 20-period windows). Above zero, a difference between two strikes grows exponentially. The dashed line is the energy, in dB.</span></figcaption>
279+ <div class="plot short"><canvas id="cLam"></canvas><canvas id="cLamO"></canvas></div>
280+ </figure>
281+
282+ <figure>
283+ <figcaption><b>What you hear</b><span>Spectrogram of the selected strike (pickup p<sub>x</sub> + p<sub>y</sub>). Broadband smear means chaos; clean lines mean regular ringing. Click a strike, or use the ← and → keys, to flip between the 12; the color scale is shared, so differences in level are real. The selected strike is also drawn bold in the energy-share plot and as the path in the bowl.</span></figcaption>
284+ <div class="strikes" role="group" aria-label="Choose a strike" id="strikeSel"></div>
285+ <div class="plot wide"><canvas id="cSpec"></canvas><canvas id="cSpecO"></canvas></div>
286+ </figure>
287+
288+ <footer>
289+ <p>Integrated in your browser with a fourth-order symplectic scheme (Yoshida), damping applied to the momenta. Time is scaled so that ω<sub>x</sub> = 1 sounds at the chosen pitch. The model is the toy behind the struck-plate study in cymbal-elastodynamics: the plate's modes are oscillators coupled through the same kind of quartic stretching energy.</p>
290+ </footer>
291+ </main>
292+ </div>
293+</div>
294+
295+<script id="worker-src" type="text/plain">
296+const C = Math.cbrt(2), Y1 = 1 / (2 - C), Y0 = -C / (2 - C);
297+const COEF = [Y1, Y0, Y1];
298+
299+function hmaxFor(E) { return Math.min(0.02, 0.1 / Math.sqrt(1 + 4 * E)); }
300+
301+// one strike: state s = [x, y, px, py]; optional tangent d for the Lyapunov exponent
302+function strike(p, angle, withTangent, audioRate) {
303+ const w2sq = p.w2 * p.w2;
304+ const P = Math.sqrt(2 * p.E0);
305+ let x = 0, y = 0, px = P * Math.cos(angle), py = P * Math.sin(angle);
306+ let dx = 0.6, dy = 0.3, dpx = -0.5, dpy = 0.55;
307+ { const n = Math.hypot(dx, dy, dpx, dpy); dx /= n; dy /= n; dpx /= n; dpy /= n; }
308+ const nOut = p.nOut, dtOut = p.dtOut;
309+ const nsub = Math.max(1, Math.ceil(dtOut / hmaxFor(p.E0)));
310+ const h = dtOut / nsub;
311+ const damp = Math.exp(-p.gamma * h);
312+ const K = p.decim;
313+ const nD = Math.floor(nOut / K);
314+ const sig = audioRate ? new Float32Array(nOut) : null;
315+ const xs = audioRate ? new Float32Array(nOut) : null;
316+ const ys = audioRate ? new Float32Array(nOut) : null;
317+ const share = new Float32Array(nD);
318+ const energy = new Float32Array(nD);
319+ const logsum = withTangent ? new Float64Array(nD) : null;
320+ let L = 0;
321+ for (let i = 0; i < nOut; i++) {
322+ for (let s = 0; s < nsub; s++) {
323+ for (let c = 0; c < 3; c++) {
324+ const hh = COEF[c] * h;
325+ let fx = -x - 2 * x * y * y, fy = -w2sq * y - 2 * x * x * y;
326+ px += 0.5 * hh * fx; py += 0.5 * hh * fy;
327+ if (withTangent) {
328+ const fxx = -1 - 2 * y * y, fyy = -w2sq - 2 * x * x, fxy = -4 * x * y;
329+ dpx += 0.5 * hh * (fxx * dx + fxy * dy); dpy += 0.5 * hh * (fxy * dx + fyy * dy);
330+ dx += hh * dpx; dy += hh * dpy;
331+ }
332+ x += hh * px; y += hh * py;
333+ fx = -x - 2 * x * y * y; fy = -w2sq * y - 2 * x * x * y;
334+ px += 0.5 * hh * fx; py += 0.5 * hh * fy;
335+ if (withTangent) {
336+ const fxx = -1 - 2 * y * y, fyy = -w2sq - 2 * x * x, fxy = -4 * x * y;
337+ dpx += 0.5 * hh * (fxx * dx + fxy * dy); dpy += 0.5 * hh * (fxy * dx + fyy * dy);
338+ }
339+ }
340+ px *= damp; py *= damp;
341+ if (withTangent) { dpx *= damp; dpy *= damp; }
342+ }
343+ if (withTangent) {
344+ const n = Math.hypot(dx, dy, dpx, dpy);
345+ L += Math.log(n); dx /= n; dy /= n; dpx /= n; dpy /= n;
346+ }
347+ if (audioRate) { sig[i] = px + py; xs[i] = x; ys[i] = y; }
348+ if ((i + 1) % K === 0) {
349+ const j = (i + 1) / K - 1;
350+ if (j < nD) {
351+ const ex = 0.5 * (px * px + x * x), ey = 0.5 * (py * py + w2sq * y * y);
352+ share[j] = ex / (ex + ey + 1e-300);
353+ energy[j] = ex + ey + x * x * y * y;
354+ if (withTangent) logsum[j] = L;
355+ }
356+ }
357+ }
358+ return { sig, xs, ys, share, energy, logsum };
359+}
360+
361+function section(p) {
362+ // undamped orbits at energy E0 started on the section, plus the strike orbit
363+ const w2sq = p.w2 * p.w2, E = p.E0;
364+ const h = hmaxFor(E), T = 1200;
365+ const n = Math.round(T / h);
366+ const ics = [];
367+ const xm = Math.sqrt(2 * E), pm = Math.sqrt(2 * E);
368+ const g = 7;
369+ for (let a = 1; a <= g; a++) for (let b = 1; b <= g; b++) {
370+ const X = -xm + 2 * xm * a / (g + 1), Pp = -pm + 2 * pm * b / (g + 1);
371+ const rest = E - 0.5 * Pp * Pp - 0.5 * X * X;
372+ if (rest > 1e-6 * E) ics.push([X, 0, Pp, Math.sqrt(2 * rest)]);
373+ }
374+ const P = Math.sqrt(2 * E);
375+ ics.push([0, 0, P * Math.cos(p.angle), P * Math.sin(p.angle)]);
376+ const out = [];
377+ for (const ic of ics) {
378+ let [x, y, px, py] = ic;
379+ let dx = 0.6, dy = 0.3, dpx = -0.5, dpy = 0.55, L = 0;
380+ { const nn = Math.hypot(dx, dy, dpx, dpy); dx /= nn; dy /= nn; dpx /= nn; dpy /= nn; }
381+ const pts = [];
382+ for (let i = 0; i < n; i++) {
383+ const yOld = y, xOld = x, pxOld = px;
384+ for (let c = 0; c < 3; c++) {
385+ const hh = COEF[c] * h;
386+ let fx = -x - 2 * x * y * y, fy = -w2sq * y - 2 * x * x * y;
387+ let fxx = -1 - 2 * y * y, fyy = -w2sq - 2 * x * x, fxy = -4 * x * y;
388+ px += 0.5 * hh * fx; py += 0.5 * hh * fy;
389+ dpx += 0.5 * hh * (fxx * dx + fxy * dy); dpy += 0.5 * hh * (fxy * dx + fyy * dy);
390+ dx += hh * dpx; dy += hh * dpy;
391+ x += hh * px; y += hh * py;
392+ fx = -x - 2 * x * y * y; fy = -w2sq * y - 2 * x * x * y;
393+ fxx = -1 - 2 * y * y; fyy = -w2sq - 2 * x * x; fxy = -4 * x * y;
394+ px += 0.5 * hh * fx; py += 0.5 * hh * fy;
395+ dpx += 0.5 * hh * (fxx * dx + fxy * dy); dpy += 0.5 * hh * (fxy * dx + fyy * dy);
396+ }
397+ if (i % 10 === 9) { const nn = Math.hypot(dx, dy, dpx, dpy); L += Math.log(nn); dx /= nn; dy /= nn; dpx /= nn; dpy /= nn; }
398+ if (yOld < 0 && y >= 0 && py > 0) {
399+ const f = -yOld / (y - yOld);
400+ pts.push((1 - f) * xOld + f * x, (1 - f) * pxOld + f * px);
401+ }
402+ }
403+ out.push({ pts: new Float32Array(pts), lam: L / T });
404+ }
405+ return out;
406+}
407+
408+let gen = 0;
409+onmessage = (ev) => {
410+ const m = ev.data;
411+ gen = m.gen;
412+ const p = m.p;
413+ const a = strike(p, p.angle, true, true);
414+ const b = strike(p, p.angle + p.delta, true, true);
415+ postMessage({ kind: "strikes", gen: m.gen, a, b },
416+ [a.sig.buffer, a.xs.buffer, a.ys.buffer, a.share.buffer, a.energy.buffer, a.logsum.buffer,
417+ b.sig.buffer, b.xs.buffer, b.ys.buffer, b.share.buffer, b.energy.buffer, b.logsum.buffer]);
418+ for (let k = 2; k < 12; k++) {
419+ const r = strike(p, p.angle + k * p.delta, false, true);
420+ postMessage({ kind: "member", gen: m.gen, k, share: r.share, sig: r.sig, xs: r.xs, ys: r.ys },
421+ [r.share.buffer, r.sig.buffer, r.xs.buffer, r.ys.buffer]);
422+ }
423+ postMessage({ kind: "section", gen: m.gen, sec: section(p) });
424+};
425+</script>
426+
427+<script>
428+(() => {
429+ const SR = 44100;
430+ const $ = (id) => document.getElementById(id);
431+ const css = (name) => getComputedStyle(document.documentElement).getPropertyValue(name).trim();
432+
433+ const PRESETS = {
434+ gentle: { E0: 0.5, w2: 1.618, ang: 40, dl: -10, t60: 4, f1: 220, dur: 4 },
435+ edge: { E0: 3, w2: 1.618, ang: 40, dl: -10, t60: 4, f1: 220, dur: 4 },
436+ strong: { E0: 20, w2: 1.618, ang: 40, dl: -10, t60: 4, f1: 220, dur: 4 },
437+ equal: { E0: 4, w2: 1.0, ang: 40, dl: -10, t60: 4, f1: 220, dur: 4 },
438+ };
439+
440+ // slider <-> value maps
441+ const eFromSlider = (v) => Math.pow(10, -1 + 3 * v / 1000);
442+ const sliderFromE = (E) => Math.round((Math.log10(E) + 1) / 3 * 1000);
443+ const fmt = (v, d = 2) => Number(v).toFixed(d);
444+
445+ function readParams() {
446+ const E0 = eFromSlider(+$("E0").value);
447+ const w2 = +$("w2").value, angDeg = +$("ang").value, dl = +$("dl").value;
448+ const t60 = +$("t60").value, f1 = +$("f1").value, dur = +$("dur").value;
449+ $("E0v").textContent = E0 < 10 ? fmt(E0, 2) : fmt(E0, 1);
450+ $("w2v").textContent = fmt(w2, 3) + (Math.abs(w2 - 1.618) < 0.002 ? " (golden)" : "");
451+ $("angv").textContent = fmt(angDeg, 1) + "°";
452+ $("dlv").textContent = "1e" + dl + " rad";
453+ $("t60v").textContent = fmt(t60, 1) + " s";
454+ $("f1v").textContent = Math.round(f1) + " Hz (y: " + Math.round(f1 * w2) + " Hz)";
455+ $("durv").textContent = fmt(dur, 1) + " s";
456+ const tau = 1 / (2 * Math.PI * f1); // seconds per time unit
457+ const nOut = Math.round(dur * SR);
458+ return {
459+ E0, w2, angle: angDeg * Math.PI / 180, delta: Math.pow(10, dl), tau,
460+ // momenta damped by exp(-gamma dt): energy decays at rate gamma, amplitude at gamma / 2
461+ gamma: 2 * Math.log(1000) / (t60 / tau), dtOut: 1 / SR / tau, nOut, decim: Math.round(SR / 400), dur,
462+ };
463+ }
464+
465+ function setControls(p) {
466+ $("E0").value = sliderFromE(p.E0); $("w2").value = p.w2; $("ang").value = p.ang; $("dl").value = p.dl;
467+ $("t60").value = p.t60; $("f1").value = p.f1; $("dur").value = p.dur;
468+ }
469+
470+ // ---------------- worker ----------------
471+ const workerURL = URL.createObjectURL(new Blob([$("worker-src").textContent], { type: "text/javascript" }));
472+ let worker = null, gen = 0, params = null;
473+ let A = null, B = null, members = [], sec = null;
474+
475+ function run() {
476+ stopAudio();
477+ if (smOn) { smOn = false; cancelAnimationFrame(smRaf); $("slow").setAttribute("aria-pressed", "false"); $("slow").textContent = "Slow motion"; }
478+ params = readParams();
479+ gen++;
480+ if (worker) worker.terminate();
481+ worker = new Worker(workerURL);
482+ worker.onmessage = onResult;
483+ A = B = null; members = []; sec = null;
484+ $("status").textContent = "Computing…";
485+ worker.postMessage({ gen, p: params });
486+ drawAll();
487+ }
488+
489+ function onResult(ev) {
490+ const m = ev.data;
491+ if (m.gen !== gen) return;
492+ if (m.kind === "strikes") {
493+ A = m.a; B = m.b;
494+ members = [];
495+ prepareAudio();
496+ updateStats();
497+ $("status").textContent = "Ready. Strikes differ by 1e" + $("dl").value + " rad in direction.";
498+ } else if (m.kind === "member") {
499+ members[m.k] = { share: m.share, sig: m.sig, xs: m.xs, ys: m.ys, audio: scaleAudio(m.sig) };
500+ updateStats();
501+ } else if (m.kind === "section") {
502+ sec = m.sec;
503+ }
504+ renderSelector();
505+ drawAll();
506+ }
507+
508+ // ---------------- strike selector ----------------
509+ const NS = 12;
510+ let sel = 0;
511+ const strikeAt = (k) => (k === 0 ? A : k === 1 ? B : members[k] || null);
512+ const strikeColor = (k) => (k === 0 ? css("--x") : k === 1 ? css("--y") : css("--ink"));
513+ const strikeName = (k) => (k === 0 ? "A" : k === 1 ? "B" : String(k + 1));
514+ function renderSelector() {
515+ const box = $("strikeSel");
516+ if (!box.childElementCount) {
517+ for (let k = 0; k < NS; k++) {
518+ const b = document.createElement("button");
519+ b.type = "button"; b.className = "sk"; b.id = "sk" + k;
520+ b.innerHTML = '<span class="sw"></span>' + strikeName(k);
521+ b.addEventListener("click", () => selectStrike(k));
522+ box.appendChild(b);
523+ }
524+ const pb = document.createElement("button");
525+ pb.type = "button"; pb.className = "btn"; pb.id = "playSel"; pb.textContent = "Play this strike";
526+ pb.addEventListener("click", () => play(sel, null));
527+ box.appendChild(pb);
528+ }
529+ for (let k = 0; k < NS; k++) {
530+ const b = $("sk" + k);
531+ b.disabled = !strikeAt(k);
532+ b.setAttribute("aria-pressed", String(k === sel));
533+ b.querySelector(".sw").style.background = strikeColor(k);
534+ b.title = strikeAt(k) ? "Strike " + strikeName(k) + ": direction + " + k + " × 1e" + $("dl").value + " rad" : "computing…";
535+ }
536+ $("playSel").disabled = !strikeAt(sel);
537+ }
538+ function selectStrike(k) {
539+ if (!strikeAt(k)) return;
540+ sel = k;
541+ renderSelector();
542+ drawPath(); drawShare(); drawSpec();
543+ }
544+ document.addEventListener("keydown", (e) => {
545+ if (e.target && (e.target.tagName === "INPUT" || e.target.tagName === "TEXTAREA")) return;
546+ if (e.key !== "ArrowRight" && e.key !== "ArrowLeft") return;
547+ const dir = e.key === "ArrowRight" ? 1 : -1;
548+ for (let i = 1; i <= NS; i++) {
549+ const k = (sel + dir * i + NS) % NS;
550+ if (strikeAt(k)) { selectStrike(k); e.preventDefault(); return; }
551+ }
552+ });
553+
554+ // ---------------- derived quantities ----------------
555+ function localLyap(logsum) {
556+ // per second, over windows of 20 periods of x (20 * 2 pi time units)
557+ const dtD = params.decim * params.dtOut; // time units per decimated sample
558+ const w = Math.max(1, Math.round(20 * 2 * Math.PI / dtD));
559+ const out = new Float32Array(logsum.length);
560+ for (let i = 0; i < logsum.length; i++) {
561+ const j = Math.max(0, i - w);
562+ out[i] = i > j ? (logsum[i] - logsum[j]) / ((i - j) * dtD) / params.tau : NaN;
563+ }
564+ return out;
565+ }
566+
567+ function chaosEnd(lam) {
568+ // last time the local exponent is clearly positive (0.02 per time unit)
569+ const thr = 0.02 / params.tau;
570+ let last = -1;
571+ for (let i = 0; i < lam.length; i++) if (lam[i] > thr) last = i;
572+ return last < 0 ? null : (last + 1) * params.decim / SR;
573+ }
574+
575+ function updateStats() {
576+ if (!A) return;
577+ const lam = localLyap(A.logsum);
578+ A.lam = lam;
579+ const end = chaosEnd(lam);
580+ $("sChaos").textContent = end === null ? "never chaotic" : fmt(end, 2) + " s";
581+ let peak = 0;
582+ for (const v of lam) if (v > peak) peak = v;
583+ $("sLam").textContent = peak > 0.02 / params.tau ? fmt(peak, 0) + " /s" : "≈ 0";
584+ const fin = (s) => { let t = 0, n = 0; for (let i = Math.floor(s.length * 0.85); i < s.length; i++) { t += s[i]; n++; } return t / n; };
585+ const fa = fin(A.share), fb = fin(B.share);
586+ $("sA").textContent = fmt(fa, 3);
587+ $("sB").textContent = fmt(fb, 3);
588+ const all = [fa, fb];
589+ for (let k = 2; k < 12; k++) if (members[k]) all.push(fin(members[k].share));
590+ const mean = all.reduce((s, v) => s + v, 0) / all.length;
591+ const sd = Math.sqrt(all.reduce((s, v) => s + (v - mean) ** 2, 0) / all.length);
592+ let lo = Math.min(...all), hi = Math.max(...all);
593+ $("sSpread").textContent = all.length < 12 ? "computing…" : fmt(lo, 2) + " to " + fmt(hi, 2);
594+ }
595+
596+ // ---------------- canvas helpers ----------------
597+ function setup(canvas) {
598+ const r = canvas.getBoundingClientRect();
599+ const dpr = window.devicePixelRatio || 1;
600+ const w = Math.max(1, Math.round(r.width * dpr)), h = Math.max(1, Math.round(r.height * dpr));
601+ if (canvas.width !== w || canvas.height !== h) { canvas.width = w; canvas.height = h; }
602+ const ctx = canvas.getContext("2d");
603+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
604+ ctx.clearRect(0, 0, r.width, r.height);
605+ return { ctx, W: r.width, H: r.height };
606+ }
607+
608+ const PAD = { l: 44, r: 12, t: 10, b: 26 };
609+
610+ function axesTime(ctx, W, H, ymin, ymax, yticks, ylabelFmt, dur) {
611+ const x0 = PAD.l, x1 = W - PAD.r, y0 = PAD.t, y1 = H - PAD.b;
612+ ctx.strokeStyle = css("--grid"); ctx.lineWidth = 1;
613+ ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono");
614+ ctx.textAlign = "right"; ctx.textBaseline = "middle";
615+ for (const v of yticks) {
616+ const y = y1 - (v - ymin) / (ymax - ymin) * (y1 - y0);
617+ ctx.beginPath(); ctx.moveTo(x0, y); ctx.lineTo(x1, y); ctx.stroke();
618+ ctx.fillText(ylabelFmt(v), x0 - 6, y);
619+ }
620+ ctx.textAlign = "center"; ctx.textBaseline = "top";
621+ const step = dur <= 2 ? 0.25 : dur <= 5 ? 0.5 : 1;
622+ for (let t = 0; t <= dur + 1e-9; t += step) {
623+ const x = x0 + t / dur * (x1 - x0);
624+ ctx.beginPath(); ctx.moveTo(x, y0); ctx.lineTo(x, y1); ctx.stroke();
625+ ctx.fillText(fmt(t, step < 1 ? (step < 0.5 ? 2 : 1) : 0) + (t === 0 ? " s" : ""), x, y1 + 6);
626+ }
627+ return { x0, x1, y0, y1, X: (t) => x0 + t / dur * (x1 - x0), Y: (v) => y1 - (v - ymin) / (ymax - ymin) * (y1 - y0) };
628+ }
629+
630+ function chaosSpan(ctx, ax) {
631+ if (!A || !A.lam) return;
632+ const end = chaosEnd(A.lam);
633+ if (end === null) return;
634+ ctx.fillStyle = css("--chaos-wash");
635+ ctx.fillRect(ax.x0, ax.y0, ax.X(end) - ax.x0, ax.y1 - ax.y0);
636+ }
637+
638+ function seriesTime(ctx, ax, arr, color, width, dash) {
639+ if (!arr) return;
640+ const n = arr.length, dt = params.decim / SR;
641+ ctx.strokeStyle = color; ctx.lineWidth = width; ctx.setLineDash(dash || []);
642+ ctx.beginPath();
643+ let started = false;
644+ for (let i = 0; i < n; i++) {
645+ const v = arr[i];
646+ if (!Number.isFinite(v)) { started = false; continue; }
647+ const x = ax.X((i + 1) * dt), y = ax.Y(v);
648+ if (!started) { ctx.moveTo(x, y); started = true; } else ctx.lineTo(x, y);
649+ }
650+ ctx.stroke(); ctx.setLineDash([]);
651+ }
652+
653+ // ---------------- plots ----------------
654+ function drawShare() {
655+ const { ctx, W, H } = setup($("cShare"));
656+ const ax = axesTime(ctx, W, H, 0, 1, [0, 0.25, 0.5, 0.75, 1], (v) => fmt(v, 2), params.dur);
657+ chaosSpan(ctx, ax);
658+ for (let k = 2; k < 12; k++) if (members[k]) seriesTime(ctx, ax, members[k].share, css("--faint"), 1);
659+ if (B) seriesTime(ctx, ax, B.share, css("--y"), 1.4);
660+ if (A) seriesTime(ctx, ax, A.share, css("--x"), 1.6);
661+ const S = strikeAt(sel);
662+ if (S) seriesTime(ctx, ax, S.share, strikeColor(sel), 3);
663+ }
664+
665+ function drawLam() {
666+ const { ctx, W, H } = setup($("cLam"));
667+ if (!A || !A.lam) { axesTime(ctx, W, H, 0, 1, [0, 1], () => "", params.dur); return; }
668+ let mx = 0;
669+ for (const v of A.lam) if (v > mx) mx = v;
670+ mx = Math.max(mx * 1.1, 0.05 / params.tau);
671+ const nice = niceStep(mx / 3);
672+ const ticks = [];
673+ for (let v = 0; v <= mx; v += nice) ticks.push(v);
674+ const ax = axesTime(ctx, W, H, -0.15 * mx, mx, ticks, (v) => fmt(v, 0), params.dur);
675+ chaosSpan(ctx, ax);
676+ ctx.strokeStyle = css("--rule"); ctx.lineWidth = 1;
677+ ctx.beginPath(); ctx.moveTo(ax.x0, ax.Y(0)); ctx.lineTo(ax.x1, ax.Y(0)); ctx.stroke();
678+ // energy in dB mapped onto the same vertical span (0 dB at top, -80 dB at bottom)
679+ const e = A.energy, e0 = e[0] || 1, eDb = new Float32Array(e.length);
680+ for (let i = 0; i < e.length; i++) {
681+ const db = 10 * Math.log10(e[i] / e0);
682+ eDb[i] = -0.15 * mx + (Math.max(db, -80) + 80) / 80 * 1.15 * mx;
683+ }
684+ seriesTime(ctx, ax, eDb, css("--muted"), 1, [4, 4]);
685+ seriesTime(ctx, ax, A.lam, css("--x"), 1.6);
686+ }
687+
688+ function niceStep(v) {
689+ const p = Math.pow(10, Math.floor(Math.log10(v)));
690+ const m = v / p;
691+ return (m < 1.5 ? 1 : m < 3.5 ? 2 : m < 7.5 ? 5 : 10) * p;
692+ }
693+
694+ let pathScale = 1;
695+ function drawPath() {
696+ const { ctx, W, H } = setup($("cPath"));
697+ if (!A) return;
698+ let m = 0;
699+ for (let i = 0; i < A.xs.length; i++) m = Math.max(m, Math.abs(A.xs[i]), Math.abs(A.ys[i]));
700+ pathScale = (Math.min(W, H) / 2 - 14) / (m || 1);
701+ drawBowl(ctx, W, H, pathScale);
702+ ctx.strokeStyle = css("--grid"); ctx.lineWidth = 1;
703+ ctx.beginPath(); ctx.moveTo(W / 2, 8); ctx.lineTo(W / 2, H - 8); ctx.moveTo(8, H / 2); ctx.lineTo(W - 8, H / 2); ctx.stroke();
704+ ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono");
705+ ctx.textAlign = "right"; ctx.textBaseline = "top"; ctx.fillText("x", W - 10, H / 2 + 4);
706+ ctx.textAlign = "left"; ctx.fillText("y", W / 2 + 6, 8);
707+ const P = strikeAt(sel) || A;
708+ ctx.strokeStyle = strikeColor(strikeAt(sel) ? sel : 0); ctx.globalAlpha = 0.22; ctx.lineWidth = 0.7;
709+ ctx.beginPath();
710+ const n = P.xs.length, stride = n > 200000 ? 2 : 1;
711+ ctx.moveTo(W / 2 + P.xs[0] * pathScale, H / 2 - P.ys[0] * pathScale);
712+ for (let i = stride; i < n; i += stride) ctx.lineTo(W / 2 + P.xs[i] * pathScale, H / 2 - P.ys[i] * pathScale);
713+ ctx.stroke(); ctx.globalAlpha = 1;
714+ ctx.fillStyle = strikeColor(strikeAt(sel) ? sel : 0); ctx.font = "600 12px " + css("--sans"); ctx.textAlign = "left"; ctx.textBaseline = "bottom";
715+ ctx.fillText("strike " + strikeName(strikeAt(sel) ? sel : 0), 10, H - 8);
716+ }
717+
718+ // potential V(x, y) as bands at fractions of the strike energy
719+ function drawBowl(ctx, W, H, scale) {
720+ const dpr = window.devicePixelRatio || 1;
721+ const nw = Math.max(1, Math.round(W)), nh = Math.max(1, Math.round(H));
722+ const off = document.createElement("canvas");
723+ off.width = nw; off.height = nh;
724+ const octx = off.getContext("2d");
725+ const img = octx.createImageData(nw, nh);
726+ const E = params.E0, w2sq = params.w2 * params.w2;
727+ const levels = [0.1, 0.25, 0.5, 1.0];
728+ const panel = hexToRgb(css("--panel")), ink = hexToRgb(css("--ink")), xcol = hexToRgb(css("--x"));
729+ for (let j = 0; j < nh; j++) {
730+ const y = (H / 2 - j) / scale;
731+ for (let i = 0; i < nw; i++) {
732+ const x = (i - W / 2) / scale;
733+ const v = (0.5 * x * x + 0.5 * w2sq * y * y + x * x * y * y) / E;
734+ let band = 0;
735+ while (band < levels.length && v > levels[band]) band++;
736+ const c = band >= levels.length ? mix(panel, ink, 0.10) : mix(panel, xcol, 0.04 + 0.035 * (levels.length - 1 - band));
737+ const o = (j * nw + i) * 4;
738+ img.data[o] = c[0]; img.data[o + 1] = c[1]; img.data[o + 2] = c[2]; img.data[o + 3] = 255;
739+ }
740+ }
741+ octx.putImageData(img, 0, 0);
742+ ctx.drawImage(off, 0, 0, W, H);
743+ }
744+
745+ // ---------------- physical setup ----------------
746+ let smIdx = 0, smOn = false, smRaf = 0, smLast = 0;
747+ function zigzag(ctx, x0, y0, x1, y1, coils, amp) {
748+ const dx = x1 - x0, dy = y1 - y0, L = Math.hypot(dx, dy) || 1;
749+ const ux = dx / L, uy = dy / L, nx = -uy, ny = ux;
750+ const lead = Math.min(10, L * 0.15);
751+ ctx.beginPath();
752+ ctx.moveTo(x0, y0);
753+ ctx.lineTo(x0 + ux * lead, y0 + uy * lead);
754+ const n = coils * 2, seg = (L - 2 * lead) / n;
755+ for (let k = 1; k < n; k++) {
756+ const s = lead + k * seg, side = k % 2 ? 1 : -1;
757+ ctx.lineTo(x0 + ux * s + nx * amp * side, y0 + uy * s + ny * amp * side);
758+ }
759+ ctx.lineTo(x1 - ux * lead, y1 - uy * lead);
760+ ctx.lineTo(x1, y1);
761+ ctx.stroke();
762+ }
763+ function drawSetup() {
764+ const { ctx, W, H } = setup($("cSetup"));
765+ const fx0 = 22, fx1 = W - 22, fy0 = 18, fy1 = H - 18;
766+ const cx = (fx0 + fx1) / 2, cy = (fy0 + fy1) / 2;
767+ const ink = css("--ink"), muted = css("--muted");
768+ // walls with hatching
769+ ctx.strokeStyle = muted; ctx.lineWidth = 1.5;
770+ ctx.strokeRect(fx0, fy0, fx1 - fx0, fy1 - fy0);
771+ ctx.lineWidth = 1; ctx.strokeStyle = css("--rule");
772+ for (let t = fx0; t < fx1; t += 9) { ctx.beginPath(); ctx.moveTo(t, fy0); ctx.lineTo(t - 6, fy0 - 6); ctx.moveTo(t, fy1); ctx.lineTo(t + 6, fy1 + 6); ctx.stroke(); }
773+ for (let t = fy0; t < fy1; t += 9) { ctx.beginPath(); ctx.moveTo(fx0, t); ctx.lineTo(fx0 - 6, t + 6); ctx.moveTo(fx1, t); ctx.lineTo(fx1 + 6, t - 6); ctx.stroke(); }
774+ // position
775+ let x = 0, y = 0, i = 0;
776+ if (A && smOn) { i = Math.min(A.xs.length - 1, Math.floor(smIdx)); x = A.xs[i]; y = A.ys[i]; }
777+ let m = 0;
778+ if (A) for (let k = 0; k < A.xs.length; k += 7) m = Math.max(m, Math.abs(A.xs[k]), Math.abs(A.ys[k]));
779+ const sc = 0.62 * Math.min(cx - fx0, cy - fy0) / (m || 1);
780+ const px = cx + x * sc, py = cy - y * sc, r = 13;
781+ const w2sq = params ? params.w2 * params.w2 : 1;
782+ const kx = 1 + 2 * y * y, ky = w2sq + 2 * x * x;
783+ const lw = (k) => Math.min(5.5, 1 + 1.3 * Math.log(1 + k));
784+ // trail
785+ if (A && smOn && i > 1) {
786+ const span = Math.round(2 * Math.PI * 1.5 / params.dtOut);
787+ ctx.strokeStyle = css("--x"); ctx.globalAlpha = 0.25; ctx.lineWidth = 1;
788+ ctx.beginPath();
789+ const i0 = Math.max(0, i - span);
790+ ctx.moveTo(cx + A.xs[i0] * sc, cy - A.ys[i0] * sc);
791+ for (let k = i0 + 1; k <= i; k++) ctx.lineTo(cx + A.xs[k] * sc, cy - A.ys[k] * sc);
792+ ctx.stroke(); ctx.globalAlpha = 1;
793+ }
794+ // springs: x springs to rollers on the side walls, y springs to rollers on the top and bottom walls
795+ ctx.lineJoin = "round";
796+ ctx.strokeStyle = css("--x"); ctx.lineWidth = lw(kx);
797+ zigzag(ctx, fx0 + 5, py, px - r, py, 7, 6);
798+ zigzag(ctx, px + r, py, fx1 - 5, py, 7, 6);
799+ ctx.strokeStyle = css("--y"); ctx.lineWidth = lw(ky);
800+ zigzag(ctx, px, fy0 + 5, px, py - r, 5, 6);
801+ zigzag(ctx, px, py + r, px, fy1 - 5, 5, 6);
802+ // rollers
803+ ctx.fillStyle = css("--panel"); ctx.strokeStyle = muted; ctx.lineWidth = 1.5;
804+ for (const [rx, ry] of [[fx0 + 5, py], [fx1 - 5, py], [px, fy0 + 5], [px, fy1 - 5]]) {
805+ ctx.beginPath(); ctx.arc(rx, ry, 4.5, 0, 2 * Math.PI); ctx.fill(); ctx.stroke();
806+ }
807+ // mass
808+ ctx.fillStyle = ink; ctx.beginPath(); ctx.arc(px, py, r, 0, 2 * Math.PI); ctx.fill();
809+ ctx.fillStyle = css("--panel"); ctx.font = "600 12px " + css("--sans"); ctx.textAlign = "center"; ctx.textBaseline = "middle";
810+ ctx.fillText("m", px, py + 0.5);
811+ // strike arrow at rest
812+ if (!smOn && params) {
813+ const a = params.angle, L = 0.42 * Math.min(cx - fx0, cy - fy0);
814+ const ex = cx + Math.cos(a) * L, ey = cy - Math.sin(a) * L;
815+ ctx.strokeStyle = css("--chaos"); ctx.fillStyle = css("--chaos"); ctx.lineWidth = 3;
816+ ctx.beginPath(); ctx.moveTo(cx + Math.cos(a) * (r + 3), cy - Math.sin(a) * (r + 3)); ctx.lineTo(ex, ey); ctx.stroke();
817+ const ha = 0.45, hl = 11;
818+ ctx.beginPath(); ctx.moveTo(ex + Math.cos(a) * 3, ey - Math.sin(a) * 3);
819+ ctx.lineTo(ex - hl * Math.cos(a - ha), ey + hl * Math.sin(a - ha));
820+ ctx.lineTo(ex - hl * Math.cos(a + ha), ey + hl * Math.sin(a + ha)); ctx.closePath(); ctx.fill();
821+ ctx.fillStyle = ink; ctx.font = "600 12px " + css("--sans"); ctx.textAlign = "left"; ctx.textBaseline = "bottom";
822+ ctx.fillText("strike, E = " + (params.E0 < 10 ? fmt(params.E0, 2) : fmt(params.E0, 1)), ex + 8, ey - 2);
823+ }
824+ // stiffness readouts
825+ ctx.font = "11px " + css("--mono"); ctx.textBaseline = "top";
826+ ctx.fillStyle = css("--x"); ctx.textAlign = "left";
827+ ctx.fillText("kx = " + fmt(kx, 2), fx0 + 10, fy0 + 8);
828+ ctx.fillStyle = css("--y"); ctx.textAlign = "right";
829+ ctx.fillText("ky = " + fmt(ky, 2), fx1 - 10, fy0 + 8);
830+ if (smOn && A) {
831+ ctx.fillStyle = muted; ctx.textAlign = "left"; ctx.textBaseline = "bottom";
832+ const e = A.energy[Math.min(A.energy.length - 1, Math.floor(i / params.decim))] || 0;
833+ ctx.fillText("t = " + fmt(i / SR * 1000, 1) + " ms of sound, energy " + fmt(e, 2), fx0 + 10, fy1 - 6);
834+ }
835+ }
836+ function smTick(now) {
837+ if (!smOn || !A) return;
838+ const dtReal = smLast ? Math.min(0.1, (now - smLast) / 1000) : 0;
839+ smLast = now;
840+ smIdx += 3.0 * dtReal / params.dtOut; // 3 time units per second: one x period in about 2 s
841+ if (smIdx >= A.xs.length - 1) { stopSlow(); return; }
842+ drawSetup();
843+ smRaf = requestAnimationFrame(smTick);
844+ }
845+ function startSlow() {
846+ if (!A) return;
847+ smOn = true; smIdx = 0; smLast = 0;
848+ $("slow").setAttribute("aria-pressed", "true"); $("slow").textContent = "Stop slow motion";
849+ smRaf = requestAnimationFrame(smTick);
850+ }
851+ function stopSlow() {
852+ smOn = false; cancelAnimationFrame(smRaf);
853+ $("slow").setAttribute("aria-pressed", "false"); $("slow").textContent = "Slow motion";
854+ drawSetup();
855+ }
856+
857+ function drawSection() {
858+ const { ctx, W, H } = setup($("cSec"));
859+ ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono");
860+ if (!sec) { ctx.textAlign = "center"; ctx.textBaseline = "middle"; ctx.fillText("mapping phase space…", W / 2, H / 2); return; }
861+ const E = params.E0, xm = Math.sqrt(2 * E), pm = Math.sqrt(2 * E);
862+ const sx = (W / 2 - 12) / xm, sy = (H / 2 - 12) / pm;
863+ ctx.strokeStyle = css("--grid"); ctx.lineWidth = 1;
864+ ctx.beginPath(); ctx.ellipse(W / 2, H / 2, xm * sx, pm * sy, 0, 0, 2 * Math.PI); ctx.stroke();
865+ ctx.textAlign = "right"; ctx.textBaseline = "top"; ctx.fillText("x", W - 10, H / 2 + 4);
866+ ctx.textAlign = "left"; ctx.fillText("pₓ", W / 2 + 6, 8);
867+ const pal = [css("--x"), css("--y"), "#3f8f6b", "#8a6bd1", "#c07a2b", "#2f8fa3"];
868+ const thr = 0.02;
869+ sec.forEach((o, i) => {
870+ const chaotic = o.lam > thr;
871+ const isStrike = i === sec.length - 1;
872+ ctx.fillStyle = chaotic ? css("--ink") : pal[i % pal.length];
873+ ctx.globalAlpha = chaotic ? 0.35 : 0.9;
874+ const s = chaotic ? 1 : 1.3;
875+ if (isStrike) return;
876+ for (let k = 0; k < o.pts.length; k += 2) ctx.fillRect(W / 2 + o.pts[k] * sx - s / 2, H / 2 - o.pts[k + 1] * sy - s / 2, s, s);
877+ });
878+ ctx.globalAlpha = 1;
879+ // strike start: x = 0, px = sqrt(2E) cos(angle)
880+ const p0 = Math.sqrt(2 * E) * Math.cos(params.angle);
881+ ctx.strokeStyle = css("--chaos"); ctx.lineWidth = 2.5;
882+ ctx.beginPath(); ctx.arc(W / 2, H / 2 - p0 * sy, 7, 0, 2 * Math.PI); ctx.stroke();
883+ const st = sec[sec.length - 1];
884+ ctx.fillStyle = css("--muted"); ctx.textAlign = "left"; ctx.textBaseline = "bottom";
885+ ctx.fillText(st.lam > thr ? "strike starts in the chaotic sea" : "strike starts on a regular orbit", 10, H - 8);
886+ }
887+
888+ // ---------------- spectrogram ----------------
889+ let specFmax = 0;
890+ function fft(re, im) {
891+ const n = re.length;
892+ for (let i = 1, j = 0; i < n; i++) {
893+ let bit = n >> 1;
894+ for (; j & bit; bit >>= 1) j ^= bit;
895+ j ^= bit;
896+ if (i < j) { [re[i], re[j]] = [re[j], re[i]]; [im[i], im[j]] = [im[j], im[i]]; }
897+ }
898+ for (let len = 2; len <= n; len <<= 1) {
899+ const ang = -2 * Math.PI / len, wr = Math.cos(ang), wi = Math.sin(ang);
900+ for (let i = 0; i < n; i += len) {
901+ let cr = 1, ci = 0;
902+ for (let k = 0; k < len / 2; k++) {
903+ const a = i + k, b = a + len / 2;
904+ const tr = re[b] * cr - im[b] * ci, ti = re[b] * ci + im[b] * cr;
905+ re[b] = re[a] - tr; im[b] = im[a] - ti; re[a] += tr; im[a] += ti;
906+ const ncr = cr * wr - ci * wi; ci = cr * wi + ci * wr; cr = ncr;
907+ }
908+ }
909+ }
910+ }
911+ function spectrogramOf(S) {
912+ if (S.spec) return S.spec;
913+ const sig = S.sig, N = 2048, hop = 256;
914+ const frames = Math.max(1, Math.floor((sig.length - N) / hop));
915+ const f1 = 1 / (2 * Math.PI * params.tau);
916+ specFmax = Math.min(16000, Math.max(2000, 10 * f1 * Math.max(1, params.w2)));
917+ const bins = Math.floor(specFmax / (SR / N));
918+ const win = new Float32Array(N);
919+ for (let i = 0; i < N; i++) win[i] = 0.5 - 0.5 * Math.cos(2 * Math.PI * i / (N - 1));
920+ const mag = new Float32Array(frames * bins);
921+ const re = new Float64Array(N), im = new Float64Array(N);
922+ let mx = -Infinity;
923+ for (let f = 0; f < frames; f++) {
924+ for (let i = 0; i < N; i++) { re[i] = sig[f * hop + i] * win[i]; im[i] = 0; }
925+ fft(re, im);
926+ for (let b = 0; b < bins; b++) {
927+ const v = 10 * Math.log10(re[b] * re[b] + im[b] * im[b] + 1e-20);
928+ mag[f * bins + b] = v; if (v > mx) mx = v;
929+ }
930+ }
931+ S.spec = { mag, frames, bins, mx };
932+ return S.spec;
933+ }
934+ function drawSpec() {
935+ const { ctx, W, H } = setup($("cSpec"));
936+ const ax = axesTime(ctx, W, H, 0, 1, [], () => "", params.dur);
937+ const S = strikeAt(sel) || A;
938+ if (!S || !S.sig) return;
939+ const mx = spectrogramOf(A).mx;
940+ const { mag, frames, bins } = spectrogramOf(S);
941+ const off = document.createElement("canvas");
942+ off.width = frames; off.height = bins;
943+ const octx = off.getContext("2d");
944+ const img = octx.createImageData(frames, bins);
945+ const ink = hexToRgb(css("--ink")), bg = hexToRgb(css("--panel")), hot = hexToRgb(css("--chaos"));
946+ for (let f = 0; f < frames; f++) for (let b = 0; b < bins; b++) {
947+ const v = Math.min(1, Math.max(0, (mag[f * bins + b] - (mx - 70)) / 70));
948+ const c = v < 0.75 ? mix(bg, ink, v / 0.75) : mix(ink, hot, (v - 0.75) / 0.25);
949+ const o = ((bins - 1 - b) * frames + f) * 4;
950+ img.data[o] = c[0]; img.data[o + 1] = c[1]; img.data[o + 2] = c[2]; img.data[o + 3] = 255;
951+ }
952+ octx.putImageData(img, 0, 0);
953+ const tEnd = (frames * 256 + 1024) / SR;
954+ ctx.imageSmoothingEnabled = true;
955+ ctx.drawImage(off, ax.x0 + 1024 / SR / params.dur * (ax.x1 - ax.x0), ax.y0, (tEnd - 1024 / SR) / params.dur * (ax.x1 - ax.x0), ax.y1 - ax.y0);
956+ // frequency labels
957+ ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono"); ctx.textAlign = "right"; ctx.textBaseline = "middle";
958+ const step = niceStep(specFmax / 4);
959+ for (let f = 0; f <= specFmax; f += step) {
960+ const y = ax.y1 - f / specFmax * (ax.y1 - ax.y0);
961+ ctx.fillText(f >= 1000 ? fmt(f / 1000, f % 1000 ? 1 : 0) + "k" : String(f), ax.x0 - 6, y);
962+ }
963+ const ks = strikeAt(sel) ? sel : 0;
964+ ctx.fillStyle = strikeColor(ks); ctx.font = "600 12px " + css("--sans"); ctx.textAlign = "right"; ctx.textBaseline = "top";
965+ ctx.fillText("strike " + strikeName(ks), ax.x1 - 6, ax.y0 + 4);
966+ }
967+ function hexToRgb(h) {
968+ h = h.replace("#", "");
969+ if (h.length === 3) h = h.split("").map((c) => c + c).join("");
970+ const n = parseInt(h, 16);
971+ return [(n >> 16) & 255, (n >> 8) & 255, n & 255];
972+ }
973+ function mix(a, b, t) { return [a[0] + (b[0] - a[0]) * t, a[1] + (b[1] - a[1]) * t, a[2] + (b[2] - a[2]) * t]; }
974+
975+ function drawAll() {
976+ if (!params) return;
977+ drawSetup(); drawPath(); drawSection(); drawShare(); drawLam(); drawSpec(); drawOverlays(null, null);
978+ }
979+
980+ // ---------------- audio ----------------
981+ let actx = null, bufs = [], src = null, playStart = 0, playing = null, raf = 0, audioGain = 1;
982+ function scaleAudio(s) {
983+ const out = new Float32Array(s.length), fade = Math.round(0.01 * SR);
984+ for (let i = 0; i < s.length; i++) {
985+ const tail = s.length - 1 - i;
986+ out[i] = Math.max(-1, Math.min(1, s[i] * audioGain * (tail < fade ? tail / fade : 1)));
987+ }
988+ return out;
989+ }
990+ function prepareAudio() {
991+ let m = 0;
992+ for (const s of [A.sig, B.sig]) for (let i = 0; i < s.length; i++) m = Math.max(m, Math.abs(s[i]));
993+ audioGain = 0.8 / (m || 1); // one gain for every strike, so loudness differences are real
994+ A.audio = scaleAudio(A.sig); B.audio = scaleAudio(B.sig);
995+ bufs = [];
996+ }
997+ function ensureCtx() {
998+ if (!actx) actx = new (window.AudioContext || window.webkitAudioContext)();
999+ if (actx.state === "suspended") actx.resume();
1000+ }
1001+ function toBuffer(data) {
1002+ const b = actx.createBuffer(1, data.length, SR);
1003+ b.copyToChannel(data, 0);
1004+ return b;
1005+ }
1006+ function stopAudio() {
1007+ if (src) { try { src.onended = null; src.stop(); } catch (e) {} src = null; }
1008+ playing = null; cancelAnimationFrame(raf);
1009+ for (const id of ["playA", "playB", "playAB", "playSel"]) if ($(id)) $(id).setAttribute("aria-pressed", "false");
1010+ drawOverlays(null, null);
1011+ }
1012+ function play(k, then) {
1013+ const S = strikeAt(k);
1014+ if (!S || !S.audio) return;
1015+ ensureCtx();
1016+ stopAudio();
1017+ if (!bufs[k]) bufs[k] = toBuffer(S.audio);
1018+ src = actx.createBufferSource();
1019+ src.buffer = bufs[k];
1020+ src.connect(actx.destination);
1021+ const next = then === null || then === undefined ? null : then;
1022+ src.onended = () => {
1023+ src = null;
1024+ if (next !== null) setTimeout(() => play(next, null), 350);
1025+ else stopAudio();
1026+ };
1027+ playStart = actx.currentTime;
1028+ playing = k;
1029+ src.start();
1030+ if (k === 0) $("playA").setAttribute("aria-pressed", "true");
1031+ else if (k === 1) $("playB").setAttribute("aria-pressed", "true");
1032+ if (k === sel && $("playSel")) $("playSel").setAttribute("aria-pressed", "true");
1033+ if (next !== null) $("playAB").setAttribute("aria-pressed", "true");
1034+ tick();
1035+ }
1036+ function tick() {
1037+ if (playing === null) return;
1038+ const t = actx.currentTime - playStart;
1039+ drawOverlays(t, playing);
1040+ raf = requestAnimationFrame(tick);
1041+ }
1042+
1043+ function drawOverlays(t, which) {
1044+ for (const [id, base] of [["cShareO", "cShare"], ["cLamO", "cLam"], ["cSpecO", "cSpec"]]) {
1045+ const { ctx, W } = setup($(id));
1046+ if (t === null) continue;
1047+ const H = $(base).getBoundingClientRect().height;
1048+ const x = PAD.l + Math.min(t, params.dur) / params.dur * (W - PAD.r - PAD.l);
1049+ ctx.strokeStyle = strikeColor(which); ctx.lineWidth = 2;
1050+ ctx.beginPath(); ctx.moveTo(x, PAD.t); ctx.lineTo(x, H - PAD.b); ctx.stroke();
1051+ }
1052+ const { ctx, W, H } = setup($("cPathO"));
1053+ if (t === null || !A) return;
1054+ const S = strikeAt(which);
1055+ if (!S || !S.xs) return;
1056+ const i1 = Math.min(S.xs.length - 1, Math.floor(t * SR));
1057+ const i0 = Math.max(0, i1 - Math.round(0.03 * SR));
1058+ ctx.strokeStyle = strikeColor(which); ctx.lineWidth = 1.6;
1059+ ctx.beginPath();
1060+ ctx.moveTo(W / 2 + S.xs[i0] * pathScale, H / 2 - S.ys[i0] * pathScale);
1061+ for (let i = i0 + 1; i <= i1; i++) ctx.lineTo(W / 2 + S.xs[i] * pathScale, H / 2 - S.ys[i] * pathScale);
1062+ ctx.stroke();
1063+ ctx.fillStyle = css("--chaos");
1064+ ctx.beginPath(); ctx.arc(W / 2 + S.xs[i1] * pathScale, H / 2 - S.ys[i1] * pathScale, 4, 0, 2 * Math.PI); ctx.fill();
1065+ }
1066+
1067+ // ---------------- wiring ----------------
1068+ $("playA").addEventListener("click", () => play(0, null));
1069+ $("playB").addEventListener("click", () => play(1, null));
1070+ $("playAB").addEventListener("click", () => play(0, 1));
1071+ $("stop").addEventListener("click", stopAudio);
1072+ $("slow").addEventListener("click", () => (smOn ? stopSlow() : startSlow()));
1073+ let timer = 0;
1074+ for (const id of ["E0", "w2", "ang", "dl", "t60", "f1", "dur"]) {
1075+ $(id).addEventListener("input", () => { readParams(); clearTimeout(timer); timer = setTimeout(run, 250); });
1076+ }
1077+ document.querySelectorAll("[data-preset]").forEach((b) => b.addEventListener("click", () => {
1078+ setControls(PRESETS[b.dataset.preset]); run();
1079+ }));
1080+ let rt = 0;
1081+ window.addEventListener("resize", () => { clearTimeout(rt); rt = setTimeout(drawAll, 120); });
1082+ try {
1083+ window.matchMedia("(prefers-color-scheme: dark)").addEventListener("change", drawAll);
1084+ new MutationObserver(drawAll).observe(document.documentElement, { attributes: true, attributeFilter: ["data-theme"] });
1085+ } catch (e) {}
1086+
1087+ setControls(PRESETS.strong);
1088+ run();
1089+})();
1090+</script>