Triangle strike simulator: MATLAB model and interactive web app
44 changed files+8375−0
.github/workflows/deploy.ymladded+69−0View file
@@ -0,0 +1,69 @@
1+name: Deploy to Pages
2+
3+# Builds webapp/ and publishes it as this repository's site.
4+#
5+# Two notes on the vault:
6+# - there is no cache service, so setup-node runs without `cache: npm`;
7+# - a mochi site is served at /<collection>/<repo>/site/ or at the root of an
8+# origin of its own. The Vite build uses a relative base (`base: './'`), so
9+# it works at either without being told the deploy path, and configure-pages
10+# is therefore only here for parity with the deploy step.
11+
12+on:
13+ push:
14+ branches: [main]
15+ workflow_dispatch:
16+
17+permissions:
18+ contents: read
19+ pages: write
20+ id-token: write
21+
22+concurrency:
23+ group: pages
24+ cancel-in-progress: false
25+
26+jobs:
27+ build:
28+ runs-on: ubuntu-latest
29+ steps:
30+ - uses: actions/checkout@v4
31+
32+ - uses: actions/setup-node@v4
33+ with:
34+ node-version: 22
35+
36+ - name: Install
37+ working-directory: webapp
38+ run: npm ci
39+
40+ # Runs matlab/triangle_validate.m through numbl: the straight-rod
41+ # closed forms, the eigensolve identities, reproducibility of the
42+ # seeded roll, and the pinned end-to-end diagnostics. Fails the build
43+ # on any mismatch.
44+ - name: Validate the model
45+ working-directory: webapp
46+ run: npm run selftest
47+
48+ - name: Build
49+ working-directory: webapp
50+ run: npm run build
51+
52+ - name: Setup Pages
53+ uses: actions/configure-pages@v5
54+
55+ - name: Upload artifact
56+ uses: actions/upload-pages-artifact@v3
57+ with:
58+ path: webapp/dist
59+
60+ deploy:
61+ needs: build
62+ runs-on: ubuntu-latest
63+ environment:
64+ name: github-pages
65+ url: ${{ steps.deployment.outputs.page_url }}
66+ steps:
67+ - name: Deploy to Pages
68+ id: deployment
69+ uses: actions/deploy-pages@v4
README.mdadded+45−0View file
@@ -0,0 +1,45 @@
1+# Triangle
2+
3+A physically-based simulation of an orchestral triangle strike, sampled as acoustic pressure at a point in space. The instrument is a steel rod bent into an open equilateral triangle; a Hertz beater contact excites the bending, stretching, and torsional modes of the bent free-free rod, computed by a small finite element eigensolve, and each element of the rod radiates as a point dipole to a listener point.
4+
5+There is no microphone model: no directivity, no capsule response, no gain staging. The output is p(r, t) in pascals.
6+
7+## Contents
8+
9+| Path | What it is |
10+|---|---|
11+| `matlab/` | The algorithm in MATLAB syntax, run by [numbl](https://numbl.org) |
12+| `webapp/` | Interactive browser app that executes the MATLAB code client-side |
13+
14+The MATLAB code runs from the command line with `cd matlab && numbl run triangle_demo.m`, or `triangle_validate.m` to run the checks. The web app is a Vite project; see `webapp/README.md`.
15+
16+## The model
17+
18+The signal chain, with each stage in its own `.m` file:
19+
20+```
21+strike or roll gesture -> Hertz beater contacts triangle_strikes.m
22+ -> bending / stretching / torsional modes of the
23+ open bent rod (finite elements + eigensolve) triangle_basis.m
24+ -> coherent sum of element dipoles triangle_radiate.m
25+ -> free-field pressure at the listener
26+ -> anti-aliased decimation -> PCM triangle_output.m
27+```
28+
29+Three points are worth stating plainly, because they are where the model commits itself.
30+
31+**The modes come from an honest discretization of the bent geometry.** The rod centerline is an open polyline (two closed corners, a gap at the third), discretized into about sixty Euler-Bernoulli beam elements with consistent mass and rotary inertia (`triangle_fem.m`). Because the geometry is planar, the six degrees of freedom per node decouple exactly into two families: in-plane bending plus stretching, and out-of-plane bending plus torsion. The bends are where the families earn their keep: a corner converts out-of-plane bending in one leg into torsion in the next, and it is this conversion, not any hand-tuned inharmonicity, that produces the dense non-overlapping cluster of partials a triangle is heard as. Each family is a generalized eigenproblem, reduced through the Cholesky factor of the mass matrix and diagonalized by singular value decomposition (the reduced matrix is symmetric positive semidefinite, and svd keeps the near-degenerate pairs of this geometry orthonormal where a nonsymmetric eigensolver does not). Shear deformation is neglected, which leaves the modes near 20 kHz a few percent sharp for an 8 mm rod. The rod is fully free-free: the suspension loop is taken as massless and lossless, so the lowest modes ring longer here than on a hand-held instrument, where the loop and the hand take their energy.
32+
33+**The strike is an event, not a simulation.** Each beater contact is a Hertz half-sine force pulse whose contact time couples the beater mass to an estimate of the rod mass engaged during the contact, closed by a short fixed-point iteration (`triangle_strikes.m`). What this buys is closed-form synthesis: between strikes every mode is a damped exponential, nothing is integrated, and a five-second ring renders in about two seconds. What it gives up is the millisecond-scale detail of the contact itself; unlike the sibling `woodblock` repository, which integrates its mallet contact, the force pulse here is prescribed. A roll is a seeded point process of such events alternating between two contact points, reproducible from its seed.
34+
35+**Radiation is a coherent sum of element dipoles.** A thin rod moving transversally radiates as a line of point dipoles (dipole volume twice the displaced volume, the added mass coefficient of a cylinder). Per mode, the element dipoles are summed with their relative retardations at the mode's centre frequency, so the interference across the frame, and its dependence on where the listener stands, is kept. Radiation damping is estimated from the incoherent sum of element dipole powers, fair above the rod's coincidence frequency (about 1.9 kHz here) and an overestimate below it. Modes above `max_mode_hz` are truncated, and their collective off-resonance response is reinstated as a broadband click from the driving-point impedance of an infinite beam; that is the least principled term of the model, and `direct_gain` turns it off.
36+
37+## Validation
38+
39+`matlab/triangle_validate.m` (also `npm run selftest` in `webapp/`) checks the element routines against the closed-form free-free bar (bending, axial, and torsional frequencies, and the rigid-body count), the mass-orthonormality and stiffness-diagonalization identities of the eigensolve, the Hertz `v^(-1/5)` contact-time scaling in the light-beater limit, the far-field `-6 dB` per doubling of distance along a fixed ray, exact reproducibility of the seeded roll, and pinned end-to-end diagnostics.
40+
41+As in the sibling `maraca` repository, there is no separate Python reference implementation: the MATLAB code is the reference, and the pinned numbers were recorded from it when the model was built.
42+
43+## Defaults
44+
45+A 15.2 cm (6 inch) steel triangle of 8 mm rod with a 15 mm corner opening, struck once at 1.5 m/s on the middle of the base with a 22 g steel beater, 20 degrees off the plane. That yields 40 modes between 164 Hz and 21.6 kHz (19 in-plane, 21 out-of-plane), a 57 µs contact, and a peak of 7.2 Pa (111 dB SPL) at the 45 cm listener before normalization, ringing through a five-second tail.
matlab/triangle_basis.madded+161−0View file
@@ -0,0 +1,161 @@
1+function basis = triangle_basis(P)
2+% TRIANGLE_BASIS Mode set of the open triangle: FEM assembly and eigensolve.
3+%
4+% The planar geometry splits the frame into two exactly decoupled families
5+% (see triangle_fem.m): in-plane modes (bending + stretching) and
6+% out-of-plane modes (bending + torsion, converted into each other at the
7+% two closed corners). Each family gives a generalized eigenproblem
8+% K phi = om^2 M phi, reduced to symmetric standard form through the
9+% Cholesky factor of the consistent mass matrix and diagonalized by svd
10+% rather than eig: the reduced matrix is symmetric positive semidefinite,
11+% so its singular vectors are its eigenvectors, and svd keeps them
12+% orthonormal even for the near-degenerate pairs this geometry produces
13+% (eig's nonsymmetric path returns complex, non-orthogonal vectors there).
14+% The eigenvectors then come back mass-normalized: phi' M phi = I. The
15+% three rigid-body modes of each family sit at zero frequency and are
16+% dropped along with anything below min_mode_hz; modes above max_mode_hz
17+% are truncated (the broadband part they carry is reinstated by the direct
18+% term in triangle_radiate.m).
19+%
20+% For radiation and excitation the basis stores, per retained mode, the
21+% displacement at every node (NUx, NUy, NUz; out-of-plane modes have only
22+% NUz, in-plane modes only NUx, NUy) and the transverse displacement at
23+% every element center (EUpx, EUpy, EUpz), with the axial component removed:
24+% axial motion of a thin rod displaces almost no air. Element centers,
25+% tangents, and lengths are kept alongside so the radiation stage can treat
26+% each element as a point dipole.
27+%
28+% Damping: zeta = 0.5 (loss_factor + loss_slope om) plus radiation damping
29+% from the incoherent sum of the element dipole powers,
30+% zeta_rad = rho0 om^3 sum_e (Vdip_e^2 |u_perp,e|^2) / (24 pi c^3),
31+% Vdip_e = 2 pi a^2 l_e (the factor 2 is the added-mass coefficient of a
32+% transversally moving cylinder). The incoherent sum ignores the
33+% cancellation between anti-phase sections, which matters only below the
34+% coincidence frequency (about 1.9 kHz for an 8 mm steel rod); above it the
35+% flexural waves are supersonic and the estimate is fair.
36+
37+geo = triangle_geometry(P);
38+nn = geo.nn;
39+nel = nn - 1;
40+E = P.youngs_modulus;
41+rho = P.density;
42+nu = P.poisson_ratio;
43+d = P.rod_diameter_m;
44+c = P.speed_of_sound;
45+rho0 = P.air_density;
46+nodes = geo.nodes;
47+
48+% ---- element geometry ----
49+ex = zeros(nel, 1);
50+ey = zeros(nel, 1);
51+tx = zeros(nel, 1);
52+ty = zeros(nel, 1);
53+el = zeros(nel, 1);
54+for e = 1:nel
55+ dx = nodes(1, e + 1) - nodes(1, e);
56+ dy = nodes(2, e + 1) - nodes(2, e);
57+ el(e) = sqrt(dx^2 + dy^2);
58+ tx(e) = dx / el(e);
59+ ty(e) = dy / el(e);
60+ ex(e) = (nodes(1, e) + nodes(1, e + 1)) / 2;
61+ ey(e) = (nodes(2, e) + nodes(2, e + 1)) / 2;
62+end
63+
64+% ---- eigensolve per family, collect retained modes ----
65+fs = [];
66+fams = [];
67+NUx = [];
68+NUy = [];
69+NUz = [];
70+EUpx = [];
71+EUpy = [];
72+EUpz = [];
73+omLo = 2 * pi * P.min_mode_hz;
74+for family = 0:1
75+ [K, M] = triangle_fem(nodes, family, E, rho, nu, d);
76+ Rc = chol(M);
77+ Am = (Rc') \ (K / Rc);
78+ Am = (Am + Am') / 2;
79+ [Uv, Sv, ~] = svd(Am);
80+ om2 = flipud(diag(Sv)); % ascending
81+ Vv = fliplr(Uv);
82+ om = sqrt(max(0, om2)) * P.stiffness_trim;
83+ keep = find(om > omLo & om < 2 * pi * P.max_mode_hz);
84+ if length(keep) == 0
85+ continue;
86+ end
87+ Phi = Rc \ Vv(:, keep); % mass-normalized shapes, 3 nn x nk
88+ nk = length(keep);
89+ if family == 0
90+ ux = Phi(1:3:3 * nn, :);
91+ uy = Phi(2:3:3 * nn, :);
92+ uz = zeros(nn, nk);
93+ eux = (ux(1:nn - 1, :) + ux(2:nn, :)) / 2;
94+ euy = (uy(1:nn - 1, :) + uy(2:nn, :)) / 2;
95+ ax = eux .* (tx * ones(1, nk)) + euy .* (ty * ones(1, nk));
96+ upx = eux - ax .* (tx * ones(1, nk));
97+ upy = euy - ax .* (ty * ones(1, nk));
98+ upz = zeros(nel, nk);
99+ else
100+ ux = zeros(nn, nk);
101+ uy = zeros(nn, nk);
102+ uz = Phi(1:3:3 * nn, :);
103+ upx = zeros(nel, nk);
104+ upy = zeros(nel, nk);
105+ upz = (uz(1:nn - 1, :) + uz(2:nn, :)) / 2;
106+ end
107+ fs = [fs; om(keep) / (2 * pi)];
108+ fams = [fams; family * ones(nk, 1)];
109+ NUx = [NUx, ux];
110+ NUy = [NUy, uy];
111+ NUz = [NUz, uz];
112+ EUpx = [EUpx, upx];
113+ EUpy = [EUpy, upy];
114+ EUpz = [EUpz, upz];
115+end
116+n = length(fs);
117+if n == 0
118+ error('No rod modes between min_mode_hz and max_mode_hz -- check geometry/material.');
119+end
120+
121+[fs, order] = sort(fs);
122+fs = reshape(fs, n, 1);
123+fams = reshape(fams(order), n, 1);
124+NUx = NUx(:, order);
125+NUy = NUy(:, order);
126+NUz = NUz(:, order);
127+EUpx = EUpx(:, order);
128+EUpy = EUpy(:, order);
129+EUpz = EUpz(:, order);
130+
131+% ---- damping: material + radiation ----
132+omega = 2 * pi * fs;
133+a = d / 2;
134+Vdip = 2 * pi * a^2 * el; % nel x 1 dipole volume per element
135+up2 = EUpx.^2 + EUpy.^2 + EUpz.^2; % nel x n
136+S2 = (Vdip.^2).' * up2; % 1 x n
137+zetaRad = rho0 * omega.^3 .* reshape(S2, n, 1) / (24 * pi * c^3);
138+zeta = 0.5 * (P.loss_factor + P.loss_slope * omega) + P.radiation_loss * zetaRad;
139+
140+basis.n = n;
141+basis.f = fs;
142+basis.omega = omega;
143+basis.zeta = zeta;
144+basis.zetaRad = zetaRad;
145+basis.family = fams;
146+basis.nodeX = nodes(1, :);
147+basis.nodeY = nodes(2, :);
148+basis.sfrac = geo.sfrac;
149+basis.totalLen = geo.totalLen;
150+basis.NUx = NUx;
151+basis.NUy = NUy;
152+basis.NUz = NUz;
153+basis.elemX = ex;
154+basis.elemY = ey;
155+basis.elemTanX = tx;
156+basis.elemTanY = ty;
157+basis.elemLen = el;
158+basis.EUpx = EUpx;
159+basis.EUpy = EUpy;
160+basis.EUpz = EUpz;
161+end
matlab/triangle_defaults.madded+60−0View file
@@ -0,0 +1,60 @@
1+function P = triangle_defaults()
2+% TRIANGLE_DEFAULTS Default parameters, mirroring the web app's params.ts.
3+% The output duration is not a parameter: it is
4+% (strike_count - 1) / strike_rate_hz + tail_s.
5+
6+% Rod geometry (equilateral, open at the bottom-left corner)
7+P.side_length_m = 0.152;
8+P.gap_m = 0.015;
9+P.rod_diameter_m = 0.008;
10+
11+% Rod material (steel)
12+P.density = 7850.0;
13+P.youngs_modulus = 2.0e11;
14+P.poisson_ratio = 0.30;
15+P.stiffness_trim = 1.00;
16+
17+% Discretization / mode retention
18+P.n_elements = 60;
19+P.min_mode_hz = 20.0;
20+P.max_mode_hz = 22000.0;
21+
22+% Damping
23+P.loss_factor = 6.0e-5;
24+P.loss_slope = 4.0e-10;
25+P.radiation_loss = 1.0;
26+
27+% Beater (steel rod tip)
28+P.beater_radius_m = 0.0045;
29+P.beater_mass_kg = 0.022;
30+P.beater_youngs = 2.0e11;
31+P.restitution = 0.60;
32+
33+% Strike / roll gesture
34+P.strike_velocity = 1.5;
35+P.strike_pos = 0.17; % arclength fraction: middle of the base
36+P.strike_angle_deg = 20.0; % 0 = in the plane of the triangle, 90 = normal to it
37+P.strike_count = 1;
38+P.strike_rate_hz = 12.0;
39+P.roll_pos2 = 0.40; % second contact point of a roll (the right leg)
40+P.timing_jitter_ms = 6.0;
41+P.velocity_spread = 0.15;
42+P.seed = 1;
43+
44+% Radiation / listener
45+P.direct_gain = 1.0;
46+P.listener_x_m = 0.20;
47+P.listener_y_m = 0.10;
48+P.listener_z_m = 0.40;
49+P.speed_of_sound = 343.0;
50+P.air_density = 1.204;
51+
52+% Rendering / output
53+P.sample_rate = 48000;
54+P.oversample = 2;
55+P.tail_s = 5.0;
56+P.pre_delay_ms = 3.0;
57+P.normalize = 1;
58+P.peak_dbfs = -3.0;
59+P.full_scale_pa = 2.0;
60+end
matlab/triangle_demo.madded+21−0View file
@@ -0,0 +1,21 @@
1+% Render a triangle strike with default parameters and print diagnostics.
2+% Run with: numbl run triangle_demo.m
3+
4+P = triangle_defaults();
5+tic;
6+[y, basis, dinfo] = triangle_render(P);
7+elapsed = toc;
8+
9+fprintf('modes kept : %d (%d in-plane, %d out-of-plane)\n', basis.n, dinfo.modesIn, dinfo.modesOut);
10+fprintf('mode freq : ');
11+fprintf('%.0f ', basis.f);
12+fprintf('Hz\n');
13+fprintf('strikes : %d (mean v %.2f m/s, mean contact %.0f us, peak force %.0f N)\n', ...
14+ dinfo.strikeCount, dinfo.meanV, dinfo.meanTauUs, dinfo.peakForceN);
15+fprintf('peak press.: %.4f Pa (%.1f dB SPL) at %.0f cm\n', ...
16+ dinfo.peakPa, dinfo.peakSplDb, dinfo.listenerDistM * 100);
17+fprintf('length : %.2f s @ %d Hz\n', length(y) / P.sample_rate, P.sample_rate);
18+fprintf('render time: %.2f s\n', elapsed);
19+if dinfo.listenerInside
20+ fprintf('WARNING : listener is at or inside the rod\n');
21+end
matlab/triangle_fem.madded+97−0View file
@@ -0,0 +1,97 @@
1+function [K, M] = triangle_fem(nodes, family, E, rho, nu, d)
2+% TRIANGLE_FEM Free-free frame system for one mode family of a bent rod.
3+%
4+% The rod centerline is the open polyline `nodes` (2 x nn, metres, in the
5+% x-y plane), discretized into nn - 1 two-node elements. Because the geometry
6+% is planar, the six degrees of freedom per node decouple exactly into two
7+% independent families of three, and each family is assembled separately:
8+% family 0 (in-plane): ux, uy, theta_z (in-plane bending + stretching)
9+% family 1 (out-of-plane): uz, theta_x, theta_y (out-of-plane bending + torsion)
10+% The bends are where the families earn their keep: within one straight leg
11+% out-of-plane bending and torsion are independent, but a change of direction
12+% rotates the (torsion, bending-rotation) pair, so the corner converts one
13+% into the other. That conversion is what the per-element rotation T applies,
14+% and it is what makes a triangle sound different from a straight rod.
15+%
16+% Elements are Euler-Bernoulli (cubic Hermite bending, linear axial and
17+% torsion interpolation) with consistent mass including rotary inertia.
18+% Shear deformation is neglected, which leaves the highest retained modes a
19+% few percent sharp for a rod this thick; see the README.
20+%
21+% Sign conventions: in-plane bending uses theta_z = +dv/dx, giving the
22+% textbook stiffness block; out-of-plane bending uses theta_y = -dw/dx, so
23+% the odd coupling terms flip sign (the similarity transform Sf below).
24+%
25+% Returns K and M, each 3 nn x 3 nn, symmetric; M is positive definite, so
26+% chol(M) reduces the generalized eigenproblem to a symmetric standard one.
27+
28+nn = size(nodes, 2);
29+ndof = 3 * nn;
30+K = zeros(ndof);
31+M = zeros(ndof);
32+
33+Aa = pi * d^2 / 4; % cross-section area
34+Iz = pi * d^4 / 64; % bending second moment (circular)
35+Jp = pi * d^4 / 32; % polar moment = torsion constant (circular)
36+Gs = E / (2 * (1 + nu));
37+
38+for e = 1:nn - 1
39+ dx = nodes(1, e + 1) - nodes(1, e);
40+ dy = nodes(2, e + 1) - nodes(2, e);
41+ l = sqrt(dx^2 + dy^2);
42+ c = dx / l;
43+ s = dy / l;
44+
45+ % bending block on the (transverse, rotation) pairs of both nodes
46+ Kb = (E * Iz / l^3) * ...
47+ [12, 6 * l, -12, 6 * l; ...
48+ 6 * l, 4 * l^2, -6 * l, 2 * l^2; ...
49+ -12, -6 * l, 12, -6 * l; ...
50+ 6 * l, 2 * l^2, -6 * l, 4 * l^2];
51+ Mb = (rho * Aa * l / 420) * ...
52+ [156, 22 * l, 54, -13 * l; ...
53+ 22 * l, 4 * l^2, 13 * l, -3 * l^2; ...
54+ 54, 13 * l, 156, -22 * l; ...
55+ -13 * l, -3 * l^2, -22 * l, 4 * l^2] ...
56+ + (rho * Iz / (30 * l)) * ...
57+ [36, 3 * l, -36, 3 * l; ...
58+ 3 * l, 4 * l^2, -3 * l, -l^2; ...
59+ -36, -3 * l, 36, -3 * l; ...
60+ 3 * l, -l^2, -3 * l, 4 * l^2];
61+
62+ Ke = zeros(6, 6);
63+ Me = zeros(6, 6);
64+ if family == 0
65+ % local order (u1, v1, t1, u2, v2, t2): axial on u, bending on (v, t)
66+ Ka2 = (E * Aa / l) * [1, -1; -1, 1];
67+ Ma2 = (rho * Aa * l / 6) * [2, 1; 1, 2];
68+ ia = [1 4];
69+ ib = [2 3 5 6];
70+ Ke(ia, ia) = Ka2;
71+ Me(ia, ia) = Ma2;
72+ Ke(ib, ib) = Kb;
73+ Me(ib, ib) = Mb;
74+ T = [c, s, 0; -s, c, 0; 0, 0, 1];
75+ else
76+ % local order (w1, tx1, ty1, w2, tx2, ty2): torsion on tx (about the leg
77+ % axis), bending on (w, ty) with theta_y = -dw/dx
78+ Kt2 = (Gs * Jp / l) * [1, -1; -1, 1];
79+ Mt2 = (rho * Jp * l / 6) * [2, 1; 1, 2];
80+ Sf = diag([1, -1, 1, -1]);
81+ ia = [2 5];
82+ ib = [1 3 4 6];
83+ Ke(ia, ia) = Kt2;
84+ Me(ia, ia) = Mt2;
85+ Ke(ib, ib) = Sf * Kb * Sf;
86+ Me(ib, ib) = Sf * Mb * Sf;
87+ T = [1, 0, 0; 0, c, s; 0, -s, c];
88+ end
89+
90+ TT = blkdiag(T, T);
91+ Kg = TT' * Ke * TT;
92+ Mg = TT' * Me * TT;
93+ dofs = [3 * e - 2, 3 * e - 1, 3 * e, 3 * e + 1, 3 * e + 2, 3 * e + 3];
94+ K(dofs, dofs) = K(dofs, dofs) + Kg;
95+ M(dofs, dofs) = M(dofs, dofs) + Mg;
96+end
97+end
matlab/triangle_geometry.madded+63−0View file
@@ -0,0 +1,63 @@
1+function G = triangle_geometry(P)
2+% TRIANGLE_GEOMETRY Node polyline of the open triangle centerline.
3+%
4+% An equilateral triangle of side side_length_m in the x-y plane, apex up:
5+% base from V2 (bottom left) to V3 (bottom right) along y = 0, apex V1 above
6+% the origin. The rod is open at the bottom-left corner: gap_m of arclength
7+% is removed there, half from the base and half from the left leg, so the
8+% centerline runs A -> V3 -> V1 -> B counterclockwise, where A and B sit
9+% gap_m / 2 from V2 along their legs. Total rod length is 3 L - gap.
10+%
11+% The polyline is discretized into about n_elements equal-length elements,
12+% with the two closed corners landing exactly on nodes. Arclength fraction 0
13+% is the free end on the base, 1 the free end on the left leg; strike_pos and
14+% roll_pos2 are expressed in this coordinate.
15+%
16+% Returns G with:
17+% nodes 2 x nn node coordinates (m)
18+% sfrac 1 x nn arclength fraction of each node
19+% totalLen rod arclength (m)
20+% nn, ne node count and per-leg element counts (1 x 3)
21+
22+L = P.side_length_m;
23+g = min(P.gap_m, 0.8 * L);
24+V1 = [0; L * sqrt(3) / 2];
25+V2 = [-L / 2; 0];
26+V3 = [L / 2; 0];
27+e1 = (V3 - V2) / L;
28+e3 = (V2 - V1) / L;
29+A = V2 + e1 * (g / 2);
30+B = V2 - e3 * (g / 2);
31+
32+starts = [A, V3, V1];
33+stops = [V3, V1, B];
34+lens = [L - g / 2, L, L - g / 2];
35+total = lens(1) + lens(2) + lens(3);
36+
37+ne = zeros(1, 3);
38+for k = 1:3
39+ ne(k) = max(2, round(P.n_elements * lens(k) / total));
40+end
41+nn = ne(1) + ne(2) + ne(3) + 1;
42+
43+nodes = zeros(2, nn);
44+sarc = zeros(1, nn);
45+nodes(:, 1) = A;
46+idx = 1;
47+acc = 0;
48+for k = 1:3
49+ for j = 1:ne(k)
50+ idx = idx + 1;
51+ frac = j / ne(k);
52+ nodes(:, idx) = starts(:, k) + (stops(:, k) - starts(:, k)) * frac;
53+ sarc(idx) = acc + lens(k) * frac;
54+ end
55+ acc = acc + lens(k);
56+end
57+
58+G.nodes = nodes;
59+G.sfrac = sarc / total;
60+G.totalLen = total;
61+G.nn = nn;
62+G.ne = ne;
63+end
matlab/triangle_output.madded+65−0View file
@@ -0,0 +1,65 @@
1+function [y, peakPa] = triangle_output(P, p, oversample)
2+% TRIANGLE_OUTPUT Anti-aliased decimation and output conditioning.
3+% Kaiser-windowed-sinc FIR lowpass at 0.45 * fs_lo, decimate by oversample,
4+% pre-roll, tail fade, then either normalize to peak_dbfs or divide by
5+% full_scale_pa. Identical in structure to the woodblock and maraca output
6+% stages.
7+
8+peakPa = max(abs(p));
9+R = oversample;
10+T = length(p);
11+
12+if R > 1
13+ taps = 12 * R + 1;
14+ M2 = (taps - 1) / 2;
15+ fcNorm = 0.45 / R; % cycles per input sample
16+ m = (0:taps - 1) - M2;
17+ h = sin(2 * pi * fcNorm * m) ./ (pi * m);
18+ h(M2 + 1) = 2 * fcNorm; % sinc value at m = 0
19+ r = m / M2;
20+ beta = 8;
21+ h = h .* (tr_i0(beta * sqrt(max(0, 1 - r.^2))) / tr_i0(beta));
22+ h = h / sum(h); % unity DC gain
23+ pf = conv(p, h);
24+ outLen = floor(T / R);
25+ y = reshape(pf((0:outLen - 1) * R + M2 + 1), 1, outLen); % centered pick
26+else
27+ y = p;
28+end
29+
30+% Pre-roll, then truncate to the nominal duration
31+padN = floor(P.pre_delay_ms * 1e-3 * P.sample_rate);
32+K = max(1, round(P.strike_count));
33+duration = (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
34+total = floor(duration * P.sample_rate);
35+if padN > 0
36+ y = [zeros(1, padN) y];
37+end
38+y = y(1:min(length(y), total));
39+
40+% Squared fade on the tail
41+fade = min(2048, floor(length(y) / 8));
42+if fade > 1
43+ y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;
44+end
45+
46+if P.normalize ~= 0
47+ pk = max(abs(y));
48+ if pk > 0
49+ y = y / pk * 10^(P.peak_dbfs / 20);
50+ end
51+else
52+ y = y / P.full_scale_pa;
53+end
54+end
55+
56+function v = tr_i0(x)
57+% Modified Bessel function of the first kind, order zero (series expansion).
58+v = ones(size(x));
59+term = ones(size(x));
60+hx = x / 2;
61+for k = 1:40
62+ term = term .* (hx / k).^2;
63+ v = v + term;
64+end
65+end
matlab/triangle_radiate.madded+181−0View file
@@ -0,0 +1,181 @@
1+function [p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, fsHi)
2+% TRIANGLE_RADIATE Modal synthesis and free-field pressure at the listener.
3+%
4+% Two contributions:
5+%
6+% 1. Modal ringing. A strike with impulse J at arclength position s gives
7+% every mode a modal-velocity jump J S(omega tau) (phi(s) . dhat), where
8+% S(x) = cos(x/2) / (1 - (x/pi)^2) is the spectrum of the half-sine
9+% contact pulse and dhat the strike direction (strike_angle_deg between
10+% the leg's in-plane normal and the plane normal, so one strike excites
11+% both mode families). The rod radiates through its transverse motion:
12+% each element is a point dipole of moment Vdip u_perp qdot (Vdip from
13+% triangle_basis.m), and per mode the element dipoles are summed
14+% coherently with their relative retardations,
15+% G = sum_e (i om rho0 / 4 pi) Vdip_e (u_e . rhat_e)
16+% (i k + 1/r_e) exp(-i k (r_e - rcen)) / r_e,
17+% evaluated at the mode's centre frequency (modes are narrow-band). The
18+% common carrier exp(-i k rcen) is stripped, so relative phase across the
19+% triangle is kept while the overall propagation delay is not double
20+% counted. This sum is where the interference across the frame lives:
21+% moving the listener re-runs only this stage.
22+%
23+% 2. Direct strike clicks. Modes above max_mode_hz are truncated; their
24+% collective off-resonance response is the immediate local motion of the
25+% rod around the contact. It is reinstated from the classic driving-point
26+% impedance of an infinite Euler beam, Z = 2 rho A c_B (1 + i) with
27+% c_B = sqrt(om kappa cL), evaluated at the pulse's own centre frequency
28+% om* = pi / tau: the contact force divided by |Z| is the local velocity,
29+% which radiates as a dipole of length half a flexural wavelength along
30+% the strike direction. Scaled by direct_gain. This is the least
31+% principled term of the model, and direct_gain = 0 turns it off.
32+%
33+% The synthesis is exact in time: between strikes every mode evolves as
34+% exp(lambda t), evaluated in chunks over each inter-strike segment (the
35+% chunking only bounds memory; nothing is integrated).
36+%
37+% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.
38+
39+c = P.speed_of_sound;
40+rho0 = P.air_density;
41+K = max(1, round(P.strike_count));
42+duration = (K - 1) / max(1e-3, P.strike_rate_hz) + P.tail_s;
43+T = floor(duration * fsHi);
44+N = basis.n;
45+I = length(tstr);
46+nel = length(basis.elemLen);
47+nn = length(basis.sfrac);
48+
49+lx = P.listener_x_m;
50+ly = P.listener_y_m;
51+lz = P.listener_z_m;
52+
53+% ---- element-to-listener geometry ----
54+rex = lx - basis.elemX; % nel x 1
55+rey = ly - basis.elemY;
56+rez = lz * ones(nel, 1);
57+re = sqrt(rex.^2 + rey.^2 + rez.^2);
58+listenerInside = min(re) < 0.02;
59+re = max(re, 0.01);
60+rhx = rex ./ re;
61+rhy = rey ./ re;
62+rhz = rez ./ re;
63+cx = sum(basis.elemX .* basis.elemLen) / sum(basis.elemLen);
64+cy = sum(basis.elemY .* basis.elemLen) / sum(basis.elemLen);
65+rcen = sqrt((lx - cx)^2 + (ly - cy)^2 + lz^2);
66+
67+% ---- per-mode radiation transfer and complex pole ----
68+a = P.rod_diameter_m / 2;
69+Vdip = 2 * pi * a^2 * basis.elemLen; % nel x 1
70+G = zeros(N, 1) + 0i;
71+lam = zeros(N, 1) + 0i;
72+phv = zeros(N, 1) + 0i; % modal-velocity phasor per unit jump
73+for m = 1:N
74+ om = basis.omega(m);
75+ zt = basis.zeta(m);
76+ k = om / c;
77+ proj = basis.EUpx(:, m) .* rhx + basis.EUpy(:, m) .* rhy + basis.EUpz(:, m) .* rhz;
78+ ph = exp(-1i * k * (re - rcen));
79+ G(m) = (1i * om * rho0 / (4 * pi)) * sum(Vdip .* proj .* (1i * k + 1 ./ re) .* ph ./ re);
80+ omd = om * sqrt(max(1e-12, 1 - zt^2));
81+ lam(m) = -zt * om + 1i * omd;
82+ phv(m) = 1 + 1i * zt * om / omd; % qdot(t) = Re[phv e^{lam t}] per unit jump
83+end
84+
85+% ---- per-strike modal weights: A(n, i) modal-velocity-jump phasors ----
86+A = zeros(N, I) + 0i;
87+dirx = zeros(1, I);
88+diry = zeros(1, I);
89+dirz = zeros(1, I);
90+hitx = zeros(1, I);
91+hity = zeros(1, I);
92+ang = P.strike_angle_deg * pi / 180;
93+for i = 1:I
94+ s = min(0.999, max(0.001, sstr(i)));
95+ idx = sum(basis.sfrac <= s);
96+ idx = min(max(idx, 1), nn - 1);
97+ w = (s - basis.sfrac(idx)) / (basis.sfrac(idx + 1) - basis.sfrac(idx));
98+ uxh = basis.NUx(idx, :) + w * (basis.NUx(idx + 1, :) - basis.NUx(idx, :));
99+ uyh = basis.NUy(idx, :) + w * (basis.NUy(idx + 1, :) - basis.NUy(idx, :));
100+ uzh = basis.NUz(idx, :) + w * (basis.NUz(idx + 1, :) - basis.NUz(idx, :));
101+ % strike direction: the leg's inward in-plane normal, tilted out of plane
102+ ex = basis.elemTanX(idx);
103+ ey = basis.elemTanY(idx);
104+ dirx(i) = cos(ang) * (-ey);
105+ diry(i) = cos(ang) * ex;
106+ dirz(i) = sin(ang);
107+ hitx(i) = basis.nodeX(idx) + w * (basis.nodeX(idx + 1) - basis.nodeX(idx));
108+ hity(i) = basis.nodeY(idx) + w * (basis.nodeY(idx + 1) - basis.nodeY(idx));
109+ proj = uxh * dirx(i) + uyh * diry(i) + uzh * dirz(i); % 1 x N
110+ x = basis.omega * taustr(i); % N x 1
111+ den = 1 - (x / pi).^2;
112+ S = cos(x / 2) ./ den;
113+ sing = abs(den) < 1e-6;
114+ S(sing) = pi / 4; % limit at x -> pi
115+ A(:, i) = phv .* (jstr(i) * S) .* reshape(proj, N, 1);
116+end
117+
118+% ---- segment synthesis: exact damped-exponential evolution, chunked ----
119+ks = min(T, max(1, round(tstr * fsHi) + 1));
120+p = zeros(1, T);
121+s = zeros(N, 1) + 0i;
122+Gt = G.';
123+dt = 1 / fsHi;
124+CH = 32768; % chunk length: bounds the exp matrix
125+prevK = ks(1);
126+i = 1;
127+while i <= I
128+ s = s + A(:, i);
129+ while i < I && ks(i + 1) == ks(i)
130+ i = i + 1;
131+ s = s + A(:, i);
132+ end
133+ if i < I
134+ nextK = ks(i + 1);
135+ else
136+ nextK = T + 1;
137+ end
138+ pos = prevK;
139+ while pos < nextK
140+ len = min(CH, nextK - pos);
141+ seg = real(Gt * ((s * ones(1, len)) .* exp(lam * ((0:len - 1) * dt))));
142+ p(pos:pos + len - 1) = p(pos:pos + len - 1) + seg;
143+ s = s .* exp(lam * (len * dt));
144+ pos = pos + len;
145+ end
146+ prevK = nextK;
147+ i = i + 1;
148+end
149+
150+% ---- direct strike clicks (driving-point velocity radiating as a dipole) ----
151+Aarea = pi * (2 * a)^2 / 4;
152+kapcL = (a / 2) * sqrt(P.youngs_modulus / P.density); % kappa * cL
153+directPeak = 0;
154+for i = 1:I
155+ rsx = lx - hitx(i);
156+ rsy = ly - hity(i);
157+ rs = max(0.01, sqrt(rsx^2 + rsy^2 + lz^2));
158+ cosg = (rsx * dirx(i) + rsy * diry(i) + lz * dirz(i)) / rs;
159+ nP = max(2, round(taustr(i) * fsHi));
160+ k0 = ks(i) + round((rs - rcen) / c * fsHi);
161+ if k0 < 1 || k0 + nP > T
162+ continue;
163+ end
164+ omStar = pi / taustr(i);
165+ cB = sqrt(omStar * kapcL);
166+ Zb = 2 * sqrt(2) * P.density * Aarea * cB; % |2 rho A cB (1 + i)|
167+ lamB = 2 * pi * cB / omStar; % flexural wavelength at om*
168+ mu = 2 * pi * a^2 * (lamB / 2) / Zb; % dipole moment per unit force
169+ amp = P.direct_gain * fstr(i) * mu * rho0 * omStar^2 * cosg / (4 * pi * rs * c);
170+ pulse = amp * sin(pi * (0:nP) / nP);
171+ p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;
172+ if abs(amp) > directPeak
173+ directPeak = abs(amp);
174+ end
175+end
176+
177+radinfo.gainAbs = abs(G);
178+radinfo.directPeakPa = directPeak;
179+radinfo.listenerInside = listenerInside;
180+radinfo.listenerDistM = rcen;
181+end
matlab/triangle_render.madded+27−0View file
@@ -0,0 +1,27 @@
1+function [y, basis, dinfo] = triangle_render(P)
2+% TRIANGLE_RENDER Full render: basis + strikes -> radiate -> output.
3+% Convenience composition of the stage functions; the web app calls the
4+% stages separately so it can cache intermediates (the eigensolve above all).
5+
6+basis = triangle_basis(P);
7+fsHi = P.sample_rate * P.oversample;
8+[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P);
9+[p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, fsHi);
10+[y, peakPa] = triangle_output(P, p, P.oversample);
11+
12+dinfo.strikeCount = sinfo.count;
13+dinfo.meanV = sinfo.meanV;
14+dinfo.maxV = sinfo.maxV;
15+dinfo.meanTauUs = sinfo.meanTauUs;
16+dinfo.peakForceN = sinfo.peakForceN;
17+dinfo.durationS = sinfo.durationS;
18+dinfo.modesIn = sum(basis.family == 0);
19+dinfo.modesOut = sum(basis.family == 1);
20+dinfo.f1 = basis.f(1);
21+dinfo.peakPa = peakPa;
22+dinfo.peakSplDb = 20 * log10(peakPa / 20e-6 + 1e-30);
23+dinfo.directPeakPa = radinfo.directPeakPa;
24+dinfo.listenerInside = radinfo.listenerInside;
25+dinfo.listenerDistM = radinfo.listenerDistM;
26+dinfo.sinfo = sinfo;
27+end
matlab/triangle_rng.madded+12−0View file
@@ -0,0 +1,12 @@
1+function [r, state] = triangle_rng(state, n)
2+% TRIANGLE_RNG Portable seeded uniform generator (32-bit LCG).
3+% Returns n uniforms in (0, 1) as a column vector, plus the advanced state.
4+% Deliberately not MATLAB's rand: the same seed must produce the identical
5+% roll pattern in the numbl CLI, Node, and the browser, and every double
6+% here stays below 2^53 so the arithmetic is exact everywhere.
7+r = zeros(n, 1);
8+for i = 1:n
9+ state = mod(state * 1664525 + 1013904223, 4294967296);
10+ r(i) = (state + 0.5) / 4294967296;
11+end
12+end
matlab/triangle_strikes.madded+114−0View file
@@ -0,0 +1,114 @@
1+function [tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P)
2+% TRIANGLE_STRIKES Beater strike events: a single strike, or a seeded roll.
3+%
4+% strike_count = 1 gives one deterministic strike of strike_velocity at
5+% strike_pos at t = 0; no random numbers are drawn. strike_count > 1 gives a
6+% roll: nominally periodic hits at strike_rate_hz, alternating between
7+% strike_pos and roll_pos2 (the beater rattling in a corner strikes the two
8+% legs that form it in turn), with Gaussian timing jitter (timing_jitter_ms)
9+% and lognormal velocity spread (velocity_spread). The first hit stays at
10+% t = 0 so the onset is fixed. All randomness comes from triangle_rng seeded
11+% with P.seed, so a given seed reproduces the identical roll in the CLI,
12+% Node, and the browser.
13+%
14+% Each strike is a Hertz contact of the beater tip (a steel sphere of radius
15+% beater_radius_m) on the rod:
16+% k_h = (4/3) E* sqrt(R_eff), R_eff = sqrt(r_b R''),
17+% 1/R'' = 1/r_b + 2/d (sphere on cylinder: R_eff is the
18+% geometric mean of the two relative
19+% principal curvatures)
20+% tau = 3.218 (m*^2 / (k_h^2 v))^{1/5} (Hertz contact time)
21+% J = (1 + restitution) m* v (momentum transferred)
22+% Fpk = pi J / (2 tau) (half-sine pulse peak)
23+% The reduced mass m* pairs the beater with the rod mass engaged during the
24+% contact. The engaged length is estimated from the flexural group velocity
25+% at the pulse's own bandwidth, l = 2 sqrt(pi kappa cL tau), which depends
26+% on tau, so a short fixed-point iteration closes the loop. This is an
27+% estimate, not a contact simulation: what it captures is that a heavy
28+% beater lengthens the contact and darkens the strike, while a light beater
29+% approaches the free-sphere limit and tau ~ v^(-1/5).
30+%
31+% Returns row vectors sorted by time: tstr (s), vstr (m/s), jstr (N s),
32+% taustr (s), fstr (peak N), sstr (arclength fraction of the contact), and
33+% sinfo with summary statistics.
34+
35+K = max(1, round(P.strike_count));
36+rate = max(1e-3, P.strike_rate_hz);
37+duration = (K - 1) / rate + P.tail_s;
38+
39+E = P.youngs_modulus;
40+rho = P.density;
41+d = P.rod_diameter_m;
42+nub = 0.30; % beater Poisson ratio, fixed
43+Estar = 1 / ((1 - nub^2) / P.beater_youngs + (1 - P.poisson_ratio^2) / E);
44+rb = P.beater_radius_m;
45+Rpp = 1 / (1 / rb + 2 / d);
46+Reff = sqrt(rb * Rpp);
47+kh = (4 / 3) * Estar * sqrt(Reff);
48+kapcL = (d / 4) * sqrt(E / rho); % kappa * cL: flexural dispersion constant
49+Aarea = pi * d^2 / 4;
50+
51+state = mod(floor(abs(P.seed)), 4294967296);
52+z1 = zeros(K, 1);
53+z2 = zeros(K, 1);
54+if K > 1
55+ [u, state] = triangle_rng(state, 4 * K);
56+ U = reshape(u, K, 4);
57+ z1 = sqrt(-2 * log(U(:, 1))) .* cos(2 * pi * U(:, 2));
58+ z2 = sqrt(-2 * log(U(:, 3))) .* cos(2 * pi * U(:, 4));
59+end
60+
61+tstr = zeros(1, K);
62+vstr = zeros(1, K);
63+jstr = zeros(1, K);
64+taustr = zeros(1, K);
65+fstr = zeros(1, K);
66+sstr = zeros(1, K);
67+for i = 1:K
68+ t = (i - 1) / rate;
69+ if i > 1
70+ t = max(0, t + P.timing_jitter_ms * 1e-3 * z1(i));
71+ end
72+ v = P.strike_velocity;
73+ if K > 1
74+ v = v * exp(P.velocity_spread * z2(i));
75+ end
76+ v = max(0.01, v);
77+
78+ % engaged-mass fixed point for the Hertz contact time
79+ tau = 30e-6;
80+ mstar = P.beater_mass_kg;
81+ for it = 1:4
82+ leng = 2 * sqrt(pi * kapcL * tau);
83+ mloc = rho * Aarea * leng;
84+ mstar = P.beater_mass_kg * mloc / (P.beater_mass_kg + mloc);
85+ tau = 3.218 * (mstar^2 / (kh^2 * v))^(1 / 5);
86+ end
87+
88+ tstr(i) = t;
89+ vstr(i) = v;
90+ jstr(i) = (1 + P.restitution) * mstar * v;
91+ taustr(i) = tau;
92+ fstr(i) = pi * jstr(i) / (2 * tau);
93+ if mod(i, 2) == 1
94+ sstr(i) = P.strike_pos;
95+ else
96+ sstr(i) = P.roll_pos2;
97+ end
98+end
99+
100+[tstr, order] = sort(tstr);
101+tstr = reshape(tstr, 1, K);
102+vstr = reshape(vstr(order), 1, K);
103+jstr = reshape(jstr(order), 1, K);
104+taustr = reshape(taustr(order), 1, K);
105+fstr = reshape(fstr(order), 1, K);
106+sstr = reshape(sstr(order), 1, K);
107+
108+sinfo.count = K;
109+sinfo.meanV = sum(vstr) / K;
110+sinfo.maxV = max(vstr);
111+sinfo.meanTauUs = sum(taustr) / K * 1e6;
112+sinfo.peakForceN = max(fstr);
113+sinfo.durationS = duration;
114+end
matlab/triangle_validate.madded+152−0View file
@@ -0,0 +1,152 @@
1+% Validate the triangle model. There is no external reference implementation,
2+% so the checks are (a) closed-form identities the physics must satisfy,
3+% (b) structural identities of the FEM eigensolve, (c) exact reproducibility
4+% of the seeded roll process, and (d) pinned end-to-end diagnostics from the
5+% default configuration, recorded when the model was built and treated as a
6+% regression reference since.
7+% Run with: numbl run triangle_validate.m
8+
9+fails = 0;
10+
11+% ---- 1. straight free-free rod against closed forms ----
12+% The same element routines that build the triangle, assembled on a straight
13+% polyline, must reproduce the textbook free-free bar: bending frequencies
14+% (beta L)^2 kappa cL / (2 pi L^2), the first axial mode cL / 2L, the first
15+% torsional mode cT / 2L, and exactly three rigid-body modes per family.
16+E = 2e11; rho = 7850; nu = 0.30; d = 0.008; L = 0.44;
17+nodes = [linspace(0, L, 61); zeros(1, 61)];
18+cL = sqrt(E / rho);
19+kap = d / 4;
20+blTab = [4.73004074 7.85320462 10.99560784];
21+fbClosed = blTab.^2 * kap * cL / (2 * pi * L^2);
22+faClosed = cL / (2 * L);
23+ftClosed = sqrt(E / (2 * (1 + nu)) / rho) / (2 * L);
24+f0 = tv_eigfreqs(nodes, 0, E, rho, nu, d);
25+f1 = tv_eigfreqs(nodes, 1, E, rho, nu, d);
26+% the rigid-body zeros come back as numerical noise well below the first
27+% elastic mode (185 Hz here); 5 Hz separates them cleanly
28+nRigid0 = sum(f0 < 5);
29+nRigid1 = sum(f1 < 5);
30+e0 = f0(f0 > 1);
31+e1 = f1(f1 > 1);
32+errB = max(abs(e0(1:3).' - fbClosed) ./ fbClosed);
33+% rotary inertia lowers the FEM bending values slightly below the closed
34+% form, which neglects it; 0.5% covers that plus discretization
35+errA = abs(min(abs(e0 - faClosed)) / faClosed);
36+errT = abs(min(abs(e1 - ftClosed)) / ftClosed);
37+ok = nRigid0 == 3 && nRigid1 == 3 && errB < 5e-3 && errA < 2e-3 && errT < 2e-3;
38+fails = fails + (~ok);
39+fprintf('straight rod: bend err %.1e, axial %.1e, torsion %.1e, rigid %d/%d %s\n', ...
40+ errB, errA, errT, nRigid0, nRigid1, tv_verdict(ok));
41+
42+% ---- 2. eigensolve identities on the triangle itself ----
43+% The Cholesky reduction must hand back mass-orthonormal shapes that
44+% diagonalize the stiffness: phi' M phi = I and phi' K phi = diag(om^2).
45+P = triangle_defaults();
46+geo = triangle_geometry(P);
47+[K, M] = triangle_fem(geo.nodes, 0, P.youngs_modulus, P.density, P.poisson_ratio, P.rod_diameter_m);
48+Rc = chol(M);
49+Am = (Rc') \ (K / Rc);
50+Am = (Am + Am') / 2;
51+[Uv, Sv, ~] = svd(Am);
52+om2 = flipud(diag(Sv));
53+Vv = fliplr(Uv);
54+Phi = Rc \ Vv;
55+eOrth = max(max(abs(Phi' * M * Phi - eye(length(om2)))));
56+KP = Phi' * K * Phi;
57+eDiag = max(max(abs(KP - diag(om2)))) / max(om2);
58+ok = eOrth < 1e-8 && eDiag < 1e-10;
59+fails = fails + (~ok);
60+fprintf('eigensolve: orthonormality %.1e, diagonalization %.1e %s\n', eOrth, eDiag, tv_verdict(ok));
61+
62+% ---- 3. Hertz contact-time scaling tau ~ v^(-1/5) ----
63+% With a beater much lighter than the engaged rod mass the reduced mass is a
64+% constant and the free-sphere scaling must emerge.
65+Q = P;
66+Q.beater_mass_kg = 5e-4;
67+[~, ~, ~, tau1, ~, ~, ~] = triangle_strikes(Q);
68+Q.strike_velocity = 2 * P.strike_velocity;
69+[~, ~, ~, tau2, ~, ~, ~] = triangle_strikes(Q);
70+ratio = tau2(1) / tau1(1);
71+ok = abs(ratio - 2^(-1 / 5)) < 0.03;
72+fails = fails + (~ok);
73+fprintf('hertz scaling: tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), tv_verdict(ok));
74+
75+% ---- 4. seeded roll is exactly reproducible ----
76+Q = P;
77+Q.strike_count = 12;
78+[t1r, v1r, j1r, ~, ~, s1r, g1] = triangle_strikes(Q);
79+[t2r, v2r, j2r, ~, ~, s2r, g2] = triangle_strikes(Q);
80+same = max(abs(t1r - t2r)) == 0 && max(abs(v1r - v2r)) == 0 && max(abs(s1r - s2r)) == 0;
81+ok = same ...
82+ && g1.count == 12 ...
83+ && abs(g1.meanV - 1.501211020) < 1e-6 ...
84+ && abs(sum(j1r) - 0.383069278968) < 1e-9 ...
85+ && abs(t1r(2) - 0.087732781124) < 1e-9;
86+fails = fails + (~ok);
87+fprintf('roll (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\n', ...
88+ Q.seed, g1.count, g1.meanV, sum(j1r), tv_verdict(ok));
89+
90+% ---- 5. far-field falloff: -6 dB per doubling of distance ----
91+% RMS over the render, with the two listeners on the same ray from the
92+% origin so every element dipole keeps its directivity factor: off the ray
93+% the in-plane dipoles fall faster than 1/r, which is physics rather than
94+% error. Modes are capped at 8 kHz so 2 m is beyond the Rayleigh distance
95+% of the frame; at 20 kHz it is not, and Fresnel-zone interference shifts
96+% the ratio by a few tenths of a dB, again physics rather than error.
97+Q = P;
98+Q.normalize = 0;
99+Q.tail_s = 1.0;
100+Q.max_mode_hz = 8000;
101+Q.listener_x_m = 0.6; Q.listener_y_m = 0.5; Q.listener_z_m = 1.8;
102+[y1, ~, ~] = triangle_render(Q);
103+Q.listener_x_m = 1.2; Q.listener_y_m = 1.0; Q.listener_z_m = 3.6;
104+[y2, ~, ~] = triangle_render(Q);
105+drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
106+ok = abs(drop + 6.02) < 0.2;
107+fails = fails + (~ok);
108+fprintf('falloff 2 m -> 4 m: %.2f dB rms (target -6.02) %s\n', drop, tv_verdict(ok));
109+
110+% ---- 6. pinned end-to-end diagnostics at the defaults ----
111+[y, basis, dinfo] = triangle_render(P);
112+ok = abs(dinfo.peakPa - 7.217688) < 0.02 ...
113+ && abs(basis.f(1) - 163.9200) < 0.1 ...
114+ && basis.n == 40;
115+fails = fails + (~ok);
116+fprintf('pinned render: peak %.4f Pa (target %.4f), mode1 %.2f Hz (target %.2f), %d modes %s\n', ...
117+ dinfo.peakPa, 7.217688, basis.f(1), 163.9200, basis.n, tv_verdict(ok));
118+
119+% ---- 7. the ring decays ----
120+% A triangle rings for seconds, so the tail is not silent; but the energy in
121+% the last half second must sit well below the first half second's.
122+fs = P.sample_rate;
123+n = length(y);
124+head = y(1:round(0.5 * fs));
125+tail = y(n - round(0.5 * fs):n);
126+dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
127+ok = dropE < -25;
128+fails = fails + (~ok);
129+fprintf('ring decay: last/first half second %.1f dB %s\n', dropE, tv_verdict(ok));
130+
131+if fails > 0
132+ fprintf('%d CHECK(S) FAILED\n', fails);
133+else
134+ fprintf('all checks passed\n');
135+end
136+
137+function f = tv_eigfreqs(nodes, family, E, rho, nu, d)
138+[K, M] = triangle_fem(nodes, family, E, rho, nu, d);
139+Rc = chol(M);
140+Am = (Rc') \ (K / Rc);
141+Am = (Am + Am') / 2;
142+[~, Sv, ~] = svd(Am);
143+f = sort(sqrt(diag(Sv))) / (2 * pi);
144+end
145+
146+function s = tv_verdict(ok)
147+if ok
148+ s = 'ok';
149+else
150+ s = 'MISMATCH';
151+end
152+end
webapp/.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+node_modules/
2+dist/
webapp/README.mdadded+46−0View file
@@ -0,0 +1,46 @@
1+# Triangle Simulator — web app
2+
3+Interactive browser implementation of the physical model described in `../README.md`. The algorithm itself is written in MATLAB syntax, in `../matlab/*.m`, and executed client-side by [numbl](https://numbl.org) inside a Web Worker; the TypeScript here is only the UI, the plots, the WAV encoder, and the marshalling of parameters and tensors in and out of the numbl runtime.
4+
5+The signal chain: a Hertz beater contact (or a seeded roll of them), the bending, stretching, and torsional modes of the open bent rod from a finite element eigensolve, a coherent sum of element dipoles to the listener point, and Kaiser-windowed anti-aliased decimation to 48 kHz.
6+
7+## Usage
8+
9+```
10+npm install
11+npm run dev # development server
12+npm run build # production build in dist/
13+npm run selftest # run ../matlab/triangle_validate.m through numbl in Node
14+```
15+
16+The self-test runs the closed-form and reproducibility checks and the pinned diagnostics. The MATLAB code can also be run directly: `cd ../matlab && numbl run triangle_demo.m` (or `triangle_validate.m`).
17+
18+## Why it renders fast
19+
20+Nothing is integrated. Between strikes every mode is a damped exponential, so the synthesis evaluates `exp(lambda t)` over each inter-strike segment in closed form (chunked only to bound memory) and adds each strike's modal velocity jump as it lands. The cost that remains is the eigensolve of the two finite element families, about a second, and that is exactly the piece the worker caches: striking somewhere else, changing the beater, or moving the microphone reuses the cached modes.
21+
22+## Reading and editing the code in the app
23+
24+The Code button opens the MATLAB source in a CodeMirror editor with MATLAB syntax highlighting. Every `.m` file in `../matlab` is listed, tagged by whether the app executes it or it is a CLI entry point. The sources are editable: an edit re-renders after the usual debounce, so changing the physics in the editor is heard and plotted the same way moving a slider is. Removing the bending correction sign flip in `triangle_fem.m`, for instance, decouples torsion from out-of-plane bending and audibly thins the cluster of partials.
25+
26+Edits live in memory only. They are deliberately not persisted, so a reload restores the code the app shipped with; "Revert file" and "Revert all" undo them sooner, and a dot marks the files that differ from the shipped source. A fault in edited code is reported with its file and line (numbl skips a file it cannot parse and then fails at the call site, so the worker captures the parse warning and reports that instead), the last successful render stays on screen, and reverting recovers.
27+
28+## How rendering is tiered
29+
30+Changing any parameter triggers a debounced fast preview (oversample 1); the Render button runs the full-quality pass. While the current result is a preview, the Play and Download buttons and the status badge say so. The worker runs each stage as a separate numbl execution and keys its caches by tier: rod geometry, material, damping, and discretization parameters invalidate the mode basis (the eigensolve); the strike and roll parameters, the beater, and the listener only re-run the strike process and the radiation; pre-delay, normalization, and bit depth only re-run the output stage. Intermediate tensors stay inside the worker as numbl RuntimeValues. Because the eigensolve does not depend on the synthesis rate, previews reuse the full-quality mode basis.
31+
32+The New seed button redraws the roll timing and velocities without touching any other parameter (a single strike is deterministic, so the button is disabled there). The Loop button repeats the render continuously, swapping in each new render as parameters change. Parameters persist in localStorage; "Reset to defaults" restores the values in `src/params.ts` (which `../matlab/triangle_defaults.m` mirrors).
33+
34+## Layout
35+
36+- `../matlab/` — the algorithm: `triangle_geometry.m`, `triangle_fem.m`, `triangle_basis.m`, `triangle_strikes.m`, `triangle_radiate.m`, `triangle_output.m`, `triangle_rng.m`, plus `triangle_render.m` / `triangle_demo.m` / `triangle_validate.m` for CLI use
37+- `src/params.ts` — parameter schema, defaults, tier assignment, persistence
38+- `src/mfiles.ts` — the `.m` sources, bundled at build time and hashed for cache keying
39+- `src/worker/renderWorker.ts` — off-thread numbl execution with tiered caching
40+- `src/lib/outline.ts` — the open-triangle centerline, mirrored from `triangle_geometry.m` for the scene and the mode-shape plot
41+- `src/components/` — 3D scene (three.js), parameter panel, code panel and editor, plots (strikes, waveform, spectrum, spectrogram, mode map, mode shapes)
42+- `src/selftest.ts` — Node harness for the MATLAB validation script
43+
44+## A note on the seeded generator
45+
46+The roll process draws from a 32-bit LCG in `triangle_rng.m` rather than MATLAB's `rand`, so a given seed produces the identical pattern in the numbl CLI, in Node, and in the browser. That is what lets the self-test pin exact impulse sums and hit times.
webapp/index.htmladded+13−0View file
@@ -0,0 +1,13 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
6+ <link rel="icon" href="data:image/svg+xml,<svg xmlns=%22http://www.w3.org/2000/svg%22 viewBox=%220 0 100 100%22><text y=%22.9em%22 font-size=%2290%22>🔺</text></svg>" />
7+ <title>Triangle Simulator</title>
8+ </head>
9+ <body>
10+ <div id="root"></div>
11+ <script type="module" src="/src/main.tsx"></script>
12+ </body>
13+</html>
webapp/package-lock.jsonadded+4230−0View file
This diff is 4,235 lines long and is not shown.
webapp/package.jsonadded+34−0View file
@@ -0,0 +1,34 @@
1+{
2+ "name": "triangle-webapp",
3+ "private": true,
4+ "version": "0.1.0",
5+ "type": "module",
6+ "scripts": {
7+ "dev": "vite",
8+ "build": "tsc --noEmit && vite build",
9+ "preview": "vite preview",
10+ "selftest": "tsx src/selftest.ts"
11+ },
12+ "dependencies": {
13+ "@codemirror/commands": "^6.11.0",
14+ "@codemirror/language": "^6.12.4",
15+ "@codemirror/legacy-modes": "^6.5.3",
16+ "@codemirror/state": "^6.7.1",
17+ "@codemirror/theme-one-dark": "^6.1.3",
18+ "@codemirror/view": "^6.43.9",
19+ "codemirror": "^6.0.2",
20+ "numbl": "^0.4.19",
21+ "react": "^18.3.1",
22+ "react-dom": "^18.3.1",
23+ "three": "^0.168.0"
24+ },
25+ "devDependencies": {
26+ "@types/react": "^18.3.5",
27+ "@types/react-dom": "^18.3.0",
28+ "@types/three": "^0.168.0",
29+ "@vitejs/plugin-react": "^4.3.1",
30+ "tsx": "^4.19.0",
31+ "typescript": "^5.5.4",
32+ "vite": "^5.4.3"
33+ }
34+}
webapp/src/App.tsxadded+280−0View file
@@ -0,0 +1,280 @@
1+import { useCallback, useEffect, useRef, useState } from 'react';
2+import { DEFAULTS, loadParams, saveParams, type TriangleParams } from './params';
3+import { DEFAULT_MFILES, type MFile } from './mfiles';
4+import { useRenderer } from './hooks/useRenderer';
5+import { encodeWav } from './lib/wav';
6+import ParamPanel from './components/ParamPanel';
7+import CodePanel from './components/CodePanel';
8+import Scene3D from './components/Scene3D';
9+import DiagnosticsPanel from './components/DiagnosticsPanel';
10+import WaveformPlot from './components/plots/WaveformPlot';
11+import SpectrumPlot from './components/plots/SpectrumPlot';
12+import SpectrogramPlot from './components/plots/SpectrogramPlot';
13+import StrikePlot from './components/plots/StrikePlot';
14+import ModeMap from './components/plots/ModeMap';
15+import ModeShapes from './components/plots/ModeShapes';
16+
17+export default function App() {
18+ const [params, setParams] = useState<TriangleParams>(() => loadParams());
19+ // Editable copy of the MATLAB implementation. Deliberately not persisted:
20+ // a reload restores the code the app shipped with.
21+ const [mfiles, setMfiles] = useState<MFile[]>(() => DEFAULT_MFILES.map((f) => ({ ...f })));
22+ const [codeOpen, setCodeOpen] = useState(false);
23+ const { result, error, busy, stage, progress, renderFull } = useRenderer(params, mfiles);
24+ const audioCtxRef = useRef<AudioContext | null>(null);
25+
26+ useEffect(() => saveParams(params), [params]);
27+
28+ const onChange = useCallback((key: keyof TriangleParams, value: number | boolean) => {
29+ setParams((p) => ({ ...p, [key]: value }));
30+ }, []);
31+
32+ const reset = useCallback(() => setParams({ ...DEFAULTS }), []);
33+
34+ // A new seed redraws the roll timing and velocities without touching
35+ // anything else. For a single strike nothing is random, so the button is
36+ // disabled there.
37+ const reseed = useCallback(
38+ () => setParams((p) => ({ ...p, seed: 1 + Math.floor(Math.random() * 9999) })),
39+ [],
40+ );
41+
42+ const onChangeFile = useCallback((name: string, source: string) => {
43+ setMfiles((fs) => fs.map((f) => (f.name === name ? { ...f, source } : f)));
44+ }, []);
45+
46+ const onResetFile = useCallback((name: string) => {
47+ setMfiles((fs) =>
48+ fs.map((f) =>
49+ f.name === name
50+ ? { ...f, source: DEFAULT_MFILES.find((d) => d.name === name)?.source ?? f.source }
51+ : f,
52+ ),
53+ );
54+ }, []);
55+
56+ const onResetAllFiles = useCallback(() => setMfiles(DEFAULT_MFILES.map((f) => ({ ...f }))), []);
57+ const codeEdited = mfiles.some(
58+ (f) => f.source !== DEFAULT_MFILES.find((d) => d.name === f.name)?.source,
59+ );
60+
61+ const [looping, setLooping] = useState(false);
62+ const loopSrcRef = useRef<AudioBufferSourceNode | null>(null);
63+
64+ const getAudioCtx = useCallback((sampleRate: number) => {
65+ let ctx = audioCtxRef.current;
66+ if (!ctx || ctx.state === 'closed') {
67+ ctx = new AudioContext({ sampleRate });
68+ audioCtxRef.current = ctx;
69+ }
70+ void ctx.resume();
71+ return ctx;
72+ }, []);
73+
74+ const makeBuffer = useCallback((ctx: AudioContext, r: NonNullable<typeof result>) => {
75+ const buf = ctx.createBuffer(1, r.y.length, r.sampleRate);
76+ buf.copyToChannel(new Float32Array(r.y), 0);
77+ return buf;
78+ }, []);
79+
80+ const play = useCallback(() => {
81+ if (!result) return;
82+ const ctx = getAudioCtx(result.sampleRate);
83+ const src = ctx.createBufferSource();
84+ src.buffer = makeBuffer(ctx, result);
85+ src.connect(ctx.destination);
86+ src.start();
87+ }, [result, getAudioCtx, makeBuffer]);
88+
89+ // Loop mode: repeat the strike or roll continuously, swapping in each new
90+ // render as parameters change so edits are audible without pressing Play.
91+ useEffect(() => {
92+ if (!looping || !result) return;
93+ const ctx = getAudioCtx(result.sampleRate);
94+ const src = ctx.createBufferSource();
95+ src.buffer = makeBuffer(ctx, result);
96+ src.loop = true;
97+ src.connect(ctx.destination);
98+ src.start();
99+ loopSrcRef.current = src;
100+ return () => {
101+ loopSrcRef.current = null;
102+ try {
103+ src.stop();
104+ } catch {
105+ // already stopped
106+ }
107+ src.disconnect();
108+ };
109+ }, [looping, result, getAudioCtx, makeBuffer]);
110+
111+ const download = useCallback(() => {
112+ if (!result) return;
113+ const blob = encodeWav(result.y, result.sampleRate, params.bit_depth === 16 ? 16 : 24);
114+ const url = URL.createObjectURL(blob);
115+ const a = document.createElement('a');
116+ a.href = url;
117+ a.download = result.quality === 'full' ? 'triangle.wav' : 'triangle_preview.wav';
118+ a.click();
119+ setTimeout(() => URL.revokeObjectURL(url), 5000);
120+ }, [result, params.bit_depth]);
121+
122+ return (
123+ <div className="app">
124+ <div className="topbar">
125+ <h1>Triangle Simulator</h1>
126+ <button className="primary" onClick={renderFull} disabled={busy}>
127+ Render
128+ </button>
129+ <button
130+ className={`play ${result?.quality === 'preview' ? 'preview' : ''}`}
131+ onClick={play}
132+ disabled={!result}
133+ title={
134+ result?.quality === 'preview'
135+ ? 'Playing the fast preview render; press Render for full quality'
136+ : 'Play the rendered strike'
137+ }
138+ >
139+ ▶ Play{result?.quality === 'preview' ? ' (preview)' : ''}
140+ </button>
141+ <button
142+ className={`loop ${looping ? 'active' : ''}`}
143+ onClick={() => setLooping((v) => !v)}
144+ disabled={!result && !looping}
145+ title="Repeat the strike continuously so parameter changes are audible immediately"
146+ >
147+ ⟳ Loop{looping ? ' on' : ''}
148+ </button>
149+ <button
150+ onClick={reseed}
151+ disabled={params.strike_count < 2}
152+ title="Draw a new roll pattern from a new random seed (rolls only)"
153+ >
154+ ⚄ New seed
155+ </button>
156+ <button onClick={download} disabled={!result}>
157+ Download .wav{result?.quality === 'preview' ? ' (preview)' : ''}
158+ </button>
159+ <button onClick={reset}>Reset to defaults</button>
160+ <button
161+ className={codeOpen ? 'code active' : 'code'}
162+ onClick={() => setCodeOpen((v) => !v)}
163+ title="View and edit the MATLAB source that numbl executes"
164+ >
165+ {'{ }'} Code{codeEdited ? ' •' : ''}
166+ </button>
167+ <div className="status">
168+ {error && <span className="badge err">code error</span>}
169+ {result && (
170+ <span className={`badge ${result.quality}`}>
171+ {result.quality === 'full' ? 'full quality' : 'preview quality — press Render'}
172+ </span>
173+ )}
174+ {busy ? (
175+ <>
176+ <span>{stage}…</span>
177+ <div className="progressbar">
178+ <div style={{ width: `${Math.round(progress * 100)}%` }} />
179+ </div>
180+ </>
181+ ) : (
182+ <span>ready</span>
183+ )}
184+ </div>
185+ </div>
186+ <div className="main">
187+ <div className="sidebar">
188+ <ParamPanel params={params} onChange={onChange} />
189+ </div>
190+ <div className="content">
191+ {error && (
192+ <div className="card wide error-card">
193+ <h2>MATLAB error</h2>
194+ <div>
195+ <strong>
196+ {error.file ? `${error.file}${error.line ? `:${error.line}` : ''}` : 'error'}
197+ </strong>{' '}
198+ {error.message}
199+ </div>
200+ <div className="code-hint">
201+ Showing the last successful render. Open Code to fix the source, or revert it.
202+ </div>
203+ </div>
204+ )}
205+ <div className="card wide">
206+ <h2>Geometry — rod, beater, microphone</h2>
207+ <Scene3D params={params} />
208+ </div>
209+ {result && (
210+ <>
211+ <div className="card">
212+ <h2>Diagnostics</h2>
213+ <DiagnosticsPanel result={result} />
214+ </div>
215+ <div className="card">
216+ <h2>Strikes</h2>
217+ <StrikePlot
218+ strikeT={result.strikeT}
219+ strikeV={result.strikeV}
220+ strikeS={result.strikeS}
221+ durationS={result.diag.durationS}
222+ />
223+ </div>
224+ <div className="card">
225+ <h2>Waveform</h2>
226+ <WaveformPlot y={result.y} sampleRate={result.sampleRate} />
227+ </div>
228+ <div className="card">
229+ <h2>Spectrum</h2>
230+ <SpectrumPlot
231+ y={result.y}
232+ sampleRate={result.sampleRate}
233+ modeF={result.modeF}
234+ modeFamily={result.modeFamily}
235+ />
236+ </div>
237+ <div className="card">
238+ <h2>Spectrogram</h2>
239+ <SpectrogramPlot y={result.y} sampleRate={result.sampleRate} />
240+ </div>
241+ <div className="card">
242+ <h2>Rod modes</h2>
243+ <ModeMap
244+ modeF={result.modeF}
245+ modeZeta={result.modeZeta}
246+ modeZetaRad={result.modeZetaRad}
247+ modeFamily={result.modeFamily}
248+ />
249+ </div>
250+ <div className="card wide">
251+ <h2>Mode shapes</h2>
252+ <ModeShapes
253+ nodeX={result.nodeX}
254+ nodeY={result.nodeY}
255+ shapeUx={result.shapeUx}
256+ shapeUy={result.shapeUy}
257+ shapeUz={result.shapeUz}
258+ modeF={result.modeF}
259+ modeFamily={result.modeFamily}
260+ />
261+ </div>
262+ </>
263+ )}
264+ </div>
265+ </div>
266+ {codeOpen && (
267+ <CodePanel
268+ mfiles={mfiles}
269+ onChangeFile={onChangeFile}
270+ onResetFile={onResetFile}
271+ onResetAll={onResetAllFiles}
272+ onClose={() => setCodeOpen(false)}
273+ onRenderFull={renderFull}
274+ error={error}
275+ busy={busy}
276+ />
277+ )}
278+ </div>
279+ );
280+}
webapp/src/components/CodeEditor.tsxadded+91−0View file
@@ -0,0 +1,91 @@
1+// CodeMirror 6 editor for one .m file. The parent gives each instance a key
2+// of the file name, so switching files remounts rather than syncing docs.
3+
4+import { useEffect, useRef } from 'react';
5+import { EditorState } from '@codemirror/state';
6+import { EditorView, keymap, lineNumbers, highlightActiveLine } from '@codemirror/view';
7+import { defaultKeymap, history, historyKeymap, indentWithTab } from '@codemirror/commands';
8+import { StreamLanguage, syntaxHighlighting, HighlightStyle } from '@codemirror/language';
9+import type { StreamParser } from '@codemirror/language';
10+import { octave } from '@codemirror/legacy-modes/mode/octave';
11+import { tags as t } from '@lezer/highlight';
12+
13+// The octave mode's keyword list is missing these MATLAB keywords, so promote
14+// them after it has tokenized the word.
15+const EXTRA_KEYWORDS = new Set(['function', 'break', 'parfor', 'spmd']);
16+
17+const matlabMode: StreamParser<unknown> = {
18+ ...octave,
19+ token(stream, state) {
20+ const start = stream.pos;
21+ const tok = octave.token(stream, state);
22+ if (EXTRA_KEYWORDS.has(stream.string.slice(start, stream.pos))) return 'keyword';
23+ return tok;
24+ },
25+};
26+
27+const matlabHighlight = HighlightStyle.define([
28+ { tag: t.keyword, color: '#e8a33d' },
29+ { tag: t.comment, color: '#6b7a86', fontStyle: 'italic' },
30+ { tag: t.string, color: '#a5d6a7' },
31+ { tag: t.number, color: '#4dd0e1' },
32+ { tag: t.operator, color: '#c9a0dc' },
33+ { tag: t.variableName, color: '#d6dae2' },
34+ { tag: t.standard(t.variableName), color: '#82b1ff' }, // builtins: zeros, sqrt, ...
35+ { tag: t.function(t.variableName), color: '#82b1ff' },
36+ { tag: t.atom, color: '#4dd0e1' },
37+ { tag: t.meta, color: '#c9a0dc' },
38+ { tag: t.invalid, color: '#e57373' },
39+]);
40+
41+const editorTheme = EditorView.theme(
42+ {
43+ '&': { color: '#d6dae2', backgroundColor: '#181b20', height: '100%' },
44+ '.cm-content': { fontFamily: 'ui-monospace, SFMono-Regular, Menlo, monospace', fontSize: '12.5px' },
45+ '.cm-gutters': { backgroundColor: '#181b20', color: '#4d5560', border: 'none' },
46+ '.cm-activeLine': { backgroundColor: '#20242b' },
47+ '.cm-cursor': { borderLeftColor: '#e8a33d' },
48+ '&.cm-focused': { outline: 'none' },
49+ '.cm-scroller': { overflow: 'auto' },
50+ '.cm-selectionBackground, ::selection': { backgroundColor: '#31465a' },
51+ '&.cm-focused .cm-selectionBackground': { backgroundColor: '#31465a' },
52+ },
53+ { dark: true },
54+);
55+
56+interface Props {
57+ source: string;
58+ onChange: (source: string) => void;
59+}
60+
61+export default function CodeEditor({ source, onChange }: Props) {
62+ const hostRef = useRef<HTMLDivElement | null>(null);
63+ const onChangeRef = useRef(onChange);
64+ onChangeRef.current = onChange;
65+
66+ useEffect(() => {
67+ const host = hostRef.current!;
68+ const view = new EditorView({
69+ parent: host,
70+ state: EditorState.create({
71+ doc: source,
72+ extensions: [
73+ lineNumbers(),
74+ highlightActiveLine(),
75+ history(),
76+ keymap.of([...defaultKeymap, ...historyKeymap, indentWithTab]),
77+ StreamLanguage.define(matlabMode),
78+ syntaxHighlighting(matlabHighlight),
79+ editorTheme,
80+ EditorView.updateListener.of((u) => {
81+ if (u.docChanged) onChangeRef.current(u.state.doc.toString());
82+ }),
83+ ],
84+ }),
85+ });
86+ return () => view.destroy();
87+ // eslint-disable-next-line react-hooks/exhaustive-deps
88+ }, []);
89+
90+ return <div className="code-editor" ref={hostRef} />;
91+}
webapp/src/components/CodePanel.tsxadded+109−0View file
@@ -0,0 +1,109 @@
1+// Overlay for reading and editing the MATLAB implementation. Edits take
2+// effect on the next debounced render and live only as long as the page.
3+
4+import { useState } from 'react';
5+import { DEFAULT_MFILES, EXECUTED_BY_APP, type MFile } from '../mfiles';
6+import type { RenderError } from '../hooks/useRenderer';
7+import CodeEditor from './CodeEditor';
8+
9+interface Props {
10+ mfiles: MFile[];
11+ onChangeFile: (name: string, source: string) => void;
12+ onResetFile: (name: string) => void;
13+ onResetAll: () => void;
14+ onClose: () => void;
15+ onRenderFull: () => void;
16+ error: RenderError | null;
17+ busy: boolean;
18+}
19+
20+function defaultSourceOf(name: string): string {
21+ return DEFAULT_MFILES.find((f) => f.name === name)?.source ?? '';
22+}
23+
24+export default function CodePanel({
25+ mfiles,
26+ onChangeFile,
27+ onResetFile,
28+ onResetAll,
29+ onClose,
30+ onRenderFull,
31+ error,
32+ busy,
33+}: Props) {
34+ const [selected, setSelected] = useState(
35+ mfiles.find((f) => EXECUTED_BY_APP.has(f.name))?.name ?? mfiles[0]?.name ?? '',
36+ );
37+ const file = mfiles.find((f) => f.name === selected) ?? mfiles[0];
38+ const dirty = (f: MFile) => f.source !== defaultSourceOf(f.name);
39+ const anyDirty = mfiles.some(dirty);
40+
41+ return (
42+ <div className="modal-backdrop" onClick={onClose}>
43+ <div className="code-panel" onClick={(e) => e.stopPropagation()}>
44+ <div className="code-header">
45+ <h2>MATLAB source</h2>
46+ <span className="code-note">
47+ Executed by numbl. Edits re-render as you type and are discarded on reload.
48+ </span>
49+ <button onClick={onRenderFull} disabled={busy}>
50+ Render full quality
51+ </button>
52+ <button onClick={onResetAll} disabled={!anyDirty}>
53+ Revert all
54+ </button>
55+ <button onClick={onClose}>Close</button>
56+ </div>
57+ <div className="code-body">
58+ <div className="code-files">
59+ {mfiles.map((f) => (
60+ <button
61+ key={f.name}
62+ className={f.name === file?.name ? 'active' : ''}
63+ onClick={() => setSelected(f.name)}
64+ >
65+ <span className="fname">
66+ {f.name}
67+ {dirty(f) ? <span className="dot" title="edited" /> : null}
68+ </span>
69+ <span className="ftag">{EXECUTED_BY_APP.has(f.name) ? 'in app' : 'CLI'}</span>
70+ </button>
71+ ))}
72+ </div>
73+ <div className="code-main">
74+ {file && (
75+ <>
76+ <div className="code-toolbar">
77+ <span>{file.name}</span>
78+ {!EXECUTED_BY_APP.has(file.name) && (
79+ <span className="code-hint">
80+ entry point for the numbl CLI; the app does not run this file
81+ </span>
82+ )}
83+ <button
84+ onClick={() => onResetFile(file.name)}
85+ disabled={!dirty(file)}
86+ style={{ marginLeft: 'auto' }}
87+ >
88+ Revert file
89+ </button>
90+ </div>
91+ <CodeEditor
92+ key={file.name}
93+ source={file.source}
94+ onChange={(src) => onChangeFile(file.name, src)}
95+ />
96+ </>
97+ )}
98+ </div>
99+ </div>
100+ {error && (
101+ <div className="code-error">
102+ <strong>{error.file ? `${error.file}${error.line ? `:${error.line}` : ''}` : 'error'}</strong>{' '}
103+ {error.message}
104+ </div>
105+ )}
106+ </div>
107+ </div>
108+ );
109+}
webapp/src/components/DiagnosticsPanel.tsxadded+65−0View file
@@ -0,0 +1,65 @@
1+import type { RenderResponse } from '../worker/renderWorker';
2+
3+export default function DiagnosticsPanel({ result }: { result: RenderResponse }) {
4+ const d = result.diag;
5+ const modes = Array.from(result.modeF);
6+ const shown = modes.slice(0, 10);
7+ return (
8+ <div className="diag">
9+ <table>
10+ <tbody>
11+ <tr>
12+ <th>Strikes</th>
13+ <td>
14+ {d.strikeCount === 1
15+ ? `single, ${d.meanV.toFixed(2)} m/s`
16+ : `${d.strikeCount} · mean ${d.meanV.toFixed(2)} m/s, max ${d.maxV.toFixed(2)} m/s`}
17+ </td>
18+ </tr>
19+ <tr>
20+ <th>Contact</th>
21+ <td>
22+ {d.meanTauUs.toFixed(0)} µs mean (Hertz) · peak force {d.peakForceN.toFixed(0)} N
23+ </td>
24+ </tr>
25+ <tr>
26+ <th>Rod modes</th>
27+ <td className="partials">
28+ {shown.map((f) => f.toFixed(0)).join(', ')}
29+ {modes.length > shown.length ? `, … (${modes.length} total)` : ''} Hz
30+ </td>
31+ </tr>
32+ <tr>
33+ <th>Families</th>
34+ <td>
35+ {d.modesIn} in-plane, {d.modesOut} out-of-plane
36+ </td>
37+ </tr>
38+ <tr>
39+ <th>Listener</th>
40+ <td>{(d.listenerDistM * 100).toFixed(0)} cm from the frame</td>
41+ </tr>
42+ <tr>
43+ <th>Peak pressure</th>
44+ <td>
45+ {d.peakPa.toFixed(3)} Pa ({d.peakSplDb.toFixed(1)} dB SPL)
46+ </td>
47+ </tr>
48+ <tr>
49+ <th>Length</th>
50+ <td>
51+ {d.durationS.toFixed(2)} s · {result.quality} · oversample {d.oversample}
52+ </td>
53+ </tr>
54+ </tbody>
55+ </table>
56+ <div className="cachenote">
57+ eigensolve: {d.solveCached ? 'cached' : `${(d.solveMs / 1000).toFixed(2)} s`} · radiation:{' '}
58+ {d.radiateCached ? 'cached' : `${(d.radiateMs / 1000).toFixed(2)} s`}
59+ </div>
60+ {d.warnings.map((w, i) => (
61+ <div key={i} className="warning">⚠ {w}</div>
62+ ))}
63+ </div>
64+ );
65+}
webapp/src/components/ParamPanel.tsxadded+123−0View file
@@ -0,0 +1,123 @@
1+import { PARAM_GROUPS, type ParamSpec, type TriangleParams } from '../params';
2+
3+interface Props {
4+ params: TriangleParams;
5+ onChange: (key: keyof TriangleParams, value: number | boolean) => void;
6+}
7+
8+const SLIDER_STEPS = 1000;
9+
10+function sliderFromValue(spec: ParamSpec, v: number): number {
11+ if (spec.log) {
12+ const t = (Math.log(v) - Math.log(spec.min)) / (Math.log(spec.max) - Math.log(spec.min));
13+ return Math.round(t * SLIDER_STEPS);
14+ }
15+ return Math.round(((v - spec.min) / (spec.max - spec.min)) * SLIDER_STEPS);
16+}
17+
18+function valueFromSlider(spec: ParamSpec, s: number): number {
19+ const t = s / SLIDER_STEPS;
20+ let v: number;
21+ if (spec.log) {
22+ v = Math.exp(Math.log(spec.min) + t * (Math.log(spec.max) - Math.log(spec.min)));
23+ // keep three significant digits so the displayed number stays readable
24+ const mag = 10 ** Math.floor(Math.log10(Math.abs(v)));
25+ v = Math.round((v / mag) * 100) / 100 * mag;
26+ } else {
27+ v = spec.min + t * (spec.max - spec.min);
28+ if (spec.step) v = Math.round(v / spec.step) * spec.step;
29+ }
30+ if (spec.integer) v = Math.round(v);
31+ return Math.min(spec.max, Math.max(spec.min, v));
32+}
33+
34+function formatValue(v: number): string {
35+ if (v === 0) return '0';
36+ const a = Math.abs(v);
37+ if (a >= 1e5 || a < 1e-3) return v.toExponential(2);
38+ return String(Math.round(v * 1e6) / 1e6);
39+}
40+
41+function ParamRow({ spec, params, onChange }: Props & { spec: ParamSpec }) {
42+ const v = params[spec.key] as number;
43+ return (
44+ <div className="param-row" title={spec.help ?? ''}>
45+ <div className="param-label">
46+ <span className="name">
47+ {spec.label}
48+ {spec.unit ? <span style={{ color: 'var(--muted)' }}> ({spec.unit})</span> : null}
49+ </span>
50+ </div>
51+ <div className="param-controls">
52+ <input
53+ type="range"
54+ min={0}
55+ max={SLIDER_STEPS}
56+ value={sliderFromValue(spec, v)}
57+ onChange={(e) => onChange(spec.key, valueFromSlider(spec, Number(e.target.value)))}
58+ />
59+ <input
60+ type="number"
61+ value={formatValue(v)}
62+ step={spec.step ?? 'any'}
63+ onChange={(e) => {
64+ const nv = Number(e.target.value);
65+ if (Number.isFinite(nv)) onChange(spec.key, spec.integer ? Math.round(nv) : nv);
66+ }}
67+ />
68+ </div>
69+ </div>
70+ );
71+}
72+
73+export default function ParamPanel({ params, onChange }: Props) {
74+ return (
75+ <div>
76+ {PARAM_GROUPS.map((g, gi) => (
77+ <details key={g.name} className="group" open={gi < 2}>
78+ <summary>{g.name}</summary>
79+ <div className="group-body">
80+ {g.specs.map((spec) => (
81+ <ParamRow key={spec.key} spec={spec} params={params} onChange={onChange} />
82+ ))}
83+ {g.name === 'Output' && (
84+ <>
85+ <div className="param-row">
86+ <div className="param-label">
87+ <span className="name">Normalize peak</span>
88+ </div>
89+ <div className="param-controls">
90+ <input
91+ type="checkbox"
92+ checked={params.normalize}
93+ onChange={(e) => onChange('normalize', e.target.checked)}
94+ />
95+ <span className="param-help">
96+ {params.normalize
97+ ? `scaled to ${params.peak_dbfs} dBFS`
98+ : `absolute: ${params.full_scale_pa} Pa = 0 dBFS`}
99+ </span>
100+ </div>
101+ </div>
102+ <div className="param-row">
103+ <div className="param-label">
104+ <span className="name">WAV bit depth</span>
105+ </div>
106+ <div className="param-controls">
107+ <select
108+ value={params.bit_depth}
109+ onChange={(e) => onChange('bit_depth', Number(e.target.value))}
110+ >
111+ <option value={16}>16-bit</option>
112+ <option value={24}>24-bit</option>
113+ </select>
114+ </div>
115+ </div>
116+ </>
117+ )}
118+ </div>
119+ </details>
120+ ))}
121+ </div>
122+ );
123+}
webapp/src/components/Scene3D.tsxadded+256−0View file
@@ -0,0 +1,256 @@
1+// 3D view of the geometry: the open triangle hanging in a vertical plane,
2+// the beater at the strike point (and the second roll contact), and the
3+// listener. Model coordinates put the triangle in the x-y plane with the
4+// apex up and z normal to the plane; three.js is y-up too, so the mapping
5+// is the identity.
6+
7+import { useEffect, useRef } from 'react';
8+import * as THREE from 'three';
9+import { OrbitControls } from 'three/examples/jsm/controls/OrbitControls.js';
10+import type { TriangleParams } from '../params';
11+import { triangleOutline } from '../lib/outline';
12+
13+function makeLabel(text: string, color: string, scale = 1): THREE.Sprite {
14+ const canvas = document.createElement('canvas');
15+ canvas.width = 320;
16+ canvas.height = 64;
17+ const ctx = canvas.getContext('2d')!;
18+ ctx.font = '600 36px system-ui, sans-serif';
19+ ctx.fillStyle = color;
20+ ctx.textAlign = 'center';
21+ ctx.textBaseline = 'middle';
22+ ctx.fillText(text, 160, 32);
23+ const tex = new THREE.CanvasTexture(canvas);
24+ const mat = new THREE.SpriteMaterial({ map: tex, depthTest: false, transparent: true });
25+ const sprite = new THREE.Sprite(mat);
26+ sprite.scale.set(0.075 * scale, 0.015 * scale, 1);
27+ return sprite;
28+}
29+
30+interface Props {
31+ params: TriangleParams;
32+}
33+
34+export default function Scene3D({ params }: Props) {
35+ const mountRef = useRef<HTMLDivElement | null>(null);
36+ const sceneRef = useRef<{
37+ scene: THREE.Scene;
38+ group: THREE.Group | null;
39+ renderer: THREE.WebGLRenderer;
40+ camera: THREE.PerspectiveCamera;
41+ } | null>(null);
42+
43+ // one-time scene setup
44+ useEffect(() => {
45+ const mount = mountRef.current!;
46+ const renderer = new THREE.WebGLRenderer({ antialias: true });
47+ renderer.setPixelRatio(window.devicePixelRatio);
48+ renderer.setClearColor(0x14161a);
49+ mount.appendChild(renderer.domElement);
50+
51+ const scene = new THREE.Scene();
52+ const camera = new THREE.PerspectiveCamera(45, 1, 0.005, 20);
53+ camera.position.set(0.16, 0.16, 0.34);
54+
55+ const controls = new OrbitControls(camera, renderer.domElement);
56+ controls.target.set(0, 0.07, 0);
57+ controls.enableDamping = true;
58+
59+ scene.add(new THREE.AmbientLight(0xffffff, 0.55));
60+ const dir = new THREE.DirectionalLight(0xffffff, 1.1);
61+ dir.position.set(0.5, 1.0, 0.6);
62+ scene.add(dir);
63+ const dir2 = new THREE.DirectionalLight(0x8899ff, 0.3);
64+ dir2.position.set(-0.6, 0.3, -0.5);
65+ scene.add(dir2);
66+
67+ const grid = new THREE.GridHelper(0.8, 32, 0x3a4048, 0x262b33);
68+ grid.position.y = -0.06;
69+ scene.add(grid);
70+
71+ // model axes: x red (along the base), y green (up), z blue (plane normal)
72+ const axisLen = 0.06;
73+ const axes: [THREE.Vector3, number, string][] = [
74+ [new THREE.Vector3(axisLen, 0, 0), 0xe57373, '+x'],
75+ [new THREE.Vector3(0, axisLen, 0), 0x81c784, '+y'],
76+ [new THREE.Vector3(0, 0, axisLen), 0x64b5f6, '+z'],
77+ ];
78+ for (const [v, color, label] of axes) {
79+ const geo = new THREE.BufferGeometry().setFromPoints([new THREE.Vector3(0, 0, 0), v]);
80+ scene.add(new THREE.Line(geo, new THREE.LineBasicMaterial({ color })));
81+ const s = makeLabel(label, '#' + color.toString(16).padStart(6, '0'));
82+ s.position.copy(v).multiplyScalar(1.15);
83+ scene.add(s);
84+ }
85+
86+ sceneRef.current = { scene, group: null, renderer, camera };
87+
88+ let raf = 0;
89+ const animate = () => {
90+ raf = requestAnimationFrame(animate);
91+ controls.update();
92+ renderer.render(scene, camera);
93+ };
94+ animate();
95+
96+ const resize = () => {
97+ const w = mount.clientWidth;
98+ const h = mount.clientHeight;
99+ renderer.setSize(w, h, false);
100+ camera.aspect = w / h;
101+ camera.updateProjectionMatrix();
102+ };
103+ resize();
104+ const ro = new ResizeObserver(resize);
105+ ro.observe(mount);
106+
107+ return () => {
108+ cancelAnimationFrame(raf);
109+ ro.disconnect();
110+ controls.dispose();
111+ renderer.dispose();
112+ mount.removeChild(renderer.domElement);
113+ sceneRef.current = null;
114+ };
115+ }, []);
116+
117+ // rebuild the geometry group whenever relevant parameters change
118+ useEffect(() => {
119+ const ctx = sceneRef.current;
120+ if (!ctx) return;
121+ if (ctx.group) {
122+ ctx.scene.remove(ctx.group);
123+ ctx.group.traverse((o) => {
124+ const mesh = o as THREE.Mesh;
125+ if (mesh.geometry) mesh.geometry.dispose();
126+ });
127+ }
128+
129+ const g = new THREE.Group();
130+ const outline = triangleOutline(params.side_length_m, params.gap_m);
131+ const a = params.rod_diameter_m / 2;
132+ const steel = new THREE.MeshStandardMaterial({
133+ color: 0xb8bec8,
134+ roughness: 0.35,
135+ metalness: 0.85,
136+ });
137+
138+ // the rod: one cylinder per leg, spheres closing the two corners and
139+ // capping the free ends at the gap
140+ for (let k = 0; k < 3; k++) {
141+ const p0 = outline.pts[k];
142+ const p1 = outline.pts[k + 1];
143+ const len = Math.hypot(p1.x - p0.x, p1.y - p0.y);
144+ const cyl = new THREE.Mesh(new THREE.CylinderGeometry(a, a, len, 20), steel);
145+ cyl.position.set((p0.x + p1.x) / 2, (p0.y + p1.y) / 2, 0);
146+ const dirv = new THREE.Vector3(p1.x - p0.x, p1.y - p0.y, 0).normalize();
147+ cyl.quaternion.setFromUnitVectors(new THREE.Vector3(0, 1, 0), dirv);
148+ g.add(cyl);
149+ }
150+ for (const p of outline.pts) {
151+ const cap = new THREE.Mesh(new THREE.SphereGeometry(a, 14, 10), steel);
152+ cap.position.set(p.x, p.y, 0);
153+ g.add(cap);
154+ }
155+
156+ // suspension loop at the apex
157+ const apexY = (params.side_length_m * Math.sqrt(3)) / 2;
158+ const loop = new THREE.Mesh(
159+ new THREE.TorusGeometry(0.011, 0.0016, 8, 24),
160+ new THREE.MeshStandardMaterial({ color: 0x8d6e63, roughness: 0.8 }),
161+ );
162+ loop.position.set(0, apexY + 0.013, 0);
163+ g.add(loop);
164+ const cordGeo = new THREE.BufferGeometry().setFromPoints([
165+ new THREE.Vector3(0, apexY + 0.022, 0),
166+ new THREE.Vector3(0, apexY + 0.08, 0),
167+ ]);
168+ g.add(new THREE.Line(cordGeo, new THREE.LineBasicMaterial({ color: 0x6b5a50 })));
169+
170+ // beater at the strike point, pointing along the strike direction
171+ const addBeater = (frac: number, primary: boolean) => {
172+ const pt = outline.at(frac);
173+ const ang = (params.strike_angle_deg * Math.PI) / 180;
174+ // inward in-plane normal of a counterclockwise outline, tilted out of plane
175+ const dirv = new THREE.Vector3(
176+ Math.cos(ang) * -pt.ty,
177+ Math.cos(ang) * pt.tx,
178+ Math.sin(ang),
179+ );
180+ const tip = new THREE.Vector3(pt.x, pt.y, 0).addScaledVector(
181+ dirv,
182+ -(a + params.beater_radius_m + 0.001),
183+ );
184+ const color = primary ? 0xe8a33d : 0x9a7a3d;
185+ const head = new THREE.Mesh(
186+ new THREE.SphereGeometry(params.beater_radius_m, 16, 12),
187+ new THREE.MeshStandardMaterial({ color, roughness: 0.4, metalness: 0.6 }),
188+ );
189+ head.position.copy(tip);
190+ g.add(head);
191+ const shaftLen = 0.07;
192+ const shaft = new THREE.Mesh(
193+ new THREE.CylinderGeometry(0.0022, 0.0022, shaftLen, 10),
194+ new THREE.MeshStandardMaterial({ color: 0x777d88, roughness: 0.5, metalness: 0.7 }),
195+ );
196+ shaft.quaternion.setFromUnitVectors(new THREE.Vector3(0, 1, 0), dirv.clone().negate());
197+ shaft.position.copy(tip).addScaledVector(dirv, -(params.beater_radius_m + shaftLen / 2));
198+ g.add(shaft);
199+ if (primary) {
200+ const lbl = makeLabel('strike', '#e8a33d');
201+ lbl.position.copy(tip).add(new THREE.Vector3(0, -0.02, 0.01));
202+ g.add(lbl);
203+ }
204+ };
205+ addBeater(params.strike_pos, true);
206+ if (params.strike_count > 1) addBeater(params.roll_pos2, false);
207+
208+ // listener: microphone marker + line from the frame centroid
209+ const lis = new THREE.Vector3(params.listener_x_m, params.listener_y_m, params.listener_z_m);
210+ const micHead = new THREE.Mesh(
211+ new THREE.SphereGeometry(0.016, 20, 16),
212+ new THREE.MeshStandardMaterial({ color: 0x4dd0e1, emissive: 0x0c3b42, roughness: 0.35 }),
213+ );
214+ micHead.position.copy(lis);
215+ g.add(micHead);
216+ const micBody = new THREE.Mesh(
217+ new THREE.CylinderGeometry(0.008, 0.008, 0.05, 12),
218+ new THREE.MeshStandardMaterial({ color: 0x455a64, roughness: 0.6 }),
219+ );
220+ micBody.position.copy(lis).add(new THREE.Vector3(0, -0.033, 0));
221+ g.add(micBody);
222+ const centroid = new THREE.Vector3(0, apexY / 2.2, 0);
223+ const lineGeo = new THREE.BufferGeometry().setFromPoints([centroid, lis]);
224+ const line = new THREE.Line(
225+ lineGeo,
226+ new THREE.LineDashedMaterial({ color: 0x4dd0e1, dashSize: 0.012, gapSize: 0.009 }),
227+ );
228+ line.computeLineDistances();
229+ g.add(line);
230+ const dist = centroid.distanceTo(lis);
231+ const micLabel = makeLabel(`mic ${(dist * 100).toFixed(0)} cm`, '#4dd0e1');
232+ micLabel.position.copy(lis).add(new THREE.Vector3(0, 0.045, 0));
233+ g.add(micLabel);
234+
235+ ctx.scene.add(g);
236+ ctx.group = g;
237+ }, [
238+ params.side_length_m,
239+ params.gap_m,
240+ params.rod_diameter_m,
241+ params.strike_pos,
242+ params.strike_angle_deg,
243+ params.strike_count,
244+ params.roll_pos2,
245+ params.beater_radius_m,
246+ params.listener_x_m,
247+ params.listener_y_m,
248+ params.listener_z_m,
249+ ]);
250+
251+ return (
252+ <div className="scene3d" ref={mountRef}>
253+ <div className="scene-hint">drag to orbit · scroll to zoom · right-drag to pan</div>
254+ </div>
255+ );
256+}
webapp/src/components/plots/ModeMap.tsxadded+97−0View file
@@ -0,0 +1,97 @@
1+// The retained rod modes: frequency against ring time T60 = 6.91 / (zeta om),
2+// colored by family. In-plane and out-of-plane modes interleave across the
3+// whole band; the marker area shows how much of each mode's damping is
4+// radiation (a large marker is a mode the air takes energy from, a small one
5+// rings on and is heard mostly by its persistence).
6+
7+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
8+
9+interface Props {
10+ modeF: Float32Array;
11+ modeZeta: Float32Array;
12+ modeZetaRad: Float32Array;
13+ modeFamily: Float32Array;
14+}
15+
16+const F_LO = 100;
17+const F_HI = 24000;
18+
19+export default function ModeMap({ modeF, modeZeta, modeZetaRad, modeFamily }: Props) {
20+ const ref = useCanvas(
21+ (ctx, w, h) => {
22+ drawFrame(ctx, w, h);
23+ const padL = 40;
24+ const padB = 20;
25+ const pw = w - padL - 10;
26+ const ph = h - padB - 14;
27+ const top = 12;
28+ const n = modeF.length;
29+ if (n === 0) return;
30+
31+ const t60 = new Float64Array(n);
32+ let tMin = Infinity;
33+ let tMax = 0;
34+ for (let i = 0; i < n; i++) {
35+ t60[i] = 6.91 / (Math.max(1e-9, modeZeta[i]) * 2 * Math.PI * modeF[i]);
36+ tMin = Math.min(tMin, t60[i]);
37+ tMax = Math.max(tMax, t60[i]);
38+ }
39+ tMin = Math.max(1e-3, tMin / 1.5);
40+ tMax = tMax * 1.5;
41+
42+ const fx = (f: number) =>
43+ padL + ((Math.log10(f) - Math.log10(F_LO)) / (Math.log10(F_HI) - Math.log10(F_LO))) * pw;
44+ const ty = (t: number) =>
45+ top + ph - ((Math.log10(t) - Math.log10(tMin)) / (Math.log10(tMax) - Math.log10(tMin))) * ph;
46+
47+ ctx.font = '10px system-ui';
48+ ctx.textAlign = 'center';
49+ for (const f of [100, 200, 500, 1000, 2000, 5000, 10000, 20000]) {
50+ const x = fx(f);
51+ ctx.strokeStyle = PLOT_COLORS.grid;
52+ ctx.beginPath();
53+ ctx.moveTo(x, top);
54+ ctx.lineTo(x, top + ph);
55+ ctx.stroke();
56+ ctx.fillStyle = PLOT_COLORS.axis;
57+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, x, h - 7);
58+ }
59+ ctx.textAlign = 'right';
60+ for (let e = Math.ceil(Math.log10(tMin)); e <= Math.floor(Math.log10(tMax)); e++) {
61+ const t = 10 ** e;
62+ const y = ty(t);
63+ ctx.strokeStyle = PLOT_COLORS.grid;
64+ ctx.beginPath();
65+ ctx.moveTo(padL, y);
66+ ctx.lineTo(padL + pw, y);
67+ ctx.stroke();
68+ ctx.fillStyle = PLOT_COLORS.axis;
69+ ctx.fillText(t >= 1 ? `${t}s` : `${t * 1000}ms`, padL - 4, y + 3);
70+ }
71+
72+ for (let i = 0; i < n; i++) {
73+ const x = fx(Math.min(F_HI, Math.max(F_LO, modeF[i])));
74+ const y = ty(t60[i]);
75+ const radShare = Math.min(1, modeZetaRad[i] / Math.max(1e-12, modeZeta[i]));
76+ const r = 2 + 5 * Math.sqrt(radShare);
77+ ctx.fillStyle = modeFamily[i] === 0 ? PLOT_COLORS.trace : PLOT_COLORS.accent;
78+ ctx.globalAlpha = 0.85;
79+ ctx.beginPath();
80+ ctx.arc(x, y, r, 0, 2 * Math.PI);
81+ ctx.fill();
82+ ctx.globalAlpha = 1;
83+ }
84+
85+ ctx.fillStyle = PLOT_COLORS.trace;
86+ ctx.textAlign = 'left';
87+ ctx.fillText('in-plane', padL + 6, top + 2);
88+ ctx.fillStyle = PLOT_COLORS.accent;
89+ ctx.fillText('out-of-plane', padL + 56, top + 2);
90+ ctx.fillStyle = PLOT_COLORS.axis;
91+ ctx.fillText(`· ${n} modes · marker = radiation share of damping`, padL + 128, top + 2);
92+ },
93+ [modeF, modeZeta, modeZetaRad, modeFamily],
94+ 210,
95+ );
96+ return <canvas ref={ref} />;
97+}
webapp/src/components/plots/ModeShapes.tsxadded+122−0View file
@@ -0,0 +1,122 @@
1+// Small multiples of the lowest mode shapes on the triangle outline.
2+// In-plane modes are drawn as the outline displaced within the plane;
3+// out-of-plane modes keep the outline in place and color it by the normal
4+// displacement (blue toward the viewer, red away), since that motion is
5+// perpendicular to the drawing. The gap at the bottom-left corner is where
6+// the free ends whip.
7+
8+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
9+
10+interface Props {
11+ nodeX: Float32Array;
12+ nodeY: Float32Array;
13+ shapeUx: Float32Array; // nn x n, column-major
14+ shapeUy: Float32Array;
15+ shapeUz: Float32Array;
16+ modeF: Float32Array;
17+ modeFamily: Float32Array;
18+}
19+
20+const SHOWN = 12;
21+const CELL_H = 110;
22+
23+export default function ModeShapes({
24+ nodeX,
25+ nodeY,
26+ shapeUx,
27+ shapeUy,
28+ shapeUz,
29+ modeF,
30+ modeFamily,
31+}: Props) {
32+ const nShown = Math.min(SHOWN, modeF.length);
33+ const rows = Math.ceil(nShown / 4);
34+ const ref = useCanvas(
35+ (ctx, w, h) => {
36+ drawFrame(ctx, w, h);
37+ const nn = nodeX.length;
38+ if (nn === 0 || nShown === 0) return;
39+ const cols = 4;
40+ const cw = w / cols;
41+ const ch = h / rows;
42+
43+ // model bounding box, common to all cells
44+ let xMin = Infinity, xMax = -Infinity, yMin = Infinity, yMax = -Infinity;
45+ for (let i = 0; i < nn; i++) {
46+ xMin = Math.min(xMin, nodeX[i]);
47+ xMax = Math.max(xMax, nodeX[i]);
48+ yMin = Math.min(yMin, nodeY[i]);
49+ yMax = Math.max(yMax, nodeY[i]);
50+ }
51+ const span = Math.max(xMax - xMin, yMax - yMin);
52+ const scale = (Math.min(cw, ch) - 34) / span;
53+ const defl = span * 0.1; // displacement drawn as 10% of the frame size
54+
55+ for (let m = 0; m < nShown; m++) {
56+ const cx = (m % cols) * cw + cw / 2;
57+ const cy = Math.floor(m / cols) * ch + ch / 2 + 4;
58+ const px = (x: number) => cx + (x - (xMin + xMax) / 2) * scale;
59+ const py = (y: number) => cy - (y - (yMin + yMax) / 2) * scale;
60+ const off = m * nn;
61+
62+ let uMax = 1e-12;
63+ for (let i = 0; i < nn; i++) {
64+ uMax = Math.max(
65+ uMax,
66+ Math.hypot(shapeUx[off + i], shapeUy[off + i], shapeUz[off + i]),
67+ );
68+ }
69+
70+ // undeformed outline, faint
71+ ctx.strokeStyle = 'rgba(139, 147, 161, 0.25)';
72+ ctx.lineWidth = 1;
73+ ctx.beginPath();
74+ for (let i = 0; i < nn; i++) {
75+ const x = px(nodeX[i]);
76+ const y = py(nodeY[i]);
77+ if (i === 0) ctx.moveTo(x, y);
78+ else ctx.lineTo(x, y);
79+ }
80+ ctx.stroke();
81+
82+ if (modeFamily[m] === 0) {
83+ // in-plane: displaced outline
84+ ctx.strokeStyle = PLOT_COLORS.trace;
85+ ctx.lineWidth = 1.6;
86+ ctx.beginPath();
87+ for (let i = 0; i < nn; i++) {
88+ const x = px(nodeX[i] + (shapeUx[off + i] / uMax) * defl);
89+ const y = py(nodeY[i] + (shapeUy[off + i] / uMax) * defl);
90+ if (i === 0) ctx.moveTo(x, y);
91+ else ctx.lineTo(x, y);
92+ }
93+ ctx.stroke();
94+ } else {
95+ // out-of-plane: color the outline by normal displacement
96+ ctx.lineWidth = 2.2;
97+ for (let i = 0; i < nn - 1; i++) {
98+ const t = (shapeUz[off + i] + shapeUz[off + i + 1]) / (2 * uMax);
99+ const mag = Math.min(1, Math.abs(t));
100+ ctx.strokeStyle =
101+ t >= 0
102+ ? `rgba(100, 181, 246, ${0.25 + 0.75 * mag})`
103+ : `rgba(229, 115, 115, ${0.25 + 0.75 * mag})`;
104+ ctx.beginPath();
105+ ctx.moveTo(px(nodeX[i]), py(nodeY[i]));
106+ ctx.lineTo(px(nodeX[i + 1]), py(nodeY[i + 1]));
107+ ctx.stroke();
108+ }
109+ }
110+
111+ ctx.font = '10px system-ui';
112+ ctx.textAlign = 'center';
113+ ctx.fillStyle = modeFamily[m] === 0 ? PLOT_COLORS.trace : PLOT_COLORS.accent;
114+ const fam = modeFamily[m] === 0 ? 'in' : 'out';
115+ ctx.fillText(`${modeF[m].toFixed(0)} Hz · ${fam}`, cx, cy + (ch / 2) - 8);
116+ }
117+ },
118+ [nodeX, nodeY, shapeUx, shapeUy, shapeUz, modeF, modeFamily, nShown, rows],
119+ rows * CELL_H,
120+ );
121+ return <canvas ref={ref} />;
122+}
webapp/src/components/plots/SpectrogramPlot.tsxadded+64−0View file
@@ -0,0 +1,64 @@
1+// Time-frequency view of the rendered strike; per-partial decay is visible as
2+// horizontal streaks of different lengths.
3+
4+import { stft } from '../../lib/fft';
5+import { useCanvas, drawFrame, PLOT_COLORS, colormap } from './canvasUtils';
6+
7+interface Props {
8+ y: Float32Array;
9+ sampleRate: number;
10+}
11+
12+export default function SpectrogramPlot({ y, sampleRate }: Props) {
13+ const ref = useCanvas(
14+ (ctx, w, h) => {
15+ drawFrame(ctx, w, h);
16+ const padL = 34;
17+ const padB = 16;
18+ const pw = w - padL - 6;
19+ const ph = h - padB - 6;
20+
21+ const { frames, bins, db, floorDb } = stft(y, 512, 128);
22+
23+ const img = ctx.createImageData(frames, bins);
24+ for (let f = 0; f < frames; f++) {
25+ for (let b = 0; b < bins; b++) {
26+ const v = db[f * bins + b];
27+ const t = Math.min(1, Math.max(0, (v - floorDb) / -floorDb));
28+ const [r, g, bl] = colormap(t);
29+ // flip vertically: low frequency at the bottom
30+ const row = bins - 1 - b;
31+ const idx = (row * frames + f) * 4;
32+ img.data[idx] = r;
33+ img.data[idx + 1] = g;
34+ img.data[idx + 2] = bl;
35+ img.data[idx + 3] = 255;
36+ }
37+ }
38+ const off = new OffscreenCanvas(frames, bins);
39+ const octx = off.getContext('2d')!;
40+ octx.putImageData(img, 0, 0);
41+ ctx.imageSmoothingEnabled = true;
42+ ctx.drawImage(off, padL, 6, pw, ph);
43+
44+ // axes
45+ ctx.fillStyle = PLOT_COLORS.axis;
46+ ctx.font = '10px system-ui';
47+ ctx.textAlign = 'right';
48+ const nyq = sampleRate / 2;
49+ for (let kf = 0; kf <= nyq; kf += 4000) {
50+ const yy = 6 + ph - (kf / nyq) * ph;
51+ ctx.fillText(kf >= 1000 ? `${kf / 1000}k` : `${kf}`, padL - 4, yy + 3);
52+ }
53+ ctx.textAlign = 'center';
54+ const durS = y.length / sampleRate;
55+ for (let t = 0; t <= durS + 1e-9; t += 0.1) {
56+ const x = padL + (t / durS) * pw;
57+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
58+ }
59+ },
60+ [y, sampleRate],
61+ 200,
62+ );
63+ return <canvas ref={ref} />;
64+}
webapp/src/components/plots/SpectrumPlot.tsxadded+108−0View file
@@ -0,0 +1,108 @@
1+// Magnitude spectrum of the rendered strike, log-frequency axis, with the
2+// retained rod modes overlaid. A triangle's spectrum is a sparse comb of
3+// inharmonic partials with almost nothing between them; the broadband floor
4+// near the attack is the direct click standing in for the truncated modes.
5+
6+import { magnitudeSpectrum } from '../../lib/fft';
7+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
8+
9+interface Props {
10+ y: Float32Array;
11+ sampleRate: number;
12+ modeF: Float32Array;
13+ modeFamily: Float32Array;
14+}
15+
16+const F_LO = 100;
17+const F_HI = 24000;
18+const DB_RANGE = 80;
19+
20+export default function SpectrumPlot({ y, sampleRate, modeF, modeFamily }: Props) {
21+ const ref = useCanvas(
22+ (ctx, w, h) => {
23+ drawFrame(ctx, w, h);
24+ const padL = 34;
25+ const padB = 18;
26+ const pw = w - padL - 6;
27+ const ph = h - padB - 8;
28+ const fx = (f: number) =>
29+ padL + ((Math.log10(f) - Math.log10(F_LO)) / (Math.log10(F_HI) - Math.log10(F_LO))) * pw;
30+
31+ ctx.font = '10px system-ui';
32+ ctx.textAlign = 'center';
33+ for (const f of [100, 200, 500, 1000, 2000, 5000, 10000, 20000]) {
34+ const x = fx(f);
35+ ctx.strokeStyle = PLOT_COLORS.grid;
36+ ctx.beginPath();
37+ ctx.moveTo(x, 8);
38+ ctx.lineTo(x, 8 + ph);
39+ ctx.stroke();
40+ ctx.fillStyle = PLOT_COLORS.axis;
41+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, x, h - 5);
42+ }
43+ ctx.textAlign = 'right';
44+ for (let db = 0; db >= -DB_RANGE; db -= 20) {
45+ const yy = 8 + (-db / DB_RANGE) * ph;
46+ ctx.strokeStyle = PLOT_COLORS.grid;
47+ ctx.beginPath();
48+ ctx.moveTo(padL, yy);
49+ ctx.lineTo(padL + pw, yy);
50+ ctx.stroke();
51+ ctx.fillStyle = PLOT_COLORS.axis;
52+ ctx.fillText(`${db}`, padL - 4, yy + 3);
53+ }
54+
55+ const { mag, size } = magnitudeSpectrum(y, 32768);
56+ let peak = 0;
57+ for (let i = 1; i < mag.length; i++) peak = Math.max(peak, mag[i]);
58+ if (peak <= 0) return;
59+
60+ // overlays first so the trace draws on top
61+ ctx.setLineDash([4, 3]);
62+ for (let i = 0; i < modeF.length; i++) {
63+ const f = modeF[i];
64+ if (f < F_LO || f > F_HI) continue;
65+ ctx.strokeStyle = modeFamily[i] === 0 ? 'rgba(77, 208, 225, 0.45)' : 'rgba(232, 163, 61, 0.45)';
66+ ctx.beginPath();
67+ ctx.moveTo(fx(f), 8);
68+ ctx.lineTo(fx(f), 8 + ph);
69+ ctx.stroke();
70+ }
71+ ctx.setLineDash([]);
72+
73+ ctx.strokeStyle = PLOT_COLORS.trace;
74+ ctx.beginPath();
75+ let started = false;
76+ for (let px = 0; px <= pw; px++) {
77+ const f = 10 ** (Math.log10(F_LO) + (px / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
78+ const bin = Math.min(mag.length - 1, Math.round((f / sampleRate) * size));
79+ const fNext =
80+ 10 ** (Math.log10(F_LO) + ((px + 1) / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
81+ const binNext = Math.min(
82+ mag.length - 1,
83+ Math.max(bin + 1, Math.round((fNext / sampleRate) * size)),
84+ );
85+ let m = 0;
86+ for (let b = bin; b < binNext; b++) m = Math.max(m, mag[b]);
87+ const db = 20 * Math.log10(m / peak + 1e-12);
88+ const yy = 8 + Math.min(ph, (-db / DB_RANGE) * ph);
89+ if (!started) {
90+ ctx.moveTo(padL + px, yy);
91+ started = true;
92+ } else {
93+ ctx.lineTo(padL + px, yy);
94+ }
95+ }
96+ ctx.stroke();
97+
98+ ctx.textAlign = 'left';
99+ ctx.fillStyle = 'rgba(77, 208, 225, 0.9)';
100+ ctx.fillText('in-plane modes', padL + 4, 18);
101+ ctx.fillStyle = 'rgba(232, 163, 61, 0.9)';
102+ ctx.fillText('out-of-plane modes', padL + 84, 18);
103+ },
104+ [y, sampleRate, modeF, modeFamily],
105+ 200,
106+ );
107+ return <canvas ref={ref} />;
108+}
webapp/src/components/plots/StrikePlot.tsxadded+77−0View file
@@ -0,0 +1,77 @@
1+// The strike events: each hit as a stem at its time and velocity, colored by
2+// which contact point it lands on. For a single strike this is one stem; for
3+// a roll the jitter and velocity spread are what keep it from sounding like
4+// a machine.
5+
6+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
7+
8+interface Props {
9+ strikeT: Float32Array;
10+ strikeV: Float32Array;
11+ strikeS: Float32Array;
12+ durationS: number;
13+}
14+
15+export default function StrikePlot({ strikeT, strikeV, strikeS, durationS }: Props) {
16+ const ref = useCanvas(
17+ (ctx, w, h) => {
18+ drawFrame(ctx, w, h);
19+ const padL = 40;
20+ const padB = 16;
21+ const pw = w - padL - 8;
22+ const ph = h - padB - 12;
23+ const top = 10;
24+
25+ let vMax = 0;
26+ for (let i = 0; i < strikeV.length; i++) vMax = Math.max(vMax, strikeV[i]);
27+ if (vMax <= 0) vMax = 1;
28+
29+ const tx = (t: number) => padL + (t / durationS) * pw;
30+
31+ ctx.font = '10px system-ui';
32+ ctx.textAlign = 'center';
33+ const tStep = durationS > 4 ? 1.0 : durationS > 1.5 ? 0.5 : 0.1;
34+ for (let t = 0; t <= durationS + 1e-9; t += tStep) {
35+ const x = tx(t);
36+ ctx.strokeStyle = PLOT_COLORS.grid;
37+ ctx.beginPath();
38+ ctx.moveTo(x, top);
39+ ctx.lineTo(x, top + ph);
40+ ctx.stroke();
41+ ctx.fillStyle = PLOT_COLORS.axis;
42+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
43+ }
44+ ctx.textAlign = 'right';
45+ ctx.fillStyle = PLOT_COLORS.axis;
46+ ctx.fillText(`${vMax.toFixed(1)} m/s`, padL - 4, top + 4);
47+ ctx.fillText('0', padL - 4, top + ph);
48+
49+ const s0 = strikeS.length > 0 ? strikeS[0] : 0;
50+ for (let i = 0; i < strikeT.length; i++) {
51+ const x = tx(strikeT[i]);
52+ const y = top + ph - (strikeV[i] / vMax) * ph * 0.92;
53+ const alt = Math.abs(strikeS[i] - s0) > 1e-9;
54+ ctx.strokeStyle = alt ? 'rgba(232, 163, 61, 0.85)' : 'rgba(77, 208, 225, 0.85)';
55+ ctx.beginPath();
56+ ctx.moveTo(x, top + ph);
57+ ctx.lineTo(x, y);
58+ ctx.stroke();
59+ ctx.fillStyle = alt ? PLOT_COLORS.accent : PLOT_COLORS.trace;
60+ ctx.beginPath();
61+ ctx.arc(x, y, 2.5, 0, 2 * Math.PI);
62+ ctx.fill();
63+ }
64+
65+ ctx.fillStyle = PLOT_COLORS.axis;
66+ ctx.textAlign = 'left';
67+ const note =
68+ strikeT.length > 1
69+ ? `${strikeT.length} strikes · cyan = strike position, amber = second contact`
70+ : `single strike at t = 0`;
71+ ctx.fillText(note, padL + 6, top + 4);
72+ },
73+ [strikeT, strikeV, strikeS, durationS],
74+ 150,
75+ );
76+ return <canvas ref={ref} />;
77+}
webapp/src/components/plots/WaveformPlot.tsxadded+64−0View file
@@ -0,0 +1,64 @@
1+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
2+
3+interface Props {
4+ y: Float32Array;
5+ sampleRate: number;
6+}
7+
8+export default function WaveformPlot({ y, sampleRate }: Props) {
9+ const ref = useCanvas(
10+ (ctx, w, h) => {
11+ drawFrame(ctx, w, h);
12+ const padL = 34;
13+ const padB = 16;
14+ const pw = w - padL - 6;
15+ const ph = h - padB - 6;
16+ const midY = 6 + ph / 2;
17+
18+ ctx.strokeStyle = PLOT_COLORS.grid;
19+ ctx.beginPath();
20+ ctx.moveTo(padL, midY);
21+ ctx.lineTo(padL + pw, midY);
22+ ctx.stroke();
23+
24+ // time gridlines every 0.1 s
25+ const durS = y.length / sampleRate;
26+ ctx.fillStyle = PLOT_COLORS.axis;
27+ ctx.font = '10px system-ui';
28+ ctx.textAlign = 'center';
29+ for (let t = 0; t <= durS + 1e-9; t += 0.1) {
30+ const x = padL + (t / durS) * pw;
31+ ctx.strokeStyle = PLOT_COLORS.grid;
32+ ctx.beginPath();
33+ ctx.moveTo(x, 6);
34+ ctx.lineTo(x, 6 + ph);
35+ ctx.stroke();
36+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
37+ }
38+ ctx.textAlign = 'right';
39+ ctx.fillText('+1', padL - 4, 12);
40+ ctx.fillText('-1', padL - 4, 6 + ph);
41+
42+ // min/max envelope per pixel column
43+ ctx.strokeStyle = PLOT_COLORS.trace;
44+ ctx.beginPath();
45+ const perPx = y.length / pw;
46+ for (let px = 0; px < pw; px++) {
47+ const i0 = Math.floor(px * perPx);
48+ const i1 = Math.min(y.length, Math.max(i0 + 1, Math.floor((px + 1) * perPx)));
49+ let mn = Infinity, mx = -Infinity;
50+ for (let i = i0; i < i1; i++) {
51+ if (y[i] < mn) mn = y[i];
52+ if (y[i] > mx) mx = y[i];
53+ }
54+ if (!Number.isFinite(mn)) continue;
55+ ctx.moveTo(padL + px + 0.5, midY - mx * (ph / 2));
56+ ctx.lineTo(padL + px + 0.5, midY - mn * (ph / 2) + 0.5);
57+ }
58+ ctx.stroke();
59+ },
60+ [y, sampleRate],
61+ 170,
62+ );
63+ return <canvas ref={ref} />;
64+}
webapp/src/components/plots/canvasUtils.tsadded+66−0View file
@@ -0,0 +1,66 @@
1+import { useEffect, useRef } from 'react';
2+
3+// Canvas hook: resizes to element width at devicePixelRatio and calls draw.
4+export function useCanvas(
5+ draw: (ctx: CanvasRenderingContext2D, w: number, h: number) => void,
6+ deps: unknown[],
7+ height: number,
8+) {
9+ const ref = useRef<HTMLCanvasElement | null>(null);
10+
11+ useEffect(() => {
12+ const canvas = ref.current;
13+ if (!canvas) return;
14+ const render = () => {
15+ const dpr = window.devicePixelRatio || 1;
16+ const w = canvas.clientWidth;
17+ canvas.width = Math.max(1, Math.round(w * dpr));
18+ canvas.height = Math.max(1, Math.round(height * dpr));
19+ canvas.style.height = `${height}px`;
20+ const ctx = canvas.getContext('2d')!;
21+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
22+ ctx.clearRect(0, 0, w, height);
23+ draw(ctx, w, height);
24+ };
25+ render();
26+ const ro = new ResizeObserver(render);
27+ ro.observe(canvas);
28+ return () => ro.disconnect();
29+ // eslint-disable-next-line react-hooks/exhaustive-deps
30+ }, deps);
31+
32+ return ref;
33+}
34+
35+export const PLOT_COLORS = {
36+ bg: '#181b20',
37+ grid: '#2a2f38',
38+ axis: '#8b93a1',
39+ trace: '#4dd0e1',
40+ accent: '#e8a33d',
41+ good: '#6fbf73',
42+ warn: '#e57373',
43+};
44+
45+export function drawFrame(ctx: CanvasRenderingContext2D, w: number, h: number) {
46+ ctx.fillStyle = PLOT_COLORS.bg;
47+ ctx.fillRect(0, 0, w, h);
48+}
49+
50+// Viridis-like colormap for the spectrogram, t in [0,1].
51+export function colormap(t: number): [number, number, number] {
52+ const stops: [number, number, number][] = [
53+ [13, 8, 60], [84, 39, 143], [140, 41, 129], [196, 66, 90],
54+ [237, 105, 37], [251, 155, 6], [247, 209, 61], [252, 255, 164],
55+ ];
56+ const x = Math.min(0.9999, Math.max(0, t)) * (stops.length - 1);
57+ const i = Math.floor(x);
58+ const f = x - i;
59+ const a = stops[i];
60+ const b = stops[i + 1];
61+ return [
62+ Math.round(a[0] + f * (b[0] - a[0])),
63+ Math.round(a[1] + f * (b[1] - a[1])),
64+ Math.round(a[2] + f * (b[2] - a[2])),
65+ ];
66+}
webapp/src/hooks/useRenderer.tsadded+114−0View file
@@ -0,0 +1,114 @@
1+// Hook owning the render worker: debounced auto-preview whenever the
2+// parameters or the .m sources change, explicit full-quality render, and a
3+// latest-wins request queue.
4+
5+import { useCallback, useEffect, useRef, useState } from 'react';
6+import type { TriangleParams } from '../params';
7+import type { MFile } from '../mfiles';
8+import type { ErrorMessage, ProgressMessage, RenderResponse } from '../worker/renderWorker';
9+
10+export interface RenderError {
11+ message: string;
12+ file?: string;
13+ line?: number;
14+}
15+
16+export interface RenderState {
17+ result: RenderResponse | null;
18+ error: RenderError | null;
19+ busy: boolean;
20+ stage: string;
21+ progress: number; // 0..1 of the current stage
22+}
23+
24+interface Request {
25+ params: TriangleParams;
26+ mfiles: MFile[];
27+ quality: 'preview' | 'full';
28+}
29+
30+export function useRenderer(params: TriangleParams, mfiles: MFile[]) {
31+ const workerRef = useRef<Worker | null>(null);
32+ const nextId = useRef(1);
33+ const inFlight = useRef<number | null>(null);
34+ const queued = useRef<Request | null>(null);
35+ const sendRef = useRef<(r: Request) => void>(() => {});
36+ const [state, setState] = useState<RenderState>({
37+ result: null,
38+ error: null,
39+ busy: false,
40+ stage: '',
41+ progress: 0,
42+ });
43+
44+ const send = useCallback((req: Request) => {
45+ const worker = workerRef.current;
46+ if (!worker) return;
47+ if (inFlight.current !== null) {
48+ queued.current = req; // latest wins
49+ return;
50+ }
51+ const id = nextId.current++;
52+ inFlight.current = id;
53+ setState((s) => ({ ...s, busy: true, stage: 'starting', progress: 0 }));
54+ worker.postMessage({ type: 'render', id, ...req });
55+ }, []);
56+ sendRef.current = send;
57+
58+ useEffect(() => {
59+ const worker = new Worker(new URL('../worker/renderWorker.ts', import.meta.url), {
60+ type: 'module',
61+ });
62+ workerRef.current = worker;
63+
64+ const finish = (update: (s: RenderState) => RenderState) => {
65+ inFlight.current = null;
66+ const next = queued.current;
67+ queued.current = null;
68+ setState((s) => {
69+ const done = update(s);
70+ return next ? { ...done, busy: true, stage: 'queued', progress: 0 } : done;
71+ });
72+ if (next) sendRef.current(next);
73+ };
74+
75+ worker.onmessage = (ev: MessageEvent<RenderResponse | ProgressMessage | ErrorMessage>) => {
76+ const msg = ev.data;
77+ if (msg.type === 'progress') {
78+ setState((s) => ({ ...s, stage: msg.stage, progress: msg.frac }));
79+ return;
80+ }
81+ if (msg.type === 'error') {
82+ // Keep the last good result on screen; the code panel shows the fault.
83+ finish((s) => ({
84+ ...s,
85+ error: { message: msg.message, file: msg.file, line: msg.line },
86+ busy: false,
87+ stage: '',
88+ progress: 0,
89+ }));
90+ return;
91+ }
92+ finish(() => ({ result: msg, error: null, busy: false, stage: '', progress: 0 }));
93+ };
94+ return () => {
95+ worker.terminate();
96+ workerRef.current = null;
97+ inFlight.current = null;
98+ queued.current = null;
99+ };
100+ }, []);
101+
102+ const renderFull = useCallback(
103+ () => send({ params, mfiles, quality: 'full' }),
104+ [params, mfiles, send],
105+ );
106+
107+ // Debounced auto-preview on any change to the parameters or the code.
108+ useEffect(() => {
109+ const t = setTimeout(() => send({ params, mfiles, quality: 'preview' }), 250);
110+ return () => clearTimeout(t);
111+ }, [params, mfiles, send]);
112+
113+ return { ...state, renderFull };
114+}
webapp/src/lib/fft.tsadded+91−0View file
@@ -0,0 +1,91 @@
1+// Minimal iterative radix-2 complex FFT, in-place on separate re/im arrays.
2+
3+export function fft(re: Float64Array, im: Float64Array): void {
4+ const n = re.length;
5+ if ((n & (n - 1)) !== 0) throw new Error('FFT length must be a power of 2');
6+ // bit-reversal permutation
7+ for (let i = 1, j = 0; i < n; i++) {
8+ let bit = n >> 1;
9+ for (; j & bit; bit >>= 1) j ^= bit;
10+ j ^= bit;
11+ if (i < j) {
12+ const tr = re[i]; re[i] = re[j]; re[j] = tr;
13+ const ti = im[i]; im[i] = im[j]; im[j] = ti;
14+ }
15+ }
16+ for (let len = 2; len <= n; len <<= 1) {
17+ const ang = (-2 * Math.PI) / len;
18+ const wr = Math.cos(ang);
19+ const wi = Math.sin(ang);
20+ for (let i = 0; i < n; i += len) {
21+ let cwr = 1, cwi = 0;
22+ for (let k = 0; k < len / 2; k++) {
23+ const ur = re[i + k], ui = im[i + k];
24+ const vr = re[i + k + len / 2] * cwr - im[i + k + len / 2] * cwi;
25+ const vi = re[i + k + len / 2] * cwi + im[i + k + len / 2] * cwr;
26+ re[i + k] = ur + vr;
27+ im[i + k] = ui + vi;
28+ re[i + k + len / 2] = ur - vr;
29+ im[i + k + len / 2] = ui - vi;
30+ const nwr = cwr * wr - cwi * wi;
31+ cwi = cwr * wi + cwi * wr;
32+ cwr = nwr;
33+ }
34+ }
35+ }
36+}
37+
38+// Magnitude spectrum of a real signal with a Hann window, zero-padded to the
39+// next power of two (at least minSize). Returns bins 0..N/2.
40+export function magnitudeSpectrum(x: Float32Array | Float64Array, minSize = 16384): {
41+ mag: Float64Array;
42+ size: number;
43+} {
44+ let size = minSize;
45+ while (size < x.length) size <<= 1;
46+ const re = new Float64Array(size);
47+ const im = new Float64Array(size);
48+ const n = x.length;
49+ for (let i = 0; i < n; i++) {
50+ const w = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (n - 1)));
51+ re[i] = x[i] * w;
52+ }
53+ fft(re, im);
54+ const half = size / 2;
55+ const mag = new Float64Array(half + 1);
56+ for (let i = 0; i <= half; i++) mag[i] = Math.hypot(re[i], im[i]);
57+ return { mag, size };
58+}
59+
60+// Simple STFT magnitude in dB for the spectrogram, frames x bins.
61+export function stft(
62+ x: Float32Array,
63+ win = 512,
64+ hop = 128,
65+): { frames: number; bins: number; db: Float32Array; floorDb: number } {
66+ const nFrames = Math.max(1, Math.floor((x.length - win) / hop) + 1);
67+ const bins = win / 2 + 1;
68+ const db = new Float32Array(nFrames * bins);
69+ const re = new Float64Array(win);
70+ const im = new Float64Array(win);
71+ const hann = new Float64Array(win);
72+ for (let i = 0; i < win; i++) hann[i] = 0.5 * (1 - Math.cos((2 * Math.PI * i) / (win - 1)));
73+ let maxDb = -Infinity;
74+ for (let f = 0; f < nFrames; f++) {
75+ const off = f * hop;
76+ for (let i = 0; i < win; i++) {
77+ re[i] = (off + i < x.length ? x[off + i] : 0) * hann[i];
78+ im[i] = 0;
79+ }
80+ fft(re, im);
81+ for (let b = 0; b < bins; b++) {
82+ const m = Math.hypot(re[b], im[b]);
83+ const v = 20 * Math.log10(m + 1e-12);
84+ db[f * bins + b] = v;
85+ if (v > maxDb) maxDb = v;
86+ }
87+ }
88+ // normalize so 0 dB is the frame-bin maximum
89+ for (let i = 0; i < db.length; i++) db[i] -= maxDb;
90+ return { frames: nFrames, bins, db, floorDb: -90 };
91+}
webapp/src/lib/outline.tsadded+60−0View file
@@ -0,0 +1,60 @@
1+// Mirror of ../matlab/triangle_geometry.m: the open-triangle centerline as a
2+// polyline in the x-y plane, apex up, gap at the bottom-left corner, running
3+// counterclockwise from the free end on the base (arclength fraction 0) to
4+// the free end on the left leg (fraction 1). Used by the 3D scene and the
5+// mode-shape plot; the physics reads the same geometry from the MATLAB side.
6+
7+export interface OutlinePoint {
8+ x: number;
9+ y: number;
10+ tx: number; // unit tangent along increasing arclength
11+ ty: number;
12+}
13+
14+export interface Outline {
15+ pts: { x: number; y: number }[]; // A, V3, V1, B
16+ segLens: number[];
17+ totalLen: number;
18+ at(frac: number): OutlinePoint;
19+}
20+
21+export function triangleOutline(sideLength: number, gapM: number): Outline {
22+ const L = sideLength;
23+ const g = Math.min(gapM, 0.8 * L);
24+ const V1 = { x: 0, y: (L * Math.sqrt(3)) / 2 };
25+ const V2 = { x: -L / 2, y: 0 };
26+ const V3 = { x: L / 2, y: 0 };
27+ const e1 = { x: 1, y: 0 };
28+ const e3 = { x: (V2.x - V1.x) / L, y: (V2.y - V1.y) / L };
29+ const A = { x: V2.x + (e1.x * g) / 2, y: V2.y + (e1.y * g) / 2 };
30+ const B = { x: V2.x - (e3.x * g) / 2, y: V2.y - (e3.y * g) / 2 };
31+ const pts = [A, V3, V1, B];
32+ const segLens = [L - g / 2, L, L - g / 2];
33+ const totalLen = segLens[0] + segLens[1] + segLens[2];
34+ return {
35+ pts,
36+ segLens,
37+ totalLen,
38+ at(frac: number): OutlinePoint {
39+ let s = Math.min(0.999, Math.max(0.001, frac)) * totalLen;
40+ for (let k = 0; k < 3; k++) {
41+ if (s <= segLens[k] || k === 2) {
42+ const a = pts[k];
43+ const b = pts[k + 1];
44+ const t = Math.min(1, s / segLens[k]);
45+ const dx = b.x - a.x;
46+ const dy = b.y - a.y;
47+ const dl = Math.hypot(dx, dy);
48+ return {
49+ x: a.x + t * dx,
50+ y: a.y + t * dy,
51+ tx: dx / dl,
52+ ty: dy / dl,
53+ };
54+ }
55+ s -= segLens[k];
56+ }
57+ throw new Error('unreachable');
58+ },
59+ };
60+}
webapp/src/lib/wav.tsadded+40−0View file
@@ -0,0 +1,40 @@
1+// Mono PCM WAV encoding (16- or 24-bit).
2+
3+export function encodeWav(y: Float32Array, sampleRate: number, bitDepth: 16 | 24): Blob {
4+ const bytesPerSample = bitDepth / 8;
5+ const dataSize = y.length * bytesPerSample;
6+ const buf = new ArrayBuffer(44 + dataSize);
7+ const dv = new DataView(buf);
8+ const writeStr = (off: number, s: string) => {
9+ for (let i = 0; i < s.length; i++) dv.setUint8(off + i, s.charCodeAt(i));
10+ };
11+ writeStr(0, 'RIFF');
12+ dv.setUint32(4, 36 + dataSize, true);
13+ writeStr(8, 'WAVE');
14+ writeStr(12, 'fmt ');
15+ dv.setUint32(16, 16, true);
16+ dv.setUint16(20, 1, true); // PCM
17+ dv.setUint16(22, 1, true); // mono
18+ dv.setUint32(24, sampleRate, true);
19+ dv.setUint32(28, sampleRate * bytesPerSample, true);
20+ dv.setUint16(32, bytesPerSample, true);
21+ dv.setUint16(34, bitDepth, true);
22+ writeStr(36, 'data');
23+ dv.setUint32(40, dataSize, true);
24+
25+ let off = 44;
26+ for (let i = 0; i < y.length; i++) {
27+ const v = Math.max(-1, Math.min(1, y[i]));
28+ if (bitDepth === 16) {
29+ dv.setInt16(off, Math.round(v * 32767), true);
30+ off += 2;
31+ } else {
32+ const iv = Math.round(v * 8388607);
33+ dv.setUint8(off, iv & 0xff);
34+ dv.setUint8(off + 1, (iv >> 8) & 0xff);
35+ dv.setUint8(off + 2, (iv >> 16) & 0xff);
36+ off += 3;
37+ }
38+ }
39+ return new Blob([buf], { type: 'audio/wav' });
40+}
webapp/src/main.tsxadded+10−0View file
@@ -0,0 +1,10 @@
1+import React from 'react';
2+import { createRoot } from 'react-dom/client';
3+import App from './App';
4+import './style.css';
5+
6+createRoot(document.getElementById('root')!).render(
7+ <React.StrictMode>
8+ <App />
9+ </React.StrictMode>,
10+);
webapp/src/mfiles.tsadded+44−0View file
@@ -0,0 +1,44 @@
1+// The MATLAB implementation, loaded from ../matlab at build time. The app
2+// keeps an editable in-memory copy of these sources (see App.tsx); edits feed
3+// straight into the render worker and are deliberately not persisted, so a
4+// reload restores the shipped code.
5+
6+export interface MFile {
7+ name: string;
8+ source: string;
9+}
10+
11+const modules = import.meta.glob('../../matlab/*.m', {
12+ eager: true,
13+ query: '?raw',
14+ import: 'default',
15+}) as Record<string, string>;
16+
17+export const DEFAULT_MFILES: MFile[] = Object.entries(modules)
18+ .map(([path, source]) => ({ name: path.split('/').pop()!, source }))
19+ .sort((a, b) => a.name.localeCompare(b.name));
20+
21+// Files the app itself executes. The others are entry points for running the
22+// model from the numbl CLI, and editing them changes nothing in the browser.
23+export const EXECUTED_BY_APP = new Set([
24+ 'triangle_geometry.m',
25+ 'triangle_fem.m',
26+ 'triangle_basis.m',
27+ 'triangle_strikes.m',
28+ 'triangle_rng.m',
29+ 'triangle_radiate.m',
30+ 'triangle_output.m',
31+]);
32+
33+// Cheap content hash (FNV-1a) so the worker can key its caches on the code.
34+export function hashMFiles(files: MFile[]): string {
35+ let h = 0x811c9dc5;
36+ for (const f of files) {
37+ const s = f.name + ' ' + f.source + '';
38+ for (let i = 0; i < s.length; i++) {
39+ h ^= s.charCodeAt(i);
40+ h = Math.imul(h, 0x01000193);
41+ }
42+ }
43+ return (h >>> 0).toString(36);
44+}
webapp/src/params.tsadded+259−0View file
@@ -0,0 +1,259 @@
1+// Parameter definitions for the triangle model, mirroring triangle_defaults.m.
2+// Each parameter is assigned a tier that determines how much recomputation a
3+// change requires:
4+// solve -> FEM assembly + eigensolve + strikes + radiation + output
5+// radiate -> strikes + radiation + output only (the cached mode basis is
6+// reused; the eigensolve is the expensive part of this model)
7+// output -> output stage only
8+//
9+// Note that the render duration is not a parameter: it is
10+// (strike_count - 1) / strike_rate_hz + tail_s.
11+
12+export interface TriangleParams {
13+ // Rod geometry
14+ side_length_m: number;
15+ gap_m: number;
16+ rod_diameter_m: number;
17+ // Rod material
18+ density: number;
19+ youngs_modulus: number;
20+ poisson_ratio: number;
21+ stiffness_trim: number;
22+ // Discretization / mode retention
23+ n_elements: number;
24+ min_mode_hz: number;
25+ max_mode_hz: number;
26+ // Damping
27+ loss_factor: number;
28+ loss_slope: number;
29+ radiation_loss: number;
30+ // Beater
31+ beater_radius_m: number;
32+ beater_mass_kg: number;
33+ beater_youngs: number;
34+ restitution: number;
35+ // Strike / roll gesture
36+ strike_velocity: number;
37+ strike_pos: number;
38+ strike_angle_deg: number;
39+ strike_count: number;
40+ strike_rate_hz: number;
41+ roll_pos2: number;
42+ timing_jitter_ms: number;
43+ velocity_spread: number;
44+ seed: number;
45+ // Radiation / listener
46+ direct_gain: number;
47+ listener_x_m: number;
48+ listener_y_m: number;
49+ listener_z_m: number;
50+ speed_of_sound: number;
51+ air_density: number;
52+ // Rendering
53+ sample_rate: number;
54+ oversample: number;
55+ tail_s: number;
56+ pre_delay_ms: number;
57+ normalize: boolean;
58+ peak_dbfs: number;
59+ full_scale_pa: number;
60+ bit_depth: number;
61+}
62+
63+export const DEFAULTS: TriangleParams = {
64+ side_length_m: 0.152,
65+ gap_m: 0.015,
66+ rod_diameter_m: 0.008,
67+ density: 7850.0,
68+ youngs_modulus: 2.0e11,
69+ poisson_ratio: 0.3,
70+ stiffness_trim: 1.0,
71+ n_elements: 60,
72+ min_mode_hz: 20.0,
73+ max_mode_hz: 22000.0,
74+ loss_factor: 6.0e-5,
75+ loss_slope: 4.0e-10,
76+ radiation_loss: 1.0,
77+ beater_radius_m: 0.0045,
78+ beater_mass_kg: 0.022,
79+ beater_youngs: 2.0e11,
80+ restitution: 0.6,
81+ strike_velocity: 1.5,
82+ strike_pos: 0.17,
83+ strike_angle_deg: 20.0,
84+ strike_count: 1,
85+ strike_rate_hz: 12.0,
86+ roll_pos2: 0.4,
87+ timing_jitter_ms: 6.0,
88+ velocity_spread: 0.15,
89+ seed: 1,
90+ direct_gain: 1.0,
91+ listener_x_m: 0.2,
92+ listener_y_m: 0.1,
93+ listener_z_m: 0.4,
94+ speed_of_sound: 343.0,
95+ air_density: 1.204,
96+ sample_rate: 48000,
97+ oversample: 2,
98+ tail_s: 5.0,
99+ pre_delay_ms: 3.0,
100+ normalize: true,
101+ peak_dbfs: -3.0,
102+ full_scale_pa: 2.0,
103+ bit_depth: 24,
104+};
105+
106+export type Tier = 'solve' | 'radiate' | 'output';
107+
108+// Keys whose change invalidates the cached mode basis (the eigensolve). The
109+// strike process, the listener, and everything downstream are deliberately
110+// absent: striking somewhere else or moving the microphone reuses the modes.
111+export const SOLVE_KEYS: (keyof TriangleParams)[] = [
112+ 'side_length_m', 'gap_m', 'rod_diameter_m',
113+ 'density', 'youngs_modulus', 'poisson_ratio', 'stiffness_trim',
114+ 'n_elements', 'min_mode_hz', 'max_mode_hz',
115+ 'loss_factor', 'loss_slope', 'radiation_loss',
116+ 'speed_of_sound', 'air_density',
117+];
118+
119+export const RADIATE_KEYS: (keyof TriangleParams)[] = [
120+ 'beater_radius_m', 'beater_mass_kg', 'beater_youngs', 'restitution',
121+ 'strike_velocity', 'strike_pos', 'strike_angle_deg',
122+ 'strike_count', 'strike_rate_hz', 'roll_pos2',
123+ 'timing_jitter_ms', 'velocity_spread', 'seed',
124+ 'direct_gain', 'listener_x_m', 'listener_y_m', 'listener_z_m',
125+ 'sample_rate', 'tail_s',
126+];
127+
128+export const OUTPUT_KEYS: (keyof TriangleParams)[] = [
129+ 'pre_delay_ms', 'normalize', 'peak_dbfs', 'full_scale_pa',
130+];
131+
132+export interface ParamSpec {
133+ key: keyof TriangleParams;
134+ label: string;
135+ unit?: string;
136+ min: number;
137+ max: number;
138+ step?: number;
139+ log?: boolean;
140+ integer?: boolean;
141+ help?: string;
142+}
143+
144+export interface ParamGroup {
145+ name: string;
146+ specs: ParamSpec[];
147+}
148+
149+export const PARAM_GROUPS: ParamGroup[] = [
150+ {
151+ name: 'Strike',
152+ specs: [
153+ { key: 'strike_velocity', label: 'Strike velocity', unit: 'm/s', min: 0.1, max: 6.0, step: 0.05, help: 'Beater speed at contact. Dynamics and timbre together: a faster strike gives a shorter Hertz contact and a brighter attack.' },
154+ { key: 'strike_pos', label: 'Strike position', min: 0.01, max: 0.99, step: 0.005, help: 'Arclength fraction along the rod, 0 at the free end on the base. ~0.17 is the middle of the base; a node of a mode silences that mode.' },
155+ { key: 'strike_angle_deg', label: 'Strike angle', unit: '°', min: 0, max: 90, step: 1, help: '0 strikes in the plane of the triangle, 90 normal to it. The two mode families are excited in this proportion.' },
156+ ],
157+ },
158+ {
159+ name: 'Roll',
160+ specs: [
161+ { key: 'strike_count', label: 'Strikes', min: 1, max: 64, step: 1, integer: true, help: '1 is a single strike; more gives a roll, the beater rattling in a corner.' },
162+ { key: 'strike_rate_hz', label: 'Roll rate', unit: 'Hz', min: 2, max: 30, step: 0.5, help: 'Hits per second during a roll.' },
163+ { key: 'roll_pos2', label: 'Second contact', min: 0.01, max: 0.99, step: 0.005, help: 'A roll alternates between the strike position and this one (the other leg of the corner).' },
164+ { key: 'timing_jitter_ms', label: 'Timing jitter', unit: 'ms', min: 0, max: 20, step: 0.5, help: 'Gaussian jitter on each hit after the first. 0 sounds mechanical.' },
165+ { key: 'velocity_spread', label: 'Velocity spread', min: 0, max: 0.6, step: 0.01, help: 'Lognormal spread of hit velocities during a roll.' },
166+ { key: 'seed', label: 'Random seed', min: 1, max: 9999, step: 1, integer: true, help: 'Any change redraws the roll timing and velocities; the same seed always reproduces them.' },
167+ ],
168+ },
169+ {
170+ name: 'Beater',
171+ specs: [
172+ { key: 'beater_radius_m', label: 'Tip radius', unit: 'm', min: 0.001, max: 0.012, step: 0.0005, help: 'Radius of the beater tip; enters the Hertz contact stiffness.' },
173+ { key: 'beater_mass_kg', label: 'Beater mass', unit: 'kg', min: 0.002, max: 0.12, log: true, help: 'A heavy beater lengthens the contact and darkens the strike; a light one brightens it.' },
174+ { key: 'beater_youngs', label: 'Beater modulus', unit: 'Pa', min: 1e9, max: 4e11, log: true, help: 'Steel ~2e11. Lower it toward wood (~1e10) for a soft, dark contact.' },
175+ { key: 'restitution', label: 'Restitution', min: 0.05, max: 0.95, step: 0.01, help: 'Coefficient of restitution; scales the momentum each strike transfers.' },
176+ ],
177+ },
178+ {
179+ name: 'Rod geometry',
180+ specs: [
181+ { key: 'side_length_m', label: 'Side length', unit: 'm', min: 0.08, max: 0.3, step: 0.001, help: '0.152 is a 6-inch triangle. Mode frequencies scale roughly as 1/L².' },
182+ { key: 'gap_m', label: 'Corner gap', unit: 'm', min: 0.004, max: 0.05, step: 0.001, help: 'Arclength removed at the open corner. The opening is what frees the ends and detunes the near-degenerate pairs.' },
183+ { key: 'rod_diameter_m', label: 'Rod diameter', unit: 'm', min: 0.004, max: 0.016, step: 0.0005, help: 'Bending frequencies scale as d; a thicker rod is also a louder radiator.' },
184+ ],
185+ },
186+ {
187+ name: 'Rod material',
188+ specs: [
189+ { key: 'density', label: 'Density', unit: 'kg/m³', min: 2000, max: 20000, log: true, help: 'Steel 7850, brass ~8500, aluminum 2700.' },
190+ { key: 'youngs_modulus', label: "Young's modulus", unit: 'Pa', min: 5e10, max: 4e11, log: true, help: 'Steel 2e11. Frequencies scale as sqrt(E).' },
191+ { key: 'poisson_ratio', label: 'Poisson ratio', min: 0.0, max: 0.45, step: 0.01, help: 'Sets the shear modulus and with it the torsional mode family.' },
192+ { key: 'stiffness_trim', label: 'Stiffness trim', min: 0.4, max: 2.0, step: 0.01, help: 'Overall pitch trim for departures of a real triangle from the ideal geometry.' },
193+ ],
194+ },
195+ {
196+ name: 'Damping',
197+ specs: [
198+ { key: 'loss_factor', label: 'Loss factor η₀', min: 1e-6, max: 2e-3, log: true, help: 'Frequency-independent material loss. Steel is tiny, which is why a triangle rings for seconds.' },
199+ { key: 'loss_slope', label: 'Loss slope η₁', unit: 's', min: 1e-11, max: 1e-8, log: true, help: 'Extra high-frequency damping; controls how fast the shimmer fades into the ring.' },
200+ { key: 'radiation_loss', label: 'Radiation loss', min: 0.0, max: 3.0, step: 0.05, help: 'Scales the acoustic radiation damping estimated from the element dipole powers.' },
201+ ],
202+ },
203+ {
204+ name: 'Listener',
205+ specs: [
206+ { key: 'listener_x_m', label: 'Listener x', unit: 'm', min: -2.0, max: 2.0, step: 0.005, help: 'In the plane of the triangle. Radiation-only: no re-solve needed.' },
207+ { key: 'listener_y_m', label: 'Listener y', unit: 'm', min: -2.0, max: 2.0, step: 0.005 },
208+ { key: 'listener_z_m', label: 'Listener z', unit: 'm', min: -2.0, max: 2.0, step: 0.005, help: 'Normal to the plane: where the out-of-plane modes radiate best.' },
209+ { key: 'direct_gain', label: 'Direct click gain', min: 0.0, max: 3.0, step: 0.05, help: 'Broadband local response at the contact, standing in for the truncated high modes. 0 leaves only the resolved modes.' },
210+ ],
211+ },
212+ {
213+ name: 'Output',
214+ specs: [
215+ { key: 'tail_s', label: 'Tail', unit: 's', min: 0.5, max: 12.0, step: 0.1, help: 'Ring time rendered after the last strike. A triangle needs seconds.' },
216+ { key: 'pre_delay_ms', label: 'Pre-delay', unit: 'ms', min: 0, max: 50, step: 0.5 },
217+ { key: 'peak_dbfs', label: 'Peak level', unit: 'dBFS', min: -24, max: 0, step: 0.5, help: 'Used when normalize is on.' },
218+ { key: 'full_scale_pa', label: 'Full scale', unit: 'Pa', min: 0.1, max: 20, step: 0.1, help: 'Used when normalize is off: this many pascals maps to 0 dBFS.' },
219+ ],
220+ },
221+ {
222+ name: 'Advanced / quality',
223+ specs: [
224+ { key: 'n_elements', label: 'Elements', min: 24, max: 120, step: 1, integer: true, help: 'Finite elements along the rod. More converges the high modes; the eigensolve cost grows as its cube.' },
225+ { key: 'min_mode_hz', label: 'Min mode freq', unit: 'Hz', min: 20, max: 500, step: 10 },
226+ { key: 'max_mode_hz', label: 'Max mode freq', unit: 'Hz', min: 2000, max: 24000, step: 500 },
227+ { key: 'oversample', label: 'Oversample', min: 1, max: 4, step: 1, integer: true, help: 'Synthesis rate multiple; 2 for render, 1 for preview.' },
228+ { key: 'speed_of_sound', label: 'Speed of sound', unit: 'm/s', min: 300, max: 400, step: 1 },
229+ { key: 'air_density', label: 'Air density', unit: 'kg/m³', min: 0.9, max: 1.5, step: 0.01 },
230+ ],
231+ },
232+];
233+
234+const STORAGE_KEY = 'triangle-webapp-params-v1';
235+
236+export function loadParams(): TriangleParams {
237+ try {
238+ const raw = localStorage.getItem(STORAGE_KEY);
239+ if (!raw) return { ...DEFAULTS };
240+ const parsed = JSON.parse(raw);
241+ const out = { ...DEFAULTS };
242+ for (const k of Object.keys(DEFAULTS) as (keyof TriangleParams)[]) {
243+ if (typeof parsed[k] === typeof DEFAULTS[k]) {
244+ (out as Record<string, unknown>)[k] = parsed[k];
245+ }
246+ }
247+ return out;
248+ } catch {
249+ return { ...DEFAULTS };
250+ }
251+}
252+
253+export function saveParams(p: TriangleParams): void {
254+ try {
255+ localStorage.setItem(STORAGE_KEY, JSON.stringify(p));
256+ } catch {
257+ // storage unavailable; persistence is a convenience only
258+ }
259+}
webapp/src/selftest.tsadded+25−0View file
@@ -0,0 +1,25 @@
1+// Runs the MATLAB validation script (../matlab/triangle_validate.m) through
2+// numbl in Node: closed-form identities, exact reproducibility of the seeded
3+// roll process, and the pinned end-to-end diagnostics.
4+// Run with: npm run selftest
5+
6+import { readFileSync, readdirSync } from 'fs';
7+import { join, dirname } from 'path';
8+import { fileURLToPath } from 'url';
9+import { executeCode } from 'numbl';
10+
11+const matlabDir = join(dirname(fileURLToPath(import.meta.url)), '..', '..', 'matlab');
12+const workspaceFiles = readdirSync(matlabDir)
13+ .filter((f) => f.endsWith('.m'))
14+ .map((f) => ({ name: f, source: readFileSync(join(matlabDir, f), 'utf8') }));
15+
16+const main = workspaceFiles.find((f) => f.name === 'triangle_validate.m')!;
17+const t0 = performance.now();
18+const result = executeCode(main.source, { displayResults: false }, workspaceFiles, 'triangle_validate.m');
19+const text = result.output.join('');
20+process.stdout.write(text);
21+console.log(`total time : ${((performance.now() - t0) / 1000).toFixed(2)} s`);
22+
23+if (text.includes('MISMATCH') || text.includes('FAILED')) {
24+ process.exit(1);
25+}
webapp/src/style.cssadded+358−0View file
@@ -0,0 +1,358 @@
1+:root {
2+ --bg: #14161a;
3+ --panel: #1d2026;
4+ --panel2: #23272f;
5+ --border: #31363f;
6+ --text: #d6dae2;
7+ --muted: #8b93a1;
8+ --accent: #e8a33d;
9+ --accent2: #4dd0e1;
10+ --good: #6fbf73;
11+ --warn: #e57373;
12+ font-family: system-ui, -apple-system, 'Segoe UI', Roboto, sans-serif;
13+}
14+
15+* { box-sizing: border-box; }
16+
17+html, body, #root {
18+ margin: 0;
19+ height: 100%;
20+ background: var(--bg);
21+ color: var(--text);
22+}
23+
24+.app {
25+ display: flex;
26+ flex-direction: column;
27+ height: 100vh;
28+}
29+
30+.topbar {
31+ display: flex;
32+ align-items: center;
33+ gap: 10px;
34+ padding: 8px 14px;
35+ background: var(--panel);
36+ border-bottom: 1px solid var(--border);
37+ flex-wrap: wrap;
38+}
39+
40+.topbar h1 {
41+ font-size: 16px;
42+ margin: 0 12px 0 0;
43+ font-weight: 600;
44+ color: var(--accent);
45+ white-space: nowrap;
46+}
47+
48+button {
49+ background: var(--panel2);
50+ color: var(--text);
51+ border: 1px solid var(--border);
52+ border-radius: 6px;
53+ padding: 6px 14px;
54+ font-size: 13px;
55+ cursor: pointer;
56+}
57+button:hover { border-color: var(--muted); }
58+button:disabled { opacity: 0.45; cursor: default; }
59+button.primary {
60+ background: var(--accent);
61+ color: #1a1204;
62+ border-color: var(--accent);
63+ font-weight: 600;
64+}
65+button.play {
66+ background: var(--good);
67+ color: #06210a;
68+ border-color: var(--good);
69+ font-weight: 600;
70+}
71+button.play.preview {
72+ background: transparent;
73+ color: var(--accent);
74+ border-color: var(--accent);
75+}
76+button.loop.active {
77+ background: var(--accent2);
78+ color: #062a2f;
79+ border-color: var(--accent2);
80+ font-weight: 600;
81+}
82+button.code { font-family: ui-monospace, SFMono-Regular, Menlo, monospace; }
83+button.code.active { border-color: var(--accent2); color: var(--accent2); }
84+
85+.status {
86+ margin-left: auto;
87+ font-size: 12px;
88+ color: var(--muted);
89+ display: flex;
90+ align-items: center;
91+ gap: 8px;
92+ min-width: 180px;
93+ justify-content: flex-end;
94+}
95+
96+.badge {
97+ font-size: 11px;
98+ border-radius: 10px;
99+ padding: 2px 8px;
100+ border: 1px solid var(--border);
101+}
102+.badge.preview { color: var(--accent); border-color: var(--accent); }
103+.badge.full { color: var(--good); border-color: var(--good); }
104+.badge.err { color: var(--warn); border-color: var(--warn); }
105+
106+.progressbar {
107+ width: 90px;
108+ height: 6px;
109+ background: var(--panel2);
110+ border-radius: 3px;
111+ overflow: hidden;
112+}
113+.progressbar > div {
114+ height: 100%;
115+ background: var(--accent2);
116+ transition: width 0.1s linear;
117+}
118+
119+.main {
120+ display: flex;
121+ flex: 1;
122+ min-height: 0;
123+}
124+
125+.sidebar {
126+ width: 340px;
127+ min-width: 340px;
128+ overflow-y: auto;
129+ background: var(--panel);
130+ border-right: 1px solid var(--border);
131+ padding: 8px;
132+}
133+
134+.content {
135+ flex: 1;
136+ overflow-y: auto;
137+ padding: 10px;
138+ display: grid;
139+ grid-template-columns: 1fr 1fr;
140+ gap: 10px;
141+ align-content: start;
142+}
143+
144+.card {
145+ background: var(--panel);
146+ border: 1px solid var(--border);
147+ border-radius: 8px;
148+ padding: 8px 10px 10px;
149+ min-width: 0;
150+}
151+.card.wide { grid-column: 1 / -1; }
152+.card h2 {
153+ font-size: 12px;
154+ text-transform: uppercase;
155+ letter-spacing: 0.06em;
156+ color: var(--muted);
157+ margin: 0 0 6px;
158+ font-weight: 600;
159+}
160+.card canvas { display: block; width: 100%; }
161+
162+.scene3d { height: 380px; position: relative; }
163+.scene3d canvas { height: 100%; border-radius: 6px; }
164+.scene-hint {
165+ position: absolute;
166+ bottom: 10px;
167+ right: 14px;
168+ font-size: 11px;
169+ color: var(--muted);
170+ pointer-events: none;
171+}
172+
173+/* Parameter panel */
174+details.group {
175+ border: 1px solid var(--border);
176+ border-radius: 6px;
177+ margin-bottom: 6px;
178+ background: var(--panel2);
179+}
180+details.group summary {
181+ cursor: pointer;
182+ padding: 6px 10px;
183+ font-size: 13px;
184+ font-weight: 600;
185+ user-select: none;
186+}
187+details.group[open] summary { border-bottom: 1px solid var(--border); }
188+.group-body { padding: 6px 10px 8px; }
189+
190+.param-row { margin-bottom: 8px; }
191+.param-label {
192+ display: flex;
193+ justify-content: space-between;
194+ font-size: 12px;
195+ margin-bottom: 2px;
196+ align-items: baseline;
197+}
198+.param-label .name { color: var(--text); }
199+.param-controls { display: flex; gap: 8px; align-items: center; }
200+.param-controls input[type='range'] { flex: 1; accent-color: var(--accent2); }
201+.param-controls input[type='number'] {
202+ width: 88px;
203+ background: var(--bg);
204+ color: var(--text);
205+ border: 1px solid var(--border);
206+ border-radius: 4px;
207+ padding: 3px 5px;
208+ font-size: 12px;
209+}
210+.param-controls select, .param-controls input[type='checkbox'] {
211+ background: var(--bg);
212+ color: var(--text);
213+ border: 1px solid var(--border);
214+ border-radius: 4px;
215+ font-size: 12px;
216+}
217+.param-help { font-size: 11px; color: var(--muted); margin-top: 1px; }
218+
219+/* Diagnostics */
220+.diag table { border-collapse: collapse; font-size: 12px; width: 100%; }
221+.diag td, .diag th {
222+ padding: 2px 8px 2px 0;
223+ text-align: left;
224+ color: var(--text);
225+ font-weight: 400;
226+}
227+.diag th { color: var(--muted); white-space: nowrap; }
228+.diag .partials { font-variant-numeric: tabular-nums; }
229+.warning { color: var(--warn); font-size: 12px; margin-top: 4px; }
230+.cachenote { color: var(--muted); font-size: 11px; margin-top: 6px; }
231+
232+.mode-list { display: flex; gap: 6px; flex-wrap: wrap; margin-bottom: 6px; }
233+.mode-list button { padding: 3px 10px; font-size: 12px; }
234+.mode-list button.active { border-color: var(--accent2); color: var(--accent2); }
235+
236+/* Code panel */
237+.modal-backdrop {
238+ position: fixed;
239+ inset: 0;
240+ background: rgba(8, 10, 13, 0.72);
241+ display: flex;
242+ align-items: center;
243+ justify-content: center;
244+ z-index: 50;
245+}
246+
247+.code-panel {
248+ width: min(1180px, 94vw);
249+ height: min(860px, 92vh);
250+ background: var(--panel);
251+ border: 1px solid var(--border);
252+ border-radius: 10px;
253+ display: flex;
254+ flex-direction: column;
255+ overflow: hidden;
256+}
257+
258+.code-header {
259+ display: flex;
260+ align-items: center;
261+ gap: 10px;
262+ padding: 9px 12px;
263+ border-bottom: 1px solid var(--border);
264+}
265+.code-header h2 {
266+ margin: 0;
267+ font-size: 13px;
268+ text-transform: uppercase;
269+ letter-spacing: 0.06em;
270+ color: var(--accent);
271+}
272+.code-header button:last-child { margin-left: 0; }
273+.code-note {
274+ font-size: 11.5px;
275+ color: var(--muted);
276+ margin-right: auto;
277+}
278+
279+.code-body { display: flex; flex: 1; min-height: 0; }
280+
281+.code-files {
282+ width: 236px;
283+ min-width: 236px;
284+ border-right: 1px solid var(--border);
285+ overflow-y: auto;
286+ padding: 6px;
287+}
288+.code-files button {
289+ display: flex;
290+ align-items: center;
291+ justify-content: space-between;
292+ gap: 6px;
293+ width: 100%;
294+ text-align: left;
295+ background: transparent;
296+ border: 1px solid transparent;
297+ border-radius: 5px;
298+ padding: 5px 8px;
299+ margin-bottom: 2px;
300+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
301+ font-size: 12px;
302+}
303+.code-files button:hover { background: var(--panel2); }
304+.code-files button.active {
305+ background: var(--panel2);
306+ border-color: var(--accent2);
307+ color: var(--accent2);
308+}
309+.code-files .fname { display: flex; align-items: center; gap: 5px; }
310+.code-files .dot {
311+ width: 6px;
312+ height: 6px;
313+ border-radius: 50%;
314+ background: var(--accent);
315+}
316+.code-files .ftag {
317+ font-size: 9.5px;
318+ color: var(--muted);
319+ text-transform: uppercase;
320+ letter-spacing: 0.05em;
321+ font-family: system-ui, sans-serif;
322+}
323+
324+.code-main { flex: 1; min-width: 0; display: flex; flex-direction: column; }
325+.code-toolbar {
326+ display: flex;
327+ align-items: center;
328+ gap: 10px;
329+ padding: 6px 10px;
330+ border-bottom: 1px solid var(--border);
331+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
332+ font-size: 12px;
333+ color: var(--text);
334+}
335+.code-hint {
336+ font-size: 11px;
337+ color: var(--muted);
338+ font-family: system-ui, sans-serif;
339+}
340+.code-editor { flex: 1; min-height: 0; overflow: hidden; }
341+
342+.code-error {
343+ border-top: 1px solid var(--warn);
344+ background: #2a1a1a;
345+ color: #ffb4a9;
346+ padding: 7px 12px;
347+ font-size: 12px;
348+ font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
349+}
350+
351+.error-card { border-color: var(--warn); }
352+.error-card h2 { color: var(--warn); }
353+.error-card div { font-size: 12px; }
354+.error-card strong { font-family: ui-monospace, SFMono-Regular, Menlo, monospace; }
355+
356+@media (max-width: 1100px) {
357+ .content { grid-template-columns: 1fr; }
358+}
webapp/src/vite-env.d.tsadded+1−0View file
@@ -0,0 +1 @@
1+/// <reference types="vite/client" />
webapp/src/worker/renderWorker.tsadded+430−0View file
@@ -0,0 +1,430 @@
1+// Web Worker that renders the triangle strike by executing the MATLAB
2+// implementation in ../matlab through numbl. The stages run as separate
3+// executeCode calls so intermediates can be cached and cheap parameter
4+// changes skip the expensive ones:
5+// solve tier -> triangle_basis (FEM assembly + eigensolve; the expensive
6+// stage, cached on geometry/material/damping params)
7+// radiate tier -> triangle_strikes + triangle_radiate (cached on the
8+// strike process, the listener, and the synthesis rate)
9+// output tier -> triangle_output (always cheap, always re-run)
10+// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).
11+
12+import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
13+import { hashMFiles, type MFile } from '../mfiles';
14+import { SOLVE_KEYS, RADIATE_KEYS, type TriangleParams } from '../params';
15+
16+export interface RenderRequest {
17+ type: 'render';
18+ id: number;
19+ params: TriangleParams;
20+ quality: 'preview' | 'full';
21+ mfiles: MFile[];
22+}
23+
24+export interface RenderResponse {
25+ type: 'result';
26+ id: number;
27+ quality: 'preview' | 'full';
28+ y: Float32Array; // final signal, [-1, 1]
29+ sampleRate: number;
30+ strikeT: Float32Array; // strike times, s
31+ strikeV: Float32Array; // strike velocities, m/s
32+ strikeS: Float32Array; // strike arclength fractions
33+ modeF: Float32Array; // mode frequencies, Hz
34+ modeZeta: Float32Array; // damping ratio of each mode
35+ modeZetaRad: Float32Array; // radiation part of the damping
36+ modeFamily: Float32Array; // 0 = in-plane, 1 = out-of-plane
37+ nodeX: Float32Array; // rod centerline nodes, m
38+ nodeY: Float32Array;
39+ nodeS: Float32Array; // arclength fraction of each node
40+ shapeUx: Float32Array; // node displacements per mode, nn x n col-major
41+ shapeUy: Float32Array;
42+ shapeUz: Float32Array;
43+ diag: {
44+ modesKept: number;
45+ modesIn: number;
46+ modesOut: number;
47+ strikeCount: number;
48+ meanV: number;
49+ maxV: number;
50+ meanTauUs: number;
51+ peakForceN: number;
52+ durationS: number;
53+ peakPa: number;
54+ peakSplDb: number;
55+ directPeakPa: number;
56+ listenerDistM: number;
57+ warnings: string[];
58+ solveMs: number; // 0 when the cached eigensolve was reused
59+ radiateMs: number; // 0 when the cached radiation was reused
60+ solveCached: boolean;
61+ radiateCached: boolean;
62+ oversample: number;
63+ };
64+}
65+
66+export interface ProgressMessage {
67+ type: 'progress';
68+ id: number;
69+ stage: 'solve' | 'radiate' | 'output';
70+ frac: number;
71+}
72+
73+// Raised by a fault in the .m code (the user can edit it), so it is reported
74+// rather than left to stall the request.
75+export interface ErrorMessage {
76+ type: 'error';
77+ id: number;
78+ message: string;
79+ file?: string;
80+ line?: number;
81+}
82+
83+// Basis fields that cross from the solve stage into the radiate stage as
84+// numbl RuntimeValues, plus their Float32 copies for the plots.
85+interface SolveEntry {
86+ vars: Record<string, RuntimeValue>;
87+ modeF: Float32Array;
88+ modeZeta: Float32Array;
89+ modeZetaRad: Float32Array;
90+ modeFamily: Float32Array;
91+ nodeX: Float32Array;
92+ nodeY: Float32Array;
93+ nodeS: Float32Array;
94+ shapeUx: Float32Array;
95+ shapeUy: Float32Array;
96+ shapeUz: Float32Array;
97+ modesIn: number;
98+ modesOut: number;
99+}
100+
101+interface RadiateEntry {
102+ p: RuntimeValue;
103+ strikeT: Float32Array;
104+ strikeV: Float32Array;
105+ strikeS: Float32Array;
106+ strikeCount: number;
107+ meanV: number;
108+ maxV: number;
109+ meanTauUs: number;
110+ peakForceN: number;
111+ durationS: number;
112+ directPeakPa: number;
113+ listenerDistM: number;
114+ listenerInside: boolean;
115+ fsHi: number;
116+ oversample: number;
117+}
118+
119+// Small LRU caches. The solve cache does not depend on the oversample rate,
120+// so a preview render right after a full one reuses the same eigensolve.
121+const solveCache = new Map<string, SolveEntry>();
122+const radiateCache = new Map<string, RadiateEntry>();
123+const SOLVE_CACHE_MAX = 2;
124+const RADIATE_CACHE_MAX = 3;
125+
126+function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
127+ const v = cache.get(key);
128+ if (v !== undefined) {
129+ cache.delete(key);
130+ cache.set(key, v);
131+ }
132+ return v;
133+}
134+
135+function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
136+ cache.delete(key);
137+ cache.set(key, value);
138+ while (cache.size > max) {
139+ const oldest = cache.keys().next().value as string;
140+ cache.delete(oldest);
141+ }
142+}
143+
144+function keyFor(params: TriangleParams, keys: readonly (keyof TriangleParams)[]): string {
145+ return keys.map((k) => String(params[k])).join('|');
146+}
147+
148+// MATLAB source lines assigning every parameter onto a struct P.
149+function pLiteral(params: TriangleParams): string {
150+ const lines: string[] = [];
151+ for (const [k, v] of Object.entries(params)) {
152+ if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
153+ else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
154+ }
155+ return lines.join('\n');
156+}
157+
158+function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
159+ return executeCode(
160+ source,
161+ { initialVariableValues: inputs, displayResults: false },
162+ mfiles,
163+ 'stage.m',
164+ );
165+}
166+
167+function asTensor(v: RuntimeValue): RuntimeTensor {
168+ const t = v as RuntimeTensor;
169+ if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
170+ return t;
171+}
172+
173+function asScalar(v: RuntimeValue): number {
174+ if (typeof v === 'number') return v;
175+ const t = v as RuntimeTensor;
176+ if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
177+ if (typeof v === 'boolean') return v ? 1 : 0;
178+ throw new Error('expected scalar result from numbl stage');
179+}
180+
181+function toF32(v: RuntimeValue): Float32Array {
182+ const t = asTensor(v);
183+ const out = new Float32Array(t.data.length);
184+ for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
185+ return out;
186+}
187+
188+// The basis fields handed from the solve stage into the radiate stage, in the
189+// order the reassembly script expects.
190+const BASIS_FIELDS = [
191+ 'f', 'omega', 'zeta', 'zetaRad', 'family',
192+ 'nodeX', 'nodeY', 'sfrac', 'totalLen',
193+ 'NUx', 'NUy', 'NUz',
194+ 'elemX', 'elemY', 'elemTanX', 'elemTanY', 'elemLen',
195+ 'EUpx', 'EUpy', 'EUpz',
196+];
197+
198+const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
199+
200+self.onmessage = (ev: MessageEvent<RenderRequest>) => {
201+ const msg = ev.data;
202+ if (msg.type !== 'render') return;
203+
204+ // numbl reports a file it could not parse through console.warn and then
205+ // skips it, so the failure surfaces at the call site as "Undefined
206+ // function". Capture the warning to report the actual syntax error instead.
207+ const warnings: string[] = [];
208+ const origWarn = console.warn;
209+ console.warn = (...args: unknown[]) => {
210+ warnings.push(args.map(String).join(' '));
211+ };
212+
213+ try {
214+ render(msg);
215+ } catch (err) {
216+ // Most likely a fault in the (user-editable) .m code. Report it so the
217+ // request completes and the message reaches the code panel.
218+ const e = err as { message?: string; file?: string; line?: number };
219+ let m: ErrorMessage = {
220+ type: 'error',
221+ id: msg.id,
222+ message: e?.message ? String(e.message) : String(err),
223+ file: typeof e?.file === 'string' ? e.file : undefined,
224+ line: typeof e?.line === 'number' ? e.line : undefined,
225+ };
226+ for (const w of warnings) {
227+ const hit = SKIP_WARNING.exec(w);
228+ if (hit) {
229+ m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
230+ break;
231+ }
232+ }
233+ (self as unknown as Worker).postMessage(m);
234+ } finally {
235+ console.warn = origWarn;
236+ }
237+};
238+
239+function render(msg: RenderRequest): void {
240+ const { id, params, mfiles } = msg;
241+ let quality = msg.quality;
242+
243+ // The .m sources are part of the cache identity: editing the code must
244+ // invalidate everything computed from the previous version.
245+ const codeKey = hashMFiles(mfiles);
246+ const solveKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);
247+
248+ const fullOversample = params.oversample;
249+ const radiateParamsKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);
250+ if (quality === 'preview' && radiateCache.has(radiateParamsKey + `|os=${fullOversample}`)) {
251+ quality = 'full'; // full-rate radiation already cached
252+ }
253+ const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample;
254+ const fsHi = params.sample_rate * oversample;
255+
256+ const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
257+
258+ // ---- solve stage: FEM assembly and eigensolve ----
259+ let solveMs = 0;
260+ let solveCached = true;
261+ let solveEntry = lruGet(solveCache, solveKey);
262+ if (!solveEntry) {
263+ solveCached = false;
264+ post({ type: 'progress', id, stage: 'solve', frac: 0.15 });
265+ const t0 = performance.now();
266+ const r = runStage(mfiles, `
267+${pLiteral(params)}
268+basis = triangle_basis(P);
269+${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}
270+modesIn = sum(basis.family == 0);
271+modesOut = sum(basis.family == 1);
272+`);
273+ solveMs = performance.now() - t0;
274+
275+ const vars: Record<string, RuntimeValue> = {};
276+ for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];
277+ solveEntry = {
278+ vars,
279+ modeF: toF32(r.variableValues.b_f),
280+ modeZeta: toF32(r.variableValues.b_zeta),
281+ modeZetaRad: toF32(r.variableValues.b_zetaRad),
282+ modeFamily: toF32(r.variableValues.b_family),
283+ nodeX: toF32(r.variableValues.b_nodeX),
284+ nodeY: toF32(r.variableValues.b_nodeY),
285+ nodeS: toF32(r.variableValues.b_sfrac),
286+ shapeUx: toF32(r.variableValues.b_NUx),
287+ shapeUy: toF32(r.variableValues.b_NUy),
288+ shapeUz: toF32(r.variableValues.b_NUz),
289+ modesIn: asScalar(r.variableValues.modesIn),
290+ modesOut: asScalar(r.variableValues.modesOut),
291+ };
292+ lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
293+ }
294+
295+ // ---- radiate stage: strike process + modal synthesis ----
296+ const radiateKey = radiateParamsKey + `|os=${oversample}`;
297+ let radiateMs = 0;
298+ let radiateCached = true;
299+ let radEntry = lruGet(radiateCache, radiateKey);
300+ if (!radEntry) {
301+ radiateCached = false;
302+ post({ type: 'progress', id, stage: 'radiate', frac: 0.55 });
303+ const t0 = performance.now();
304+ const r = runStage(mfiles, `
305+${pLiteral(params)}
306+${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}
307+basis.n = length(basis.f);
308+[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P);
309+[p, radinfo] = triangle_radiate(P, basis, tstr, vstr, jstr, taustr, fstr, sstr, ${fsHi});
310+strikeCount = sinfo.count;
311+meanV = sinfo.meanV;
312+maxV = sinfo.maxV;
313+meanTauUs = sinfo.meanTauUs;
314+peakForceN = sinfo.peakForceN;
315+durationS = sinfo.durationS;
316+directPeakPa = radinfo.directPeakPa;
317+listenerDistM = radinfo.listenerDistM;
318+listenerInside = radinfo.listenerInside;
319+`, solveEntry.vars);
320+ radiateMs = performance.now() - t0;
321+ radEntry = {
322+ p: r.variableValues.p,
323+ strikeT: toF32(r.variableValues.tstr),
324+ strikeV: toF32(r.variableValues.vstr),
325+ strikeS: toF32(r.variableValues.sstr),
326+ strikeCount: asScalar(r.variableValues.strikeCount),
327+ meanV: asScalar(r.variableValues.meanV),
328+ maxV: asScalar(r.variableValues.maxV),
329+ meanTauUs: asScalar(r.variableValues.meanTauUs),
330+ peakForceN: asScalar(r.variableValues.peakForceN),
331+ durationS: asScalar(r.variableValues.durationS),
332+ directPeakPa: asScalar(r.variableValues.directPeakPa),
333+ listenerDistM: asScalar(r.variableValues.listenerDistM),
334+ listenerInside: asScalar(r.variableValues.listenerInside) !== 0,
335+ fsHi,
336+ oversample,
337+ };
338+ lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
339+ }
340+
341+ // ---- output stage ----
342+ post({ type: 'progress', id, stage: 'output', frac: 0.9 });
343+ const rOut = runStage(mfiles, `
344+${pLiteral(params)}
345+[y, peakPa] = triangle_output(P, p, ${radEntry.oversample});
346+`, { p: radEntry.p });
347+
348+ const y32 = toF32(rOut.variableValues.y);
349+ const peakPa = asScalar(rOut.variableValues.peakPa);
350+
351+ const warnings: string[] = [];
352+ if (radEntry.listenerInside) {
353+ warnings.push(
354+ 'the listener is within 2 cm of the rod; the dipole field is evaluated at a clamped minimum distance',
355+ );
356+ }
357+ if (!params.normalize) {
358+ let peak = 0;
359+ for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
360+ if (peak > 1.0) {
361+ warnings.push(
362+ `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,
363+ );
364+ }
365+ }
366+ const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];
367+ if (topMode < 0.4 * params.max_mode_hz) {
368+ warnings.push(
369+ `the highest retained mode is only ${topMode.toFixed(0)} Hz; the shimmer above it falls entirely on the direct click`,
370+ );
371+ }
372+
373+ const resp: RenderResponse = {
374+ type: 'result',
375+ id,
376+ quality,
377+ y: y32,
378+ sampleRate: params.sample_rate,
379+ strikeT: radEntry.strikeT.slice(),
380+ strikeV: radEntry.strikeV.slice(),
381+ strikeS: radEntry.strikeS.slice(),
382+ modeF: solveEntry.modeF.slice(),
383+ modeZeta: solveEntry.modeZeta.slice(),
384+ modeZetaRad: solveEntry.modeZetaRad.slice(),
385+ modeFamily: solveEntry.modeFamily.slice(),
386+ nodeX: solveEntry.nodeX.slice(),
387+ nodeY: solveEntry.nodeY.slice(),
388+ nodeS: solveEntry.nodeS.slice(),
389+ shapeUx: solveEntry.shapeUx.slice(),
390+ shapeUy: solveEntry.shapeUy.slice(),
391+ shapeUz: solveEntry.shapeUz.slice(),
392+ diag: {
393+ modesKept: solveEntry.modeF.length,
394+ modesIn: solveEntry.modesIn,
395+ modesOut: solveEntry.modesOut,
396+ strikeCount: radEntry.strikeCount,
397+ meanV: radEntry.meanV,
398+ maxV: radEntry.maxV,
399+ meanTauUs: radEntry.meanTauUs,
400+ peakForceN: radEntry.peakForceN,
401+ durationS: radEntry.durationS,
402+ peakPa,
403+ peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
404+ directPeakPa: radEntry.directPeakPa,
405+ listenerDistM: radEntry.listenerDistM,
406+ warnings,
407+ solveMs,
408+ radiateMs,
409+ solveCached,
410+ radiateCached,
411+ oversample: radEntry.oversample,
412+ },
413+ };
414+ (self as unknown as Worker).postMessage(resp, [
415+ resp.y.buffer,
416+ resp.strikeT.buffer,
417+ resp.strikeV.buffer,
418+ resp.strikeS.buffer,
419+ resp.modeF.buffer,
420+ resp.modeZeta.buffer,
421+ resp.modeZetaRad.buffer,
422+ resp.modeFamily.buffer,
423+ resp.nodeX.buffer,
424+ resp.nodeY.buffer,
425+ resp.nodeS.buffer,
426+ resp.shapeUx.buffer,
427+ resp.shapeUy.buffer,
428+ resp.shapeUz.buffer,
429+ ]);
430+}
webapp/tsconfig.jsonadded+20−0View file
@@ -0,0 +1,20 @@
1+{
2+ "compilerOptions": {
3+ "target": "ES2022",
4+ "useDefineForClassFields": true,
5+ "lib": ["ES2022", "DOM", "DOM.Iterable", "WebWorker"],
6+ "module": "ESNext",
7+ "skipLibCheck": true,
8+ "moduleResolution": "bundler",
9+ "resolveJsonModule": true,
10+ "isolatedModules": true,
11+ "noEmit": true,
12+ "jsx": "react-jsx",
13+ "strict": true,
14+ "noUnusedLocals": true,
15+ "noUnusedParameters": true,
16+ "noFallthroughCasesInSwitch": true
17+ },
18+ "include": ["src"],
19+ "exclude": ["src/selftest.ts"]
20+}
webapp/vite.config.tsadded+9−0View file
@@ -0,0 +1,9 @@
1+import { defineConfig } from 'vite';
2+import react from '@vitejs/plugin-react';
3+
4+export default defineConfig({
5+ plugins: [react()],
6+ base: './',
7+ worker: { format: 'es' },
8+ server: { fs: { allow: ['..'] } }, // the .m sources live in ../matlab
9+});