Struck FPUT chain: model, experiments, figures, web app
10 changed files+1295−0
.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+__pycache__/
2+*.npz
README.mdadded+44−0View file
@@ -0,0 +1,44 @@
1+# fput-chain
2+
3+A struck Fermi-Pasta-Ulam-Tsingou chain: the many-mode companion of the two-oscillator toy in `../chaos-toy-examples`, and a simple stand-in for the struck plate of `cymbal-elastodynamics/flat-plate-2`.
4+
5+## The model
6+
7+$N$ unit masses on a line between two fixed walls, joined by springs with force law $f(\Delta) = \Delta + \beta \Delta^3$ ($\beta = 1$), and optionally held by on-site springs of stiffness $\mu^2$. The normal modes are $\sin(jk\pi/(N+1))$ with $\omega_k^2 = \mu^2 + 4\sin^2(k\pi/(2(N+1)))$: for $\mu = 0$ the low modes are nearly harmonic (like a string), and $\mu > 0$ detunes them (like a plate). The stiffness of a stretched spring is $1 + 3\beta\Delta^2$, so stretching couples the modes, as stretching the midplane does in a plate. Damping is friction $c_0$ on every mass plus dashpots $c_1$ between neighbours; mode $k$ decays (in amplitude) at $(c_0 + c_1(\omega_k^2 - \mu^2))/2$, the chain's version of the plate's $\sigma_0 + \sigma_1 k^2$ damping. The strike gives a few masses near one end a smooth velocity profile with kinetic energy $E$.
8+
9+`model.py` integrates ensembles of strikes with a fourth-order symplectic scheme (Yoshida) and the exact propagator of the damping, integrates the tangent dynamics for the largest Lyapunov exponent, and computes mode energies and the effective fraction of modes $e^S/N$ ($S$ the entropy of the mode shares; 1 for equal shares, about 0.7 for random equipartition, $1/N$ for a single mode).
10+
11+## Results (2026-10-05; $N = 32$, strike on mass 8, width 2 masses)
12+
13+**Equipartition and chaos vs strike energy, undamped** (`results/scan.png`, 4000 time units, 6 strikes per energy):
14+
15+| E | 0.5 | 0.8 | 2 | 5 | 13 | 32 | 50 |
16+| --- | --- | --- | --- | --- | --- | --- | --- |
17+| $\mu = 0$: Lyapunov exponent | 0.0007 | 0.0015 | 0.006 | 0.021 | 0.050 | 0.10 | 0.13 |
18+| $\mu = 0$: fraction of modes sharing the energy | 0.18 | 0.50 | 0.67 | 0.67 | 0.64 | 0.64 | 0.67 |
19+| $\mu = 1$: Lyapunov exponent | 0 | 0 | 0.0001 | 0.0002 | 0.0009 | 0.013 | 0.037 |
20+| $\mu = 1$: fraction of modes sharing the energy | 0.17 | 0.17 | 0.18 | 0.18 | 0.18 | 0.25 | 0.43 |
21+
22+Below $E \approx 0.5$ the harmonic chain keeps the strike's energy in the modes it started in, the classic FPUT observation. From $E \approx 1$ to 2 it spreads to equipartition, and the Lyapunov exponent grows roughly in proportion to $E$. The detuned chain ($\mu = 1$) holds the energy in its initial modes up to $E \approx 30$, some twenty times more: detuning removes the near-resonances through which the modes exchange energy (`results/cascade.png` shows the mode energies over time). The same effect, much weaker, separated the equal and detuned pairs of the two-oscillator toy, and the plate's spectrum is strongly inharmonic.
23+
24+**Damped strikes and the frozen spectrum** (`results/damped_mu0.png`, `results/damped_mu1.png`; lowest mode decaying at 4e-4, highest at 3e-3 per time unit; 16 strikes whose strengths differ by $10^{-10}$ times the member index):
25+
26+- $\mu = 0$, $E_0 = 0.5$: regular; the energy stays in the lowest four modes and all 16 strikes end identical.
27+- $\mu = 0$, $E_0 = 5$: a real cascade, energy spread over about 70% of the modes for the first 1000 time units, and a positive local Lyapunov exponent (about 0.015) for a few hundred. Yet all 16 strikes end identical. The chaotic phase amplifies differences by about $e^{10}$, which takes $10^{-10}$ to $10^{-6}$, still invisible.
28+- $\mu = 0$, $E_0 = 30$: local exponent up to 0.08 for about 600 time units, amplification far beyond $e^{23} \approx 10^{10}$, and the 16 strikes end with different spectra: the share of energy in each of the main modes differs by about 0.57 decades (about 6 dB, standard deviation across strikes).
29+- $\mu = 1$: no strike up to $E_0 = 30$ ends differently from its neighbours.
30+
31+So an energy cascade is not the same as an audible fingerprint. Whether a tiny difference between strikes grows to a visible one depends on the amplification accumulated over the chaotic phase, $\exp(\int \lambda\, dt)$, against the size of the difference; a strike can share its energy among all modes and still be perfectly repeatable. The damping enters directly: the web app shows that with the highest mode damped twice as fast (ring time 0.15 s instead of 0.3 s at 110 Hz), the strong strike at $E_0 = 30$ no longer separates at all.
32+
33+## Web application
34+
35+`web/index.html` (single page, no build step; published as a private artifact at https://claude.ai/artifact/C6oaNiG2zCRX9VR7Wn4Wfx) simulates 8 strikes in the browser and plays strikes A and B. Controls: strike energy, number of masses, on-site springs, strike position and width, the A/B difference, ring times of the lowest and highest modes, the pitch of the lowest mode, and length. Panels: the chain with displacements drawn sideways and springs drawn by stiffness (with slow motion), the energy of strike A per mode over time, the effective fraction of modes for all 8 strikes, their final spectra, the local Lyapunov exponent with the energy, and a spectrogram of any of the 8 strikes (a selector above it, or the arrow keys, flips between strikes on a shared color scale, highlights the chosen strike in the two plots above, and plays it). Presets: Gentle ($E = 0.5$), Cascade only ($E = 5$), Strong ($E = 30$), Strong, detuned ($\mu = 1$, $E = 30$). Locally: `cd web && python -m http.server`.
36+
37+Time in the app is scaled so the lowest mode sounds at the chosen pitch, so a ring time in seconds means different damping in the chain's own time units for different $\mu$: with $\mu = 1$ the lowest mode is at $\omega \approx 1$ instead of 0.095, and one second of sound is about ten times fewer time units.
38+
39+## Running
40+
41+Experiments run on the office workstation (`~/work/fput-chain`, Python 3.9 with numpy and scipy), plots locally:
42+
43+ WORKERS=32 python experiments.py all # results/*.npz
44+ python plots.py # results/*.png
experiments.pyadded+113−0View file
@@ -0,0 +1,113 @@
1+"""Experiments with the struck FPUT-beta chain. Each writes results/<name>.npz; plots.py draws them.
2+
3+ python experiments.py scan # Lyapunov exponent and equipartition vs strike energy (undamped)
4+ python experiments.py cascade # mode energies over time at three strike energies (undamped)
5+ python experiments.py damped # damped strikes: ensembles of 16 nearby strikes, frozen outcome
6+ python experiments.py all
7+
8+N = 32, beta = 1, strike centered on mass 8 with width 2 masses. mu = 0 (harmonic at low k) or mu = 1
9+(on-site springs: an inharmonic spectrum, closer to the plate).
10+"""
11+
12+from __future__ import annotations
13+
14+import os
15+import sys
16+from multiprocessing import Pool
17+
18+os.environ.setdefault("OMP_NUM_THREADS", "1")
19+
20+import numpy as np # noqa: E402
21+
22+from model import Chain, run, spectral_entropy, strike # noqa: E402
23+
24+HERE = os.path.dirname(os.path.abspath(__file__))
25+RES = os.path.join(HERE, "results")
26+N, CENTER, WIDTH = 32, 8, 2.0
27+
28+
29+def dt_for(E, mu):
30+ return min(0.1, 0.25 / np.sqrt(4 + mu**2 + 12 * E))
31+
32+
33+def _scan_one(args):
34+ mu, E = args
35+ ch = Chain(N=N, mu=mu)
36+ M = 6
37+ u0, v0 = strike(ch, E, CENTER, WIDTH, M=M, scale=1 + 1e-6 * np.arange(M))
38+ dt = dt_for(E, mu)
39+ r = run(ch, u0, v0, 4000.0, dt, lyap=True, record_every=int(round(20 / dt)))
40+ neff = spectral_entropy(r["Ek"]) # (T, M)
41+ return mu, E, r["lyap"], r["t"], neff
42+
43+
44+def scan(pool):
45+ energies = np.geomspace(0.05, 50, 16)
46+ out = {}
47+ for mu, E, lam, t, neff in pool.imap_unordered(_scan_one, [(mu, E) for mu in (0.0, 1.0) for E in energies]):
48+ i = int(np.argmin(abs(energies - E)))
49+ out[f"mu{mu:g}_E{i}_lyap"], out[f"mu{mu:g}_E{i}_neff"], out[f"mu{mu:g}_t"] = lam, neff, t
50+ print(f" scan mu {mu:g} E {E:.3g}: lambda {np.median(lam):.4f}, n_eff/N final {np.median(neff[-1]):.2f}", flush=True)
51+ np.savez(os.path.join(RES, "scan.npz"), energies=energies, **out)
52+
53+
54+def _cascade_one(args):
55+ mu, E = args
56+ ch = Chain(N=N, mu=mu)
57+ u0, v0 = strike(ch, E, CENTER, WIDTH)
58+ dt = dt_for(E, mu)
59+ r = run(ch, u0, v0, 3000.0, dt, record_every=int(round(2 / dt)))
60+ return mu, E, r["t"], r["Ek"][:, 0, :], ch.omega
61+
62+
63+def cascade(pool):
64+ out = {}
65+ jobs = [(mu, E) for mu in (0.0, 1.0) for E in (0.5, 2.0, 10.0)]
66+ for mu, E, t, Ek, om in pool.imap_unordered(_cascade_one, jobs):
67+ k = f"mu{mu:g}_E{E:g}"
68+ out[k + "_t"], out[k + "_Ek"], out[k + "_omega"] = t, Ek, om
69+ print(f" cascade mu {mu:g} E {E:g}", flush=True)
70+ np.savez(os.path.join(RES, "cascade.npz"), **out)
71+
72+
73+def damping_for(ch, sigma_low, sigma_high):
74+ """c0, c1 so that the lowest mode decays (amplitude) at sigma_low and the highest at sigma_high."""
75+ om = ch.omega
76+ a1, aN = om[0] ** 2 - ch.mu**2, om[-1] ** 2 - ch.mu**2
77+ c1 = 2 * (sigma_high - sigma_low) / (aN - a1)
78+ c0 = 2 * sigma_low - c1 * a1
79+ return c0, c1
80+
81+
82+def _damped_one(args):
83+ mu, E0 = args
84+ base = Chain(N=N, mu=mu)
85+ c0, c1 = damping_for(base, 4e-4, 3e-3)
86+ ch = Chain(N=N, mu=mu, c0=c0, c1=c1)
87+ M = 16
88+ u0, v0 = strike(ch, E0, CENTER, WIDTH, M=M, scale=1 + 1e-10 * np.arange(M))
89+ dt = dt_for(E0, mu)
90+ r = run(ch, u0, v0, 8000.0, dt, lyap=True, record_every=int(round(5 / dt)))
91+ return mu, E0, r["t"], r["Ek"], r["E"], r["logsum"], c0, c1
92+
93+
94+def damped(pool):
95+ out = {}
96+ jobs = [(mu, E0) for mu in (0.0, 1.0) for E0 in (0.5, 5.0, 30.0)]
97+ for mu, E0, t, Ek, E, logsum, c0, c1 in pool.imap_unordered(_damped_one, jobs):
98+ k = f"mu{mu:g}_E{E0:g}"
99+ out[k + "_t"], out[k + "_Ek"], out[k + "_E"], out[k + "_logsum"] = t, Ek, E, logsum
100+ out[k + "_c"] = np.array([c0, c1])
101+ print(f" damped mu {mu:g} E0 {E0:g}", flush=True)
102+ np.savez(os.path.join(RES, "damped.npz"), **out)
103+
104+
105+if __name__ == "__main__":
106+ os.makedirs(RES, exist_ok=True)
107+ which = sys.argv[1:] or ["all"]
108+ if "all" in which:
109+ which = ["scan", "cascade", "damped"]
110+ with Pool(int(os.environ.get("WORKERS", "16"))) as pool:
111+ for w in which:
112+ print(w, flush=True)
113+ dict(scan=scan, cascade=cascade, damped=damped)[w](pool)
model.pyadded+152−0View file
@@ -0,0 +1,152 @@
1+"""The struck FPUT-beta chain, vectorized over an ensemble of strikes.
2+
3+N unit masses between two fixed walls (u_0 = u_{N+1} = 0), joined by springs with force law
4+ f(Delta) = Delta + beta Delta^3 (Fermi-Pasta-Ulam-Tsingou, beta model)
5+and optionally held by on-site springs of stiffness mu^2 (which makes the spectrum inharmonic). The
6+energy is
7+ H = sum_j v_j^2 / 2 + sum_{j=0}^{N} (Delta_j^2 / 2 + beta Delta_j^4 / 4) + mu^2 sum_j u_j^2 / 2,
8+with Delta_j = u_{j+1} - u_j. The linear normal modes are sin(j k pi / (N + 1)), k = 1..N, with
9+ omega_k^2 = mu^2 + 4 sin^2(k pi / (2 (N + 1))).
10+
11+Damping: friction c0 on every mass and dashpots c1 between neighbours,
12+ dv_j/dt += -c0 v_j + c1 (v_{j+1} - 2 v_j + v_{j-1}),
13+so that mode k decays (in amplitude) at (c0 + c1 (omega_k^2 - mu^2)) / 2. The dashpot term is the chain's version
14+of the plate's sigma1 gradient damping.
15+
16+Integration: 4th-order symplectic composition (Yoshida) of leapfrog for the conservative part; after each
17+step the velocities are multiplied by the exact propagator exp(-dt D) of the damping (an N x N matrix,
18+precomputed). The tangent dynamics (largest Lyapunov exponent) uses the same splitting, linearized.
19+"""
20+
21+from __future__ import annotations
22+
23+from dataclasses import dataclass
24+
25+import numpy as np
26+from scipy.linalg import expm
27+
28+_C = 2 ** (1 / 3)
29+YOSHIDA = (1 / (2 - _C), -_C / (2 - _C), 1 / (2 - _C))
30+
31+
32+@dataclass(frozen=True)
33+class Chain:
34+ N: int = 32
35+ beta: float = 1.0
36+ mu: float = 0.0
37+ c0: float = 0.0
38+ c1: float = 0.0
39+
40+ @property
41+ def k(self):
42+ return np.arange(1, self.N + 1)
43+
44+ @property
45+ def omega(self):
46+ return np.sqrt(self.mu**2 + 4 * np.sin(self.k * np.pi / (2 * (self.N + 1))) ** 2)
47+
48+ @property
49+ def modes(self):
50+ """Orthonormal mode matrix S[j, k] = sqrt(2 / (N + 1)) sin(j k pi / (N + 1)), j, k = 1..N."""
51+ j = np.arange(1, self.N + 1)
52+ return np.sqrt(2 / (self.N + 1)) * np.sin(np.outer(j, self.k) * np.pi / (self.N + 1))
53+
54+ def damping_matrix(self):
55+ N = self.N
56+ L = -2 * np.eye(N) + np.eye(N, k=1) + np.eye(N, k=-1)
57+ return self.c0 * np.eye(N) - self.c1 * L
58+
59+ def force(self, u):
60+ """u: (M, N) -> force (M, N)."""
61+ up = np.pad(u, ((0, 0), (1, 1)))
62+ d = np.diff(up, axis=1) # (M, N + 1) stretches
63+ f = d + self.beta * d**3
64+ return f[:, 1:] - f[:, :-1] - self.mu**2 * u
65+
66+ def tangent_force(self, u, du):
67+ up = np.pad(u, ((0, 0), (1, 1)))
68+ dup = np.pad(du, ((0, 0), (1, 1)))
69+ d = np.diff(up, axis=1)
70+ dd = np.diff(dup, axis=1)
71+ g = (1 + 3 * self.beta * d**2) * dd
72+ return g[:, 1:] - g[:, :-1] - self.mu**2 * du
73+
74+ def energy(self, u, v):
75+ up = np.pad(u, ((0, 0), (1, 1)))
76+ d = np.diff(up, axis=1)
77+ return 0.5 * (v**2).sum(1) + (0.5 * d**2 + 0.25 * self.beta * d**4).sum(1) + 0.5 * self.mu**2 * (u**2).sum(1)
78+
79+ def mode_energies(self, u, v):
80+ S = self.modes
81+ q, p = u @ S, v @ S
82+ return 0.5 * (p**2 + (self.omega * q) ** 2)
83+
84+
85+def strike(chain: Chain, E, center, width, M=1, scale=None):
86+ """Initial state of a strike: at rest at u = 0, velocity profile exp(-((j - center) / width)^2 / 2)
87+ normalized to kinetic energy E. scale (M,) multiplies the velocity of each member (for ensembles)."""
88+ j = np.arange(1, chain.N + 1)
89+ prof = np.exp(-0.5 * ((j - center) / width) ** 2)
90+ prof *= np.sqrt(2 * E / (prof**2).sum())
91+ v = np.tile(prof, (M, 1))
92+ if scale is not None:
93+ v *= np.asarray(scale)[:, None]
94+ return np.zeros_like(v), v
95+
96+
97+def run(chain: Chain, u0, v0, t_end, dt, lyap=False, record_every=None, renorm_every=10):
98+ """Integrate an ensemble (M, N). Returns final state, the finite-time Lyapunov exponent per member
99+ (if lyap), and records (t, mode energies, total energy, Lyapunov log-sum) every record_every steps."""
100+ u, v = np.array(u0, float), np.array(v0, float)
101+ M = u.shape[0]
102+ D = chain.damping_matrix()
103+ damp = expm(-dt * D) if (chain.c0 or chain.c1) else None
104+ n = int(round(t_end / dt))
105+ if lyap:
106+ rng = np.random.default_rng(0)
107+ du = rng.standard_normal(u.shape)
108+ dv = rng.standard_normal(v.shape)
109+ nrm = np.sqrt((du**2).sum(1) + (dv**2).sum(1))[:, None]
110+ du /= nrm
111+ dv /= nrm
112+ logsum = np.zeros(M)
113+ rec = dict(t=[], Ek=[], E=[], logsum=[])
114+ for i in range(n):
115+ for c in YOSHIDA:
116+ h = c * dt
117+ v = v + 0.5 * h * chain.force(u)
118+ if lyap:
119+ dv = dv + 0.5 * h * chain.tangent_force(u, du)
120+ du = du + h * dv
121+ u = u + h * v
122+ v = v + 0.5 * h * chain.force(u)
123+ if lyap:
124+ dv = dv + 0.5 * h * chain.tangent_force(u, du)
125+ if damp is not None:
126+ v = v @ damp.T
127+ if lyap:
128+ dv = dv @ damp.T
129+ if lyap and (i + 1) % renorm_every == 0:
130+ nrm = np.sqrt((du**2).sum(1) + (dv**2).sum(1))
131+ logsum += np.log(nrm)
132+ du /= nrm[:, None]
133+ dv /= nrm[:, None]
134+ if record_every and (i + 1) % record_every == 0:
135+ rec["t"].append((i + 1) * dt)
136+ rec["Ek"].append(chain.mode_energies(u, v))
137+ rec["E"].append(chain.energy(u, v))
138+ rec["logsum"].append(logsum.copy())
139+ out = dict(u=u, v=v)
140+ if lyap:
141+ out["lyap"] = logsum / (n * dt)
142+ if record_every:
143+ out.update({k: np.array(x) for k, x in rec.items()})
144+ return out
145+
146+
147+def spectral_entropy(Ek):
148+ """Normalized effective number of modes exp(S) / N, S = -sum e ln e, e = E_k / sum E_k (last axis).
149+ 1 at equipartition, 1/N when all energy is in one mode."""
150+ e = Ek / Ek.sum(-1, keepdims=True)
151+ S = -(e * np.log(np.where(e > 0, e, 1))).sum(-1)
152+ return np.exp(S) / Ek.shape[-1]
plots.pyadded+109−0View file
@@ -0,0 +1,109 @@
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+from model import spectral_entropy # noqa: E402
14+
15+HERE = os.path.dirname(os.path.abspath(__file__))
16+RES = os.path.join(HERE, "results")
17+
18+
19+def scan():
20+ d = np.load(os.path.join(RES, "scan.npz"))
21+ energies = d["energies"]
22+ fig, ax = plt.subplots(1, 3, figsize=(17, 4.6))
23+ rows = []
24+ for mu, col in ((0.0, "C0"), (1.0, "C3")):
25+ lam = np.array([np.median(d[f"mu{mu:g}_E{i}_lyap"]) for i in range(len(energies))])
26+ nfin = np.array([np.median(d[f"mu{mu:g}_E{i}_neff"][-1]) for i in range(len(energies))])
27+ ax[0].loglog(energies, lam, "o-", color=col, label=f"mu = {mu:g}")
28+ ax[1].semilogx(energies, nfin, "o-", color=col, label=f"mu = {mu:g}")
29+ for i in range(0, len(energies), 3):
30+ ne = np.median(d[f"mu{mu:g}_E{i}_neff"], axis=1)
31+ t = 20.0 * (np.arange(len(ne)) + 1) # records every 20 time units
32+ ax[2].semilogx(t, ne, color=col, alpha=0.3 + 0.7 * i / len(energies), lw=1)
33+ rows.append((mu, lam, nfin))
34+ ax[0].set(xlabel="strike energy E (N = 32, beta = 1)", ylabel="Lyapunov exponent (4000 time units)",
35+ title="Rate of separation")
36+ ax[1].set(xlabel="strike energy E", ylabel="effective fraction of modes at t = 4000",
37+ title="Equipartition (about 0.7 for random equipartition)")
38+ ax[2].set(xlabel="time", ylabel="effective fraction of modes", title="Spreading over time (darker: more energy)")
39+ for a in ax:
40+ a.grid(True, which="both", alpha=0.3)
41+ ax[0].legend()
42+ ax[1].legend()
43+ fig.tight_layout()
44+ fig.savefig(os.path.join(RES, "scan.png"), dpi=100)
45+ return energies, rows
46+
47+
48+def cascade():
49+ d = np.load(os.path.join(RES, "cascade.npz"))
50+ Es = (0.5, 2.0, 10.0)
51+ fig, ax = plt.subplots(2, 3, figsize=(17, 7.5), sharex=True)
52+ for i, mu in enumerate((0.0, 1.0)):
53+ for j, E in enumerate(Es):
54+ k = f"mu{mu:g}_E{E:g}"
55+ t, Ek, om = d[k + "_t"], d[k + "_Ek"], d[k + "_omega"]
56+ frac = Ek / Ek.sum(1, keepdims=True)
57+ im = ax[i, j].pcolormesh(t, np.arange(1, Ek.shape[1] + 1), np.log10(frac.T + 1e-12),
58+ vmin=-6, vmax=0, cmap="magma", shading="auto", rasterized=True)
59+ ax[i, j].set(title=f"mu = {mu:g}, E = {E:g}", ylabel="mode k" if j == 0 else None)
60+ if i == 1:
61+ ax[i, j].set_xlabel("time")
62+ fig.colorbar(im, ax=ax, label="log10 share of energy in mode k", shrink=0.8)
63+ fig.savefig(os.path.join(RES, "cascade.png"), dpi=100, bbox_inches="tight")
64+
65+
66+def damped(mu):
67+ d = np.load(os.path.join(RES, "damped.npz"))
68+ E0s = (0.5, 5.0, 30.0)
69+ fig, ax = plt.subplots(4, 3, figsize=(17, 14))
70+ summary = []
71+ for j, E0 in enumerate(E0s):
72+ k = f"mu{mu:g}_E{E0:g}"
73+ t, Ek, E, logsum = d[k + "_t"], d[k + "_Ek"], d[k + "_E"], d[k + "_logsum"]
74+ neff = spectral_entropy(Ek) # (T, M)
75+ ax[0, j].semilogy(t, E[:, 0], "k")
76+ ax[0, j].set(title=f"mu = {mu:g}, strike energy E0 = {E0:g}", ylabel="energy")
77+ ax[1, j].plot(t, neff, lw=0.6)
78+ ax[1, j].set(ylabel="effective fraction of modes\n(16 strikes differing by 1e-10)", ylim=(0, 1))
79+ frac = Ek[-1] / Ek[-1].sum(1, keepdims=True) # (M, N)
80+ ks = np.arange(1, frac.shape[1] + 1)
81+ for m in range(frac.shape[0]):
82+ ax[2, j].semilogy(ks, frac[m], lw=0.7)
83+ ax[2, j].set(xlabel="mode k", ylabel="final share of energy", ylim=(1e-8, 1))
84+ w = max(1, int(round(100 / (t[1] - t[0]))))
85+ loc = (logsum[w:, 0] - logsum[:-w, 0]) / (t[w:] - t[:-w])
86+ ax[3, j].plot(t[w:], loc, "k", lw=0.8)
87+ ax[3, j].axhline(0, color="0.6", lw=0.5)
88+ ax[3, j].set(xlabel="time", ylabel="local Lyapunov exponent (100-unit windows)")
89+ top = frac.mean(0) > 1e-3
90+ spread = np.log10(frac[:, top]).std(0)
91+ summary.append((E0, float(np.median(spread)), float(spread.max()), int(top.sum())))
92+ fig.tight_layout()
93+ fig.savefig(os.path.join(RES, f"damped_mu{mu:g}.png"), dpi=95)
94+ return summary
95+
96+
97+if __name__ == "__main__":
98+ if os.path.exists(os.path.join(RES, "scan.npz")):
99+ energies, rows = scan()
100+ for mu, lam, nfin in rows:
101+ print(f"mu {mu:g}: " + " ".join(f"E {E:.2g}: lam {l:.4f} n {n:.2f}" for E, l, n in zip(energies, lam, nfin)))
102+ if os.path.exists(os.path.join(RES, "cascade.npz")):
103+ cascade()
104+ if os.path.exists(os.path.join(RES, "damped.npz")):
105+ for mu in (0.0, 1.0):
106+ for E0, med, mx, n in damped(mu):
107+ print(f"damped mu {mu:g}, E0 {E0:g}: spread of final log10 mode shares across strikes over the {n} modes "
108+ f"holding > 0.1%: median {med:.3f}, max {mx:.3f} (decades)")
109+ print("done")
results/cascade.pngadded+0−0View file
Binary file not shown.
results/damped_mu0.pngadded+0−0View file
Binary file not shown.
results/damped_mu1.pngadded+0−0View file
Binary file not shown.
results/scan.pngadded+0−0View file
Binary file not shown.
web/index.htmladded+875−0View file
@@ -0,0 +1,875 @@
1+<meta charset="utf-8">
2+<title>Struck FPUT Chain</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; --panel: #16202c; --ink: #e3e9f0; --muted: #93a3b5; --rule: #2a3747; --grid: #1f2a37;
27+ --x: #6fa8e8; --y: #d57bae; --chaos: #e8b931; --chaos-wash: rgba(232, 185, 49, 0.14);
28+ --faint: rgba(227, 233, 240, 0.16); --focus: #6fa8e8;
29+ }
30+}
31+:root[data-theme="dark"] {
32+ color-scheme: dark;
33+ --paper: #0f161f; --panel: #16202c; --ink: #e3e9f0; --muted: #93a3b5; --rule: #2a3747; --grid: #1f2a37;
34+ --x: #6fa8e8; --y: #d57bae; --chaos: #e8b931; --chaos-wash: rgba(232, 185, 49, 0.14);
35+ --faint: rgba(227, 233, 240, 0.16); --focus: #6fa8e8;
36+}
37+* { box-sizing: border-box; }
38+body { background: var(--paper); color: var(--ink); font-family: var(--sans); font-size: 15px; line-height: 1.5; padding-inline: 16px; padding-block: 20px 48px; }
39+.wrap { max-width: 1280px; margin: 0 auto; display: grid; gap: 20px; }
40+header { display: grid; gap: 6px; }
41+h1 { font-size: clamp(26px, 4vw, 38px); font-weight: 750; font-stretch: 78%; letter-spacing: -0.01em; line-height: 1.05; margin: 0; text-wrap: balance; }
42+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%; }
43+.lede { max-width: 74ch; color: var(--muted); margin: 0; }
44+.layout { display: grid; grid-template-columns: 300px minmax(0, 1fr); gap: 20px; align-items: start; }
45+@media (max-width: 900px) { .layout { grid-template-columns: 1fr; } }
46+aside { background: var(--panel); border: 1px solid var(--rule); border-radius: 6px; padding: 16px; display: grid; gap: 16px; position: sticky; top: calc(env(safe-area-inset-top, 0px) + 12px); }
47+@media (max-width: 900px) { aside { position: static; } }
48+.label { font-size: 11px; font-weight: 650; letter-spacing: 0.08em; text-transform: uppercase; color: var(--muted); }
49+.presets { display: flex; flex-wrap: wrap; gap: 6px; }
50+.chip { font: inherit; font-size: 13px; padding: 4px 10px; border-radius: 999px; border: 1px solid var(--rule); background: var(--paper); color: var(--ink); cursor: pointer; }
51+.chip:hover { border-color: var(--muted); }
52+.field { display: grid; gap: 4px; }
53+.field .row { display: flex; justify-content: space-between; align-items: baseline; gap: 8px; }
54+.field label { font-size: 14px; font-weight: 550; }
55+.val { font-family: var(--mono); font-size: 13px; font-variant-numeric: tabular-nums; color: var(--ink); white-space: nowrap; }
56+.hint { font-size: 12px; color: var(--muted); margin: 0; }
57+input[type="range"] { width: 100%; accent-color: var(--x); }
58+.player { display: grid; gap: 8px; }
59+.btns { display: grid; grid-template-columns: 1fr 1fr; gap: 8px; }
60+.btn { font: inherit; font-weight: 600; font-size: 14px; padding: 9px 10px; border-radius: 6px; border: 1px solid var(--rule); background: var(--paper); color: var(--ink); cursor: pointer; display: flex; align-items: center; justify-content: center; gap: 8px; }
61+.btn:hover { border-color: var(--muted); }
62+.btn .sw { width: 10px; height: 10px; border-radius: 2px; flex: none; }
63+.btn[aria-pressed="true"] { background: var(--ink); color: var(--paper); border-color: var(--ink); }
64+button:focus-visible, input:focus-visible { outline: 2px solid var(--focus); outline-offset: 2px; }
65+.stats { display: grid; grid-template-columns: auto 1fr; gap: 4px 12px; font-size: 13px; margin: 0; }
66+.stats dt { color: var(--muted); }
67+.stats dd { margin: 0; font-family: var(--mono); font-variant-numeric: tabular-nums; text-align: right; }
68+.status { font-size: 12px; color: var(--muted); min-height: 1.5em; font-family: var(--mono); }
69+main { display: grid; gap: 16px; min-width: 0; }
70+.pair { display: grid; grid-template-columns: 1fr 1fr; gap: 16px; }
71+@media (max-width: 640px) { .pair { grid-template-columns: 1fr; } }
72+figure { margin: 0; display: grid; gap: 6px; min-width: 0; }
73+figcaption { display: grid; gap: 2px; }
74+figcaption b { font-weight: 650; font-size: 15px; }
75+figcaption span { font-size: 13px; color: var(--muted); max-width: 80ch; }
76+.plot { position: relative; width: 100%; background: var(--panel); border: 1px solid var(--rule); border-radius: 6px; }
77+.plot canvas { position: absolute; inset: 0; width: 100%; height: 100%; display: block; }
78+.chainplot { aspect-ratio: 3.4 / 1; min-height: 190px; }
79+.wide { aspect-ratio: 3.2 / 1; min-height: 170px; }
80+.half { aspect-ratio: 1.6 / 1; min-height: 170px; }
81+.short { aspect-ratio: 4.6 / 1; min-height: 140px; }
82+.setuptext { display: grid; gap: 10px; font-size: 14px; grid-template-columns: minmax(0, 1fr) auto; align-items: start; }
83+@media (max-width: 640px) { .setuptext { grid-template-columns: 1fr; } }
84+.setuptext p { margin: 0; max-width: 80ch; }
85+.k { font-family: var(--mono); font-size: 13px; white-space: nowrap; }
86+.strikes { display: flex; flex-wrap: wrap; gap: 6px; align-items: center; }
87+.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; }
88+.sk .sw { width: 8px; height: 8px; border-radius: 2px; }
89+.sk[aria-pressed="true"] { background: var(--ink); color: var(--paper); border-color: var(--ink); }
90+.sk:disabled { opacity: 0.4; cursor: default; }
91+.strikes .btn { padding: 4px 12px; font-size: 13px; margin-left: auto; }
92+.legend { display: flex; flex-wrap: wrap; gap: 14px; font-size: 12px; color: var(--muted); }
93+.legend i { display: inline-block; width: 18px; height: 3px; vertical-align: middle; margin-right: 6px; border-radius: 2px; }
94+footer { color: var(--muted); font-size: 13px; max-width: 80ch; }
95+footer p { margin: 0 0 8px; }
96+</style>
97+
98+<div class="wrap">
99+ <header>
100+ <h1>Struck FPUT Chain
101+ <span class="eq">f(Δ) = Δ + βΔ³ between neighbours, β = 1, walls at both ends</span>
102+ </h1>
103+ <p class="lede">A row of masses joined by springs that stiffen as they stretch, hit near one end and left to ring down: the many-mode version of the coupled-oscillator toy. A gentle strike keeps its energy in the few modes it started in. A strong one spreads it over all of them, chaotically, until damping takes the energy back below the threshold, and the chain then rings on with whatever spectrum the chaos left behind. Strikes A and B differ by a tiny change in strength.</p>
104+ </header>
105+
106+ <div class="layout">
107+ <aside aria-label="Controls">
108+ <div class="player">
109+ <div class="label">Listen</div>
110+ <div class="btns">
111+ <button class="btn" id="playA" type="button"><span class="sw" style="background: var(--x)"></span>Strike A</button>
112+ <button class="btn" id="playB" type="button"><span class="sw" style="background: var(--y)"></span>Strike B</button>
113+ <button class="btn" id="playAB" type="button">A, then B</button>
114+ <button class="btn" id="stop" type="button">Stop</button>
115+ </div>
116+ <div class="status" id="status">Computing…</div>
117+ </div>
118+
119+ <div class="field">
120+ <div class="label">Presets</div>
121+ <div class="presets">
122+ <button class="chip" type="button" data-preset="gentle">Gentle</button>
123+ <button class="chip" type="button" data-preset="cascade">Cascade only</button>
124+ <button class="chip" type="button" data-preset="strong">Strong</button>
125+ <button class="chip" type="button" data-preset="detuned">Strong, detuned</button>
126+ </div>
127+ </div>
128+
129+ <div class="field">
130+ <div class="row"><label for="E0">Strike energy E</label><span class="val" id="E0v"></span></div>
131+ <input id="E0" type="range" min="0" max="1000" step="1">
132+ <p class="hint">Around E ≈ 5 energy spreads over most modes for a while, but the spreading is too brief to amplify a tiny difference between strikes; by E ≈ 30 it is (try a larger A/B difference at E ≈ 5).</p>
133+ </div>
134+ <div class="field">
135+ <div class="row"><label for="N">Masses</label><span class="val" id="Nv"></span></div>
136+ <input id="N" type="range" min="8" max="64" step="1">
137+ </div>
138+ <div class="field">
139+ <div class="row"><label for="mu">On-site springs μ</label><span class="val" id="muv"></span></div>
140+ <input id="mu" type="range" min="0" max="2" step="0.05">
141+ <p class="hint">Springs from each mass to the frame detune the modes (as in a plate) and hold the cascade back: at μ = 1 even E = 30 stays put.</p>
142+ </div>
143+ <div class="field">
144+ <div class="row"><label for="pos">Strike position</label><span class="val" id="posv"></span></div>
145+ <input id="pos" type="range" min="0.05" max="0.5" step="0.01">
146+ </div>
147+ <div class="field">
148+ <div class="row"><label for="wid">Strike width</label><span class="val" id="widv"></span></div>
149+ <input id="wid" type="range" min="0.6" max="6" step="0.1">
150+ </div>
151+ <div class="field">
152+ <div class="row"><label for="dl">Difference between A and B</label><span class="val" id="dlv"></span></div>
153+ <input id="dl" type="range" min="-14" max="-2" step="0.5">
154+ </div>
155+ <div class="field">
156+ <div class="row"><label for="tlo">Ring time, lowest mode</label><span class="val" id="tlov"></span></div>
157+ <input id="tlo" type="range" min="0.5" max="10" step="0.1">
158+ </div>
159+ <div class="field">
160+ <div class="row"><label for="thi">Ring time, highest mode</label><span class="val" id="thiv"></span></div>
161+ <input id="thi" type="range" min="0.05" max="5" step="0.05">
162+ <p class="hint">Friction on the masses and dashpots between them: high modes die faster.</p>
163+ </div>
164+ <div class="field">
165+ <div class="row"><label for="f1">Pitch of the lowest mode</label><span class="val" id="f1v"></span></div>
166+ <input id="f1" type="range" min="40" max="440" step="1">
167+ </div>
168+ <div class="field">
169+ <div class="row"><label for="dur">Length</label><span class="val" id="durv"></span></div>
170+ <input id="dur" type="range" min="1" max="6" step="0.5">
171+ </div>
172+
173+ <dl class="stats" aria-live="polite">
174+ <dt>Chaotic phase ends</dt><dd id="sChaos">–</dd>
175+ <dt>Peak separation rate</dt><dd id="sLam">–</dd>
176+ <dt>Modes sharing the energy, peak</dt><dd id="sNeff">–</dd>
177+ <dt>Final spectra, A vs 8 strikes</dt><dd id="sSpread">–</dd>
178+ </dl>
179+ </aside>
180+
181+ <main>
182+ <figure>
183+ <figcaption><b>Physical setup</b><span>Masses on a line between two walls, each joined to the next by a spring and a dashpot. The masses move along the line; here each mass's displacement is drawn sideways so the waves are visible.</span></figcaption>
184+ <div class="plot chainplot"><canvas id="cChain"></canvas></div>
185+ <div class="setuptext">
186+ <p>The springs push back with f(Δ) = Δ + βΔ³ when stretched by Δ, so each one stiffens as it stretches: its stiffness is <span class="k">1 + 3βΔ²</span> (line thickness). That is what couples the modes, exactly as stretching the midplane couples the modes of a plate. The strike gives the marked masses a sudden velocity; the sound is the velocity of the mass at the pickup ring.</p>
187+ <button class="btn" id="slow" type="button" aria-pressed="false">Slow motion</button>
188+ </div>
189+ </figure>
190+
191+ <figure>
192+ <figcaption><b>Where the energy is</b><span>Share of strike A's energy in each normal mode over time (log scale, lowest mode at the bottom). A cascade shows as energy climbing to the higher modes.</span></figcaption>
193+ <div class="plot wide"><canvas id="cModes"></canvas><canvas id="cModesO"></canvas></div>
194+ </figure>
195+
196+ <div class="pair">
197+ <figure>
198+ <figcaption><b>How many modes share it</b><span>Effective fraction of modes holding the energy (1 for equal shares; about 0.7 for random equipartition), for 8 strikes. Shaded: the chaotic phase of strike A.</span></figcaption>
199+ <div class="plot half"><canvas id="cNeff"></canvas><canvas id="cNeffO"></canvas></div>
200+ </figure>
201+ <figure>
202+ <figcaption><b>The spectrum it rings with</b><span>Each strike's final share of energy per mode. Where the strikes differ, the fingerprint is chaos's doing.</span></figcaption>
203+ <div class="plot half"><canvas id="cFinal"></canvas></div>
204+ </figure>
205+ </div>
206+ <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>6 more strikes</span></div>
207+
208+ <figure>
209+ <figcaption><b>How fast nearby strikes separate</b><span>Local Lyapunov exponent of strike A (per second). The dashed line is the energy in dB.</span></figcaption>
210+ <div class="plot short"><canvas id="cLam"></canvas><canvas id="cLamO"></canvas></div>
211+ </figure>
212+
213+ <figure>
214+ <figcaption><b>What you hear</b><span>Spectrogram of the selected strike at the pickup. Click a strike, or use the ← and → keys, to flip between them; the color scale is shared, so differences in level are real. The selected strike is also drawn bold in the two plots above.</span></figcaption>
215+ <div class="strikes" role="group" aria-label="Choose a strike" id="strikeSel"></div>
216+ <div class="plot wide"><canvas id="cSpec"></canvas><canvas id="cSpecO"></canvas></div>
217+ </figure>
218+
219+ <footer>
220+ <p>Integrated in your browser with a fourth-order symplectic scheme (Yoshida); the damping is applied once per audio sample through its exact propagator, which the chain's sine modes diagonalize. Time is scaled so that the lowest mode sounds at the chosen pitch. Fermi, Pasta, Ulam and Tsingou (1955) expected such a chain to share its energy among all modes and found that at low energy it does not; above a threshold it does, and the sharing is chaotic.</p>
221+ </footer>
222+ </main>
223+ </div>
224+</div>
225+
226+<script id="worker-src" type="text/plain">
227+const C = Math.cbrt(2), YS = [1 / (2 - C), -C / (2 - C), 1 / (2 - C)];
228+
229+function setupChain(p) {
230+ const N = p.N, S = new Float64Array(N * N), om = new Float64Array(N), lap = new Float64Array(N);
231+ const c = Math.sqrt(2 / (N + 1));
232+ for (let j = 0; j < N; j++) for (let k = 0; k < N; k++) S[j * N + k] = c * Math.sin((j + 1) * (k + 1) * Math.PI / (N + 1));
233+ for (let k = 0; k < N; k++) { lap[k] = 4 * Math.sin((k + 1) * Math.PI / (2 * (N + 1))) ** 2; om[k] = Math.sqrt(p.mu * p.mu + lap[k]); }
234+ // damping: amplitude decay of mode k = (c0 + c1 lap_k) / 2; propagator over one output sample
235+ const P = new Float64Array(N * N);
236+ const ex = new Float64Array(N);
237+ for (let k = 0; k < N; k++) ex[k] = Math.exp(-p.dtOut * (p.c0 + p.c1 * lap[k]));
238+ for (let i = 0; i < N; i++) for (let j = 0; j < N; j++) {
239+ let s = 0;
240+ for (let k = 0; k < N; k++) s += S[i * N + k] * ex[k] * S[j * N + k];
241+ P[i * N + j] = s;
242+ }
243+ return { S, om, P };
244+}
245+
246+function strikeRun(p, ch, scale, opts) {
247+ const N = p.N, mu2 = p.mu * p.mu, beta = 1;
248+ const u = new Float64Array(N), v = new Float64Array(N), f = new Float64Array(N);
249+ const du = new Float64Array(N), dv = new Float64Array(N), g = new Float64Array(N);
250+ const tmp = new Float64Array(N);
251+ // strike profile
252+ let n2 = 0;
253+ for (let j = 0; j < N; j++) { const w = Math.exp(-0.5 * ((j + 1 - p.center) / p.width) ** 2); v[j] = w; n2 += w * w; }
254+ const amp = Math.sqrt(2 * p.E0 / n2) * scale;
255+ for (let j = 0; j < N; j++) v[j] *= amp;
256+ if (opts.tangent) { for (let j = 0; j < N; j++) { du[j] = Math.sin(1.7 * j + 0.3); dv[j] = Math.cos(2.3 * j + 0.1); } normalize(du, dv); }
257+ const nsub = Math.max(1, Math.ceil(p.dtOut / p.hmax)), h0 = p.dtOut / nsub;
258+ const nOut = p.nOut, K = p.decim, nD = Math.floor(nOut / K);
259+ const sig = opts.audio ? new Float32Array(nOut) : null;
260+ const snapN = opts.snap ? Math.min(nOut, p.snapN) : 0;
261+ const snap = opts.snap ? new Float32Array(snapN * N) : null;
262+ const shares = opts.modes ? new Float32Array(nD * N) : null;
263+ const neff = new Float32Array(nD);
264+ const energy = new Float32Array(nD);
265+ const logsum = opts.tangent ? new Float64Array(nD) : null;
266+ const finalShare = new Float32Array(N);
267+ let L = 0;
268+ const pk = p.pickup - 1;
269+
270+ function force(uu, out) {
271+ let dPrev = uu[0], fPrev = dPrev + beta * dPrev * dPrev * dPrev; // left wall: stretch u_1 - 0
272+ for (let j = 0; j < N; j++) {
273+ const d = (j < N - 1 ? uu[j + 1] : 0) - uu[j];
274+ const fn = d + beta * d * d * d;
275+ out[j] = fn - fPrev - mu2 * uu[j];
276+ fPrev = fn;
277+ }
278+ }
279+ function tforce(uu, dd, out) {
280+ let sPrev = (1 + 3 * beta * uu[0] * uu[0]) * dd[0];
281+ for (let j = 0; j < N; j++) {
282+ const d = (j < N - 1 ? uu[j + 1] : 0) - uu[j];
283+ const dl = (j < N - 1 ? dd[j + 1] : 0) - dd[j];
284+ const sn = (1 + 3 * beta * d * d) * dl;
285+ out[j] = sn - sPrev - mu2 * dd[j];
286+ sPrev = sn;
287+ }
288+ }
289+ function normalize(a, b) {
290+ let s = 0;
291+ for (let j = 0; j < N; j++) s += a[j] * a[j] + b[j] * b[j];
292+ const n = Math.sqrt(s);
293+ for (let j = 0; j < N; j++) { a[j] /= n; b[j] /= n; }
294+ return n;
295+ }
296+ function modeShares(out) {
297+ let tot = 0;
298+ for (let k = 0; k < N; k++) {
299+ let q = 0, pp = 0;
300+ for (let j = 0; j < N; j++) { const s = ch.S[j * N + k]; q += s * u[j]; pp += s * v[j]; }
301+ const e = 0.5 * (pp * pp + ch.om[k] * ch.om[k] * q * q);
302+ out[k] = e; tot += e;
303+ }
304+ for (let k = 0; k < N; k++) out[k] /= (tot || 1);
305+ return tot;
306+ }
307+ const es = new Float64Array(N);
308+ force(u, f);
309+ if (opts.tangent) tforce(u, du, g);
310+ for (let i = 0; i < nOut; i++) {
311+ for (let s = 0; s < nsub; s++) {
312+ for (let c = 0; c < 3; c++) {
313+ const h = YS[c] * h0;
314+ for (let j = 0; j < N; j++) v[j] += 0.5 * h * f[j];
315+ if (opts.tangent) { for (let j = 0; j < N; j++) { dv[j] += 0.5 * h * g[j]; du[j] += h * dv[j]; } }
316+ for (let j = 0; j < N; j++) u[j] += h * v[j];
317+ force(u, f);
318+ for (let j = 0; j < N; j++) v[j] += 0.5 * h * f[j];
319+ if (opts.tangent) { tforce(u, du, g); for (let j = 0; j < N; j++) dv[j] += 0.5 * h * g[j]; }
320+ }
321+ }
322+ // damping over one output sample (exact propagator)
323+ for (let a = 0; a < N; a++) { let s = 0; const row = a * N; for (let b = 0; b < N; b++) s += ch.P[row + b] * v[b]; tmp[a] = s; }
324+ v.set(tmp);
325+ if (opts.tangent) {
326+ for (let a = 0; a < N; a++) { let s = 0; const row = a * N; for (let b = 0; b < N; b++) s += ch.P[row + b] * dv[b]; tmp[a] = s; }
327+ dv.set(tmp);
328+ L += Math.log(normalize(du, dv));
329+ tforce(u, du, g);
330+ }
331+ if (sig) sig[i] = v[pk];
332+ if (i < snapN) for (let j = 0; j < N; j++) snap[i * N + j] = u[j];
333+ if ((i + 1) % K === 0) {
334+ const jd = (i + 1) / K - 1;
335+ if (jd < nD) {
336+ const tot = modeShares(es);
337+ let S = 0;
338+ for (let k = 0; k < N; k++) if (es[k] > 0) S -= es[k] * Math.log(es[k]);
339+ neff[jd] = Math.exp(S) / N;
340+ energy[jd] = tot;
341+ if (shares) for (let k = 0; k < N; k++) shares[jd * N + k] = es[k];
342+ if (logsum) logsum[jd] = L;
343+ }
344+ }
345+ }
346+ modeShares(es);
347+ for (let k = 0; k < N; k++) finalShare[k] = es[k];
348+ return { sig, snap, shares, neff, energy, logsum, finalShare };
349+}
350+
351+onmessage = (ev) => {
352+ const { gen, p } = ev.data;
353+ const ch = setupChain(p);
354+ postMessage({ kind: "progress", gen, text: "Computing strike A…" });
355+ const a = strikeRun(p, ch, 1, { audio: true, snap: true, modes: true, tangent: true });
356+ postMessage({ kind: "progress", gen, text: "Computing strike B…" });
357+ const b = strikeRun(p, ch, 1 + p.delta, { audio: true, snap: false, modes: false, tangent: false });
358+ postMessage({ kind: "strikes", gen, a, b, om: ch.om });
359+ for (let k = 2; k < 8; k++) {
360+ postMessage({ kind: "progress", gen, text: "Computing strike " + (k + 1) + " of 8…" });
361+ const r = strikeRun(p, ch, 1 + k * p.delta, { audio: true, snap: false, modes: false, tangent: false });
362+ postMessage({ kind: "member", gen, k, neff: r.neff, finalShare: r.finalShare, sig: r.sig }, [r.sig.buffer]);
363+ }
364+ postMessage({ kind: "done", gen });
365+};
366+</script>
367+
368+<script>
369+(() => {
370+ const SR = 44100;
371+ const $ = (id) => document.getElementById(id);
372+ const css = (name) => getComputedStyle(document.documentElement).getPropertyValue(name).trim();
373+ const fmt = (v, d = 2) => Number(v).toFixed(d);
374+
375+ const PRESETS = {
376+ gentle: { E0: 0.5, N: 32, mu: 0, pos: 0.25, wid: 2, dl: -10, tlo: 3, thi: 0.3, f1: 110, dur: 4 },
377+ cascade: { E0: 5, N: 32, mu: 0, pos: 0.25, wid: 2, dl: -10, tlo: 3, thi: 0.3, f1: 110, dur: 4 },
378+ strong: { E0: 30, N: 32, mu: 0, pos: 0.25, wid: 2, dl: -10, tlo: 3, thi: 0.3, f1: 110, dur: 4 },
379+ detuned: { E0: 30, N: 32, mu: 1, pos: 0.25, wid: 2, dl: -10, tlo: 3, thi: 0.3, f1: 110, dur: 4 },
380+ };
381+ const eFromSlider = (v) => Math.pow(10, -1.3 + 3.3 * v / 1000);
382+ const sliderFromE = (E) => Math.round((Math.log10(E) + 1.3) / 3.3 * 1000);
383+
384+ function readParams() {
385+ const E0 = eFromSlider(+$("E0").value), N = +$("N").value, mu = +$("mu").value;
386+ const pos = +$("pos").value, wid = +$("wid").value, dl = +$("dl").value;
387+ const tlo = +$("tlo").value, thi = Math.min(+$("thi").value, tlo), f1 = +$("f1").value, dur = +$("dur").value;
388+ const center = Math.max(1, Math.round(pos * (N + 1)));
389+ $("E0v").textContent = (E0 < 10 ? fmt(E0, 2) : fmt(E0, 1)) + " (" + fmt(E0 / N, 3) + " per mass)";
390+ $("Nv").textContent = String(N);
391+ $("muv").textContent = fmt(mu, 2);
392+ $("posv").textContent = "mass " + center + " of " + N;
393+ $("widv").textContent = fmt(wid, 1) + " masses";
394+ $("dlv").textContent = "1e" + dl + " in strength";
395+ $("tlov").textContent = fmt(tlo, 1) + " s";
396+ $("thiv").textContent = fmt(thi, 2) + " s";
397+ const lap1 = 4 * Math.sin(Math.PI / (2 * (N + 1))) ** 2, lapN = 4 * Math.sin(N * Math.PI / (2 * (N + 1))) ** 2;
398+ const om1 = Math.sqrt(mu * mu + lap1), omN = Math.sqrt(mu * mu + lapN);
399+ $("f1v").textContent = Math.round(f1) + " Hz (top: " + Math.round(f1 * omN / om1) + " Hz)";
400+ $("durv").textContent = fmt(dur, 1) + " s";
401+ const tau = om1 / (2 * Math.PI * f1); // seconds per time unit
402+ // amplitude decay rates per time unit, then c0, c1 from (c0 + c1 lap) / 2
403+ const sLo = Math.log(1000) / tlo * tau, sHi = Math.log(1000) / thi * tau;
404+ const c1 = 2 * (sHi - sLo) / (lapN - lap1), c0 = Math.max(0, 2 * sLo - c1 * lap1);
405+ const nOut = Math.round(dur * SR);
406+ return {
407+ E0, N, mu, center, width: wid, delta: Math.pow(10, dl), tau, c0, c1, dur, f1, omN, om1,
408+ dtOut: 1 / SR / tau, nOut, decim: Math.round(SR / 400), pickup: Math.max(1, Math.round(0.38 * (N + 1))),
409+ hmax: Math.min(0.2, 0.45 / Math.sqrt(4 + mu * mu + 12 * E0)), snapN: Math.round(0.12 * SR),
410+ };
411+ }
412+ function setControls(p) {
413+ $("E0").value = sliderFromE(p.E0); $("N").value = p.N; $("mu").value = p.mu; $("pos").value = p.pos; $("wid").value = p.wid;
414+ $("dl").value = p.dl; $("tlo").value = p.tlo; $("thi").value = p.thi; $("f1").value = p.f1; $("dur").value = p.dur;
415+ }
416+
417+ // ---------------- worker ----------------
418+ const workerURL = URL.createObjectURL(new Blob([$("worker-src").textContent], { type: "text/javascript" }));
419+ let worker = null, gen = 0, params = null, A = null, B = null, OM = null, members = [];
420+ let sel = 0, audioGain = 1;
421+ const strikeAt = (k) => (k === 0 ? A : k === 1 ? B : members[k] || null);
422+ const strikeColor = (k) => (k === 0 ? css("--x") : k === 1 ? css("--y") : css("--ink"));
423+ const strikeName = (k) => (k === 0 ? "A" : k === 1 ? "B" : String(k + 1));
424+
425+ function run() {
426+ stopAudio(); stopSlow(true);
427+ params = readParams();
428+ gen++;
429+ if (worker) worker.terminate();
430+ worker = new Worker(workerURL);
431+ worker.onmessage = onResult;
432+ A = B = null; members = []; bufs = [];
433+ $("status").textContent = "Computing…";
434+ worker.postMessage({ gen, p: params });
435+ drawAll();
436+ }
437+ function onResult(ev) {
438+ const m = ev.data;
439+ if (m.gen !== gen) return;
440+ if (m.kind === "progress") { $("status").textContent = m.text; return; }
441+ if (m.kind === "strikes") {
442+ A = m.a; B = m.b; OM = m.om;
443+ A.lam = localLyap(A.logsum);
444+ prepareAudio();
445+ } else if (m.kind === "member") {
446+ members[m.k] = m;
447+ m.audio = scaleAudio(m.sig);
448+ } else if (m.kind === "done") {
449+ $("status").textContent = "Ready. Strikes differ by 1e" + $("dl").value + " in strength.";
450+ }
451+ updateStats();
452+ renderSelector();
453+ drawAll();
454+ }
455+
456+ function renderSelector() {
457+ const box = $("strikeSel");
458+ if (!box.childElementCount) {
459+ for (let k = 0; k < 8; k++) {
460+ const b = document.createElement("button");
461+ b.type = "button"; b.className = "sk"; b.dataset.k = k; b.id = "sk" + k;
462+ b.innerHTML = '<span class="sw"></span>' + strikeName(k);
463+ b.addEventListener("click", () => selectStrike(k));
464+ box.appendChild(b);
465+ }
466+ const pb = document.createElement("button");
467+ pb.type = "button"; pb.className = "btn"; pb.id = "playSel"; pb.textContent = "Play this strike";
468+ pb.addEventListener("click", () => play(sel, null));
469+ box.appendChild(pb);
470+ }
471+ for (let k = 0; k < 8; k++) {
472+ const b = $("sk" + k);
473+ b.disabled = !strikeAt(k);
474+ b.setAttribute("aria-pressed", String(k === sel));
475+ b.querySelector(".sw").style.background = strikeColor(k);
476+ b.title = strikeAt(k) ? "Strike " + strikeName(k) + ": strength × (1 + " + k + " × 1e" + $("dl").value + ")" : "computing…";
477+ }
478+ $("playSel").disabled = !strikeAt(sel);
479+ }
480+ function selectStrike(k) {
481+ if (!strikeAt(k)) return;
482+ sel = k;
483+ renderSelector();
484+ drawNeff(); drawFinal(); drawSpec();
485+ }
486+ document.addEventListener("keydown", (e) => {
487+ if (e.target && (e.target.tagName === "INPUT" || e.target.tagName === "TEXTAREA")) return;
488+ if (e.key !== "ArrowRight" && e.key !== "ArrowLeft") return;
489+ const dir = e.key === "ArrowRight" ? 1 : -1;
490+ for (let i = 1; i <= 8; i++) {
491+ const k = (sel + dir * i + 8) % 8;
492+ if (strikeAt(k)) { selectStrike(k); e.preventDefault(); return; }
493+ }
494+ });
495+
496+ function localLyap(logsum) {
497+ const dtD = params.decim * params.dtOut;
498+ const w = Math.max(1, Math.round(200 / dtD)); // windows of 200 time units
499+ const out = new Float32Array(logsum.length);
500+ for (let i = 0; i < logsum.length; i++) {
501+ const j = Math.max(0, i - w);
502+ out[i] = i > j ? (logsum[i] - logsum[j]) / ((i - j) * dtD) / params.tau : NaN;
503+ }
504+ return out;
505+ }
506+ function chaosEnd(lam) {
507+ const thr = 0.005 / params.tau;
508+ let last = -1;
509+ for (let i = 0; i < lam.length; i++) if (lam[i] > thr) last = i;
510+ return last < 0 ? null : (last + 1) * params.decim / SR;
511+ }
512+ function updateStats() {
513+ if (!A) { for (const id of ["sChaos", "sLam", "sNeff", "sSpread"]) $(id).textContent = "–"; return; }
514+ const end = chaosEnd(A.lam);
515+ $("sChaos").textContent = end === null ? "never chaotic" : end > params.dur - 0.15 ? "still chaotic at the end" : fmt(end, 2) + " s";
516+ let peak = 0; for (const v of A.lam) if (v > peak) peak = v;
517+ $("sLam").textContent = peak > 0.005 / params.tau ? fmt(peak, 0) + " /s" : "≈ 0";
518+ let nmax = 0; for (const v of A.neff) nmax = Math.max(nmax, v);
519+ $("sNeff").textContent = fmt(nmax * 100, 0) + "% of " + params.N;
520+ // spread of final spectra: median over modes holding > 0.1% of the std of log10 share (decades)
521+ const all = [A.finalShare, B.finalShare];
522+ for (let k = 2; k < 8; k++) if (members[k]) all.push(members[k].finalShare);
523+ if (all.length < 8) { $("sSpread").textContent = "computing…"; return; }
524+ const N = params.N, sd = [];
525+ for (let k = 0; k < N; k++) {
526+ const vals = all.map((s) => Math.log10(Math.max(s[k], 1e-12)));
527+ const meanShare = all.reduce((s, x) => s + x[k], 0) / all.length;
528+ if (meanShare < 1e-3) continue;
529+ const m = vals.reduce((s, v) => s + v, 0) / vals.length;
530+ sd.push(Math.sqrt(vals.reduce((s, v) => s + (v - m) ** 2, 0) / vals.length));
531+ }
532+ sd.sort((a, b) => a - b);
533+ const med = sd.length ? sd[Math.floor(sd.length / 2)] : 0;
534+ $("sSpread").textContent = med < 0.01 ? "identical" : "differ by " + fmt(10 * med, 1) + " dB";
535+ }
536+
537+ // ---------------- canvas helpers ----------------
538+ function setup(canvas) {
539+ const r = canvas.getBoundingClientRect();
540+ const dpr = window.devicePixelRatio || 1;
541+ const w = Math.max(1, Math.round(r.width * dpr)), h = Math.max(1, Math.round(r.height * dpr));
542+ if (canvas.width !== w || canvas.height !== h) { canvas.width = w; canvas.height = h; }
543+ const ctx = canvas.getContext("2d");
544+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
545+ ctx.clearRect(0, 0, r.width, r.height);
546+ return { ctx, W: r.width, H: r.height };
547+ }
548+ const PAD = { l: 44, r: 12, t: 10, b: 26 };
549+ function niceStep(v) { const p = Math.pow(10, Math.floor(Math.log10(v))), m = v / p; return (m < 1.5 ? 1 : m < 3.5 ? 2 : m < 7.5 ? 5 : 10) * p; }
550+ function axesTime(ctx, W, H, ymin, ymax, yticks, yfmt, dur) {
551+ const x0 = PAD.l, x1 = W - PAD.r, y0 = PAD.t, y1 = H - PAD.b;
552+ ctx.strokeStyle = css("--grid"); ctx.lineWidth = 1;
553+ ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono");
554+ ctx.textAlign = "right"; ctx.textBaseline = "middle";
555+ for (const v of yticks) {
556+ const y = y1 - (v - ymin) / (ymax - ymin) * (y1 - y0);
557+ ctx.beginPath(); ctx.moveTo(x0, y); ctx.lineTo(x1, y); ctx.stroke();
558+ ctx.fillText(yfmt(v), x0 - 6, y);
559+ }
560+ ctx.textAlign = "center"; ctx.textBaseline = "top";
561+ const step = dur <= 2 ? 0.25 : dur <= 5 ? 0.5 : 1;
562+ for (let t = 0; t <= dur + 1e-9; t += step) {
563+ const x = x0 + t / dur * (x1 - x0);
564+ ctx.beginPath(); ctx.moveTo(x, y0); ctx.lineTo(x, y1); ctx.stroke();
565+ ctx.fillText(fmt(t, step < 1 ? (step < 0.5 ? 2 : 1) : 0) + (t === 0 ? " s" : ""), x, y1 + 6);
566+ }
567+ return { x0, x1, y0, y1, X: (t) => x0 + t / dur * (x1 - x0), Y: (v) => y1 - (v - ymin) / (ymax - ymin) * (y1 - y0) };
568+ }
569+ function chaosSpan(ctx, ax) {
570+ if (!A || !A.lam) return;
571+ const end = chaosEnd(A.lam);
572+ if (end === null) return;
573+ ctx.fillStyle = css("--chaos-wash");
574+ ctx.fillRect(ax.x0, ax.y0, ax.X(end) - ax.x0, ax.y1 - ax.y0);
575+ }
576+ function seriesTime(ctx, ax, arr, color, width, dash) {
577+ if (!arr) return;
578+ const dt = params.decim / SR;
579+ ctx.strokeStyle = color; ctx.lineWidth = width; ctx.setLineDash(dash || []);
580+ ctx.beginPath();
581+ let started = false;
582+ for (let i = 0; i < arr.length; i++) {
583+ const v = arr[i];
584+ if (!Number.isFinite(v)) { started = false; continue; }
585+ const x = ax.X((i + 1) * dt), y = ax.Y(v);
586+ if (!started) { ctx.moveTo(x, y); started = true; } else ctx.lineTo(x, y);
587+ }
588+ ctx.stroke(); ctx.setLineDash([]);
589+ }
590+ function hexToRgb(h) { h = h.replace("#", ""); if (h.length === 3) h = h.split("").map((c) => c + c).join(""); const n = parseInt(h, 16); return [(n >> 16) & 255, (n >> 8) & 255, n & 255]; }
591+ 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]; }
592+ function heatColor(v, bg, ink, hot) { return v < 0.72 ? mix(bg, ink, v / 0.72) : mix(ink, hot, (v - 0.72) / 0.28); }
593+
594+ // ---------------- chain (physical setup) ----------------
595+ let smOn = false, smIdx = 0, smRaf = 0, smLast = 0;
596+ function zigzag(ctx, x0, y0, x1, y1, coils, amp) {
597+ const dx = x1 - x0, dy = y1 - y0, L = Math.hypot(dx, dy) || 1, ux = dx / L, uy = dy / L, nx = -uy, ny = ux;
598+ const lead = Math.min(3, L * 0.15);
599+ ctx.beginPath(); ctx.moveTo(x0, y0); ctx.lineTo(x0 + ux * lead, y0 + uy * lead);
600+ const n = coils * 2, seg = (L - 2 * lead) / n;
601+ for (let k = 1; k < n; k++) { const s = lead + k * seg, side = k % 2 ? 1 : -1; ctx.lineTo(x0 + ux * s + nx * amp * side, y0 + uy * s + ny * amp * side); }
602+ ctx.lineTo(x1 - ux * lead, y1 - uy * lead); ctx.lineTo(x1, y1); ctx.stroke();
603+ }
604+ function drawChain() {
605+ const { ctx, W, H } = setup($("cChain"));
606+ if (!params) return;
607+ const N = params.N, x0 = 26, x1 = W - 26, cy = H / 2 + 6;
608+ const gap = (x1 - x0) / (N + 1);
609+ const ink = css("--ink"), muted = css("--muted");
610+ // walls
611+ ctx.strokeStyle = muted; ctx.lineWidth = 2;
612+ ctx.beginPath(); ctx.moveTo(x0, 16); ctx.lineTo(x0, H - 16); ctx.moveTo(x1, 16); ctx.lineTo(x1, H - 16); ctx.stroke();
613+ ctx.lineWidth = 1; ctx.strokeStyle = css("--rule");
614+ for (let y = 16; y < H - 16; y += 8) { ctx.beginPath(); ctx.moveTo(x0, y); ctx.lineTo(x0 - 6, y + 6); ctx.moveTo(x1, y); ctx.lineTo(x1 + 6, y - 6); ctx.stroke(); }
615+ // displacement
616+ const u = new Float64Array(N);
617+ let umax = 0;
618+ if (A && A.snap) {
619+ const nS = A.snap.length / N;
620+ for (let i = 0; i < nS * N; i += 3) umax = Math.max(umax, Math.abs(A.snap[i]));
621+ if (smOn) { const i = Math.min(nS - 1, Math.floor(smIdx)); for (let j = 0; j < N; j++) u[j] = A.snap[i * N + j]; }
622+ }
623+ const sc = umax > 0 ? (H / 2 - 26) / umax : 1;
624+ const px = (j) => x0 + (j + 1) * gap, py = (j) => cy - u[j] * sc;
625+ // springs, thickness by stiffness 1 + 3 beta Delta^2
626+ ctx.strokeStyle = css("--x"); ctx.lineJoin = "round";
627+ for (let j = -1; j < N; j++) {
628+ const ua = j >= 0 ? u[j] : 0, ub = j < N - 1 ? u[j + 1] : 0;
629+ const d = ub - ua, kst = 1 + 3 * d * d;
630+ ctx.lineWidth = Math.min(5, 0.8 + 1.4 * Math.log(kst) + 0.6);
631+ const xa = j >= 0 ? px(j) : x0, ya = j >= 0 ? py(j) : cy, xb = j < N - 1 ? px(j + 1) : x1, yb = j < N - 1 ? py(j + 1) : cy;
632+ zigzag(ctx, xa, ya, xb, yb, Math.max(2, Math.round(gap / 7)), Math.min(4, gap / 6));
633+ }
634+ // masses
635+ const r = Math.max(3, Math.min(8, gap * 0.28));
636+ for (let j = 0; j < N; j++) {
637+ ctx.fillStyle = ink; ctx.beginPath(); ctx.arc(px(j), py(j), r, 0, 2 * Math.PI); ctx.fill();
638+ }
639+ // pickup ring
640+ const pk = params.pickup - 1;
641+ ctx.strokeStyle = css("--y"); ctx.lineWidth = 2;
642+ ctx.beginPath(); ctx.arc(px(pk), py(pk), r + 5, 0, 2 * Math.PI); ctx.stroke();
643+ ctx.fillStyle = css("--y"); ctx.font = "11px " + css("--mono"); ctx.textAlign = "center"; ctx.textBaseline = "top";
644+ ctx.fillText("pickup", px(pk), H - 18);
645+ // strike arrows (velocity profile) when at rest
646+ if (!smOn) {
647+ ctx.strokeStyle = css("--chaos"); ctx.fillStyle = css("--chaos"); ctx.lineWidth = 2.5;
648+ for (let j = 0; j < N; j++) {
649+ const w = Math.exp(-0.5 * ((j + 1 - params.center) / params.width) ** 2);
650+ if (w < 0.08) continue;
651+ const L = 34 * w, X = px(j), Y0 = cy - r - 4;
652+ ctx.beginPath(); ctx.moveTo(X, Y0); ctx.lineTo(X, Y0 - L); ctx.stroke();
653+ ctx.beginPath(); ctx.moveTo(X, Y0 - L - 6); ctx.lineTo(X - 4, Y0 - L + 1); ctx.lineTo(X + 4, Y0 - L + 1); ctx.closePath(); ctx.fill();
654+ }
655+ ctx.fillStyle = ink; ctx.font = "600 12px " + css("--sans"); ctx.textAlign = "left"; ctx.textBaseline = "top";
656+ ctx.fillText("strike, E = " + (params.E0 < 10 ? fmt(params.E0, 2) : fmt(params.E0, 1)), px(Math.max(0, params.center - 1)) + 12, 10);
657+ } else if (A) {
658+ ctx.fillStyle = muted; ctx.font = "11px " + css("--mono"); ctx.textAlign = "left"; ctx.textBaseline = "top";
659+ ctx.fillText("t = " + fmt(smIdx / SR * 1000, 2) + " ms of sound (displacement drawn sideways, magnified)", x0 + 10, 8);
660+ }
661+ }
662+ function smTick(now) {
663+ if (!smOn || !A) return;
664+ const dtReal = smLast ? Math.min(0.1, (now - smLast) / 1000) : 0;
665+ smLast = now;
666+ smIdx += 2.5 * dtReal / params.dtOut; // 2.5 time units per second
667+ if (smIdx >= A.snap.length / params.N - 1) { stopSlow(); return; }
668+ drawChain();
669+ smRaf = requestAnimationFrame(smTick);
670+ }
671+ function startSlow() { if (!A) return; smOn = true; smIdx = 0; smLast = 0; $("slow").setAttribute("aria-pressed", "true"); $("slow").textContent = "Stop slow motion"; smRaf = requestAnimationFrame(smTick); }
672+ function stopSlow(silent) { smOn = false; cancelAnimationFrame(smRaf); $("slow").setAttribute("aria-pressed", "false"); $("slow").textContent = "Slow motion"; if (!silent) drawChain(); }
673+
674+ // ---------------- plots ----------------
675+ function drawModes() {
676+ const { ctx, W, H } = setup($("cModes"));
677+ const ax = axesTime(ctx, W, H, 0, 1, [], () => "", params.dur);
678+ if (!A || !A.shares) return;
679+ const N = params.N, nD = A.neff.length;
680+ const off = document.createElement("canvas");
681+ off.width = nD; off.height = N;
682+ const octx = off.getContext("2d"), img = octx.createImageData(nD, N);
683+ const bg = hexToRgb(css("--panel")), ink = hexToRgb(css("--ink")), hot = hexToRgb(css("--chaos"));
684+ for (let i = 0; i < nD; i++) for (let k = 0; k < N; k++) {
685+ const v = Math.min(1, Math.max(0, (Math.log10(Math.max(A.shares[i * N + k], 1e-12)) + 6) / 6));
686+ const c = heatColor(v, bg, ink, hot);
687+ const o = ((N - 1 - k) * nD + i) * 4;
688+ img.data[o] = c[0]; img.data[o + 1] = c[1]; img.data[o + 2] = c[2]; img.data[o + 3] = 255;
689+ }
690+ octx.putImageData(img, 0, 0);
691+ ctx.imageSmoothingEnabled = false;
692+ ctx.drawImage(off, ax.x0, ax.y0, ax.x1 - ax.x0, ax.y1 - ax.y0);
693+ ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono"); ctx.textAlign = "right"; ctx.textBaseline = "middle";
694+ const step = N <= 16 ? 4 : N <= 40 ? 8 : 16;
695+ for (let k = step; k <= N; k += step) ctx.fillText("k " + k, ax.x0 - 6, ax.y1 - (k - 0.5) / N * (ax.y1 - ax.y0));
696+ ctx.fillText("k 1", ax.x0 - 6, ax.y1 - 0.5 / N * (ax.y1 - ax.y0));
697+ const end = chaosEnd(A.lam);
698+ if (end !== null) { ctx.strokeStyle = css("--chaos"); ctx.setLineDash([5, 4]); ctx.lineWidth = 1.5; ctx.beginPath(); ctx.moveTo(ax.X(end), ax.y0); ctx.lineTo(ax.X(end), ax.y1); ctx.stroke(); ctx.setLineDash([]); }
699+ }
700+ function drawNeff() {
701+ const { ctx, W, H } = setup($("cNeff"));
702+ const ax = axesTime(ctx, W, H, 0, 1, [0, 0.25, 0.5, 0.75, 1], (v) => fmt(v, 2), params.dur);
703+ chaosSpan(ctx, ax);
704+ for (let k = 2; k < 8; k++) if (members[k]) seriesTime(ctx, ax, members[k].neff, css("--faint"), 1);
705+ if (B) seriesTime(ctx, ax, B.neff, css("--y"), 1.4);
706+ if (A) seriesTime(ctx, ax, A.neff, css("--x"), 1.6);
707+ const S = strikeAt(sel);
708+ if (S) seriesTime(ctx, ax, S.neff, strikeColor(sel), 3);
709+ }
710+ function drawFinal() {
711+ const { ctx, W, H } = setup($("cFinal"));
712+ if (!params) return;
713+ const N = params.N, x0 = PAD.l, x1 = W - PAD.r, y0 = PAD.t, y1 = H - PAD.b;
714+ const lo = -8, hi = 0;
715+ const X = (k) => x0 + (k - 1) / Math.max(1, N - 1) * (x1 - x0), Y = (v) => y1 - (v - lo) / (hi - lo) * (y1 - y0);
716+ ctx.strokeStyle = css("--grid"); ctx.lineWidth = 1; ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono");
717+ ctx.textAlign = "right"; ctx.textBaseline = "middle";
718+ for (let v = lo; v <= hi; v += 2) { ctx.beginPath(); ctx.moveTo(x0, Y(v)); ctx.lineTo(x1, Y(v)); ctx.stroke(); ctx.fillText("1e" + v, x0 - 6, Y(v)); }
719+ ctx.textAlign = "center"; ctx.textBaseline = "top";
720+ const step = N <= 16 ? 4 : N <= 40 ? 8 : 16;
721+ for (let k = 1; k <= N; k += (k === 1 ? step - 1 : step)) ctx.fillText("k " + k, X(k), y1 + 6);
722+ const line = (s, color, w) => {
723+ if (!s) return;
724+ ctx.strokeStyle = color; ctx.lineWidth = w; ctx.beginPath();
725+ for (let k = 1; k <= N; k++) { const v = Math.max(lo, Math.log10(Math.max(s[k - 1], 1e-12))); if (k === 1) ctx.moveTo(X(k), Y(v)); else ctx.lineTo(X(k), Y(v)); }
726+ ctx.stroke();
727+ };
728+ for (let k = 2; k < 8; k++) if (members[k]) line(members[k].finalShare, css("--faint"), 1);
729+ if (B) line(B.finalShare, css("--y"), 1.4);
730+ if (A) line(A.finalShare, css("--x"), 1.6);
731+ const S = strikeAt(sel);
732+ if (S) line(S.finalShare, strikeColor(sel), 3);
733+ }
734+ function drawLam() {
735+ const { ctx, W, H } = setup($("cLam"));
736+ if (!A || !A.lam) { axesTime(ctx, W, H, 0, 1, [0, 1], () => "", params.dur); return; }
737+ let mx = 0; for (const v of A.lam) if (v > mx) mx = v;
738+ mx = Math.max(mx * 1.1, 0.02 / params.tau);
739+ const st = niceStep(mx / 3), ticks = [];
740+ for (let v = 0; v <= mx; v += st) ticks.push(v);
741+ const ax = axesTime(ctx, W, H, -0.15 * mx, mx, ticks, (v) => fmt(v, 0), params.dur);
742+ chaosSpan(ctx, ax);
743+ ctx.strokeStyle = css("--rule"); ctx.beginPath(); ctx.moveTo(ax.x0, ax.Y(0)); ctx.lineTo(ax.x1, ax.Y(0)); ctx.stroke();
744+ const e = A.energy, e0 = e[0] || 1, eDb = new Float32Array(e.length);
745+ for (let i = 0; i < e.length; i++) eDb[i] = -0.15 * mx + (Math.max(10 * Math.log10(e[i] / e0), -80) + 80) / 80 * 1.15 * mx;
746+ seriesTime(ctx, ax, eDb, css("--muted"), 1, [4, 4]);
747+ seriesTime(ctx, ax, A.lam, css("--x"), 1.6);
748+ }
749+
750+ // ---------------- spectrogram ----------------
751+ let specFmax = 0;
752+ function fft(re, im) {
753+ const n = re.length;
754+ for (let i = 1, j = 0; i < n; i++) { let bit = n >> 1; for (; j & bit; bit >>= 1) j ^= bit; j ^= bit; if (i < j) { [re[i], re[j]] = [re[j], re[i]]; [im[i], im[j]] = [im[j], im[i]]; } }
755+ for (let len = 2; len <= n; len <<= 1) {
756+ const ang = -2 * Math.PI / len, wr = Math.cos(ang), wi = Math.sin(ang);
757+ for (let i = 0; i < n; i += len) { let cr = 1, ci = 0; for (let k = 0; k < len / 2; k++) { const a = i + k, b = a + len / 2; const tr = re[b] * cr - im[b] * ci, ti = re[b] * ci + im[b] * cr; re[b] = re[a] - tr; im[b] = im[a] - ti; re[a] += tr; im[a] += ti; const ncr = cr * wr - ci * wi; ci = cr * wi + ci * wr; cr = ncr; } }
758+ }
759+ }
760+ function spectrogramOf(S) {
761+ if (S.spec) return S.spec;
762+ const sig = S.sig, N = 2048, hop = 256, frames = Math.max(1, Math.floor((sig.length - N) / hop));
763+ const bins = Math.floor(specFmax / (SR / N));
764+ const win = new Float32Array(N); for (let i = 0; i < N; i++) win[i] = 0.5 - 0.5 * Math.cos(2 * Math.PI * i / (N - 1));
765+ const mag = new Float32Array(frames * bins), re = new Float64Array(N), im = new Float64Array(N);
766+ let mx = -Infinity;
767+ for (let f = 0; f < frames; f++) {
768+ for (let i = 0; i < N; i++) { re[i] = sig[f * hop + i] * win[i]; im[i] = 0; }
769+ fft(re, im);
770+ for (let b = 0; b < bins; b++) { const v = 10 * Math.log10(re[b] * re[b] + im[b] * im[b] + 1e-20); mag[f * bins + b] = v; if (v > mx) mx = v; }
771+ }
772+ S.spec = { mag, frames, bins, mx };
773+ return S.spec;
774+ }
775+ function drawSpec() {
776+ const { ctx, W, H } = setup($("cSpec"));
777+ const ax = axesTime(ctx, W, H, 0, 1, [], () => "", params.dur);
778+ const S = strikeAt(sel) || A;
779+ if (!S || !S.sig) return;
780+ const ref = spectrogramOf(A).mx;
781+ const { mag, frames, bins } = spectrogramOf(S);
782+ const mx = ref;
783+ const off = document.createElement("canvas"); off.width = frames; off.height = bins;
784+ const octx = off.getContext("2d"), img = octx.createImageData(frames, bins);
785+ const bg = hexToRgb(css("--panel")), ink = hexToRgb(css("--ink")), hot = hexToRgb(css("--chaos"));
786+ for (let f = 0; f < frames; f++) for (let b = 0; b < bins; b++) {
787+ const v = Math.min(1, Math.max(0, (mag[f * bins + b] - (mx - 70)) / 70));
788+ const c = heatColor(v, bg, ink, hot), o = ((bins - 1 - b) * frames + f) * 4;
789+ img.data[o] = c[0]; img.data[o + 1] = c[1]; img.data[o + 2] = c[2]; img.data[o + 3] = 255;
790+ }
791+ octx.putImageData(img, 0, 0);
792+ const t0 = 1024 / SR, t1 = (frames * 256 + 1024) / SR;
793+ ctx.drawImage(off, ax.X(t0), ax.y0, ax.X(t1) - ax.X(t0), ax.y1 - ax.y0);
794+ ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono"); ctx.textAlign = "right"; ctx.textBaseline = "middle";
795+ const step = niceStep(specFmax / 4);
796+ for (let f = 0; f <= specFmax; f += step) ctx.fillText(f >= 1000 ? fmt(f / 1000, f % 1000 ? 1 : 0) + "k" : String(f), ax.x0 - 6, ax.y1 - f / specFmax * (ax.y1 - ax.y0));
797+ ctx.fillStyle = strikeColor(sel); ctx.font = "600 12px " + css("--sans"); ctx.textAlign = "right"; ctx.textBaseline = "top";
798+ ctx.fillText("strike " + strikeName(sel), ax.x1 - 6, ax.y0 + 4);
799+ }
800+
801+ function drawAll() {
802+ if (!params) return;
803+ drawChain(); drawModes(); drawNeff(); drawFinal(); drawLam(); drawSpec(); drawOverlays(null, null);
804+ }
805+
806+ // ---------------- audio ----------------
807+ let actx = null, bufs = [], src = null, playStart = 0, playing = null, raf = 0;
808+ function scaleAudio(s) {
809+ const fade = Math.round(0.01 * SR), o = new Float32Array(s.length);
810+ for (let i = 0; i < s.length; i++) { const tail = s.length - 1 - i; o[i] = Math.max(-1, Math.min(1, s[i] * audioGain * (tail < fade ? tail / fade : 1))); }
811+ return o;
812+ }
813+ function prepareAudio() {
814+ let m = 0;
815+ for (const s of [A.sig, B.sig]) for (let i = 0; i < s.length; i++) m = Math.max(m, Math.abs(s[i]));
816+ audioGain = 0.8 / (m || 1); // one gain for every strike, so loudness differences are real
817+ A.audio = scaleAudio(A.sig); B.audio = scaleAudio(B.sig); bufs = [];
818+ specFmax = Math.min(16000, Math.max(1500, 1.6 * params.f1 * params.omN / params.om1));
819+ }
820+ function ensureCtx() { if (!actx) actx = new (window.AudioContext || window.webkitAudioContext)(); if (actx.state === "suspended") actx.resume(); }
821+ function toBuffer(d) { const b = actx.createBuffer(1, d.length, SR); b.copyToChannel(d, 0); return b; }
822+ function stopAudio() {
823+ if (src) { try { src.onended = null; src.stop(); } catch (e) {} src = null; }
824+ playing = null; cancelAnimationFrame(raf);
825+ for (const id of ["playA", "playB", "playAB", "playSel"]) if ($(id)) $(id).setAttribute("aria-pressed", "false");
826+ drawOverlays(null, null);
827+ }
828+ function play(k, then) {
829+ const S = strikeAt(k);
830+ if (!S || !S.audio) return;
831+ ensureCtx(); stopAudio();
832+ if (!bufs[k]) bufs[k] = toBuffer(S.audio);
833+ src = actx.createBufferSource(); src.buffer = bufs[k]; src.connect(actx.destination);
834+ src.onended = () => { src = null; if (then !== null && then !== undefined) setTimeout(() => play(then, null), 350); else stopAudio(); };
835+ playStart = actx.currentTime; playing = k; src.start();
836+ if (k === 0) $("playA").setAttribute("aria-pressed", "true");
837+ else if (k === 1) $("playB").setAttribute("aria-pressed", "true");
838+ if (k === sel && $("playSel")) $("playSel").setAttribute("aria-pressed", "true");
839+ if (then !== null && then !== undefined) $("playAB").setAttribute("aria-pressed", "true");
840+ tick();
841+ }
842+ function tick() { if (playing === null) return; drawOverlays(actx.currentTime - playStart, playing); raf = requestAnimationFrame(tick); }
843+ function drawOverlays(t, which) {
844+ for (const [id, base] of [["cModesO", "cModes"], ["cNeffO", "cNeff"], ["cLamO", "cLam"], ["cSpecO", "cSpec"]]) {
845+ const { ctx, W } = setup($(id));
846+ if (t === null) continue;
847+ const H = $(base).getBoundingClientRect().height;
848+ const x = PAD.l + Math.min(t, params.dur) / params.dur * (W - PAD.r - PAD.l);
849+ ctx.strokeStyle = strikeColor(which); ctx.lineWidth = 2;
850+ ctx.beginPath(); ctx.moveTo(x, PAD.t); ctx.lineTo(x, H - PAD.b); ctx.stroke();
851+ }
852+ }
853+
854+ // ---------------- wiring ----------------
855+ $("playA").addEventListener("click", () => play(0, null));
856+ $("playB").addEventListener("click", () => play(1, null));
857+ $("playAB").addEventListener("click", () => play(0, 1));
858+ $("stop").addEventListener("click", stopAudio);
859+ $("slow").addEventListener("click", () => (smOn ? stopSlow() : startSlow()));
860+ let timer = 0;
861+ for (const id of ["E0", "N", "mu", "pos", "wid", "dl", "tlo", "thi", "f1", "dur"]) {
862+ $(id).addEventListener("input", () => { readParams(); clearTimeout(timer); timer = setTimeout(run, 300); });
863+ }
864+ document.querySelectorAll("[data-preset]").forEach((b) => b.addEventListener("click", () => { setControls(PRESETS[b.dataset.preset]); run(); }));
865+ let rt = 0;
866+ window.addEventListener("resize", () => { clearTimeout(rt); rt = setTimeout(drawAll, 120); });
867+ try {
868+ window.matchMedia("(prefers-color-scheme: dark)").addEventListener("change", drawAll);
869+ new MutationObserver(drawAll).observe(document.documentElement, { attributes: true, attributeFilter: ["data-theme"] });
870+ } catch (e) {}
871+
872+ setControls(PRESETS.strong);
873+ run();
874+})();
875+</script>