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