1// The design: a MATLAB-syntax script, run by numbl, that fills in a struct
2// named `design` with the material, the bounding box, the outline, the
3// thickness field, the stress-free profile, the stand, and the losses.
4// The outline, the thickness, and the profile are function handles
5// @(x, y, r, theta), vectorized over grid matrices (or plain numbers).
6// compileDesign runs the script and checks the struct; evaluatePlate runs it
7// again with a short trailer that samples the fields on meshgrid matrices of
8// the requested spacing, inside numbl, and reads the matrices back.
9//
10// Lengths are meters throughout. A script is self-contained: nothing is in
11// scope but numbl's own functions and whatever the script defines itself,
12// including ordinary MATLAB functions written at its bottom (numbl hoists
13// them, so the trailer may still follow them).
15import { executeCode, type RuntimeValue, type RuntimeTensor, type RuntimeStruct } from "numbl";
16import { RUNTIME_FILES } from "./mruntime.ts";
17import { LOSS_FREQS, constantLoss, defaultLoss, lossFromTable } from "./loss.ts";
18import {
19 type Grid,
20 type Loss,
21 type Material,
22 type Pad,
23 type Plate,
24 type Stand,
25 insideNode,
26 makeGrid,
27 nodeWeight,
28 nodeX,
29 nodeY,
30 withCellMask,
31} from "./params.ts";
33/** A compiled design: the checked contents of the `design` struct, plus the script that produced it. */
34export interface Design {
35 name: string;
36 material: Material;
37 Lx: number;
38 Ly: number;
39 /** Whether the script gave an outline; without one the plate is the whole bounding box. */
40 hasInside: boolean;
41 /** The thickness when it is a plain number; null for a function. */
42 thicknessConst: number | null;
43 /** The profile when it is a plain number (or absent); null for a function. */
44 profileConst: number | null;
45 stand: Stand | null;
46 /** The losses when they do not need the script (a number, a table, or the default); null for a function handle. */
47 lossConst: Loss | null;
48 source: string;
49}
51/** Defaults of the pad of a stand, a typical felt washer. */
52export const DEFAULT_PAD: Pad = { innerRadius: 7e-3, outerRadius: 20e-3, thickness: 12e-3, E: 0.5e6, lossFactor: 0.3 };
53export const DEFAULT_PRELOAD = 1e-3;
55/** Message of an error thrown by the script or by the checks. */
56export class DesignError extends Error {}
58const FILE_NAME = "design.m";
60/**
61 * Runs a design script and returns its workspace variables. `userLines` is
62 * the length of the part of the source the user wrote, so that an error in
63 * the appended trailer is not reported at a line the user cannot see.
64 */
65export function runDesignScript(source: string, userLines = Infinity): Record<string, RuntimeValue> {
66 try {
67 return executeCode(source, { displayResults: false, onOutput: () => {} }, RUNTIME_FILES, FILE_NAME).variableValues;
68 } catch (e) {
69 throw new DesignError(describeError(e, Math.min(source.split("\n").length, userLines)));
70 }
71}
73function describeError(e: unknown, userLines: number): string {
74 if (!(e instanceof Error)) return String(e);
75 const { line, file } = e as { line?: unknown; file?: unknown };
76 const inScript = file === undefined || file === null || file === FILE_NAME;
77 const where = inScript && typeof line === "number" && line >= 1 && line <= userLines ? `line ${line}: ` : "";
78 return `${where}${e.message}`;
79}
81// ---------- Reading runtime values ----------
83function isStruct(v: RuntimeValue | undefined): v is RuntimeStruct {
84 return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "struct";
85}
86function isTensor(v: RuntimeValue | undefined): v is RuntimeTensor {
87 return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "tensor";
88}
89function isFunction(v: RuntimeValue | undefined): boolean {
90 return typeof v === "object" && v !== null && (v as { kind?: string }).kind === "function";
91}
92function isEmpty(v: RuntimeValue | undefined): boolean {
93 return v === undefined || (isTensor(v) && v.data.length === 0);
94}
96/** A scalar number from a runtime value, or null when it is not one. */
97function scalar(v: RuntimeValue | undefined): number | null {
98 if (typeof v === "number") return v;
99 if (typeof v === "boolean") return v ? 1 : 0;
100 if (isTensor(v) && v.data.length === 1 && (!v.imag || v.imag[0] === 0)) return v.data[0];
101 return null;
102}
104function num(v: RuntimeValue | undefined, what: string, min = -Infinity): number {
105 const s = scalar(v);
106 if (s === null || !isFinite(s)) throw new DesignError(`${what} must be a finite number`);
107 if (s < min) throw new DesignError(`${what} must be at least ${min}`);
108 return s;
109}
111function bool(v: RuntimeValue | undefined, what: string, dflt: boolean): boolean {
112 if (v === undefined || isEmpty(v)) return dflt;
113 const s = scalar(v);
114 if (s === null) throw new DesignError(`${what} must be true or false`);
115 return s !== 0;
116}
118function field(s: RuntimeStruct, name: string): RuntimeValue | undefined {
119 return s.fields.get(name);
120}
122function text(v: RuntimeValue | undefined): string | null {
123 if (typeof v === "string") return v;
124 if (typeof v === "object" && v !== null && (v as { kind?: string }).kind === "char") return (v as { value: string }).value;
125 return null;
126}
128// ---------- Compilation ----------
130/**
131 * Runs the script and checks the `design` struct it leaves in the
132 * workspace. The field functions are not called here; evaluatePlate samples
133 * them on the grid.
134 */
135export function compileDesign(source: string): Design {
136 const vars = runDesignScript(source);
137 const d = vars.design;
138 if (d === undefined) throw new DesignError("the script must fill in a struct named design (design.material, design.bounds, design.thickness, ...)");
139 if (!isStruct(d)) throw new DesignError("design must be a struct");
141 const nameV = field(d, "name");
142 const name = (nameV !== undefined ? text(nameV) : null)?.trim() || "Untitled design";
144 const m = field(d, "material");
145 if (!isStruct(m)) throw new DesignError("design.material must be a struct with fields E, nu, rho");
146 const material: Material = { E: num(field(m, "E"), "design.material.E", 1e6), nu: num(field(m, "nu"), "design.material.nu", 0), rho: num(field(m, "rho"), "design.material.rho", 1) };
147 if (material.nu >= 0.5) throw new DesignError("design.material.nu must be below 0.5");
149 const b = field(d, "bounds");
150 if (!isStruct(b)) throw new DesignError("design.bounds must be a struct with fields width, height");
151 const Lx = num(field(b, "width"), "design.bounds.width", 1e-3), Ly = num(field(b, "height"), "design.bounds.height", 1e-3);
153 const inside = field(d, "inside");
154 let hasInside = false;
155 if (!isEmpty(inside)) {
156 if (!isFunction(inside)) throw new DesignError("design.inside must be a function handle @(x, y, r, theta) returning true inside the plate, for example @(x, y, r, theta) r <= D/2");
157 hasInside = true;
158 }
160 const th = field(d, "thickness");
161 if (th === undefined) throw new DesignError("design.thickness is required (a number in meters, or a function handle @(x, y, r, theta))");
162 let thicknessConst: number | null = null;
163 if (!isFunction(th)) {
164 thicknessConst = scalar(th);
165 if (thicknessConst === null || !isFinite(thicknessConst)) throw new DesignError("design.thickness must be a number or a function handle @(x, y, r, theta)");
166 if (thicknessConst <= 0) throw new DesignError("design.thickness must be positive");
167 }
169 const pr = field(d, "profile");
170 let profileConst: number | null = 0;
171 if (pr !== undefined && !isEmpty(pr)) {
172 if (isFunction(pr)) profileConst = null;
173 else {
174 profileConst = scalar(pr);
175 if (profileConst === null || !isFinite(profileConst)) throw new DesignError("design.profile must be a number or a function handle @(x, y, r, theta)");
176 }
177 }
179 const stand = readStand(field(d, "stand"));
180 const lossConst = readLoss(field(d, "loss"));
182 return { name, material, Lx, Ly, hasInside, thicknessConst, profileConst, stand, lossConst, source };
183}
185function readStand(st: RuntimeValue | undefined): Stand | null {
186 if (st === undefined || isEmpty(st) || st === false) return null;
187 if (!isStruct(st)) throw new DesignError("design.stand must be a struct (design.stand.pad, design.stand.topPad, design.stand.preload), or [] for a free plate");
188 const opt = (s: RuntimeStruct, k: string, what: string, dflt: number, min = -Infinity) => (field(s, k) === undefined || isEmpty(field(s, k)) ? dflt : num(field(s, k), what, min));
189 const padV = field(st, "pad");
190 let pad: Pad = { ...DEFAULT_PAD };
191 if (padV !== undefined && !isEmpty(padV)) {
192 if (!isStruct(padV)) throw new DesignError("design.stand.pad must be a struct with fields innerRadius, outerRadius, thickness, E, lossFactor");
193 pad = {
194 innerRadius: opt(padV, "innerRadius", "design.stand.pad.innerRadius", DEFAULT_PAD.innerRadius, 0),
195 outerRadius: opt(padV, "outerRadius", "design.stand.pad.outerRadius", DEFAULT_PAD.outerRadius, 0),
196 thickness: opt(padV, "thickness", "design.stand.pad.thickness", DEFAULT_PAD.thickness, 1e-5),
197 E: opt(padV, "E", "design.stand.pad.E", DEFAULT_PAD.E, 1),
198 lossFactor: opt(padV, "lossFactor", "design.stand.pad.lossFactor", DEFAULT_PAD.lossFactor, 0),
199 };
200 if (pad.outerRadius <= pad.innerRadius) throw new DesignError("design.stand.pad.outerRadius must exceed innerRadius");
201 }
202 return {
203 x: opt(st, "x", "design.stand.x", 0),
204 y: opt(st, "y", "design.stand.y", 0),
205 pad,
206 topPad: bool(field(st, "topPad"), "design.stand.topPad", true),
207 preload: opt(st, "preload", "design.stand.preload", DEFAULT_PRELOAD, 0),
208 };
209}
211/** The losses when they can be read without the script: null when T60 is a function handle. */
212function readLoss(lo: RuntimeValue | undefined): Loss | null {
213 if (lo === undefined || isEmpty(lo)) return defaultLoss();
214 if (!isStruct(lo)) throw new DesignError("design.loss must be a struct with the field T60: a decay time in seconds, a table [f T60; ...], or a function handle @(f)");
215 const t = field(lo, "T60");
216 if (t === undefined || isEmpty(t)) return defaultLoss();
217 if (isFunction(t)) return null;
218 const s = scalar(t);
219 if (s !== null) {
220 if (!(s > 0) || !isFinite(s)) throw new DesignError("design.loss.T60 must be a positive decay time in seconds");
221 return constantLoss(s);
222 }
223 if (isTensor(t) && t.shape.length === 2 && t.shape[1] === 2 && t.shape[0] >= 1) {
224 const rows = t.shape[0];
225 const fs: number[] = [], ts: number[] = [];
226 for (let r = 0; r < rows; r++) {
227 const f = t.data[r], t60 = t.data[rows + r]; // column-major
228 if (!(f > 0) || !(t60 > 0) || !isFinite(f) || !isFinite(t60)) throw new DesignError("design.loss.T60 table rows must be [frequency T60] with positive values");
229 fs.push(f); ts.push(t60);
230 }
231 return lossFromTable(fs, ts);
232 }
233 throw new DesignError("design.loss.T60 must be a number, a two-column table [f T60; ...], or a function handle @(f)");
234}
236// ---------- Evaluation on a grid ----------
238/** Subtracts the weighted least-squares plane a + b x + c y from f over the plate nodes. */
239export function removePlane(f: Float64Array, g: Grid): void {
240 let sw = 0, sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0, sf = 0, sxf = 0, syf = 0;
241 for (let j = 0; j < g.ny; j++) {
242 for (let i = 0; i < g.nx; i++) {
243 const wgt = nodeWeight(g, i, j);
244 if (wgt === 0) continue;
245 const x = nodeX(g, i), y = nodeY(g, j), v = f[j * g.nx + i];
246 sw += wgt; sx += wgt * x; sy += wgt * y;
247 sxx += wgt * x * x; sxy += wgt * x * y; syy += wgt * y * y;
248 sf += wgt * v; sxf += wgt * x * v; syf += wgt * y * v;
249 }
250 }
251 if (sw === 0) return;
252 const det = sw * (sxx * syy - sxy * sxy) - sx * (sx * syy - sxy * sy) + sy * (sx * sxy - sxx * sy);
253 if (Math.abs(det) < 1e-300) return;
254 const a = (sf * (sxx * syy - sxy * sxy) - sx * (sxf * syy - sxy * syf) + sy * (sxf * sxy - sxx * syf)) / det;
255 const b = (sw * (sxf * syy - sxy * syf) - sf * (sx * syy - sxy * sy) + sy * (sx * syf - sxf * sy)) / det;
256 const c = (sw * (sxx * syf - sxf * sxy) - sx * (sx * syf - sxf * sy) + sf * (sx * sxy - sxx * sy)) / det;
257 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);
258}
260/** The trailer appended to the script to sample the fields on the grid (all names start with cymbalsim_). */
261function trailer(d: Design, g: Grid): string {
262 const fmt = (v: number) => v.toPrecision(17);
263 return [
264 "",
265 "% ---- appended by the simulator: sample the fields on the grid ----",
266 `cymbalsim_nx = ${g.nx}; cymbalsim_ny = ${g.ny};`,
267 `cymbalsim_hx = ${fmt(g.hx)}; cymbalsim_hy = ${fmt(g.hy)};`,
268 `[cymbalsim_X, cymbalsim_Y] = meshgrid(((0:cymbalsim_nx-1) - (cymbalsim_nx-1)/2) * cymbalsim_hx, ((0:cymbalsim_ny-1) - (cymbalsim_ny-1)/2) * cymbalsim_hy);`,
269 "cymbalsim_R = hypot(cymbalsim_X, cymbalsim_Y); cymbalsim_T = atan2(cymbalsim_Y, cymbalsim_X);",
270 "cymbalsim_H = cymbalsim_field(design.thickness, 'thickness', cymbalsim_X, cymbalsim_Y, cymbalsim_R, cymbalsim_T);",
271 d.profileConst === null
272 ? "cymbalsim_Z = cymbalsim_field(design.profile, 'profile', cymbalsim_X, cymbalsim_Y, cymbalsim_R, cymbalsim_T);"
273 : `cymbalsim_Z = ${fmt(d.profileConst)} * ones(cymbalsim_ny, cymbalsim_nx);`,
274 d.hasInside
275 ? [
276 `[cymbalsim_XC, cymbalsim_YC] = meshgrid(((0:cymbalsim_nx-2) - (cymbalsim_nx-2)/2) * cymbalsim_hx, ((0:cymbalsim_ny-2) - (cymbalsim_ny-2)/2) * cymbalsim_hy);`,
277 "cymbalsim_IN = cymbalsim_field(design.inside, 'inside', cymbalsim_XC, cymbalsim_YC, hypot(cymbalsim_XC, cymbalsim_YC), atan2(cymbalsim_YC, cymbalsim_XC));",
278 ].join("\n")
279 : "",
280 d.lossConst === null
281 ? [`cymbalsim_F = [${Array.from(LOSS_FREQS, fmt).join(" ")}];`, "cymbalsim_T60 = cymbalsim_t60(design.loss.T60, cymbalsim_F);"].join("\n")
282 : "",
283 "",
284 ].join("\n");
285}
287/** A real [rows, cols] matrix from the workspace as a Float64Array in column-major order. */
288function matrix(vars: Record<string, RuntimeValue>, name: string, rows: number, cols: number, what: string): Float64Array {
289 const v = vars[name];
290 if (!isTensor(v) || v.shape.length !== 2 || v.shape[0] !== rows || v.shape[1] !== cols) throw new DesignError(`${what} did not produce a ${rows} x ${cols} matrix`);
291 if (v.imag && v.imag.some((x) => x !== 0)) throw new DesignError(`${what} returned complex values`);
292 return v.data;
293}
295/**
296 * Samples the design on a grid of about the requested spacing: the mask from
297 * the outline at the cell centers, the thickness at the nodes (which must be
298 * positive on the plate), the profile at the nodes with its mean plane
299 * removed, and the decay times when they are a function. Only the slopes of
300 * the profile enter the physics, so its mean is immaterial, but its tilt is
301 * not: the von Karman strains are not invariant under rotation, and a tilted
302 * stress-free shape would add a spurious membrane coupling.
303 */
304export function evaluatePlate(d: Design, spacing: number): Plate {
305 const full = makeGrid(d.Lx, d.Ly, spacing);
306 const { nx, ny } = full;
307 const vars = runDesignScript(d.source + trailer(d, full), d.source.split("\n").length);
308 // meshgrid matrices are [ny, nx] in column-major order: entry (i, j) sits at i * ny + j.
309 const H = matrix(vars, "cymbalsim_H", ny, nx, "thickness");
310 const Z = matrix(vars, "cymbalsim_Z", ny, nx, "profile");
311 let g = full;
312 if (d.hasInside) {
313 const IN = matrix(vars, "cymbalsim_IN", ny - 1, nx - 1, "inside");
314 g = withCellMask(full, (ci, cj) => IN[ci * (ny - 1) + cj] !== 0);
315 }
316 let loss = d.lossConst;
317 if (loss === null) {
318 const T = matrix(vars, "cymbalsim_T60", 1, LOSS_FREQS.length, "design.loss.T60");
319 for (let i = 0; i < T.length; i++) if (!(T[i] > 0) || !isFinite(T[i])) throw new DesignError(`design.loss.T60 must be positive; it is ${T[i]} at ${LOSS_FREQS[i].toFixed(0)} Hz`);
320 loss = { freqs: LOSS_FREQS, t60: Float64Array.from(T) };
321 }
322 const h = new Float64Array(g.n), w0 = new Float64Array(g.n);
323 let any = false;
324 for (let j = 0; j < ny; j++) {
325 for (let i = 0; i < nx; i++) {
326 if (!insideNode(g, i, j)) continue;
327 any = true;
328 const k = j * nx + i, m = i * ny + j;
329 const x = nodeX(g, i), y = nodeY(g, j);
330 const at = `(${x.toFixed(4)}, ${y.toFixed(4)})`;
331 if (!isFinite(H[m])) throw new DesignError(`thickness is ${H[m]} at ${at}`);
332 if (H[m] <= 0) throw new DesignError(`thickness must be positive; it is ${H[m]} at ${at}`);
333 if (!isFinite(Z[m])) throw new DesignError(`profile is ${Z[m]} at ${at}`);
334 h[k] = H[m];
335 w0[k] = Z[m];
336 }
337 }
338 if (!any) throw new DesignError("the outline contains no grid cells; enlarge it or reduce the grid spacing");
339 removePlane(w0, g);
340 for (let j = 0; j < ny; j++) for (let i = 0; i < nx; i++) if (!insideNode(g, i, j)) w0[j * nx + i] = 0;
341 return { name: d.name, Lx: d.Lx, Ly: d.Ly, grid: g, mat: d.material, h, w0, stand: d.stand, loss };
342}
344/** Whether the plate is flat (no profile), in which case bending and in-plane motion decouple in the linear model. */
345export function isFlat(pl: Plate): boolean {
346 for (let k = 0; k < pl.w0.length; k++) if (pl.w0[k] !== 0) return false;
347 return true;
348}
350/** Peak-to-peak height and the largest slope of the profile over the plate nodes. */
351export function profileStats(pl: Plate): { peakToPeak: number; maxSlope: number } {
352 const g = pl.grid, w0 = pl.w0;
353 let lo = Infinity, hi = -Infinity, maxSlope = 0;
354 const at = (i: number, j: number) => w0[j * g.nx + i];
355 const slope = (i: number, j: number, di: number, dj: number, hg: number) => {
356 const hasP = insideNode(g, i + di, j + dj), hasM = insideNode(g, i - di, j - dj);
357 if (hasP && hasM) return (at(i + di, j + dj) - at(i - di, j - dj)) / (2 * hg);
358 if (hasP) return (at(i + di, j + dj) - at(i, j)) / hg;
359 if (hasM) return (at(i, j) - at(i - di, j - dj)) / hg;
360 return 0;
361 };
362 for (let j = 0; j < g.ny; j++) {
363 for (let i = 0; i < g.nx; i++) {
364 if (!insideNode(g, i, j)) continue;
365 const v = at(i, j);
366 lo = Math.min(lo, v); hi = Math.max(hi, v);
367 maxSlope = Math.max(maxSlope, Math.hypot(slope(i, j, 1, 0, g.hx), slope(i, j, 0, 1, g.hy)));
368 }
369 }
370 return { peakToPeak: isFinite(hi - lo) ? hi - lo : 0, maxSlope };
371}