1// The cymbal design: a JavaScript script that returns an object describing
2// the material, the bounding box, the outline, the thickness field, the
3// stress-free profile, the mount, and the losses. The script runs in the page
4// with a set of helper functions in scope (see HELPERS) and returns a plain
5// object; compileDesign checks it and evaluatePlate samples the fields on a
6// grid of the requested spacing.
7//
8// Lengths are meters throughout; `inch` and `mm` are available as constants.
9// Every field function receives (x, y, r, theta): the point relative to the
10// center of the bounding box, its distance from the center, and its angle
11// from the x axis in radians.
13import {
14 type Grid,
15 type InsideFn,
16 type Loss,
17 type Material,
18 type Mount,
19 type Plate,
20 insideNode,
21 makeGrid,
22 nodeWeight,
23 nodeX,
24 nodeY,
25} from "./params.ts";
27export type PointFn = (x: number, y: number, r: number, theta: number) => number;
28export type OutlineFn = (x: number, y: number, r: number, theta: number) => boolean;
30/** A compiled design: the script's object, checked and with defaults filled in. */
31export interface Design {
32 name: string;
33 material: Material;
34 Lx: number;
35 Ly: number;
36 inside: OutlineFn | null;
37 thickness: PointFn;
38 profile: PointFn;
39 mount: Mount | null;
40 loss: Loss;
41}
43export const DEFAULT_LOSS: Loss = { sigma0: 1.15, sigma1: 5e-4, sigmaIn: 20 };
45// ---------- Helpers available to the script ----------
47/** Deterministic uniform generator on [0, 1) (mulberry32). */
48function mulberry32(seed: number): () => number {
49 let a = Math.floor(seed) >>> 0;
50 return () => {
51 a = (a + 0x6d2b79f5) >>> 0;
52 let t = a;
53 t = Math.imul(t ^ (t >>> 15), t | 1);
54 t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
55 return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
56 };
57}
59const NOISE_TERMS = 64;
61/**
62 * A smooth random field of unit RMS with the covariance exp(-d^2 / length^2):
63 * a sum of plane waves whose wavevectors are Gaussian with standard deviation
64 * sqrt(2) / length per component. The same seed gives the same field on any grid.
65 */
66export function noise2d(seed: number, length: number): (x: number, y: number) => number {
67 const rnd = mulberry32(seed);
68 const gauss = () => Math.sqrt(-2 * Math.log(1 - rnd())) * Math.cos(2 * Math.PI * rnd());
69 const L = Math.max(length, 1e-4);
70 const terms = Array.from({ length: NOISE_TERMS }, () => ({ kx: (gauss() * Math.SQRT2) / L, ky: (gauss() * Math.SQRT2) / L, ph: 2 * Math.PI * rnd() }));
71 const norm = 1 / Math.sqrt(NOISE_TERMS / 2);
72 return (x, y) => {
73 let v = 0;
74 for (const t of terms) v += Math.cos(t.kx * x + t.ky * y + t.ph);
75 return v * norm;
76 };
77}
79export const HELPERS = {
80 inch: 0.0254,
81 mm: 1e-3,
82 /** Outline of a circle of the given diameter, with an optional center hole of diameter `hole`. */
83 circle: (diameter: number, hole = 0): OutlineFn => (_x, _y, r) => r <= diameter / 2 && r >= hole / 2,
84 /** Outline of the ellipse with the given width and height, with an optional center hole. */
85 ellipse: (width: number, height: number, hole = 0): OutlineFn => (x, y, r) => (x / (width / 2)) ** 2 + (y / (height / 2)) ** 2 <= 1 && r >= hole / 2,
86 /** Outline of the rectangle of the given width and height (centered), with an optional center hole. */
87 rectangle: (width: number, height: number, hole = 0): OutlineFn => (x, y, r) => Math.abs(x) <= width / 2 && Math.abs(y) <= height / 2 && r >= hole / 2,
88 clamp: (x: number, lo: number, hi: number): number => Math.min(hi, Math.max(lo, x)),
89 lerp: (a: number, b: number, t: number): number => a + (b - a) * t,
90 /** 0 below e0, 1 above e1, a smooth cubic between. */
91 smoothstep: (e0: number, e1: number, x: number): number => {
92 const t = Math.min(1, Math.max(0, (x - e0) / (e1 - e0)));
93 return t * t * (3 - 2 * t);
94 },
95 /** Gaussian bump of the given height, falling to 1/e at `radius` from the center; `d` is the distance from the center. */
96 bump: (d: number, radius: number, height: number): number => height * Math.exp(-(d * d) / (radius * radius)),
97 /** Cosine dome of the given height and base radius (zero slope at the rim and the apex); `d` is the distance from the center. */
98 dome: (d: number, radius: number, height: number): number => (d >= radius ? 0 : 0.5 * height * (1 + Math.cos((Math.PI * d) / radius))),
99 /** Parabolic bow of the given sag across a span L centered on 0; `s` is the coordinate along the span; the offset in -1..1 moves the crest toward one end. */
100 bow: (s: number, L: number, sag: number, offset = 0): number => {
101 const t = (2 * s) / L;
102 return sag * (1 - t * t) * (1 + offset * t);
103 },
104 noise2d,
105};
107export type Helpers = typeof HELPERS;
109// ---------- Compilation ----------
111/** Message of an error thrown by the script or by the checks. */
112export class DesignError extends Error {}
114function num(v: unknown, what: string, min = -Infinity): number {
115 if (typeof v !== "number" || !isFinite(v)) throw new DesignError(`${what} must be a finite number`);
116 if (v < min) throw new DesignError(`${what} must be at least ${min}`);
117 return v;
118}
120function fieldFn(v: unknown, what: string): PointFn {
121 if (typeof v === "number") {
122 if (!isFinite(v)) throw new DesignError(`${what} must be a finite number or a function`);
123 return () => v;
124 }
125 if (typeof v === "function") return v as PointFn;
126 throw new DesignError(`${what} must be a number or a function (x, y, r, theta) => number`);
127}
129/**
130 * Runs the script and checks the object it returns. The script body has the
131 * helpers in scope as free variables and must `return` the design object.
132 */
133export function compileDesign(source: string): Design {
134 let raw: unknown;
135 try {
136 const fn = new Function(...Object.keys(HELPERS), source) as (...args: unknown[]) => unknown;
137 raw = fn(...Object.values(HELPERS));
138 } catch (e) {
139 throw new DesignError(`script error: ${e instanceof Error ? e.message : String(e)}`);
140 }
141 if (!raw || typeof raw !== "object") throw new DesignError("the script must return an object (end it with `return { ... };`)");
142 const o = raw as Record<string, unknown>;
143 const name = typeof o.name === "string" && o.name.trim() ? o.name.trim() : "Untitled design";
145 const m = o.material as Record<string, unknown> | undefined;
146 if (!m || typeof m !== "object") throw new DesignError("material must be an object { E, nu, rho }");
147 const material: Material = { E: num(m.E, "material.E", 1e6), nu: num(m.nu, "material.nu", 0), rho: num(m.rho, "material.rho", 1) };
148 if (material.nu >= 0.5) throw new DesignError("material.nu must be below 0.5");
150 const b = o.bounds as Record<string, unknown> | undefined;
151 if (!b || typeof b !== "object") throw new DesignError("bounds must be an object { width, height }");
152 const Lx = num(b.width, "bounds.width", 1e-3), Ly = num(b.height, "bounds.height", 1e-3);
154 let inside: OutlineFn | null = null;
155 if (o.inside !== undefined && o.inside !== null) {
156 if (typeof o.inside !== "function") throw new DesignError("inside must be a function (x, y, r, theta) => boolean");
157 inside = o.inside as OutlineFn;
158 }
160 if (o.thickness === undefined) throw new DesignError("thickness is required (a number in meters, or a function)");
161 const thickness = fieldFn(o.thickness, "thickness");
162 const profile = o.profile === undefined || o.profile === null ? () => 0 : fieldFn(o.profile, "profile");
164 let mount: Mount | null = null;
165 if (o.mount !== undefined && o.mount !== null && o.mount !== false) {
166 const mo = o.mount as Record<string, unknown>;
167 if (typeof mo !== "object") throw new DesignError("mount must be an object { x, y, radius, stiffness, damping } or null");
168 mount = {
169 x: mo.x === undefined ? 0 : num(mo.x, "mount.x"),
170 y: mo.y === undefined ? 0 : num(mo.y, "mount.y"),
171 radius: num(mo.radius, "mount.radius", 0),
172 stiffness: num(mo.stiffness, "mount.stiffness", 0),
173 damping: mo.damping === undefined ? 0 : num(mo.damping, "mount.damping", 0),
174 };
175 }
177 const lo = (o.loss ?? {}) as Record<string, unknown>;
178 if (typeof lo !== "object") throw new DesignError("loss must be an object { sigma0, sigma1, sigmaIn }");
179 const loss: Loss = {
180 sigma0: lo.sigma0 === undefined ? DEFAULT_LOSS.sigma0 : num(lo.sigma0, "loss.sigma0", 0),
181 sigma1: lo.sigma1 === undefined ? DEFAULT_LOSS.sigma1 : num(lo.sigma1, "loss.sigma1", 0),
182 sigmaIn: lo.sigmaIn === undefined ? DEFAULT_LOSS.sigmaIn : num(lo.sigmaIn, "loss.sigmaIn", 0),
183 };
185 return { name, material, Lx, Ly, inside, thickness, profile, mount, loss };
186}
188// ---------- Evaluation on a grid ----------
190/** Subtracts the weighted least-squares plane a + b x + c y from f over the plate nodes. */
191export function removePlane(f: Float64Array, g: Grid): void {
192 let sw = 0, sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0, sf = 0, sxf = 0, syf = 0;
193 for (let j = 0; j < g.ny; j++) {
194 for (let i = 0; i < g.nx; i++) {
195 const wgt = nodeWeight(g, i, j);
196 if (wgt === 0) continue;
197 const x = nodeX(g, i), y = nodeY(g, j), v = f[j * g.nx + i];
198 sw += wgt; sx += wgt * x; sy += wgt * y;
199 sxx += wgt * x * x; sxy += wgt * x * y; syy += wgt * y * y;
200 sf += wgt * v; sxf += wgt * x * v; syf += wgt * y * v;
201 }
202 }
203 if (sw === 0) return;
204 // Normal equations of the 3-parameter fit, solved by Cramer's rule.
205 const det = sw * (sxx * syy - sxy * sxy) - sx * (sx * syy - sxy * sy) + sy * (sx * sxy - sxx * sy);
206 if (Math.abs(det) < 1e-300) return;
207 const a = (sf * (sxx * syy - sxy * sxy) - sx * (sxf * syy - sxy * syf) + sy * (sxf * sxy - sxx * syf)) / det;
208 const b = (sw * (sxf * syy - sxy * syf) - sf * (sx * syy - sxy * sy) + sy * (sx * syf - sxf * sy)) / det;
209 const c = (sw * (sxx * syf - sxf * sxy) - sx * (sx * syf - sxf * sy) + sf * (sx * sxy - sxx * sy)) / det;
210 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) f[j * g.nx + i] -= a + b * nodeX(g, i) + c * nodeY(g, j);
211}
213function callField(fn: PointFn, x: number, y: number, what: string): number {
214 let v: unknown;
215 try {
216 v = fn(x, y, Math.hypot(x, y), Math.atan2(y, x));
217 } catch (e) {
218 throw new DesignError(`${what} threw at (${x.toFixed(4)}, ${y.toFixed(4)}): ${e instanceof Error ? e.message : String(e)}`);
219 }
220 if (typeof v !== "number" || !isFinite(v)) throw new DesignError(`${what} returned ${String(v)} at (${x.toFixed(4)}, ${y.toFixed(4)})`);
221 return v;
222}
224/**
225 * Samples the design on a grid of about the requested spacing: the mask from
226 * the outline, the thickness at the nodes (which must be positive on the
227 * plate), and the profile at the nodes with its mean plane removed. Only the
228 * slopes of the profile enter the physics, so its mean is immaterial, but its
229 * tilt is not: the von Karman strains are not invariant under rotation, and a
230 * tilted stress-free shape would add a spurious membrane coupling.
231 */
232export function evaluatePlate(d: Design, spacing: number): Plate {
233 let inside: InsideFn | undefined;
234 if (d.inside) {
235 const fn = d.inside;
236 inside = (x, y) => {
237 try {
238 return !!fn(x, y, Math.hypot(x, y), Math.atan2(y, x));
239 } catch (e) {
240 throw new DesignError(`inside threw at (${x.toFixed(4)}, ${y.toFixed(4)}): ${e instanceof Error ? e.message : String(e)}`);
241 }
242 };
243 }
244 const g = makeGrid(d.Lx, d.Ly, spacing, inside);
245 const h = new Float64Array(g.n), w0 = new Float64Array(g.n);
246 let any = false;
247 for (let j = 0; j < g.ny; j++) {
248 for (let i = 0; i < g.nx; i++) {
249 if (!insideNode(g, i, j)) continue;
250 any = true;
251 const x = nodeX(g, i), y = nodeY(g, j), k = j * g.nx + i;
252 h[k] = callField(d.thickness, x, y, "thickness");
253 if (h[k] <= 0) throw new DesignError(`thickness must be positive; it is ${h[k]} at (${x.toFixed(4)}, ${y.toFixed(4)})`);
254 w0[k] = callField(d.profile, x, y, "profile");
255 }
256 }
257 if (!any) throw new DesignError("the outline contains no grid cells; enlarge it or reduce the grid spacing");
258 removePlane(w0, g);
259 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) if (!insideNode(g, i, j)) w0[j * g.nx + i] = 0;
260 return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, mount: d.mount, loss: d.loss };
261}
263/** Peak-to-peak height and the largest slope of the profile over the plate nodes. */
264export function profileStats(pl: Plate): { peakToPeak: number; maxSlope: number } {
265 const g = pl.grid, w0 = pl.w0;
266 let lo = Infinity, hi = -Infinity, maxSlope = 0;
267 const at = (i: number, j: number) => w0[j * g.nx + i];
268 const slope = (i: number, j: number, di: number, dj: number, hg: number) => {
269 const hasP = insideNode(g, i + di, j + dj), hasM = insideNode(g, i - di, j - dj);
270 if (hasP && hasM) return (at(i + di, j + dj) - at(i - di, j - dj)) / (2 * hg);
271 if (hasP) return (at(i + di, j + dj) - at(i, j)) / hg;
272 if (hasM) return (at(i, j) - at(i - di, j - dj)) / hg;
273 return 0;
274 };
275 for (let j = 0; j < g.ny; j++) {
276 for (let i = 0; i < g.nx; i++) {
277 if (!insideNode(g, i, j)) continue;
278 const v = at(i, j);
279 lo = Math.min(lo, v); hi = Math.max(hi, v);
280 maxSlope = Math.max(maxSlope, Math.hypot(slope(i, j, 1, 0, g.hx), slope(i, j, 0, 1, g.hy)));
281 }
282 }
283 return { peakToPeak: isFinite(hi - lo) ? hi - lo : 0, maxSlope };
284}