/ proof-of-concept / toy-instruments
Sign in
proof-of-concept / toy-instruments
15fbb00
Branches
main
toy-instruments / src / instruments / bead-cymbal / model.ts
445 lines · 16.3 KBCodeBlameHistory
15fbb00Cymbal: replace felt washers with a minimal rigid hub (pinned beads around the hole), no support dampingJeremy Magland 1import type { BeadSystem, Bending, Excitation } from '../../sim/types';
3import { dampingFromT60 } from '../bead-string/model';
5// A cymbal: a thin bronze plate, round, free at the rim, with a hole in the
15fbb00Cymbal: replace felt washers with a minimal rigid hub (pinned beads around the hole), no support dampingJeremy Magland 6// centre; the beads around the hole are held fixed (a small rigid hub), so the
7// rest of the plate rings like a free plate, with no added damping. Beads on a triangular lattice in
c696b9aCymbal: plate bending via per-triangle shape operators, felt centre support, stick and mallet with angle and friction, far-field microphoneJeremy Magland 8// the x-z plane, lifted to the profile height y(r) (flat by default; optional
9// bow and bell). In-plane springs on every edge; plate bending from the
10// dihedral angles across interior edges, per triangle (see Bending in sim/types). Vibration is mostly along y.
12export type CymbalExcitation = 'stick' | 'mallet' | 'strike' | 'force';
14export interface CymbalParams {
15 // geometry
16 diameter: number; // m
17 rings: number; // beads from the centre to the rim
18 holeRadius: number; // m
19 bowHeight: number; // m, dome rise from rim to centre (0 = flat)
20 bellDiameter: number; // m
21 bellHeight: number; // m (0 = no bell)
22 // material
23 thickness: number; // m
24 youngs: number; // Pa
25 density: number; // kg/m^3
26 soften: number; // in-plane stiffness divided by this (1 = real bronze)
15fbb00Cymbal: replace felt washers with a minimal rigid hub (pinned beads around the hole), no support dampingJeremy Magland 27 // centre support: beads within this radius are held fixed (a rigid hub)
28 hubRadius: number; // m
30 t60Low: number;
31 fLow: number;
32 t60High: number;
33 fHigh: number;
34 // excitation
35 excitation: CymbalExcitation;
36 excR: number; // fraction of radius
37 excAngle: number; // degrees, azimuth in the x-z plane
38 hitAngle: number; // degrees from the surface normal, tilted toward the rim
39 excWidth: number; // m, Gaussian sigma of the patch (0 = one bead)
40 stickMass: number;
41 stickSpeed: number;
42 stickK: number; // N/m^1.5 (Hertz)
43 malletMass: number;
44 malletSpeed: number;
45 malletK: number; // N/m^p
46 malletP: number;
47 friction: number;
48 strikeVel: number;
49 forcePeak: number;
50 forceMs: number;
51 // microphone
52 micDistance: number; // m
53 micElevation: number; // degrees above the plane
54 micAzimuth: number; // degrees
55 // simulation
56 duration: number;
57 fs: number;
58 cfl: number;
59 frameRate: number;
62export function profileHeight(p: CymbalParams, r: number): number {
63 const R = p.diameter / 2;
64 let y = p.bowHeight * (1 - (r / R) ** 2);
65 const rb = p.bellDiameter / 2;
66 if (p.bellHeight > 0 && rb > 0 && r < rb) y += p.bellHeight * 0.5 * (1 + Math.cos((Math.PI * r) / rb));
67 return y;
70export interface CymbalMesh {
71 a: number;
72 pos: [number, number][]; // (x, z)
73 tris: [number, number, number][];
74 edges: { i: number; j: number; opp: number[] }[]; // opposite vertices of the 1 or 2 triangles
77export function cymbalMesh(R: number, rings: number, hole: number): CymbalMesh {
78 const a = R / (rings + 0.5);
79 const at = (i: number, j: number): [number, number] => [a * (i + 0.5 * j), a * (Math.sqrt(3) / 2) * j];
80 const index = new Map<string, number>();
81 const pos: [number, number][] = [];
82 const M = rings + 2;
83 for (let j = -M; j <= M; j++)
84 for (let i = -M - Math.abs(j); i <= M + Math.abs(j); i++) {
85 const r = Math.hypot(...at(i, j));
86 if (r <= R - 0.5 * a + 1e-9 && !(r < hole)) {
87 index.set(`${i},${j}`, pos.length);
88 pos.push(at(i, j));
89 }
90 }
91 const id = (i: number, j: number) => index.get(`${i},${j}`);
92 const tris: [number, number, number][] = [];
93 for (const k of index.keys()) {
94 const [i, j] = k.split(',').map(Number);
95 const t1 = [id(i, j), id(i + 1, j), id(i, j + 1)];
96 const t2 = [id(i + 1, j), id(i + 1, j + 1), id(i, j + 1)];
97 for (const t of [t1, t2]) if (t.every((v) => v !== undefined)) tris.push(t as [number, number, number]);
98 }
99 const em = new Map<string, { i: number; j: number; opp: number[] }>();
100 for (const [x, y, z] of tris)
101 for (const [u, v, o] of [[x, y, z], [y, z, x], [z, x, y]]) {
102 const key = u < v ? `${u},${v}` : `${v},${u}`;
103 let e = em.get(key);
104 if (!e) em.set(key, (e = { i: Math.min(u, v), j: Math.max(u, v), opp: [] }));
105 e.opp.push(o);
106 }
107 // scale so the triangulated area equals the disk's (minus the hole): the
108 // outermost row of beads otherwise leaves the plate about a/2 short
109 let meshArea = 0;
110 for (const [u, v, w] of tris) {
111 const [p, q, r] = [pos[u], pos[v], pos[w]];
112 meshArea += 0.5 * Math.abs((q[0] - p[0]) * (r[1] - p[1]) - (r[0] - p[0]) * (q[1] - p[1]));
113 }
114 const f = Math.sqrt((Math.PI * (R * R - hole * hole)) / meshArea);
115 for (const p of pos) {
116 p[0] *= f;
117 p[1] *= f;
118 }
119 return { a: a * f, pos, tris, edges: [...em.values()] };
122export interface CymbalDerived {
123 mesh: CymbalMesh;
124 n: number;
125 rest: Float64Array; // 3n
126 area: Float64Array; // n, lumped area per bead
127 normal: Float64Array; // 3n, upward unit normals
128 mass: Float64Array;
129 k: number; // interior spring stiffness
130 D: number; // bending stiffness
15fbb00Cymbal: replace felt washers with a minimal rigid hub (pinned beads around the hole), no support dampingJeremy Magland 131 support: number[]; // beads held fixed at the hub
133 sigma1: number;
134 waveSpeedP: number; // in-plane P waves (after softening)
135 bendingCoef: number; // sqrt(D / rho t): bending waves omega = coef * k^2
136 warnings: string[];
139// Triangular lattice, spacing a, Poisson ratio 1/3:
140// in-plane: spring stiffness k = sqrt(3)/2 E t (an edge with one triangle gets half)
141// bending: D = E t^3 / (12 (1 - nu^2)), through per-triangle shape operators
142export function deriveCymbal(p: CymbalParams): CymbalDerived {
143 const R = p.diameter / 2;
144 const mesh = cymbalMesh(R, Math.max(2, Math.round(p.rings)), p.holeRadius);
145 const n = mesh.pos.length;
146 const rest = new Float64Array(3 * n);
147 mesh.pos.forEach(([x, z], i) => {
148 rest[3 * i] = x;
149 rest[3 * i + 1] = profileHeight(p, Math.hypot(x, z));
150 rest[3 * i + 2] = z;
151 });
152 const area = new Float64Array(n);
153 const normal = new Float64Array(3 * n);
154 for (const [u, v, w] of mesh.tris) {
155 const e1 = [0, 1, 2].map((c) => rest[3 * v + c] - rest[3 * u + c]);
156 const e2 = [0, 1, 2].map((c) => rest[3 * w + c] - rest[3 * u + c]);
157 let nx = e1[1] * e2[2] - e1[2] * e2[1];
158 let ny = e1[2] * e2[0] - e1[0] * e2[2];
159 let nz = e1[0] * e2[1] - e1[1] * e2[0];
160 const A2 = Math.hypot(nx, ny, nz);
161 if (ny < 0) [nx, ny, nz] = [-nx, -ny, -nz];
162 for (const b of [u, v, w]) {
163 area[b] += A2 / 6;
164 normal[3 * b] += nx;
165 normal[3 * b + 1] += ny;
166 normal[3 * b + 2] += nz;
167 }
168 }
169 for (let i = 0; i < n; i++) {
170 const l = Math.hypot(normal[3 * i], normal[3 * i + 1], normal[3 * i + 2]) || 1;
171 for (let c = 0; c < 3; c++) normal[3 * i + c] /= l;
172 }
173 const rhoT = p.density * p.thickness;
174 const mass = area.map((A) => rhoT * A);
175 const Et = (p.youngs * p.thickness) / Math.max(p.soften, 1e-9);
176 const k = (Math.sqrt(3) / 2) * Et;
177 const nu = 1 / 3;
178 const D = (p.youngs * p.thickness ** 3) / (12 * (1 - nu * nu));
179 const support: number[] = [];
180 mesh.pos.forEach(([x, z], i) => {
15fbb00Cymbal: replace felt washers with a minimal rigid hub (pinned beads around the hole), no support dampingJeremy Magland 181 if (Math.hypot(x, z) <= p.hubRadius + 1e-9) support.push(i);
183 const warnings: string[] = [];
185 warnings.push(
186 support.length
187 ? 'Fewer than three beads are held at the hub: the plate can tilt freely about them and will swing or spin after an off-centre hit. Increase the hub radius.'
188 : 'No beads are held at the hub: the plate is free and will drift away after a hit. Increase the hub radius.',
189 );
191 warnings.push(
192 `In-plane stiffness is softened ${p.soften}×. A flat plate bends the same at small amplitude; curved shapes lose curvature stiffening and big hits exaggerate the nonlinear coupling.`,
193 );
194 const { sigma0, sigma1, note } = dampingFromT60(p);
195 if (note) warnings.push(note);
196 const lambdaPlus2Mu = (3 * Math.sqrt(3) * k) / 4;
197 return {
198 mesh,
199 n,
200 rest,
201 area,
202 normal,
203 mass,
204 k,
205 D,
206 support,
207 sigma0,
208 sigma1,
209 waveSpeedP: Math.sqrt(lambdaPlus2Mu / rhoT),
210 bendingCoef: Math.sqrt(D / rhoT),
211 warnings,
212 };
215// highest frequency whose bending wavelength spans six lattice spacings (where
216// the lattice is still a fair plate)
217export function bendingLimitHz(d: CymbalDerived): number {
218 const kb = (2 * Math.PI) / (6 * d.mesh.a);
219 return (d.bendingCoef * kb * kb) / (2 * Math.PI);
15fbb00Cymbal: replace felt washers with a minimal rigid hub (pinned beads around the hole), no support dampingJeremy Magland 222export function buildCymbal(d: CymbalDerived): BeadSystem {
224 const si: number[] = [];
225 const sj: number[] = [];
226 const sk: number[] = [];
227 const sL0: number[] = [];
228 const sc: number[] = [];
229 const len = (i: number, j: number) =>
230 Math.hypot(rest[3 * j] - rest[3 * i], rest[3 * j + 1] - rest[3 * i + 1], rest[3 * j + 2] - rest[3 * i + 2]);
231 for (const e of mesh.edges) {
232 si.push(e.i);
233 sj.push(e.j);
234 sk.push((d.k / 2) * e.opp.length);
235 sL0.push(len(e.i, e.j));
236 sc.push(0);
237 }
239 for (const i of d.support) pinned[i] = 1;
241 return {
242 n,
243 rest,
244 mass: d.mass,
246 anchors: [],
248 sj: Int32Array.from(sj),
249 sk: Float64Array.from(sk),
250 sL0: Float64Array.from(sL0),
251 sc: Float64Array.from(sc),
252 sigma0: d.sigma0,
253 bending,
254 };
257export function hitPoint(p: CymbalParams): [number, number] {
258 const r = (p.diameter / 2) * p.excR;
259 const th = (p.excAngle * Math.PI) / 180;
260 return [r * Math.cos(th), r * Math.sin(th)];
263function nearest(d: CymbalDerived, x: number, z: number): number {
264 let best = 0;
265 let bd = Infinity;
266 d.mesh.pos.forEach(([px, pz], i) => {
267 const q = (px - x) ** 2 + (pz - z) ** 2;
268 if (q < bd) [bd, best] = [q, i];
269 });
270 return best;
273// patch weights (sum 1), the inward normal there and the hit direction
274// (tilted from the normal toward the rim by hitAngle)
275export function hitGeometry(p: CymbalParams, d: CymbalDerived) {
276 const [cx, cz] = hitPoint(p);
277 const w = new Float64Array(d.n);
278 if (p.excWidth > 0)
279 d.mesh.pos.forEach(([x, z], i) => {
280 const r2 = ((x - cx) ** 2 + (z - cz) ** 2) / (p.excWidth * p.excWidth);
281 w[i] = r2 < 9 ? Math.exp(-0.5 * r2) : 0;
282 });
283 if (!w.some((v) => v > 0)) w[nearest(d, cx, cz)] = 1;
284 const sum = w.reduce((s, v) => s + v, 0);
285 for (let i = 0; i < d.n; i++) w[i] /= sum;
286 const up = [0, 0, 0];
287 for (let i = 0; i < d.n; i++) for (let c = 0; c < 3; c++) up[c] += w[i] * d.normal[3 * i + c];
288 const ul = Math.hypot(...up);
289 const N = up.map((v) => v / ul);
290 // outward tangent: radial direction minus its normal part
291 const rl = Math.hypot(cx, cz) || 1;
292 const er = [cx / rl, 0, cz / rl];
293 const dot = er[0] * N[0] + er[2] * N[2];
294 let t = [er[0] - dot * N[0], -dot * N[1], er[2] - dot * N[2]];
295 const tl = Math.hypot(...t) || 1;
296 t = t.map((v) => v / tl);
297 const al = (p.hitAngle * Math.PI) / 180;
298 const inward: [number, number, number] = [-N[0], -N[1], -N[2]];
299 const dir: [number, number, number] = [0, 1, 2].map((c) => Math.cos(al) * inward[c] + Math.sin(al) * t[c]) as [number, number, number];
300 return { w, inward, dir };
303export function buildCymbalExcitation(p: CymbalParams, d: CymbalDerived): Excitation {
304 const n = d.n;
305 const u0 = new Float64Array(3 * n);
306 const v0 = new Float64Array(3 * n);
307 const { w, inward, dir } = hitGeometry(p, d);
308 const wmax = Math.max(...w);
309 if (p.excitation === 'strike') {
310 for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) v0[3 * i + c] = (p.strikeVel * w[i] * dir[c]) / wmax;
311 return { u0, v0 };
312 }
313 if (p.excitation === 'force') {
314 const shape = new Float64Array(3 * n);
315 for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) shape[3 * i + c] = w[i] * dir[c];
316 return { u0, v0, force: { shape, F0: p.forcePeak, tau: p.forceMs / 1000 } };
317 }
318 const stick = p.excitation === 'stick';
319 const speed = stick ? p.stickSpeed : p.malletSpeed;
320 return {
321 u0,
322 v0,
323 mallet: {
324 weights: w,
325 normal: inward,
326 vel0: [dir[0] * speed, dir[1] * speed, dir[2] * speed],
327 mass: stick ? p.stickMass : p.malletMass,
328 K: stick ? p.stickK : p.malletK,
329 p: stick ? 1.5 : p.malletP,
330 mu: p.friction,
331 reach: Math.max(2.5 * p.excWidth, 0.75 * d.mesh.a),
332 },
333 };
336const RHO_AIR = 1.2;
337const C_AIR = 343;
339export function micDirection(p: CymbalParams): [number, number, number] {
340 const el = (p.micElevation * Math.PI) / 180;
341 const az = (p.micAzimuth * Math.PI) / 180;
342 return [Math.cos(el) * Math.cos(az), Math.sin(el), Math.cos(el) * Math.sin(az)];
345// Rayleigh integral: p(t) = rho / (2 pi) sum_i A_i a_n,i(t - R_i / c) / R_i,
346// with a_n the acceleration along the normal on the microphone's side.
347export function buildMic(p: CymbalParams, d: CymbalDerived): { weights: Float64Array; delay: Float64Array } {
348 const [mx, my, mz] = micDirection(p).map((v) => v * p.micDistance);
349 const side = p.micElevation >= 0 ? 1 : -1;
350 const weights = new Float64Array(3 * d.n);
351 const delay = new Float64Array(d.n);
352 for (let i = 0; i < d.n; i++) {
353 const R = Math.hypot(mx - d.rest[3 * i], my - d.rest[3 * i + 1], mz - d.rest[3 * i + 2]);
354 delay[i] = (R / C_AIR) * p.fs;
355 const g = (side * RHO_AIR * d.area[i]) / (2 * Math.PI * R);
356 for (let c = 0; c < 3; c++) weights[3 * i + c] = g * d.normal[3 * i + c];
357 }
358 return { weights, delay };
361// Hinges on interior edges, oriented so a positive change folds the surface
362// toward its upward normal; per-triangle plate-energy matrices, condensed over
363// free edges (the missing angle takes the value that leaves no moment there).
364function buildBending(d: CymbalDerived, nu: number): Bending {
365 const { rest, mesh, normal } = d;
366 const key = (i: number, j: number) => (i < j ? `${i},${j}` : `${j},${i}`);
367 const hingeOf = new Map<string, number>();
368 const w0: number[] = [], w1: number[] = [], e0: number[] = [], e1: number[] = [];
369 for (const e of mesh.edges)
370 if (e.opp.length === 2) {
371 hingeOf.set(key(e.i, e.j), w0.length);
372 w0.push(e.opp[0]);
373 w1.push(e.opp[1]);
374 e0.push(e.i);
375 e1.push(e.j);
376 }
377 const nh = w0.length;
378 const g = new Float64Array(13);
379 const theta0 = new Float64Array(nh);
380 const sign = new Float64Array(nh);
381 const lift = 1e-4 * mesh.a;
382 for (let e = 0; e < nh; e++) {
383 theta0[e] = hingeAngleGrad(rest, w0[e], w1[e], e0[e], e1[e], g);
384 const Y = rest.slice();
385 for (const w of [w0[e], w1[e]]) for (let c = 0; c < 3; c++) Y[3 * w + c] += lift * normal[3 * w + c];
386 sign[e] = Math.sign(hingeAngleGrad(Y, w0[e], w1[e], e0[e], e1[e], g) - theta0[e]) || 1;
387 }
388 const nt = mesh.tris.length;
389 const tri = new Int32Array(3 * nt).fill(-1);
390 const M = new Float64Array(9 * nt);
391 const P = (i: number) => [rest[3 * i], rest[3 * i + 1], rest[3 * i + 2]];
392 mesh.tris.forEach(([u, v, w], t) => {
393 const ed = [[v, w], [w, u], [u, v]];
394 const vec = ed.map(([i, j]) => P(j).map((x, c) => x - P(i)[c]));
395 const len = vec.map((q) => Math.hypot(...q));
396 const tu = vec.map((q, k) => q.map((x) => x / len[k]));
397 const [x1, x2] = [vec[2], vec[1].map((x) => -x)]; // u->v and u->w
398 const A = 0.5 * Math.hypot(x1[1] * x2[2] - x1[2] * x2[1], x1[2] * x2[0] - x1[0] * x2[2], x1[0] * x2[1] - x1[1] * x2[0]);
399 const beta = len.map((l) => l / (2 * A));
400 const full = new Float64Array(9);
401 for (let a = 0; a < 3; a++)
402 for (let b = 0; b < 3; b++) {
403 const c2 = (tu[a][0] * tu[b][0] + tu[a][1] * tu[b][1] + tu[a][2] * tu[b][2]) ** 2;
404 full[3 * a + b] = A * d.D * beta[a] * beta[b] * (nu + (1 - nu) * c2);
405 }
406 const have = ed.map(([i, j]) => hingeOf.get(key(i, j)) ?? -1);
407 have.forEach((h, a) => (tri[3 * t + a] = h));
408 const K = [0, 1, 2].filter((a) => have[a] >= 0);
409 const B = [0, 1, 2].filter((a) => have[a] < 0);
410 // condense: M' = M_KK - M_KB M_BB^-1 M_BK
411 const inv = invertSmall(B.map((a) => B.map((b) => full[3 * a + b])));
412 for (const a of K)
413 for (const b of K) {
414 let m = full[3 * a + b];
415 for (let x = 0; x < B.length; x++)
416 for (let y = 0; y < B.length; y++) m -= full[3 * a + B[x]] * inv[x][y] * full[3 * B[y] + b];
417 M[9 * t + 3 * a + b] = m;
418 }
419 });
420 return {
421 w0: Int32Array.from(w0),
422 w1: Int32Array.from(w1),
423 e0: Int32Array.from(e0),
424 e1: Int32Array.from(e1),
425 sign,
426 theta0,
427 tri,
428 M,
429 sigma1: d.sigma1,
430 };
433function invertSmall(m: number[][]): number[][] {
434 if (m.length === 0) return [];
435 if (m.length === 1) return [[1 / m[0][0]]];
436 if (m.length === 2) {
437 const det = m[0][0] * m[1][1] - m[0][1] * m[1][0];
438 return [
439 [m[1][1] / det, -m[0][1] / det],
440 [-m[1][0] / det, m[0][0] / det],
441 ];
442 }
443 // all three edges free: an isolated triangle has no bending energy
444 return m.map((r) => r.map(() => 0));
↑↓ move↵openescclose