1// Simulation parameters (strike, listener, numerics), the grid with its mask,
2// the plate as evaluated on the grid (material, thickness field, stress-free
3// shape, mount), and the derived quantities the kernel needs. All values are
4// SI. The cymbal design itself is produced by the script in design.ts.
6export interface SimParams {
7 // Strike
8 strikeFracX: number; // strike point as a fraction of the bounding-box width from the left edge
9 strikeFracY: number; // strike point as a fraction of the bounding-box height from the bottom edge
10 strikeVelocity: number; // stick tip speed at impact [m/s]
11 stickMass: number; // effective mass of the stick [kg]
12 contactStiffness: number; // Hertz contact stiffness K [N/m^alpha]
13 contactExponent: number; // Hertz exponent alpha (1.5 for a rounded tip)
14 contactRadius: number; // radius of the contact patch [m]
15 // Listener
16 listenerDist: number; // [m]
17 listenerAzimuth: number; // [deg], rotation about the plate's vertical axis (0 = on the normal)
18 listenerElevation: number; // [deg]
19 // Numerics and output
20 spacing: number; // requested grid spacing [m]
21 inPlaneMass: number; // factor by which the in-plane inertia is scaled; 1 is physical, larger lengthens the time step
22 rigidProjection: boolean; // remove the six rigid-body modes every RB_INTERVAL steps (free plates only)
23 fs: number; // audio sample rate [Hz]
24 fps: number; // movie frame rate, also the batch rate [Hz]
25 duration: number; // recorded duration after the impact [s]
26}
28export const defaultParams: SimParams = {
29 strikeFracX: 0.7,
30 strikeFracY: 0.5,
31 strikeVelocity: 3,
32 stickMass: 0.03,
33 contactStiffness: 5e8,
34 contactExponent: 1.5,
35 contactRadius: 5e-3,
36 listenerDist: 1.0,
37 listenerAzimuth: 0,
38 listenerElevation: 0,
39 spacing: 0.008,
40 inPlaneMass: 1,
41 rigidProjection: false,
42 fs: 48000,
43 fps: 120,
44 duration: 4,
45};
47export const SPEED_OF_SOUND = 343; // [m/s]
48export const RHO_AIR = 1.2; // [kg/m^3]
49export const CFL = 0.9; // fraction of the stability limit used for the time step
50// The ring buffer of the retarded-time listener sum holds one batch of samples
51// plus this margin, which bounds the spread of the delays across the plate.
52export const HIST_MARGIN = 128; // [samples]
53export const RB_INTERVAL = 32; // steps between rigid-body projections
55// ---------- Grid ----------
57export interface Grid {
58 nx: number;
59 ny: number;
60 hx: number;
61 hy: number;
62 n: number;
63 /** 1 for nodes that belong to the plate, 0 otherwise; absent for a full rectangle. */
64 mask?: Uint8Array;
65}
67/** True when the point (x, y), relative to the center of the bounding box, lies in the plate. */
68export type InsideFn = (x: number, y: number) => boolean;
70/**
71 * The grid covering the bounding box Lx by Ly with about the requested
72 * spacing, restricted to the outline `inside` when one is given.
73 */
74export function makeGrid(Lx: number, Ly: number, spacing: number, inside?: InsideFn): Grid {
75 const nx = Math.max(6, Math.round(Lx / spacing) + 1);
76 const ny = Math.max(6, Math.round(Ly / spacing) + 1);
77 const g: Grid = { nx, ny, hx: Lx / (nx - 1), hy: Ly / (ny - 1), n: nx * ny };
78 return inside ? withMask(g, inside) : g;
79}
81/**
82 * Restricts a grid to an outline. The edge is a staircase of whole cells: a
83 * cell belongs to the plate when its center lies inside the outline, and a
84 * node when it is a corner of such a cell. Deciding by cell centers rather
85 * than by nodes keeps the covered area unbiased; a node-based mask loses
86 * about half a cell all around the edge, which raises every frequency by a
87 * few percent. A node with no cell of its own would carry no mass, and the
88 * trailing loop removes any such node. A hole in the outline (a cymbal's
89 * center hole) is handled by the same rule.
90 */
91export function withMask(g: Grid, inside: InsideFn): Grid {
92 return withCellMask(g, (ci, cj) => inside(nodeX(g, ci) + g.hx / 2, nodeY(g, cj) + g.hy / 2));
93}
95/** The same restriction with the membership of each cell (ci, cj) given directly. */
96export function withCellMask(g: Grid, cellInside: (ci: number, cj: number) => boolean): Grid {
97 const mask = new Uint8Array(g.n);
98 for (let cj = 0; cj < g.ny - 1; cj++) {
99 for (let ci = 0; ci < g.nx - 1; ci++) {
100 if (cellInside(ci, cj)) {
101 mask[cj * g.nx + ci] = 1; mask[cj * g.nx + ci + 1] = 1;
102 mask[(cj + 1) * g.nx + ci] = 1; mask[(cj + 1) * g.nx + ci + 1] = 1;
103 }
104 }
105 }
106 const gm: Grid = { ...g, mask };
107 for (let changed = true; changed; ) {
108 changed = false;
109 for (let j = 0; j < g.ny; j++) {
110 for (let i = 0; i < g.nx; i++) {
111 if (mask[j * g.nx + i] && nodeWeight(gm, i, j) === 0) {
112 mask[j * g.nx + i] = 0;
113 changed = true;
114 }
115 }
116 }
117 }
118 return gm;
119}
121/** The ellipse inscribed in the bounding box (kept for the tests of the discretization). */
122export function withEllipseMask(g: Grid): Grid {
123 const a = ((g.nx - 1) * g.hx) / 2, b = ((g.ny - 1) * g.hy) / 2;
124 return withMask(g, (x, y) => (x / a) ** 2 + (y / b) ** 2 <= 1);
125}
127/** Node coordinates relative to the center of the bounding box. */
128export function nodeX(g: Grid, i: number): number {
129 return (i - (g.nx - 1) / 2) * g.hx;
130}
131export function nodeY(g: Grid, j: number): number {
132 return (j - (g.ny - 1) / 2) * g.hy;
133}
135/** True when node (i, j) exists and belongs to the plate. */
136export function insideNode(g: Grid, i: number, j: number): boolean {
137 if (i < 0 || i >= g.nx || j < 0 || j >= g.ny) return false;
138 return g.mask ? g.mask[j * g.nx + i] !== 0 : true;
139}
141/** True when cell (ci, cj), with nodes (ci..ci+1, cj..cj+1), lies within the plate. */
142export function insideCell(g: Grid, ci: number, cj: number): boolean {
143 return insideNode(g, ci, cj) && insideNode(g, ci + 1, cj) && insideNode(g, ci + 1, cj + 1) && insideNode(g, ci, cj + 1);
144}
146/**
147 * Quadrature and lumped-mass weight of a node: the fraction of its four
148 * surrounding cells that belong to the plate. On a full rectangle this is the
149 * trapezoidal rule (1 inside, 1/2 on an edge, 1/4 at a corner). Zero outside.
150 */
151export function nodeWeight(g: Grid, i: number, j: number): number {
152 if (!insideNode(g, i, j)) return 0;
153 let c = 0;
154 if (insideCell(g, i, j)) c++;
155 if (insideCell(g, i - 1, j)) c++;
156 if (insideCell(g, i, j - 1)) c++;
157 if (insideCell(g, i - 1, j - 1)) c++;
158 return c / 4;
159}
161/** Per-node weights as a flat array. */
162export function nodeWeights(g: Grid): Float32Array<ArrayBuffer> {
163 const w = new Float32Array(g.n);
164 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) w[j * g.nx + i] = nodeWeight(g, i, j);
165 return w;
166}
168/** Per-cell membership (1 or 0) as a flat array of (nx - 1) * (ny - 1). */
169export function cellMask(g: Grid): Float32Array<ArrayBuffer> {
170 const c = new Float32Array((g.nx - 1) * (g.ny - 1));
171 for (let cj = 0; cj < g.ny - 1; cj++) for (let ci = 0; ci < g.nx - 1; ci++) c[cj * (g.nx - 1) + ci] = insideCell(g, ci, cj) ? 1 : 0;
172 return c;
173}
175/** Number of nodes in the plate. */
176export function countInside(g: Grid): number {
177 if (!g.mask) return g.n;
178 let c = 0;
179 for (let k = 0; k < g.n; k++) c += g.mask[k];
180 return c;
181}
183// ---------- Material and plate ----------
185export interface Material {
186 E: number; // Young's modulus [Pa]
187 nu: number; // Poisson ratio
188 rho: number; // density [kg/m^3]
189}
191/** Losses of the plate (see THEORY.md). */
192export interface Loss {
193 sigma0: number; // frequency-independent loss of all motion [1/s]
194 sigma1: number; // loss proportional to the squared wavenumber, transverse motion [m^2/s]
195 sigmaIn: number; // additional frequency-independent loss of the in-plane motion [1/s]
196}
198/**
199 * The stand: the plate nodes within `radius` of (x, y) are tied to the ground
200 * by a spring of total stiffness `stiffness` and a dashpot of total rate
201 * `damping`, both spread over the nodes in proportion to their weight and
202 * acting on all three displacement components.
203 */
204export interface Mount {
205 x: number; // [m] from the center of the bounding box
206 y: number;
207 radius: number; // [m]
208 stiffness: number; // [N/m]
209 damping: number; // [N s/m]
210}
212/** A design evaluated on a grid: everything the simulation needs to know about the plate. */
213export interface Plate {
214 name: string;
215 Lx: number; // bounding box [m]
216 Ly: number;
217 grid: Grid;
218 mat: Material;
219 /** Thickness at every node [m]; zero outside the plate. */
220 h: Float64Array;
221 /** Stress-free shape at every node [m], mean plane removed; zero outside the plate. */
222 w0: Float64Array;
223 mount: Mount | null;
224 loss: Loss;
225}
227/** Bending stiffness [N m] of a plate of thickness h. */
228export function bendingStiffness(m: Material, h: number): number {
229 return (m.E * h ** 3) / (12 * (1 - m.nu * m.nu));
230}
231/** Membrane stiffness E h / (1 - nu^2) [N/m]. */
232export function membraneStiffness(m: Material, h: number): number {
233 return (m.E * h) / (1 - m.nu * m.nu);
234}
235/** Membrane shear stiffness E h / (2 (1 + nu)) [N/m]. */
236export function shearStiffness(m: Material, h: number): number {
237 return (m.E * h) / (2 * (1 + m.nu));
238}
240/** Mass of node (i, j) [kg]: density, local thickness, cell area, and the node's weight. */
241export function nodeMass(pl: Plate, i: number, j: number): number {
242 const g = pl.grid;
243 return pl.mat.rho * pl.h[j * g.nx + i] * g.hx * g.hy * nodeWeight(g, i, j);
244}
246/** Thickness of a cell: the mean over its four corners (zero for a cell outside the plate). */
247export function cellThickness(pl: Plate, ci: number, cj: number): number {
248 const g = pl.grid;
249 if (!insideCell(g, ci, cj)) return 0;
250 const k = cj * g.nx + ci;
251 return 0.25 * (pl.h[k] + pl.h[k + 1] + pl.h[k + g.nx] + pl.h[k + g.nx + 1]);
252}
254/** Smallest, mean, and largest thickness over the plate nodes. */
255export function thicknessRange(pl: Plate): { min: number; mean: number; max: number } {
256 const g = pl.grid;
257 let lo = Infinity, hi = 0, sum = 0, sw = 0;
258 for (let j = 0; j < g.ny; j++) {
259 for (let i = 0; i < g.nx; i++) {
260 const wgt = nodeWeight(g, i, j);
261 if (wgt === 0) continue;
262 const h = pl.h[j * g.nx + i];
263 lo = Math.min(lo, h); hi = Math.max(hi, h); sum += wgt * h; sw += wgt;
264 }
265 }
266 return { min: isFinite(lo) ? lo : 0, mean: sw > 0 ? sum / sw : 0, max: hi };
267}
269/** Total mass [kg] and the moments of inertia about the in-plane axes through the center of the bounding box. */
270export function massProperties(pl: Plate): { mass: number; Ix: number; Iy: number } {
271 const g = pl.grid;
272 let mass = 0, Ix = 0, Iy = 0;
273 for (let j = 0; j < g.ny; j++) {
274 for (let i = 0; i < g.nx; i++) {
275 const mk = nodeMass(pl, i, j);
276 mass += mk;
277 Ix += mk * nodeY(g, j) ** 2;
278 Iy += mk * nodeX(g, i) ** 2;
279 }
280 }
281 return { mass, Ix, Iy };
282}
284/**
285 * 60 dB decay time [s] of a bending mode of angular frequency omega under the
286 * loss model sigma0 + sigma1 k^2, using the flat-plate dispersion k^2 = omega / kappa
287 * for a plate of thickness h.
288 */
289export function t60(m: Material, h: number, loss: Loss, omega: number): number {
290 const kappa = Math.sqrt(bendingStiffness(m, h) / (m.rho * h));
291 const sigma = loss.sigma0 + loss.sigma1 * (omega / kappa);
292 return (3 * Math.LN10) / sigma;
293}
295/**
296 * Frequency [Hz] at which the three-point second difference lowers a bending
297 * frequency by the relative amount `err` on a grid of spacing hmax, for a
298 * plate of thickness h. A plane bending wave along a grid axis with
299 * wavenumber k has the discrete frequency ratio (2 - 2 cos kh) / (kh)^2.
300 */
301export function dispersionFrequency(m: Material, h: number, g: Grid, err: number): number {
302 const hg = Math.max(g.hx, g.hy);
303 let lo = 0, hi = Math.PI;
304 for (let it = 0; it < 60; it++) {
305 const x = (lo + hi) / 2;
306 if ((2 - 2 * Math.cos(x)) / (x * x) > 1 - err) lo = x; else hi = x;
307 }
308 const k = lo / hg;
309 const kappa = Math.sqrt(bendingStiffness(m, h) / (m.rho * h));
310 return (kappa * k * k) / (2 * Math.PI);
311}
313// ---------- Geometry of the excitation and the mount ----------
315/** Strike point in plate coordinates. */
316export function strikePoint(p: SimParams, pl: Plate): [number, number] {
317 return [(p.strikeFracX - 0.5) * pl.Lx, (p.strikeFracY - 0.5) * pl.Ly];
318}
320/** The plate node nearest to a point, by index; -1 when the plate has no nodes. */
321export function nearestNode(g: Grid, x: number, y: number): number {
322 let nearest = -1, best = Infinity;
323 for (let j = 0; j < g.ny; j++) {
324 for (let i = 0; i < g.nx; i++) {
325 if (nodeWeight(g, i, j) === 0) continue;
326 const d2 = (nodeX(g, i) - x) ** 2 + (nodeY(g, j) - y) ** 2;
327 if (d2 < best) { best = d2; nearest = j * g.nx + i; }
328 }
329 }
330 return nearest;
331}
333export interface Patch {
334 nodes: number[];
335 /** The plate node nearest the center; always a member. */
336 nearest: number;
337 /** Total mass of the member nodes [kg]. */
338 mass: number;
339 /** Total node weight of the members. */
340 weight: number;
341 /** Mass-weighted center [m]. */
342 center: [number, number];
343 /** Index box enclosing the members (i0, i1, j0, j1, inclusive). */
344 box: [number, number, number, number];
345}
347/**
348 * The plate nodes within `radius` of (x, y), plus the nearest one. Used for
349 * the drumstick contact patch and for the mount. The contact force is spread
350 * over the patch in proportion to the node masses, and the contact coordinate
351 * is their mass-weighted mean displacement, so the two are conjugate and the
352 * contact has a potential energy. The mass-weighted center of the patch is
353 * where the inertia-relief reaction to the force is taken.
354 */
355export function patchAround(pl: Plate, x: number, y: number, radius: number): Patch {
356 const g = pl.grid;
357 const nearest = nearestNode(g, x, y);
358 const nodes: number[] = [];
359 let mass = 0, weight = 0, cx = 0, cy = 0;
360 const box: [number, number, number, number] = [g.nx, -1, g.ny, -1];
361 for (let j = 0; j < g.ny; j++) {
362 for (let i = 0; i < g.nx; i++) {
363 const k = j * g.nx + i;
364 const wgt = nodeWeight(g, i, j);
365 if (wgt === 0) continue;
366 const d2 = (nodeX(g, i) - x) ** 2 + (nodeY(g, j) - y) ** 2;
367 if (d2 <= radius * radius || k === nearest) {
368 nodes.push(k);
369 const mk = nodeMass(pl, i, j);
370 mass += mk; weight += wgt; cx += mk * nodeX(g, i); cy += mk * nodeY(g, j);
371 box[0] = Math.min(box[0], i); box[1] = Math.max(box[1], i);
372 box[2] = Math.min(box[2], j); box[3] = Math.max(box[3], j);
373 }
374 }
375 }
376 return { nodes, nearest, mass, weight, center: mass > 0 ? [cx / mass, cy / mass] : [x, y], box };
377}
379export function contactPatch(p: SimParams, pl: Plate): Patch {
380 const [sx, sy] = strikePoint(p, pl);
381 return patchAround(pl, sx, sy, p.contactRadius);
382}
384/** The mount patch, or null when the plate is free. */
385export function mountPatch(pl: Plate): Patch | null {
386 if (!pl.mount) return null;
387 return patchAround(pl, pl.mount.x, pl.mount.y, pl.mount.radius);
388}