1// The crash in miniature: a Fermi-Pasta-Ulam chain, masses on springs whose
2// force carries a quadratic (geometric) correction, started in its lowest
3// mode. At small amplitude the energy stays where it was put; at large
4// amplitude the nonlinearity couples the modes and energy migrates up the
5// ladder. The shell's bending-stretching coupling does to its nodal-diameter
6// modes what this chain does to its sine modes.
8import { useMemo, useState } from 'react';
9import { useCanvas, drawFrame, PLOT_COLORS, colormap } from '../../components/plots/canvasUtils';
10import { Figure, FigControls, FigSlider, usePlayer } from '../figparts';
12const HEIGHT = 260;
13const N = 24; // interior masses
14const ALPHA = 0.25; // quadratic coupling
15const TSTEPS = 30000;
16const DT = 0.06;
17const FRAMES = 240; // energy snapshots across the run
19function simulate(amp: number) {
20 const u = new Float64Array(N + 2);
21 const v = new Float64Array(N + 2);
22 for (let i = 1; i <= N; i++) u[i] = amp * Math.sin((Math.PI * i) / (N + 1));
24 const wk = new Float64Array(N);
25 for (let k = 1; k <= N; k++) wk[k - 1] = 2 * Math.sin((Math.PI * k) / (2 * (N + 1)));
27 const energies = new Float32Array(FRAMES * N);
28 const a = new Float64Array(N + 2);
29 const accel = (uu: Float64Array) => {
30 for (let i = 1; i <= N; i++) {
31 const dr = uu[i + 1] - uu[i];
32 const dl = uu[i] - uu[i - 1];
33 a[i] = dr - dl + ALPHA * (dr * dr - dl * dl);
34 }
35 };
37 const snap = (frame: number) => {
38 for (let k = 1; k <= N; k++) {
39 let q = 0, p = 0;
40 for (let i = 1; i <= N; i++) {
41 const s = Math.sin((Math.PI * k * i) / (N + 1));
42 q += u[i] * s;
43 p += v[i] * s;
44 }
45 const c = Math.sqrt(2 / (N + 1));
46 q *= c;
47 p *= c;
48 energies[frame * N + (k - 1)] = 0.5 * (p * p + wk[k - 1] * wk[k - 1] * q * q);
49 }
50 };
52 const perFrame = Math.floor(TSTEPS / FRAMES);
53 // leapfrog: cheap and symplectic, so the ladder is not an artifact of decay
54 accel(u);
55 for (let f = 0; f < FRAMES; f++) {
56 snap(f);
57 for (let s = 0; s < perFrame; s++) {
58 for (let i = 1; i <= N; i++) v[i] += 0.5 * DT * a[i];
59 for (let i = 1; i <= N; i++) u[i] += DT * v[i];
60 accel(u);
61 for (let i = 1; i <= N; i++) v[i] += 0.5 * DT * a[i];
62 }
63 }
64 return energies;
65}
67// Sonification: the same chain, retuned so mode 1 sits at 220 Hz (its top
68// mode then lands near 3.5 kHz), lightly damped, and integrated at audio
69// rate; the output is the acceleration of one off-center mass. At small
70// amplitude it is a plain 220 Hz tone; at large amplitude the cascade is
71// audible as a brightening after the onset.
72const AUDIO_SR = 24000;
73const AUDIO_F1 = 220;
75export function renderChainAudio(amp: number): Float32Array<ArrayBuffer> {
76 const dur = 2.5;
77 const n = Math.round(AUDIO_SR * dur);
78 const u = new Float64Array(N + 2);
79 const v = new Float64Array(N + 2);
80 const a = new Float64Array(N + 2);
81 for (let i = 1; i <= N; i++) u[i] = amp * Math.sin((Math.PI * i) / (N + 1));
82 const w1 = 2 * Math.sin(Math.PI / (2 * (N + 1)));
83 const tauPerSec = (2 * Math.PI * AUDIO_F1) / w1; // natural time per real second
84 const sub = 2;
85 const dtau = tauPerSec / AUDIO_SR / sub;
86 const cdamp = (2 * 2.76) / tauPerSec; // amplitude T60 near 2.5 s
88 const accel = () => {
89 for (let i = 1; i <= N; i++) {
90 const dr = u[i + 1] - u[i];
91 const dl = u[i] - u[i - 1];
92 a[i] = dr - dl + ALPHA * (dr * dr - dl * dl) - cdamp * v[i];
93 }
94 };
95 const iOut = 5;
96 const y = new Float32Array(n);
97 accel();
98 for (let k = 0; k < n; k++) {
99 for (let s = 0; s < sub; s++) {
100 for (let i = 1; i <= N; i++) v[i] += 0.5 * dtau * a[i];
101 for (let i = 1; i <= N; i++) u[i] += dtau * v[i];
102 accel();
103 for (let i = 1; i <= N; i++) v[i] += 0.5 * dtau * a[i];
104 }
105 y[k] = a[iOut];
106 }
107 let pk = 0;
108 for (let k = 0; k < n; k++) pk = Math.max(pk, Math.abs(y[k]));
109 const fadeIn = Math.round(0.004 * AUDIO_SR);
110 const fadeOut = Math.round(0.12 * AUDIO_SR);
111 for (let k = 0; k < n; k++) {
112 let g = 0.45 / (pk || 1);
113 if (k < fadeIn) g *= k / fadeIn;
114 if (k > n - fadeOut) g *= (n - k) / fadeOut;
115 y[k] *= g;
116 }
117 return y;
118}
120export default function CascadeFig() {
121 const [amp, setAmp] = useState(4);
122 const energies = useMemo(() => simulate(amp), [amp]);
123 const play = usePlayer();
125 const ref = useCanvas(
126 (ctx, w, h) => {
127 drawFrame(ctx, w, h);
128 const padL = 70;
129 const padR = 14;
130 const padT = 10;
131 const padB = 28;
132 const pw = w - padL - padR;
133 const ph = h - padT - padB;
135 let e0 = 0;
136 for (let k = 0; k < N; k++) e0 += energies[k];
137 const floorDb = -55;
139 // paint at data resolution, then stretch (putImageData ignores the
140 // dpr transform, so go through an offscreen canvas as the plots do)
141 const off = document.createElement('canvas');
142 off.width = FRAMES;
143 off.height = N;
144 const octx = off.getContext('2d')!;
145 const img = octx.createImageData(FRAMES, N);
146 for (let mode = 0; mode < N; mode++) {
147 const row = N - 1 - mode; // mode 1 at the bottom
148 for (let frame = 0; frame < FRAMES; frame++) {
149 const e = energies[frame * N + mode] / (e0 || 1);
150 const db = 10 * Math.log10(Math.max(1e-9, e));
151 const t = Math.min(1, Math.max(0, (db - floorDb) / -floorDb));
152 const [r, g, b] = colormap(t);
153 const o = 4 * (row * FRAMES + frame);
154 img.data[o] = r;
155 img.data[o + 1] = g;
156 img.data[o + 2] = b;
157 img.data[o + 3] = 255;
158 }
159 }
160 octx.putImageData(img, 0, 0);
161 ctx.imageSmoothingEnabled = false;
162 ctx.drawImage(off, padL, padT, pw, ph);
164 ctx.fillStyle = PLOT_COLORS.axis;
165 ctx.font = '10px system-ui';
166 ctx.textAlign = 'right';
167 for (const k of [1, 6, 12, 18, 24]) {
168 const yy = padT + ph - ((k - 0.5) / N) * ph;
169 ctx.fillText(`${k}`, padL - 6, yy + 3);
170 }
171 ctx.save();
172 ctx.translate(14, padT + ph / 2);
173 ctx.rotate(-Math.PI / 2);
174 ctx.textAlign = 'center';
175 ctx.fillText('mode number (energy, dB re start)', 0, 0);
176 ctx.restore();
177 ctx.textAlign = 'center';
178 ctx.fillText('time →', padL + pw / 2, h - 8);
179 ctx.textAlign = 'left';
180 ctx.fillText('all energy starts in mode 1 (bottom row)', padL + 6, h - 8);
181 },
182 [energies],
183 HEIGHT,
184 );
186 return (
187 <Figure
188 caption={
189 <>
190 Modal energy (color, log scale) against time for a chain of {N} masses whose spring
191 forces carry a quadratic geometric term, started with all of its energy in mode 1.
192 Below an amplitude near 1 the system is effectively linear and the energy stays in
193 the bottom row for the whole run. Push the amplitude up and the modes cease to be
194 independent: energy climbs the ladder, and at the top of the range the spectrum
195 fills in. This, plus damping, is the crash: the strike puts energy in at the bottom,
196 and the nonlinearity spends the next several hundred milliseconds carrying it
197 upward, which is why a crash brightens after the attack rather than at it. The play
198 button integrates the same chain retuned so its lowest mode sits at 220 Hz and
199 lightly damped: quiet strikes give a plain tone, and somewhere past an amplitude of
200 2 the brightening-after-onset becomes plainly audible.
201 </>
202 }
203 >
204 <canvas ref={ref} />
205 <FigControls>
206 <FigSlider
207 label="strike amplitude"
208 value={amp}
209 min={0.2}
210 max={8}
211 step={0.1}
212 format={(v) => v.toFixed(1)}
213 onChange={setAmp}
214 />
215 <button onClick={() => play(renderChainAudio(amp), AUDIO_SR)}>
216 ▶ Hear this chain (mode 1 tuned to 220 Hz)
217 </button>
218 </FigControls>
219 </Figure>
220 );
221}