4a8abcfString of beads: 3D bead-spring solver, control panel, waveform, spectrogram, wav, displacement viewerJeremy Magland 1import type { Anchor, BeadSystem, Excitation } from '../../sim/types';
2import { tridiagEigenvalues } from '../../sim/tridiag';
4export type EndType = 'clamped' | 'sliding' | 'free';
5export type ExcitationType = 'pluck' | 'strike' | 'hammer';
7export interface StringParams {
8 // geometry
9 n: number; // beads
10 spacing: number; // m, bead spacing at equilibrium (also bead-to-wall)
11 // properties
12 mass: number; // kg per bead
13 k: number; // spring stiffness, N/m
14 tension: number; // N, equilibrium tension (sets rest length L0 = spacing - T/k)
15 leftEnd: EndType;
16 rightEnd: EndType;
17 pins: string; // 1-based bead numbers, comma separated
18 // damping (0 = off)
19 t60Low: number; // s
20 fLow: number; // Hz
21 t60High: number; // s
22 fHigh: number; // Hz
23 // excitation
24 excitation: ExcitationType;
25 excBead: number; // 1-based
26 excWidth: number; // beads (Gaussian sigma); 0 = one bead
27 pluckAmp: number; // m
28 strikeVel: number; // m/s
29 hammerForce: number; // N, peak total
30 hammerMs: number; // ms
31 excAngle: number; // degrees from y toward z
32 // pickup
33 pickupMode: 'bead' | 'sum';
34 pickupBead: number; // 1-based
35 pickupAxis: 'y' | 'z' | 'x';
36 // simulation
37 duration: number;
38 fs: number;
39 cfl: number;
40 frameRate: number;
41}
43export function parsePins(s: string, n: number): number[] {
44 return s
45 .split(/[\s,;]+/)
46 .filter((t) => t.length > 0)
47 .map((t) => Math.round(Number(t)) - 1)
48 .filter((i) => Number.isFinite(i) && i >= 0 && i < n);
49}
51export interface Derived {
52 restLength: number;
53 prestretch: number; // (L - L0) / L0
54 length: number; // wall to wall (m)
55 linearDensity: number; // kg/m
56 waveSpeedT: number; // transverse, continuum limit (m/s)
57 waveSpeedL: number; // longitudinal, continuum limit (m/s)
58 transverseHz: number[]; // linear modes (lowest first)
59 longitudinalHz: number[];
60 sigma0: number;
61 sigma1: number;
62 dashpot: number;
63 warnings: string[];
64}
66// Linear damping sigma(omega) = sigma0 + sigma1 omega^2 for transverse modes,
67// fitted to two decay times T60 = 3 ln(10) / sigma.
68export function dampingFromT60(p: StringParams): { sigma0: number; sigma1: number; note?: string } {
69 if (!(p.t60Low > 0) && !(p.t60High > 0)) return { sigma0: 0, sigma1: 0 };
70 const ln1000 = 3 * Math.LN10;
71 if (!(p.t60High > 0) || !(p.fHigh > p.fLow)) return { sigma0: ln1000 / p.t60Low, sigma1: 0 };
72 if (!(p.t60Low > 0)) {
73 const w = 2 * Math.PI * p.fHigh;
74 return { sigma0: 0, sigma1: ln1000 / p.t60High / (w * w) };
75 }
76 const s1 = ln1000 / p.t60Low;
77 const s2 = ln1000 / p.t60High;
78 const w1 = 2 * Math.PI * p.fLow;
79 const w2 = 2 * Math.PI * p.fHigh;
80 let sigma1 = (s2 - s1) / (w2 * w2 - w1 * w1);
81 let sigma0 = s1 - sigma1 * w1 * w1;
82 let note: string | undefined;
83 if (sigma1 < 0) {
84 sigma1 = 0;
85 sigma0 = s1;
86 note = 'High-frequency T60 is longer than low-frequency T60; using frequency-independent damping.';
87 } else if (sigma0 < 0) {
88 sigma0 = 0;
89 sigma1 = s1 / (w1 * w1);
90 note = 'Decay times imply negative mass damping; clamped to zero.';
91 }
92 return { sigma0, sigma1, note };
93}
95// Linearized modes about the straight equilibrium. Transverse: each spring acts
96// with stiffness T/L (zero to a sliding ring, which follows the bead).
97// Longitudinal: stiffness k (to clamped or sliding walls; nothing to a free end).
98function linearModes(p: StringParams, count: number): { trans: number[]; long: number[] } {
99 const pinned = new Set(parsePins(p.pins, p.n));
100 const idx: number[] = [];
101 for (let i = 0; i < p.n; i++) if (!pinned.has(i)) idx.push(i);
102 const build = (kSpring: number, wallLeft: number, wallRight: number) => {
103 const d = new Float64Array(idx.length);
104 const e = new Float64Array(Math.max(0, idx.length - 1));
105 idx.forEach((i, r) => {
106 let diag = 0;
107 diag += i === 0 ? wallLeft : kSpring;
108 diag += i === p.n - 1 ? wallRight : kSpring;
109 d[r] = diag / p.mass;
110 if (r < idx.length - 1 && idx[r + 1] === i + 1) e[r] = -kSpring / p.mass;
111 });
112 return tridiagEigenvalues(d, e, count)
113 .filter((l) => l > 1e-9 * Math.max(...d, 1))
114 .map((l) => Math.sqrt(l) / (2 * Math.PI));
115 };
116 const kt = p.tension / p.spacing;
117 const tWall = (end: EndType) => (end === 'clamped' ? kt : 0);
118 const lWall = (end: EndType) => (end === 'free' ? 0 : p.k);
119 return {
120 trans: build(kt, tWall(p.leftEnd), tWall(p.rightEnd)),
121 long: build(p.k, lWall(p.leftEnd), lWall(p.rightEnd)),
122 };
123}
125export function derive(p: StringParams): Derived {
126 const warnings: string[] = [];
127 const restLength = p.spacing - p.tension / p.k;
128 if (!(restLength > 0)) warnings.push('Tension too high for this stiffness: rest length L0 = L - T/k is not positive.');
129 if (p.tension > 0 && (p.leftEnd === 'free' || p.rightEnd === 'free'))
130 warnings.push('A free end cannot hold tension: the string is not in equilibrium and will contract. Use sliding ends (tension kept, free transverse motion) or set T = 0.');
131 if (p.tension < 0) warnings.push('Negative tension (compression): the straight string is unstable to buckling.');
132 const linearDensity = p.mass / p.spacing;
133 const { sigma0, sigma1, note } = dampingFromT60(p);
134 if (note) warnings.push(note);
135 const kt = p.tension > 0 ? p.tension / p.spacing : p.k;
136 const modes = linearModes(p, 24);
137 return {
138 restLength,
139 prestretch: p.spacing / restLength - 1,
140 length: (p.n + 1) * p.spacing,
141 linearDensity,
142 waveSpeedT: Math.sqrt(Math.max(p.tension, 0) / linearDensity),
143 waveSpeedL: Math.sqrt((p.k * p.spacing) / linearDensity),
144 transverseHz: modes.trans,
145 longitudinalHz: modes.long,
146 sigma0,
147 sigma1,
148 dashpot: 2 * sigma1 * kt,
149 warnings,
150 };
151}
153export function buildSystem(p: StringParams, d: Derived): BeadSystem {
154 const n = p.n;
155 const rest = new Float64Array(3 * n);
156 for (let i = 0; i < n; i++) rest[3 * i] = (i + 1) * p.spacing;
157 const mass = new Float64Array(n).fill(p.mass);
158 const pinned = new Uint8Array(n);
159 for (const i of parsePins(p.pins, n)) pinned[i] = 1;
160 const anchors: Anchor[] = [];
161 const si: number[] = [];
162 const sj: number[] = [];
163 const addEnd = (end: EndType, x: number, bead: number) => {
164 if (end === 'free') return;
165 anchors.push({ pos: [x, 0, 0], slideBead: end === 'sliding' ? bead : -1 });
166 const node = n + anchors.length - 1;
167 // anchor on the left of bead 0 is si; keep orientation left -> right
168 if (bead === 0) {
169 si.push(node);
170 sj.push(0);
171 } else {
172 si.push(bead);
173 sj.push(node);
174 }
175 };
176 addEnd(p.leftEnd, 0, 0);
177 for (let i = 0; i < n - 1; i++) {
178 si.push(i);
179 sj.push(i + 1);
180 }
181 addEnd(p.rightEnd, (n + 1) * p.spacing, n - 1);
182 const ns = si.length;
183 return {
184 n,
185 rest,
186 mass,
187 pinned,
188 anchors,
189 si: Int32Array.from(si),
190 sj: Int32Array.from(sj),
191 sk: new Float64Array(ns).fill(p.k),
192 sL0: new Float64Array(ns).fill(d.restLength),
193 sc: new Float64Array(ns).fill(d.dashpot),
194 sigma0: d.sigma0,
195 };
196}
198function gaussianWeights(n: number, center: number, width: number): Float64Array {
199 const w = new Float64Array(n);
200 if (!(width > 0)) {
201 w[center] = 1;
202 return w;
203 }
204 for (let i = 0; i < n; i++) w[i] = Math.exp(-0.5 * ((i - center) / width) ** 2);
205 return w;
206}
208const clampIndex = (i: number, n: number) => Math.min(n - 1, Math.max(0, Math.round(i) - 1));
210export function buildExcitation(p: StringParams): Excitation {
211 const n = p.n;
212 const u0 = new Float64Array(3 * n);
213 const v0 = new Float64Array(3 * n);
214 const th = (p.excAngle * Math.PI) / 180;
215 const dy = Math.cos(th);
216 const dz = Math.sin(th);
217 const c = clampIndex(p.excBead, n);
218 if (p.excitation === 'pluck') {
219 // triangle from the nearest support on each side (pinned bead or wall) to the plucked bead
220 const pins = new Set(parsePins(p.pins, n));
221 let left = -1;
222 for (let i = c - 1; i >= 0; i--) if (pins.has(i)) { left = i; break; }
223 let right = n;
224 for (let i = c + 1; i < n; i++) if (pins.has(i)) { right = i; break; }
225 for (let i = left + 1; i < right; i++) {
226 const s = i <= c ? (i - left) / (c - left) : (right - i) / (right - c);
227 u0[3 * i + 1] = p.pluckAmp * s * dy;
228 u0[3 * i + 2] = p.pluckAmp * s * dz;
229 }
230 return { u0, v0 };
231 }
232 const w = gaussianWeights(n, c, p.excWidth);
233 if (p.excitation === 'strike') {
234 for (let i = 0; i < n; i++) {
235 v0[3 * i + 1] = p.strikeVel * w[i] * dy;
236 v0[3 * i + 2] = p.strikeVel * w[i] * dz;
237 }
238 return { u0, v0 };
239 }
240 let sum = 0;
241 for (let i = 0; i < n; i++) sum += w[i];
242 const shape = new Float64Array(3 * n);
243 for (let i = 0; i < n; i++) {
244 shape[3 * i + 1] = (w[i] / sum) * dy;
245 shape[3 * i + 2] = (w[i] / sum) * dz;
246 }
247 return { u0, v0, force: { shape, F0: p.hammerForce, tau: p.hammerMs / 1000 } };
248}
250export function buildPickup(p: StringParams): Float64Array {
251 const w = new Float64Array(3 * p.n);
252 const axis = p.pickupAxis === 'x' ? 0 : p.pickupAxis === 'y' ? 1 : 2;
253 if (p.pickupMode === 'sum') for (let i = 0; i < p.n; i++) w[3 * i + axis] = 1;
254 else w[3 * clampIndex(p.pickupBead, p.n) + axis] = 1;
255 return w;
256}