1import type { Anchor, BeadSystem, Excitation } from '../../sim/types';
2import { dampingFromT60 } from '../bead-string/model';
4// A round drumhead: beads on a triangular (hex) lattice in the x-z plane,
5// cut to a circle, held by a rim of fixed lattice points just outside it.
6// Transverse vibration is along y, as on the string.
8export type DrumExcitation = 'mallet' | 'strike' | 'force' | 'press';
10export interface DrumParams {
11 // geometry
12 diameter: number; // m, rim to rim
13 rings: number; // beads from centre to edge (lattice spacing a = R / (rings + 1))
14 // membrane properties (continuum values; the lattice is chosen to match them)
15 density: number; // kg/m^2, surface density
16 tension: number; // N/m, isotropic membrane tension
17 inPlane: number; // N/m, in-plane stiffness E t
18 // damping
19 t60Low: number;
20 fLow: number;
21 t60High: number;
22 fHigh: number;
23 // excitation
24 excitation: DrumExcitation;
25 excR: number; // fraction of radius, 0 = centre
26 excAngle: number; // degrees, in the x-z plane from +x
27 excWidth: number; // m, Gaussian sigma of the contact patch (0 = one bead)
28 malletMass: number; // kg
29 malletSpeed: number; // m/s
30 feltK: number; // N/m^p
31 feltP: number;
32 strikeVel: number; // m/s
33 forcePeak: number; // N
34 forceMs: number; // ms
35 pressAmp: number; // m
36 // pickup
37 pickupMode: 'bead' | 'sum';
38 pickupR: number;
39 pickupAngle: number;
40 pickupAxis: 'y' | 'x' | 'z';
41 // simulation
42 duration: number;
43 fs: number;
44 cfl: number;
45 frameRate: number;
46}
48export interface Lattice {
49 a: number; // spacing (m)
50 pos: [number, number][]; // bead positions (x, z)
51 rim: [number, number][]; // anchor positions (x, z)
52 edges: [number, number][]; // node pairs; node >= pos.length is rim[node - n]
53}
55const key = (i: number, j: number) => `${i},${j}`;
57export function hexLattice(R: number, rings: number): Lattice {
58 const a = R / (rings + 1);
59 const at = (i: number, j: number): [number, number] => [a * (i + 0.5 * j), a * (Math.sqrt(3) / 2) * j];
60 const inside = (i: number, j: number) => Math.hypot(...at(i, j)) <= R - 0.5 * a + 1e-9;
61 const index = new Map<string, number>();
62 const pos: [number, number][] = [];
63 const M = rings + 3;
64 for (let j = -M; j <= M; j++)
65 for (let i = -M - Math.abs(j); i <= M + Math.abs(j); i++)
66 if (inside(i, j)) {
67 index.set(key(i, j), pos.length);
68 pos.push(at(i, j));
69 }
70 const n = pos.length;
71 const rimIndex = new Map<string, number>();
72 const rim: [number, number][] = [];
73 const edges: [number, number][] = [];
74 const nb = [[1, 0], [0, 1], [-1, 1]];
75 const all = [[1, 0], [0, 1], [-1, 1], [-1, 0], [0, -1], [1, -1]];
76 for (const [k, b] of index) {
77 const [i, j] = k.split(',').map(Number);
78 // bead-bead springs: each pair once
79 for (const [di, dj] of nb) {
80 const o = index.get(key(i + di, j + dj));
81 if (o !== undefined) edges.push([b, o]);
82 }
83 // bead-rim springs: every missing neighbour is a rim point
84 for (const [di, dj] of all) {
85 const kk = key(i + di, j + dj);
86 if (index.has(kk)) continue;
87 let r = rimIndex.get(kk);
88 if (r === undefined) {
89 r = rim.length;
90 rimIndex.set(kk, r);
91 rim.push(at(i + di, j + dj));
92 }
93 edges.push([b, n + r]);
94 }
95 }
96 return { a, pos, rim, edges };
97}
99export interface DrumDerived {
100 lattice: Lattice;
101 n: number;
102 beadMass: number;
103 springK: number;
104 springTension: number;
105 restLength: number;
106 prestretch: number;
107 waveSpeedT: number;
108 waveSpeedP: number; // in-plane longitudinal (P) wave speed
109 idealF01: number; // ideal clamped membrane fundamental
110 sigma0: number;
111 sigma1: number;
112 dashpot: number;
113 warnings: string[];
114}
116// Triangular spring lattice, spacing a, spring stiffness k, spring tension f:
117// membrane tension tau = sqrt(3) f / a, 2D Young's modulus E t = 2 k / sqrt(3)
118// (Poisson ratio 1/3), area per bead sqrt(3)/2 a^2.
119export function deriveDrum(p: DrumParams): DrumDerived {
120 const R = p.diameter / 2;
121 const lattice = hexLattice(R, Math.max(1, Math.round(p.rings)));
122 const a = lattice.a;
123 const beadMass = p.density * (Math.sqrt(3) / 2) * a * a;
124 const springK = (Math.sqrt(3) / 2) * p.inPlane;
125 const springTension = (p.tension * a) / Math.sqrt(3);
126 const restLength = a - springTension / springK;
127 const warnings: string[] = [];
128 if (!(restLength > 0)) warnings.push('Tension too high for this in-plane stiffness: the springs would need a negative rest length.');
129 const { sigma0, sigma1, note } = dampingFromT60(p);
130 if (note) warnings.push(note);
131 const waveSpeedT = Math.sqrt(Math.max(p.tension, 0) / p.density);
132 const lambdaPlus2Mu = (3 * Math.sqrt(3) * springK) / 4;
133 return {
134 lattice,
135 n: lattice.pos.length,
136 beadMass,
137 springK,
138 springTension,
139 restLength,
140 prestretch: a / restLength - 1,
141 waveSpeedT,
142 waveSpeedP: Math.sqrt(lambdaPlus2Mu / p.density),
143 idealF01: (2.404826 * waveSpeedT) / (2 * Math.PI * R),
144 sigma0,
145 sigma1,
146 // each spring's transverse stiffness is f / a, so this damps transverse modes by sigma1 omega^2
147 dashpot: 2 * sigma1 * (springTension > 0 ? springTension / a : springK),
148 warnings,
149 };
150}
152export function buildDrum(d: DrumDerived): BeadSystem {
153 const { pos, rim, edges } = d.lattice;
154 const n = pos.length;
155 const rest = new Float64Array(3 * n);
156 pos.forEach(([x, z], i) => {
157 rest[3 * i] = x;
158 rest[3 * i + 2] = z;
159 });
160 const anchors: Anchor[] = rim.map(([x, z]) => ({ pos: [x, 0, z], slideBead: -1 }));
161 const ns = edges.length;
162 return {
163 n,
164 rest,
165 mass: new Float64Array(n).fill(d.beadMass),
166 pinned: new Uint8Array(n),
167 anchors,
168 si: Int32Array.from(edges.map((e) => e[0])),
169 sj: Int32Array.from(edges.map((e) => e[1])),
170 sk: new Float64Array(ns).fill(d.springK),
171 sL0: new Float64Array(ns).fill(d.restLength),
172 sc: new Float64Array(ns).fill(d.dashpot),
173 sigma0: d.sigma0,
174 };
175}
177const polar = (R: number, frac: number, deg: number): [number, number] => {
178 const th = (deg * Math.PI) / 180;
179 return [R * frac * Math.cos(th), R * frac * Math.sin(th)];
180};
182export function nearestBead(d: DrumDerived, x: number, z: number): number {
183 let best = 0;
184 let bd = Infinity;
185 d.lattice.pos.forEach(([px, pz], i) => {
186 const q = (px - x) ** 2 + (pz - z) ** 2;
187 if (q < bd) {
188 bd = q;
189 best = i;
190 }
191 });
192 return best;
193}
195// Gaussian patch around the excitation point; width 0 = the nearest bead
196function patch(p: DrumParams, d: DrumDerived): Float64Array {
197 const [cx, cz] = polar(p.diameter / 2, p.excR, p.excAngle);
198 const n = d.n;
199 const w = new Float64Array(n);
200 if (!(p.excWidth > 0)) {
201 w[nearestBead(d, cx, cz)] = 1;
202 return w;
203 }
204 d.lattice.pos.forEach(([x, z], i) => {
205 const r2 = ((x - cx) ** 2 + (z - cz) ** 2) / (p.excWidth * p.excWidth);
206 w[i] = r2 < 9 ? Math.exp(-0.5 * r2) : 0;
207 });
208 if (!w.some((v) => v > 0)) w[nearestBead(d, cx, cz)] = 1;
209 return w;
210}
212export function excitationPoint(p: DrumParams): [number, number] {
213 return polar(p.diameter / 2, p.excR, p.excAngle);
214}
216export function buildDrumExcitation(p: DrumParams, d: DrumDerived): Excitation {
217 const n = d.n;
218 const u0 = new Float64Array(3 * n);
219 const v0 = new Float64Array(3 * n);
220 const w = patch(p, d);
221 let max = 0;
222 let sum = 0;
223 for (let i = 0; i < n; i++) {
224 max = Math.max(max, w[i]);
225 sum += w[i];
226 }
227 // everything acts downward (-y), onto the head from above
228 if (p.excitation === 'press') {
229 for (let i = 0; i < n; i++) u0[3 * i + 1] = (-p.pressAmp * w[i]) / max;
230 return { u0, v0 };
231 }
232 if (p.excitation === 'strike') {
233 for (let i = 0; i < n; i++) v0[3 * i + 1] = (-p.strikeVel * w[i]) / max;
234 return { u0, v0 };
235 }
236 const weights = w.map((v) => v / sum);
237 if (p.excitation === 'force') {
238 const shape = new Float64Array(3 * n);
239 for (let i = 0; i < n; i++) shape[3 * i + 1] = -weights[i];
240 return { u0, v0, force: { shape, F0: p.forcePeak, tau: p.forceMs / 1000 } };
241 }
242 return {
243 u0,
244 v0,
245 mallet: { weights, dir: [0, -1, 0], mass: p.malletMass, K: p.feltK, p: p.feltP, v0: p.malletSpeed },
246 };
247}
249export function buildDrumPickup(p: DrumParams, d: DrumDerived): Float64Array {
250 const w = new Float64Array(3 * d.n);
251 const axis = p.pickupAxis === 'x' ? 0 : p.pickupAxis === 'y' ? 1 : 2;
252 if (p.pickupMode === 'sum') for (let i = 0; i < d.n; i++) w[3 * i + axis] = 1;
253 else w[3 * nearestBead(d, ...polar(p.diameter / 2, p.pickupR, p.pickupAngle)) + axis] = 1;
254 return w;
255}