Gong wire: FPUT chain with bending stiffness and a dome, experiments, web app
7 changed files+1249−0
.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+__pycache__/
2+*.npz
README.mdadded+43−0View file
@@ -0,0 +1,43 @@
1+# gong-wire
2+
3+The third chaos toy, after `../chaos-toy-examples` (two coupled oscillators) and `../fput-chain` (the struck FPUT chain). The FPUT chain cascades among its 32 nearly harmonic modes, which span only about 4 kHz when the lowest is at 100 Hz; a gong or a cymbal spreads energy upward through a dense, inharmonic spectrum reaching well above 10 kHz, over tens of milliseconds. This project adds the two ingredients that make the difference, bending stiffness and curvature, to the FPUT chain.
4+
5+## The model (`wire.py`)
6+
7+$N$ unit beads on a wire between two pinned ends move sideways, $y_j$. The wire may be built as a shallow arch of height $H_d$, rest shape $z_j = 4 H_d x(1-x)$. The energy is
8+
9+$$\sum_j \left[\tfrac12 T D_j^2 + \beta \left(z'_j D_j + \tfrac12 D_j^2\right)^2\right] + \tfrac12 K_b \sum_j (y_{j+1} - 2y_j + y_{j-1})^2, \qquad D_j = y_{j+1} - y_j,$$
10+
11+with $z'_j$ the slope of link $j$ at rest. Each term has a role:
12+
13+- *Stretching* ($\beta$). A sideways offset $D$ between neighbours stretches the link by about $D^2/2$. With $H_d = K_b = 0$ the model is exactly the FPUT-$\beta$ chain, read as the sideways motion of a beaded string; this is what passes energy between modes, as the midplane stretching does in a plate.
14+- *Bending* ($K_b$). Resistance to kinks makes the frequencies grow like the square of the wavenumber at short wavelengths, as in a plate: with $K_b = 4$ and 64 beads the modes span 100 Hz to about 17 kHz instead of 100 Hz to 4 kHz.
15+- *Dome* ($H_d$). With a curved rest shape the link stretch has a part linear in $D$, which stiffens the wire (as curvature stiffens a shell) and adds a quadratic coupling ($z' D^3$ in the energy), through which two modes can drive a third at the sum of their frequencies, the three-wave coupling of gongs and cymbals.
16+
17+With a dome the normal modes are no longer sines and are computed numerically. Damping is friction plus dashpots between neighbours; with bending, the dashpot decay rate grows roughly in proportion to frequency, as in a metal plate. `run` uses leapfrog time stepping with a Crank-Nicolson step for the damping. Checks: the force equals minus the energy gradient to 1e-10, the stiffness matrix matches the linearized force to 1e-7, and energy is conserved to about 1e-4 in the exploration runs.
18+
19+## Results (2026-10-05)
20+
21+**Exploration, undamped** (`experiments.py explore`, `results/explore.png`; 64 beads, strike centered on bead 14 with width 3, lowest mode at 100 Hz). Share of the energy in modes above 4 kHz at 10, 50 and 300 ms:
22+
23+| bending | dome | E = 10 | E = 50 |
24+| --- | --- | --- | --- |
25+| 0 | 0 | 3%, 20%, 19% (top mode 4.1 kHz) | 12%, 12%, 11% |
26+| 4 | 0 | 0.6%, 6%, 16% | 15%, 48%, 74% |
27+| 4 | 3 | 1%, 9%, 21% | 24%, 44%, 63% |
28+| 4 | 8 | 1%, 8%, 22% | 15%, 42%, 61% |
29+
30+At $E = 2$ almost nothing reaches 4 kHz in any configuration. With bending, a strong strike's energy climbs into the 4 to 17 kHz band over the first 50 to 100 ms: a bloom rather than the instant mixing of the plain chain, which has no high band to climb into. The dome speeds up the early spreading (at $E = 10$ the share above 1 kHz after 10 ms rises from 24% without a dome to 34 to 43% with one), but its effect is smaller than that of bending.
31+
32+**Bloom without chaos.** In the web app (damped: lowest mode ringing 3 s, highest 0.3 s) the cascade and the chaos separate clearly. At $E = 20$ about 10% of the energy reaches above 4 kHz within 20 ms, yet the local Lyapunov exponent stays near zero and strikes differing by $10^{-6}$ end identical. At $E = 60$ a perturbation grows only by about $e^{4.6}$ within 200 ms; at $E = 200$, 59% of the energy reaches above 4 kHz by 35 ms, the motion is chaotic for the first 0.11 s (peak rate 334 per second, total growth about $e^{14}$), and 8 strikes differing by $10^{-6}$ end with spectra that differ by 4.5 dB. As in the FPUT chain, a cascade is not the same as sensitivity: the bloom is reproducible from strike to strike over a wide range of strengths, and only very hard strikes leave a different spectrum each time. Bending detunes the modes, which presumably raises the threshold for chaos much as the on-site springs did in the FPUT chain.
33+
34+## Web application
35+
36+`web/index.html` (single page; published as a private artifact at https://claude.ai/artifact/RqEkYpuX4wqHXT7Rcfs3ZU) runs 8 strikes in parallel Web Workers. Controls: strike energy, bending stiffness, dome height, number of beads, strike position and width, the difference between strikes, ring times of the lowest and highest modes, the pitch of the lowest mode (held so the highest mode stays below 20 kHz), and length. Panels: the domed wire with its beads (with slow motion), strike A's energy per mode labelled by pitch, the share of each strike's energy above 1 and 4 kHz, the final spectra against pitch, the local Lyapunov exponent with the energy, and a spectrogram of any strike (selector or arrow keys). Presets: Soft strike ($E = 3$), Gong ($E = 20$), Hard strike ($E = 200$), String (FPUT, no bending or dome). Locally: `cd web && python -m http.server`.
37+
38+## Running
39+
40+Experiments run on the office workstation (`~/work/gong-wire`, Python 3.9 with numpy and scipy), plots locally:
41+
42+ WORKERS=18 python experiments.py explore # results/explore.npz
43+ python plots.py # results/explore.png
experiments.pyadded+68−0View file
@@ -0,0 +1,68 @@
1+"""Exploration of the struck wire (wire.py). Writes results/<name>.npz; plots.py draws them.
2+
3+ python experiments.py explore # undamped: bending x dome x strike energy; high-band share over time
4+
5+Time is reported in seconds of sound with the lowest mode at F_LOW Hz, as in the web app.
6+"""
7+
8+from __future__ import annotations
9+
10+import os
11+import sys
12+from multiprocessing import Pool
13+
14+os.environ.setdefault("OMP_NUM_THREADS", "1")
15+
16+import numpy as np # noqa: E402
17+
18+from wire import Wire, run, strike # noqa: E402
19+
20+HERE = os.path.dirname(os.path.abspath(__file__))
21+RES = os.path.join(HERE, "results")
22+N, CENTER, WIDTH, F_LOW = 64, 14, 3.0, 100.0
23+
24+
25+def dt_for(w: Wire, E):
26+ om, _ = w.modes()
27+ return 0.4 / np.sqrt(om[-1] ** 2 + 12 * E)
28+
29+
30+def tau_for(w: Wire):
31+ om, _ = w.modes()
32+ return om[0] / (2 * np.pi * F_LOW) # seconds of sound per time unit
33+
34+
35+def _explore_one(args):
36+ Kb, Hd, E = args
37+ w = Wire(N=N, Kb=Kb, Hd=Hd)
38+ tau = tau_for(w)
39+ dt = dt_for(w, E)
40+ t_end = 0.3 / tau
41+ rec = max(1, int(round(0.002 / tau / dt)))
42+ y0, v0 = strike(w, E, CENTER, WIDTH, M=2, scale=[1.0, 1 + 1e-10])
43+ r = run(w, y0, v0, t_end, dt, record_every=rec)
44+ f = r["omega"] / (2 * np.pi * tau) # mode frequencies in Hz of sound
45+ return Kb, Hd, E, r["t"] * tau, r["Ek"][:, 0, :], f, r["E"][:, 0]
46+
47+
48+def explore(pool):
49+ jobs = [(Kb, Hd, E) for Kb in (0.0, 4.0) for Hd in (0.0, 3.0, 8.0) for E in (2.0, 10.0, 50.0)]
50+ out = {}
51+ for Kb, Hd, E, t, Ek, f, Etot in pool.imap_unordered(_explore_one, jobs):
52+ k = f"Kb{Kb:g}_Hd{Hd:g}_E{E:g}"
53+ out[k + "_t"], out[k + "_Ek"], out[k + "_f"], out[k + "_E"] = t, Ek, f, Etot
54+ share = Ek / Ek.sum(1, keepdims=True)
55+ hi1, hi4 = share[:, f > 1000].sum(1), share[:, f > 4000].sum(1)
56+ print(f" Kb {Kb:g} Hd {Hd:g} E {E:g}: top {f[-1]:.0f} Hz; share above 1 kHz at 10/50/300 ms "
57+ f"{hi1[4]:.2f} {hi1[24]:.2f} {hi1[-1]:.2f}; above 4 kHz {hi4[4]:.3f} {hi4[24]:.3f} {hi4[-1]:.3f}; "
58+ f"energy drift {abs(Etot[-1] / Etot[0] - 1):.1e}", flush=True)
59+ np.savez(os.path.join(RES, "explore.npz"), **out)
60+
61+
62+if __name__ == "__main__":
63+ os.makedirs(RES, exist_ok=True)
64+ which = sys.argv[1:] or ["explore"]
65+ with Pool(int(os.environ.get("WORKERS", "18"))) as pool:
66+ for w in which:
67+ print(w, flush=True)
68+ dict(explore=explore)[w](pool)
plots.pyadded+46−0View file
@@ -0,0 +1,46 @@
1+"""Figure from results/explore.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+
16+
17+def explore():
18+ d = np.load(os.path.join(RES, "explore.npz"))
19+ cases = [(0.0, 0.0), (4.0, 0.0), (4.0, 3.0), (4.0, 8.0)]
20+ Es = (2.0, 10.0, 50.0)
21+ fig, ax = plt.subplots(2, 4, figsize=(20, 8.5))
22+ for j, (Kb, Hd) in enumerate(cases):
23+ k = f"Kb{Kb:g}_Hd{Hd:g}_E50"
24+ t, Ek, f = d[k + "_t"], d[k + "_Ek"], d[k + "_f"]
25+ share = Ek / Ek.sum(1, keepdims=True)
26+ ax[0, j].pcolormesh(t * 1000, f, np.log10(share.T + 1e-12), vmin=-6, vmax=0, cmap="magma", shading="auto",
27+ rasterized=True)
28+ ax[0, j].set(yscale="log", ylim=(f[0], f[-1]), title=f"bending {Kb:g}, dome {Hd:g}, E = 50",
29+ xlabel="time (ms of sound)", ylabel="mode pitch (Hz)" if j == 0 else None)
30+ for E, col in zip(Es, ("C0", "C1", "C3")):
31+ k = f"Kb{Kb:g}_Hd{Hd:g}_E{E:g}"
32+ t, Ek, f = d[k + "_t"], d[k + "_Ek"], d[k + "_f"]
33+ share = Ek / Ek.sum(1, keepdims=True)
34+ ax[1, j].plot(t * 1000, share[:, f > 4000].sum(1), color=col, label=f"E = {E:g}, above 4 kHz")
35+ ax[1, j].plot(t * 1000, share[:, f > 1000].sum(1), color=col, lw=0.8, ls="--", label=f"E = {E:g}, above 1 kHz")
36+ ax[1, j].set(ylim=(0, 1), xlabel="time (ms of sound)", ylabel="share of energy" if j == 0 else None)
37+ ax[1, j].grid(alpha=0.3)
38+ ax[1, 0].legend(fontsize=7)
39+ fig.suptitle("Undamped, 64 beads, lowest mode at 100 Hz: energy share per mode (top, log color, E = 50) and above 1 / 4 kHz (bottom)")
40+ fig.tight_layout()
41+ fig.savefig(os.path.join(RES, "explore.png"), dpi=90)
42+
43+
44+if __name__ == "__main__":
45+ explore()
46+ print("wrote results/explore.png")
results/explore.pngadded+0−0View file
Binary file not shown.
web/index.htmladded+922−0View file
@@ -0,0 +1,922 @@
1+<meta charset="utf-8">
2+<title>Gong Wire</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+ --brass: #b07a22;
18+ --chaos-wash: rgba(232, 185, 49, 0.18);
19+ --faint: rgba(23, 35, 52, 0.16);
20+ --focus: #2563a6;
21+ --sans: "Archivo", "Helvetica Neue", Arial, sans-serif;
22+ --mono: "IBM Plex Mono", ui-monospace, "SFMono-Regular", Menlo, monospace;
23+}
24+@media (prefers-color-scheme: dark) {
25+ :root:not([data-theme="light"]) {
26+ color-scheme: dark;
27+ --paper: #0f161f; --panel: #16202c; --ink: #e3e9f0; --muted: #93a3b5; --rule: #2a3747; --grid: #1f2a37;
28+ --x: #6fa8e8; --y: #d57bae; --chaos: #e8b931; --brass: #d8a64f; --chaos-wash: rgba(232, 185, 49, 0.14);
29+ --faint: rgba(227, 233, 240, 0.16); --focus: #6fa8e8;
30+ }
31+}
32+:root[data-theme="dark"] {
33+ color-scheme: dark;
34+ --paper: #0f161f; --panel: #16202c; --ink: #e3e9f0; --muted: #93a3b5; --rule: #2a3747; --grid: #1f2a37;
35+ --x: #6fa8e8; --y: #d57bae; --chaos: #e8b931; --brass: #d8a64f; --chaos-wash: rgba(232, 185, 49, 0.14);
36+ --faint: rgba(227, 233, 240, 0.16); --focus: #6fa8e8;
37+}
38+* { box-sizing: border-box; }
39+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; }
40+.wrap { max-width: 1280px; margin: 0 auto; display: grid; gap: 20px; }
41+header { display: grid; gap: 6px; }
42+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; }
43+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%; }
44+.lede { max-width: 76ch; color: var(--muted); margin: 0; }
45+.layout { display: grid; grid-template-columns: 300px minmax(0, 1fr); gap: 20px; align-items: start; }
46+@media (max-width: 900px) { .layout { grid-template-columns: 1fr; } }
47+aside { background: var(--panel); border: 1px solid var(--rule); border-radius: 6px; padding: 16px; display: grid; grid-template-columns: minmax(0, 1fr); gap: 16px; position: sticky; top: calc(env(safe-area-inset-top, 0px) + 12px); }
48+@media (max-width: 900px) { aside { position: static; } }
49+.label { font-size: 11px; font-weight: 650; letter-spacing: 0.08em; text-transform: uppercase; color: var(--muted); }
50+.presets { display: flex; flex-wrap: wrap; gap: 6px; }
51+.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; }
52+.chip:hover { border-color: var(--muted); }
53+.chip[aria-pressed="true"] { background: var(--ink); color: var(--paper); border-color: var(--ink); }
54+input[type="range"]:disabled { opacity: 0.4; }
55+.field { display: grid; gap: 4px; }
56+.field .row { display: flex; justify-content: space-between; align-items: baseline; gap: 8px; }
57+.field label { font-size: 14px; font-weight: 550; }
58+.val { font-family: var(--mono); font-size: 13px; font-variant-numeric: tabular-nums; color: var(--ink); white-space: nowrap; }
59+.hint { font-size: 12px; color: var(--muted); margin: 0; }
60+input[type="range"] { width: 100%; accent-color: var(--x); }
61+.player { display: grid; gap: 8px; }
62+.btns { display: grid; grid-template-columns: 1fr 1fr; gap: 8px; }
63+.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; }
64+.btn:hover { border-color: var(--muted); }
65+.btn .sw { width: 10px; height: 10px; border-radius: 2px; flex: none; }
66+.btn[aria-pressed="true"] { background: var(--ink); color: var(--paper); border-color: var(--ink); }
67+button:focus-visible, input:focus-visible { outline: 2px solid var(--focus); outline-offset: 2px; }
68+.stats { display: grid; grid-template-columns: auto 1fr; gap: 4px 12px; font-size: 13px; margin: 0; }
69+.stats dt { color: var(--muted); }
70+.stats dd { margin: 0; font-family: var(--mono); font-variant-numeric: tabular-nums; text-align: right; }
71+.status { font-size: 12px; color: var(--muted); min-height: 1.5em; font-family: var(--mono); }
72+main { display: grid; gap: 16px; min-width: 0; }
73+.pair { display: grid; grid-template-columns: 1fr 1fr; gap: 16px; }
74+@media (max-width: 640px) { .pair { grid-template-columns: 1fr; } }
75+figure { margin: 0; display: grid; gap: 6px; min-width: 0; }
76+figcaption { display: grid; gap: 2px; }
77+figcaption b { font-weight: 650; font-size: 15px; }
78+figcaption span { font-size: 13px; color: var(--muted); max-width: 82ch; }
79+.plot { position: relative; width: 100%; background: var(--panel); border: 1px solid var(--rule); border-radius: 6px; }
80+.plot canvas { position: absolute; inset: 0; width: 100%; height: 100%; display: block; }
81+.chainplot { aspect-ratio: 3.4 / 1; min-height: 190px; }
82+.wide { aspect-ratio: 3.2 / 1; min-height: 170px; }
83+.half { aspect-ratio: 1.6 / 1; min-height: 170px; }
84+.short { aspect-ratio: 4.6 / 1; min-height: 140px; }
85+.setuptext { display: grid; gap: 10px; font-size: 14px; grid-template-columns: minmax(0, 1fr) auto; align-items: start; }
86+@media (max-width: 640px) { .setuptext { grid-template-columns: 1fr; } }
87+.setuptext p { margin: 0; max-width: 84ch; }
88+.strikes { display: flex; flex-wrap: wrap; gap: 6px; align-items: center; }
89+.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; }
90+.sk .sw { width: 8px; height: 8px; border-radius: 2px; }
91+.sk[aria-pressed="true"] { background: var(--ink); color: var(--paper); border-color: var(--ink); }
92+.sk:disabled { opacity: 0.4; cursor: default; }
93+.strikes .btn { padding: 4px 12px; font-size: 13px; margin-left: auto; }
94+.legend { display: flex; flex-wrap: wrap; gap: 14px; font-size: 12px; color: var(--muted); }
95+.legend i { display: inline-block; width: 18px; height: 3px; vertical-align: middle; margin-right: 6px; border-radius: 2px; }
96+footer { color: var(--muted); font-size: 13px; max-width: 82ch; }
97+footer p { margin: 0 0 8px; }
98+/* keep the control column inside its panel whatever the readouts say */
99+aside > * { min-width: 0; }
100+aside .field .row { flex-wrap: wrap; }
101+aside .val { white-space: normal; text-align: right; overflow-wrap: anywhere; }
102+</style>
103+
104+<div class="wrap">
105+ <header>
106+ <h1>Gong Wire
107+ <span class="eq">beads on a stiff, domed wire: stretching + bending + curvature, struck once</span>
108+ </h1>
109+ <p class="lede">The FPUT chain made into a small gong. Beads on a wire move sideways; the wire resists stretching (which couples the modes), resists bending (which spreads the modes from about 100 Hz to 17 kHz, as in a plate), and can be built as a dome (which adds the three-wave coupling of curved gongs and cymbals). A soft strike rings in its low modes. A hard one sends energy climbing into the high frequencies over tens of milliseconds, the bloom of a gong or a crash, until damping, strongest at high frequencies, pulls it back down.</p>
110+ </header>
111+
112+ <div class="layout">
113+ <aside aria-label="Controls">
114+ <div class="player">
115+ <div class="label">Listen</div>
116+ <div class="btns">
117+ <button class="btn" id="playA" type="button"><span class="sw" style="background: var(--x)"></span>Strike A</button>
118+ <button class="btn" id="playB" type="button"><span class="sw" style="background: var(--y)"></span>Strike B</button>
119+ <button class="btn" id="playAB" type="button">A, then B</button>
120+ <button class="btn" id="stop" type="button">Stop</button>
121+ </div>
122+ <div class="status" id="status">Computing…</div>
123+ </div>
124+
125+ <div class="field">
126+ <div class="row"><label for="lis">Listen at</label><span class="val" id="lisv"></span></div>
127+ <div class="presets" role="group" aria-label="Listener">
128+ <button class="chip" type="button" id="lisOne" aria-pressed="true">One bead</button>
129+ <button class="chip" type="button" id="lisSum" aria-pressed="false">All beads (sum)</button>
130+ </div>
131+ <input id="lis" type="range" min="0.02" max="0.98" step="0.01" value="0.38">
132+ <p class="hint">The sound is the velocity of the bead in the ring, or the sum of all their velocities. Changing it recomputes the strikes.</p>
133+ </div>
134+
135+ <div class="field">
136+ <div class="label">Presets</div>
137+ <div class="presets">
138+ <button class="chip" type="button" data-preset="soft">Soft strike</button>
139+ <button class="chip" type="button" data-preset="gong">Gong</button>
140+ <button class="chip" type="button" data-preset="hard">Hard strike</button>
141+ <button class="chip" type="button" data-preset="string">String (FPUT)</button>
142+ </div>
143+ </div>
144+
145+ <div class="field">
146+ <div class="row"><label for="E0">Strike energy E</label><span class="val" id="E0v"></span></div>
147+ <input id="E0" type="range" min="0" max="1000" step="1">
148+ <p class="hint">With bending on, energy climbs above 4 kHz from about E ≈ 10; the strikes only become sensitive to tiny differences (chaotic) from about E ≈ 100.</p>
149+ </div>
150+ <div class="field">
151+ <div class="row"><label for="Kb">Bending stiffness</label><span class="val" id="Kbv"></span></div>
152+ <input id="Kb" type="range" min="0" max="8" step="0.1">
153+ <p class="hint">0 is the plain FPUT string (modes up to about 4 kHz); more bending spreads them to plate-like spacing.</p>
154+ </div>
155+ <div class="field">
156+ <div class="row"><label for="Hd">Dome height</label><span class="val" id="Hdv"></span></div>
157+ <input id="Hd" type="range" min="0" max="12" step="0.25">
158+ <p class="hint">Curvature: stiffer, and a quadratic coupling that starts the spreading sooner.</p>
159+ </div>
160+ <div class="field">
161+ <div class="row"><label for="N">Beads</label><span class="val" id="Nv"></span></div>
162+ <input id="N" type="range" min="24" max="96" step="1">
163+ </div>
164+ <div class="field">
165+ <div class="row"><label for="pos">Strike position</label><span class="val" id="posv"></span></div>
166+ <input id="pos" type="range" min="0.05" max="0.5" step="0.01">
167+ </div>
168+ <div class="field">
169+ <div class="row"><label for="wid">Strike width</label><span class="val" id="widv"></span></div>
170+ <input id="wid" type="range" min="1" max="8" step="0.1">
171+ </div>
172+ <div class="field">
173+ <div class="row"><label for="dl">Difference between strikes</label><span class="val" id="dlv"></span></div>
174+ <input id="dl" type="range" min="-14" max="-2" step="0.5">
175+ </div>
176+ <div class="field">
177+ <div class="row"><label for="tlo">Ring time, lowest mode</label><span class="val" id="tlov"></span></div>
178+ <input id="tlo" type="range" min="0.5" max="10" step="0.1">
179+ </div>
180+ <div class="field">
181+ <div class="row"><label for="thi">Ring time, highest mode</label><span class="val" id="thiv"></span></div>
182+ <input id="thi" type="range" min="0.05" max="5" step="0.05">
183+ </div>
184+ <div class="field">
185+ <div class="row"><label for="f1">Pitch of the lowest mode</label><span class="val" id="f1v"></span></div>
186+ <input id="f1" type="range" min="40" max="300" step="1">
187+ <p class="hint" id="f1note" hidden>Held lower than the slider so the highest mode stays under 20 kHz.</p>
188+ </div>
189+ <div class="field">
190+ <div class="row"><label for="dur">Length</label><span class="val" id="durv"></span></div>
191+ <input id="dur" type="range" min="1" max="5" step="0.5">
192+ </div>
193+
194+ <dl class="stats" aria-live="polite">
195+ <dt>Energy above 4 kHz, peak</dt><dd id="sBloom">–</dd>
196+ <dt>… reached at</dt><dd id="sBloomT">–</dd>
197+ <dt>Chaotic phase ends</dt><dd id="sChaos">–</dd>
198+ <dt>Peak separation rate</dt><dd id="sLam">–</dd>
199+ <dt>Final spectra, 8 strikes</dt><dd id="sSpread">–</dd>
200+ </dl>
201+ </aside>
202+
203+ <main>
204+ <figure>
205+ <figcaption><b>Physical setup</b><span>Beads on a wire between two pinned ends, built as a shallow dome (dashed). The beads move sideways, drawn magnified; link thickness shows how hard each link is being stretched.</span></figcaption>
206+ <div class="plot chainplot"><canvas id="cChain"></canvas></div>
207+ <div class="setuptext">
208+ <p>Three things hold each bead. The links pull it back toward its neighbours and, when stretched, pull harder: this stretching is what passes energy between modes, the same mechanism as in the FPUT chain or the midplane of a plate. Bending stiffness resists kinks, so short waves vibrate much faster than long ones, as in a plate or cymbal (frequency grows with the square of the wavenumber). The dome makes even a small sideways motion stretch the links, which stiffens the wire and adds a coupling that lets two modes feed a third at the sum of their frequencies. The strike gives a few beads a sudden velocity; the sound is the velocity of the bead in the listener ring, or the sum over all beads (Listen at).</p>
209+ <button class="btn" id="slow" type="button" aria-pressed="false">Slow motion</button>
210+ </div>
211+ </figure>
212+
213+ <figure>
214+ <figcaption><b>Where the energy is</b><span>Share of strike A's energy in each normal mode over time, log scale, labelled by the mode's pitch. A cascade shows as energy climbing toward the top.</span></figcaption>
215+ <div class="plot wide"><canvas id="cModes"></canvas><canvas id="cModesO"></canvas></div>
216+ </figure>
217+
218+ <div class="pair">
219+ <figure>
220+ <figcaption><b>How high it reaches</b><span>Share of each strike's energy in modes above 4 kHz (solid) and above 1 kHz (thin), for all 8 strikes. Shaded: the chaotic phase of strike A.</span></figcaption>
221+ <div class="plot half"><canvas id="cHi"></canvas><canvas id="cHiO"></canvas></div>
222+ </figure>
223+ <figure>
224+ <figcaption><b>The spectrum it rings with</b><span>Each strike's final share of energy per mode, against the mode's pitch.</span></figcaption>
225+ <div class="plot half"><canvas id="cFinal"></canvas></div>
226+ </figure>
227+ </div>
228+ <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>
229+
230+ <figure>
231+ <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>
232+ <div class="plot short"><canvas id="cLam"></canvas><canvas id="cLamO"></canvas></div>
233+ </figure>
234+
235+ <figure>
236+ <figcaption><b>What you hear</b><span>Spectrogram of the selected strike at the listener. 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>
237+ <div class="strikes" role="group" aria-label="Choose a strike" id="strikeSel"></div>
238+ <div class="plot wide"><canvas id="cSpec"></canvas><canvas id="cSpecO"></canvas></div>
239+ </figure>
240+
241+ <footer>
242+ <p>Computed in your browser, the 8 strikes in parallel: leapfrog (Störmer–Verlet) time stepping with the damping applied by a Crank–Nicolson step; normal modes found numerically, since the dome makes them differ from sines. Time is scaled so that the lowest mode sounds at the chosen pitch, which is held low enough that the highest mode stays below 20 kHz. Damping is friction on the beads plus dashpots between neighbours; with bending, that makes the decay rate grow roughly in proportion to frequency, as in a metal plate.</p>
243+ </footer>
244+ </main>
245+ </div>
246+</div>
247+
248+<script id="worker-src" type="text/plain">
249+// ---------------- the wire ----------------
250+function linkSlopes(N, Hd) {
251+ const zp = new Float64Array(N + 1);
252+ for (let j = 0; j <= N; j++) {
253+ const x0 = j / (N + 1), x1 = (j + 1) / (N + 1);
254+ zp[j] = 4 * Hd * (x1 * (1 - x1) - x0 * (1 - x0));
255+ }
256+ return zp;
257+}
258+function stiffness(p, zp) {
259+ const N = p.N, K = new Float64Array(N * N);
260+ for (let j = 0; j <= N; j++) {
261+ const kl = p.T + 2 * p.beta * zp[j] * zp[j], a = j - 1, b = j;
262+ if (a >= 0) K[a * N + a] += kl;
263+ if (b < N) K[b * N + b] += kl;
264+ if (a >= 0 && b < N) { K[a * N + b] -= kl; K[b * N + a] -= kl; }
265+ }
266+ if (p.Kb) {
267+ // Kb L^2, L the Dirichlet Laplacian (pinned ends)
268+ const L = (i, j) => (i === j ? -2 : Math.abs(i - j) === 1 ? 1 : 0);
269+ for (let i = 0; i < N; i++) for (let j = Math.max(0, i - 2); j <= Math.min(N - 1, i + 2); j++) {
270+ let s = 0;
271+ for (let m = Math.max(0, i - 1); m <= Math.min(N - 1, i + 1); m++) s += L(i, m) * L(m, j);
272+ K[i * N + j] += p.Kb * s;
273+ }
274+ }
275+ return K;
276+}
277+function jacobiEigen(Ain, n) {
278+ // cyclic Jacobi: returns eigenvalues (ascending) and eigenvectors as columns V[i * n + k]
279+ const A = Float64Array.from(Ain), V = new Float64Array(n * n);
280+ for (let i = 0; i < n; i++) V[i * n + i] = 1;
281+ for (let sweep = 0; sweep < 60; sweep++) {
282+ let off = 0;
283+ for (let i = 0; i < n; i++) for (let j = i + 1; j < n; j++) off += A[i * n + j] * A[i * n + j];
284+ if (off < 1e-22) break;
285+ for (let p = 0; p < n - 1; p++) for (let q = p + 1; q < n; q++) {
286+ const apq = A[p * n + q];
287+ if (Math.abs(apq) < 1e-300) continue;
288+ const app = A[p * n + p], aqq = A[q * n + q];
289+ const th = (aqq - app) / (2 * apq);
290+ const t = Math.sign(th || 1) / (Math.abs(th) + Math.sqrt(th * th + 1));
291+ const c = 1 / Math.sqrt(t * t + 1), s = t * c;
292+ for (let k = 0; k < n; k++) {
293+ const akp = A[k * n + p], akq = A[k * n + q];
294+ A[k * n + p] = c * akp - s * akq; A[k * n + q] = s * akp + c * akq;
295+ }
296+ for (let k = 0; k < n; k++) {
297+ const apk = A[p * n + k], aqk = A[q * n + k];
298+ A[p * n + k] = c * apk - s * aqk; A[q * n + k] = s * apk + c * aqk;
299+ }
300+ for (let k = 0; k < n; k++) {
301+ const vkp = V[k * n + p], vkq = V[k * n + q];
302+ V[k * n + p] = c * vkp - s * vkq; V[k * n + q] = s * vkp + c * vkq;
303+ }
304+ }
305+ }
306+ const idx = Array.from({ length: n }, (_, i) => i).sort((a, b) => A[a * n + a] - A[b * n + b]);
307+ const vals = new Float64Array(n), vecs = new Float64Array(n * n);
308+ idx.forEach((k, m) => { vals[m] = A[k * n + k]; for (let i = 0; i < n; i++) vecs[i * n + m] = V[i * n + k]; });
309+ return { vals, vecs };
310+}
311+function setupWire(p) {
312+ const zp = linkSlopes(p.N, p.Hd);
313+ const { vals, vecs } = jacobiEigen(stiffness(p, zp), p.N);
314+ const om = new Float64Array(p.N);
315+ for (let k = 0; k < p.N; k++) om[k] = Math.sqrt(Math.max(vals[k], 0));
316+ return { zp, om, V: vecs };
317+}
318+
319+function strikeRun(p, w, scale, opts) {
320+ const N = p.N, T = p.T, beta = p.beta, Kb = p.Kb, zp = w.zp;
321+ const y = new Float64Array(N), v = new Float64Array(N), f = new Float64Array(N);
322+ const dy = new Float64Array(N), dv = new Float64Array(N), g = new Float64Array(N);
323+ const D = new Float64Array(N + 1), fl = new Float64Array(N + 1), c = new Float64Array(N + 2);
324+ let n2 = 0;
325+ for (let j = 0; j < N; j++) { const s = Math.exp(-0.5 * ((j + 1 - p.center) / p.width) ** 2); v[j] = s; n2 += s * s; }
326+ const amp = Math.sqrt(2 * p.E0 / n2) * scale;
327+ for (let j = 0; j < N; j++) v[j] *= amp;
328+ if (opts.tangent) { for (let j = 0; j < N; j++) { dy[j] = Math.sin(1.7 * j + 0.3); dv[j] = Math.cos(2.3 * j + 0.1); } normalize(dy, dv); }
329+
330+ function bend(yy, out, k) {
331+ // out -= Kb L (L yy); curvature c_j for j = 1..N (pinned ends: c_0 = c_{N+1} = 0)
332+ for (let j = 0; j < N; j++) c[j + 1] = (j > 0 ? yy[j - 1] : 0) - 2 * yy[j] + (j < N - 1 ? yy[j + 1] : 0);
333+ c[0] = 0; c[N + 1] = 0;
334+ for (let j = 0; j < N; j++) out[j] -= k * (c[j] - 2 * c[j + 1] + c[j + 2]);
335+ }
336+ function force(yy, out) {
337+ for (let j = 0; j <= N; j++) {
338+ const d = (j < N ? yy[j] : 0) - (j > 0 ? yy[j - 1] : 0);
339+ const e = zp[j] * d + 0.5 * d * d;
340+ D[j] = d; fl[j] = T * d + 2 * beta * e * (zp[j] + d);
341+ }
342+ for (let j = 0; j < N; j++) out[j] = fl[j + 1] - fl[j];
343+ if (Kb) bend(yy, out, Kb);
344+ }
345+ function tforce(yy, dd, out) {
346+ // linearized force on dd at yy (uses D from the last force(yy) call)
347+ for (let j = 0; j <= N; j++) {
348+ const dd_ = (j < N ? dd[j] : 0) - (j > 0 ? dd[j - 1] : 0);
349+ const s = zp[j] + D[j];
350+ fl[j] = (T + 2 * beta * (s * s + zp[j] * D[j] + 0.5 * D[j] * D[j])) * dd_;
351+ }
352+ for (let j = 0; j < N; j++) out[j] = fl[j + 1] - fl[j];
353+ if (Kb) bend(dd, out, Kb);
354+ }
355+ function normalize(a, b) { let s = 0; for (let j = 0; j < N; j++) s += a[j] * a[j] + b[j] * b[j]; const n = Math.sqrt(s); for (let j = 0; j < N; j++) { a[j] /= n; b[j] /= n; } return n; }
356+
357+ // Crank-Nicolson damping, (I + h/2 Dm) v+ = (I - h/2 Dm) v, Dm = c0 I - c1 L: tridiagonal (Thomas)
358+ const nsub = Math.max(1, Math.ceil(p.dtOut / p.hmax)), h = p.dtOut / nsub;
359+ const a = 0.5 * h, diag = 1 + a * (p.c0 + 2 * p.c1), offd = -a * p.c1, rdiag = 1 - a * (p.c0 + 2 * p.c1), roff = a * p.c1;
360+ const cp = new Float64Array(N), rhs = new Float64Array(N);
361+ cp[0] = offd / diag;
362+ for (let j = 1; j < N; j++) cp[j] = offd / (diag - offd * cp[j - 1]);
363+ function damp(x) {
364+ for (let j = 0; j < N; j++) rhs[j] = rdiag * x[j] + roff * ((j > 0 ? x[j - 1] : 0) + (j < N - 1 ? x[j + 1] : 0));
365+ let prev = rhs[0] / diag; x[0] = prev;
366+ for (let j = 1; j < N; j++) { prev = (rhs[j] - offd * prev) / (diag - offd * cp[j - 1]); x[j] = prev; }
367+ for (let j = N - 2; j >= 0; j--) x[j] -= cp[j] * x[j + 1];
368+ }
369+
370+ const nOut = p.nOut, K = p.decim, nD = Math.floor(nOut / K);
371+ const sig = new Float32Array(nOut);
372+ const snapN = opts.snap ? Math.min(nOut, p.snapN) : 0, snap = opts.snap ? new Float32Array(snapN * N) : null;
373+ const shares = opts.modes ? new Float32Array(nD * N) : null;
374+ const hi1 = new Float32Array(nD), hi4 = new Float32Array(nD), energy = new Float32Array(nD);
375+ const logsum = opts.tangent ? new Float64Array(nD) : null;
376+ const finalShare = new Float32Array(N), es = new Float64Array(N);
377+ const pk = p.pickup - 1;
378+ let L = 0;
379+ function modeShares() {
380+ let tot = 0;
381+ for (let k = 0; k < N; k++) {
382+ let q = 0, pp = 0;
383+ for (let j = 0; j < N; j++) { const s = w.V[j * N + k]; q += s * y[j]; pp += s * v[j]; }
384+ const e = 0.5 * (pp * pp + w.om[k] * w.om[k] * q * q);
385+ es[k] = e; tot += e;
386+ }
387+ for (let k = 0; k < N; k++) es[k] /= (tot || 1);
388+ return tot;
389+ }
390+ force(y, f);
391+ if (opts.tangent) tforce(y, dy, g);
392+ for (let i = 0; i < nOut; i++) {
393+ for (let s = 0; s < nsub; s++) {
394+ for (let j = 0; j < N; j++) { v[j] += 0.5 * h * f[j]; y[j] += h * v[j]; }
395+ if (opts.tangent) for (let j = 0; j < N; j++) { dv[j] += 0.5 * h * g[j]; dy[j] += h * dv[j]; }
396+ force(y, f);
397+ for (let j = 0; j < N; j++) v[j] += 0.5 * h * f[j];
398+ if (opts.tangent) { tforce(y, dy, g); for (let j = 0; j < N; j++) dv[j] += 0.5 * h * g[j]; }
399+ damp(v);
400+ if (opts.tangent) damp(dv);
401+ }
402+ if (!Number.isFinite(v[pk])) return { blowup: true };
403+ if (opts.tangent) { L += Math.log(normalize(dy, dv)); tforce(y, dy, g); } // g must follow the rescaled dy
404+ if (p.listenSum) { let sv = 0; for (let j = 0; j < N; j++) sv += v[j]; sig[i] = sv; } else sig[i] = v[pk];
405+ if (i < snapN) for (let j = 0; j < N; j++) snap[i * N + j] = y[j];
406+ if ((i + 1) % K === 0) {
407+ const jd = (i + 1) / K - 1;
408+ if (jd < nD) {
409+ energy[jd] = modeShares();
410+ let a1 = 0, a4 = 0;
411+ for (let k = 0; k < N; k++) { if (p.fk[k] > 1000) a1 += es[k]; if (p.fk[k] > 4000) a4 += es[k]; }
412+ hi1[jd] = a1; hi4[jd] = a4;
413+ if (shares) for (let k = 0; k < N; k++) shares[jd * N + k] = es[k];
414+ if (logsum) logsum[jd] = L;
415+ }
416+ }
417+ }
418+ modeShares();
419+ for (let k = 0; k < N; k++) finalShare[k] = es[k];
420+ return { sig, snap, shares, hi1, hi4, energy, logsum, finalShare };
421+}
422+
423+onmessage = (ev) => {
424+ const { gen, p, k } = ev.data;
425+ const w = setupWire(p);
426+ let r = strikeRun(p, w, 1 + k * p.delta, { snap: k <= 1, modes: k === 0, tangent: k === 0 });
427+ if (r.blowup) { const p2 = Object.assign({}, p, { hmax: p.hmax / 3 }); r = strikeRun(p2, w, 1 + k * p.delta, { snap: k <= 1, modes: k === 0, tangent: k === 0 }); }
428+ const tr = [r.sig.buffer, r.hi1.buffer, r.hi4.buffer, r.energy.buffer, r.finalShare.buffer];
429+ if (r.snap) tr.push(r.snap.buffer);
430+ if (r.shares) tr.push(r.shares.buffer);
431+ postMessage({ gen, k, r, om: w.om }, tr);
432+};
433+</script>
434+
435+<script>
436+(() => {
437+ const SR = 44100, NS = 8;
438+ const $ = (id) => document.getElementById(id);
439+ const css = (name) => getComputedStyle(document.documentElement).getPropertyValue(name).trim();
440+ const fmt = (v, d = 2) => Number(v).toFixed(d);
441+
442+ const PRESETS = {
443+ soft: { E0: 3, Kb: 4, Hd: 3, N: 64, pos: 0.22, wid: 3, dl: -6, tlo: 3, thi: 0.3, f1: 100, dur: 3 },
444+ gong: { E0: 20, Kb: 4, Hd: 3, N: 64, pos: 0.22, wid: 3, dl: -6, tlo: 3, thi: 0.3, f1: 100, dur: 3 },
445+ hard: { E0: 200, Kb: 4, Hd: 3, N: 64, pos: 0.22, wid: 3, dl: -6, tlo: 3, thi: 0.3, f1: 100, dur: 2 },
446+ string: { E0: 20, Kb: 0, Hd: 0, N: 64, pos: 0.22, wid: 3, dl: -6, tlo: 3, thi: 0.3, f1: 100, dur: 3 },
447+ };
448+ const eFromSlider = (v) => Math.pow(10, -0.3 + 2.9 * v / 1000);
449+ const sliderFromE = (E) => Math.round((Math.log10(E) + 0.3) / 2.9 * 1000);
450+
451+ // modes of the linear wire, for the pitch mapping (also computed in the workers)
452+ function wireModes(N, Kb, Hd) {
453+ const src = $("worker-src").textContent;
454+ if (!wireModes.fn) wireModes.fn = new Function(src.replace(/onmessage =[\s\S]*$/, "") + "; return setupWire;")();
455+ return wireModes.fn({ N, T: 1, beta: 1, Kb, Hd });
456+ }
457+
458+ function readParams() {
459+ const E0 = eFromSlider(+$("E0").value), Kb = +$("Kb").value, Hd = +$("Hd").value, N = +$("N").value;
460+ const pos = +$("pos").value, wid = +$("wid").value, dl = +$("dl").value;
461+ const tlo = +$("tlo").value, thi = Math.min(+$("thi").value, tlo), dur = +$("dur").value;
462+ const w = wireModes(N, Kb, Hd), om1 = w.om[0], omN = w.om[N - 1];
463+ const f1max = 19000 / (omN / om1);
464+ const f1 = Math.min(+$("f1").value, f1max);
465+ const center = Math.max(1, Math.round(pos * (N + 1)));
466+ $("E0v").textContent = E0 < 10 ? fmt(E0, 2) : fmt(E0, 1);
467+ $("Kbv").textContent = fmt(Kb, 1);
468+ $("Hdv").textContent = fmt(Hd, 2);
469+ $("Nv").textContent = String(N);
470+ { const sum = $("lisSum").getAttribute("aria-pressed") === "true", j = Math.max(1, Math.min(N, Math.round(+$("lis").value * (N + 1))));
471+ $("lisv").textContent = sum ? "all " + N + " beads" : "bead " + j + " of " + N; $("lis").disabled = sum; }
472+ $("posv").textContent = "bead " + center + " of " + N;
473+ $("widv").textContent = fmt(wid, 1) + " beads";
474+ $("dlv").textContent = "1e" + dl + " in strength";
475+ $("tlov").textContent = fmt(tlo, 1) + " s";
476+ $("thiv").textContent = fmt(thi, 2) + " s";
477+ const ftop = f1 * omN / om1;
478+ $("f1v").textContent = Math.round(f1) + " Hz · top " + (ftop >= 1000 ? fmt(ftop / 1000, 1) + "k" : Math.round(ftop));
479+ $("f1note").hidden = !(+$("f1").value > f1max);
480+ $("durv").textContent = fmt(dur, 1) + " s";
481+ const tau = om1 / (2 * Math.PI * f1);
482+ // damping from the ring times: amplitude decay (c0 + c1 lambda) / 2, lambda the Laplacian eigenvalue
483+ const lap1 = 4 * Math.sin(Math.PI / (2 * (N + 1))) ** 2, lapN = 4 * Math.sin(N * Math.PI / (2 * (N + 1))) ** 2;
484+ const sLo = Math.log(1000) / tlo * tau, sHi = Math.log(1000) / thi * tau;
485+ const c1 = Math.max(0, 2 * (sHi - sLo) / (lapN - lap1)), c0 = Math.max(0, 2 * sLo - c1 * lap1);
486+ const fk = Array.from(w.om, (o) => f1 * o / om1);
487+ return {
488+ E0, Kb, Hd, N, T: 1, beta: 1, center, width: wid, delta: Math.pow(10, dl), tau, c0, c1, dur, f1, fk, om1, omN,
489+ dtOut: 1 / SR / tau, nOut: Math.round(dur * SR), decim: Math.round(SR / 400), pickup: Math.max(1, Math.min(N, Math.round(+$("lis").value * (N + 1)))), listenSum: $("lisSum").getAttribute("aria-pressed") === "true",
490+ hmax: 0.7 / Math.sqrt(omN * omN + 12 * E0 / Math.sqrt(N / 32)), snapN: Math.round(0.12 * SR),
491+ };
492+ }
493+ function setControls(p) {
494+ $("E0").value = sliderFromE(p.E0); $("Kb").value = p.Kb; $("Hd").value = p.Hd; $("N").value = p.N; $("pos").value = p.pos;
495+ $("wid").value = p.wid; $("dl").value = p.dl; $("tlo").value = p.tlo; $("thi").value = p.thi; $("f1").value = p.f1; $("dur").value = p.dur;
496+ }
497+
498+ // ---------------- workers (one strike each, in parallel) ----------------
499+ const workerURL = URL.createObjectURL(new Blob([$("worker-src").textContent], { type: "text/javascript" }));
500+ const NW = Math.max(2, Math.min(NS, (navigator.hardwareConcurrency || 4) - 1));
501+ let workers = [], gen = 0, params = null, S = [], queue = [], done = 0, t0run = 0;
502+
503+ function run() {
504+ stopAudio(); stopSlow(true);
505+ params = readParams();
506+ gen++;
507+ for (const w of workers) w.terminate();
508+ workers = [];
509+ S = new Array(NS).fill(null); done = 0; bufs = [];
510+ queue = Array.from({ length: NS }, (_, k) => k);
511+ t0run = performance.now();
512+ for (let i = 0; i < Math.min(NW, NS); i++) {
513+ const w = new Worker(workerURL);
514+ w.onmessage = onResult;
515+ workers.push(w);
516+ next(w);
517+ }
518+ $("status").textContent = "Computing 8 strikes in parallel…";
519+ renderSelector(); updateStats(); drawAll();
520+ }
521+ function next(w) {
522+ if (!queue.length) return;
523+ const k = queue.shift();
524+ w.postMessage({ gen, p: params, k });
525+ }
526+ function onResult(ev) {
527+ const m = ev.data;
528+ if (m.gen !== gen) return;
529+ next(ev.target);
530+ const r = m.r;
531+ S[m.k] = r;
532+ done++;
533+ if (m.k === 0) { r.lam = localLyap(r.logsum); audioGain = 0.8 / (maxAbs(r.sig) || 1); }
534+ if (S[0]) for (let k = 0; k < NS; k++) if (S[k] && !S[k].audio) S[k].audio = scaleAudio(S[k].sig);
535+ $("status").textContent = done < NS ? "Computed " + done + " of 8 strikes…" : "Ready (" + fmt((performance.now() - t0run) / 1000, 1) + " s). Strikes differ by 1e" + $("dl").value + " in strength.";
536+ renderSelector(); updateStats(); drawAll();
537+ }
538+ const maxAbs = (a) => { let m = 0; for (let i = 0; i < a.length; i++) m = Math.max(m, Math.abs(a[i])); return m; };
539+
540+ function localLyap(logsum) {
541+ const dtD = params.decim * params.dtOut;
542+ const w = Math.max(1, Math.round(0.05 / (params.decim / SR))); // 50 ms windows
543+ const out = new Float32Array(logsum.length);
544+ for (let i = 0; i < logsum.length; i++) { const j = Math.max(0, i - w); out[i] = i > j ? (logsum[i] - logsum[j]) / ((i - j) * dtD) / params.tau : NaN; }
545+ return out;
546+ }
547+ function chaosEnd(lam) {
548+ const thr = 40; // per second
549+ let last = -1;
550+ for (let i = 0; i < lam.length; i++) if (lam[i] > thr) last = i;
551+ return last < 0 ? null : (last + 1) * params.decim / SR;
552+ }
553+ function updateStats() {
554+ const A = S[0];
555+ if (!A) { for (const id of ["sBloom", "sBloomT", "sChaos", "sLam", "sSpread"]) $(id).textContent = "–"; return; }
556+ let pk = 0, pi = 0;
557+ for (let i = 0; i < A.hi4.length; i++) if (A.hi4[i] > pk) { pk = A.hi4[i]; pi = i; }
558+ $("sBloom").textContent = fmt(100 * pk, pk < 0.01 ? 2 : 0) + "%";
559+ $("sBloomT").textContent = fmt((pi + 1) * params.decim / SR * 1000, 0) + " ms";
560+ const end = chaosEnd(A.lam);
561+ $("sChaos").textContent = end === null ? "never chaotic" : end > params.dur - 0.15 ? "still chaotic at the end" : fmt(end, 2) + " s";
562+ let peak = 0; for (const v of A.lam) if (v > peak) peak = v;
563+ $("sLam").textContent = peak > 40 ? fmt(peak, 0) + " /s" : "≈ 0";
564+ if (S.some((x) => !x)) { $("sSpread").textContent = "computing…"; return; }
565+ const N = params.N, sd = [];
566+ for (let k = 0; k < N; k++) {
567+ const mean = S.reduce((s, x) => s + x.finalShare[k], 0) / NS;
568+ if (mean < 1e-3) continue;
569+ const vals = S.map((x) => Math.log10(Math.max(x.finalShare[k], 1e-12)));
570+ const m = vals.reduce((s, v) => s + v, 0) / NS;
571+ sd.push(Math.sqrt(vals.reduce((s, v) => s + (v - m) ** 2, 0) / NS));
572+ }
573+ sd.sort((a, b) => a - b);
574+ const med = sd.length ? sd[Math.floor(sd.length / 2)] : 0;
575+ $("sSpread").textContent = med < 0.01 ? "identical" : "differ by " + fmt(10 * med, 1) + " dB";
576+ }
577+
578+ // ---------------- selector ----------------
579+ let sel = 0;
580+ const strikeColor = (k) => (k === 0 ? css("--x") : k === 1 ? css("--y") : css("--ink"));
581+ const strikeName = (k) => (k === 0 ? "A" : k === 1 ? "B" : String(k + 1));
582+ function renderSelector() {
583+ const box = $("strikeSel");
584+ if (!box.childElementCount) {
585+ for (let k = 0; k < NS; k++) {
586+ const b = document.createElement("button");
587+ b.type = "button"; b.className = "sk"; b.id = "sk" + k;
588+ b.innerHTML = '<span class="sw"></span>' + strikeName(k);
589+ b.addEventListener("click", () => selectStrike(k));
590+ box.appendChild(b);
591+ }
592+ const pb = document.createElement("button");
593+ pb.type = "button"; pb.className = "btn"; pb.id = "playSel"; pb.textContent = "Play this strike";
594+ pb.addEventListener("click", () => play(sel, null));
595+ box.appendChild(pb);
596+ }
597+ for (let k = 0; k < NS; k++) {
598+ const b = $("sk" + k);
599+ b.disabled = !S[k];
600+ b.setAttribute("aria-pressed", String(k === sel));
601+ b.querySelector(".sw").style.background = strikeColor(k);
602+ }
603+ $("playSel").disabled = !S[sel];
604+ }
605+ function selectStrike(k) { if (!S[k]) return; sel = k; renderSelector(); drawHi(); drawFinal(); drawSpec(); }
606+ document.addEventListener("keydown", (e) => {
607+ if (e.target && (e.target.tagName === "INPUT" || e.target.tagName === "TEXTAREA")) return;
608+ if (e.key !== "ArrowRight" && e.key !== "ArrowLeft") return;
609+ const dir = e.key === "ArrowRight" ? 1 : -1;
610+ for (let i = 1; i <= NS; i++) { const k = (sel + dir * i + NS) % NS; if (S[k]) { selectStrike(k); e.preventDefault(); return; } }
611+ });
612+
613+ // ---------------- canvas helpers ----------------
614+ function setup(canvas) {
615+ const r = canvas.getBoundingClientRect(), dpr = window.devicePixelRatio || 1;
616+ const w = Math.max(1, Math.round(r.width * dpr)), h = Math.max(1, Math.round(r.height * dpr));
617+ if (canvas.width !== w || canvas.height !== h) { canvas.width = w; canvas.height = h; }
618+ const ctx = canvas.getContext("2d");
619+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0); ctx.clearRect(0, 0, r.width, r.height);
620+ return { ctx, W: r.width, H: r.height };
621+ }
622+ const PAD = { l: 48, r: 12, t: 10, b: 26 };
623+ 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; }
624+ const hz = (f) => (f >= 1000 ? fmt(f / 1000, f >= 10000 || f % 1000 === 0 ? 0 : 1) + "k" : String(Math.round(f)));
625+ function axesTime(ctx, W, H, ymin, ymax, yticks, yfmt, dur) {
626+ const x0 = PAD.l, x1 = W - PAD.r, y0 = PAD.t, y1 = H - PAD.b;
627+ ctx.strokeStyle = css("--grid"); ctx.lineWidth = 1; ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono");
628+ ctx.textAlign = "right"; ctx.textBaseline = "middle";
629+ for (const v of yticks) { const y = y1 - (v - ymin) / (ymax - ymin) * (y1 - y0); ctx.beginPath(); ctx.moveTo(x0, y); ctx.lineTo(x1, y); ctx.stroke(); ctx.fillText(yfmt(v), x0 - 6, y); }
630+ ctx.textAlign = "center"; ctx.textBaseline = "top";
631+ const step = dur <= 2 ? 0.25 : dur <= 5 ? 0.5 : 1;
632+ for (let t = 0; t <= dur + 1e-9; t += step) { const x = x0 + t / dur * (x1 - x0); ctx.beginPath(); ctx.moveTo(x, y0); ctx.lineTo(x, y1); ctx.stroke(); ctx.fillText(fmt(t, step < 1 ? (step < 0.5 ? 2 : 1) : 0) + (t === 0 ? " s" : ""), x, y1 + 6); }
633+ return { x0, x1, y0, y1, X: (t) => x0 + t / dur * (x1 - x0), Y: (v) => y1 - (v - ymin) / (ymax - ymin) * (y1 - y0) };
634+ }
635+ function chaosSpan(ctx, ax) {
636+ if (!S[0] || !S[0].lam) return;
637+ const end = chaosEnd(S[0].lam);
638+ if (end === null) return;
639+ ctx.fillStyle = css("--chaos-wash"); ctx.fillRect(ax.x0, ax.y0, ax.X(end) - ax.x0, ax.y1 - ax.y0);
640+ }
641+ function seriesTime(ctx, ax, arr, color, width, dash) {
642+ if (!arr) return;
643+ const dt = params.decim / SR;
644+ ctx.strokeStyle = color; ctx.lineWidth = width; ctx.setLineDash(dash || []); ctx.beginPath();
645+ let started = false;
646+ for (let i = 0; i < arr.length; i++) { const v = arr[i]; if (!Number.isFinite(v)) { started = false; continue; } const x = ax.X((i + 1) * dt), y = ax.Y(v); if (!started) { ctx.moveTo(x, y); started = true; } else ctx.lineTo(x, y); }
647+ ctx.stroke(); ctx.setLineDash([]);
648+ }
649+ 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]; }
650+ 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]; }
651+ function heatColor(v, bg, ink, hot) { return v < 0.72 ? mix(bg, ink, v / 0.72) : mix(ink, hot, (v - 0.72) / 0.28); }
652+
653+ // ---------------- the wire (physical setup) ----------------
654+ let smOn = false, smIdx = 0, smRaf = 0, smLast = 0;
655+ function drawChain() {
656+ const { ctx, W, H } = setup($("cChain"));
657+ if (!params) return;
658+ const N = params.N, x0 = 26, x1 = W - 26, base = H - 34;
659+ const gap = (x1 - x0) / (N + 1), ink = css("--ink"), muted = css("--muted");
660+ // supports
661+ ctx.fillStyle = muted;
662+ for (const X of [x0, x1]) { ctx.beginPath(); ctx.moveTo(X, base); ctx.lineTo(X - 8, base + 14); ctx.lineTo(X + 8, base + 14); ctx.closePath(); ctx.fill(); }
663+ // dome rest shape (drawn up to 40% of the height for the tallest dome)
664+ const domePx = (H - 70) * 0.55 * Math.min(1, params.Hd / 12);
665+ const zx = (j) => 4 * domePx * (j / (N + 1)) * (1 - j / (N + 1));
666+ ctx.strokeStyle = css("--rule"); ctx.setLineDash([4, 4]); ctx.lineWidth = 1;
667+ ctx.beginPath(); for (let j = 0; j <= N + 1; j++) { const X = x0 + j * gap, Y = base - zx(j); if (j === 0) ctx.moveTo(X, Y); else ctx.lineTo(X, Y); } ctx.stroke(); ctx.setLineDash([]);
668+ // displacement
669+ const A = S[0], y = new Float64Array(N);
670+ let ym = 0;
671+ if (A && A.snap) {
672+ const nS = A.snap.length / N;
673+ for (let i = 0; i < nS * N; i += 3) ym = Math.max(ym, Math.abs(A.snap[i]));
674+ if (smOn) { const i = Math.min(nS - 1, Math.floor(smIdx)); for (let j = 0; j < N; j++) y[j] = A.snap[i * N + j]; }
675+ }
676+ const sc = ym > 0 ? 34 / ym : 1;
677+ const px = (j) => x0 + (j + 1) * gap, py = (j) => base - zx(j + 1) - y[j] * sc;
678+ // links, thickness by stretch stiffness
679+ const zp = (j) => 4 * params.Hd * (((j + 1) / (N + 1)) * (1 - (j + 1) / (N + 1)) - (j / (N + 1)) * (1 - j / (N + 1)));
680+ ctx.strokeStyle = css("--brass"); ctx.lineCap = "round";
681+ for (let j = 0; j <= N; j++) {
682+ const ya = j > 0 ? y[j - 1] : 0, yb = j < N ? y[j] : 0, d = yb - ya, s = zp(j) + d;
683+ const k = 1 + 2 * (s * s + zp(j) * d + 0.5 * d * d);
684+ ctx.lineWidth = Math.min(6, 1 + 1.6 * Math.log(k) + 0.5);
685+ const xa = j > 0 ? px(j - 1) : x0, yA = j > 0 ? py(j - 1) : base, xb = j < N ? px(j) : x1, yB = j < N ? py(j) : base;
686+ ctx.beginPath(); ctx.moveTo(xa, yA); ctx.lineTo(xb, yB); ctx.stroke();
687+ }
688+ const r = Math.max(2.5, Math.min(6, gap * 0.3));
689+ ctx.fillStyle = ink;
690+ for (let j = 0; j < N; j++) { ctx.beginPath(); ctx.arc(px(j), py(j), r, 0, 2 * Math.PI); ctx.fill(); }
691+ const pk = params.pickup - 1;
692+ ctx.strokeStyle = css("--y"); ctx.fillStyle = css("--y"); ctx.font = "11px " + css("--mono"); ctx.textAlign = "center"; ctx.textBaseline = "top";
693+ if (params.listenSum) {
694+ ctx.lineWidth = 1.2;
695+ for (let j = 0; j < N; j++) { ctx.beginPath(); ctx.arc(px(j), py(j), r + 3, 0, 2 * Math.PI); ctx.stroke(); }
696+ ctx.fillText("listening to all beads (sum)", (x0 + x1) / 2, H - 16);
697+ } else {
698+ ctx.lineWidth = 2; ctx.beginPath(); ctx.arc(px(pk), py(pk), r + 5, 0, 2 * Math.PI); ctx.stroke();
699+ ctx.fillText("listener", px(pk), H - 16);
700+ }
701+ if (!smOn) {
702+ ctx.strokeStyle = css("--chaos"); ctx.fillStyle = css("--chaos"); ctx.lineWidth = 2.5;
703+ let topY = H;
704+ for (let j = 0; j < N; j++) {
705+ const s = Math.exp(-0.5 * ((j + 1 - params.center) / params.width) ** 2);
706+ if (s < 0.08) continue;
707+ const L = 30 * s, X = px(j), Y0 = py(j) - r - 4;
708+ ctx.beginPath(); ctx.moveTo(X, Y0); ctx.lineTo(X, Y0 - L); ctx.stroke();
709+ 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();
710+ topY = Math.min(topY, Y0 - L - 8);
711+ }
712+ ctx.fillStyle = ink; ctx.font = "600 12px " + css("--sans"); ctx.textAlign = "left"; ctx.textBaseline = "bottom";
713+ ctx.fillText("strike, E = " + (params.E0 < 10 ? fmt(params.E0, 2) : fmt(params.E0, 1)), px(params.center - 1) + 10, Math.max(14, topY + 4));
714+ } else if (A) {
715+ ctx.fillStyle = muted; ctx.font = "11px " + css("--mono"); ctx.textAlign = "left"; ctx.textBaseline = "top";
716+ ctx.fillText("t = " + fmt(smIdx / SR * 1000, 2) + " ms of sound (sideways motion magnified)", x0 + 10, 8);
717+ }
718+ }
719+ function smTick(now) {
720+ const A = S[0];
721+ if (!smOn || !A) return;
722+ const dtReal = smLast ? Math.min(0.1, (now - smLast) / 1000) : 0;
723+ smLast = now;
724+ smIdx += 0.004 * dtReal * SR; // 4 ms of sound per second
725+ if (smIdx >= A.snap.length / params.N - 1) { stopSlow(); return; }
726+ drawChain();
727+ smRaf = requestAnimationFrame(smTick);
728+ }
729+ function startSlow() { if (!S[0]) return; smOn = true; smIdx = 0; smLast = 0; $("slow").setAttribute("aria-pressed", "true"); $("slow").textContent = "Stop slow motion"; smRaf = requestAnimationFrame(smTick); }
730+ function stopSlow(silent) { smOn = false; cancelAnimationFrame(smRaf); $("slow").setAttribute("aria-pressed", "false"); $("slow").textContent = "Slow motion"; if (!silent) drawChain(); }
731+
732+ // ---------------- plots ----------------
733+ function drawModes() {
734+ const { ctx, W, H } = setup($("cModes"));
735+ const ax = axesTime(ctx, W, H, 0, 1, [], () => "", params.dur);
736+ const A = S[0];
737+ if (!A || !A.shares) return;
738+ const N = params.N, nD = A.hi4.length, off = document.createElement("canvas");
739+ off.width = nD; off.height = N;
740+ const octx = off.getContext("2d"), img = octx.createImageData(nD, N);
741+ const bg = hexToRgb(css("--panel")), ink = hexToRgb(css("--ink")), hot = hexToRgb(css("--chaos"));
742+ for (let i = 0; i < nD; i++) for (let k = 0; k < N; k++) {
743+ const v = Math.min(1, Math.max(0, (Math.log10(Math.max(A.shares[i * N + k], 1e-12)) + 6) / 6));
744+ const c = heatColor(v, bg, ink, hot), o = ((N - 1 - k) * nD + i) * 4;
745+ img.data[o] = c[0]; img.data[o + 1] = c[1]; img.data[o + 2] = c[2]; img.data[o + 3] = 255;
746+ }
747+ octx.putImageData(img, 0, 0);
748+ ctx.imageSmoothingEnabled = false;
749+ ctx.drawImage(off, ax.x0, ax.y0, ax.x1 - ax.x0, ax.y1 - ax.y0);
750+ ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono"); ctx.textAlign = "right"; ctx.textBaseline = "middle";
751+ const marks = [100, 200, 500, 1000, 2000, 4000, 8000, 16000];
752+ let lastY = 1e9;
753+ for (const f of marks) {
754+ let k = -1; for (let i = 0; i < N; i++) if (params.fk[i] <= f) k = i;
755+ if (k < 0 || params.fk[N - 1] < f * 0.8) continue;
756+ const y = ax.y1 - (k + 0.5) / N * (ax.y1 - ax.y0);
757+ if (lastY - y < 13) continue;
758+ ctx.fillText(hz(params.fk[k]), ax.x0 - 6, y); lastY = y;
759+ }
760+ }
761+ function drawHi() {
762+ const { ctx, W, H } = setup($("cHi"));
763+ const ax = axesTime(ctx, W, H, 0, 1, [0, 0.25, 0.5, 0.75, 1], (v) => fmt(v * 100, 0) + "%", params.dur);
764+ chaosSpan(ctx, ax);
765+ for (let k = NS - 1; k >= 0; k--) {
766+ const x = S[k]; if (!x) continue;
767+ const col = k === 0 ? css("--x") : k === 1 ? css("--y") : css("--faint");
768+ seriesTime(ctx, ax, x.hi1, col, 0.8);
769+ seriesTime(ctx, ax, x.hi4, col, k <= 1 ? 1.6 : 1);
770+ }
771+ if (S[sel]) { seriesTime(ctx, ax, S[sel].hi4, strikeColor(sel), 3); seriesTime(ctx, ax, S[sel].hi1, strikeColor(sel), 1.4, [4, 3]); }
772+ }
773+ function drawFinal() {
774+ const { ctx, W, H } = setup($("cFinal"));
775+ if (!params) return;
776+ const N = params.N, x0 = PAD.l, x1 = W - PAD.r, y0 = PAD.t, y1 = H - PAD.b, lo = -8, hi = 0;
777+ const lf0 = Math.log10(params.fk[0]), lf1 = Math.log10(params.fk[N - 1]);
778+ const X = (f) => x0 + (Math.log10(f) - lf0) / Math.max(1e-9, lf1 - lf0) * (x1 - x0), Y = (v) => y1 - (v - lo) / (hi - lo) * (y1 - y0);
779+ ctx.strokeStyle = css("--grid"); ctx.lineWidth = 1; ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono");
780+ ctx.textAlign = "right"; ctx.textBaseline = "middle";
781+ 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)); }
782+ ctx.textAlign = "center"; ctx.textBaseline = "top";
783+ for (const f of [100, 200, 500, 1000, 2000, 5000, 10000, 20000]) if (f >= params.fk[0] * 0.95 && f <= params.fk[N - 1] * 1.05) { ctx.beginPath(); ctx.moveTo(X(f), y0); ctx.lineTo(X(f), y1); ctx.stroke(); ctx.fillText(hz(f), X(f), y1 + 6); }
784+ const line = (s, color, w) => {
785+ if (!s) return;
786+ ctx.strokeStyle = color; ctx.lineWidth = w; ctx.beginPath();
787+ for (let k = 0; k < N; k++) { const v = Math.max(lo, Math.log10(Math.max(s[k], 1e-12))); if (k === 0) ctx.moveTo(X(params.fk[k]), Y(v)); else ctx.lineTo(X(params.fk[k]), Y(v)); }
788+ ctx.stroke();
789+ };
790+ for (let k = 2; k < NS; k++) if (S[k]) line(S[k].finalShare, css("--faint"), 1);
791+ if (S[1]) line(S[1].finalShare, css("--y"), 1.4);
792+ if (S[0]) line(S[0].finalShare, css("--x"), 1.6);
793+ if (S[sel]) line(S[sel].finalShare, strikeColor(sel), 3);
794+ }
795+ function drawLam() {
796+ const { ctx, W, H } = setup($("cLam"));
797+ const A = S[0];
798+ if (!A || !A.lam) { axesTime(ctx, W, H, 0, 1, [0, 1], () => "", params.dur); return; }
799+ let mx = 0; for (const v of A.lam) if (Number.isFinite(v) && v > mx) mx = v;
800+ mx = Math.max(mx * 1.1, 60);
801+ const st = niceStep(mx / 3), ticks = [];
802+ for (let v = 0; v <= mx; v += st) ticks.push(v);
803+ const ax = axesTime(ctx, W, H, -0.15 * mx, mx, ticks, (v) => fmt(v, 0), params.dur);
804+ chaosSpan(ctx, ax);
805+ ctx.strokeStyle = css("--rule"); ctx.beginPath(); ctx.moveTo(ax.x0, ax.Y(0)); ctx.lineTo(ax.x1, ax.Y(0)); ctx.stroke();
806+ const e = A.energy, e0 = e[0] || 1, eDb = new Float32Array(e.length);
807+ 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;
808+ seriesTime(ctx, ax, eDb, css("--muted"), 1, [4, 4]);
809+ seriesTime(ctx, ax, A.lam, css("--x"), 1.6);
810+ }
811+
812+ // ---------------- spectrogram ----------------
813+ function fft(re, im) {
814+ const n = re.length;
815+ 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]]; } }
816+ for (let len = 2; len <= n; len <<= 1) {
817+ const ang = -2 * Math.PI / len, wr = Math.cos(ang), wi = Math.sin(ang);
818+ 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; } }
819+ }
820+ }
821+ const specFmax = () => Math.min(SR / 2, Math.max(2000, 1.1 * params.fk[params.N - 1]));
822+ function spectrogramOf(x) {
823+ if (x.spec) return x.spec;
824+ const sig = x.sig, N = 1024, hop = 256, frames = Math.max(1, Math.floor((sig.length - N) / hop));
825+ const bins = Math.floor(specFmax() / (SR / N));
826+ 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));
827+ const mag = new Float32Array(frames * bins), re = new Float64Array(N), im = new Float64Array(N);
828+ let mx = -Infinity;
829+ for (let f = 0; f < frames; f++) {
830+ for (let i = 0; i < N; i++) { re[i] = sig[f * hop + i] * win[i]; im[i] = 0; }
831+ fft(re, im);
832+ 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; }
833+ }
834+ x.spec = { mag, frames, bins, mx };
835+ return x.spec;
836+ }
837+ function drawSpec() {
838+ const { ctx, W, H } = setup($("cSpec"));
839+ const ax = axesTime(ctx, W, H, 0, 1, [], () => "", params.dur);
840+ const X = S[sel] || S[0];
841+ if (!X || !S[0]) return;
842+ const ref = spectrogramOf(S[0]).mx, { mag, frames, bins } = spectrogramOf(X);
843+ const off = document.createElement("canvas"); off.width = frames; off.height = bins;
844+ const octx = off.getContext("2d"), img = octx.createImageData(frames, bins);
845+ const bg = hexToRgb(css("--panel")), ink = hexToRgb(css("--ink")), hot = hexToRgb(css("--chaos"));
846+ for (let f = 0; f < frames; f++) for (let b = 0; b < bins; b++) {
847+ const v = Math.min(1, Math.max(0, (mag[f * bins + b] - (ref - 80)) / 80));
848+ const c = heatColor(v, bg, ink, hot), o = ((bins - 1 - b) * frames + f) * 4;
849+ img.data[o] = c[0]; img.data[o + 1] = c[1]; img.data[o + 2] = c[2]; img.data[o + 3] = 255;
850+ }
851+ octx.putImageData(img, 0, 0);
852+ const t0 = 512 / SR, t1 = (frames * 256 + 512) / SR;
853+ ctx.drawImage(off, ax.X(t0), ax.y0, ax.X(t1) - ax.X(t0), ax.y1 - ax.y0);
854+ ctx.fillStyle = css("--muted"); ctx.font = "11px " + css("--mono"); ctx.textAlign = "right"; ctx.textBaseline = "middle";
855+ const fm = specFmax(), step = niceStep(fm / 4);
856+ for (let f = 0; f <= fm; f += step) ctx.fillText(hz(f), ax.x0 - 6, ax.y1 - f / fm * (ax.y1 - ax.y0));
857+ const ks = S[sel] ? sel : 0;
858+ ctx.fillStyle = strikeColor(ks); ctx.font = "600 12px " + css("--sans"); ctx.textAlign = "right"; ctx.textBaseline = "top";
859+ ctx.fillText("strike " + strikeName(ks), ax.x1 - 6, ax.y0 + 4);
860+ }
861+
862+ function drawAll() { if (!params) return; drawChain(); drawModes(); drawHi(); drawFinal(); drawLam(); drawSpec(); drawOverlays(null, null); }
863+
864+ // ---------------- audio ----------------
865+ let actx = null, bufs = [], src = null, playStart = 0, playing = null, raf = 0, audioGain = 1;
866+ function scaleAudio(s) { const o = new Float32Array(s.length), fade = Math.round(0.01 * SR); 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))); } return o; }
867+ function ensureCtx() { if (!actx) actx = new (window.AudioContext || window.webkitAudioContext)(); if (actx.state === "suspended") actx.resume(); }
868+ function toBuffer(d) { const b = actx.createBuffer(1, d.length, SR); b.copyToChannel(d, 0); return b; }
869+ function stopAudio() {
870+ if (src) { try { src.onended = null; src.stop(); } catch (e) {} src = null; }
871+ playing = null; cancelAnimationFrame(raf);
872+ for (const id of ["playA", "playB", "playAB", "playSel"]) if ($(id)) $(id).setAttribute("aria-pressed", "false");
873+ drawOverlays(null, null);
874+ }
875+ function play(k, then) {
876+ const x = S[k];
877+ if (!x || !x.audio) return;
878+ ensureCtx(); stopAudio();
879+ if (!bufs[k]) bufs[k] = toBuffer(x.audio);
880+ src = actx.createBufferSource(); src.buffer = bufs[k]; src.connect(actx.destination);
881+ const nx = then === null || then === undefined ? null : then;
882+ src.onended = () => { src = null; if (nx !== null) setTimeout(() => play(nx, null), 350); else stopAudio(); };
883+ playStart = actx.currentTime; playing = k; src.start();
884+ if (k === 0) $("playA").setAttribute("aria-pressed", "true"); else if (k === 1) $("playB").setAttribute("aria-pressed", "true");
885+ if (k === sel && $("playSel")) $("playSel").setAttribute("aria-pressed", "true");
886+ if (nx !== null) $("playAB").setAttribute("aria-pressed", "true");
887+ tick();
888+ }
889+ function tick() { if (playing === null) return; drawOverlays(actx.currentTime - playStart, playing); raf = requestAnimationFrame(tick); }
890+ function drawOverlays(t, which) {
891+ for (const [id, base] of [["cModesO", "cModes"], ["cHiO", "cHi"], ["cLamO", "cLam"], ["cSpecO", "cSpec"]]) {
892+ const { ctx, W } = setup($(id));
893+ if (t === null) continue;
894+ const H = $(base).getBoundingClientRect().height, x = PAD.l + Math.min(t, params.dur) / params.dur * (W - PAD.r - PAD.l);
895+ ctx.strokeStyle = strikeColor(which); ctx.lineWidth = 2; ctx.beginPath(); ctx.moveTo(x, PAD.t); ctx.lineTo(x, H - PAD.b); ctx.stroke();
896+ }
897+ }
898+
899+ // ---------------- wiring ----------------
900+ $("playA").addEventListener("click", () => play(0, null));
901+ $("playB").addEventListener("click", () => play(1, null));
902+ $("playAB").addEventListener("click", () => play(0, 1));
903+ $("stop").addEventListener("click", stopAudio);
904+ for (const [id, other] of [["lisOne", "lisSum"], ["lisSum", "lisOne"]]) $(id).addEventListener("click", () => {
905+ if ($(id).getAttribute("aria-pressed") === "true") return;
906+ $(id).setAttribute("aria-pressed", "true"); $(other).setAttribute("aria-pressed", "false"); readParams(); run();
907+ });
908+ $("slow").addEventListener("click", () => (smOn ? stopSlow() : startSlow()));
909+ let timer = 0;
910+ for (const id of ["lis", "E0", "Kb", "Hd", "N", "pos", "wid", "dl", "tlo", "thi", "f1", "dur"]) $(id).addEventListener("input", () => { readParams(); clearTimeout(timer); timer = setTimeout(run, 350); });
911+ document.querySelectorAll("[data-preset]").forEach((b) => b.addEventListener("click", () => { setControls(PRESETS[b.dataset.preset]); run(); }));
912+ let rt = 0;
913+ window.addEventListener("resize", () => { clearTimeout(rt); rt = setTimeout(drawAll, 120); });
914+ try {
915+ window.matchMedia("(prefers-color-scheme: dark)").addEventListener("change", drawAll);
916+ new MutationObserver(drawAll).observe(document.documentElement, { attributes: true, attributeFilter: ["data-theme"] });
917+ } catch (e) {}
918+
919+ setControls(PRESETS.gong);
920+ run();
921+})();
922+</script>
wire.pyadded+168−0View file
@@ -0,0 +1,168 @@
1+"""A struck beaded wire: the FPUT chain generalized toward a cymbal or gong.
2+
3+N unit masses on a wire between two pinned ends move sideways, with displacements y_j (y_0 = y_{N+1} = 0).
4+The wire may be built as a shallow arch (a "dome") of height Hd: rest shape z_j = 4 Hd x (1 - x), x = j / (N + 1).
5+Energy:
6+ links: sum_j [ T D_j^2 / 2 + beta (z'_j D_j + D_j^2 / 2)^2 ], D_j = y_{j+1} - y_j, z'_j = z_{j+1} - z_j
7+ bending: sum_j Kb (y_{j+1} - 2 y_j + y_{j-1})^2 / 2 (pinned ends: y_{-1} = -y_1, y_{N+2} = -y_N)
8+ on-site: mu^2 sum_j y_j^2 / 2
9+With Hd = 0 and Kb = 0 this is the FPUT-beta chain (sideways motion of a pre-tensioned beaded string; the link
10+stretch is about D^2 / 2, so the stretching energy is quartic in D). The bending term makes the spectrum
11+dispersive, omega ~ k^2 at large Kb, like a plate. The dome makes the link stretch linear in D as well
12+(z' D), which adds stiffness (as curvature does in a shell) and a quadratic nonlinearity (z' D^3), the
13+three-wave coupling of gongs and cymbals.
14+
15+The linear part is K = T L + 2 beta diag-ish(z'^2) + Kb B + mu^2 I; with a dome its normal modes are not sines, so
16+they are computed numerically (eigh). Damping: friction c0 and neighbour dashpots c1 (Laplacian of the
17+velocity), applied per step with Crank-Nicolson (tridiagonal). Time stepping: Stormer-Verlet (leapfrog).
18+"""
19+
20+from __future__ import annotations
21+
22+from dataclasses import dataclass
23+
24+import numpy as np
25+from scipy.linalg import eigh, solve_banded
26+
27+
28+@dataclass(frozen=True)
29+class Wire:
30+ N: int = 64
31+ T: float = 1.0
32+ beta: float = 1.0
33+ Kb: float = 0.0
34+ Hd: float = 0.0
35+ mu: float = 0.0
36+ c0: float = 0.0
37+ c1: float = 0.0
38+
39+ @property
40+ def zp(self):
41+ x = np.arange(self.N + 2) / (self.N + 1)
42+ z = 4 * self.Hd * x * (1 - x)
43+ return np.diff(z) # (N + 1,) slopes of the links
44+
45+ def force(self, y):
46+ """y: (M, N) -> force (M, N)."""
47+ M, N = y.shape
48+ yp = np.pad(y, ((0, 0), (1, 1)))
49+ D = np.diff(yp, axis=1) # (M, N + 1)
50+ zp = self.zp
51+ e = zp * D + 0.5 * D * D # link stretch
52+ fl = self.T * D + 2 * self.beta * e * (zp + D) # dV/dD
53+ f = fl[:, 1:] - fl[:, :-1]
54+ if self.Kb:
55+ # bending: curvature c_j = y_{j+1} - 2 y_j + y_{j-1}, j = 0..N+1 with pinned ghosts
56+ yg = np.concatenate([-y[:, :1], np.zeros((M, 1)), y, np.zeros((M, 1)), -y[:, -1:]], axis=1) # j = -1..N+2
57+ c = yg[:, 2:] - 2 * yg[:, 1:-1] + yg[:, :-2] # j = 0..N+1
58+ # dU/dy_j = Kb (c_{j-1} - 2 c_j + c_{j+1}), j = 1..N
59+ f -= self.Kb * (c[:, :-2] - 2 * c[:, 1:-1] + c[:, 2:])
60+ f -= self.mu**2 * y
61+ return f
62+
63+ def stiffness(self):
64+ """Linear stiffness matrix K (N x N) at y = 0."""
65+ N = self.N
66+ kl = self.T + 2 * self.beta * self.zp**2 # link stiffness at rest
67+ K = np.zeros((N, N))
68+ for j in range(N + 1): # link j joins y_j and y_{j+1} (1-based; 0 and N+1 fixed)
69+ a, b = j - 1, j
70+ if a >= 0:
71+ K[a, a] += kl[j]
72+ if b < N:
73+ K[b, b] += kl[j]
74+ if a >= 0 and b < N:
75+ K[a, b] -= kl[j]
76+ K[b, a] -= kl[j]
77+ if self.Kb:
78+ L = -2 * np.eye(N) + np.eye(N, k=1) + np.eye(N, k=-1) # Dirichlet Laplacian; pinned bending = Kb L^2
79+ K += self.Kb * L @ L
80+ K += self.mu**2 * np.eye(N)
81+ return K
82+
83+ def modes(self):
84+ w2, V = eigh(self.stiffness())
85+ return np.sqrt(np.maximum(w2, 0)), V
86+
87+ def energy(self, y, v):
88+ M, N = y.shape
89+ yp = np.pad(y, ((0, 0), (1, 1)))
90+ D = np.diff(yp, axis=1)
91+ e = self.zp * D + 0.5 * D * D
92+ U = (0.5 * self.T * D * D + self.beta * e * e).sum(1)
93+ if self.Kb:
94+ yg = np.concatenate([-y[:, :1], np.zeros((M, 1)), y, np.zeros((M, 1)), -y[:, -1:]], axis=1)
95+ c = yg[:, 2:] - 2 * yg[:, 1:-1] + yg[:, :-2]
96+ U += 0.5 * self.Kb * (c[:, 1:-1] ** 2).sum(1)
97+ U += 0.5 * self.mu**2 * (y * y).sum(1)
98+ return 0.5 * (v * v).sum(1) + U
99+
100+
101+def strike(w: Wire, E, center, width, M=1, scale=None):
102+ j = np.arange(1, w.N + 1)
103+ prof = np.exp(-0.5 * ((j - center) / width) ** 2)
104+ prof *= np.sqrt(2 * E / (prof**2).sum())
105+ v = np.tile(prof, (M, 1))
106+ if scale is not None:
107+ v *= np.asarray(scale)[:, None]
108+ return np.zeros_like(v), v
109+
110+
111+def run(w: Wire, y0, v0, t_end, dt, record_every=100, lyap=False):
112+ """Stormer-Verlet (kick-drift-kick), damping by Crank-Nicolson on v after each step."""
113+ y, v = np.array(y0, float), np.array(v0, float)
114+ M, N = y.shape
115+ om, V = w.modes()
116+ damp = bool(w.c0 or w.c1)
117+ if damp:
118+ # (I + dt/2 Dm) v+ = (I - dt/2 Dm) v, Dm = c0 I - c1 L
119+ a = 0.5 * dt
120+ ab = np.zeros((3, N))
121+ ab[0, 1:] = -a * w.c1
122+ ab[1, :] = 1 + a * (w.c0 + 2 * w.c1)
123+ ab[2, :-1] = -a * w.c1
124+
125+ def apply_damp(x):
126+ rhs = (1 - a * (w.c0 + 2 * w.c1)) * x
127+ rhs[:, 1:] += a * w.c1 * x[:, :-1]
128+ rhs[:, :-1] += a * w.c1 * x[:, 1:]
129+ return solve_banded((1, 1), ab, rhs.T).T
130+ if lyap:
131+ rng = np.random.default_rng(0)
132+ dy, dv = rng.standard_normal(y.shape), rng.standard_normal(v.shape)
133+ nrm = np.sqrt((dy**2).sum(1) + (dv**2).sum(1))[:, None]
134+ dy /= nrm
135+ dv /= nrm
136+ eps = 1e-7
137+ logsum = np.zeros(M)
138+ n = int(round(t_end / dt))
139+ rec = dict(t=[], Ek=[], E=[], logsum=[])
140+ f = w.force(y)
141+ for i in range(n):
142+ v = v + 0.5 * dt * f
143+ y = y + dt * v
144+ f = w.force(y)
145+ v = v + 0.5 * dt * f
146+ if lyap:
147+ # tangent by finite differences of the force (keeps the model in one place)
148+ g = (w.force(y + eps * dy) - f) / eps
149+ dv = dv + dt * g
150+ dy = dy + dt * dv
151+ if damp:
152+ v = apply_damp(v)
153+ if lyap:
154+ dv = apply_damp(dv)
155+ if lyap and (i + 1) % 10 == 0:
156+ nrm = np.sqrt((dy**2).sum(1) + (dv**2).sum(1))
157+ logsum += np.log(nrm)
158+ dy /= nrm[:, None]
159+ dv /= nrm[:, None]
160+ if (i + 1) % record_every == 0:
161+ q, p = y @ V, v @ V
162+ rec["t"].append((i + 1) * dt)
163+ rec["Ek"].append(0.5 * (p * p + (om * q) ** 2))
164+ rec["E"].append(w.energy(y, v))
165+ rec["logsum"].append(logsum.copy())
166+ out = {k: np.array(x) for k, x in rec.items()}
167+ out.update(y=y, v=v, omega=om)
168+ return out