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