Cymbal simulator: shell-of-revolution modal model with a nonlinear energy cascade, numbl MATLAB reference, and interactive web app
48 changed files+9187−0
.github/workflows/deploy.ymladded+70−0View file
@@ -0,0 +1,70 @@
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/cymbal_validate.m through numbl: the free circular plate
41+ # closed forms, the eigensolve and radiation identities, energy
42+ # conservation of the cascade, reproducibility of the seeded stroke
43+ # pattern, and the pinned end-to-end diagnostics. Fails the build on any
44+ # mismatch.
45+ - name: Validate the model
46+ working-directory: webapp
47+ run: npm run selftest
48+
49+ - name: Build
50+ working-directory: webapp
51+ run: npm run build
52+
53+ - name: Setup Pages
54+ uses: actions/configure-pages@v5
55+
56+ - name: Upload artifact
57+ uses: actions/upload-pages-artifact@v3
58+ with:
59+ path: webapp/dist
60+
61+ deploy:
62+ needs: build
63+ runs-on: ubuntu-latest
64+ environment:
65+ name: github-pages
66+ url: ${{ steps.deployment.outputs.page_url }}
67+ steps:
68+ - name: Deploy to Pages
69+ id: deployment
70+ uses: actions/deploy-pages@v4
README.mdadded+49−0View file
@@ -0,0 +1,49 @@
1+# Cymbal
2+
3+A physically-based simulation of a struck cymbal, sampled as acoustic pressure at a point in space. The instrument is a bronze shell of revolution with a raised bell and a tapered bow; a Hertz stick contact drives its bending and stretching modes, computed by a one-dimensional finite element eigensolve per nodal-diameter order, and the geometric nonlinearity that turns a hard stroke into a crash is carried as a cascade of energy up the modal frequency ladder. The plate radiates to a listener point as a sheet of normal dipoles.
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 cymbal_demo.m`, or `cymbal_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+meridian profile: bell, bow, thickness taper cymbal_profile.m
22+ -> thin-shell elements for one Fourier order cymbal_shell.m
23+ -> one eigensolve per nodal-diameter order cymbal_basis.m
24+ + radiated power -> radiation damping
25+stroke gesture -> Hertz stick contacts cymbal_strikes.m
26+ -> modal energies, and their nonlinear flow
27+ up the frequency ladder cymbal_cascade.m
28+ -> oscillator bank + dipole-sheet transfer cymbal_radiate.m
29+ -> free-field pressure at the listener
30+ -> anti-aliased decimation -> PCM cymbal_output.m
31+```
32+
33+Three points are worth stating plainly, because they are where the model commits itself.
34+
35+**The modes come from a shell of revolution, one Fourier order at a time.** A cymbal is axisymmetric, so its displacement field separates exactly in the azimuth: every mode carries `n` nodal diameters and varies as `cos(n theta)`, and the two-dimensional eigenproblem breaks into one one-dimensional problem per `n`. Each is discretized along the meridian with two-node Love-Kirchhoff shell elements, linear in the meridional and circumferential displacements and cubic Hermite in the normal one, integrated by Gauss quadrature against the tabulated profile rather than frozen into conical frusta (`cymbal_shell.m`). Rotary inertia is kept, which shifts the retained frequencies by well under a percent but is what makes the mass matrix positive definite, so each order reduces through the Cholesky factor of the mass matrix and is diagonalized by svd. The curvature of the bow is what earns this formulation its keep: on a flat disc, bending and stretching are independent, and it is the meridian curvature that couples them and pulls the partials of a cymbal away from the Chladni frequencies of a plate. Shear deformation is neglected, and the mounting is taken as free-free with a point dashpot at the hole for the felt washer, which reaches the `n = 0` and `n = 1` modes and almost nothing else. Only the `cos` member of each degenerate pair is retained: that is exact for the model as posed, since every contact sits on the meridian `theta = 0` where the `sin` member has a node. A real cymbal is not perfectly round, and the resulting slow beating of each split pair is outside this model.
36+
37+**The nonlinearity is carried as energy, not as a coupled integration.** A cymbal struck hard is not a linear resonator. Its transverse displacement is comparable to the plate thickness, the von Karman stretching of the mid-surface couples the modes, and energy injected low by the stick migrates upward over tens of milliseconds, which is what the ear hears as a crash opening into shimmer after its attack. Integrating that directly means a dense cubic coupling tensor over hundreds of modes, so what is integrated instead is the energy of each mode (`cymbal_cascade.m`): each mode has an up-target a fixed frequency ratio above it and a down-target the same ratio below, and energy performs a biased random walk between them at a rate `cascade_gain * om * (w/t)^2`, proportional to the mode's own squared amplitude-to-thickness ratio because that is the order of the geometric nonlinearity. Steps are taken in log frequency rather than in mode index on purpose, so the swell does not depend on how densely the band happens to be resolved. Because the flux out of a mode goes as the square of its energy, a nearly empty mode passes nothing on, the front of the cascade advances rather than the whole band draining at once, and the swell stops on its own. What this buys is the instrument's central behaviour: the same cymbal pings when tapped and crashes when hit, with no velocity layers anywhere in the model. What it gives up is phase. Amplitude comes from the cascade and phase from a linear impulse response carried alongside it, so a stroke re-phases the modes it finds quiet and a coherent attack survives, but the intermodulation products of a real cascade do not. `cascade_gain` is the one fitted number in the model; setting it to zero leaves an exactly linear plate, which sounds like a small gong.
38+
39+**Radiation is a dipole sheet, computed two ways because two regimes are needed.** Each surface element carries a dipole along the local normal of volume `dS * ldip`, where `ldip` is half the shorter of the flexural and acoustic wavelengths, capped at the plate radius: the distance over which the air short-circuits from one face to the other. That length is the one phenomenological quantity in the radiation model, and below the plate's coincidence frequency, about 18 kHz for 1 mm bronze, it is the flexural wavelength that sets it, which is why a cymbal is such a poor radiator in the low register and improves steadily as the cascade fills the top. For the listener, who at half a metre from a 40 cm plate is nowhere near the far field, the element dipoles are summed coherently over the azimuth with their true distances and the near-field term kept (`cymbal_radiate.m`). For the radiated power, which sets the radiation damping and is defined only in the far field, the azimuthal integral of the same kernel collapses onto Bessel functions and the power integral reduces to one quadrature in `cos(theta)` (`cymbal_basis.m`); the two agree to a few parts in a thousand, which the validation checks. Summing element powers incoherently, as the sibling `triangle` does for a rod, would not do here: for a continuous surface that estimate halves when the elements do. The multipole structure comes out of the azimuthal sum on its own. A mode with `n` nodal diameters is a `2^n`-pole, so the high-order modes are weak until their wavelength is short, and on the symmetry axis every one of them cancels exactly.
40+
41+## Validation
42+
43+`matlab/cymbal_validate.m` (also `npm run selftest` in `webapp/`) checks the shell elements in their flat limit against the closed-form free circular plate, whose frequency equation is the vanishing of a two-by-two determinant in `J_n` and `I_n` under zero radial moment and zero Kelvin-Kirchhoff shear, together with the rigid-body count each Fourier order must have (two at `n = 0`, two at `n = 1`, none above); the mass-orthonormality and stiffness-diagonalization identities of the eigensolve; the Bessel far-field reduction against brute-force summation over a finely discretized dipole sheet; the exact cancellation of every mode with a nodal diameter on the symmetry axis; the Hertz `v^(-1/5)` contact-time scaling in the light-stick limit; exact energy conservation of the cascade scheme with the damping off and the top of the ladder reflecting, and its collapse onto closed-form damped exponentials when `cascade_gain` is zero; the far-field `-6 dB` per doubling of distance; exact reproducibility of the seeded stroke pattern; and pinned end-to-end diagnostics.
44+
45+As in the sibling `maraca`, `triangle`, `recorder` and `bass` repositories, 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.
46+
47+## Defaults
48+
49+A 16-inch B20 bronze crash-ride: 20.3 cm radius, a 4.5 cm bell rising 1.7 cm above a bow that itself rises 2 cm from the rim, and a thickness tapering from 1.15 mm at the bell to 0.75 mm at the edge, which comes to 1.05 kg. Solving to 48 nodal diameters retains 303 modes between 49 Hz and 8.0 kHz. The default gesture is a single crash on the edge at 3 m/s with a 30 g nylon-tipped stick held 40 degrees off the surface normal: an 88 microsecond Hertz contact, a 291 N peak force, a peak plate displacement of 150 micrometres or 0.16 thicknesses, and a peak of 13.4 Pa (116 dB SPL) at the 58 cm listener before normalization. About 6% of the struck energy leaves the top of the retained band, which is the shimmer above 8 kHz the model does not carry; raising `max_mode_hz` and `n_diameters` together recovers it at the cost of the eigensolves. Setting `strike_count` to 8 gives a repeated figure alternating between the edge and the bow, and dropping `strike_velocity` to a few tenths of a metre per second turns the same instrument back into a ping, because the cascade has nothing left to work with.
matlab/cymbal_at.madded+35−0View file
@@ -0,0 +1,35 @@
1+function [uh, wh, thk, drds, dzds, cs, ct] = cymbal_at(basis, sfrac)
2+% CYMBAL_AT Mode shapes and local geometry at an arclength fraction.
3+%
4+% sfrac addresses the meridian from the hole (0) to the rim (1); it is the
5+% coordinate the strike positions are expressed in. The meridional shape U is
6+% linear within an element, so it is interpolated linearly; the normal shape W
7+% is cubic Hermite, and is interpolated with its own shape functions from the
8+% nodal (W, dW/ds) pair, which matters because at the top of the retained band
9+% there are only about four nodes per flexural wavelength.
10+%
11+% Returns row vectors uh, wh (1 x n modes) and the local thickness, meridian
12+% direction cosines and principal curvatures at that point.
13+
14+nn = length(basis.snode);
15+S = basis.S;
16+s = min(0.999999, max(0, sfrac)) * S;
17+idx = sum(basis.snode <= s);
18+idx = min(max(idx, 1), nn - 1);
19+L = basis.snode(idx + 1) - basis.snode(idx);
20+xi = (s - basis.snode(idx)) / L;
21+
22+uh = basis.NU(idx, :) + xi * (basis.NU(idx + 1, :) - basis.NU(idx, :));
23+h1 = 1 - 3 * xi^2 + 2 * xi^3;
24+h2 = L * (xi - 2 * xi^2 + xi^3);
25+h3 = 3 * xi^2 - 2 * xi^3;
26+h4 = L * (-xi^2 + xi^3);
27+wh = h1 * basis.NW(idx, :) + h2 * basis.NT(idx, :) ...
28+ + h3 * basis.NW(idx + 1, :) + h4 * basis.NT(idx + 1, :);
29+
30+thk = basis.thknode(idx) + xi * (basis.thknode(idx + 1) - basis.thknode(idx));
31+drds = basis.drdsnode(idx) + xi * (basis.drdsnode(idx + 1) - basis.drdsnode(idx));
32+dzds = basis.dzdsnode(idx) + xi * (basis.dzdsnode(idx + 1) - basis.dzdsnode(idx));
33+cs = basis.csnode(idx) + xi * (basis.csnode(idx + 1) - basis.csnode(idx));
34+ct = basis.ctnode(idx) + xi * (basis.ctnode(idx + 1) - basis.ctnode(idx));
35+end
matlab/cymbal_basis.madded+241−0View file
@@ -0,0 +1,241 @@
1+function basis = cymbal_basis(P)
2+% CYMBAL_BASIS Mode set of the cymbal: one shell eigensolve per nodal-diameter
3+% order, plus damping and the radiation data the later stages need.
4+%
5+% Because the plate is a shell of revolution, the modes separate exactly in
6+% the azimuth (see cymbal_shell.m), so the two-dimensional eigenproblem breaks
7+% into n_diameters + 1 one-dimensional ones. Each is reduced to symmetric
8+% standard form through the Cholesky factor of the consistent mass matrix and
9+% diagonalized by svd rather than eig: the reduced matrix is symmetric
10+% positive semidefinite, so its singular vectors are its eigenvectors, and svd
11+% keeps them orthonormal for the closely spaced pairs this geometry produces.
12+% The eigenvectors come back mass-normalized, phi' M phi = I. Rigid-body modes
13+% (two at n = 0, two at n = 1, none above) sit at zero frequency and are
14+% dropped with everything below min_mode_hz; modes above max_mode_hz are
15+% truncated, and the broadband part they carry is reinstated by the direct
16+% term in cymbal_radiate.m.
17+%
18+% Only the cos(n theta) member of each degenerate pair is retained. That is
19+% exact for the model as posed: every contact sits on the meridian theta = 0,
20+% where the sin member has a node. A real cymbal is not perfectly round, which
21+% splits each pair into a slowly beating doublet; that splitting is outside
22+% this model, and the README says so.
23+%
24+% Damping has three parts.
25+% internal damp_rate / om + 0.5 (loss_factor + loss_slope om). The first
26+% term is a frequency-independent decay RATE, which is what a
27+% resistive path out of the plate gives and what keeps the low
28+% modes of a cymbal ringing for tens of seconds; the other two are
29+% the usual constant and frequency-proportional loss factors, and
30+% they are what makes the shimmer die well before the hum does;
31+% mount a point dashpot of mount_loss N s/m on the axial motion at the
32+% inner edge, standing in for the felt washer and the stand. With
33+% mass-normalized shapes its modal damping is c phi_axial^2 / 2 om,
34+% which falls off steeply with n because the shapes of high
35+% nodal-diameter modes vanish at the hole. This is what kills the
36+% low n = 0 gong modes a mounted cymbal does not sustain;
37+% radiation computed from the actual radiated power, below.
38+%
39+% Radiated power is evaluated in the far field, where it is defined, rather
40+% than by summing element powers as the sibling triangle does (an incoherent
41+% element sum is mesh dependent for a continuous surface: halving the element
42+% size halves the estimate). The plate is treated as a sheet of normal dipoles
43+% of volume dS * ldip, and for a cos(n theta) source on a surface of revolution
44+% the azimuthal integral of the far-field kernel is a Bessel function, so the
45+% far field is exactly
46+% p(th_o, ph_o) = -(om^2 rho0 / (4 pi c R)) e^{-ikR} A_n(th_o) cos(n ph_o)
47+% A_n = ldip sum_e r_e L_e W_e e^{i k z_e cos th_o}
48+% [ 2 pi i^n J_n(al) r'_e cos th_o
49+% - pi (i^{n+1} J_{n+1}(al) + i^{n-1} J_{n-1}(al)) z'_e sin th_o ]
50+% al = k r_e sin th_o
51+% and the power integral over the sphere collapses to one quadrature in
52+% cos th_o. The cos(n ph_o) factor is the model's multipole structure made
53+% explicit: a mode with n nodal diameters is a 2^n-pole, so high-n modes
54+% radiate very little at low frequency, and none of them radiate on the axis.
55+%
56+% ldip is the one phenomenological length in the radiation model: the distance
57+% over which the air short-circuits from one face of the plate to the other.
58+% It is taken as half the shorter of the flexural and acoustic wavelengths,
59+% capped at the plate radius, and scaled by dipole_scale. Below the plate's
60+% coincidence frequency (about 18 kHz for 1 mm bronze) it is the flexural
61+% wavelength that sets it, which is why a cymbal is such a poor radiator in
62+% the low register and gets steadily better as the cascade fills the top.
63+
64+G = cymbal_profile(P);
65+E = P.youngs_modulus;
66+nu = P.poisson_ratio;
67+rho = P.density;
68+c = P.speed_of_sound;
69+rho0 = P.air_density;
70+nn = G.nn;
71+nel = nn - 1;
72+Lref = G.S / nel;
73+nmax = max(0, round(P.n_diameters));
74+
75+% ---- radiation sample points: two per element, Hermite-evaluated ----
76+nsub = 2;
77+nsurf = nel * nsub;
78+surfS = zeros(nsurf, 1);
79+surfLen = zeros(nsurf, 1);
80+surfEl = zeros(nsurf, 1);
81+surfXi = zeros(nsurf, 1);
82+for e = 1:nel
83+ L = G.snode(e + 1) - G.snode(e);
84+ for j = 1:nsub
85+ idx = (e - 1) * nsub + j;
86+ xi = (j - 0.5) / nsub;
87+ surfS(idx) = G.snode(e) + xi * L;
88+ surfLen(idx) = L / nsub;
89+ surfEl(idx) = e;
90+ surfXi(idx) = xi;
91+ end
92+end
93+surfR = reshape(interp1(G.s, G.r, surfS), nsurf, 1);
94+surfZ = reshape(interp1(G.s, G.z, surfS), nsurf, 1);
95+surfDr = reshape(interp1(G.s, G.drds, surfS), nsurf, 1);
96+surfDz = reshape(interp1(G.s, G.dzds, surfS), nsurf, 1);
97+
98+% ---- eigensolve per circumferential order ----
99+fs = [];
100+nd = [];
101+NUall = [];
102+NWall = [];
103+NTall = [];
104+SWall = [];
105+omLo = 2 * pi * P.min_mode_hz;
106+omHi = 2 * pi * P.max_mode_hz;
107+for n = 0:nmax
108+ [K, M] = cymbal_shell(G, n, E, nu, rho, Lref);
109+ Rc = chol(M);
110+ Am = (Rc') \ (K / Rc);
111+ Am = (Am + Am') / 2;
112+ [Uv, Sv, ~] = svd(Am);
113+ om2 = flipud(diag(Sv));
114+ Vv = fliplr(Uv);
115+ om = sqrt(max(0, om2)) * P.stiffness_trim;
116+ keep = find(om > omLo & om < omHi);
117+ if length(keep) == 0
118+ continue;
119+ end
120+ Phi = Rc \ Vv(:, keep);
121+ nk = length(keep);
122+ Uz = Phi(1:4:4 * nn, :); % meridional, nn x nk
123+ Wz = Phi(3:4:4 * nn, :); % normal
124+ Tz = Phi(4:4:4 * nn, :) * Lref; % dW/ds, unscaled back to physical
125+ % normal displacement at the radiation sample points (Hermite)
126+ SW = zeros(nsurf, nk);
127+ for idx = 1:nsurf
128+ e = surfEl(idx);
129+ xi = surfXi(idx);
130+ L = G.snode(e + 1) - G.snode(e);
131+ h1 = 1 - 3 * xi^2 + 2 * xi^3;
132+ h2 = L * (xi - 2 * xi^2 + xi^3);
133+ h3 = 3 * xi^2 - 2 * xi^3;
134+ h4 = L * (-xi^2 + xi^3);
135+ SW(idx, :) = h1 * Wz(e, :) + h2 * Tz(e, :) + h3 * Wz(e + 1, :) + h4 * Tz(e + 1, :);
136+ end
137+ fs = [fs; om(keep) / (2 * pi)];
138+ nd = [nd; n * ones(nk, 1)];
139+ NUall = [NUall, Uz];
140+ NWall = [NWall, Wz];
141+ NTall = [NTall, Tz];
142+ SWall = [SWall, SW];
143+end
144+N = length(fs);
145+if N == 0
146+ error('No cymbal modes between min_mode_hz and max_mode_hz -- check geometry/material.');
147+end
148+[fs, order] = sort(fs);
149+fs = reshape(fs, N, 1);
150+nd = reshape(nd(order), N, 1);
151+NUall = NUall(:, order);
152+NWall = NWall(:, order);
153+NTall = NTall(:, order);
154+SWall = SWall(:, order);
155+omega = 2 * pi * fs;
156+
157+% ---- the acoustic short-circuit length, per mode ----
158+tMean = 2 * pi * trapz(G.s, G.r .* G.thk) / G.area;
159+cL = sqrt(E / (rho * (1 - nu^2)));
160+kap = tMean / sqrt(12);
161+kFlex = sqrt(omega / (kap * cL));
162+kAc = omega / c;
163+ldip = P.dipole_scale * min(P.radius_m, pi ./ max(kFlex, kAc));
164+
165+% ---- radiated power, and the radiation damping that follows from it ----
166+nMu = 24;
167+mu = reshape(((1:nMu) - 0.5) / nMu * 2 - 1, 1, nMu); % cos(theta_o) midpoints
168+dmu = 2 / nMu;
169+sth = sqrt(max(0, 1 - mu.^2));
170+zetaRad = zeros(N, 1);
171+sigmaEff = zeros(N, 1);
172+for m = 1:N
173+ om = omega(m);
174+ k = om / c;
175+ nn_d = nd(m);
176+ alpha = (k * surfR) * sth; % nsurf x nMu
177+ Jn = besselj(nn_d, alpha);
178+ Jp1 = besselj(nn_d + 1, alpha);
179+ if nn_d == 0
180+ Jm1 = -Jp1;
181+ else
182+ Jm1 = besselj(nn_d - 1, alpha);
183+ end
184+ % i^(n+1) J_{n+1} + i^(n-1) J_{n-1} = i^n * i * (J_{n+1} - J_{n-1})
185+ ipow = 1i^nn_d;
186+ Tk = ipow * (2 * pi * Jn .* (surfDr * mu) ...
187+ - 1i * pi * (Jp1 - Jm1) .* (surfDz * sth));
188+ carrier = exp(1i * (k * surfZ) * mu);
189+ wgt = ldip(m) * surfR .* surfLen .* SWall(:, m);
190+ An = sum((wgt * ones(1, nMu)) .* Tk .* carrier, 1); % 1 x nMu
191+ Kn = sum(abs(An).^2) * dmu;
192+ if nn_d == 0
193+ piN = 2 * pi;
194+ else
195+ piN = pi;
196+ end
197+ zetaRad(m) = om^3 * rho0 * piN * Kn / (32 * pi^2 * c^3);
198+ % radiation efficiency, for the diagnostics only
199+ vv = sum(surfR .* surfLen .* SWall(:, m).^2) * piN;
200+ sigmaEff(m) = 2 * zetaRad(m) * om / max(1e-30, rho0 * c * vv);
201+end
202+
203+% ---- mount: a point dashpot on the axial motion at the inner edge ----
204+axHole = NUall(1, :).' * G.dzdsnode(1) + NWall(1, :).' * G.drdsnode(1);
205+zetaMount = P.mount_loss * axHole.^2 ./ (2 * omega);
206+
207+zeta = P.damp_rate ./ omega ...
208+ + 0.5 * (P.loss_factor + P.loss_slope * omega) ...
209+ + P.radiation_loss * zetaRad + zetaMount;
210+
211+basis.n = N;
212+basis.f = fs;
213+basis.omega = omega;
214+basis.zeta = zeta;
215+basis.zetaRad = zetaRad;
216+basis.zetaMount = zetaMount;
217+basis.sigma = sigmaEff;
218+basis.ndia = nd;
219+basis.ldip = ldip;
220+basis.NU = NUall;
221+basis.NW = NWall;
222+basis.NT = NTall;
223+basis.rnode = reshape(G.rnode, nn, 1);
224+basis.znode = reshape(G.znode, nn, 1);
225+basis.snode = reshape(G.snode, nn, 1);
226+basis.drdsnode = reshape(G.drdsnode, nn, 1);
227+basis.dzdsnode = reshape(G.dzdsnode, nn, 1);
228+basis.thknode = reshape(G.thknode, nn, 1);
229+basis.csnode = reshape(interp1(G.s, G.cs, G.snode), nn, 1);
230+basis.ctnode = reshape(interp1(G.s, G.ct, G.snode), nn, 1);
231+basis.surfR = surfR;
232+basis.surfZ = surfZ;
233+basis.surfDr = surfDr;
234+basis.surfDz = surfDz;
235+basis.surfLen = surfLen;
236+basis.SW = SWall;
237+basis.S = G.S;
238+basis.area = G.area;
239+basis.mass = G.mass;
240+basis.thkMean = tMean;
241+end
matlab/cymbal_cascade.madded+157−0View file
@@ -0,0 +1,157 @@
1+function [env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, taustr, sstr, duration)
2+% CYMBAL_CASCADE Modal excitation and the nonlinear flow of energy up the
3+% frequency ladder.
4+%
5+% This is the stage that makes the instrument a cymbal rather than a gong or a
6+% bell. A cymbal struck hard is not a linear resonator: the transverse
7+% displacement is comparable to the plate thickness, so the von Karman
8+% stretching of the mid-surface couples the modes, and energy injected in the
9+% low modes by the stick migrates upward over tens of milliseconds. That is
10+% what the ear hears as a crash swelling into shimmer after its attack, and it
11+% is why the same cymbal struck softly just pings: the coupling strength goes
12+% as the square of the amplitude, so a quiet stroke stays linear.
13+%
14+% Integrating the von Karman equations is out of reach here (the cubic
15+% coupling is a dense N^3 tensor, and the sibling repositories render in
16+% seconds). What is integrated instead is the ENERGY of each mode. Each mode
17+% is given an up-target, the lowest retained mode at least cascade_ratio times
18+% its frequency, and a down-target the same ratio below, and energy performs a
19+% biased random walk between them:
20+% dE_m/dt = -(1 + beta) Gam_m E_m + (arrivals from below and above)
21+% - 2 zeta_m om_m E_m
22+% Gam_m = cascade_gain om_m (w_m / t_mean)^2 = cascade_gain E_m / (om_m M t^2)
23+% w_m^2 = E_m / (om_m^2 M)
24+% with beta = cascade_back < 1 the downward-to-upward ratio, so the net flux is
25+% toward high frequency. Steps are taken in log frequency rather than in mode
26+% index on purpose: four-wave interactions move energy by a fraction of an
27+% octave at a time, not by one mode, and a ladder defined on the index would
28+% make the swell depend on how densely the band happens to be resolved. The
29+% rate is proportional to om_m because the higher a mode sits the faster its
30+% interactions run, and to the square of that mode's own
31+% amplitude-to-thickness ratio because that is the order of the geometric
32+% nonlinearity. Taking the amplitude locally rather than as a plate-wide rms
33+% is what keeps the cascade honest at the top of the band: the flux out of a
34+% mode goes as the SQUARE of its energy, so a nearly empty mode passes nothing
35+% on, the front of the cascade advances rather than the whole band draining at
36+% once, and the swell stops on its own as the energy thins out.
37+%
38+% The scheme conserves energy exactly when the damping is off and the top of
39+% the ladder reflects (cascade_leak = 0), which the validation checks; with
40+% cascade_leak = 1 the flux out of the top is lost, standing in for the
41+% truncated modes above max_mode_hz. cascade_gain is the one fitted number in
42+% the model: it is a dimensionless coupling constant, set so that a default
43+% crash takes about the 40 ms a real one takes to open up.
44+%
45+% What this buys is a cymbal that responds to dynamics the way the instrument
46+% does, at the cost of one control-rate loop rather than a full nonlinear
47+% integration. What it gives up is phase: energy moves between modes but the
48+% intermodulation products of a real cascade, and the split-pair beating of an
49+% out-of-round cymbal, are not reproduced. cascade_gain = 0 turns the whole
50+% mechanism off and leaves an exactly linear plate.
51+%
52+% Each contact adds ½ (dv_m)^2 to the energy of every mode, where
53+% dv_m = J S(om_m tau) (phi_m(x_c) . dhat)
54+% with S(x) = cos(x/2) / (1 - (x/pi)^2) the spectrum of the half-sine contact
55+% pulse and dhat the stroke direction (strike_angle_deg away from the surface
56+% normal, in the meridian plane). Energies add incoherently, which is right
57+% for a plate already ringing and exact for the first stroke into a silent one.
58+%
59+% Note that only the MAGNITUDE of the excitation survives into the energy
60+% state. The signed jumps DV and the frame each contact lands in are returned
61+% alongside, because cymbal_radiate.m needs them for the oscillator phases:
62+% amplitude comes from the cascade, phase from the linear impulse response.
63+%
64+% Returns env (N x nf+1, modal VELOCITY amplitude at each frame boundary),
65+% tframe (1 x nf+1), DV (N x number of contacts), kf (the frame index of each
66+% contact) and cinfo diagnostics.
67+
68+N = basis.n;
69+omega = basis.omega;
70+zeta = basis.zeta;
71+h = max(1e-4, P.cascade_step_ms * 1e-3);
72+nf = max(1, ceil(duration / h));
73+tframe = (0:nf) * h;
74+
75+% ---- modal velocity jump of each contact ----
76+I = length(tstr);
77+DV = zeros(N, I);
78+ang = P.strike_angle_deg * pi / 180;
79+for i = 1:I
80+ [uh, wh] = cymbal_at(basis, sstr(i));
81+ proj = -cos(ang) * wh + sin(ang) * uh; % 1 x N, along the stroke
82+ x = omega * taustr(i);
83+ den = 1 - (x / pi).^2;
84+ S = cos(x / 2) ./ den;
85+ sing = abs(den) < 1e-6;
86+ S(sing) = pi / 4; % limit at x -> pi
87+ DV(:, i) = jstr(i) * S .* reshape(proj, N, 1);
88+end
89+kf = min(nf, max(1, floor(tstr / h) + 1)); % frame each contact lands in
90+
91+decay = exp(-2 * zeta .* omega * h);
92+gmax = 0.5 / (h * (1 + P.cascade_back));
93+Mtot = basis.mass;
94+teff = basis.thkMean;
95+
96+% ---- the log-frequency ladder ----
97+% upIdx(m) is the lowest mode at least cascade_ratio times f(m); modes with no
98+% such partner point at N + 1, the sink standing in for the truncated band.
99+% dnIdx(m) is the highest mode at most f(m) / cascade_ratio, or m itself at the
100+% bottom, which reflects.
101+ratio = max(1.02, P.cascade_ratio);
102+f = basis.f;
103+upIdx = zeros(N, 1);
104+dnIdx = zeros(N, 1);
105+for m = 1:N
106+ ju = find(f > f(m) * ratio);
107+ if length(ju) == 0
108+ upIdx(m) = N + 1;
109+ else
110+ upIdx(m) = ju(1);
111+ end
112+ jd = find(f < f(m) / ratio);
113+ if length(jd) == 0
114+ dnIdx(m) = m;
115+ else
116+ dnIdx(m) = jd(length(jd));
117+ end
118+end
119+
120+env = zeros(N, nf + 1);
121+Ee = zeros(N, 1);
122+peakStrength = 0;
123+peakWrms = 0;
124+lost = 0;
125+for f = 1:nf
126+ hit = find(kf == f);
127+ for q = 1:length(hit)
128+ Ee = Ee + 0.5 * DV(:, hit(q)).^2;
129+ end
130+ env(:, f) = sqrt(2 * Ee);
131+
132+ wrms2 = sum(Ee ./ omega.^2) / Mtot;
133+ strength = wrms2 / teff^2;
134+ if strength > peakStrength
135+ peakStrength = strength;
136+ peakWrms = sqrt(wrms2);
137+ end
138+ gam = min(gmax, P.cascade_gain * Ee ./ (omega * Mtot * teff^2));
139+ up = gam .* Ee * h;
140+ dn = P.cascade_back * gam .* Ee * h;
141+ arrU = accumarray(upIdx, up, [N + 1, 1]);
142+ arrD = accumarray(dnIdx, dn, [N, 1]);
143+ leak = arrU(N + 1);
144+ En = Ee - up - dn + arrU(1:N) + arrD;
145+ En(N) = En(N) + (1 - P.cascade_leak) * leak;
146+ lost = lost + P.cascade_leak * leak;
147+ Ee = En .* decay;
148+end
149+env(:, nf + 1) = sqrt(2 * Ee);
150+
151+cinfo.frames = nf;
152+cinfo.stepS = h;
153+cinfo.peakStrength = peakStrength;
154+cinfo.peakWrmsUm = peakWrms * 1e6;
155+cinfo.cascadeLostJ = lost;
156+cinfo.injectedJ = sum(sum(0.5 * DV.^2));
157+end
matlab/cymbal_defaults.madded+76−0View file
@@ -0,0 +1,76 @@
1+function P = cymbal_defaults()
2+% CYMBAL_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+% Plate geometry (a 16-inch crash-ride)
7+P.radius_m = 0.2032;
8+P.hole_radius_m = 0.0064;
9+P.bell_radius_m = 0.0450;
10+P.bell_height_m = 0.0170;
11+P.bow_height_m = 0.0200;
12+P.thickness_bell_m = 0.00115;
13+P.thickness_rim_m = 0.00075;
14+
15+% Material (B20 bell bronze)
16+P.density = 8700.0;
17+P.youngs_modulus = 1.05e11;
18+P.poisson_ratio = 0.33;
19+P.stiffness_trim = 1.00;
20+
21+% Discretization / mode retention
22+P.n_elements = 28;
23+P.n_diameters = 48; % highest circumferential order solved for
24+P.min_mode_hz = 40.0;
25+P.max_mode_hz = 8000.0;
26+
27+% Damping
28+P.damp_rate = 0.35; % s^-1, frequency-independent decay rate
29+P.loss_factor = 1.5e-4;
30+P.loss_slope = 2.0e-9;
31+P.radiation_loss = 1.0;
32+P.mount_loss = 3.0; % N s/m dashpot at the inner edge (felt and stand)
33+
34+% Stick (nylon-tip drumstick)
35+P.tip_radius_m = 0.0050;
36+P.stick_mass_kg = 0.030;
37+P.stick_youngs = 5.0e9;
38+P.restitution = 0.55;
39+
40+% Strike gesture
41+P.strike_velocity = 3.0;
42+P.strike_r = 0.94; % 0 at the hole, 1 at the rim
43+P.strike_angle_deg = 40.0; % 0 = normal to the surface, 90 = tangential
44+P.strike_count = 1;
45+P.strike_rate_hz = 6.0;
46+P.roll_r2 = 0.80; % second contact radius when strike_count > 1
47+P.timing_jitter_ms = 8.0;
48+P.velocity_spread = 0.18;
49+P.seed = 1;
50+
51+% Nonlinear cascade
52+P.cascade_gain = 800.0;
53+P.cascade_back = 0.45; % downward / upward transfer ratio
54+P.cascade_ratio = 1.19; % frequency ratio of one transfer step (2^0.25)
55+P.cascade_leak = 1.0; % 1 = energy leaving the top mode is lost
56+P.cascade_step_ms = 0.5;
57+
58+% Radiation / listener
59+P.direct_gain = 1.0;
60+P.dipole_scale = 1.0;
61+P.n_azimuth = 96; % azimuthal samples in the listener transfer
62+P.listener_x_m = 0.30;
63+P.listener_y_m = 0.00;
64+P.listener_z_m = 0.50;
65+P.speed_of_sound = 343.0;
66+P.air_density = 1.204;
67+
68+% Rendering / output
69+P.sample_rate = 48000;
70+P.oversample = 2;
71+P.tail_s = 3.0;
72+P.pre_delay_ms = 3.0;
73+P.normalize = 1;
74+P.peak_dbfs = -3.0;
75+P.full_scale_pa = 5.0;
76+end
matlab/cymbal_demo.madded+26−0View file
@@ -0,0 +1,26 @@
1+% Render a cymbal crash with default parameters and print diagnostics.
2+% Set P.strike_count = 8 for a repeated ride figure, or P.cascade_gain = 0
3+% to hear the same plate as an exactly linear resonator.
4+% Run with: numbl run cymbal_demo.m
5+
6+P = cymbal_defaults();
7+tic;
8+[y, basis, dinfo] = cymbal_render(P);
9+elapsed = toc;
10+
11+fprintf('modes kept : %d up to n = %d nodal diameters (%.0f Hz to %.0f Hz)\n', ...
12+ dinfo.modesKept, dinfo.maxN, basis.f(1), basis.f(end));
13+fprintf('lowest : ');
14+fprintf('%.0f ', basis.f(1:min(10, basis.n)));
15+fprintf('Hz\n');
16+fprintf('strikes : %d at %.1f Hz (mean v %.2f m/s, mean contact %.0f us, peak force %.0f N)\n', ...
17+ dinfo.strikeCount, P.strike_rate_hz, dinfo.meanV, dinfo.meanTauUs, dinfo.peakForceN);
18+fprintf('cascade : peak w_rms %.1f um = %.2f thicknesses, %.1f%% of the injected energy left the top\n', ...
19+ dinfo.peakWrmsUm, sqrt(dinfo.peakStrength), 100 * dinfo.cascadeLostJ / max(1e-30, dinfo.injectedJ));
20+fprintf('peak press.: %.4f Pa (%.1f dB SPL) at %.0f cm\n', ...
21+ dinfo.peakPa, dinfo.peakSplDb, dinfo.listenerDistM * 100);
22+fprintf('length : %.2f s @ %d Hz\n', length(y) / P.sample_rate, P.sample_rate);
23+fprintf('render time: %.2f s\n', elapsed);
24+if dinfo.listenerInside
25+ fprintf('WARNING : listener is at or inside the plate\n');
26+end
matlab/cymbal_output.madded+65−0View file
@@ -0,0 +1,65 @@
1+function [y, peakPa] = cymbal_output(P, p, oversample)
2+% CYMBAL_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 triangle 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 .* (cy_i0(beta * sqrt(max(0, 1 - r.^2))) / cy_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 = cy_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/cymbal_point.madded+10−0View file
@@ -0,0 +1,10 @@
1+function [r, z] = cymbal_point(basis, sfrac)
2+% CYMBAL_POINT Meridian coordinates (r, z) at an arclength fraction.
3+nn = length(basis.snode);
4+s = min(0.999999, max(0, sfrac)) * basis.S;
5+idx = sum(basis.snode <= s);
6+idx = min(max(idx, 1), nn - 1);
7+xi = (s - basis.snode(idx)) / (basis.snode(idx + 1) - basis.snode(idx));
8+r = basis.rnode(idx) + xi * (basis.rnode(idx + 1) - basis.rnode(idx));
9+z = basis.znode(idx) + xi * (basis.znode(idx + 1) - basis.znode(idx));
10+end
matlab/cymbal_profile.madded+115−0View file
@@ -0,0 +1,115 @@
1+function G = cymbal_profile(P)
2+% CYMBAL_PROFILE Meridian of the cymbal as a shell of revolution.
3+%
4+% The cymbal is a surface of revolution about the z axis: a central mounting
5+% hole of radius hole_radius_m, a raised bell out to bell_radius_m, and a bow
6+% falling from the bell to the rim, which sits at z = 0 and r = radius_m.
7+% Arclength s runs from the hole (s = 0) to the rim (s = S), so r increases
8+% with s and the parameter strike_r addresses the plate as s / S.
9+%
10+% The bow is a parabola in the rim-referred coordinate u = (a - r) / (a - rb),
11+% z = bow_height_m u^2, which is flat at the rim and steepest where it meets
12+% the bell. The bell is the quartic in r that matches the bow's height, slope
13+% and curvature at r = rb and comes to a horizontal apex at r = 0. Matching
14+% the curvature matters: the finite element uses the meridian curvature and
15+% its arclength derivative, so a profile that is only C1 would put a delta
16+% function in the bending strain at the bell/bow join.
17+%
18+% Everything the element routine needs is tabulated on a fine grid and
19+% interpolated onto the Gauss points:
20+% r, z meridian coordinates (m)
21+% s arclength from the hole (m)
22+% thk shell thickness (m), tapering linearly from bell to rim
23+% drds dr/ds = cos(meridian angle)
24+% dzds dz/ds
25+% cs 1 / R_s, meridian curvature, positive where the surface curves
26+% away from the outward normal (a sphere of radius R gives 1/R)
27+% ct 1 / R_theta = -dz/ds / r, the circumferential curvature
28+% csp, ctp d/ds of those two
29+%
30+% Returns G with those fields plus S (total arclength), snode (node
31+% arclengths), and the node-sampled geometry.
32+
33+a = P.radius_m;
34+rh = max(1e-4, min(P.hole_radius_m, 0.5 * a));
35+rb = min(max(P.bell_radius_m, 2 * rh), 0.8 * a);
36+hb = P.bell_height_m;
37+hw = P.bow_height_m;
38+
39+M = 2001;
40+r = linspace(rh, a, M);
41+
42+% ---- bow: z = hw u^2 with u = (a - r) / (a - rb) ----
43+den = a - rb;
44+u = (a - r) / den;
45+zb = hw * u.^2;
46+zb_r = -2 * hw * u / den;
47+zb_rr = (2 * hw / den^2) * ones(1, M);
48+
49+% ---- bell: quartic matching z, z' and z'' at rb, horizontal apex at r = 0 ----
50+mb = -2 * hw / den; % dz/dr of the bow at r = rb
51+kb = 2 * hw / den^2; % d2z/dr2 of the bow at r = rb
52+A0 = hw + hb; % apex height (z at r = 0)
53+Amat = [rb^2, rb^3, rb^4; 2 * rb, 3 * rb^2, 4 * rb^3; 2, 6 * rb, 12 * rb^2];
54+coef = Amat \ [-hb; mb; kb];
55+B2 = coef(1); C3 = coef(2); E4 = coef(3);
56+zc = A0 + B2 * r.^2 + C3 * r.^3 + E4 * r.^4;
57+zc_r = 2 * B2 * r + 3 * C3 * r.^2 + 4 * E4 * r.^3;
58+zc_rr = 2 * B2 + 6 * C3 * r + 12 * E4 * r.^2;
59+
60+inBell = r < rb;
61+z = zb;
62+z_r = zb_r;
63+z_rr = zb_rr;
64+z(inBell) = zc(inBell);
65+z_r(inBell) = zc_r(inBell);
66+z_rr(inBell) = zc_rr(inBell);
67+
68+% ---- arclength and the differential geometry of the meridian ----
69+q = sqrt(1 + z_r.^2);
70+sarc = cumtrapz(r, q);
71+S = sarc(M);
72+drds = 1 ./ q;
73+dzds = z_r ./ q;
74+d2zds2 = z_rr ./ q.^4;
75+cs = -d2zds2 ./ drds;
76+ct = -dzds ./ r;
77+ctp = -d2zds2 ./ r + dzds .* drds ./ r.^2;
78+
79+% cs involves the third derivative of the profile; the fine grid is dense
80+% enough that a central difference is accurate to ~1e-8 relative
81+csp = zeros(1, M);
82+csp(2:M - 1) = (cs(3:M) - cs(1:M - 2)) ./ (sarc(3:M) - sarc(1:M - 2));
83+csp(1) = csp(2);
84+csp(M) = csp(M - 1);
85+
86+% ---- thickness: constant through the bell, linear taper across the bow ----
87+taper = min(1, max(0, (a - r) / den));
88+thk = P.thickness_rim_m + (P.thickness_bell_m - P.thickness_rim_m) * taper;
89+
90+% ---- nodes, equally spaced in arclength ----
91+nn = max(6, round(P.n_elements)) + 1;
92+snode = linspace(0, S, nn);
93+
94+G.r = r;
95+G.z = z;
96+G.s = sarc;
97+G.thk = thk;
98+G.drds = drds;
99+G.dzds = dzds;
100+G.cs = cs;
101+G.ct = ct;
102+G.csp = csp;
103+G.ctp = ctp;
104+G.S = S;
105+G.M = M;
106+G.nn = nn;
107+G.snode = snode;
108+G.rnode = interp1(sarc, r, snode);
109+G.znode = interp1(sarc, z, snode);
110+G.drdsnode = interp1(sarc, drds, snode);
111+G.dzdsnode = interp1(sarc, dzds, snode);
112+G.thknode = interp1(sarc, thk, snode);
113+G.area = 2 * pi * trapz(sarc, r);
114+G.mass = 2 * pi * P.density * trapz(sarc, r .* thk);
115+end
matlab/cymbal_radiate.madded+187−0View file
@@ -0,0 +1,187 @@
1+function [p, radinfo] = cymbal_radiate(P, basis, env, tframe, DV, kf, tstr, taustr, fstr, sstr, fsHi)
2+% CYMBAL_RADIATE Free-field pressure at the listener from the modal envelopes.
3+%
4+% Two contributions:
5+%
6+% 1. The mode bank. Each retained mode is a normal-dipole sheet on the plate:
7+% a surface element of area dS carries a dipole of volume dS * ldip along
8+% the local normal (ldip from cymbal_basis.m), and the element dipoles are
9+% summed coherently with their relative retardations at the mode's centre
10+% frequency,
11+% G = sum_e (i om rho0 / 4 pi) V_e (u_e . rhat_e) (i k + 1/r_e)
12+% exp(-i k (r_e - rcen)) / r_e,
13+% over n_azimuth azimuths and two sample points per meridian element. The
14+% common carrier exp(-i k rcen) is stripped, so relative phase across the
15+% plate is kept while the overall propagation delay is not double counted.
16+% Unlike the far-field form used for the radiated power, this keeps the
17+% 1/r^2 near-field term and the true distances, which matters: a listener
18+% half a metre from a 40 cm plate is nowhere near the far field. Only the
19+% normal component of the surface motion radiates; the in-plane component
20+% of a thin plate displaces almost no air.
21+%
22+% Because the azimuthal shape is exactly cos(n theta), the sum over
23+% azimuths carries the multipole cancellation of the mode: on the axis of
24+% symmetry every mode with n >= 1 cancels exactly, and off axis the high-n
25+% modes are weak until their wavelength is short. Moving the listener
26+% re-runs only this stage.
27+%
28+% The synthesis itself is an oscillator bank rather than the closed-form
29+% damped exponentials of the sibling repositories: the cascade makes every
30+% amplitude time-varying, so each mode is a fixed-frequency oscillator
31+% whose envelope is the cascade solution, interpolated linearly across each
32+% control frame. Within a frame the two matrix products below are exact,
33+% and the per-frame phase advance is a unit-modulus multiply, so nothing
34+% drifts.
35+%
36+% Amplitude comes from the cascade; PHASE comes from a linear impulse
37+% response carried alongside it. Each oscillator keeps a unit phasor and the
38+% magnitude the mode would have had with no cascade at all; a contact adds
39+% its signed modal velocity jump to that complex state and the phasor is
40+% renormalized. A mode that was quiet is therefore re-phased by the stroke,
41+% which is what makes the attack a coherent click, while a mode already
42+% ringing loudly keeps the phase it had, which is also right. What is lost
43+% is the intermodulation the real cascade would produce, since the energy
44+% it moves into a mode arrives with whatever phase that mode already had.
45+%
46+% 2. Direct strike clicks. Modes above max_mode_hz are truncated; their
47+% collective off-resonance response is the immediate local motion of the
48+% plate under the tip. It is reinstated from the driving-point impedance of
49+% an infinite plate, which for a plate is a pure resistance
50+% Z = 8 sqrt(D rho t): the contact force divided by Z is the local
51+% velocity, which radiates as a dipole over about half a flexural
52+% wavelength, evaluated at the pulse's own centre frequency om* = pi / tau.
53+% Scaled by direct_gain. This is the least principled term of the model,
54+% and direct_gain = 0 turns it off.
55+%
56+% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.
57+
58+c = P.speed_of_sound;
59+rho0 = P.air_density;
60+N = basis.n;
61+nsurf = length(basis.surfR);
62+nth = max(8, round(P.n_azimuth));
63+lx = P.listener_x_m;
64+ly = P.listener_y_m;
65+lz = P.listener_z_m;
66+
67+% ---- surface element geometry, meridian samples x azimuths ----
68+th = reshape(((0:nth - 1) + 0.5) / nth * 2 * pi, 1, nth);
69+dth = 2 * pi / nth;
70+ct = cos(th);
71+st = sin(th);
72+ex = basis.surfR * ct; % nsurf x nth
73+ey = basis.surfR * st;
74+ez = basis.surfZ * ones(1, nth);
75+nx = (-basis.surfDz) * ct;
76+ny = (-basis.surfDz) * st;
77+nz = basis.surfDr * ones(1, nth);
78+dSa = (basis.surfR .* basis.surfLen) * (dth * ones(1, nth));
79+
80+dx = lx - ex;
81+dy = ly - ey;
82+dz = lz - ez;
83+re = sqrt(dx.^2 + dy.^2 + dz.^2);
84+listenerInside = min(min(re)) < 0.02;
85+re = max(re, 0.01);
86+cosg = (nx .* dx + ny .* dy + nz .* dz) ./ re; % nsurf x nth
87+
88+cx = 0;
89+cy = 0;
90+cz = sum(basis.surfR .* basis.surfLen .* basis.surfZ) / sum(basis.surfR .* basis.surfLen);
91+rcen = sqrt((lx - cx)^2 + (ly - cy)^2 + (lz - cz)^2);
92+
93+% ---- per-mode listener transfer ----
94+G = zeros(N, 1) + 0i;
95+for m = 1:N
96+ om = basis.omega(m);
97+ k = om / c;
98+ azi = cos(basis.ndia(m) * th); % 1 x nth
99+ amp = (basis.ldip(m) * dSa .* (basis.SW(:, m) * ones(1, nth))) .* (ones(nsurf, 1) * azi);
100+ ker = amp .* cosg .* (1i * k + 1 ./ re) .* exp(-1i * k * (re - rcen)) ./ re;
101+ G(m) = (1i * om * rho0 / (4 * pi)) * sum(sum(ker));
102+end
103+
104+% ---- oscillator bank driven by the cascade envelopes ----
105+T = floor(tframe(end) * fsHi);
106+p = zeros(1, T);
107+nf = length(tframe) - 1;
108+L = max(1, round((tframe(2) - tframe(1)) * fsHi));
109+omd = basis.omega .* sqrt(max(1e-12, 1 - basis.zeta.^2));
110+Eb = exp(1i * omd * ((0:L - 1) / fsHi)); % N x L, hoisted
111+stepPh = exp(1i * omd * (L / fsHi)); % per-frame advance
112+decLin = exp(-basis.zeta .* basis.omega * (L / fsHi));
113+ramp = (0:L - 1) / L;
114+sPh = ones(N, 1) + 0i;
115+aLin = zeros(N, 1);
116+pos = 1;
117+for f = 1:nf
118+ if pos > T
119+ break;
120+ end
121+ hit = find(kf == f);
122+ for q = 1:length(hit)
123+ cc = aLin .* sPh + DV(:, hit(q));
124+ aLin = abs(cc);
125+ sPh = cc ./ max(1e-300, aLin);
126+ end
127+ len = min(L, T - pos + 1);
128+ A0 = env(:, f);
129+ A1 = env(:, f + 1);
130+ c1 = (G .* A0 .* sPh).';
131+ c2 = (G .* (A1 - A0) .* sPh).';
132+ seg = real((c1 * Eb(:, 1:len)) + (c2 * Eb(:, 1:len)) .* ramp(1:len));
133+ p(pos:pos + len - 1) = seg;
134+ sPh = sPh .* stepPh;
135+ aLin = aLin .* decLin;
136+ pos = pos + L;
137+end
138+
139+% ---- direct strike clicks ----
140+E = P.youngs_modulus;
141+nu = P.poisson_ratio;
142+rho = P.density;
143+cL = sqrt(E / (rho * (1 - nu^2)));
144+I = length(tstr);
145+ks = min(T, max(1, round(tstr * fsHi) + 1));
146+ang = P.strike_angle_deg * pi / 180;
147+directPeak = 0;
148+for i = 1:I
149+ [~, ~, thk, drds, dzds] = cymbal_at(basis, sstr(i));
150+ [rs, zs] = cymbal_point(basis, sstr(i));
151+ dvx = cos(ang) * dzds + sin(ang) * drds; % stroke direction at theta = 0
152+ dvz = -cos(ang) * drds + sin(ang) * dzds;
153+ rsx = lx - rs;
154+ rsy = ly;
155+ rsz = lz - zs;
156+ rr = max(0.01, sqrt(rsx^2 + rsy^2 + rsz^2));
157+ cg = (rsx * dvx + rsz * dvz) / rr;
158+ nP = max(2, round(taustr(i) * fsHi));
159+ % the click carries the true path delay relative to the carrier that was
160+ % stripped from G, which is negative when the contact is nearer the
161+ % listener than the plate centroid
162+ k0 = max(1, ks(i) + round((rr - rcen) / c * fsHi));
163+ if k0 + nP > T
164+ continue;
165+ end
166+ D = E * thk^3 / (12 * (1 - nu^2));
167+ Zp = 8 * sqrt(D * rho * thk);
168+ omStar = pi / taustr(i);
169+ kap = thk / sqrt(12);
170+ kB = sqrt(omStar / (kap * cL));
171+ lamB = 2 * pi / kB;
172+ Aeff = pi * (lamB / 4)^2; % patch driven by the contact
173+ ldir = pi / max(kB, omStar / c);
174+ mu = Aeff * ldir / Zp; % dipole volume per unit force
175+ amp = P.direct_gain * fstr(i) * mu * rho0 * omStar^2 * cg / (4 * pi * rr * c);
176+ pulse = amp * sin(pi * (0:nP) / nP);
177+ p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;
178+ if abs(amp) > directPeak
179+ directPeak = abs(amp);
180+ end
181+end
182+
183+radinfo.gainAbs = abs(G);
184+radinfo.directPeakPa = directPeak;
185+radinfo.listenerInside = listenerInside;
186+radinfo.listenerDistM = rcen;
187+end
matlab/cymbal_render.madded+33−0View file
@@ -0,0 +1,33 @@
1+function [y, basis, dinfo] = cymbal_render(P)
2+% CYMBAL_RENDER Full render: basis + strikes -> cascade -> radiate -> output.
3+% Convenience composition of the stage functions; the web app calls the stages
4+% separately so it can cache intermediates (the shell eigensolves above all).
5+
6+basis = cymbal_basis(P);
7+fsHi = P.sample_rate * P.oversample;
8+[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = cymbal_strikes(P, basis);
9+[env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, taustr, sstr, sinfo.durationS);
10+[p, radinfo] = cymbal_radiate(P, basis, env, tframe, DV, kf, tstr, taustr, fstr, sstr, fsHi);
11+[y, peakPa] = cymbal_output(P, p, P.oversample);
12+
13+dinfo.strikeCount = sinfo.count;
14+dinfo.meanV = sinfo.meanV;
15+dinfo.maxV = sinfo.maxV;
16+dinfo.meanTauUs = sinfo.meanTauUs;
17+dinfo.peakForceN = sinfo.peakForceN;
18+dinfo.durationS = sinfo.durationS;
19+dinfo.modesKept = basis.n;
20+dinfo.maxN = max(basis.ndia);
21+dinfo.f1 = basis.f(1);
22+dinfo.peakPa = peakPa;
23+dinfo.peakSplDb = 20 * log10(peakPa / 20e-6 + 1e-30);
24+dinfo.directPeakPa = radinfo.directPeakPa;
25+dinfo.listenerInside = radinfo.listenerInside;
26+dinfo.listenerDistM = radinfo.listenerDistM;
27+dinfo.peakStrength = cinfo.peakStrength;
28+dinfo.peakWrmsUm = cinfo.peakWrmsUm;
29+dinfo.injectedJ = cinfo.injectedJ;
30+dinfo.cascadeLostJ = cinfo.cascadeLostJ;
31+dinfo.cascadeFrames = cinfo.frames;
32+dinfo.sinfo = sinfo;
33+end
matlab/cymbal_rng.madded+12−0View file
@@ -0,0 +1,12 @@
1+function [r, state] = cymbal_rng(state, n)
2+% CYMBAL_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+% strike 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/cymbal_shell.madded+133−0View file
@@ -0,0 +1,133 @@
1+function [K, M] = cymbal_shell(G, n, E, nu, rho, Lref)
2+% CYMBAL_SHELL Free-free thin-shell system for one circumferential order n.
3+%
4+% The cymbal is a shell of revolution, so its modes separate exactly in the
5+% azimuth: every mode has n nodal diameters and is
6+% u(s, th) = U(s) cos(n th) (meridional, along the surface toward the rim)
7+% v(s, th) = V(s) sin(n th) (circumferential)
8+% w(s, th) = W(s) cos(n th) (along the outward normal)
9+% with a partner rotated by pi / 2n. A strike on the meridian th = 0 excites
10+% only the cos member of each pair, which is why the model solves one member
11+% per (n, k) and says so in the README. The 1-D energy functional below is the
12+% th integral of the Love-Kirchhoff shell energy; every quadratic term carries
13+% cos^2 or sin^2, so the integral contributes the same constant factor (pi, or
14+% 2 pi at n = 0) to K and M and drops out of the eigenvalues. At n = 0 the V
15+% field decouples and gives the axisymmetric torsional modes, which have no
16+% normal motion and therefore radiate nothing.
17+%
18+% Strains (Kraus, with the curvatures signed so that a sphere of radius R has
19+% 1 / R_s = 1 / R_theta = +1 / R for the outward normal):
20+% eps_s = U' + W cs
21+% eps_th = (n V + U r') / r + W ct
22+% gamma = -n U / r + V' - V r' / r (engineering shear)
23+% beta_s = U cs - W' (rotation of the normal)
24+% beta_th = n W / r + V ct
25+% kap_s = beta_s' = U' cs + U cs' - W''
26+% kap_th = (n beta_th + beta_s r') / r
27+% kap_sth = -n beta_s / r + beta_th' - beta_th r' / r (the full twist,
28+% factor 2 included)
29+% and the energy densities are the standard C = E t / (1 - nu^2) membrane and
30+% D = E t^3 / (12 (1 - nu^2)) bending forms, integrated with weight r ds.
31+%
32+% Elements have two nodes and four degrees of freedom each: U and V are linear
33+% (C0), W is cubic Hermite (C1) with dW/ds as its second freedom, so beta_s and
34+% its derivative are both well defined. Four-point Gauss quadrature per
35+% element evaluates the tabulated geometry at the quadrature points rather
36+% than freezing it per element, so a curved meridian is integrated rather than
37+% approximated by conical frusta.
38+%
39+% Rotary inertia rho t^3 / 12 (beta_s^2 + beta_th^2) is kept. It shifts the
40+% retained frequencies by well under a percent, but it is what gives the
41+% dW/ds freedoms mass and so makes M positive definite, which is what lets the
42+% eigensolve reduce through chol(M). The dW/ds freedoms are then scaled by
43+% Lref (a congruence, so the eigenvalues are untouched) to bring them onto the
44+% same numerical scale as the displacements.
45+%
46+% Returns K and M, each 4 nn x 4 nn, symmetric, M positive definite. The
47+% returned system is in the SCALED freedoms; multiply an eigenvector by
48+% diag(1, 1, 1, Lref) per node to recover physical dW/ds.
49+
50+nn = G.nn;
51+ndof = 4 * nn;
52+K = zeros(ndof);
53+M = zeros(ndof);
54+
55+% four-point Gauss on [0, 1]
56+gx = [0.5 - 0.5 * 0.8611363115940526, 0.5 - 0.5 * 0.3399810435848563, ...
57+ 0.5 + 0.5 * 0.3399810435848563, 0.5 + 0.5 * 0.8611363115940526];
58+gw = [0.3478548451374538, 0.6521451548625461, ...
59+ 0.6521451548625461, 0.3478548451374538] / 2;
60+
61+for e = 1:nn - 1
62+ s0 = G.snode(e);
63+ s1 = G.snode(e + 1);
64+ L = s1 - s0;
65+ Ke = zeros(8);
66+ Me = zeros(8);
67+ for g = 1:4
68+ xi = gx(g);
69+ sg = s0 + xi * L;
70+ rr = interp1(G.s, G.r, sg);
71+ tt = interp1(G.s, G.thk, sg);
72+ dr = interp1(G.s, G.drds, sg);
73+ cs = interp1(G.s, G.cs, sg);
74+ ct = interp1(G.s, G.ct, sg);
75+ csp = interp1(G.s, G.csp, sg);
76+ ctp = interp1(G.s, G.ctp, sg);
77+
78+ % shape functions, DOF order [U1 V1 W1 W1' U2 V2 W2 W2']
79+ NU = [1 - xi, 0, 0, 0, xi, 0, 0, 0];
80+ NUp = [-1 / L, 0, 0, 0, 1 / L, 0, 0, 0];
81+ NV = [0, 1 - xi, 0, 0, 0, xi, 0, 0];
82+ NVp = [0, -1 / L, 0, 0, 0, 1 / L, 0, 0];
83+ h1 = 1 - 3 * xi^2 + 2 * xi^3;
84+ h2 = xi - 2 * xi^2 + xi^3;
85+ h3 = 3 * xi^2 - 2 * xi^3;
86+ h4 = -xi^2 + xi^3;
87+ d1 = -6 * xi + 6 * xi^2;
88+ d2 = 1 - 4 * xi + 3 * xi^2;
89+ d3 = 6 * xi - 6 * xi^2;
90+ d4 = -2 * xi + 3 * xi^2;
91+ e1 = -6 + 12 * xi;
92+ e2 = -4 + 6 * xi;
93+ e3 = 6 - 12 * xi;
94+ e4 = -2 + 6 * xi;
95+ NW = [0, 0, h1, L * h2, 0, 0, h3, L * h4];
96+ NWp = [0, 0, d1 / L, d2, 0, 0, d3 / L, d4];
97+ NWpp = [0, 0, e1 / L^2, e2 / L, 0, 0, e3 / L^2, e4 / L];
98+
99+ Bes = NUp + cs * NW;
100+ Bet = (n * NV + dr * NU) / rr + ct * NW;
101+ Bg = -n * NU / rr + NVp - dr * NV / rr;
102+ Bbs = cs * NU - NWp;
103+ Bbt = n * NW / rr + ct * NV;
104+ Bks = csp * NU + cs * NUp - NWpp;
105+ Bkt = (n * Bbt + dr * Bbs) / rr;
106+ Bbtp = n * (NWp / rr - dr * NW / rr^2) + ctp * NV + ct * NVp;
107+ Bkst = -n * Bbs / rr + Bbtp - dr * Bbt / rr;
108+
109+ Cm = E * tt / (1 - nu^2);
110+ Db = E * tt^3 / (12 * (1 - nu^2));
111+ wq = gw(g) * L * rr;
112+
113+ Ke = Ke + wq * (Cm * (Bes' * Bes + Bet' * Bet ...
114+ + nu * (Bes' * Bet + Bet' * Bes) ...
115+ + (1 - nu) / 2 * (Bg' * Bg)) ...
116+ + Db * (Bks' * Bks + Bkt' * Bkt ...
117+ + nu * (Bks' * Bkt + Bkt' * Bks) ...
118+ + (1 - nu) / 2 * (Bkst' * Bkst)));
119+ Me = Me + wq * (rho * tt * (NU' * NU + NV' * NV + NW' * NW) ...
120+ + rho * tt^3 / 12 * (Bbs' * Bbs + Bbt' * Bbt));
121+ end
122+ d = [4 * e - 3, 4 * e - 2, 4 * e - 1, 4 * e, 4 * e + 1, 4 * e + 2, 4 * e + 3, 4 * e + 4];
123+ K(d, d) = K(d, d) + Ke;
124+ M(d, d) = M(d, d) + Me;
125+end
126+
127+% scale the dW/ds freedoms so every freedom is numerically comparable
128+sc = ones(ndof, 1);
129+sc(4:4:ndof) = Lref;
130+Sd = sc * sc';
131+K = (K .* Sd + (K .* Sd)') / 2;
132+M = (M .* Sd + (M .* Sd)') / 2;
133+end
matlab/cymbal_strikes.madded+119−0View file
@@ -0,0 +1,119 @@
1+function [tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = cymbal_strikes(P, basis)
2+% CYMBAL_STRIKES Stick contacts: a single strike, or a seeded repeated figure.
3+%
4+% strike_count = 1 gives one deterministic strike of strike_velocity at
5+% strike_r at t = 0; no random numbers are drawn. strike_count > 1 repeats it
6+% at strike_rate_hz, alternating between strike_r and roll_r2 (a ride pattern
7+% moving between the bow and the edge), with Gaussian timing jitter
8+% (timing_jitter_ms) and lognormal velocity spread (velocity_spread). The
9+% first hit stays at t = 0 so the onset is fixed. All randomness comes from
10+% cymbal_rng seeded with P.seed, so a given seed reproduces the identical
11+% pattern in the CLI, in Node, and in the browser. Every contact sits on the
12+% meridian theta = 0; see cymbal_basis.m for why.
13+%
14+% Each contact is a Hertz impact of the stick tip on the plate:
15+% E* = 1 / ((1 - nu_t^2) / E_t + (1 - nu_p^2) / E_p)
16+% R_eff= 1 / sqrt((1/r_tip + 1/R_s)(1/r_tip + 1/R_theta)) (sphere on a
17+% doubly curved surface, at the local principal curvatures)
18+% k_h = (4/3) E* sqrt(R_eff)
19+% tau = 3.218 (m*^2 / (k_h^2 v))^{1/5} (Hertz contact time)
20+% J = (1 + restitution) m* v
21+% Fpk = pi J / (2 tau) (half-sine peak)
22+%
23+% The reduced mass m* pairs the stick with the plate mass engaged during the
24+% contact. Here the sibling triangle's flexural-group-velocity estimate is
25+% replaced by something exact: the driving-point impedance of an infinite
26+% plate is a pure resistance, Z = 8 sqrt(D rho t), independent of frequency,
27+% so the mass the stick sees during a pulse of duration tau is simply Z tau. A
28+% short fixed-point iteration closes the loop, since tau depends on m* which
29+% depends on tau. What this captures is that a heavy stick, a soft tip or a
30+% slow stroke all lengthen the contact and darken the strike; what it leaves
31+% out is the millisecond detail of the contact itself, which is prescribed
32+% rather than simulated.
33+%
34+% Returns row vectors sorted by time: tstr (s), vstr (m/s), jstr (N s),
35+% taustr (s), fstr (peak N), sstr (arclength fraction of the contact), and
36+% sinfo with summary statistics.
37+
38+K = max(1, round(P.strike_count));
39+rate = max(1e-3, P.strike_rate_hz);
40+duration = (K - 1) / rate + P.tail_s;
41+
42+E = P.youngs_modulus;
43+nu = P.poisson_ratio;
44+rho = P.density;
45+nut = 0.35; % stick tip Poisson ratio, fixed
46+Estar = 1 / ((1 - nut^2) / P.stick_youngs + (1 - nu^2) / E);
47+rt = max(1e-4, P.tip_radius_m);
48+
49+state = mod(floor(abs(P.seed)), 4294967296);
50+z1 = zeros(K, 1);
51+z2 = zeros(K, 1);
52+if K > 1
53+ [u, state] = cymbal_rng(state, 4 * K);
54+ U = reshape(u, K, 4);
55+ z1 = sqrt(-2 * log(U(:, 1))) .* cos(2 * pi * U(:, 2));
56+ z2 = sqrt(-2 * log(U(:, 3))) .* cos(2 * pi * U(:, 4));
57+end
58+
59+tstr = zeros(1, K);
60+vstr = zeros(1, K);
61+jstr = zeros(1, K);
62+taustr = zeros(1, K);
63+fstr = zeros(1, K);
64+sstr = zeros(1, K);
65+for i = 1:K
66+ t = (i - 1) / rate;
67+ if i > 1
68+ t = max(0, t + P.timing_jitter_ms * 1e-3 * z1(i));
69+ end
70+ v = P.strike_velocity;
71+ if K > 1
72+ v = v * exp(P.velocity_spread * z2(i));
73+ end
74+ v = max(0.01, v);
75+ if mod(i, 2) == 1
76+ sfrac = P.strike_r;
77+ else
78+ sfrac = P.roll_r2;
79+ end
80+ sfrac = min(0.999, max(0.001, sfrac));
81+
82+ % local plate properties at the contact
83+ [~, ~, thk, ~, ~, cs, ct] = cymbal_at(basis, sfrac);
84+ Rrel = 1 / sqrt(max(1e-9, (1 / rt + cs) * (1 / rt + ct)));
85+ kh = (4 / 3) * Estar * sqrt(Rrel);
86+ D = E * thk^3 / (12 * (1 - nu^2));
87+ Zp = 8 * sqrt(D * rho * thk); % driving-point resistance of the plate
88+
89+ tau = 1e-4;
90+ mstar = P.stick_mass_kg;
91+ for it = 1:6
92+ mloc = Zp * tau;
93+ mstar = P.stick_mass_kg * mloc / (P.stick_mass_kg + mloc);
94+ tau = 3.218 * (mstar^2 / (kh^2 * v))^(1 / 5);
95+ end
96+
97+ tstr(i) = t;
98+ vstr(i) = v;
99+ jstr(i) = (1 + P.restitution) * mstar * v;
100+ taustr(i) = tau;
101+ fstr(i) = pi * jstr(i) / (2 * tau);
102+ sstr(i) = sfrac;
103+end
104+
105+[tstr, order] = sort(tstr);
106+tstr = reshape(tstr, 1, K);
107+vstr = reshape(vstr(order), 1, K);
108+jstr = reshape(jstr(order), 1, K);
109+taustr = reshape(taustr(order), 1, K);
110+fstr = reshape(fstr(order), 1, K);
111+sstr = reshape(sstr(order), 1, K);
112+
113+sinfo.count = K;
114+sinfo.meanV = sum(vstr) / K;
115+sinfo.maxV = max(vstr);
116+sinfo.meanTauUs = sum(taustr) / K * 1e6;
117+sinfo.peakForceN = max(fstr);
118+sinfo.durationS = duration;
119+end
matlab/cymbal_validate.madded+288−0View file
@@ -0,0 +1,288 @@
1+% Validate the cymbal model. There is no external reference implementation, so
2+% the checks are (a) closed-form identities the physics must satisfy, (b)
3+% structural identities of the shell eigensolve and of the radiation integral,
4+% (c) exact reproducibility of the seeded strike process, and (d) pinned
5+% end-to-end diagnostics from the default configuration, recorded when the
6+% model was built and treated as a regression reference since.
7+% Run with: numbl run cymbal_validate.m
8+
9+fails = 0;
10+
11+% ---- 1. flat free circular plate against the closed form ----
12+% The same shell routine that builds the cymbal, run on a flat meridian of
13+% constant thickness with a pinhole at the centre, must reproduce the
14+% completely free circular plate: the roots of the 2 x 2 determinant built
15+% from J_n and I_n under vanishing radial moment and Kelvin-Kirchhoff shear,
16+% and exactly two rigid-body modes at n = 0 (axial translation and rotation
17+% about the axis), two at n = 1 (translation and tilt) and none above.
18+P = cymbal_defaults();
19+Q = P;
20+Q.bow_height_m = 0; Q.bell_height_m = 0;
21+Q.thickness_bell_m = 0.001; Q.thickness_rim_m = 0.001;
22+Q.hole_radius_m = 0.002; Q.radius_m = 0.2; Q.n_elements = 40;
23+Gf = cymbal_profile(Q);
24+Dp = Q.youngs_modulus * 0.001^3 / (12 * (1 - Q.poisson_ratio^2));
25+scal = sqrt(Dp / (Q.density * 0.001)) / Q.radius_m^2 / (2 * pi);
26+Lr = Gf.S / (Gf.nn - 1);
27+errP = 0;
28+rigid = zeros(1, 5);
29+for n = 0:4
30+ f = cv_freqs(Gf, n, Q, Lr);
31+ rigid(n + 1) = sum(f < 3);
32+ fe = f(f > 3);
33+ lam = cv_plate_roots(n, Q.poisson_ratio);
34+ fc = lam.^2 * scal;
35+ nc = min(4, min(length(fe), length(fc)));
36+ errP = max(errP, max(abs(fe(1:nc).' - fc(1:nc)) ./ fc(1:nc)));
37+end
38+ok = errP < 3e-3 && rigid(1) == 2 && rigid(2) == 2 && sum(rigid(3:5)) == 0;
39+fails = fails + (~ok);
40+fprintf('flat plate : max freq err %.1e, rigid modes %d %d %d %d %d %s\n', ...
41+ errP, rigid(1), rigid(2), rigid(3), rigid(4), rigid(5), cv_verdict(ok));
42+
43+% ---- 2. eigensolve identities on the cymbal itself ----
44+% The Cholesky reduction must hand back mass-orthonormal shapes that
45+% diagonalize the stiffness: phi' M phi = I and phi' K phi = diag(om^2).
46+G = cymbal_profile(P);
47+Lref = G.S / (G.nn - 1);
48+[K, M] = cymbal_shell(G, 5, P.youngs_modulus, P.poisson_ratio, P.density, Lref);
49+Rc = chol(M);
50+Am = (Rc') \ (K / Rc);
51+Am = (Am + Am') / 2;
52+[Uv, Sv, ~] = svd(Am);
53+om2 = flipud(diag(Sv));
54+Phi = Rc \ fliplr(Uv);
55+eOrth = max(max(abs(Phi' * M * Phi - eye(length(om2)))));
56+eDiag = max(max(abs(Phi' * K * Phi - diag(om2)))) / max(om2);
57+ok = eOrth < 1e-8 && eDiag < 1e-10;
58+fails = fails + (~ok);
59+fprintf('eigensolve : orthonormality %.1e, diagonalization %.1e %s\n', eOrth, eDiag, cv_verdict(ok));
60+
61+% ---- 3. the Bessel far field against a direct sum over the dipole sheet ----
62+% cymbal_basis reduces the azimuthal integral of the far-field kernel to
63+% Bessel functions in order to get the radiated power; that reduction must
64+% agree with brute force summation over a finely discretized dipole sheet.
65+Q = P;
66+Q.n_diameters = 4; Q.max_mode_hz = 1200;
67+b4 = cymbal_basis(Q);
68+errF = 0;
69+for m = [1 3 5 9]
70+ [pd, pb] = cv_farfield(Q, b4, m, 200, 720, 0.6, 0.4);
71+ errF = max(errF, abs(pb / pd - 1));
72+end
73+ok = errF < 0.01;
74+fails = fails + (~ok);
75+fprintf('far field : Bessel vs direct dipole sum, max err %.2e %s\n', errF, cv_verdict(ok));
76+
77+% ---- 4. on the axis, every mode with a nodal diameter cancels exactly ----
78+% A cos(n theta) source on a surface of revolution has no monopole, dipole or
79+% any lower multipole than 2^n, and on the symmetry axis all of them vanish.
80+% This is the sharpest available check that the azimuthal sum in the listener
81+% transfer is assembled correctly.
82+Q = P;
83+Q.n_diameters = 6; Q.max_mode_hz = 1500;
84+Q.listener_x_m = 0; Q.listener_y_m = 0; Q.listener_z_m = 0.6;
85+b6 = cymbal_basis(Q);
86+[tstr, ~, jstr, taustr, fstr, sstr, sinf] = cymbal_strikes(Q, b6);
87+[env, tfr, DV, kf] = cymbal_cascade(Q, b6, tstr, jstr, taustr, sstr, 0.05);
88+[~, ri] = cymbal_radiate(Q, b6, env, tfr, DV, kf, tstr, taustr, fstr, sstr, Q.sample_rate);
89+g0 = max(ri.gainAbs(b6.ndia == 0));
90+g1 = max(ri.gainAbs(b6.ndia >= 1));
91+ok = g1 / g0 < 1e-10;
92+fails = fails + (~ok);
93+fprintf('on axis : |G| of n>=1 modes is %.1e of the n=0 modes %s\n', g1 / g0, cv_verdict(ok));
94+
95+% ---- 5. Hertz contact-time scaling tau ~ v^(-1/5) ----
96+% With a stick much lighter than the plate mass engaged during the contact the
97+% reduced mass is constant and the free-sphere scaling must emerge.
98+Q = P;
99+Q.strike_count = 1;
100+Q.stick_mass_kg = 2e-5;
101+[~, ~, ~, tau1] = cymbal_strikes(Q, b6);
102+Q.strike_velocity = 2 * P.strike_velocity;
103+[~, ~, ~, tau2] = cymbal_strikes(Q, b6);
104+ratio = tau2(1) / tau1(1);
105+ok = abs(ratio - 2^(-1 / 5)) < 0.03;
106+fails = fails + (~ok);
107+fprintf('hertz : tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), cv_verdict(ok));
108+
109+% ---- 6. the cascade conserves energy, and switches off cleanly ----
110+% With the damping removed and the top of the ladder reflecting, the transfer
111+% scheme is a Markov generator on the modal energies and must conserve their
112+% sum to round-off. With cascade_gain = 0 every envelope must instead collapse
113+% onto the closed-form exp(-zeta om t) of an uncoupled damped oscillator.
114+Q = P;
115+Q.n_diameters = 6; Q.max_mode_hz = 1500; Q.cascade_leak = 0;
116+bz = b6;
117+bz.zeta = zeros(bz.n, 1);
118+[tstr, ~, jstr, taustr, ~, sstr] = cymbal_strikes(Q, bz);
119+[envC, ~, DVc] = cymbal_cascade(Q, bz, tstr, jstr, taustr, sstr, 0.4);
120+Etot = sum(envC.^2 / 2, 1);
121+eCons = max(abs(Etot / Etot(1) - 1));
122+Q.cascade_gain = 0;
123+[envL, tfrL] = cymbal_cascade(Q, b6, tstr, jstr, taustr, sstr, 0.4);
124+pred = envL(:, 1) * exp(-(b6.zeta .* b6.omega).' * 0 + 0);
125+pred = (envL(:, 1) * ones(1, length(tfrL))) .* exp(-(b6.zeta .* b6.omega) * tfrL);
126+sel = envL(:, 1) > 1e-12;
127+eLin = max(max(abs(envL(sel, :) ./ pred(sel, :) - 1)));
128+ok = eCons < 1e-10 && eLin < 1e-9;
129+fails = fails + (~ok);
130+fprintf('cascade : energy drift %.1e, linear limit err %.1e %s\n', eCons, eLin, cv_verdict(ok));
131+
132+% ---- 7. far-field falloff: -6 dB per doubling of distance ----
133+% RMS over the render, with the two listeners on the same ray from the origin
134+% so every element dipole keeps its directivity factor.
135+Q = P;
136+Q.normalize = 0; Q.n_diameters = 12; Q.max_mode_hz = 3000;
137+Q.tail_s = 0.4; Q.oversample = 1;
138+Q.listener_x_m = 0.8; Q.listener_y_m = 0.5; Q.listener_z_m = 1.8;
139+y1 = cymbal_render(Q);
140+Q.listener_x_m = 1.6; Q.listener_y_m = 1.0; Q.listener_z_m = 3.6;
141+y2 = cymbal_render(Q);
142+drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
143+ok = abs(drop + 6.02) < 0.25;
144+fails = fails + (~ok);
145+fprintf('falloff : 2 m -> 4 m %.2f dB rms (target -6.02) %s\n', drop, cv_verdict(ok));
146+
147+% ---- 8. the seeded ride figure is exactly reproducible ----
148+Q = P;
149+Q.strike_count = 8;
150+[t1, v1, j1, ~, ~, s1, g1] = cymbal_strikes(Q, b6);
151+[t2, v2, ~, ~, ~, s2] = cymbal_strikes(Q, b6);
152+same = max(abs(t1 - t2)) == 0 && max(abs(v1 - v2)) == 0 && max(abs(s1 - s2)) == 0;
153+ok = same ...
154+ && g1.count == 8 ...
155+ && abs(g1.meanV - 3.111518886) < 1e-6 ...
156+ && abs(sum(j1) - 0.150818144) < 1e-9 ...
157+ && abs(t1(2) - 0.159427523) < 1e-9;
158+fails = fails + (~ok);
159+fprintf('ride (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\n', ...
160+ Q.seed, g1.count, g1.meanV, sum(j1), cv_verdict(ok));
161+
162+% ---- 9. pinned end-to-end diagnostics at the defaults ----
163+[y, basis, dinfo] = cymbal_render(P);
164+ok = abs(dinfo.peakPa - 13.4297) < 0.02 ...
165+ && abs(basis.f(1) - 49.0808) < 0.02 ...
166+ && basis.n == 303;
167+fails = fails + (~ok);
168+fprintf('pinned : peak %.4f Pa (target %.4f), mode1 %.3f Hz (target %.3f), %d modes %s\n', ...
169+ dinfo.peakPa, 13.4297, basis.f(1), 49.0808, basis.n, cv_verdict(ok));
170+
171+% ---- 10. the wash decays ----
172+% A cymbal rings for seconds, so the tail is not silent, but the energy in the
173+% last half second must sit well below that of the opening half second.
174+fs = P.sample_rate;
175+n = length(y);
176+head = y(1:round(0.5 * fs));
177+tail = y(n - round(0.5 * fs):n);
178+dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
179+ok = dropE < -20 && dropE > -80;
180+fails = fails + (~ok);
181+fprintf('decay : last/first half second %.1f dB %s\n', dropE, cv_verdict(ok));
182+
183+if fails > 0
184+ fprintf('%d CHECK(S) FAILED\n', fails);
185+else
186+ fprintf('all checks passed\n');
187+end
188+
189+function f = cv_freqs(G, n, P, Lref)
190+[K, M] = cymbal_shell(G, n, P.youngs_modulus, P.poisson_ratio, P.density, Lref);
191+Rc = chol(M);
192+Am = (Rc') \ (K / Rc);
193+Am = (Am + Am') / 2;
194+[~, Sv, ~] = svd(Am);
195+f = sort(sqrt(diag(Sv))) / (2 * pi);
196+end
197+
198+function lam = cv_plate_roots(n, nu)
199+% Roots of the free-edge frequency equation of a circular plate, bracketed on
200+% a coarse scan and bisected.
201+lam = [];
202+prev = cv_plate_det(n, nu, 0.30);
203+for x = 0.30:0.01:26
204+ cur = cv_plate_det(n, nu, x + 0.01);
205+ if prev * cur < 0
206+ lo = x;
207+ hi = x + 0.01;
208+ for it = 1:60
209+ mid = (lo + hi) / 2;
210+ if cv_plate_det(n, nu, lo) * cv_plate_det(n, nu, mid) < 0
211+ hi = mid;
212+ else
213+ lo = mid;
214+ end
215+ end
216+ lam = [lam, (lo + hi) / 2];
217+ end
218+ prev = cur;
219+end
220+end
221+
222+function d = cv_plate_det(n, nu, x)
223+% W = A J_n(lam r/a) + B I_n(lam r/a); the free edge sets the radial moment
224+% and the Kelvin-Kirchhoff shear to zero, and the 2 x 2 determinant vanishes.
225+J = besselj(n, x); I = besseli(n, x);
226+Jp = besselj(n - 1, x) - (n / x) * J;
227+Ip = besseli(n - 1, x) - (n / x) * I;
228+Jpp = -Jp / x + (n^2 / x^2 - 1) * J;
229+Ipp = -Ip / x + (n^2 / x^2 + 1) * I;
230+a11 = Jpp + nu * (Jp / x - n^2 * J / x^2);
231+a12 = Ipp + nu * (Ip / x - n^2 * I / x^2);
232+a21 = -Jp - (1 - nu) * (n^2 / x^2) * (Jp - J / x);
233+a22 = Ip - (1 - nu) * (n^2 / x^2) * (Ip - I / x);
234+d = a11 * a22 - a12 * a21;
235+end
236+
237+function [pd, pb] = cv_farfield(P, b, m, R, nth, tho, pho)
238+% Far-field pressure of one mode, per unit modal velocity, computed by direct
239+% summation over a theta-discretized dipole sheet and by the Bessel reduction.
240+c = P.speed_of_sound;
241+rho0 = P.air_density;
242+om = b.omega(m);
243+k = om / c;
244+n = b.ndia(m);
245+ld = b.ldip(m);
246+rh = [sin(tho) * cos(pho), sin(tho) * sin(pho), cos(tho)];
247+th = ((0:nth - 1) + 0.5) / nth * 2 * pi;
248+dth = 2 * pi / nth;
249+acc = 0;
250+for e = 1:length(b.surfR)
251+ x = b.surfR(e) * cos(th);
252+ y = b.surfR(e) * sin(th);
253+ z = b.surfZ(e) * ones(1, nth);
254+ nx = -b.surfDz(e) * cos(th);
255+ ny = -b.surfDz(e) * sin(th);
256+ nz = b.surfDr(e) * ones(1, nth);
257+ dx = R * rh(1) - x;
258+ dy = R * rh(2) - y;
259+ dz = R * rh(3) - z;
260+ rr = sqrt(dx.^2 + dy.^2 + dz.^2);
261+ proj = (nx .* dx + ny .* dy + nz .* dz) ./ rr * b.SW(e, m) .* cos(n * th);
262+ V = ld * b.surfR(e) * b.surfLen(e) * dth;
263+ acc = acc + sum((1i * om * rho0 / (4 * pi)) * V * proj .* (1i * k + 1 ./ rr) .* exp(-1i * k * rr) ./ rr);
264+end
265+pd = abs(acc);
266+mu = cos(tho);
267+sth = sin(tho);
268+al = k * b.surfR * sth;
269+Jn = besselj(n, al);
270+Jp1 = besselj(n + 1, al);
271+if n == 0
272+ Jm1 = -Jp1;
273+else
274+ Jm1 = besselj(n - 1, al);
275+end
276+ip = 1i^n;
277+Tk = ip * (2 * pi * Jn .* (b.surfDr * mu) - 1i * pi * (Jp1 - Jm1) .* (b.surfDz * sth));
278+An = sum(ld * b.surfR .* b.surfLen .* b.SW(:, m) .* Tk .* exp(1i * (k * b.surfZ) * mu));
279+pb = abs(-(om^2 * rho0 / (4 * pi * c * R)) * exp(-1i * k * R) * An * cos(n * pho));
280+end
281+
282+function s = cv_verdict(ok)
283+if ok
284+ s = 'ok';
285+else
286+ s = 'MISMATCH';
287+end
288+end
webapp/.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+node_modules/
2+dist/
webapp/README.mdadded+50−0View file
@@ -0,0 +1,50 @@
1+# Cymbal 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 stick contact (or a seeded repeated figure, with `Strokes` above 1), the modes of a bronze shell of revolution from one finite element eigensolve per nodal-diameter order, a nonlinear cascade of modal energy up the frequency ladder, an oscillator bank radiating through a coherent sum of surface 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/cymbal_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 cymbal_demo.m` (or `cymbal_validate.m`).
17+
18+## Where the time goes
19+
20+The expensive stage is the mode basis: one generalized eigensolve per nodal-diameter order, about two seconds at the defaults, and that is exactly the piece the worker caches. Striking somewhere else, changing the stick, moving the microphone or retuning the cascade all reuse the cached modes. The synthesis that follows is an oscillator bank rather than the closed-form exponentials of the sibling `triangle` repository, because the cascade makes every amplitude time-varying; the per-frame phase advance is a unit-modulus multiply and the oscillator basis is hoisted out of the loop, so a three-second render costs about the same again.
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. Zeroing the meridian curvature `cs` in `cymbal_shell.m`, for instance, decouples bending from stretching and collapses the partials back onto those of a flat disc.
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: plate geometry, material, damping, radiation scaling and discretization parameters invalidate the mode basis; the stick, the gesture, the cascade and the listener only re-run the strike process, the cascade and the radiation; pre-delay, normalization and bit depth only re-run the output stage. Intermediate tensors stay inside the worker as numbl RuntimeValues.
31+
32+The New seed button redraws the timing and velocities of a repeated figure without touching any other parameter (a single stroke 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/cymbal_defaults.m` mirrors).
33+
34+## Layout
35+
36+- `../matlab/` — the algorithm: `cymbal_profile.m`, `cymbal_shell.m`, `cymbal_basis.m`, `cymbal_at.m`, `cymbal_point.m`, `cymbal_strikes.m`, `cymbal_rng.m`, `cymbal_cascade.m`, `cymbal_radiate.m`, `cymbal_output.m`, plus `cymbal_render.m` / `cymbal_demo.m` / `cymbal_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/profile.ts` — the cymbal meridian, mirrored from `cymbal_profile.m` for the scene
41+- `src/components/` — 3D scene (three.js), parameter panel, code panel and editor, plots (strokes, cascade, waveform, spectrum, spectrogram, mode map, mode shapes)
42+- `src/selftest.ts` — Node harness for the MATLAB validation script
43+
44+## The cascade plot
45+
46+The plot the sibling repositories do not have. The cascade solve tracks the energy of every mode at the control rate, and the plot collapses that onto a log-frequency by time grid, so the upward flux is visible directly. Dropping `Cascade gain` to zero flattens it into horizontal streaks, every mode simply decaying where the stick left it, and the sound turns from a crash into a gong.
47+
48+## A note on the seeded generator
49+
50+The stroke pattern draws from a 32-bit LCG in `cymbal_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 stroke 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>Cymbal 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": "cymbal-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+289−0View file
@@ -0,0 +1,289 @@
1+import { useCallback, useEffect, useRef, useState } from 'react';
2+import { DEFAULTS, loadParams, saveParams, type CymbalParams } 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 ModeGallery from './components/plots/ModeGallery';
16+import CascadePlot from './components/plots/CascadePlot';
17+
18+export default function App() {
19+ const [params, setParams] = useState<CymbalParams>(() => loadParams());
20+ // Editable copy of the MATLAB implementation. Deliberately not persisted:
21+ // a reload restores the code the app shipped with.
22+ const [mfiles, setMfiles] = useState<MFile[]>(() => DEFAULT_MFILES.map((f) => ({ ...f })));
23+ const [codeOpen, setCodeOpen] = useState(false);
24+ const { result, error, busy, stage, progress, renderFull } = useRenderer(params, mfiles);
25+ const audioCtxRef = useRef<AudioContext | null>(null);
26+
27+ useEffect(() => saveParams(params), [params]);
28+
29+ const onChange = useCallback((key: keyof CymbalParams, value: number | boolean) => {
30+ setParams((p) => ({ ...p, [key]: value }));
31+ }, []);
32+
33+ const reset = useCallback(() => setParams({ ...DEFAULTS }), []);
34+
35+ // A new seed redraws the timing and velocities of a repeated figure without
36+ // touching anything else. A single stroke is deterministic, so the button is
37+ // disabled there.
38+ const reseed = useCallback(
39+ () => setParams((p) => ({ ...p, seed: 1 + Math.floor(Math.random() * 9999) })),
40+ [],
41+ );
42+
43+ const onChangeFile = useCallback((name: string, source: string) => {
44+ setMfiles((fs) => fs.map((f) => (f.name === name ? { ...f, source } : f)));
45+ }, []);
46+
47+ const onResetFile = useCallback((name: string) => {
48+ setMfiles((fs) =>
49+ fs.map((f) =>
50+ f.name === name
51+ ? { ...f, source: DEFAULT_MFILES.find((d) => d.name === name)?.source ?? f.source }
52+ : f,
53+ ),
54+ );
55+ }, []);
56+
57+ const onResetAllFiles = useCallback(() => setMfiles(DEFAULT_MFILES.map((f) => ({ ...f }))), []);
58+ const codeEdited = mfiles.some(
59+ (f) => f.source !== DEFAULT_MFILES.find((d) => d.name === f.name)?.source,
60+ );
61+
62+ const [looping, setLooping] = useState(false);
63+ const loopSrcRef = useRef<AudioBufferSourceNode | null>(null);
64+
65+ const getAudioCtx = useCallback((sampleRate: number) => {
66+ let ctx = audioCtxRef.current;
67+ if (!ctx || ctx.state === 'closed') {
68+ ctx = new AudioContext({ sampleRate });
69+ audioCtxRef.current = ctx;
70+ }
71+ void ctx.resume();
72+ return ctx;
73+ }, []);
74+
75+ const makeBuffer = useCallback((ctx: AudioContext, r: NonNullable<typeof result>) => {
76+ const buf = ctx.createBuffer(1, r.y.length, r.sampleRate);
77+ buf.copyToChannel(new Float32Array(r.y), 0);
78+ return buf;
79+ }, []);
80+
81+ const play = useCallback(() => {
82+ if (!result) return;
83+ const ctx = getAudioCtx(result.sampleRate);
84+ const src = ctx.createBufferSource();
85+ src.buffer = makeBuffer(ctx, result);
86+ src.connect(ctx.destination);
87+ src.start();
88+ }, [result, getAudioCtx, makeBuffer]);
89+
90+ // Loop mode: repeat the stroke continuously, swapping in each new render as
91+ // parameters change so edits are audible without pressing Play.
92+ useEffect(() => {
93+ if (!looping || !result) return;
94+ const ctx = getAudioCtx(result.sampleRate);
95+ const src = ctx.createBufferSource();
96+ src.buffer = makeBuffer(ctx, result);
97+ src.loop = true;
98+ src.connect(ctx.destination);
99+ src.start();
100+ loopSrcRef.current = src;
101+ return () => {
102+ loopSrcRef.current = null;
103+ try {
104+ src.stop();
105+ } catch {
106+ // already stopped
107+ }
108+ src.disconnect();
109+ };
110+ }, [looping, result, getAudioCtx, makeBuffer]);
111+
112+ const download = useCallback(() => {
113+ if (!result) return;
114+ const blob = encodeWav(result.y, result.sampleRate, params.bit_depth === 16 ? 16 : 24);
115+ const url = URL.createObjectURL(blob);
116+ const a = document.createElement('a');
117+ a.href = url;
118+ a.download = result.quality === 'full' ? 'cymbal.wav' : 'cymbal_preview.wav';
119+ a.click();
120+ setTimeout(() => URL.revokeObjectURL(url), 5000);
121+ }, [result, params.bit_depth]);
122+
123+ return (
124+ <div className="app">
125+ <div className="topbar">
126+ <h1>Cymbal Simulator</h1>
127+ <button className="primary" onClick={renderFull} disabled={busy}>
128+ Render
129+ </button>
130+ <button
131+ className={`play ${result?.quality === 'preview' ? 'preview' : ''}`}
132+ onClick={play}
133+ disabled={!result}
134+ title={
135+ result?.quality === 'preview'
136+ ? 'Playing the fast preview render; press Render for full quality'
137+ : 'Play the rendered stroke'
138+ }
139+ >
140+ ▶ Play{result?.quality === 'preview' ? ' (preview)' : ''}
141+ </button>
142+ <button
143+ className={`loop ${looping ? 'active' : ''}`}
144+ onClick={() => setLooping((v) => !v)}
145+ disabled={!result && !looping}
146+ title="Repeat the render continuously so parameter changes are audible immediately"
147+ >
148+ ⟳ Loop{looping ? ' on' : ''}
149+ </button>
150+ <button
151+ onClick={reseed}
152+ disabled={params.strike_count < 2}
153+ title="Draw new timing and velocities from a new random seed (repeated figures only)"
154+ >
155+ ⚄ New seed
156+ </button>
157+ <button onClick={download} disabled={!result}>
158+ Download .wav{result?.quality === 'preview' ? ' (preview)' : ''}
159+ </button>
160+ <button onClick={reset}>Reset to defaults</button>
161+ <button
162+ className={codeOpen ? 'code active' : 'code'}
163+ onClick={() => setCodeOpen((v) => !v)}
164+ title="View and edit the MATLAB source that numbl executes"
165+ >
166+ {'{ }'} Code{codeEdited ? ' •' : ''}
167+ </button>
168+ <div className="status">
169+ {error && <span className="badge err">code error</span>}
170+ {result && (
171+ <span className={`badge ${result.quality}`}>
172+ {result.quality === 'full' ? 'full quality' : 'preview quality — press Render'}
173+ </span>
174+ )}
175+ {busy ? (
176+ <>
177+ <span>{stage}…</span>
178+ <div className="progressbar">
179+ <div style={{ width: `${Math.round(progress * 100)}%` }} />
180+ </div>
181+ </>
182+ ) : (
183+ <span>ready</span>
184+ )}
185+ </div>
186+ </div>
187+ <div className="main">
188+ <div className="sidebar">
189+ <ParamPanel params={params} onChange={onChange} />
190+ </div>
191+ <div className="content">
192+ {error && (
193+ <div className="card wide error-card">
194+ <h2>MATLAB error</h2>
195+ <div>
196+ <strong>
197+ {error.file ? `${error.file}${error.line ? `:${error.line}` : ''}` : 'error'}
198+ </strong>{' '}
199+ {error.message}
200+ </div>
201+ <div className="code-hint">
202+ Showing the last successful render. Open Code to fix the source, or revert it.
203+ </div>
204+ </div>
205+ )}
206+ <div className="card wide">
207+ <h2>Geometry — plate, stick, microphone</h2>
208+ <Scene3D params={params} />
209+ </div>
210+ {result && (
211+ <>
212+ <div className="card">
213+ <h2>Diagnostics</h2>
214+ <DiagnosticsPanel result={result} />
215+ </div>
216+ <div className="card">
217+ <h2>Strokes</h2>
218+ <StrikePlot
219+ strikeT={result.strikeT}
220+ strikeV={result.strikeV}
221+ strikeS={result.strikeS}
222+ durationS={result.diag.durationS}
223+ />
224+ </div>
225+ <div className="card wide">
226+ <h2>Modal energy — the cascade</h2>
227+ <CascadePlot
228+ cascadeDb={result.cascadeDb}
229+ bands={result.cascadeBands}
230+ frames={result.cascadeFrames}
231+ fLo={result.cascadeFLo}
232+ fHi={result.cascadeFHi}
233+ tEnd={result.cascadeTEnd}
234+ />
235+ </div>
236+ <div className="card">
237+ <h2>Waveform</h2>
238+ <WaveformPlot y={result.y} sampleRate={result.sampleRate} />
239+ </div>
240+ <div className="card">
241+ <h2>Spectrum</h2>
242+ <SpectrumPlot
243+ y={result.y}
244+ sampleRate={result.sampleRate}
245+ modeF={result.modeF}
246+ modeN={result.modeN}
247+ />
248+ </div>
249+ <div className="card">
250+ <h2>Spectrogram</h2>
251+ <SpectrogramPlot y={result.y} sampleRate={result.sampleRate} />
252+ </div>
253+ <div className="card">
254+ <h2>Plate modes</h2>
255+ <ModeMap
256+ modeF={result.modeF}
257+ modeZeta={result.modeZeta}
258+ modeZetaRad={result.modeZetaRad}
259+ modeN={result.modeN}
260+ />
261+ </div>
262+ <div className="card wide">
263+ <h2>Mode shapes</h2>
264+ <ModeGallery
265+ nodeR={result.nodeR}
266+ shapeW={result.shapeW}
267+ modeF={result.modeF}
268+ modeN={result.modeN}
269+ />
270+ </div>
271+ </>
272+ )}
273+ </div>
274+ </div>
275+ {codeOpen && (
276+ <CodePanel
277+ mfiles={mfiles}
278+ onChangeFile={onChangeFile}
279+ onResetFile={onResetFile}
280+ onResetAll={onResetAllFiles}
281+ onClose={() => setCodeOpen(false)}
282+ onRenderFull={renderFull}
283+ error={error}
284+ busy={busy}
285+ />
286+ )}
287+ </div>
288+ );
289+}
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+71−0View file
@@ -0,0 +1,71 @@
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, 8);
7+ return (
8+ <div className="diag">
9+ <table>
10+ <tbody>
11+ <tr>
12+ <th>Strokes</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>Plate modes</th>
27+ <td className="partials">
28+ {d.modesKept} up to n = {d.maxN} nodal diameters · {d.plateMassKg.toFixed(2)} kg
29+ </td>
30+ </tr>
31+ <tr>
32+ <th>Lowest</th>
33+ <td className="partials">{shown.map((f) => f.toFixed(0)).join(', ')} … Hz</td>
34+ </tr>
35+ <tr>
36+ <th>Nonlinearity</th>
37+ <td>
38+ peak w<sub>rms</sub> {d.peakWrmsUm.toFixed(0)} µm ={' '}
39+ {Math.sqrt(d.peakStrength).toFixed(2)} thicknesses ·{' '}
40+ {(100 * d.cascadeLostFrac).toFixed(1)}% left the top of the band
41+ </td>
42+ </tr>
43+ <tr>
44+ <th>Listener</th>
45+ <td>{(d.listenerDistM * 100).toFixed(0)} cm from the plate</td>
46+ </tr>
47+ <tr>
48+ <th>Peak pressure</th>
49+ <td>
50+ {d.peakPa.toFixed(3)} Pa ({d.peakSplDb.toFixed(1)} dB SPL) · click{' '}
51+ {d.directPeakPa.toFixed(2)} Pa
52+ </td>
53+ </tr>
54+ <tr>
55+ <th>Length</th>
56+ <td>
57+ {d.durationS.toFixed(2)} s · {result.quality} · oversample {d.oversample}
58+ </td>
59+ </tr>
60+ </tbody>
61+ </table>
62+ <div className="cachenote">
63+ eigensolves: {d.solveCached ? 'cached' : `${(d.solveMs / 1000).toFixed(2)} s`} · cascade +
64+ radiation: {d.radiateCached ? 'cached' : `${(d.radiateMs / 1000).toFixed(2)} s`}
65+ </div>
66+ {d.warnings.map((w, i) => (
67+ <div key={i} className="warning">⚠ {w}</div>
68+ ))}
69+ </div>
70+ );
71+}
webapp/src/components/ParamPanel.tsxadded+123−0View file
@@ -0,0 +1,123 @@
1+import { PARAM_GROUPS, type ParamSpec, type CymbalParams } from '../params';
2+
3+interface Props {
4+ params: CymbalParams;
5+ onChange: (key: keyof CymbalParams, 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+312−0View file
@@ -0,0 +1,312 @@
1+// 3D view of the geometry: the cymbal on its stand, the stick at the contact
2+// point (and the second contact radius of a repeated figure), and the
3+// listener. Model coordinates put the symmetry axis along z with the rim at
4+// z = 0 and the bell rising toward +z, and the stroke on the theta = 0
5+// meridian, which is +x. three.js is y-up, so the scene rotates the model by
6+// -90 degrees about x and the model's z becomes the scene's y.
7+
8+import { useEffect, useRef } from 'react';
9+import * as THREE from 'three';
10+import { OrbitControls } from 'three/examples/jsm/controls/OrbitControls.js';
11+import type { CymbalParams } from '../params';
12+import { cymbalMeridian } from '../lib/profile';
13+
14+function makeLabel(text: string, color: string, scale = 1): THREE.Sprite {
15+ const canvas = document.createElement('canvas');
16+ canvas.width = 320;
17+ canvas.height = 64;
18+ const ctx = canvas.getContext('2d')!;
19+ ctx.font = '600 36px system-ui, sans-serif';
20+ ctx.fillStyle = color;
21+ ctx.textAlign = 'center';
22+ ctx.textBaseline = 'middle';
23+ ctx.fillText(text, 160, 32);
24+ const tex = new THREE.CanvasTexture(canvas);
25+ const mat = new THREE.SpriteMaterial({ map: tex, depthTest: false, transparent: true });
26+ const sprite = new THREE.Sprite(mat);
27+ sprite.scale.set(0.11 * scale, 0.022 * scale, 1);
28+ return sprite;
29+}
30+
31+// model (r along +x at theta = 0, z along the axis) -> scene (y up)
32+function toScene(x: number, y: number, z: number): THREE.Vector3 {
33+ return new THREE.Vector3(x, z, -y);
34+}
35+
36+interface Props {
37+ params: CymbalParams;
38+}
39+
40+export default function Scene3D({ params }: Props) {
41+ const mountRef = useRef<HTMLDivElement | null>(null);
42+ const sceneRef = useRef<{
43+ scene: THREE.Scene;
44+ group: THREE.Group | null;
45+ renderer: THREE.WebGLRenderer;
46+ camera: THREE.PerspectiveCamera;
47+ } | null>(null);
48+
49+ useEffect(() => {
50+ const mount = mountRef.current!;
51+ const renderer = new THREE.WebGLRenderer({ antialias: true });
52+ renderer.setPixelRatio(window.devicePixelRatio);
53+ renderer.setClearColor(0x14161a);
54+ mount.appendChild(renderer.domElement);
55+
56+ const scene = new THREE.Scene();
57+ const camera = new THREE.PerspectiveCamera(45, 1, 0.01, 30);
58+ camera.position.set(0.42, 0.36, 0.55);
59+
60+ const controls = new OrbitControls(camera, renderer.domElement);
61+ controls.target.set(0, 0.03, 0);
62+ controls.enableDamping = true;
63+
64+ scene.add(new THREE.AmbientLight(0xffffff, 0.5));
65+ const dir = new THREE.DirectionalLight(0xffffff, 1.2);
66+ dir.position.set(0.6, 1.0, 0.5);
67+ scene.add(dir);
68+ const dir2 = new THREE.DirectionalLight(0x8899ff, 0.35);
69+ dir2.position.set(-0.6, 0.4, -0.6);
70+ scene.add(dir2);
71+
72+ const grid = new THREE.GridHelper(1.6, 32, 0x3a4048, 0x262b33);
73+ grid.position.y = -0.36;
74+ scene.add(grid);
75+
76+ // model axes: x red (the stroke meridian), y green (the other in-plane
77+ // direction), z blue (the symmetry axis)
78+ const axisLen = 0.09;
79+ const axes: [THREE.Vector3, number, string][] = [
80+ [toScene(axisLen, 0, 0), 0xe57373, '+x'],
81+ [toScene(0, axisLen, 0), 0x81c784, '+y'],
82+ [toScene(0, 0, axisLen), 0x64b5f6, '+z axis'],
83+ ];
84+ for (const [v, color, label] of axes) {
85+ const geo = new THREE.BufferGeometry().setFromPoints([new THREE.Vector3(0, 0, 0), v]);
86+ scene.add(new THREE.Line(geo, new THREE.LineBasicMaterial({ color })));
87+ const s = makeLabel(label, '#' + color.toString(16).padStart(6, '0'));
88+ s.position.copy(v).multiplyScalar(1.2);
89+ scene.add(s);
90+ }
91+
92+ sceneRef.current = { scene, group: null, renderer, camera };
93+
94+ let raf = 0;
95+ const animate = () => {
96+ raf = requestAnimationFrame(animate);
97+ controls.update();
98+ renderer.render(scene, camera);
99+ };
100+ animate();
101+
102+ const resize = () => {
103+ const w = mount.clientWidth;
104+ const h = mount.clientHeight;
105+ renderer.setSize(w, h, false);
106+ camera.aspect = w / h;
107+ camera.updateProjectionMatrix();
108+ };
109+ resize();
110+ const ro = new ResizeObserver(resize);
111+ ro.observe(mount);
112+
113+ return () => {
114+ cancelAnimationFrame(raf);
115+ ro.disconnect();
116+ controls.dispose();
117+ renderer.dispose();
118+ mount.removeChild(renderer.domElement);
119+ sceneRef.current = null;
120+ };
121+ }, []);
122+
123+ useEffect(() => {
124+ const ctx = sceneRef.current;
125+ if (!ctx) return;
126+ if (ctx.group) {
127+ ctx.scene.remove(ctx.group);
128+ ctx.group.traverse((o) => {
129+ const mesh = o as THREE.Mesh;
130+ if (mesh.geometry) mesh.geometry.dispose();
131+ });
132+ }
133+
134+ const g = new THREE.Group();
135+ const mer = cymbalMeridian(
136+ params.radius_m,
137+ params.hole_radius_m,
138+ params.bell_radius_m,
139+ params.bell_height_m,
140+ params.bow_height_m,
141+ params.thickness_bell_m,
142+ params.thickness_rim_m,
143+ );
144+
145+ // meridian sample points, used both for the lathe and for placing the stick
146+ const NS = 64;
147+ const rs: number[] = [];
148+ for (let i = 0; i <= NS; i++) rs.push(mer.rHole + (i / NS) * (mer.rRim - mer.rHole));
149+ // arclength, so that strike_r addresses the same point the model does
150+ const arc: number[] = [0];
151+ for (let i = 1; i <= NS; i++) {
152+ arc.push(arc[i - 1] + Math.hypot(rs[i] - rs[i - 1], mer.z(rs[i]) - mer.z(rs[i - 1])));
153+ }
154+ const total = arc[NS];
155+ const atFrac = (frac: number) => {
156+ const s = Math.min(0.9999, Math.max(0, frac)) * total;
157+ let i = 0;
158+ while (i < NS - 1 && arc[i + 1] < s) i++;
159+ const t = (s - arc[i]) / Math.max(1e-12, arc[i + 1] - arc[i]);
160+ const r = rs[i] + t * (rs[i + 1] - rs[i]);
161+ const z = mer.z(rs[i]) + t * (mer.z(rs[i + 1]) - mer.z(rs[i]));
162+ const dr = rs[i + 1] - rs[i];
163+ const dz = mer.z(rs[i + 1]) - mer.z(rs[i]);
164+ const dl = Math.hypot(dr, dz) || 1;
165+ return { r, z, tr: dr / dl, tz: dz / dl };
166+ };
167+
168+ // the plate: a lathe of the meridian offset to both faces
169+ const pts: THREE.Vector2[] = [];
170+ for (let i = 0; i <= NS; i++) {
171+ const h = mer.thickness(rs[i]) / 2;
172+ pts.push(new THREE.Vector2(rs[i], mer.z(rs[i]) + h));
173+ }
174+ for (let i = NS; i >= 0; i--) {
175+ const h = mer.thickness(rs[i]) / 2;
176+ pts.push(new THREE.Vector2(rs[i], mer.z(rs[i]) - h));
177+ }
178+ pts.push(pts[0].clone());
179+ const lathe = new THREE.Mesh(
180+ new THREE.LatheGeometry(pts, 96),
181+ new THREE.MeshStandardMaterial({
182+ color: 0xc9a227,
183+ roughness: 0.32,
184+ metalness: 0.9,
185+ side: THREE.DoubleSide,
186+ }),
187+ );
188+ lathe.rotation.x = -Math.PI / 2;
189+ g.add(lathe);
190+
191+ // lathing grooves, as faint rings
192+ for (let k = 1; k < 9; k++) {
193+ const r = mer.rBell + (k / 9) * (mer.rRim - mer.rBell);
194+ const ringPts: THREE.Vector3[] = [];
195+ for (let a = 0; a <= 72; a++) {
196+ const th = (a / 72) * 2 * Math.PI;
197+ ringPts.push(toScene(r * Math.cos(th), r * Math.sin(th), mer.z(r) + 0.0008));
198+ }
199+ const ring = new THREE.Line(
200+ new THREE.BufferGeometry().setFromPoints(ringPts),
201+ new THREE.LineBasicMaterial({ color: 0x8a6d18, transparent: true, opacity: 0.5 }),
202+ );
203+ g.add(ring);
204+ }
205+
206+ // stand: rod through the hole, felt washer, tripod post
207+ const rod = new THREE.Mesh(
208+ new THREE.CylinderGeometry(0.0035, 0.0035, 0.09, 12),
209+ new THREE.MeshStandardMaterial({ color: 0x9aa0aa, roughness: 0.4, metalness: 0.8 }),
210+ );
211+ rod.position.set(0, mer.apexHeight + 0.01, 0);
212+ g.add(rod);
213+ const felt = new THREE.Mesh(
214+ new THREE.CylinderGeometry(0.012, 0.012, 0.006, 16),
215+ new THREE.MeshStandardMaterial({ color: 0x2f3a44, roughness: 0.95 }),
216+ );
217+ felt.position.set(0, mer.apexHeight - 0.004, 0);
218+ g.add(felt);
219+ const post = new THREE.Mesh(
220+ new THREE.CylinderGeometry(0.008, 0.008, 0.36, 12),
221+ new THREE.MeshStandardMaterial({ color: 0x6f757f, roughness: 0.5, metalness: 0.7 }),
222+ );
223+ post.position.set(0, mer.apexHeight - 0.19, 0);
224+ g.add(post);
225+
226+ // stick at the contact point, along the stroke direction
227+ const addStick = (frac: number, primary: boolean) => {
228+ const pt = atFrac(frac);
229+ const ang = (params.strike_angle_deg * Math.PI) / 180;
230+ // surface normal (-dz, dr) in the meridian plane, at theta = 0
231+ const nrm = toScene(-pt.tz, 0, pt.tr);
232+ const tan = toScene(pt.tr, 0, pt.tz);
233+ // stroke direction: -cos(angle) * normal + sin(angle) * tangent
234+ const dirv = nrm.clone().multiplyScalar(-Math.cos(ang)).addScaledVector(tan, Math.sin(ang));
235+ const contact = toScene(pt.r, 0, pt.z);
236+ const tip = contact.clone().addScaledVector(dirv, -(params.tip_radius_m + 0.0008));
237+ const color = primary ? 0xe8a33d : 0x9a7a3d;
238+ const head = new THREE.Mesh(
239+ new THREE.SphereGeometry(params.tip_radius_m, 16, 12),
240+ new THREE.MeshStandardMaterial({ color, roughness: 0.35 }),
241+ );
242+ head.position.copy(tip);
243+ g.add(head);
244+ const shaftLen = 0.24;
245+ const shaft = new THREE.Mesh(
246+ new THREE.CylinderGeometry(0.0055, 0.0072, shaftLen, 12),
247+ new THREE.MeshStandardMaterial({ color: 0xb08c5a, roughness: 0.75 }),
248+ );
249+ shaft.quaternion.setFromUnitVectors(new THREE.Vector3(0, 1, 0), dirv.clone().negate());
250+ shaft.position.copy(tip).addScaledVector(dirv, -(params.tip_radius_m + shaftLen / 2));
251+ g.add(shaft);
252+ if (primary) {
253+ const lbl = makeLabel('stroke', '#e8a33d');
254+ lbl.position.copy(tip).add(new THREE.Vector3(0, 0.05, 0));
255+ g.add(lbl);
256+ }
257+ };
258+ addStick(params.strike_r, true);
259+ if (params.strike_count > 1) addStick(params.roll_r2, false);
260+
261+ // listener
262+ const lis = toScene(params.listener_x_m, params.listener_y_m, params.listener_z_m);
263+ const micHead = new THREE.Mesh(
264+ new THREE.SphereGeometry(0.022, 20, 16),
265+ new THREE.MeshStandardMaterial({ color: 0x4dd0e1, emissive: 0x0c3b42, roughness: 0.35 }),
266+ );
267+ micHead.position.copy(lis);
268+ g.add(micHead);
269+ const micBody = new THREE.Mesh(
270+ new THREE.CylinderGeometry(0.011, 0.011, 0.07, 12),
271+ new THREE.MeshStandardMaterial({ color: 0x455a64, roughness: 0.6 }),
272+ );
273+ micBody.position.copy(lis).add(new THREE.Vector3(0, -0.046, 0));
274+ g.add(micBody);
275+ const centroid = toScene(0, 0, mer.z(mer.rRim * 0.6));
276+ const lineGeo = new THREE.BufferGeometry().setFromPoints([centroid, lis]);
277+ const line = new THREE.Line(
278+ lineGeo,
279+ new THREE.LineDashedMaterial({ color: 0x4dd0e1, dashSize: 0.02, gapSize: 0.014 }),
280+ );
281+ line.computeLineDistances();
282+ g.add(line);
283+ const micLabel = makeLabel(`mic ${(centroid.distanceTo(lis) * 100).toFixed(0)} cm`, '#4dd0e1');
284+ micLabel.position.copy(lis).add(new THREE.Vector3(0, 0.06, 0));
285+ g.add(micLabel);
286+
287+ ctx.scene.add(g);
288+ ctx.group = g;
289+ }, [
290+ params.radius_m,
291+ params.hole_radius_m,
292+ params.bell_radius_m,
293+ params.bell_height_m,
294+ params.bow_height_m,
295+ params.thickness_bell_m,
296+ params.thickness_rim_m,
297+ params.strike_r,
298+ params.strike_angle_deg,
299+ params.strike_count,
300+ params.roll_r2,
301+ params.tip_radius_m,
302+ params.listener_x_m,
303+ params.listener_y_m,
304+ params.listener_z_m,
305+ ]);
306+
307+ return (
308+ <div className="scene3d" ref={mountRef}>
309+ <div className="scene-hint">drag to orbit · scroll to zoom · right-drag to pan</div>
310+ </div>
311+ );
312+}
webapp/src/components/plots/CascadePlot.tsxadded+73−0View file
@@ -0,0 +1,73 @@
1+// Where the energy is, and where it goes. The cascade solve tracks the energy
2+// of every mode at the control rate; this collapses it onto a log-frequency by
3+// time grid, so the upward flux is visible directly as the bright region
4+// climbing out of the strike. Turning cascade_gain down to zero flattens it
5+// into horizontal streaks: every mode simply decaying where it was struck.
6+
7+import { useCanvas, drawFrame, PLOT_COLORS, colormap } from './canvasUtils';
8+
9+interface Props {
10+ cascadeDb: Float32Array; // nBands x nFrames, column-major
11+ bands: number;
12+ frames: number;
13+ fLo: number;
14+ fHi: number;
15+ tEnd: number;
16+}
17+
18+const DB_RANGE = 70;
19+
20+export default function CascadePlot({ cascadeDb, bands, frames, fLo, fHi, tEnd }: Props) {
21+ const ref = useCanvas(
22+ (ctx, w, h) => {
23+ drawFrame(ctx, w, h);
24+ const padL = 38;
25+ const padB = 16;
26+ const pw = w - padL - 8;
27+ const ph = h - padB - 8;
28+ if (bands < 1 || frames < 1) return;
29+
30+ let maxDb = -Infinity;
31+ for (let i = 0; i < cascadeDb.length; i++) maxDb = Math.max(maxDb, cascadeDb[i]);
32+ if (!Number.isFinite(maxDb)) return;
33+
34+ const img = ctx.createImageData(frames, bands);
35+ for (let f = 0; f < frames; f++) {
36+ for (let b = 0; b < bands; b++) {
37+ const v = cascadeDb[f * bands + b] - maxDb;
38+ const t = Math.min(1, Math.max(0, (v + DB_RANGE) / DB_RANGE));
39+ const [r, g, bl] = colormap(t);
40+ const row = bands - 1 - b;
41+ const idx = (row * frames + f) * 4;
42+ img.data[idx] = r;
43+ img.data[idx + 1] = g;
44+ img.data[idx + 2] = bl;
45+ img.data[idx + 3] = 255;
46+ }
47+ }
48+ const off = new OffscreenCanvas(frames, bands);
49+ off.getContext('2d')!.putImageData(img, 0, 0);
50+ ctx.imageSmoothingEnabled = true;
51+ ctx.drawImage(off, padL, 8, pw, ph);
52+
53+ ctx.fillStyle = PLOT_COLORS.axis;
54+ ctx.font = '10px system-ui';
55+ ctx.textAlign = 'right';
56+ const lLo = Math.log10(fLo);
57+ const lHi = Math.log10(Math.max(fLo * 1.2, fHi));
58+ for (const f of [50, 100, 200, 500, 1000, 2000, 5000, 10000]) {
59+ if (f < fLo || f > fHi) continue;
60+ const yy = 8 + ph - ((Math.log10(f) - lLo) / (lHi - lLo)) * ph;
61+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, padL - 4, yy + 3);
62+ }
63+ ctx.textAlign = 'center';
64+ const step = tEnd > 4 ? 1 : tEnd > 1.5 ? 0.5 : 0.1;
65+ for (let t = 0; t <= tEnd + 1e-9; t += step) {
66+ ctx.fillText(`${t.toFixed(1)}s`, padL + (t / tEnd) * pw, h - 4);
67+ }
68+ },
69+ [cascadeDb, bands, frames, fLo, fHi, tEnd],
70+ 200,
71+ );
72+ return <canvas ref={ref} />;
73+}
webapp/src/components/plots/ModeGallery.tsxadded+96−0View file
@@ -0,0 +1,96 @@
1+// The lowest mode shapes, drawn as the plate is seen from above: the normal
2+// displacement W(r) cos(n theta), blue toward the viewer and red away. The
3+// nodal diameters are the straight lines through the centre and the nodal
4+// circles are the rings, so a mode is read off directly as its (n, k) pair.
5+// Only the cos member of each degenerate pair is retained by the model (see
6+// cymbal_basis.m), and it is the member aligned with the stroke, drawn here at
7+// theta = 0 pointing right.
8+
9+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
10+
11+interface Props {
12+ nodeR: Float32Array;
13+ shapeW: Float32Array; // nn x n, column-major
14+ modeF: Float32Array;
15+ modeN: Float32Array;
16+}
17+
18+const SHOWN = 12;
19+const CELL = 118;
20+const COLS = 4;
21+
22+export default function ModeGallery({ nodeR, shapeW, modeF, modeN }: Props) {
23+ const nShown = Math.min(SHOWN, modeF.length);
24+ const rows = Math.max(1, Math.ceil(nShown / COLS));
25+ const ref = useCanvas(
26+ (ctx, w, h) => {
27+ drawFrame(ctx, w, h);
28+ const nn = nodeR.length;
29+ if (nn === 0 || nShown === 0) return;
30+ const cw = w / COLS;
31+ const ch = h / rows;
32+ const rMax = nodeR[nn - 1];
33+ const rMin = nodeR[0];
34+ const rad = Math.min(cw, ch) / 2 - 16;
35+
36+ for (let m = 0; m < nShown; m++) {
37+ const cx = (m % COLS) * cw + cw / 2;
38+ const cy = Math.floor(m / COLS) * ch + ch / 2 - 4;
39+ const off = m * nn;
40+ let uMax = 1e-30;
41+ for (let i = 0; i < nn; i++) uMax = Math.max(uMax, Math.abs(shapeW[off + i]));
42+ const nDia = modeN[m];
43+
44+ const size = Math.max(8, Math.round(2 * rad));
45+ const img = ctx.createImageData(size, size);
46+ for (let py = 0; py < size; py++) {
47+ for (let px = 0; px < size; px++) {
48+ const dx = (px - size / 2) / (size / 2);
49+ const dy = (py - size / 2) / (size / 2);
50+ const rr = Math.hypot(dx, dy);
51+ const idx = (py * size + px) * 4;
52+ if (rr > 1) {
53+ img.data[idx + 3] = 0;
54+ continue;
55+ }
56+ const rPhys = rMin + rr * (rMax - rMin);
57+ // linear interpolation of the nodal normal shape
58+ let j = Math.min(nn - 2, Math.max(0, Math.floor(((rPhys - rMin) / (rMax - rMin)) * (nn - 1))));
59+ const t = (rPhys - nodeR[j]) / Math.max(1e-12, nodeR[j + 1] - nodeR[j]);
60+ const wr = shapeW[off + j] + t * (shapeW[off + j + 1] - shapeW[off + j]);
61+ const th = Math.atan2(dy, dx);
62+ const v = (wr / uMax) * Math.cos(nDia * th);
63+ const a = Math.min(1, Math.abs(v));
64+ if (v >= 0) {
65+ img.data[idx] = Math.round(24 + 76 * a);
66+ img.data[idx + 1] = Math.round(30 + 151 * a);
67+ img.data[idx + 2] = Math.round(38 + 208 * a);
68+ } else {
69+ img.data[idx] = Math.round(24 + 205 * a);
70+ img.data[idx + 1] = Math.round(30 + 85 * a);
71+ img.data[idx + 2] = Math.round(38 + 77 * a);
72+ }
73+ img.data[idx + 3] = 255;
74+ }
75+ }
76+ const oc = new OffscreenCanvas(size, size);
77+ oc.getContext('2d')!.putImageData(img, 0, 0);
78+ ctx.drawImage(oc, cx - rad, cy - rad, 2 * rad, 2 * rad);
79+
80+ ctx.strokeStyle = 'rgba(139, 147, 161, 0.5)';
81+ ctx.lineWidth = 1;
82+ ctx.beginPath();
83+ ctx.arc(cx, cy, rad, 0, 2 * Math.PI);
84+ ctx.stroke();
85+
86+ ctx.font = '10px system-ui';
87+ ctx.textAlign = 'center';
88+ ctx.fillStyle = PLOT_COLORS.axis;
89+ ctx.fillText(`${modeF[m].toFixed(0)} Hz · n = ${nDia}`, cx, cy + rad + 12);
90+ }
91+ },
92+ [nodeR, shapeW, modeF, modeN, nShown, rows],
93+ rows * CELL,
94+ );
95+ return <canvas ref={ref} />;
96+}
webapp/src/components/plots/ModeMap.tsxadded+102−0View file
@@ -0,0 +1,102 @@
1+// The retained plate modes: frequency against ring time T60 = 6.91 / (zeta om),
2+// colored by nodal-diameter order and sized by how much of each mode's damping
3+// is radiation. Two things a cymbal shows here that the sibling triangle does
4+// not: the modes are dense enough that the cloud is continuous above about a
5+// kilohertz, and the radiation share is negligible almost everywhere, because
6+// a thin bronze plate below its coincidence frequency is a very poor radiator
7+// and loses its energy internally instead.
8+
9+import { useCanvas, drawFrame, PLOT_COLORS, modeColor } from './canvasUtils';
10+
11+interface Props {
12+ modeF: Float32Array;
13+ modeZeta: Float32Array;
14+ modeZetaRad: Float32Array;
15+ modeN: Float32Array;
16+}
17+
18+const F_LO = 40;
19+const F_HI = 20000;
20+
21+export default function ModeMap({ modeF, modeZeta, modeZetaRad, modeN }: Props) {
22+ let maxN = 1;
23+ for (let i = 0; i < modeN.length; i++) maxN = Math.max(maxN, modeN[i]);
24+ const ref = useCanvas(
25+ (ctx, w, h) => {
26+ drawFrame(ctx, w, h);
27+ const padL = 40;
28+ const padB = 20;
29+ const pw = w - padL - 10;
30+ const ph = h - padB - 14;
31+ const top = 12;
32+ const n = modeF.length;
33+ if (n === 0) return;
34+
35+ const t60 = new Float64Array(n);
36+ let tMin = Infinity;
37+ let tMax = 0;
38+ for (let i = 0; i < n; i++) {
39+ t60[i] = 6.91 / (Math.max(1e-9, modeZeta[i]) * 2 * Math.PI * modeF[i]);
40+ tMin = Math.min(tMin, t60[i]);
41+ tMax = Math.max(tMax, t60[i]);
42+ }
43+ tMin = Math.max(1e-3, tMin / 1.5);
44+ tMax = tMax * 1.5;
45+
46+ const fx = (f: number) =>
47+ padL + ((Math.log10(f) - Math.log10(F_LO)) / (Math.log10(F_HI) - Math.log10(F_LO))) * pw;
48+ const ty = (t: number) =>
49+ top + ph - ((Math.log10(t) - Math.log10(tMin)) / (Math.log10(tMax) - Math.log10(tMin))) * ph;
50+
51+ ctx.font = '10px system-ui';
52+ ctx.textAlign = 'center';
53+ for (const f of [50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000]) {
54+ const x = fx(f);
55+ ctx.strokeStyle = PLOT_COLORS.grid;
56+ ctx.beginPath();
57+ ctx.moveTo(x, top);
58+ ctx.lineTo(x, top + ph);
59+ ctx.stroke();
60+ ctx.fillStyle = PLOT_COLORS.axis;
61+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, x, h - 7);
62+ }
63+ ctx.textAlign = 'right';
64+ for (let e = Math.ceil(Math.log10(tMin)); e <= Math.floor(Math.log10(tMax)); e++) {
65+ const t = 10 ** e;
66+ const y = ty(t);
67+ ctx.strokeStyle = PLOT_COLORS.grid;
68+ ctx.beginPath();
69+ ctx.moveTo(padL, y);
70+ ctx.lineTo(padL + pw, y);
71+ ctx.stroke();
72+ ctx.fillStyle = PLOT_COLORS.axis;
73+ ctx.fillText(t >= 1 ? `${t}s` : `${t * 1000}ms`, padL - 4, y + 3);
74+ }
75+
76+ for (let i = 0; i < n; i++) {
77+ const x = fx(Math.min(F_HI, Math.max(F_LO, modeF[i])));
78+ const y = ty(t60[i]);
79+ const radShare = Math.min(1, modeZetaRad[i] / Math.max(1e-12, modeZeta[i]));
80+ const r = 1.6 + 5 * Math.sqrt(radShare);
81+ ctx.fillStyle = modeColor(modeN[i], maxN, 0.8);
82+ ctx.beginPath();
83+ ctx.arc(x, y, r, 0, 2 * Math.PI);
84+ ctx.fill();
85+ }
86+
87+ ctx.textAlign = 'left';
88+ ctx.fillStyle = PLOT_COLORS.axis;
89+ ctx.fillText(`${n} modes · marker = radiation share of damping · nodal diameters 0`, padL + 6, top + 2);
90+ const x0 = padL + 300;
91+ for (let k = 0; k <= maxN; k++) {
92+ ctx.fillStyle = modeColor(k, maxN, 0.95);
93+ ctx.fillRect(x0 + (k / maxN) * 60, top - 5, 60 / (maxN + 1) + 1, 8);
94+ }
95+ ctx.fillStyle = PLOT_COLORS.axis;
96+ ctx.fillText(`${maxN}`, x0 + 64, top + 2);
97+ },
98+ [modeF, modeZeta, modeZetaRad, modeN, maxN],
99+ 210,
100+ );
101+ return <canvas ref={ref} />;
102+}
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+114−0View file
@@ -0,0 +1,114 @@
1+// Magnitude spectrum of the rendered stroke, log-frequency axis, with the
2+// retained plate modes overlaid and colored by nodal-diameter order. A
3+// cymbal's spectrum is far denser than the sibling triangle's: below a
4+// kilohertz the partials are still resolvable, above it the modal overlap
5+// exceeds one and the spectrum fills in, which is what a wash is.
6+
7+import { magnitudeSpectrum } from '../../lib/fft';
8+import { useCanvas, drawFrame, PLOT_COLORS, modeColor } from './canvasUtils';
9+
10+interface Props {
11+ y: Float32Array;
12+ sampleRate: number;
13+ modeF: Float32Array;
14+ modeN: Float32Array;
15+}
16+
17+const F_LO = 60;
18+const F_HI = 24000;
19+const DB_RANGE = 80;
20+
21+export default function SpectrumPlot({ y, sampleRate, modeF, modeN }: Props) {
22+ let maxN = 1;
23+ for (let i = 0; i < modeN.length; i++) maxN = Math.max(maxN, modeN[i]);
24+ const ref = useCanvas(
25+ (ctx, w, h) => {
26+ drawFrame(ctx, w, h);
27+ const padL = 34;
28+ const padB = 18;
29+ const pw = w - padL - 6;
30+ const ph = h - padB - 8;
31+ const fx = (f: number) =>
32+ padL + ((Math.log10(f) - Math.log10(F_LO)) / (Math.log10(F_HI) - Math.log10(F_LO))) * pw;
33+
34+ ctx.font = '10px system-ui';
35+ ctx.textAlign = 'center';
36+ for (const f of [100, 200, 500, 1000, 2000, 5000, 10000, 20000]) {
37+ const x = fx(f);
38+ ctx.strokeStyle = PLOT_COLORS.grid;
39+ ctx.beginPath();
40+ ctx.moveTo(x, 8);
41+ ctx.lineTo(x, 8 + ph);
42+ ctx.stroke();
43+ ctx.fillStyle = PLOT_COLORS.axis;
44+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, x, h - 5);
45+ }
46+ ctx.textAlign = 'right';
47+ for (let db = 0; db >= -DB_RANGE; db -= 20) {
48+ const yy = 8 + (-db / DB_RANGE) * ph;
49+ ctx.strokeStyle = PLOT_COLORS.grid;
50+ ctx.beginPath();
51+ ctx.moveTo(padL, yy);
52+ ctx.lineTo(padL + pw, yy);
53+ ctx.stroke();
54+ ctx.fillStyle = PLOT_COLORS.axis;
55+ ctx.fillText(`${db}`, padL - 4, yy + 3);
56+ }
57+
58+ const { mag, size } = magnitudeSpectrum(y, 32768);
59+ let peak = 0;
60+ for (let i = 1; i < mag.length; i++) peak = Math.max(peak, mag[i]);
61+ if (peak <= 0) return;
62+
63+ // mode lines first so the trace draws on top
64+ for (let i = 0; i < modeF.length; i++) {
65+ const f = modeF[i];
66+ if (f < F_LO || f > F_HI) continue;
67+ ctx.strokeStyle = modeColor(modeN[i], maxN, 0.35);
68+ ctx.beginPath();
69+ ctx.moveTo(fx(f), 8 + ph * 0.45);
70+ ctx.lineTo(fx(f), 8 + ph);
71+ ctx.stroke();
72+ }
73+
74+ ctx.strokeStyle = PLOT_COLORS.trace;
75+ ctx.beginPath();
76+ let started = false;
77+ for (let px = 0; px <= pw; px++) {
78+ const f = 10 ** (Math.log10(F_LO) + (px / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
79+ const bin = Math.min(mag.length - 1, Math.round((f / sampleRate) * size));
80+ const fNext =
81+ 10 ** (Math.log10(F_LO) + ((px + 1) / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
82+ const binNext = Math.min(
83+ mag.length - 1,
84+ Math.max(bin + 1, Math.round((fNext / sampleRate) * size)),
85+ );
86+ let m = 0;
87+ for (let b = bin; b < binNext; b++) m = Math.max(m, mag[b]);
88+ const db = 20 * Math.log10(m / peak + 1e-12);
89+ const yy = 8 + Math.min(ph, (-db / DB_RANGE) * ph);
90+ if (!started) {
91+ ctx.moveTo(padL + px, yy);
92+ started = true;
93+ } else {
94+ ctx.lineTo(padL + px, yy);
95+ }
96+ }
97+ ctx.stroke();
98+
99+ ctx.textAlign = 'left';
100+ ctx.fillStyle = PLOT_COLORS.axis;
101+ ctx.fillText('mode lines colored by nodal diameters, 0', padL + 4, 18);
102+ const x0 = padL + 200;
103+ for (let n = 0; n <= maxN; n++) {
104+ ctx.fillStyle = modeColor(n, maxN, 0.95);
105+ ctx.fillRect(x0 + (n / maxN) * 70, 11, 70 / (maxN + 1) + 1, 8);
106+ }
107+ ctx.fillStyle = PLOT_COLORS.axis;
108+ ctx.fillText(`${maxN}`, x0 + 74, 18);
109+ },
110+ [y, sampleRate, modeF, modeN, maxN],
111+ 200,
112+ );
113+ return <canvas ref={ref} />;
114+}
webapp/src/components/plots/StrikePlot.tsxadded+77−0View file
@@ -0,0 +1,77 @@
1+// The stroke events: each contact as a stem at its time and velocity, colored
2+// by which radius it lands on. For a single crash this is one stem; for a
3+// repeated figure the jitter and velocity spread are what keep it from
4+// sounding like 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} strokes · cyan = contact radius, amber = second radius`
70+ : `single stroke 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+77−0View file
@@ -0,0 +1,77 @@
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+}
67+
68+// Color a mode by its number of nodal diameters: cyan at n = 0 through amber
69+// at the highest order retained. Used by every plot that shows modes, so the
70+// orders read the same everywhere.
71+export function modeColor(n: number, maxN: number, alpha = 1): string {
72+ const t = Math.min(1, Math.max(0, n / Math.max(1, maxN)));
73+ const r = Math.round(77 + t * (232 - 77));
74+ const g = Math.round(208 + t * (163 - 208));
75+ const b = Math.round(225 + t * (61 - 225));
76+ return `rgba(${r}, ${g}, ${b}, ${alpha})`;
77+}
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 { CymbalParams } 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: CymbalParams;
26+ mfiles: MFile[];
27+ quality: 'preview' | 'full';
28+}
29+
30+export function useRenderer(params: CymbalParams, 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/profile.tsadded+74−0View file
@@ -0,0 +1,74 @@
1+// Mirror of ../matlab/cymbal_profile.m: the cymbal meridian z(r), used by the
2+// 3D scene and the mode gallery. The physics reads the same geometry from the
3+// MATLAB side; this is only for drawing.
4+//
5+// Bow (r >= bell radius): a parabola in the rim-referred coordinate, flat at
6+// the rim. Bell (r < bell radius): the quartic that matches the bow's height,
7+// slope and curvature at the join and comes to a horizontal apex at r = 0.
8+
9+export interface Meridian {
10+ rHole: number;
11+ rBell: number;
12+ rRim: number;
13+ z(r: number): number;
14+ thickness(r: number): number;
15+ apexHeight: number;
16+}
17+
18+export function cymbalMeridian(
19+ radius: number,
20+ holeRadius: number,
21+ bellRadius: number,
22+ bellHeight: number,
23+ bowHeight: number,
24+ thicknessBell: number,
25+ thicknessRim: number,
26+): Meridian {
27+ const a = radius;
28+ const rh = Math.max(1e-4, Math.min(holeRadius, 0.5 * a));
29+ const rb = Math.min(Math.max(bellRadius, 2 * rh), 0.8 * a);
30+ const den = a - rb;
31+ const mb = (-2 * bowHeight) / den;
32+ const kb = (2 * bowHeight) / (den * den);
33+ // solve the 3x3 for the quartic coefficients (B r^2 + C r^3 + E r^4)
34+ const M = [
35+ [rb * rb, rb ** 3, rb ** 4],
36+ [2 * rb, 3 * rb * rb, 4 * rb ** 3],
37+ [2, 6 * rb, 12 * rb * rb],
38+ ];
39+ const rhs = [-bellHeight, mb, kb];
40+ const [B, C, E] = solve3(M, rhs);
41+ const A0 = bowHeight + bellHeight;
42+ return {
43+ rHole: rh,
44+ rBell: rb,
45+ rRim: a,
46+ apexHeight: A0,
47+ z(r: number): number {
48+ if (r < rb) return A0 + B * r * r + C * r ** 3 + E * r ** 4;
49+ const u = (a - r) / den;
50+ return bowHeight * u * u;
51+ },
52+ thickness(r: number): number {
53+ const taper = Math.min(1, Math.max(0, (a - r) / den));
54+ return thicknessRim + (thicknessBell - thicknessRim) * taper;
55+ },
56+ };
57+}
58+
59+function solve3(M: number[][], b: number[]): [number, number, number] {
60+ const A = M.map((row, i) => [...row, b[i]]);
61+ for (let c = 0; c < 3; c++) {
62+ let piv = c;
63+ for (let r = c + 1; r < 3; r++) if (Math.abs(A[r][c]) > Math.abs(A[piv][c])) piv = r;
64+ [A[c], A[piv]] = [A[piv], A[c]];
65+ const d = A[c][c] || 1e-30;
66+ for (let k = c; k < 4; k++) A[c][k] /= d;
67+ for (let r = 0; r < 3; r++) {
68+ if (r === c) continue;
69+ const f = A[r][c];
70+ for (let k = c; k < 4; k++) A[r][k] -= f * A[c][k];
71+ }
72+ }
73+ return [A[0][3], A[1][3], A[2][3]];
74+}
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+47−0View file
@@ -0,0 +1,47 @@
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+ 'cymbal_profile.m',
25+ 'cymbal_shell.m',
26+ 'cymbal_basis.m',
27+ 'cymbal_at.m',
28+ 'cymbal_point.m',
29+ 'cymbal_strikes.m',
30+ 'cymbal_rng.m',
31+ 'cymbal_cascade.m',
32+ 'cymbal_radiate.m',
33+ 'cymbal_output.m',
34+]);
35+
36+// Cheap content hash (FNV-1a) so the worker can key its caches on the code.
37+export function hashMFiles(files: MFile[]): string {
38+ let h = 0x811c9dc5;
39+ for (const f of files) {
40+ const s = f.name + ' ' + f.source + '';
41+ for (let i = 0; i < s.length; i++) {
42+ h ^= s.charCodeAt(i);
43+ h = Math.imul(h, 0x01000193);
44+ }
45+ }
46+ return (h >>> 0).toString(36);
47+}
webapp/src/params.tsadded+310−0View file
@@ -0,0 +1,310 @@
1+// Parameter definitions for the cymbal model, mirroring cymbal_defaults.m.
2+// Each parameter is assigned a tier that determines how much recomputation a
3+// change requires:
4+// solve -> shell eigensolves + strikes + cascade + radiation + output
5+// radiate -> strikes + cascade + radiation + output only (the cached mode
6+// basis is reused; the eigensolves are the expensive part here)
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 CymbalParams {
13+ // Plate geometry
14+ radius_m: number;
15+ hole_radius_m: number;
16+ bell_radius_m: number;
17+ bell_height_m: number;
18+ bow_height_m: number;
19+ thickness_bell_m: number;
20+ thickness_rim_m: number;
21+ // Material
22+ density: number;
23+ youngs_modulus: number;
24+ poisson_ratio: number;
25+ stiffness_trim: number;
26+ // Discretization / mode retention
27+ n_elements: number;
28+ n_diameters: number;
29+ min_mode_hz: number;
30+ max_mode_hz: number;
31+ // Damping
32+ damp_rate: number;
33+ loss_factor: number;
34+ loss_slope: number;
35+ radiation_loss: number;
36+ mount_loss: number;
37+ // Stick
38+ tip_radius_m: number;
39+ stick_mass_kg: number;
40+ stick_youngs: number;
41+ restitution: number;
42+ // Strike gesture
43+ strike_velocity: number;
44+ strike_r: number;
45+ strike_angle_deg: number;
46+ strike_count: number;
47+ strike_rate_hz: number;
48+ roll_r2: number;
49+ timing_jitter_ms: number;
50+ velocity_spread: number;
51+ seed: number;
52+ // Nonlinear cascade
53+ cascade_gain: number;
54+ cascade_back: number;
55+ cascade_ratio: number;
56+ cascade_leak: number;
57+ cascade_step_ms: number;
58+ // Radiation / listener
59+ direct_gain: number;
60+ dipole_scale: number;
61+ n_azimuth: number;
62+ listener_x_m: number;
63+ listener_y_m: number;
64+ listener_z_m: number;
65+ speed_of_sound: number;
66+ air_density: number;
67+ // Rendering
68+ sample_rate: number;
69+ oversample: number;
70+ tail_s: number;
71+ pre_delay_ms: number;
72+ normalize: boolean;
73+ peak_dbfs: number;
74+ full_scale_pa: number;
75+ bit_depth: number;
76+}
77+
78+export const DEFAULTS: CymbalParams = {
79+ radius_m: 0.2032,
80+ hole_radius_m: 0.0064,
81+ bell_radius_m: 0.045,
82+ bell_height_m: 0.017,
83+ bow_height_m: 0.02,
84+ thickness_bell_m: 0.00115,
85+ thickness_rim_m: 0.00075,
86+ density: 8700.0,
87+ youngs_modulus: 1.05e11,
88+ poisson_ratio: 0.33,
89+ stiffness_trim: 1.0,
90+ n_elements: 28,
91+ n_diameters: 48,
92+ min_mode_hz: 40.0,
93+ max_mode_hz: 8000.0,
94+ damp_rate: 0.35,
95+ loss_factor: 1.5e-4,
96+ loss_slope: 2.0e-9,
97+ radiation_loss: 1.0,
98+ mount_loss: 3.0,
99+ tip_radius_m: 0.005,
100+ stick_mass_kg: 0.03,
101+ stick_youngs: 5.0e9,
102+ restitution: 0.55,
103+ strike_velocity: 3.0,
104+ strike_r: 0.94,
105+ strike_angle_deg: 40.0,
106+ strike_count: 1,
107+ strike_rate_hz: 6.0,
108+ roll_r2: 0.8,
109+ timing_jitter_ms: 8.0,
110+ velocity_spread: 0.18,
111+ seed: 1,
112+ cascade_gain: 800.0,
113+ cascade_back: 0.45,
114+ cascade_ratio: 1.19,
115+ cascade_leak: 1.0,
116+ cascade_step_ms: 0.5,
117+ direct_gain: 1.0,
118+ dipole_scale: 1.0,
119+ n_azimuth: 96,
120+ listener_x_m: 0.3,
121+ listener_y_m: 0.0,
122+ listener_z_m: 0.5,
123+ speed_of_sound: 343.0,
124+ air_density: 1.204,
125+ sample_rate: 48000,
126+ oversample: 2,
127+ tail_s: 3.0,
128+ pre_delay_ms: 3.0,
129+ normalize: true,
130+ peak_dbfs: -3.0,
131+ full_scale_pa: 5.0,
132+ bit_depth: 24,
133+};
134+
135+export type Tier = 'solve' | 'radiate' | 'output';
136+
137+// Keys whose change invalidates the cached mode basis (the shell eigensolves).
138+// The stick, the gesture, the cascade and the listener are deliberately
139+// absent: striking somewhere else or moving the microphone reuses the modes.
140+export const SOLVE_KEYS: (keyof CymbalParams)[] = [
141+ 'radius_m', 'hole_radius_m', 'bell_radius_m', 'bell_height_m', 'bow_height_m',
142+ 'thickness_bell_m', 'thickness_rim_m',
143+ 'density', 'youngs_modulus', 'poisson_ratio', 'stiffness_trim',
144+ 'n_elements', 'n_diameters', 'min_mode_hz', 'max_mode_hz',
145+ 'damp_rate', 'loss_factor', 'loss_slope', 'radiation_loss', 'mount_loss',
146+ 'dipole_scale', 'speed_of_sound', 'air_density',
147+];
148+
149+export const RADIATE_KEYS: (keyof CymbalParams)[] = [
150+ 'tip_radius_m', 'stick_mass_kg', 'stick_youngs', 'restitution',
151+ 'strike_velocity', 'strike_r', 'strike_angle_deg',
152+ 'strike_count', 'strike_rate_hz', 'roll_r2',
153+ 'timing_jitter_ms', 'velocity_spread', 'seed',
154+ 'cascade_gain', 'cascade_back', 'cascade_ratio', 'cascade_leak', 'cascade_step_ms',
155+ 'direct_gain', 'n_azimuth',
156+ 'listener_x_m', 'listener_y_m', 'listener_z_m',
157+ 'sample_rate', 'tail_s',
158+];
159+
160+export const OUTPUT_KEYS: (keyof CymbalParams)[] = [
161+ 'pre_delay_ms', 'normalize', 'peak_dbfs', 'full_scale_pa',
162+];
163+
164+export interface ParamSpec {
165+ key: keyof CymbalParams;
166+ label: string;
167+ unit?: string;
168+ min: number;
169+ max: number;
170+ step?: number;
171+ log?: boolean;
172+ integer?: boolean;
173+ help?: string;
174+}
175+
176+export interface ParamGroup {
177+ name: string;
178+ specs: ParamSpec[];
179+}
180+
181+export const PARAM_GROUPS: ParamGroup[] = [
182+ {
183+ name: 'Stroke',
184+ specs: [
185+ { key: 'strike_velocity', label: 'Stroke velocity', unit: 'm/s', min: 0.1, max: 8.0, step: 0.05, help: 'Stick speed at contact. Dynamics and timbre together: the cascade strength goes as the square of the amplitude, so a soft stroke pings and a hard one crashes.' },
186+ { key: 'strike_r', label: 'Contact radius', min: 0.02, max: 0.995, step: 0.005, help: 'Arclength fraction from the hole (0) to the rim (1). ~0.95 is a crash on the edge, ~0.6 a ride on the bow, ~0.1 the bell.' },
187+ { key: 'strike_angle_deg', label: 'Stroke angle', unit: '°', min: 0, max: 85, step: 1, help: '0 strikes square onto the surface, 85 is a glancing swipe along it. Only the normal component drives the plate strongly.' },
188+ ],
189+ },
190+ {
191+ name: 'Nonlinear cascade',
192+ specs: [
193+ { key: 'cascade_gain', label: 'Cascade gain', min: 0, max: 4000, step: 10, help: 'Coupling constant of the energy cascade. 0 leaves an exactly linear plate, which sounds like a gong; raising it opens the crash up. This is the one fitted number in the model.' },
194+ { key: 'cascade_ratio', label: 'Step ratio', min: 1.03, max: 2.0, step: 0.01, help: 'Frequency ratio of one transfer step. Smaller means a slower, smoother climb; 1.19 is a quarter of an octave.' },
195+ { key: 'cascade_back', label: 'Backward share', min: 0, max: 0.95, step: 0.01, help: 'Downward transfers per upward one. At 1 the flux stops and energy only diffuses.' },
196+ { key: 'cascade_leak', label: 'Top-of-band leak', min: 0, max: 1, step: 0.05, help: 'Fraction of the flux leaving the highest retained mode that is lost, standing in for the truncated band above. 0 reflects it back.' },
197+ ],
198+ },
199+ {
200+ name: 'Repeats',
201+ specs: [
202+ { key: 'strike_count', label: 'Strokes', min: 1, max: 32, step: 1, integer: true, help: '1 is a single crash; more gives a repeated figure alternating between the two contact radii.' },
203+ { key: 'strike_rate_hz', label: 'Rate', unit: 'Hz', min: 1, max: 20, step: 0.5, help: 'Strokes per second.' },
204+ { key: 'roll_r2', label: 'Second radius', min: 0.02, max: 0.995, step: 0.005, help: 'Even-numbered strokes land here: a ride pattern moving between the bow and the edge.' },
205+ { key: 'timing_jitter_ms', label: 'Timing jitter', unit: 'ms', min: 0, max: 40, step: 0.5, help: 'Gaussian jitter on each stroke after the first. 0 sounds mechanical.' },
206+ { key: 'velocity_spread', label: 'Velocity spread', min: 0, max: 0.7, step: 0.01, help: 'Lognormal spread of stroke velocities.' },
207+ { key: 'seed', label: 'Random seed', min: 1, max: 9999, step: 1, integer: true, help: 'Any change redraws the timing and velocities; the same seed always reproduces them.' },
208+ ],
209+ },
210+ {
211+ name: 'Stick',
212+ specs: [
213+ { key: 'tip_radius_m', label: 'Tip radius', unit: 'm', min: 0.001, max: 0.02, step: 0.0005, help: 'Radius of the stick tip; enters the Hertz contact stiffness.' },
214+ { key: 'stick_mass_kg', label: 'Stick mass', unit: 'kg', min: 0.002, max: 0.3, log: true, help: 'A heavy stick or mallet lengthens the contact and darkens the stroke.' },
215+ { key: 'stick_youngs', label: 'Tip modulus', unit: 'Pa', min: 1e8, max: 2e11, log: true, help: 'Nylon ~5e9, wood ~1.2e10, felt much lower. A soft tip lengthens the contact and rolls the attack off.' },
216+ { key: 'restitution', label: 'Restitution', min: 0.05, max: 0.95, step: 0.01, help: 'Coefficient of restitution; scales the momentum each stroke transfers.' },
217+ ],
218+ },
219+ {
220+ name: 'Plate geometry',
221+ specs: [
222+ { key: 'radius_m', label: 'Radius', unit: 'm', min: 0.09, max: 0.32, step: 0.001, help: '0.2032 is a 16-inch cymbal. Mode frequencies scale roughly as 1/R².' },
223+ { key: 'thickness_rim_m', label: 'Rim thickness', unit: 'm', min: 0.0003, max: 0.003, step: 0.00005, help: 'A thin edge is a crash; a thick one is a ride. Frequencies scale as the thickness.' },
224+ { key: 'thickness_bell_m', label: 'Bell thickness', unit: 'm', min: 0.0004, max: 0.004, step: 0.00005, help: 'The taper from bell to rim is what spreads the partials of a cymbal apart.' },
225+ { key: 'bell_radius_m', label: 'Bell radius', unit: 'm', min: 0.015, max: 0.09, step: 0.001, help: 'A larger bell stiffens the middle and pushes the low modes up.' },
226+ { key: 'bell_height_m', label: 'Bell height', unit: 'm', min: 0.0, max: 0.05, step: 0.001, help: 'Height of the dome above the bow.' },
227+ { key: 'bow_height_m', label: 'Bow rise', unit: 'm', min: 0.0, max: 0.06, step: 0.001, help: 'How far the bow rises from rim to bell. 0 is a flat disc, which rings quite differently: the curvature couples bending to stretching.' },
228+ { key: 'hole_radius_m', label: 'Hole radius', unit: 'm', min: 0.002, max: 0.02, step: 0.0005 },
229+ ],
230+ },
231+ {
232+ name: 'Material',
233+ specs: [
234+ { key: 'density', label: 'Density', unit: 'kg/m³', min: 2000, max: 12000, log: true, help: 'B20 bronze 8700, brass ~8500, nickel silver ~8700, aluminium 2700.' },
235+ { key: 'youngs_modulus', label: "Young's modulus", unit: 'Pa', min: 3e10, max: 3e11, log: true, help: 'Bell bronze ~1.05e11. Frequencies scale as sqrt(E).' },
236+ { key: 'poisson_ratio', label: 'Poisson ratio', min: 0.0, max: 0.45, step: 0.01 },
237+ { key: 'stiffness_trim', label: 'Stiffness trim', min: 0.4, max: 2.0, step: 0.01, help: 'Overall pitch trim for the hammering and lathing a real cymbal carries and this profile does not.' },
238+ ],
239+ },
240+ {
241+ name: 'Damping',
242+ specs: [
243+ { key: 'damp_rate', label: 'Decay rate', unit: 's⁻¹', min: 0.02, max: 5.0, log: true, help: 'Frequency-independent decay rate: T60 = 6.91 / this, before the other terms. It is what makes a cymbal ring for tens of seconds.' },
244+ { key: 'loss_factor', label: 'Loss factor η₀', min: 1e-6, max: 3e-3, log: true, help: 'Constant material loss factor; gives a T60 falling as 1/f, so it shortens the shimmer before the hum.' },
245+ { key: 'loss_slope', label: 'Loss slope η₁', unit: 's', min: 1e-11, max: 1e-7, log: true, help: 'Extra high-frequency damping; controls how quickly the top of the wash fades.' },
246+ { key: 'mount_loss', label: 'Mount damping', unit: 'N s/m', min: 0, max: 200, step: 0.5, help: 'A dashpot on the axial motion at the hole, standing in for the felt washer and stand. It reaches the n = 0 and n = 1 modes and almost nothing else.' },
247+ { key: 'radiation_loss', label: 'Radiation loss', min: 0.0, max: 3.0, step: 0.05, help: 'Scales the radiation damping computed from the radiated power. Metal plates are poor radiators, so this barely matters below a few kHz.' },
248+ ],
249+ },
250+ {
251+ name: 'Listener',
252+ specs: [
253+ { key: 'listener_x_m', label: 'Listener x', unit: 'm', min: -3.0, max: 3.0, step: 0.01, help: 'In the plane of the rim. The stroke is on the +x meridian.' },
254+ { key: 'listener_y_m', label: 'Listener y', unit: 'm', min: -3.0, max: 3.0, step: 0.01 },
255+ { key: 'listener_z_m', label: 'Listener z', unit: 'm', min: -3.0, max: 3.0, step: 0.01, help: 'Along the symmetry axis. Directly above the bell every mode with a nodal diameter cancels exactly, leaving only the n = 0 modes.' },
256+ { 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.' },
257+ { key: 'dipole_scale', label: 'Dipole length scale', min: 0.1, max: 3.0, step: 0.05, help: 'Scales the acoustic short-circuit length. The one phenomenological length in the radiation model.' },
258+ ],
259+ },
260+ {
261+ name: 'Output',
262+ specs: [
263+ { key: 'tail_s', label: 'Tail', unit: 's', min: 0.3, max: 12.0, step: 0.1, help: 'Ring time rendered after the last stroke.' },
264+ { key: 'pre_delay_ms', label: 'Pre-delay', unit: 'ms', min: 0, max: 50, step: 0.5 },
265+ { key: 'peak_dbfs', label: 'Peak level', unit: 'dBFS', min: -24, max: 0, step: 0.5, help: 'Used when normalize is on.' },
266+ { key: 'full_scale_pa', label: 'Full scale', unit: 'Pa', min: 0.1, max: 50, step: 0.1, help: 'Used when normalize is off: this many pascals maps to 0 dBFS.' },
267+ ],
268+ },
269+ {
270+ name: 'Advanced / quality',
271+ specs: [
272+ { key: 'n_elements', label: 'Meridian elements', min: 12, max: 60, step: 1, integer: true, help: 'Shell elements along the meridian. More converges the high modes; each eigensolve costs its cube.' },
273+ { key: 'n_diameters', label: 'Max nodal diameters', min: 4, max: 64, step: 1, integer: true, help: 'Highest circumferential order solved for. One eigensolve per order, so this is the main cost. Too low and the band above a couple of kHz is missing its high-order modes.' },
274+ { key: 'min_mode_hz', label: 'Min mode freq', unit: 'Hz', min: 20, max: 400, step: 5 },
275+ { key: 'max_mode_hz', label: 'Max mode freq', unit: 'Hz', min: 1000, max: 16000, step: 250, help: 'Everything above is truncated; the cascade treats it as a sink and the direct click stands in for it.' },
276+ { key: 'cascade_step_ms', label: 'Cascade step', unit: 'ms', min: 0.1, max: 4.0, step: 0.1, help: 'Control-rate step of the energy solve, and the envelope resolution of the oscillator bank.' },
277+ { key: 'n_azimuth', label: 'Azimuth samples', min: 24, max: 256, step: 8, integer: true, help: 'Azimuthal samples in the listener transfer. Must comfortably exceed twice the highest nodal-diameter order.' },
278+ { key: 'oversample', label: 'Oversample', min: 1, max: 4, step: 1, integer: true, help: 'Synthesis rate multiple; 2 for render, 1 for preview. The mode bank needs none, the contact click does.' },
279+ { key: 'speed_of_sound', label: 'Speed of sound', unit: 'm/s', min: 300, max: 400, step: 1 },
280+ { key: 'air_density', label: 'Air density', unit: 'kg/m³', min: 0.9, max: 1.5, step: 0.01 },
281+ ],
282+ },
283+];
284+
285+const STORAGE_KEY = 'cymbal-webapp-params-v1';
286+
287+export function loadParams(): CymbalParams {
288+ try {
289+ const raw = localStorage.getItem(STORAGE_KEY);
290+ if (!raw) return { ...DEFAULTS };
291+ const parsed = JSON.parse(raw);
292+ const out = { ...DEFAULTS };
293+ for (const k of Object.keys(DEFAULTS) as (keyof CymbalParams)[]) {
294+ if (typeof parsed[k] === typeof DEFAULTS[k]) {
295+ (out as Record<string, unknown>)[k] = parsed[k];
296+ }
297+ }
298+ return out;
299+ } catch {
300+ return { ...DEFAULTS };
301+ }
302+}
303+
304+export function saveParams(p: CymbalParams): void {
305+ try {
306+ localStorage.setItem(STORAGE_KEY, JSON.stringify(p));
307+ } catch {
308+ // storage unavailable; persistence is a convenience only
309+ }
310+}
webapp/src/selftest.tsadded+25−0View file
@@ -0,0 +1,25 @@
1+// Runs the MATLAB validation script (../matlab/cymbal_validate.m) through
2+// numbl in Node: closed-form identities, exact reproducibility of the seeded
3+// strike 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 === 'cymbal_validate.m')!;
17+const t0 = performance.now();
18+const result = executeCode(main.source, { displayResults: false }, workspaceFiles, 'cymbal_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+481−0View file
@@ -0,0 +1,481 @@
1+// Web Worker that renders the cymbal stroke 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 -> cymbal_basis (one shell eigensolve per nodal-diameter
6+// order; the expensive stage, cached on geometry, material,
7+// damping and discretization parameters)
8+// radiate tier -> cymbal_strikes + cymbal_cascade + cymbal_radiate (cached
9+// on the stick, the gesture, the cascade, the listener and
10+// the synthesis rate)
11+// output tier -> cymbal_output (always cheap, always re-run)
12+// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).
13+
14+import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
15+import { hashMFiles, type MFile } from '../mfiles';
16+import { SOLVE_KEYS, RADIATE_KEYS, type CymbalParams } from '../params';
17+
18+export interface RenderRequest {
19+ type: 'render';
20+ id: number;
21+ params: CymbalParams;
22+ quality: 'preview' | 'full';
23+ mfiles: MFile[];
24+}
25+
26+export interface RenderResponse {
27+ type: 'result';
28+ id: number;
29+ quality: 'preview' | 'full';
30+ y: Float32Array; // final signal, [-1, 1]
31+ sampleRate: number;
32+ strikeT: Float32Array; // stroke times, s
33+ strikeV: Float32Array; // stroke velocities, m/s
34+ strikeS: Float32Array; // stroke arclength fractions
35+ modeF: Float32Array; // mode frequencies, Hz
36+ modeZeta: Float32Array; // damping ratio of each mode
37+ modeZetaRad: Float32Array; // radiation part of the damping
38+ modeSigma: Float32Array; // radiation efficiency of each mode
39+ modeN: Float32Array; // nodal diameters
40+ nodeR: Float32Array; // meridian nodes, m
41+ nodeZ: Float32Array;
42+ shapeW: Float32Array; // normal displacement per mode, nn x n col-major
43+ cascadeDb: Float32Array; // energy grid, nBands x nFrames col-major, dB
44+ cascadeBands: number;
45+ cascadeFrames: number;
46+ cascadeFLo: number;
47+ cascadeFHi: number;
48+ cascadeTEnd: number;
49+ diag: {
50+ modesKept: number;
51+ maxN: number;
52+ strikeCount: number;
53+ meanV: number;
54+ maxV: number;
55+ meanTauUs: number;
56+ peakForceN: number;
57+ durationS: number;
58+ peakPa: number;
59+ peakSplDb: number;
60+ directPeakPa: number;
61+ listenerDistM: number;
62+ peakWrmsUm: number;
63+ peakStrength: number;
64+ cascadeLostFrac: number;
65+ plateMassKg: number;
66+ warnings: string[];
67+ solveMs: number; // 0 when the cached eigensolve was reused
68+ radiateMs: number; // 0 when the cached radiation was reused
69+ solveCached: boolean;
70+ radiateCached: boolean;
71+ oversample: number;
72+ };
73+}
74+
75+export interface ProgressMessage {
76+ type: 'progress';
77+ id: number;
78+ stage: 'solve' | 'radiate' | 'output';
79+ frac: number;
80+}
81+
82+// Raised by a fault in the .m code (the user can edit it), so it is reported
83+// rather than left to stall the request.
84+export interface ErrorMessage {
85+ type: 'error';
86+ id: number;
87+ message: string;
88+ file?: string;
89+ line?: number;
90+}
91+
92+interface SolveEntry {
93+ vars: Record<string, RuntimeValue>;
94+ modeF: Float32Array;
95+ modeZeta: Float32Array;
96+ modeZetaRad: Float32Array;
97+ modeSigma: Float32Array;
98+ modeN: Float32Array;
99+ nodeR: Float32Array;
100+ nodeZ: Float32Array;
101+ shapeW: Float32Array;
102+ maxN: number;
103+ plateMassKg: number;
104+}
105+
106+interface RadiateEntry {
107+ p: RuntimeValue;
108+ strikeT: Float32Array;
109+ strikeV: Float32Array;
110+ strikeS: Float32Array;
111+ cascadeDb: Float32Array;
112+ cascadeBands: number;
113+ cascadeFrames: number;
114+ cascadeTEnd: number;
115+ strikeCount: number;
116+ meanV: number;
117+ maxV: number;
118+ meanTauUs: number;
119+ peakForceN: number;
120+ durationS: number;
121+ directPeakPa: number;
122+ listenerDistM: number;
123+ listenerInside: boolean;
124+ peakWrmsUm: number;
125+ peakStrength: number;
126+ cascadeLostFrac: number;
127+ oversample: number;
128+}
129+
130+// Small LRU caches. The solve cache does not depend on the oversample rate,
131+// so a preview render right after a full one reuses the same eigensolves.
132+const solveCache = new Map<string, SolveEntry>();
133+const radiateCache = new Map<string, RadiateEntry>();
134+const SOLVE_CACHE_MAX = 2;
135+const RADIATE_CACHE_MAX = 3;
136+
137+function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
138+ const v = cache.get(key);
139+ if (v !== undefined) {
140+ cache.delete(key);
141+ cache.set(key, v);
142+ }
143+ return v;
144+}
145+
146+function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
147+ cache.delete(key);
148+ cache.set(key, value);
149+ while (cache.size > max) {
150+ const oldest = cache.keys().next().value as string;
151+ cache.delete(oldest);
152+ }
153+}
154+
155+function keyFor(params: CymbalParams, keys: readonly (keyof CymbalParams)[]): string {
156+ return keys.map((k) => String(params[k])).join('|');
157+}
158+
159+// MATLAB source lines assigning every parameter onto a struct P.
160+function pLiteral(params: CymbalParams): string {
161+ const lines: string[] = [];
162+ for (const [k, v] of Object.entries(params)) {
163+ if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
164+ else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
165+ }
166+ return lines.join('\n');
167+}
168+
169+function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
170+ return executeCode(
171+ source,
172+ { initialVariableValues: inputs, displayResults: false },
173+ mfiles,
174+ 'stage.m',
175+ );
176+}
177+
178+function asTensor(v: RuntimeValue): RuntimeTensor {
179+ const t = v as RuntimeTensor;
180+ if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
181+ return t;
182+}
183+
184+function asScalar(v: RuntimeValue): number {
185+ if (typeof v === 'number') return v;
186+ const t = v as RuntimeTensor;
187+ if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
188+ if (typeof v === 'boolean') return v ? 1 : 0;
189+ throw new Error('expected scalar result from numbl stage');
190+}
191+
192+function toF32(v: RuntimeValue): Float32Array {
193+ const t = asTensor(v);
194+ const out = new Float32Array(t.data.length);
195+ for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
196+ return out;
197+}
198+
199+// The basis fields handed from the solve stage into the radiate stage, in the
200+// order the reassembly script expects.
201+const BASIS_FIELDS = [
202+ 'f', 'omega', 'zeta', 'zetaRad', 'zetaMount', 'sigma', 'ndia', 'ldip',
203+ 'NU', 'NW', 'NT',
204+ 'rnode', 'znode', 'snode', 'drdsnode', 'dzdsnode', 'thknode', 'csnode', 'ctnode',
205+ 'surfR', 'surfZ', 'surfDr', 'surfDz', 'surfLen', 'SW',
206+ 'S', 'area', 'mass', 'thkMean',
207+];
208+
209+const CASCADE_BANDS = 48;
210+const CASCADE_FRAMES = 240;
211+
212+const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
213+
214+self.onmessage = (ev: MessageEvent<RenderRequest>) => {
215+ const msg = ev.data;
216+ if (msg.type !== 'render') return;
217+
218+ // numbl reports a file it could not parse through console.warn and then
219+ // skips it, so the failure surfaces at the call site as "Undefined
220+ // function". Capture the warning to report the actual syntax error instead.
221+ const warnings: string[] = [];
222+ const origWarn = console.warn;
223+ console.warn = (...args: unknown[]) => {
224+ warnings.push(args.map(String).join(' '));
225+ };
226+
227+ try {
228+ render(msg);
229+ } catch (err) {
230+ const e = err as { message?: string; file?: string; line?: number };
231+ let m: ErrorMessage = {
232+ type: 'error',
233+ id: msg.id,
234+ message: e?.message ? String(e.message) : String(err),
235+ file: typeof e?.file === 'string' ? e.file : undefined,
236+ line: typeof e?.line === 'number' ? e.line : undefined,
237+ };
238+ for (const w of warnings) {
239+ const hit = SKIP_WARNING.exec(w);
240+ if (hit) {
241+ m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
242+ break;
243+ }
244+ }
245+ (self as unknown as Worker).postMessage(m);
246+ } finally {
247+ console.warn = origWarn;
248+ }
249+};
250+
251+function render(msg: RenderRequest): void {
252+ const { id, params, mfiles } = msg;
253+ let quality = msg.quality;
254+
255+ // The .m sources are part of the cache identity: editing the code must
256+ // invalidate everything computed from the previous version.
257+ const codeKey = hashMFiles(mfiles);
258+ const solveKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);
259+
260+ const fullOversample = params.oversample;
261+ const radiateParamsKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);
262+ if (quality === 'preview' && radiateCache.has(radiateParamsKey + `|os=${fullOversample}`)) {
263+ quality = 'full'; // full-rate radiation already cached
264+ }
265+ const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample;
266+ const fsHi = params.sample_rate * oversample;
267+
268+ const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
269+
270+ // ---- solve stage: one shell eigensolve per nodal-diameter order ----
271+ let solveMs = 0;
272+ let solveCached = true;
273+ let solveEntry = lruGet(solveCache, solveKey);
274+ if (!solveEntry) {
275+ solveCached = false;
276+ post({ type: 'progress', id, stage: 'solve', frac: 0.12 });
277+ const t0 = performance.now();
278+ const r = runStage(mfiles, `
279+${pLiteral(params)}
280+basis = cymbal_basis(P);
281+${BASIS_FIELDS.map((f) => `b_${f} = basis.${f};`).join('\n')}
282+maxN = max(basis.ndia);
283+`);
284+ solveMs = performance.now() - t0;
285+
286+ const vars: Record<string, RuntimeValue> = {};
287+ for (const f of BASIS_FIELDS) vars[`b_${f}`] = r.variableValues[`b_${f}`];
288+ solveEntry = {
289+ vars,
290+ modeF: toF32(r.variableValues.b_f),
291+ modeZeta: toF32(r.variableValues.b_zeta),
292+ modeZetaRad: toF32(r.variableValues.b_zetaRad),
293+ modeSigma: toF32(r.variableValues.b_sigma),
294+ modeN: toF32(r.variableValues.b_ndia),
295+ nodeR: toF32(r.variableValues.b_rnode),
296+ nodeZ: toF32(r.variableValues.b_znode),
297+ shapeW: toF32(r.variableValues.b_NW),
298+ maxN: asScalar(r.variableValues.maxN),
299+ plateMassKg: asScalar(r.variableValues.b_mass),
300+ };
301+ lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
302+ }
303+
304+ // ---- radiate stage: strikes + cascade + modal synthesis ----
305+ const radiateKey = radiateParamsKey + `|os=${oversample}`;
306+ let radiateMs = 0;
307+ let radiateCached = true;
308+ let radEntry = lruGet(radiateCache, radiateKey);
309+ if (!radEntry) {
310+ radiateCached = false;
311+ post({ type: 'progress', id, stage: 'radiate', frac: 0.5 });
312+ const t0 = performance.now();
313+ const r = runStage(mfiles, `
314+${pLiteral(params)}
315+${BASIS_FIELDS.map((f) => `basis.${f} = b_${f};`).join('\n')}
316+basis.n = length(basis.f);
317+[tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = cymbal_strikes(P, basis);
318+[env, tframe, DV, kf, cinfo] = cymbal_cascade(P, basis, tstr, jstr, taustr, sstr, sinfo.durationS);
319+[p, radinfo] = cymbal_radiate(P, basis, env, tframe, DV, kf, tstr, taustr, fstr, sstr, ${fsHi});
320+% collapse the modal energies onto a log-frequency by time grid for the plot
321+nb = ${CASCADE_BANDS};
322+nt = ${CASCADE_FRAMES};
323+fLo = max(1, min(basis.f));
324+fHi = max(fLo * 1.2, max(basis.f));
325+bandOf = min(nb, max(1, 1 + floor(nb * (log(basis.f) - log(fLo)) / (log(fHi) - log(fLo) + 1e-12))));
326+nf = size(env, 2);
327+grid = zeros(nb, nt);
328+for q = 1:nt
329+ c0 = 1 + floor((q - 1) * nf / nt);
330+ c1 = max(c0, floor(q * nf / nt));
331+ e = sum(env(:, c0:c1).^2, 2) / (c1 - c0 + 1);
332+ grid(:, q) = accumarray(bandOf, e, [nb, 1]);
333+end
334+cascadeDb = 10 * log10(grid + 1e-30);
335+cascadeTEnd = tframe(length(tframe));
336+strikeCount = sinfo.count;
337+meanV = sinfo.meanV;
338+maxV = sinfo.maxV;
339+meanTauUs = sinfo.meanTauUs;
340+peakForceN = sinfo.peakForceN;
341+durationS = sinfo.durationS;
342+directPeakPa = radinfo.directPeakPa;
343+listenerDistM = radinfo.listenerDistM;
344+listenerInside = radinfo.listenerInside;
345+peakWrmsUm = cinfo.peakWrmsUm;
346+peakStrength = cinfo.peakStrength;
347+cascadeLostFrac = cinfo.cascadeLostJ / max(1e-30, cinfo.injectedJ);
348+`, solveEntry.vars);
349+ radiateMs = performance.now() - t0;
350+ radEntry = {
351+ p: r.variableValues.p,
352+ strikeT: toF32(r.variableValues.tstr),
353+ strikeV: toF32(r.variableValues.vstr),
354+ strikeS: toF32(r.variableValues.sstr),
355+ cascadeDb: toF32(r.variableValues.cascadeDb),
356+ cascadeBands: CASCADE_BANDS,
357+ cascadeFrames: CASCADE_FRAMES,
358+ cascadeTEnd: asScalar(r.variableValues.cascadeTEnd),
359+ strikeCount: asScalar(r.variableValues.strikeCount),
360+ meanV: asScalar(r.variableValues.meanV),
361+ maxV: asScalar(r.variableValues.maxV),
362+ meanTauUs: asScalar(r.variableValues.meanTauUs),
363+ peakForceN: asScalar(r.variableValues.peakForceN),
364+ durationS: asScalar(r.variableValues.durationS),
365+ directPeakPa: asScalar(r.variableValues.directPeakPa),
366+ listenerDistM: asScalar(r.variableValues.listenerDistM),
367+ listenerInside: asScalar(r.variableValues.listenerInside) !== 0,
368+ peakWrmsUm: asScalar(r.variableValues.peakWrmsUm),
369+ peakStrength: asScalar(r.variableValues.peakStrength),
370+ cascadeLostFrac: asScalar(r.variableValues.cascadeLostFrac),
371+ oversample,
372+ };
373+ lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
374+ }
375+
376+ // ---- output stage ----
377+ post({ type: 'progress', id, stage: 'output', frac: 0.9 });
378+ const rOut = runStage(mfiles, `
379+${pLiteral(params)}
380+[y, peakPa] = cymbal_output(P, p, ${radEntry.oversample});
381+`, { p: radEntry.p });
382+
383+ const y32 = toF32(rOut.variableValues.y);
384+ const peakPa = asScalar(rOut.variableValues.peakPa);
385+
386+ const warnings: string[] = [];
387+ if (radEntry.listenerInside) {
388+ warnings.push(
389+ 'the listener is within 2 cm of the plate; the dipole field is evaluated at a clamped minimum distance',
390+ );
391+ }
392+ if (!params.normalize) {
393+ let peak = 0;
394+ for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
395+ if (peak > 1.0) {
396+ warnings.push(
397+ `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,
398+ );
399+ }
400+ }
401+ if (params.n_azimuth < 2.5 * solveEntry.maxN) {
402+ warnings.push(
403+ `only ${params.n_azimuth} azimuth samples for ${solveEntry.maxN} nodal diameters; the listener transfer of the high-order modes is aliased`,
404+ );
405+ }
406+ if (radEntry.cascadeLostFrac > 0.25) {
407+ warnings.push(
408+ `${(100 * radEntry.cascadeLostFrac).toFixed(0)}% of the struck energy left the top of the retained band; raise max_mode_hz or lower cascade_leak`,
409+ );
410+ }
411+ const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];
412+ if (topMode < 0.6 * params.max_mode_hz) {
413+ warnings.push(
414+ `the highest retained mode is only ${topMode.toFixed(0)} Hz; raise max nodal diameters to fill the band`,
415+ );
416+ }
417+
418+ const resp: RenderResponse = {
419+ type: 'result',
420+ id,
421+ quality,
422+ y: y32,
423+ sampleRate: params.sample_rate,
424+ strikeT: radEntry.strikeT.slice(),
425+ strikeV: radEntry.strikeV.slice(),
426+ strikeS: radEntry.strikeS.slice(),
427+ modeF: solveEntry.modeF.slice(),
428+ modeZeta: solveEntry.modeZeta.slice(),
429+ modeZetaRad: solveEntry.modeZetaRad.slice(),
430+ modeSigma: solveEntry.modeSigma.slice(),
431+ modeN: solveEntry.modeN.slice(),
432+ nodeR: solveEntry.nodeR.slice(),
433+ nodeZ: solveEntry.nodeZ.slice(),
434+ shapeW: solveEntry.shapeW.slice(),
435+ cascadeDb: radEntry.cascadeDb.slice(),
436+ cascadeBands: radEntry.cascadeBands,
437+ cascadeFrames: radEntry.cascadeFrames,
438+ cascadeFLo: solveEntry.modeF[0],
439+ cascadeFHi: solveEntry.modeF[solveEntry.modeF.length - 1],
440+ cascadeTEnd: radEntry.cascadeTEnd,
441+ diag: {
442+ modesKept: solveEntry.modeF.length,
443+ maxN: solveEntry.maxN,
444+ strikeCount: radEntry.strikeCount,
445+ meanV: radEntry.meanV,
446+ maxV: radEntry.maxV,
447+ meanTauUs: radEntry.meanTauUs,
448+ peakForceN: radEntry.peakForceN,
449+ durationS: radEntry.durationS,
450+ peakPa,
451+ peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
452+ directPeakPa: radEntry.directPeakPa,
453+ listenerDistM: radEntry.listenerDistM,
454+ peakWrmsUm: radEntry.peakWrmsUm,
455+ peakStrength: radEntry.peakStrength,
456+ cascadeLostFrac: radEntry.cascadeLostFrac,
457+ plateMassKg: solveEntry.plateMassKg,
458+ warnings,
459+ solveMs,
460+ radiateMs,
461+ solveCached,
462+ radiateCached,
463+ oversample: radEntry.oversample,
464+ },
465+ };
466+ (self as unknown as Worker).postMessage(resp, [
467+ resp.y.buffer,
468+ resp.strikeT.buffer,
469+ resp.strikeV.buffer,
470+ resp.strikeS.buffer,
471+ resp.modeF.buffer,
472+ resp.modeZeta.buffer,
473+ resp.modeZetaRad.buffer,
474+ resp.modeSigma.buffer,
475+ resp.modeN.buffer,
476+ resp.nodeR.buffer,
477+ resp.nodeZ.buffer,
478+ resp.shapeW.buffer,
479+ resp.cascadeDb.buffer,
480+ ]);
481+}
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+});