Beam: free-free ends in the physics (rigid motions removed by projection, modes without the rigid pair); not yet in the UI
3 changed files+125−48
src/beam/README.mdmodified+1−0View file
@@ -14,6 +14,7 @@ Live site: https://chaos-apps--listening-for-chaos.vault1-sites.magland.org/beam
1414 - 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.
1515 - 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.
1616 - 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).
17+- Free-free ends (physics in `beam.ts`, step 1 of 2, 2026-10-09; not yet in the Ends menu): no springs or joints beyond the end masses; masses at cell centres (j + 1/2) a, length N a (4 m at the defaults, against 4.5 m between two supports). Rigid motions (slide, tilt) are free, and the rails would make tension and stretching charge for a tilt, so energy and forces are evaluated at Pi u (the least-squares straight line removed) and the forces projected (`removeRigid`); the integrator re-projects u and p each step (removed energy into W). Modes leave out the two rigid ones (`Modes.count` = N - 2). Checked (Node): other ends bit-identical (48 cases); ratios at N = 64 2.755, 5.397, 8.911, 13.289 (free-free beam theory 2.757, 5.404, 8.933, 13.344); slide + tilt changes V by ~1e-15 relative; forces = -grad V to ~5e-11 with zero net force and torque; 200 T1 damped nonlinear runs keep momentum and angular momentum at round-off (1e-16, 1e-15), energy balance <= 8e-4. A kick at 1.0 m puts 22 % of its energy into rigid motion.
1718 - 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.
1819 - 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.
1920 - 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+117−43View file
@@ -21,11 +21,20 @@
2121 // quartic one, and adds linear stiffness k3 D0^2 per spring (local) or along
2222 // D0 as a whole (Kirchhoff).
2323 // Ends: "pinned" (both ends: u = 0, no moment), "clamped" (both: u = 0, zero
24-// slope) or "cantilever" (clamped at the left, free at the right). The
25-// supports are nodes -1 and N at u = 0. A pinned support has no joint term
26-// (its second difference, with the mirror ghost u[-2] = -u[0], is 0); a
27-// clamped one has a joint with ghost u[-2] = u[0], so c = 2 u[0]. A free end
28-// has no spring or joint beyond the last mass.
24+// slope), "cantilever" (clamped at the left, free at the right) or "free"
25+// (both free). The supports are nodes -1 and N at u = 0. A pinned support has
26+// no joint term (its second difference, with the mirror ghost u[-2] = -u[0],
27+// is 0); a clamped one has a joint with ghost u[-2] = u[0], so c = 2 u[0]. A
28+// free end has no spring or joint beyond its last mass.
29+// Free-free: nothing holds the beam, so sliding (all u equal) and tilting (u
30+// linear along the beam) are free motions, and the rails would make tension
31+// and stretching charge for a tilt. Every deformation is therefore measured
32+// from the best-fit rigid position: energy and forces at u' = Pi u, with Pi
33+// removing the least-squares straight line (slide and tilt), and the forces
34+// projected by Pi too, so they carry no net momentum or angular momentum.
35+// The integrator keeps u and p free of rigid parts (what it removes, from
36+// round-off or the loss, goes into W), and a strike's rigid part is removed
37+// by the caller (removeRigid): the beam is followed in its own frame.
2938 // Optional damping: linear drag on each mass proportional to its mass,
3039 // force -gamma m v = -gamma p, so every linear mode's amplitude decays as
3140 // exp(-t / tau) with tau = 2 / gamma.
@@ -40,7 +49,7 @@
4049 // State: displacements u, momenta p (kg m/s), and W, the energy (J)
4150 // dissipated so far (so E + W is conserved).
4251
43-export type Ends = "clamped" | "pinned" | "cantilever";
52+export type Ends = "clamped" | "pinned" | "cantilever" | "free";
4453 export type Stretch = "local" | "kirchhoff";
4554
4655 export interface Params {
@@ -68,7 +77,7 @@ export interface Params {
6877 // At N = REF_N the parameters come back unchanged.
6978 export const REF_N = 8;
7079 export function refine(base: Params, N: number): Params {
71- const span = (n: number) => (base.ends === "cantilever" ? n : n + 1);
80+ const span = (n: number) => (freeRight(base) ? n : n + 1);
7281 const r = span(base.N) / span(N);
7382 const P = { ...base, N, m: base.m * r, k: base.k / r, kappa: base.kappa / r ** 3, k3: base.k3 / r ** 3, a: base.a * r };
7483 // tauTop is set at the top of the band, which moves with N: keep c / r.
@@ -80,23 +89,68 @@ export function refine(base: Params, N: number): Params {
8089 export const dragRate = (P: Params) => (P.tau > 0 ? 2 / P.tau : 0);
8190 export const isDamped = (P: Params) => P.tau > 0 || P.tauTop > 0;
8291
83-const freeRight = (P: Params) => P.ends === "cantilever";
84-// Springs j = 0..nS-1 join nodes j - 1 and j.
85-export const springCount = (P: Params) => (freeRight(P) ? P.N : P.N + 1);
86-// Joints at masses 0..nJ-1 (each needs both neighbours).
87-const jointCount = (P: Params) => (freeRight(P) ? P.N - 1 : P.N);
92+const freeLeft = (P: Params) => P.ends === "free";
93+const freeRight = (P: Params) => P.ends === "cantilever" || P.ends === "free";
94+const clampedLeft = (P: Params) => P.ends === "clamped" || P.ends === "cantilever";
95+// Springs j = j0..j1-1 join nodes j - 1 and j (node -1 and N: the supports).
96+const springRange = (P: Params) => [freeLeft(P) ? 1 : 0, freeRight(P) ? P.N : P.N + 1];
97+export const springCount = (P: Params) => {
98+ const [j0, j1] = springRange(P);
99+ return j1 - j0;
100+};
101+// Joints at masses i0..i1-1 (each needs both neighbours).
102+const jointRange = (P: Params) => [freeLeft(P) ? 1 : 0, freeRight(P) ? P.N - 1 : P.N];
88103
89104 // Largest eigenvalue of L: 4 sin^2(N pi / (2 (N + 1))) between two supports;
90-// with a free end, 4 sin^2((2N - 1) pi / (2 (2N + 1))).
105+// with one free end, 4 sin^2((2N - 1) pi / (2 (2N + 1))); with two,
106+// 4 sin^2((N - 1) pi / (2N)).
91107 const laplacianMax = (P: Params) =>
92- 4 * Math.sin(freeRight(P) ? ((2 * P.N - 1) * Math.PI) / (2 * (2 * P.N + 1)) : (P.N * Math.PI) / (2 * (P.N + 1))) ** 2;
108+ 4 *
109+ Math.sin(
110+ P.ends === "free"
111+ ? ((P.N - 1) * Math.PI) / (2 * P.N)
112+ : freeRight(P)
113+ ? ((2 * P.N - 1) * Math.PI) / (2 * (2 * P.N + 1))
114+ : (P.N * Math.PI) / (2 * (P.N + 1)),
115+ ) ** 2;
93116 // Loss coefficient c (kg/s); 0 when off.
94117 const lossCoefficient = (P: Params) => (P.tauTop > 0 ? (2 * P.m) / (P.tauTop * laplacianMax(P)) : 0);
95118 // (L x) into out.
96119 function laplacian(P: Params, x: Float64Array, out: Float64Array): void {
97120 const N = P.N;
98121 for (let j = 0; j < N; j++) out[j] = 2 * x[j] - (j > 0 ? x[j - 1] : 0) - (j + 1 < N ? x[j + 1] : 0);
99- if (freeRight(P)) out[N - 1] -= x[N - 1]; // no spring beyond the free end
122+ if (freeRight(P)) out[N - 1] -= x[N - 1]; // no spring beyond a free end
123+ if (freeLeft(P)) out[0] -= x[0];
124+}
125+
126+// Rigid motions of a free-free beam: orthonormal slide e1 = 1 / sqrt(N) and
127+// tilt e2 ~ j - (N - 1) / 2. removeRigid takes their parts out of v (in
128+// place) and returns the squared norm removed.
129+export function removeRigid(v: Float64Array): number {
130+ const N = v.length, c = (N - 1) / 2;
131+ let s1 = 0, s2 = 0, n2 = 0;
132+ for (let j = 0; j < N; j++) {
133+ s1 += v[j];
134+ s2 += (j - c) * v[j];
135+ n2 += (j - c) ** 2;
136+ }
137+ const a1 = s1 / N, a2 = s2 / n2;
138+ let removed = 0;
139+ for (let j = 0; j < N; j++) {
140+ const r = a1 + a2 * (j - c);
141+ v[j] -= r;
142+ removed += r * r;
143+ }
144+ return removed;
145+}
146+// Scratch for the projected displacements (forces runs every step).
147+let rigidScratch = new Float64Array(0);
148+function deformed(P: Params, u: Float64Array): Float64Array {
149+ if (P.ends !== "free") return u;
150+ if (rigidScratch.length !== u.length) rigidScratch = new Float64Array(u.length);
151+ rigidScratch.set(u);
152+ removeRigid(rigidScratch);
153+ return rigidScratch;
100154 }
101155
102156 export interface State {
@@ -114,24 +168,28 @@ export const copyState = (s: State): State => ({ u: Float64Array.from(s.u), p: F
114168 // Displacement of node i: a mass, or a support (0).
115169 const at = (u: Float64Array, i: number) => (i >= 0 && i < u.length ? u[i] : 0);
116170
117-// Rest shape of mass j (m), with x its distance from the left support in
118-// spacings and L the span: between two supports a half sine, rise sin(pi x / L)
119-// with L = N + 1; for a cantilever 1 - cos(pi x / (2 L)) times the rise with
120-// L = N, flat at the clamp and highest at the free end.
171+// Rest shape of mass j (m), with x its distance from the left support (or
172+// end) in spacings and L the span: between two supports a half sine, rise
173+// sin(pi x / L) with L = N + 1; for a cantilever 1 - cos(pi x / (2 L)) times
174+// the rise with L = N, flat at the clamp and highest at the free end; free-
175+// free a half sine over the length, x = j + 1/2, L = N (an arch).
121176 export function restShape(P: Params, j: number): number {
122177 if (P.rise === 0) return 0;
178+ if (P.ends === "free") return P.rise * Math.sin((Math.PI * (j + 0.5)) / P.N);
123179 const x = j + 1;
124180 return freeRight(P) ? P.rise * (1 - Math.cos((Math.PI * x) / (2 * P.N))) : P.rise * Math.sin((Math.PI * x) / (P.N + 1));
125181 }
126182
127-// Rest offsets D0 of the springs (spring j joins nodes j - 1 and j), kept for
128-// the last parameters seen (forces runs every step).
183+// Rest offsets D0 of the springs (spring j joins nodes j - 1 and j; indexed
184+// by j, 0..N, unused entries 0), kept for the last parameters seen (forces
185+// runs every step).
129186 let restKey: { N: number; ends: Ends; rise: number } | null = null;
130187 let restD0 = new Float64Array(0);
131188 function restOffsets(P: Params): Float64Array {
132189 if (!restKey || restKey.N !== P.N || restKey.ends !== P.ends || restKey.rise !== P.rise) {
133190 const w0 = (i: number) => (i >= 0 && i < P.N ? restShape(P, i) : 0);
134- restD0 = Float64Array.from({ length: springCount(P) }, (_, j) => w0(j) - w0(j - 1));
191+ const [j0, j1] = springRange(P);
192+ restD0 = Float64Array.from({ length: P.N + 1 }, (_, j) => (j >= j0 && j < j1 ? w0(j) - w0(j - 1) : 0));
135193 restKey = { N: P.N, ends: P.ends, rise: P.rise };
136194 }
137195 return restD0;
@@ -139,10 +197,11 @@ function restOffsets(P: Params): Float64Array {
139197
140198 // Force on each mass, written into F. With linear set, only the part linear
141199 // in u, -K u (as with k3 = 0 for a straight beam).
142-export function forces(P: Params, u: Float64Array, F: Float64Array, linear = false): void {
200+export function forces(P: Params, u0: Float64Array, F: Float64Array, linear = false): void {
143201 const { N, k, kappa, k3 } = P;
144202 F.fill(0);
145- const nS = springCount(P);
203+ const u = deformed(P, u0);
204+ const nS = springCount(P), [j0, j1] = springRange(P), [i0, i1] = jointRange(P);
146205 const D0 = restOffsets(P);
147206 // Stretch e of each spring: d^2 + 2 D0 d, or its linear part 2 D0 d.
148207 const stretch = (j: number, d: number) => (linear ? 2 * D0[j] * d : d * d + 2 * D0[j] * d);
@@ -150,10 +209,10 @@ export function forces(P: Params, u: Float64Array, F: Float64Array, linear = fal
150209 let kK = 0;
151210 if (P.stretch === "kirchhoff") {
152211 let S = 0;
153- for (let j = 0; j < nS; j++) S += stretch(j, at(u, j) - at(u, j - 1));
212+ for (let j = j0; j < j1; j++) S += stretch(j, at(u, j) - at(u, j - 1));
154213 kK = (k3 * S) / nS;
155214 }
156- for (let j = 0; j < nS; j++) {
215+ for (let j = j0; j < j1; j++) {
157216 const d = at(u, j) - at(u, j - 1);
158217 // Spring j pulls its ends together with force dV/dd = (k + k3 e) d + k3 e D0
159218 // (local; Kirchhoff: kK for k3 e); the linear part drops k3 e d.
@@ -163,16 +222,17 @@ export function forces(P: Params, u: Float64Array, F: Float64Array, linear = fal
163222 if (j < N) F[j] -= f;
164223 }
165224 if (kappa > 0) {
166- for (let i = 0, nJ = jointCount(P); i < nJ; i++) {
225+ for (let i = i0; i < i1; i++) {
167226 const c = kappa * (at(u, i - 1) - 2 * u[i] + at(u, i + 1));
168227 if (i > 0) F[i - 1] -= c;
169228 F[i] += 2 * c;
170229 if (i + 1 < N) F[i + 1] -= c;
171230 }
172231 // Clamped supports: c = 2 u at the end mass, energy 2 kappa u^2.
173- if (P.ends !== "pinned") F[0] -= 4 * kappa * u[0];
232+ if (clampedLeft(P)) F[0] -= 4 * kappa * u[0];
174233 if (P.ends === "clamped") F[N - 1] -= 4 * kappa * u[N - 1];
175234 }
235+ if (P.ends === "free") removeRigid(F);
176236 }
177237
178238 export function kinetic(P: Params, s: State): number {
@@ -183,11 +243,11 @@ export function kinetic(P: Params, s: State): number {
183243
184244 export function potential(P: Params, s: State): number {
185245 const { N, k, kappa, k3 } = P;
186- const u = s.u;
187- const nS = springCount(P);
246+ const u = deformed(P, s.u);
247+ const nS = springCount(P), [j0, j1] = springRange(P), [i0, i1] = jointRange(P);
188248 const D0 = restOffsets(P);
189249 let V = 0, S = 0;
190- for (let j = 0; j < nS; j++) {
250+ for (let j = j0; j < j1; j++) {
191251 const d = at(u, j) - at(u, j - 1);
192252 const d2 = d ** 2;
193253 const e = d2 + 2 * D0[j] * d; // stretch from the rest shape
@@ -195,8 +255,8 @@ export function potential(P: Params, s: State): number {
195255 S += e;
196256 }
197257 if (P.stretch === "kirchhoff") V += (k3 * S * S) / (4 * nS);
198- for (let i = 0, nJ = jointCount(P); i < nJ; i++) V += (kappa * (at(u, i - 1) - 2 * u[i] + at(u, i + 1)) ** 2) / 2;
199- if (P.ends !== "pinned") V += 2 * kappa * u[0] ** 2;
258+ for (let i = i0; i < i1; i++) V += (kappa * (at(u, i - 1) - 2 * u[i] + at(u, i + 1)) ** 2) / 2;
259+ if (clampedLeft(P)) V += 2 * kappa * u[0] ** 2;
200260 if (P.ends === "clamped") V += 2 * kappa * u[N - 1] ** 2;
201261 return V;
202262 }
@@ -247,6 +307,11 @@ export class Verlet {
247307 for (let j = 0; j < N; j++) p[j] += 0.5 * h * F[j];
248308 if (this.c > 0) this.loss(s, h / 2);
249309 if (gamma > 0) this.drag(s, Math.exp(-0.5 * gamma * h));
310+ if (this.P.ends === "free") {
311+ // Keep the beam in its own frame (the loss's end dashpots, round-off).
312+ removeRigid(u);
313+ s.W += removeRigid(p) / (2 * m);
314+ }
250315 }
251316 // (I + a L) p' = (I - a L) p, a = c dt / 2m, by the Thomas algorithm.
252317 private loss(s: State, dt: number): void {
@@ -258,7 +323,8 @@ export class Verlet {
258323 laplacian(this.P, p, rhs);
259324 for (let j = 0; j < N; j++) rhs[j] = p[j] - a * rhs[j];
260325 // Matrix I + a L: diagonal 1 + 2a (1 + a at a free end), off-diagonals -a.
261- const diag = (j: number) => 1 + (freeRight(this.P) && j === N - 1 ? a : 2 * a);
326+ const fl = freeLeft(this.P), fr = freeRight(this.P);
327+ const diag = (j: number) => 1 + ((fr && j === N - 1) || (fl && j === 0) ? a : 2 * a);
262328 let b = diag(0);
263329 cp[0] = -a / b;
264330 p[0] = rhs[0] / b;
@@ -326,12 +392,15 @@ function jacobi(A: Float64Array[]): { values: Float64Array; V: Float64Array[] }
326392 }
327393
328394 // Normal modes of the beam at small amplitude (the linear part of the forces,
329-// which includes the stiffness the curvature adds), mode n = 1..N stored at
330-// index n - 1 in order of frequency, from the stiffness matrix K (assembled
395+// which includes the stiffness the curvature adds), mode n = 1..count stored
396+// at index n - 1 in order of frequency, from the stiffness matrix K (assembled
331397 // from forces(), so it matches the dynamics exactly). Shapes S[n - 1] have unit
332398 // norm and start upward at the left; mode amplitudes are Q = S u, P = S p.
399+// count = N, except free-free: N - 2, the two rigid motions (zero frequency,
400+// removed from the dynamics) left out, so mode 1 is the first bending mode.
333401 export class Modes {
334402 readonly N: number;
403+ readonly count: number;
335404 readonly m: number;
336405 readonly omega: Float64Array; // rad/s
337406 readonly S: Float64Array[]; // S[n - 1][j]
@@ -349,9 +418,11 @@ export class Modes {
349418 K.push(F.map((f) => -f)); // column i of K (K is symmetric)
350419 }
351420 const { values, V } = jacobi(K);
352- const order = Array.from(values.keys()).sort((x, y) => values[x] - values[y]);
421+ const rigid = P.ends === "free" ? 2 : 0;
422+ const order = Array.from(values.keys()).sort((x, y) => values[x] - values[y]).slice(rigid);
423+ this.count = N - rigid;
353424 this.omega = Float64Array.from(order, (i) => Math.sqrt(Math.max(0, values[i]) / m));
354- this.Kmax = values[order[N - 1]];
425+ this.Kmax = values[order[this.count - 1]];
355426 this.S = order.map((i) => {
356427 const row = Float64Array.from({ length: N }, (_, j) => V[j][i]);
357428 const first = row.find((x) => Math.abs(x) > 1e-9) ?? 1;
@@ -363,7 +434,7 @@ export class Modes {
363434 // out. Their sum is the energy without the nonlinear part of the stretching.
364435 energies(s: State, out: Float64Array): void {
365436 const m = this.m;
366- for (let n = 0; n < this.N; n++) {
437+ for (let n = 0; n < this.count; n++) {
367438 const row = this.S[n];
368439 let Q = 0, Pn = 0;
369440 for (let j = 0; j < this.N; j++) {
@@ -442,10 +513,13 @@ export function strikeMode(P: Params, modes: Modes, n: number, E: number): State
442513 return s;
443514 }
444515
445-// Positions along the beam (m from the left support): mass j (0-based) is
446-// at (j + 1) a; the beam's length is its span in spacings times a.
447-export const massPosition = (P: Params, j: number) => (j + 1) * P.a;
516+// Positions along the beam (m from the left support or end): mass j
517+// (0-based) is at (j + 1) a, free-free at (j + 1/2) a (cell centres); the
518+// beam's length is its span in spacings times a. massIndex: the fractional
519+// mass index at x.
520+export const massPosition = (P: Params, j: number) => (j + (freeLeft(P) ? 0.5 : 1)) * P.a;
448521 export const beamLength = (P: Params) => (freeRight(P) ? P.N : P.N + 1) * P.a;
522+export const massIndex = (P: Params, x: number) => x / P.a - (freeLeft(P) ? 0.5 : 1);
449523
450524 // The hammer's shape: a triangular bump 1 - |x_j - x| / w on the masses it
451525 // covers, x (m) its centre and w (m) its half-width (at least one spacing);
@@ -458,7 +532,7 @@ export function hammerShape(P: Params, x: number, w: number): Float64Array {
458532 phi[j] = Math.max(0, 1 - Math.abs(massPosition(P, j) - x) / we);
459533 if (phi[j] > 0) any = true;
460534 }
461- if (!any) phi[Math.min(P.N - 1, Math.max(0, Math.round(x / P.a) - 1))] = 1;
535+ if (!any) phi[Math.min(P.N - 1, Math.max(0, Math.round(massIndex(P, x))))] = 1;
462536 return phi;
463537 }
464538
src/beam/signals.tsmodified+7−5View file
@@ -4,7 +4,7 @@
44 // Point signals are read at a position x along the beam (signalAt), so the
55 // same listening point serves every N.
66
7-import { beamLength, type Params, type State } from "./beam";
7+import { beamLength, massIndex, type Params, type State } from "./beam";
88
99 export interface Signal {
1010 key: string;
@@ -35,12 +35,14 @@ export const signalByKey = (key: string) => SIGNALS.find((s) => s.key === key) ?
3535
3636 // A signal at x (m from the left support), linear between the two nearest
3737 // masses; a support counts as 0 (it does not move), and past the last mass
38-// of a cantilever the last mass's value holds. Exactly the mass's value when
39-// x is on a mass.
38+// at a free end that mass's value holds. Exactly the mass's value when x is
39+// on a mass.
4040 export function signalAt(sig: Signal, P: Params, s: State, F: Float64Array, x: number): number {
4141 if (!sig.perMass) return sig.value(P, s, F, 0);
42- const at = (j: number) => (j < 0 ? 0 : j >= P.N ? (P.ends === "cantilever" ? sig.value(P, s, F, P.N - 1) : 0) : sig.value(P, s, F, j));
43- const f = Math.min(Math.max(x, 0), beamLength(P)) / P.a - 1; // position in mass indices
42+ const fl = P.ends === "free", fr = P.ends === "cantilever" || fl;
43+ const at = (j: number) =>
44+ j < 0 ? (fl ? sig.value(P, s, F, 0) : 0) : j >= P.N ? (fr ? sig.value(P, s, F, P.N - 1) : 0) : sig.value(P, s, F, j);
45+ const f = massIndex(P, Math.min(Math.max(x, 0), beamLength(P))); // position in mass indices
4446 const j0 = Math.floor(f), frac = f - j0;
4547 return frac === 0 ? at(j0) : (1 - frac) * at(j0) + frac * at(j0 + 1);
4648 }