4ddcc1bCymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web appJeremy Magland 1// Mirror of ../matlab/cymbal_profile.m: the cymbal meridian z(r), used by the
2// 3D scene and the mode gallery. The physics reads the same geometry from the
3// MATLAB side; this is only for drawing.
4//
5// Bow (r >= bell radius): a parabola in the rim-referred coordinate, flat at
6// the rim. Bell (r < bell radius): the quartic that matches the bow's height,
7// slope and curvature at the join and comes to a horizontal apex at r = 0.
9export interface Meridian {
10 rHole: number;
11 rBell: number;
12 rRim: number;
13 z(r: number): number;
14 thickness(r: number): number;
15 apexHeight: number;
16}
18export function cymbalMeridian(
19 radius: number,
20 holeRadius: number,
21 bellRadius: number,
22 bellHeight: number,
23 bowHeight: number,
24 thicknessBell: number,
25 thicknessRim: number,
26): Meridian {
27 const a = radius;
28 const rh = Math.max(1e-4, Math.min(holeRadius, 0.5 * a));
29 const rb = Math.min(Math.max(bellRadius, 2 * rh), 0.8 * a);
30 const den = a - rb;
31 const mb = (-2 * bowHeight) / den;
32 const kb = (2 * bowHeight) / (den * den);
33 // solve the 3x3 for the quartic coefficients (B r^2 + C r^3 + E r^4)
34 const M = [
35 [rb * rb, rb ** 3, rb ** 4],
36 [2 * rb, 3 * rb * rb, 4 * rb ** 3],
37 [2, 6 * rb, 12 * rb * rb],
38 ];
39 const rhs = [-bellHeight, mb, kb];
40 const [B, C, E] = solve3(M, rhs);
41 const A0 = bowHeight + bellHeight;
42 return {
43 rHole: rh,
44 rBell: rb,
45 rRim: a,
46 apexHeight: A0,
47 z(r: number): number {
48 if (r < rb) return A0 + B * r * r + C * r ** 3 + E * r ** 4;
49 const u = (a - r) / den;
50 return bowHeight * u * u;
51 },
52 thickness(r: number): number {
53 const taper = Math.min(1, Math.max(0, (a - r) / den));
54 return thicknessRim + (thicknessBell - thicknessRim) * taper;
55 },
56 };
57}
59function solve3(M: number[][], b: number[]): [number, number, number] {
60 const A = M.map((row, i) => [...row, b[i]]);
61 for (let c = 0; c < 3; c++) {
62 let piv = c;
63 for (let r = c + 1; r < 3; r++) if (Math.abs(A[r][c]) > Math.abs(A[piv][c])) piv = r;
64 [A[c], A[piv]] = [A[piv], A[c]];
65 const d = A[c][c] || 1e-30;
66 for (let k = c; k < 4; k++) A[c][k] /= d;
67 for (let r = 0; r < 3; r++) {
68 if (r === c) continue;
69 const f = A[r][c];
70 for (let k = c; k < 4; k++) A[r][k] -= f * A[c][k];
71 }
72 }
73 return [A[0][3], A[1][3], A[2][3]];
74}
a136eacMake the contact point a first-class controlJeremy Magland 75
76export interface MeridianSamples {
77 r: Float64Array; // radius of each sample, m
78 z: Float64Array; // height above the rim plane, m
79 sFrac: Float64Array; // arclength fraction from the hole, the model's strike_r
80 total: number; // meridian arclength, m
81}
83// Sample the meridian and tabulate arclength, so that the model's strike_r
84// (an arclength fraction) can be converted to and from a radius. Both the 3D
85// scene and the strike picker need this; the physics reads the same geometry
86// from cymbal_profile.m.
87export function meridianSamples(mer: Meridian, n = 256): MeridianSamples {
88 const r = new Float64Array(n + 1);
89 const z = new Float64Array(n + 1);
90 const s = new Float64Array(n + 1);
91 for (let i = 0; i <= n; i++) {
92 r[i] = mer.rHole + (i / n) * (mer.rRim - mer.rHole);
93 z[i] = mer.z(r[i]);
94 }
95 for (let i = 1; i <= n; i++) s[i] = s[i - 1] + Math.hypot(r[i] - r[i - 1], z[i] - z[i - 1]);
96 const total = s[n];
97 const sFrac = new Float64Array(n + 1);
98 for (let i = 0; i <= n; i++) sFrac[i] = s[i] / total;
99 return { r, z, sFrac, total };
100}
102// strike_r -> the point on the meridian, matching cymbal_point.m
103export function pointAtFrac(m: MeridianSamples, frac: number) {
104 const f = Math.min(0.999999, Math.max(0, frac));
105 let i = 0;
106 while (i < m.sFrac.length - 2 && m.sFrac[i + 1] < f) i++;
107 const t = (f - m.sFrac[i]) / Math.max(1e-12, m.sFrac[i + 1] - m.sFrac[i]);
108 const dr = m.r[i + 1] - m.r[i];
109 const dz = m.z[i + 1] - m.z[i];
110 const dl = Math.hypot(dr, dz) || 1;
111 return { r: m.r[i] + t * dr, z: m.z[i] + t * dz, tr: dr / dl, tz: dz / dl };
112}
114// radius -> strike_r, the inverse the picker needs when the pointer moves
115export function fracAtRadius(m: MeridianSamples, radius: number): number {
116 const n = m.r.length - 1;
117 const rr = Math.min(m.r[n], Math.max(m.r[0], radius));
118 let i = 0;
119 while (i < n - 1 && m.r[i + 1] < rr) i++;
120 const t = (rr - m.r[i]) / Math.max(1e-12, m.r[i + 1] - m.r[i]);
121 return m.sFrac[i] + t * (m.sFrac[i + 1] - m.sFrac[i]);
122}