1// Convert the `sys` struct produced by a numbl script into a BeadSystem.
2// The script uses 1-based bead indices; they become 0-based here.
4import type { RuntimeValue } from "numbl";
5import { hingeAngle } from "./hinge";
6import type { BeadSystem, Hit, ListenQuantity, Method } from "./types";
8type Mat = { rows: number; cols: number; data: Float64Array };
10class SysError extends Error {}
12export const DEFAULT_T60 = 2.5; // s
14function asMat(v: RuntimeValue | undefined, name: string): Mat | null {
15 if (v === undefined) return null;
16 if (typeof v === "number") return { rows: 1, cols: 1, data: Float64Array.of(v) };
17 if (typeof v === "boolean") return { rows: 1, cols: 1, data: Float64Array.of(v ? 1 : 0) };
18 if (typeof v === "object" && v !== null && (v as { kind?: string }).kind === "tensor") {
19 const t = v as unknown as { data: Float64Array; shape: number[]; imag?: Float64Array };
20 if (t.imag) throw new SysError(`${name} must be real`);
21 const rows = t.shape[0] ?? 1;
22 const cols = t.shape.slice(1).reduce((a, b) => a * b, 1);
23 return { rows, cols, data: t.data };
24 }
25 throw new SysError(`${name} must be numeric`);
26}
28function scalar(v: RuntimeValue | undefined, name: string, def: number): number {
29 const m = asMat(v, name);
30 if (!m || m.data.length === 0) return def;
31 if (m.data.length !== 1) throw new SysError(`${name} must be a scalar`);
32 return m.data[0];
33}
35function vec3(v: RuntimeValue | undefined, name: string, def: [number, number, number]): [number, number, number] {
36 const m = asMat(v, name);
37 if (!m || m.data.length === 0) return def;
38 if (m.data.length !== 3) throw new SysError(`${name} must have 3 elements`);
39 return [m.data[0], m.data[1], m.data[2]];
40}
42function str(v: RuntimeValue | undefined, name: string, def: string): string {
43 if (v === undefined) return def;
44 if (typeof v === "string") return v;
45 if (typeof v === "object" && v !== null && (v as { kind?: string }).kind === "char") {
46 return (v as unknown as { value: string }).value;
47 }
48 throw new SysError(`${name} must be a string`);
49}
51function perBead(v: RuntimeValue | undefined, name: string, n: number, def: number): Float64Array {
52 const out = new Float64Array(n);
53 const m = asMat(v, name);
54 if (!m || m.data.length === 0) return out.fill(def);
55 if (m.data.length === 1) return out.fill(m.data[0]);
56 if (m.data.length !== n) throw new SysError(`${name} must be a scalar or have one entry per bead (${n})`);
57 out.set(m.data);
58 return out;
59}
61function indices(v: RuntimeValue | undefined, name: string, n: number): Int32Array {
62 const m = asMat(v, name);
63 if (!m) return new Int32Array(0);
64 const isLogical = (v as { _isLogical?: boolean })._isLogical || typeof v === "boolean";
65 if (isLogical) {
66 if (m.data.length !== n) throw new SysError(`${name}: a logical mask must have one entry per bead (${n})`);
67 const out: number[] = [];
68 m.data.forEach((b, i) => b && out.push(i));
69 return Int32Array.from(out);
70 }
71 const out = new Int32Array(m.data.length);
72 m.data.forEach((b, q) => {
73 if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`${name}: bead index ${b} is not in 1..${n}`);
74 out[q] = b - 1;
75 });
76 return out;
77}
79function fields(v: RuntimeValue | undefined, name: string): Map<string, RuntimeValue>[] {
80 if (v === undefined) return [];
81 const k = (v as { kind?: string }).kind;
82 if (k === "struct") return [(v as unknown as { fields: Map<string, RuntimeValue> }).fields];
83 if (k === "struct_array") {
84 return (v as unknown as { elements: { fields: Map<string, RuntimeValue> }[] }).elements.map(e => e.fields);
85 }
86 throw new SysError(`${name} must be a struct or struct array`);
87}
89const KNOWN = new Set([
90 "pos", "mass", "pinned", "springs", "bends", "hinges", "dashpots", "friction", "hits", "listen",
91 "duration", "sampleRate", "dt", "safety", "method", "frameRate", "modes", "decayLaw",
92]);
94export function parseSystem(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warnings: string[] } {
95 try {
96 return parseInner(sysValue);
97 } catch (e) {
98 if (e instanceof SysError) throw new Error(e.message);
99 throw e;
100 }
101}
103function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warnings: string[] } {
104 if (sysValue === undefined) throw new SysError("The script must define a struct named sys.");
105 const [f] = fields(sysValue, "sys");
106 if (!f) throw new SysError("sys must be a single struct");
107 const warnings: string[] = [];
108 for (const key of f.keys()) if (!KNOWN.has(key)) warnings.push(`sys.${key} is not used by the solver.`);
110 const P = asMat(f.get("pos"), "sys.pos");
111 if (!P || P.cols !== 3 || P.rows < 1) throw new SysError("sys.pos must be an n-by-3 matrix of bead positions");
112 const n = P.rows;
113 const pos0 = new Float64Array(3 * n);
114 for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) pos0[3 * i + c] = P.data[i + c * n];
116 if (!f.has("mass")) throw new SysError("sys.mass is required (scalar or one per bead)");
117 const mass = perBead(f.get("mass"), "sys.mass", n, 0);
118 for (let i = 0; i < n; i++) if (!(mass[i] > 0)) throw new SysError(`sys.mass must be positive (bead ${i + 1})`);
120 const pinned = new Uint8Array(n);
121 for (const b of indices(f.get("pinned"), "sys.pinned", n)) pinned[b] = 1;
123 const S = asMat(f.get("springs"), "sys.springs");
124 const ns = S && S.data.length ? S.rows : 0;
125 if (S && ns && S.cols !== 4 && S.cols !== 5) {
126 throw new SysError("sys.springs must have 4 or 5 columns: [i j k L0] or [i j k L0 k3]");
127 }
128 const si = new Int32Array(ns);
129 const sj = new Int32Array(ns);
130 const sk = new Float64Array(ns);
131 const sL0 = new Float64Array(ns);
132 const sk3 = new Float64Array(ns);
133 for (let s = 0; s < ns; s++) {
134 const col = (c: number) => S!.data[s + c * ns];
135 const i = col(0);
136 const j = col(1);
137 for (const b of [i, j]) {
138 if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.springs row ${s + 1}: bead index ${b} is not in 1..${n}`);
139 }
140 if (i === j) throw new SysError(`sys.springs row ${s + 1} joins bead ${i} to itself`);
141 si[s] = i - 1;
142 sj[s] = j - 1;
143 sk[s] = col(2);
144 sL0[s] = col(3);
145 sk3[s] = S!.cols === 5 ? col(4) : 0;
146 if (sk[s] < 0 || sk3[s] < 0 || sL0[s] < 0) throw new SysError(`sys.springs row ${s + 1}: k, L0 and k3 must be nonnegative`);
147 }
149 const D = asMat(f.get("dashpots"), "sys.dashpots");
150 const nd = D && D.data.length ? D.rows : 0;
151 if (D && nd && D.cols !== 3 && D.cols !== 4) throw new SysError("sys.dashpots must have 3 or 4 columns: [i j c] or [i j c ca]");
152 const di = new Int32Array(nd);
153 const dj = new Int32Array(nd);
154 const dc = new Float64Array(nd);
155 const dca = new Float64Array(nd);
156 for (let d = 0; d < nd; d++) {
157 const col = (c: number) => D!.data[d + c * nd];
158 for (const b of [col(0), col(1)]) {
159 if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.dashpots row ${d + 1}: bead index ${b} is not in 1..${n}`);
160 }
161 di[d] = col(0) - 1;
162 dj[d] = col(1) - 1;
163 dc[d] = col(2);
164 dca[d] = D!.cols === 4 ? col(3) : 0;
165 if (dca[d] < 0) throw new SysError(`sys.dashpots row ${d + 1}: ca must be nonnegative`);
166 if (dc[d] < 0) throw new SysError(`sys.dashpots row ${d + 1}: c must be nonnegative`);
167 }
168 const B = asMat(f.get("bends"), "sys.bends");
169 const nbend = B && B.data.length ? B.rows : 0;
170 if (B && nbend && B.cols !== 4 && B.cols !== 5) {
171 throw new SysError("sys.bends must have 4 or 5 columns: [i j k kb] or [i j k kb cb] (j the middle bead)");
172 }
173 const bi = new Int32Array(nbend);
174 const bj = new Int32Array(nbend);
175 const bk = new Int32Array(nbend);
176 const bkb = new Float64Array(nbend);
177 const bcb = new Float64Array(nbend);
178 for (let r = 0; r < nbend; r++) {
179 const col = (c: number) => B!.data[r + c * nbend];
180 const ids = [col(0), col(1), col(2)];
181 for (const b of ids) {
182 if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.bends row ${r + 1}: bead index ${b} is not in 1..${n}`);
183 }
184 if (new Set(ids).size !== 3) throw new SysError(`sys.bends row ${r + 1} repeats a bead`);
185 bi[r] = ids[0] - 1;
186 bj[r] = ids[1] - 1;
187 bk[r] = ids[2] - 1;
188 bkb[r] = col(3);
189 bcb[r] = B!.cols === 5 ? col(4) : 0;
190 if (bkb[r] < 0 || bcb[r] < 0) throw new SysError(`sys.bends row ${r + 1}: kb and cb must be nonnegative`);
191 }
193 const Hm = asMat(f.get("hinges"), "sys.hinges");
194 const nh = Hm && Hm.data.length ? Hm.rows : 0;
195 if (Hm && nh && Hm.cols !== 5 && Hm.cols !== 6) {
196 throw new SysError("sys.hinges must have 5 or 6 columns: [i j k l kh] or [i j k l kh ch] (j-k the shared edge, i and l the opposite beads)");
197 }
198 const hi = new Int32Array(nh);
199 const hj = new Int32Array(nh);
200 const hk = new Int32Array(nh);
201 const hl = new Int32Array(nh);
202 const hkh = new Float64Array(nh);
203 const hch = new Float64Array(nh);
204 const hth0 = new Float64Array(nh);
205 const gtmp = new Float64Array(12);
206 for (let r = 0; r < nh; r++) {
207 const col = (c: number) => Hm!.data[r + c * nh];
208 const ids = [col(0), col(1), col(2), col(3)];
209 for (const b of ids) {
210 if (!Number.isInteger(b) || b < 1 || b > n) throw new SysError(`sys.hinges row ${r + 1}: bead index ${b} is not in 1..${n}`);
211 }
212 if (new Set(ids).size !== 4) throw new SysError(`sys.hinges row ${r + 1} repeats a bead`);
213 hi[r] = ids[0] - 1;
214 hj[r] = ids[1] - 1;
215 hk[r] = ids[2] - 1;
216 hl[r] = ids[3] - 1;
217 hkh[r] = col(4);
218 hch[r] = Hm!.cols === 6 ? col(5) : 0;
219 if (hkh[r] < 0 || hch[r] < 0) throw new SysError(`sys.hinges row ${r + 1}: kh and ch must be nonnegative`);
220 hth0[r] = hingeAngle(pos0, hi[r], hj[r], hk[r], hl[r], gtmp);
221 if (!Number.isFinite(hth0[r]) || !gtmp.every(Number.isFinite)) {
222 throw new SysError(`sys.hinges row ${r + 1}: a triangle has zero area at the rest positions`);
223 }
224 }
226 let friction = perBead(f.get("friction"), "sys.friction", n, 0);
227 if (!f.has("friction") && !f.has("dashpots") && !(B && B.cols === 5) && !(Hm && Hm.cols === 6)) {
228 // Default damping: every mode falls 60 dB in DEFAULT_T60 seconds
229 // (amplitude rate sigma = ln(1000)/T60, from friction 2 sigma m per bead).
230 const sigma = Math.log(1000) / DEFAULT_T60;
231 friction = mass.map(m => 2 * sigma * m);
232 warnings.push(`No damping given: using the default, every mode falls 60 dB in ${DEFAULT_T60} s. Set sys.friction = 0 for none.`);
233 }
235 const hits: Hit[] = fields(f.get("hits"), "sys.hits").map((h, q) => {
236 const nm = `sys.hits(${q + 1})`;
237 const beads = indices(h.get("beads"), `${nm}.beads`, n);
238 if (beads.length === 0) throw new SysError(`${nm}.beads is required`);
239 const duration = scalar(h.get("duration"), `${nm}.duration`, 0);
240 const hit: Hit = {
241 t: scalar(h.get("t"), `${nm}.t`, 0),
242 beads,
243 velocity: vec3(h.get("velocity"), `${nm}.velocity`, [0, 0, 0]),
244 force: vec3(h.get("force"), `${nm}.force`, [0, 0, 0]),
245 duration,
246 };
247 // In a struct array every hit has every field, so treat empty [] as absent.
248 const given = (k: string) => (asMat(h.get(k), `${nm}.${k}`)?.data.length ?? 0) > 0;
249 if (duration > 0 && given("velocity")) throw new SysError(`${nm}: give velocity (impulse) or force with duration (pulse), not both`);
250 if (duration === 0 && given("force")) throw new SysError(`${nm}.force needs a positive duration`);
251 return hit;
252 });
253 if (hits.length === 0) warnings.push("sys.hits is empty: nothing excites the system.");
255 const [L] = fields(f.get("listen"), "sys.listen");
256 if (!L) throw new SysError("sys.listen is required, e.g. struct('beads', 10, 'dir', [0 0 1])");
257 const listenBeads = indices(L.get("beads") ?? L.get("bead"), "sys.listen.beads", n);
258 if (listenBeads.length === 0) throw new SysError("sys.listen.beads is required");
259 const dir = vec3(L.get("dir"), "sys.listen.dir", [0, 0, 1]);
260 const listenQuantity = str(L.get("quantity"), "sys.listen.quantity", "acceleration") as ListenQuantity;
261 if (!["acceleration", "velocity", "displacement"].includes(listenQuantity)) {
262 throw new SysError("sys.listen.quantity must be 'acceleration' (default), 'velocity' or 'displacement'");
263 }
264 const dn = Math.hypot(...dir);
265 if (dn === 0) throw new SysError("sys.listen.dir must be nonzero");
267 const method = str(f.get("method"), "sys.method", "verlet") as Method;
268 if (method !== "verlet" && method !== "yoshida4") throw new SysError("sys.method must be 'verlet' or 'yoshida4'");
270 const sys: BeadSystem = {
271 n, pos0, mass, pinned, si, sj, sk, sL0, sk3, bi, bj, bk, bkb, bcb, hi, hj, hk, hl, hkh, hch, hth0, di, dj, dc, dca, friction, hits, listenBeads,
272 listenDir: [dir[0] / dn, dir[1] / dn, dir[2] / dn],
273 listenQuantity,
274 duration: scalar(f.get("duration"), "sys.duration", 2),
275 sampleRate: scalar(f.get("sampleRate"), "sys.sampleRate", 48000),
276 dt: scalar(f.get("dt"), "sys.dt", 0),
277 safety: scalar(f.get("safety"), "sys.safety", 0.5),
278 method,
279 frameRate: scalar(f.get("frameRate"), "sys.frameRate", 500),
280 modes: scalar(f.get("modes"), "sys.modes", 1) !== 0,
281 };
282 if (!(sys.duration > 0)) throw new SysError("sys.duration must be positive");
283 if (!(sys.safety > 0 && sys.safety <= 1)) throw new SysError("sys.safety must be in (0, 1]");
284 if (ns === 0 && nbend === 0 && nh === 0) warnings.push("sys.springs, sys.bends and sys.hinges are empty: the beads do not interact.");
285 return { sys, warnings };
286}