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(
69 p: Pick<StringParams, 't60Low' | 'fLow' | 't60High' | 'fHigh'>,
70): { sigma0: number; sigma1: number; note?: string } {
71 if (!(p.t60Low > 0) && !(p.t60High > 0)) return { sigma0: 0, sigma1: 0 };
72 const ln1000 = 3 * Math.LN10;
73 if (!(p.t60High > 0) || !(p.fHigh > p.fLow)) return { sigma0: ln1000 / p.t60Low, sigma1: 0 };
74 if (!(p.t60Low > 0)) {
75 const w = 2 * Math.PI * p.fHigh;
76 return { sigma0: 0, sigma1: ln1000 / p.t60High / (w * w) };
77 }
78 const s1 = ln1000 / p.t60Low;
79 const s2 = ln1000 / p.t60High;
80 const w1 = 2 * Math.PI * p.fLow;
81 const w2 = 2 * Math.PI * p.fHigh;
82 let sigma1 = (s2 - s1) / (w2 * w2 - w1 * w1);
83 let sigma0 = s1 - sigma1 * w1 * w1;
84 let note: string | undefined;
85 if (sigma1 < 0) {
86 sigma1 = 0;
87 sigma0 = s1;
88 note = 'High-frequency T60 is longer than low-frequency T60; using frequency-independent damping.';
89 } else if (sigma0 < 0) {
90 sigma0 = 0;
91 sigma1 = s1 / (w1 * w1);
92 note = 'Decay times imply negative mass damping; clamped to zero.';
93 }
94 return { sigma0, sigma1, note };
95}
97// Linearized modes about the straight equilibrium. Transverse: each spring acts
98// with stiffness T/L (zero to a sliding ring, which follows the bead).
99// Longitudinal: stiffness k (to clamped or sliding walls; nothing to a free end).
100function linearModes(p: StringParams, count: number): { trans: number[]; long: number[] } {
101 const pinned = new Set(parsePins(p.pins, p.n));
102 const idx: number[] = [];
103 for (let i = 0; i < p.n; i++) if (!pinned.has(i)) idx.push(i);
104 const build = (kSpring: number, wallLeft: number, wallRight: number) => {
105 const d = new Float64Array(idx.length);
106 const e = new Float64Array(Math.max(0, idx.length - 1));
107 idx.forEach((i, r) => {
108 let diag = 0;
109 diag += i === 0 ? wallLeft : kSpring;
110 diag += i === p.n - 1 ? wallRight : kSpring;
111 d[r] = diag / p.mass;
112 if (r < idx.length - 1 && idx[r + 1] === i + 1) e[r] = -kSpring / p.mass;
113 });
114 return tridiagEigenvalues(d, e, count)
115 .filter((l) => l > 1e-9 * Math.max(...d, 1))
116 .map((l) => Math.sqrt(l) / (2 * Math.PI));
117 };
118 const kt = p.tension / p.spacing;
119 const tWall = (end: EndType) => (end === 'clamped' ? kt : 0);
120 const lWall = (end: EndType) => (end === 'free' ? 0 : p.k);
121 return {
122 trans: build(kt, tWall(p.leftEnd), tWall(p.rightEnd)),
123 long: build(p.k, lWall(p.leftEnd), lWall(p.rightEnd)),
124 };
125}
127export function derive(p: StringParams): Derived {
128 const warnings: string[] = [];
129 const restLength = p.spacing - p.tension / p.k;
130 if (!(restLength > 0)) warnings.push('Tension too high for this stiffness: rest length L0 = L - T/k is not positive.');
131 if (p.tension > 0 && (p.leftEnd === 'free' || p.rightEnd === 'free'))
132 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.');
133 if (p.tension < 0) warnings.push('Negative tension (compression): the straight string is unstable to buckling.');
134 const linearDensity = p.mass / p.spacing;
135 const { sigma0, sigma1, note } = dampingFromT60(p);
136 if (note) warnings.push(note);
137 const kt = p.tension > 0 ? p.tension / p.spacing : p.k;
138 const modes = linearModes(p, 24);
139 return {
140 restLength,
141 prestretch: p.spacing / restLength - 1,
142 length: (p.n + 1) * p.spacing,
143 linearDensity,
144 waveSpeedT: Math.sqrt(Math.max(p.tension, 0) / linearDensity),
145 waveSpeedL: Math.sqrt((p.k * p.spacing) / linearDensity),
146 transverseHz: modes.trans,
147 longitudinalHz: modes.long,
148 sigma0,
149 sigma1,
150 dashpot: 2 * sigma1 * kt,
151 warnings,
152 };
153}
155export function buildSystem(p: StringParams, d: Derived): BeadSystem {
156 const n = p.n;
157 const rest = new Float64Array(3 * n);
158 for (let i = 0; i < n; i++) rest[3 * i] = (i + 1) * p.spacing;
159 const mass = new Float64Array(n).fill(p.mass);
160 const pinned = new Uint8Array(n);
161 for (const i of parsePins(p.pins, n)) pinned[i] = 1;
162 const anchors: Anchor[] = [];
163 const si: number[] = [];
164 const sj: number[] = [];
165 const addEnd = (end: EndType, x: number, bead: number) => {
166 if (end === 'free') return;
167 anchors.push({ pos: [x, 0, 0], slideBead: end === 'sliding' ? bead : -1 });
168 const node = n + anchors.length - 1;
169 // anchor on the left of bead 0 is si; keep orientation left -> right
170 if (bead === 0) {
171 si.push(node);
172 sj.push(0);
173 } else {
174 si.push(bead);
175 sj.push(node);
176 }
177 };
178 addEnd(p.leftEnd, 0, 0);
179 for (let i = 0; i < n - 1; i++) {
180 si.push(i);
181 sj.push(i + 1);
182 }
183 addEnd(p.rightEnd, (n + 1) * p.spacing, n - 1);
184 const ns = si.length;
185 return {
186 n,
187 rest,
188 mass,
189 pinned,
190 anchors,
191 si: Int32Array.from(si),
192 sj: Int32Array.from(sj),
193 sk: new Float64Array(ns).fill(p.k),
194 sL0: new Float64Array(ns).fill(d.restLength),
195 sc: new Float64Array(ns).fill(d.dashpot),
196 sigma0: d.sigma0,
197 };
198}
200function gaussianWeights(n: number, center: number, width: number): Float64Array {
201 const w = new Float64Array(n);
202 if (!(width > 0)) {
203 w[center] = 1;
204 return w;
205 }
206 for (let i = 0; i < n; i++) w[i] = Math.exp(-0.5 * ((i - center) / width) ** 2);
207 return w;
208}
210const clampIndex = (i: number, n: number) => Math.min(n - 1, Math.max(0, Math.round(i) - 1));
212export function buildExcitation(p: StringParams): Excitation {
213 const n = p.n;
214 const u0 = new Float64Array(3 * n);
215 const v0 = new Float64Array(3 * n);
216 const th = (p.excAngle * Math.PI) / 180;
217 const dy = Math.cos(th);
218 const dz = Math.sin(th);
219 const c = clampIndex(p.excBead, n);
220 if (p.excitation === 'pluck') {
221 // triangle from the nearest support on each side (pinned bead or wall) to the plucked bead
222 const pins = new Set(parsePins(p.pins, n));
223 let left = -1;
224 for (let i = c - 1; i >= 0; i--) if (pins.has(i)) { left = i; break; }
225 let right = n;
226 for (let i = c + 1; i < n; i++) if (pins.has(i)) { right = i; break; }
227 for (let i = left + 1; i < right; i++) {
228 const s = i <= c ? (i - left) / (c - left) : (right - i) / (right - c);
229 u0[3 * i + 1] = p.pluckAmp * s * dy;
230 u0[3 * i + 2] = p.pluckAmp * s * dz;
231 }
232 return { u0, v0 };
233 }
234 const w = gaussianWeights(n, c, p.excWidth);
235 if (p.excitation === 'strike') {
236 for (let i = 0; i < n; i++) {
237 v0[3 * i + 1] = p.strikeVel * w[i] * dy;
238 v0[3 * i + 2] = p.strikeVel * w[i] * dz;
239 }
240 return { u0, v0 };
241 }
242 let sum = 0;
243 for (let i = 0; i < n; i++) sum += w[i];
244 const shape = new Float64Array(3 * n);
245 for (let i = 0; i < n; i++) {
246 shape[3 * i + 1] = (w[i] / sum) * dy;
247 shape[3 * i + 2] = (w[i] / sum) * dz;
248 }
249 return { u0, v0, force: { shape, F0: p.hammerForce, tau: p.hammerMs / 1000 } };
250}
252export function buildPickup(p: StringParams): Float64Array {
253 const w = new Float64Array(3 * p.n);
254 const axis = p.pickupAxis === 'x' ? 0 : p.pickupAxis === 'y' ? 1 : 2;
255 if (p.pickupMode === 'sum') for (let i = 0; i < p.n; i++) w[3 * i + axis] = 1;
256 else w[3 * clampIndex(p.pickupBead, p.n) + axis] = 1;
257 return w;
258}