1import { describe, expect, it } from 'vitest';
2import { buildCymbal, buildCymbalExcitation, buildMic, deriveCymbal, type CymbalParams } from '../src/instruments/bead-cymbal/model';
3import { CYMBAL_PRESETS } from '../src/instruments/bead-cymbal/presets';
4import { bendingEnergy, simulate } from '../src/sim/solver';
5import { spectrum } from '../src/ui/fft';
7// the flat plate (the tests compare with flat-plate theory)
8const flat = CYMBAL_PRESETS.find((p) => p.name.startsWith('Flat 8" toy, stick'))!.params;
9const base = (over: Partial<CymbalParams> = {}): CymbalParams => ({ ...flat, ...over });
11describe('cymbal lattice', () => {
12 it('bending energy matches the plate (nu = 1/3) for bending, doming and twisting', () => {
13 const p = base({ holeRadius: 0, rings: 12 });
14 const d = deriveCymbal(p);
15 const sys = buildCymbal(d);
16 const k = 0.5; // 1/m, small: k * a << 1
17 const R = p.diameter / 2;
18 const nu = 1 / 3;
19 const inner = (t: number) => d.mesh.tris[t].every((i) => Math.hypot(sys.rest[3 * i], sys.rest[3 * i + 2]) < 0.6 * R);
20 let area = 0;
21 d.mesh.tris.forEach((tr, t) => {
22 if (!inner(t)) return;
23 const [u, v, w] = tr.map((i) => [sys.rest[3 * i], sys.rest[3 * i + 2]]);
24 area += 0.5 * Math.abs((v[0] - u[0]) * (w[1] - u[1]) - (w[0] - u[0]) * (v[1] - u[1]));
25 });
26 const shapes: [string, (x: number, z: number) => number, number][] = [
27 ['cylinder', (x) => 0.5 * k * x * x, 0.5 * k * k],
28 ['sphere', (x, z) => 0.5 * k * (x * x + z * z), 0.5 * (4 - 2 * (1 - nu)) * k * k],
29 ['saddle', (x, z) => 0.5 * k * (x * x - z * z), (1 - nu) * k * k],
30 ['twist', (x, z) => k * x * z, (1 - nu) * k * k],
31 ];
32 for (const [name, f, density] of shapes) {
33 const X = new Float64Array(3 * (sys.n + sys.anchors.length));
34 X.set(sys.rest);
35 for (let i = 0; i < sys.n; i++) X[3 * i + 1] = f(sys.rest[3 * i], sys.rest[3 * i + 2]);
36 const ratio = bendingEnergy(sys.bending!, X, inner) / (d.D * density * area);
37 console.log(name, 'bending energy / plate', ratio.toFixed(4));
38 expect(ratio).toBeCloseTo(1, 2);
39 }
40 });
42 it('a free flat plate rings near the (2,0) mode of a free circular plate', () => {
43 // no felts, no hole, no damping; initial shape r^2 cos(2 phi), small
44 const p = base({ holeRadius: 0, rings: 9, t60Low: 0, t60High: 0, duration: 0.3, frameRate: 50 });
45 const d = deriveCymbal(p);
46 const sys = buildCymbal(d);
47 const exc = buildCymbalExcitation({ ...p, excitation: 'strike', strikeVel: 0 }, d);
48 for (let i = 0; i < sys.n; i++) {
49 const x = sys.rest[3 * i], z = sys.rest[3 * i + 2];
50 exc.u0[3 * i + 1] = 1e-7 * (x * x - z * z) / 0.01;
51 }
52 const pick = new Float64Array(3 * sys.n);
53 let edge = 0;
54 for (let i = 0; i < sys.n; i++) if (sys.rest[3 * i] > sys.rest[3 * edge]) edge = i;
55 pick[3 * edge + 1] = 1;
56 const r = simulate(sys, exc, { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: pick });
57 const s = spectrum(r.audio, r.fs);
58 let best = 1;
59 for (let k = 1; k < s.freq.length && s.freq[k] < 300; k++) if (s.db[k] > s.db[best]) best = k;
60 const R = p.diameter / 2;
61 const ideal = (5.253 / (2 * Math.PI * R * R)) * d.bendingCoef; // Leissa, nu ~ 1/3
62 console.log('free (2,0) measured', s.freq[best].toFixed(1), 'Hz, ideal', ideal.toFixed(1), 'Hz');
63 expect(Math.abs(s.freq[best] / ideal - 1)).toBeLessThan(0.05);
64 }, 60000);
66 it('energy balances with a glancing stick hit, friction, damping and rigid motion removed', () => {
67 const p = base({ rings: 7, hitAngle: 60, stickSpeed: 6, bowHeight: 0.01, bellHeight: 0.01, duration: 0.03, frameRate: 2000 });
68 const d = deriveCymbal(p);
69 const mic = buildMic(p, d);
70 const r = simulate(buildCymbal(d), buildCymbalExcitation(p, d), { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: mic.weights, pickupDelay: mic.delay });
71 const emax = Math.max(...r.energy);
72 const bal = Math.max(...r.balance.map(Math.abs));
73 const cf = r.contactForce!;
74 console.log('contact samples', cf.reduce((c, v) => c + (v > 0 ? 1 : 0), 0), 'peak N', Math.max(...cf).toFixed(0), 'balance', (bal / emax).toExponential(2));
75 expect(bal / emax).toBeLessThan(1e-2);
76 expect(cf[cf.length - 1]).toBe(0);
77 expect(r.audio.some((v) => v !== 0)).toBe(true);
78 });
80 it('removing rigid motion leaves the plate with no net momentum and keeps the hit point vibrating', () => {
81 const p = base({ rings: 7, duration: 0.05, frameRate: 2000 });
82 const d = deriveCymbal(p);
83 const sys = buildCymbal(d);
84 const mic = buildMic(p, d);
85 const r = simulate(sys, buildCymbalExcitation(p, d), { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: mic.weights, pickupDelay: mic.delay });
86 // mean vertical displacement (rigid translation) stays near zero compared with the vibration
87 const n = sys.n;
88 let meanMax = 0, ampMax = 0;
89 for (let f = 0; f < r.frameTimes.length; f++) {
90 let m = 0, mt = 0;
91 for (let i = 0; i < n; i++) {
92 m += sys.mass[i] * r.frames[f * 3 * n + 3 * i + 1];
93 mt += sys.mass[i];
94 ampMax = Math.max(ampMax, Math.abs(r.frames[f * 3 * n + 3 * i + 1]));
95 }
96 meanMax = Math.max(meanMax, Math.abs(m / mt));
97 }
98 console.log('mean y', (meanMax * 1000).toFixed(4), 'mm, max |y|', (ampMax * 1000).toFixed(3), 'mm');
99 expect(meanMax / ampMax).toBeLessThan(1e-2);
100 expect(Math.max(...r.balance.map(Math.abs)) / Math.max(...r.energy)).toBeLessThan(1e-2);
101 }, 60000);
102});