/ listening-for-chaos / chaos-apps
Sign in
listening-for-chaos / chaos-apps
Beam: mallet strike (Hertz contact, mass M and stiffness K in Parameters; shared contact routine in src/common/mallet.ts); hammer shape factored out; kick stays the default
Jeremy Magland <jmagland@flatironinstitute.org> committed commit 5b795e3f4571 parent 8ce184c Browse files
5 changed files+139−20
src/beam/README.mdmodified+1−0View file
@@ -12,6 +12,7 @@ Live site: https://chaos-apps--listening-for-chaos.vault1-sites.magland.org/beam
1212
1313 - SI units. Defaults m = 1 kg, k = 0, κ = 20 N/m, k3 = 4000 N/m^3, a = 0.5 m (drawing only), τ = 540 s, local stretching, clamped ends.
1414 - Resolution: the Parameters values are for the beam with 8 masses (REF_N); Masses (2 to 64) refines the same beam (`refine` in `beam.ts`): same length, rise, drag, strike energies; with spacing ratio r = span(8) / span(N) (span N + 1 spacings between two supports, N for a cantilever) m x r, k / r, κ / r^3, k3 / r^3, a x r, and the loss dashpots / r (τ top is converted so a wave of given length decays alike at any N). N = 8 returns the parameters unchanged (checked, all ends). The derived values show the scaled constants when N is not 8. Before this (to 2026-10-09) more masses made the beam longer at fixed spacing: mode 1 fell as ~1/N^2 (0.211 Hz at N = 8, 0.015 Hz at 32), so at fixed τ the note shrank (7.9 s at N = 8, 0.6 s at 32, 100 Hz pitch). Now (clamped, local, τ 540, τ top 50): f1 0.207 / 0.211 / 0.206 / 0.202 / 0.200 Hz at N = 4 / 8 / 16 / 32 / 64; f2/f1 2.24 / 2.59 / 2.71 / 2.75 / 2.754 and f3/f1 3.37 / 4.68 / 5.20 / 5.35 / 5.39 (beam theory 2.757, 5.404); pinned 3.998, 8.986 at 64 (4, 9); cantilever 6.265, 17.53 (6.27, 17.55); mode-1 decay ~375 s, mode 3 ~122 s; nonlinear scale 176 / 167 / 161 / 158 mJ at 8 / 16 / 32 / 64; default note 7.9 / 7.7 / 7.5 s. Cost: the step falls ~1/N^2 (0.58 ms at 8, 0.087 ms at 32 for 2.51 J), so the default note (with the Lyapunov shadow) took ~1 s to compute at N = 8, 8.2 s at 16 and 55 s at 32 in the pane at step 0.03; at 0.2, 0.3 s, 1.5 s and 7.9 s. Strike and listening are at positions (m), so they mean the same at any N. Default blow at 1.0 m, 2.51 J, clamped (mode-energy shares just after the strike, modes 1 / 2 / 3, and above mode 8; λ over 0 to 300 s with τ 540 s): N = 8 5.1 / 18.8 / 26.2 %, 0 %, λ 0.13 to 0.15 /s; N = 16 8.4 / 26.5 / 30.8 %, 2.7 %, 0.12 to 0.13; N = 32 8.4 / 26.2 / 30.6 %, 2.9 %, 0.068; N = 64 8.7 / 26.8 / 30.3 %, 2.8 %, 0.071. N = 8 differs because the bump covers a single mass there. The old blow on mass 2 (to 2026-10-09) moved up the spectrum with N: above mode 8 38 % at N = 16, 89 % at 32, 99 % at 64; λ 0.09, 0.024, 0.011.
15+- Mallet (strike "mallet at"; `../common/mallet.ts`, beam only): mass M (default 1 kg) arriving with energy E presses into the beam through a Hertz contact, force K δ^1.5 (K default 1000 N/m^1.5), δ the mallet face's depth past the weighted mean displacement under the hammer bump (weights the bump normalized to sum 1; force distributed by the same weights). The contact phase is integrated at a fine step (bounded by the rigid-wall contact time / 400 and the struck masses' contact frequency) until lift-off; drag and loss are left out during it; the note and animation start at lift-off; the readout gives contact time, energy share delivered, peak force. Scan (2.51 J at 1.0 m, M = 1 kg): K 1e2 / 1e3 / 1e4 / 1e5: contact 470 / 150 / 55 / 22 ms (N = 8 and 32 alike), energy in 78-83 / 92-100 / 95-100 / 96-100 %; at N = 32 brightness (energy-weighted mode frequency / f1) just after lift-off 5.4 / 6.7 / 8.5 / 9.3 against 9.7 for the kick (the kick is the hard-mallet limit). Light mallets bounce off (M = 0.25 kg: 40 to 64 % in). In the sound at 2.51 J the difference is small: the chaotic phase redistributes the energy over the modes within a fraction of a second (spectral centroid at 1 s of sound, N = 32: kick 4.3 kHz, mallets 4.0 to 4.8 kHz; with τ top 200 s: kick 363 Hz, mallets 300 to 400 Hz).
1516 - Masses on rails: only sideways motion u_j, so every axial force is taken by the rails and tension and stretching are defined at a free end too.
1617 - Potential: springs between neighbours (and to each support), k d^2 / 2 with d = u_j - u_{j-1}; bending κ c^2 / 2 at each mass with both neighbours, c = u_{j-1} - 2 u_j + u_{j+1} (supports at u = 0). Pinned support: no extra term (ghost u_{-2} = -u_0 gives c = 0 there). Clamped: a joint at the support with ghost u_{-2} = u_0, c = 2 u_0, energy 2 κ u_0^2. Free end: no spring or joint beyond the last mass.
1718 - Stretching: "local", k3 d^4 / 4 per spring (FPUT); or "Kirchhoff", k3 S^2 / (4 nS) with S the sum of d^2 over the nS springs, which acts as one tension k3 S / nS shared by all springs. They agree when all springs are stretched alike.
src/beam/beam.tsmodified+18−14View file
@@ -447,23 +447,27 @@ export function strikeMode(P: Params, modes: Modes, n: number, E: number): State
447447 export const massPosition = (P: Params, j: number) => (j + 1) * P.a;
448448 export const beamLength = (P: Params) => (freeRight(P) ? P.N : P.N + 1) * P.a;
449449
450-// A hammer blow at x (m) of half-width w (m, at least one spacing): velocity
451-// in a triangular bump, 1 - |x_j - x| / w on the masses it covers, scaled
452-// to give energy E. The same blow at any N (refine) when w is fixed: with
453-// w the spacing at REF_N and x on a mass there, only that mass moves at REF_N.
454-export function strikeAt(P: Params, x: number, w: number, E: number): State {
455- const s = newState(P.N), we = Math.max(w, P.a), phi = new Float64Array(P.N);
456- let S = 0;
450+// The hammer's shape: a triangular bump 1 - |x_j - x| / w on the masses it
451+// covers, x (m) its centre and w (m) its half-width (at least one spacing);
452+// at a support beyond the masses' reach, the nearest mass alone. With w the
453+// spacing at REF_N and x on a mass there, only that mass at REF_N.
454+export function hammerShape(P: Params, x: number, w: number): Float64Array {
455+ const we = Math.max(w, P.a), phi = new Float64Array(P.N);
456+ let any = false;
457457 for (let j = 0; j < P.N; j++) {
458458 phi[j] = Math.max(0, 1 - Math.abs(massPosition(P, j) - x) / we);
459- S += phi[j] * phi[j];
460- }
461- if (S === 0) {
462- // x beyond the masses' reach (at a support): the nearest mass.
463- const j = Math.min(P.N - 1, Math.max(0, Math.round(x / P.a) - 1));
464- phi[j] = 1;
465- S = 1;
459+ if (phi[j] > 0) any = true;
466460 }
461+ if (!any) phi[Math.min(P.N - 1, Math.max(0, Math.round(x / P.a) - 1))] = 1;
462+ return phi;
463+}
464+
465+// A hammer blow at x: velocity in the hammer's shape, scaled to give
466+// energy E. The same blow at any N (refine) when w is fixed.
467+export function strikeAt(P: Params, x: number, w: number, E: number): State {
468+ const s = newState(P.N), phi = hammerShape(P, x, w);
469+ let S = 0;
470+ for (let j = 0; j < P.N; j++) S += phi[j] * phi[j];
467471 const scale = Math.sqrt((2 * P.m * E) / S);
468472 for (let j = 0; j < P.N; j++) s.p[j] = phi[j] * scale;
469473 return s;
src/beam/index.htmlmodified+6−2View file
@@ -16,7 +16,8 @@
1616 <button id="strike" class="primary">Strike</button>
1717 <select id="strike-kind">
1818 <option value="mode">mode</option>
19- <option value="at" selected>at</option>
19+ <option value="at" selected>kick at</option>
20+ <option value="mallet">mallet at</option>
2021 </select>
2122 <input id="strike-index" type="number" min="1" value="2" step="1" hidden />
2223 <span id="strike-at"><input id="strike-x" type="number" min="0" value="1" step="0.05" /> <span class="unit">m</span></span>
@@ -28,6 +29,7 @@
2829 <input id="energy" type="range" min="-5" max="2" step="0.01" value="0.4" />
2930 <span id="energy-value" class="num"></span>
3031 </div>
32+ <div class="row"><span id="strike-info" class="muted"></span></div>
3133 <div class="row">
3234 <span class="muted">Masses</span>
3335 <input id="n" type="number" min="2" max="64" step="1" value="8" />
@@ -71,6 +73,8 @@
7173 <label>a <input id="p-a" type="number" step="any" /> <span class="unit">m (spacing)</span></label>
7274 <label>τ <input id="p-tau" type="number" step="any" /> <span class="unit">s (decay, 0 = off)</span></label>
7375 <label>τ top <input id="p-tautop" type="number" step="any" /> <span class="unit">s (decay of the highest modes from the loss, 0 = off)</span></label>
76+ <label>mallet M <input id="p-mallet-m" type="number" step="any" /> <span class="unit">kg</span></label>
77+ <label>mallet K <input id="p-mallet-k" type="number" step="any" /> <span class="unit">N/m^1.5 (contact stiffness)</span></label>
7478 <label class="wide">stretching <select id="p-stretch">
7579 <option value="local">local (each spring)</option>
7680 <option value="kirchhoff">Kirchhoff (one shared tension)</option>
@@ -136,7 +140,7 @@
136140
137141 <template id="help-beam">
138142 <h2>The beam</h2>
139- <p>Equal masses on vertical rails, moving sideways only, shown in real time (displacements scaled to fit). Neighbours are joined by springs (tension k), and each mass resists bending (κ, on how far it sits off the line through its neighbours). The springs also stiffen as they stretch (k3): this is the nonlinearity. Strike gives the energy on the slider as a hammer blow at a point (m from the left support; a short triangular bump of velocity, one spacing of the 8-mass beam wide on each side, so the blow is the same at any number of masses), or in the shape of one normal mode.</p>
143+ <p>Equal masses on vertical rails, moving sideways only, shown in real time (displacements scaled to fit). Neighbours are joined by springs (tension k), and each mass resists bending (κ, on how far it sits off the line through its neighbours). The springs also stiffen as they stretch (k3): this is the nonlinearity. Strike gives the energy on the slider as a hammer blow at a point (m from the left support; a short triangular bump of velocity, one spacing of the 8-mass beam wide on each side, so the blow is the same at any number of masses), with a mallet at a point, or in the shape of one normal mode. A mallet (mass M, contact stiffness K in Parameters) arrives with the strike energy and presses into the beam over the same bump through a Hertz contact (force K δ^1.5), then lifts off; the readout gives the contact time, the share of the energy that went into the beam, and the peak force. The motion is shown from lift-off.</p>
140144 <p>The ends are clamped (fixed, and held level), pinned (fixed, free to turn), or a cantilever: clamped at the left and free at the right. Bending makes the overtones inharmonic, as in a bar or a cymbal. Weak strikes ring as independent normal modes; strong ones trade energy between them. Ghost adds a second beam started a tiny distance away, drawn in orange: in chaotic motion the two separate exponentially. The space bar pauses.</p>
141145 </template>
142146 <template id="help-params">
src/beam/main.tsmodified+33−4View file
@@ -1,7 +1,7 @@
11 import "../common/style.css";
22 import "./style.css";
33 import {
4- Modes, REF_N, Verlet, beamLength, copyState, decayTimes, energy, isDamped, newState, nonlinearScale, refine, restShape, stepFor, strikeAt, strikeMode,
4+ Modes, REF_N, Verlet, beamLength, copyState, decayTimes, energy, forces, hammerShape, isDamped, newState, nonlinearScale, refine, restShape, stepFor, strikeAt, strikeMode,
55 type Ends, type Params, type State, type Stretch,
66 } from "./beam";
77 import { fit, indexColor } from "../common/canvas";
@@ -14,6 +14,7 @@ import { drawModeBars } from "../common/modebars";
1414 import { LiveModeHistory, ModeHistory } from "../common/modehistory";
1515 import { setupMute } from "../common/mute";
1616 import { PASS } from "../common/antialias";
17+import { malletStrike, type MalletParams } from "../common/mallet";
1718 import { DIAG_EVERY, END_RATIO, LYAP_WINDOW, UNDAMPED_LENGTH } from "../common/note-render";
1819 import { ScatterPlot } from "../common/scatter";
1920 import { Player } from "../common/sound";
@@ -35,6 +36,7 @@ const kindSelect = document.getElementById("strike-kind") as HTMLSelectElement;
3536 const indexInput = document.getElementById("strike-index") as HTMLInputElement;
3637 const strikeAtEl = document.getElementById("strike-at")!;
3738 const strikeXInput = document.getElementById("strike-x") as HTMLInputElement;
39+const strikeInfo = document.getElementById("strike-info")!;
3840 const speedSelect = document.getElementById("speed") as HTMLSelectElement;
3941 const pauseButton = document.getElementById("pause") as HTMLButtonElement;
4042
@@ -46,6 +48,8 @@ const pauseButton = document.getElementById("pause") as HTMLButtonElement;
4648 // refined to the chosen number of masses (see refine).
4749 let base: Params = { N: REF_N, m: 1, k: 0, kappa: 20, k3: 4000, stretch: "local", ends: "clamped", rise: 0, a: 0.5, tau: 540, tauTop: 0 };
4850 let P: Params = base;
51+// The mallet (for "mallet at" strikes): a physical object, not refined with N.
52+let mallet: MalletParams = { M: 1, K: 1000 };
4953 const SAMPLE_DT = 0.01; // s between plotted samples
5054
5155 // The beam. It sits at rest until the first Strike; changing parameters puts
@@ -100,10 +104,27 @@ document.getElementById("strike")!.addEventListener("click", () => {
100104 indexInput.value = String(strikeIndex());
101105 strikeXInput.value = String(position(strikeXInput));
102106 const mode = kindSelect.value === "mode";
103- strikeLabel = mode ? `mode${strikeIndex()}` : `x${position(strikeXInput)}m`;
107+ const x = position(strikeXInput);
108+ strikeLabel = mode ? `mode${strikeIndex()}` : kindSelect.value === "mallet" ? `mallet${mallet.M}kg${mallet.K}-x${x}m` : `x${x}m`;
104109 // The hammer's half-width is one spacing of the beam at REF_N masses, so
105110 // the blow is the same at any N (at REF_N, on a mass, only that mass moves).
106- s = mode ? strikeMode(P, modes, strikeIndex(), E) : strikeAt(P, position(strikeXInput), base.a, E);
111+ strikeInfo.textContent = "";
112+ if (mode) s = strikeMode(P, modes, strikeIndex(), E);
113+ else if (kindSelect.value === "at") s = strikeAt(P, x, base.a, E);
114+ else {
115+ // The mallet's face presses on the same bump (as a weighted mean of the
116+ // displacements under it); the note starts at lift-off.
117+ const phi = hammerShape(P, x, base.a), sum = phi.reduce((a, b) => a + b, 0);
118+ const r = malletStrike(
119+ { N: P.N, m: P.m, forces: (u, F) => forces(P, u, F), weights: phi.map((v) => v / sum), step: stepFor(P, modes, E) },
120+ mallet, E,
121+ );
122+ s = { u: r.u, p: r.p, W: 0 };
123+ const share = energy(P, s) / E;
124+ strikeInfo.textContent =
125+ `mallet: contact ${(1000 * r.contactTime).toPrecision(3)} ms, ${(100 * share).toFixed(0)}% of the energy into the beam` +
126+ `${share < 0.995 ? " (the rest leaves with the mallet)" : ""}, peak force ${r.peakForce.toPrecision(3)} N; shown from lift-off`;
127+ }
107128 lastStrike = copyState(s);
108129 verlet = new Verlet(P, s);
109130 E0 = energy(P, s);
@@ -231,10 +252,16 @@ const paramInputs = {
231252 };
232253 const stretchSelect = document.getElementById("p-stretch") as HTMLSelectElement;
233254 const paramStatus = document.getElementById("param-status")!;
255+const malletInputs = {
256+ M: document.getElementById("p-mallet-m") as HTMLInputElement,
257+ K: document.getElementById("p-mallet-k") as HTMLInputElement,
258+};
234259 function showParams() {
235260 for (const [key, el] of Object.entries(paramInputs)) el.value = String(base[key as keyof typeof paramInputs]);
236261 stretchSelect.value = P.stretch;
237262 endsSelect.value = P.ends;
263+ malletInputs.M.value = String(mallet.M);
264+ malletInputs.K.value = String(mallet.K);
238265 showDerived();
239266 }
240267 // Amplitude decay times: one for all modes (drag alone), or per mode.
@@ -496,14 +523,16 @@ nInput.addEventListener("change", () => {
496523 endsSelect.addEventListener("change", () => setParams({ ...base, ends: endsSelect.value as Ends }, P.N));
497524 document.getElementById("param-apply")!.addEventListener("click", () => {
498525 const v = Object.fromEntries(Object.entries(paramInputs).map(([key, el]) => [key, Number(el.value)]));
526+ const mv = { M: Number(malletInputs.M.value), K: Number(malletInputs.K.value) };
499527 const ok = Object.entries(v).every(
500528 ([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop", "k", "kappa", "rise"].includes(key) && x === 0)),
501- );
529+ ) && Object.values(mv).every((x) => Number.isFinite(x) && x > 0);
502530 if (!ok || !(v.k > 0 || v.kappa > 0)) {
503531 paramStatus.textContent = "parameters must be positive (τ, τ top, rise, k or κ may be 0, not both k and κ)";
504532 return;
505533 }
506534 paramStatus.textContent = "";
535+ mallet = mv;
507536 setParams({ ...base, ...v, stretch: stretchSelect.value as Stretch }, P.N);
508537 });
509538
src/common/mallet.tsadded+81−0View file
@@ -0,0 +1,81 @@
1+// A mallet strike: a mass M arriving with kinetic energy E presses into the
2+// struck system through a Hertz contact, force K delta^1.5 while delta > 0
3+// (delta: how far the mallet's face has passed the struck point), 0 once it
4+// lifts off. The contact phase is integrated here with a fine step, from
5+// first touch to lift-off, and the system's state at lift-off is returned:
6+// the note then starts from it. Some energy leaves with the rebounding
7+// mallet. Drag and other losses are left out during the contact, which is far
8+// shorter than their decay times.
9+//
10+// The struck point is a weighted mean of the displacements, w = sum c_j u_j
11+// with sum c_j = 1 (one mass, or a bump over several), and the contact force
12+// F acts on mass j as F c_j (from the contact energy K delta^2.5 / 2.5).
13+
14+export interface MalletParams {
15+ M: number; // kg, mallet mass
16+ K: number; // N/m^1.5, Hertz contact stiffness
17+}
18+
19+export interface Struck {
20+ N: number;
21+ m: number; // kg, each mass
22+ forces(u: Float64Array, F: Float64Array): void; // the system's own (conservative) forces
23+ weights: Float64Array; // c_j, summing to 1
24+ step: number; // s, the system's own integrator step at this energy
25+}
26+
27+export interface MalletResult {
28+ u: Float64Array;
29+ p: Float64Array;
30+ contactTime: number; // s, first touch to lift-off
31+ peakForce: number; // N
32+ malletEnergy: number; // J, kinetic energy the mallet takes away
33+}
34+
35+// Contact time of the same mallet on a rigid wall: 3.218 delta_max / v,
36+// delta_max = (5 M v^2 / 4 K)^(2/5). Sets the scale of the fine step.
37+export function rigidContactTime(mp: MalletParams, E: number): number {
38+ const v = Math.sqrt((2 * E) / mp.M);
39+ return (3.218 * ((5 * mp.M * v * v) / (4 * mp.K)) ** 0.4) / v;
40+}
41+
42+export function malletStrike(sys: Struck, mp: MalletParams, E: number): MalletResult {
43+ const { N, m, weights: c } = sys;
44+ const u = new Float64Array(N), p = new Float64Array(N), F = new Float64Array(N);
45+ let y = 0, vy = Math.sqrt((2 * E) / mp.M); // mallet face position and velocity
46+ // Fine step: resolve the contact on both sides, the mallet (rigid-wall
47+ // contact time) and the struck masses (effective mass 1 / sum c_j^2 / m on
48+ // the contact's stiffest, 1.5 K delta^0.5 at the rigid-wall depth).
49+ const tr = rigidContactTime(mp, E);
50+ const dMax = ((5 * mp.M * (2 * E / mp.M)) / (4 * mp.K)) ** 0.4;
51+ let inv = 0;
52+ for (let j = 0; j < N; j++) inv += (c[j] * c[j]) / m;
53+ const h = Math.min(sys.step, tr / 400, 0.05 * Math.sqrt(1 / inv / (1.5 * mp.K * Math.sqrt(dMax))));
54+ let t = 0, peak = 0, Fc = 0;
55+ // Total force on the masses and on the mallet at the current positions.
56+ const contact = () => {
57+ let w = 0;
58+ for (let j = 0; j < N; j++) w += c[j] * u[j];
59+ const d = y - w;
60+ Fc = d > 0 ? mp.K * d ** 1.5 : 0;
61+ sys.forces(u, F);
62+ for (let j = 0; j < N; j++) F[j] += Fc * c[j];
63+ return d;
64+ };
65+ contact();
66+ // Velocity Verlet over the contact; stop at lift-off (or after a long
67+ // time, for a mallet much heavier than what it strikes).
68+ for (let d = 0; t < 1000 * tr; ) {
69+ for (let j = 0; j < N; j++) p[j] += 0.5 * h * F[j];
70+ vy -= (0.5 * h * Fc) / mp.M;
71+ for (let j = 0; j < N; j++) u[j] += (h * p[j]) / m;
72+ y += h * vy;
73+ d = contact();
74+ for (let j = 0; j < N; j++) p[j] += 0.5 * h * F[j];
75+ vy -= (0.5 * h * Fc) / mp.M;
76+ t += h;
77+ peak = Math.max(peak, Fc);
78+ if (d <= 0 && t > h) break;
79+ }
80+ return { u, p, contactTime: t, peakForce: peak, malletEnergy: 0.5 * mp.M * vy * vy };
81+}
↑↓ move↵openescclose