Maraca shake simulator: MATLAB model and interactive web app
A few hundred grains in a closed spherical shell, drawn as a seeded
stochastic process of Hertz impacts, driving the shell's spherical-harmonic
modes and radiating to a listener through the exact sphere multipole
transfer. Structured after the sibling woodblock repository: the algorithm
in MATLAB syntax run by numbl, a Vite web app that executes it client-side
in a worker, and a validation script wired into the deploy workflow.
41 changed files+8031−0
.github/workflows/deploy.ymladded+68−0View file
@@ -0,0 +1,68 @@
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/maraca_validate.m through numbl: closed-form identities,
41+ # reproducibility of the seeded impact process, and the pinned
42+ # end-to-end diagnostics. Fails the build on any mismatch.
43+ - name: Validate the model
44+ working-directory: webapp
45+ run: npm run selftest
46+
47+ - name: Build
48+ working-directory: webapp
49+ run: npm run build
50+
51+ - name: Setup Pages
52+ uses: actions/configure-pages@v5
53+
54+ - name: Upload artifact
55+ uses: actions/upload-pages-artifact@v3
56+ with:
57+ path: webapp/dist
58+
59+ deploy:
60+ needs: build
61+ runs-on: ubuntu-latest
62+ environment:
63+ name: github-pages
64+ url: ${{ steps.deployment.outputs.page_url }}
65+ steps:
66+ - name: Deploy to Pages
67+ id: deployment
68+ uses: actions/deploy-pages@v4
README.mdadded+45−0View file
@@ -0,0 +1,45 @@
1+# Maraca
2+
3+A physically-based simulation of a maraca shake pattern, sampled as acoustic pressure at a point in space. A few hundred grains inside a closed spherical shell produce a few thousand stochastic grain-shell impacts across the strokes of the gesture; each impact is a Hertz contact that drives the shell's spherical-harmonic modes, and the modes radiate to a listener point through the exact multipole transfer of a sphere.
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 maraca_demo.m`, or `maraca_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+shake gesture -> stochastic impact process maraca_impacts.m
22+ -> Hertz grain contacts on the shell wall
23+ -> spherical shell modes (two membrane branches + bending) maraca_basis.m
24+ -> exact sphere multipole radiation maraca_radiate.m
25+ -> retarded free-field pressure at the listener
26+ -> anti-aliased decimation -> PCM maraca_output.m
27+```
28+
29+Three points are worth stating plainly, because they are where the model commits itself.
30+
31+**The grains are not tracked individually.** A stroke is a draw from a point process: a Poisson count of impacts whose mean is `grain_count * bounces_per_grain`, arrival times from a Gamma(3) burst after the hand reverses, Rayleigh-distributed impact speeds, and impact points scattered by a von Mises-Fisher distribution about the pole the grains are thrown at, alternating between the two poles on successive strokes. Grain sizes are lognormal. Each impact is then an independent Hertz contact, so the contact time falls as `v^(-1/5)` and harder or faster grains excite the shell higher up: timbre varies with dynamics on its own, rather than through a scripted velocity layer. What this buys is speed, a full pattern renders in well under a second, and reproducibility from a seed. What it gives up is any coupling between grains, so a shell packed tight enough for the grains to jam is outside the model.
32+
33+**The shell modes come from membrane theory of a complete sphere, with a bending correction.** Each degree `l` gives two branches from a quadratic in `Omega^2`; the `l = 0` upper root is the breathing mode, and the `l = 1` lower root is exactly zero, being rigid translation, which the hand holds. The lower branch flattens toward `Omega^2 = 1 - nu^2` at high `l`, so a flat-plate bending term is added to it, which is what spreads the high-degree modes upward. The interior air enters as a stiffness on the breathing mode alone, since modes with `l >= 1` displace no net volume.
34+
35+**Radiation uses the exact multipole transfer of a sphere**, `-i rho0 c h_l(kr) / h_l'(ka)`, evaluated at each mode's centre frequency. Because every `m` of a degree shares one frequency and one transfer, the spherical-harmonic addition theorem collapses the `m` sum to a Legendre polynomial in the angle between the impact point and the listener, so one damped resonator per retained mode carries the whole degree. Modal truncation then discards the dense high-degree modes whose collective off-resonance response is the immediate forced motion of the wall around the contact, and that broadband part is reinstated separately as a point-driven-plate term. That last term is the least principled piece of the model, and `direct_gain` turns it off.
36+
37+## Validation
38+
39+`matlab/maraca_validate.m` (also `npm run selftest` in `webapp/`) checks the spherical Hankel functions against their closed forms and against finite differences, the breathing mode against its closed form, the Hertz `v^(-1/5)` contact-time scaling, the far-field `-6 dB` per doubling of distance, exact reproducibility of the seeded impact process, and pinned end-to-end diagnostics.
40+
41+Unlike the sibling `woodblock` repository, there is no separate Python reference implementation here: the MATLAB code is the reference, and the pinned numbers were recorded from it when the model was built.
42+
43+## Defaults
44+
45+A 3.4 cm gourd shell, 320 grains at 1.3 bounces each, eight strokes at 4.5 Hz with an accent on every second stroke. That yields 3628 impacts with a 39 µs mean Hertz contact, 11 retained shell modes from 6446 Hz up, a breathing mode at 14751 Hz, and a peak of 14.9 Pa (117 dB SPL) at the 42 cm listener before normalization.
matlab/maraca_basis.madded+111−0View file
@@ -0,0 +1,111 @@
1+function basis = maraca_basis(P)
2+% MARACA_BASIS Mode set of a complete thin spherical shell.
3+%
4+% Membrane theory of the full sphere gives two branches per degree l from
5+% Omega^4 - (lam + 1 + 3 nu) Omega^2 + (lam - 2)(1 - nu^2) = 0,
6+% lam = l(l+1), Omega = omega R / c_p, c_p = sqrt(E / (rho (1 - nu^2))).
7+% Checks: l = 0 upper root is the breathing mode Omega^2 = 2(1+nu); the l = 1
8+% lower root is exactly 0 (rigid translation, held by the hand and neglected).
9+% The lower branch flattens toward Omega^2 = 1 - nu^2 at high l, so a bending
10+% term beta^2 (lam - 2)^2, beta^2 = h^2/(12 R^2), is added to it; this is the
11+% flat-plate asymptote and is what spreads the high-l modes upward in
12+% frequency. The membrane eigenvector is kept for the bending-corrected modes,
13+% an approximation that is good while the bending term is a correction.
14+%
15+% Eigenvector ratio (tangential potential V over radial W, W = 1):
16+% V = (2(1+nu) - Omega^2) / (lam (1+nu)),
17+% giving the modal mass factor mu = 1 + lam V^2 and modal mass rho h R^2 mu.
18+%
19+% The closed cavity enters as an air stiffness on the l = 0 (breathing)
20+% components only: d(omega^2) = 3 rho0 c^2 / (rho h R), scaled by
21+% cavity_coupling. Modes with l >= 1 displace no net volume, and the interior
22+% inertial loading they see is of order rho0 R / ((l+1) rho h), well under a
23+% percent here, so it is neglected.
24+%
25+% Damping: zeta = 0.5 (loss_factor + loss_slope omega) plus radiation loss
26+% rho0 c sigma_l / (2 rho h omega mu), sigma_l = Re[-i h_l(ka)/h_l'(ka)].
27+
28+R = P.shell_radius_m;
29+h = P.shell_thickness_m;
30+nu = P.poisson_ratio;
31+rho = P.density;
32+cp = sqrt(P.youngs_modulus / (rho * (1 - nu^2)));
33+beta2 = h^2 / (12 * R^2);
34+c = P.speed_of_sound;
35+rho0 = P.air_density;
36+
37+cavShift2 = P.cavity_coupling * 3 * rho0 * c^2 / (rho * h * R); % rad^2/s^2
38+
39+maxN = 2 * (P.max_mode_l + 1);
40+ll = zeros(maxN, 1);
41+ff = zeros(maxN, 1);
42+mus = zeros(maxN, 1);
43+n = 0;
44+breathingDryHz = 0;
45+breathingHz = 0;
46+for l = 0:P.max_mode_l
47+ lam = l * (l + 1);
48+ b = lam + 1 + 3 * nu;
49+ det2 = (lam - 2) * (1 - nu^2);
50+ disc = sqrt(max(0, b^2 - 4 * det2));
51+ om2lo = (b - disc) / 2;
52+ om2hi = (b + disc) / 2;
53+ if l == 0
54+ cand = om2hi; % breathing; lower root is negative
55+ elseif l == 1
56+ cand = om2hi; % lower root 0 = rigid translation
57+ else
58+ cand = [om2lo + beta2 * (lam - 2)^2, om2hi];
59+ end
60+ for ci = 1:length(cand)
61+ Om2 = cand(ci);
62+ om = sqrt(Om2) * cp / R * P.stiffness_trim;
63+ if l == 0
64+ breathingDryHz = om / (2 * pi);
65+ om = sqrt(om^2 + cavShift2); % closed-cavity air stiffness
66+ breathingHz = om / (2 * pi);
67+ end
68+ fm = om / (2 * pi);
69+ if fm > P.max_mode_hz
70+ continue;
71+ end
72+ if lam > 0
73+ V = (2 * (1 + nu) - Om2) / (lam * (1 + nu));
74+ else
75+ V = 0;
76+ end
77+ n = n + 1;
78+ ll(n) = l;
79+ ff(n) = fm;
80+ mus(n) = 1 + lam * V^2;
81+ end
82+end
83+if n == 0
84+ error('No shell modes below max_mode_hz -- check geometry/material.');
85+end
86+[ff, order] = sort(ff(1:n, 1));
87+ff = reshape(ff, n, 1);
88+lls = reshape(ll(order), n, 1);
89+muss = reshape(mus(order), n, 1);
90+
91+omega = 2 * pi * ff;
92+sigma = zeros(n, 1);
93+for m = 1:n
94+ ka = omega(m) / c * R;
95+ [hh, hhp] = maraca_hankel(lls(m), ka);
96+ sigma(m) = real(-1i * hh(lls(m) + 1) / hhp(lls(m) + 1));
97+end
98+zeta = 0.5 * (P.loss_factor + P.loss_slope * omega) ...
99+ + P.radiation_loss * rho0 * c * sigma ./ (2 * rho * h * omega .* muss);
100+
101+basis.n = n;
102+basis.l = lls;
103+basis.f = ff;
104+basis.omega = omega;
105+basis.zeta = zeta;
106+basis.mu = muss;
107+basis.massKg = rho * h * R^2 * muss;
108+basis.sigma = sigma;
109+basis.breathingDryHz = breathingDryHz;
110+basis.breathingHz = breathingHz;
111+end
matlab/maraca_defaults.madded+62−0View file
@@ -0,0 +1,62 @@
1+function P = maraca_defaults()
2+% MARACA_DEFAULTS Default parameters, mirroring the web app's params.ts.
3+% The output duration is not a parameter: it is
4+% stroke_count / stroke_rate_hz + tail_s.
5+
6+% Shell geometry (complete thin spherical shell)
7+P.shell_radius_m = 0.034;
8+P.shell_thickness_m = 0.0022;
9+
10+% Shell material (dried gourd / hardwood veneer scale)
11+P.density = 720.0;
12+P.youngs_modulus = 2.5e9;
13+P.poisson_ratio = 0.30;
14+P.stiffness_trim = 1.00;
15+P.max_mode_l = 16;
16+P.max_mode_hz = 22000.0;
17+
18+% Damping
19+P.loss_factor = 0.015;
20+P.loss_slope = 6.0e-7;
21+P.radiation_loss = 1.0;
22+
23+% Grains (seed beads)
24+P.grain_radius_m = 0.0015;
25+P.grain_density = 1300.0;
26+P.grain_youngs = 6.0e8;
27+P.grain_size_spread = 0.25;
28+P.restitution = 0.55;
29+
30+% Shake gesture / stochastic impact process
31+P.stroke_rate_hz = 4.5;
32+P.stroke_count = 8;
33+P.accent_period = 2;
34+P.accent_gain = 1.5;
35+P.grain_count = 320;
36+P.bounces_per_grain = 1.3;
37+P.shake_velocity = 2.2;
38+P.slosh_lag_ms = 14.0;
39+P.burst_width_ms = 45.0;
40+P.impact_spread_deg = 35.0;
41+P.seed = 1;
42+
43+% Closed air cavity
44+P.cavity_coupling = 1.0;
45+
46+% Radiation / listener
47+P.direct_gain = 1.0;
48+P.listener_x_m = 0.15;
49+P.listener_y_m = 0.25;
50+P.listener_z_m = 0.30;
51+P.speed_of_sound = 343.0;
52+P.air_density = 1.204;
53+
54+% Rendering / output
55+P.sample_rate = 48000;
56+P.oversample = 2;
57+P.tail_s = 0.40;
58+P.pre_delay_ms = 3.0;
59+P.normalize = 1;
60+P.peak_dbfs = -3.0;
61+P.full_scale_pa = 20.0;
62+end
matlab/maraca_demo.madded+26−0View file
@@ -0,0 +1,26 @@
1+% Render a maraca shake with default parameters and print diagnostics.
2+% Run with: numbl run maraca_demo.m
3+
4+P = maraca_defaults();
5+tic;
6+[y, basis, dinfo] = maraca_render(P);
7+elapsed = toc;
8+
9+fprintf('modes kept : %d\n', basis.n);
10+fprintf('mode l : ');
11+fprintf('%d ', basis.l);
12+fprintf('\nmode freq : ');
13+fprintf('%.0f ', basis.f);
14+fprintf('Hz\n');
15+fprintf('breathing : %.0f Hz (dry %.0f Hz, cavity shift %+.1f Hz)\n', ...
16+ dinfo.breathingHz, dinfo.breathingDryHz, dinfo.breathingHz - dinfo.breathingDryHz);
17+fprintf('grains : %d, %.2f bounces each per stroke\n', P.grain_count, P.bounces_per_grain);
18+fprintf('impacts : %d (mean v %.2f m/s, max %.2f; mean contact %.0f us)\n', ...
19+ dinfo.impactCount, dinfo.meanV, dinfo.maxV, dinfo.meanTauUs);
20+fprintf('peak press.: %.3f 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 inside/at the shell surface\n');
26+end
matlab/maraca_hankel.madded+54−0View file
@@ -0,0 +1,54 @@
1+function [h, hp] = maraca_hankel(lmax, x)
2+% MARACA_HANKEL Spherical Hankel functions of the second kind and derivatives.
3+% h(l+1) = h_l^(2)(x) = j_l(x) - i y_l(x), hp(l+1) = d/dx h_l^(2)(x), for
4+% l = 0..lmax, real x > 0. With the e^{+i omega t} convention used throughout,
5+% h^(2)(kr) is the outgoing wave.
6+%
7+% j_l uses downward (Miller) recursion, stable for x below l; y_l uses upward
8+% recursion, stable everywhere. Verified against the closed forms for l = 0, 1.
9+
10+if x <= 0
11+ error('maraca_hankel: x must be positive');
12+end
13+
14+% ---- j_l: downward recursion, normalized by j_0 = sin(x)/x ----
15+Lstart = lmax + 15 + ceil(x);
16+jj = zeros(Lstart + 2, 1); % jj(n+1) holds j_n up to a common scale
17+jj(Lstart + 2) = 0;
18+jj(Lstart + 1) = 1e-30;
19+for n = Lstart:-1:1
20+ jj(n) = ((2 * n + 1) / x) * jj(n + 1) - jj(n + 2);
21+ if abs(jj(n)) > 1e250 % renormalize to avoid overflow
22+ jj = jj * 1e-250;
23+ end
24+end
25+% The explicit second index works around a numbl slicing bug that can hand
26+% back jj(1:lmax+1) as a row vector depending on what else the script calls.
27+j = jj(1:lmax + 1, 1) * ((sin(x) / x) / jj(1));
28+
29+% ---- y_l: upward recursion ----
30+y = zeros(lmax + 1, 1);
31+y(1) = -cos(x) / x;
32+if lmax >= 1
33+ y(2) = -cos(x) / x^2 - sin(x) / x;
34+ for n = 1:lmax - 1
35+ y(n + 2) = ((2 * n + 1) / x) * y(n + 1) - y(n);
36+ end
37+end
38+
39+h = j - 1i * y;
40+
41+% ---- derivatives: h'_l = h_{l-1} - (l+1)/x h_l, h'_0 = -h_1 ----
42+hp = zeros(lmax + 1, 1) + 0i;
43+h1 = j(min(lmax + 1, 2)) - 1i * y(min(lmax + 1, 2));
44+if lmax == 0
45+ % need h_1 explicitly for h'_0
46+ y1 = -cos(x) / x^2 - sin(x) / x;
47+ j1 = sin(x) / x^2 - cos(x) / x;
48+ h1 = j1 - 1i * y1;
49+end
50+hp(1) = -h1;
51+for l = 1:lmax
52+ hp(l + 1) = h(l) - ((l + 1) / x) * h(l + 1);
53+end
54+end
matlab/maraca_impacts.madded+149−0View file
@@ -0,0 +1,149 @@
1+function [timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P)
2+% MARACA_IMPACTS Stochastic grain-shell impact process for a shake pattern.
3+%
4+% The shell holds grain_count grains. A stroke throws them at one wall, where
5+% each makes bounces_per_grain contacts on average (one arrival plus whatever
6+% rattle it has before the next reversal), so the expected number of impacts
7+% in stroke k is grain_count * bounces_per_grain * sqrt(accent_k). The grains
8+% are not tracked individually: at a few hundred grains in a shell barely
9+% wider than a few grain diameters, the arrival process is what is audible,
10+% not which grain arrived. That is the central approximation of this model,
11+% and it is what makes a few thousand impacts cost milliseconds.
12+%
13+% Each stroke k of the gesture produces a burst of independent grain impacts:
14+% count ~ Poisson with that mean (Gaussian approximation, exact enough
15+% at hundreds per stroke);
16+% times = stroke onset + slosh_lag + Gamma(3, burst_width/3) delay, the
17+% collective slosh of the grain mass arriving at the far wall;
18+% speeds ~ Rayleigh with mean shake_velocity * sqrt(accent_k);
19+% points ~ von Mises-Fisher about the pole the grains are thrown at, which
20+% alternates between +x and -x on successive strokes;
21+% size ~ lognormal spread of the grain radius (grain_size_spread).
22+%
23+% Each impact is a Hertz contact of a small sphere on the shell wall:
24+% k_h = (4/3) E* sqrt(r_g), 1/E* = (1-nu_g^2)/E_g + (1-nu_s^2)/E_s
25+% tau = 3.218 (m^2 / (k_h^2 v))^{1/5} (Hertz contact time)
26+% J = (1 + restitution) m v (momentum transferred)
27+% Fpk = pi J / (2 tau) (half-sine pulse peak)
28+% The force pulse is modeled as a half sine of duration tau carrying impulse J,
29+% so J, tau, and Fpk are mutually consistent by construction.
30+%
31+% All randomness comes from maraca_rng seeded with P.seed, so a given seed
32+% reproduces the identical process in the CLI, Node, and the browser.
33+%
34+% Returns row vectors sorted by time: timp (s), jimp (N s), tauimp (s),
35+% fimp (peak N), and ximp (3 x I unit impact positions on the shell).
36+% ginfo carries the expected-rate envelope and summary statistics for display.
37+
38+state = mod(floor(abs(P.seed)), 4294967296);
39+K = P.stroke_count;
40+lag = P.slosh_lag_ms * 1e-3;
41+theta = (P.burst_width_ms * 1e-3) / 3; % Gamma(3, theta): mean = width
42+spreadRad = P.impact_spread_deg * pi / 180;
43+kappa = 1 / max(1e-4, spreadRad^2);
44+duration = P.stroke_count / P.stroke_rate_hz + P.tail_s;
45+
46+% Hertz constants shared by all impacts (grain Poisson ratio fixed at 0.35)
47+nuG = 0.35;
48+Estar = 1 / ((1 - nuG^2) / P.grain_youngs + (1 - P.poisson_ratio^2) / P.youngs_modulus);
49+rg0 = P.grain_radius_m;
50+mg0 = (4 / 3) * pi * rg0^3 * P.grain_density;
51+
52+maxI = ceil(K * P.grain_count * P.bounces_per_grain * sqrt(P.accent_gain) * 1.5) + 64;
53+timp = zeros(1, maxI);
54+jimp = zeros(1, maxI);
55+tauimp = zeros(1, maxI);
56+fimp = zeros(1, maxI);
57+vimp = zeros(1, maxI);
58+ximp = zeros(3, maxI);
59+I = 0;
60+
61+counts = zeros(K, 1);
62+accents = zeros(K, 1);
63+for k = 1:K
64+ tk = (k - 1) / P.stroke_rate_hz;
65+ if mod(k - 1, max(1, round(P.accent_period))) == 0
66+ a = P.accent_gain;
67+ else
68+ a = 1.0;
69+ end
70+ accents(k) = a;
71+ meanCount = P.grain_count * P.bounces_per_grain * sqrt(a);
72+ [u, state] = maraca_rng(state, 2);
73+ z = sqrt(-2 * log(u(1))) * cos(2 * pi * u(2));
74+ nk = max(0, round(meanCount + sqrt(meanCount) * z));
75+ counts(k) = nk;
76+ if nk == 0
77+ continue;
78+ end
79+ [u, state] = maraca_rng(state, nk * 8);
80+ U = reshape(u, nk, 8);
81+
82+ % times: Gamma(3, theta) as the sum of three exponentials
83+ tt = tk + lag - theta * log(U(:, 1) .* U(:, 2) .* U(:, 3));
84+ % speeds: Rayleigh, mean shake_velocity * sqrt(a), capped at 2.5x the mean
85+ % (a grain cannot exceed a small multiple of the shell speed)
86+ sigV = P.shake_velocity * sqrt(a) / sqrt(pi / 2);
87+ vv = sigV * sqrt(-2 * log(U(:, 4)));
88+ vv = min(vv, 2.5 * P.shake_velocity * sqrt(a));
89+ % positions: von Mises-Fisher about the active pole (+x odd, -x even strokes)
90+ w = 1 + log(U(:, 5) * (1 - exp(-2 * kappa)) + exp(-2 * kappa)) / kappa;
91+ st = sqrt(max(0, 1 - w.^2));
92+ phi = 2 * pi * U(:, 6);
93+ pole = 1 - 2 * mod(k - 1, 2); % +1, -1, +1, ...
94+ xs = [pole * w.'; (st .* cos(phi)).'; (st .* sin(phi)).'];
95+ % grain size: lognormal scale on the radius
96+ zs = sqrt(-2 * log(U(:, 7))) .* cos(2 * pi * U(:, 8));
97+ s = exp(P.grain_size_spread * zs);
98+ mg = mg0 * s.^3;
99+ kh = (4 / 3) * Estar * sqrt(rg0 * s);
100+
101+ tau = 3.218 * (mg.^2 ./ (kh.^2 .* vv)).^(1 / 5);
102+ J = (1 + P.restitution) * mg .* vv;
103+ Fpk = pi * J ./ (2 * tau);
104+
105+ keep = find(tt < duration - 2e-3);
106+ nkeep = length(keep);
107+ idx = I + 1:I + nkeep;
108+ timp(idx) = tt(keep).';
109+ jimp(idx) = J(keep).';
110+ tauimp(idx) = tau(keep).';
111+ fimp(idx) = Fpk(keep).';
112+ vimp(idx) = vv(keep).';
113+ ximp(:, idx) = xs(:, keep);
114+ I = I + nkeep;
115+end
116+if I == 0
117+ error('The impact process produced no impacts -- check the gesture parameters.');
118+end
119+
120+[timp, order] = sort(timp(1, 1:I));
121+timp = reshape(timp, 1, I);
122+jimp = reshape(jimp(order), 1, I);
123+tauimp = reshape(tauimp(order), 1, I);
124+fimp = reshape(fimp(order), 1, I);
125+vimp = reshape(vimp(order), 1, I);
126+xs = ximp(:, 1:I);
127+ximp = xs(:, order);
128+
129+% Expected-rate envelope on a 1 ms grid, for plotting and validation
130+tg = 0:1e-3:duration;
131+rate = zeros(size(tg));
132+for k = 1:K
133+ tk = (k - 1) / P.stroke_rate_hz + lag;
134+ tau0 = max(0, tg - tk);
135+ g = (tau0.^2 / (2 * theta^3)) .* exp(-tau0 / theta); % Gamma(3,theta) pdf
136+ rate = rate + P.grain_count * P.bounces_per_grain * sqrt(accents(k)) * g;
137+end
138+
139+ginfo.count = I;
140+ginfo.counts = counts;
141+ginfo.accents = accents;
142+ginfo.meanV = sum(vimp) / I;
143+ginfo.maxV = max(vimp);
144+ginfo.vimp = vimp;
145+ginfo.meanTauUs = sum(tauimp) / I * 1e6;
146+ginfo.rateT = tg;
147+ginfo.rateHz = rate;
148+ginfo.durationS = duration;
149+end
matlab/maraca_output.madded+63−0View file
@@ -0,0 +1,63 @@
1+function [y, peakPa] = maraca_output(P, p, oversample)
2+% MARACA_OUTPUT Anti-aliased decimation and output conditioning.
3+% Kaiser-windowed-sinc FIR lowpass at 0.45 * fs_lo, decimate by oversample,
4+% pre-roll, tail fade, then either normalize to peak_dbfs or divide by
5+% full_scale_pa. Identical in structure to the woodblock output stage.
6+
7+peakPa = max(abs(p));
8+R = oversample;
9+T = length(p);
10+
11+if R > 1
12+ taps = 12 * R + 1;
13+ M2 = (taps - 1) / 2;
14+ fcNorm = 0.45 / R; % cycles per input sample
15+ m = (0:taps - 1) - M2;
16+ h = sin(2 * pi * fcNorm * m) ./ (pi * m);
17+ h(M2 + 1) = 2 * fcNorm; % sinc value at m = 0
18+ r = m / M2;
19+ beta = 8;
20+ h = h .* (mr_i0(beta * sqrt(max(0, 1 - r.^2))) / mr_i0(beta));
21+ h = h / sum(h); % unity DC gain
22+ pf = conv(p, h);
23+ outLen = floor(T / R);
24+ y = reshape(pf((0:outLen - 1) * R + M2 + 1), 1, outLen); % centered pick
25+else
26+ y = p;
27+end
28+
29+% Pre-roll, then truncate to the nominal duration
30+padN = floor(P.pre_delay_ms * 1e-3 * P.sample_rate);
31+duration = P.stroke_count / P.stroke_rate_hz + P.tail_s;
32+total = floor(duration * P.sample_rate);
33+if padN > 0
34+ y = [zeros(1, padN) y];
35+end
36+y = y(1:min(length(y), total));
37+
38+% Squared fade on the tail
39+fade = min(256, floor(length(y) / 8));
40+if fade > 1
41+ y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;
42+end
43+
44+if P.normalize ~= 0
45+ pk = max(abs(y));
46+ if pk > 0
47+ y = y / pk * 10^(P.peak_dbfs / 20);
48+ end
49+else
50+ y = y / P.full_scale_pa;
51+end
52+end
53+
54+function v = mr_i0(x)
55+% Modified Bessel function of the first kind, order zero (series expansion).
56+v = ones(size(x));
57+term = ones(size(x));
58+hx = x / 2;
59+for k = 1:40
60+ term = term .* (hx / k).^2;
61+ v = v + term;
62+end
63+end
matlab/maraca_radiate.madded+142−0View file
@@ -0,0 +1,142 @@
1+function [p, radinfo] = maraca_radiate(P, basis, timp, jimp, tauimp, fimp, ximp, fsHi)
2+% MARACA_RADIATE Modal synthesis and free-field pressure at the listener.
3+%
4+% Two contributions:
5+%
6+% 1. Modal ringing. An impact with impulse J at surface point x_i gives every
7+% mode (l, m) a velocity jump J S(omega tau) Y_lm(x_i) / m_l, where
8+% S(x) = cos(x/2) / (1 - (x/pi)^2) is the spectrum of the half-sine contact
9+% pulse (harder or faster grains excite higher modes). At mode frequency
10+% the sphere radiates each (l, m) with the exact multipole transfer
11+% -i rho0 c h_l(kr) / h_l'(ka); modes are narrow-band, so the transfer is
12+% evaluated at the mode centre frequency. Because all m of a degree share
13+% one frequency and one transfer, the addition theorem
14+% sum_m Y_lm(x_i) Y_lm(x_r) = (2l+1)/(4pi) P_l(cos gamma)
15+% collapses the m sum: only one damped resonator per retained mode is
16+% synthesized, carrying listener-projected weights. The listener-direction
17+% projection therefore lives in this stage, and moving the listener re-runs
18+% only this stage. The carrier phase e^{-ikr} is stripped from h_l(kr) and
19+% the propagation applied as a common group delay, so envelope and phase
20+% are not double counted.
21+%
22+% 2. Direct impact clicks. Modal truncation at max_mode_hz discards the dense
23+% high-l modes whose collective off-resonance response is the immediate
24+% forced motion of the wall around the contact point. That broadband
25+% component is reinstated with the classic point-driven infinite plate
26+% result p = rho0 F(t - r/c) / (2 pi rho_s h r), scaled by direct_gain,
27+% weighted (1 + cos gamma)/2 for shadowing of impacts on the far side, and
28+% advanced/delayed by -R cos(gamma)/c for the position on the sphere. This
29+% is the least principled term of the model: the plate formula ignores the
30+% shell's curvature and closure.
31+%
32+% The synthesis is exact in time: between impacts every resonator evolves as
33+% exp(lambda t), evaluated per inter-impact segment; nothing is integrated.
34+%
35+% Returns p (1 x T, pascals at fsHi) and radinfo diagnostics.
36+
37+c = P.speed_of_sound;
38+rho0 = P.air_density;
39+R = P.shell_radius_m;
40+duration = P.stroke_count / P.stroke_rate_hz + P.tail_s;
41+T = floor(duration * fsHi);
42+N = basis.n;
43+I = length(timp);
44+
45+lx = P.listener_x_m; ly = P.listener_y_m; lz = P.listener_z_m;
46+rmag = sqrt(lx^2 + ly^2 + lz^2);
47+listenerInside = rmag < R + 0.02;
48+rmag = max(rmag, R + 0.005);
49+xhat = [lx; ly; lz] / max(1e-9, sqrt(lx^2 + ly^2 + lz^2));
50+
51+% ---- listener geometry per impact: Legendre table P_l(cos gamma) ----
52+cosg = xhat.' * ximp; % 1 x I
53+lmax = max(basis.l);
54+PL = zeros(lmax + 1, I);
55+PL(1, :) = ones(1, I);
56+if lmax >= 1
57+ PL(2, :) = cosg;
58+ for l = 1:lmax - 1
59+ PL(l + 2, :) = ((2 * l + 1) * cosg .* PL(l + 1, :) - l * PL(l, :)) / (l + 1);
60+ end
61+end
62+
63+% ---- per-mode radiation transfer and complex pole ----
64+G = zeros(N, 1) + 0i;
65+lam = zeros(N, 1) + 0i;
66+phv = zeros(N, 1) + 0i; % velocity phasor per unit jump
67+for m = 1:N
68+ l = basis.l(m);
69+ om = basis.omega(m);
70+ zt = basis.zeta(m);
71+ k = om / c;
72+ [ha, hpa] = maraca_hankel(l, k * R);
73+ [hr, hpr] = maraca_hankel(l, k * rmag);
74+ hred = hr(l + 1) * exp(1i * k * rmag); % strip carrier: no double delay
75+ G(m) = -1i * rho0 * c * hred / hpa(l + 1);
76+ omd = om * sqrt(max(1e-12, 1 - zt^2));
77+ lam(m) = -zt * om + 1i * omd;
78+ phv(m) = 1 + 1i * zt * om / omd; % v(t) = Re[phv e^{lam t}]
79+end
80+
81+% ---- per-impact modal weights: A(n, i) velocity-jump phasors ----
82+x = basis.omega * tauimp; % N x 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+Cl = ((2 * basis.l + 1) / (4 * pi)) .* (1 ./ basis.massKg); % N x 1
88+A = (phv * jimp) .* S .* (Cl * ones(1, I)) .* PL(basis.l + 1, :);
89+
90+% ---- segment synthesis: exact damped-exponential evolution ----
91+ks = min(T, max(1, round(timp * fsHi) + 1)); % impact sample indices, sorted
92+p = zeros(1, T);
93+s = zeros(N, 1) + 0i;
94+Gt = G.';
95+dt = 1 / fsHi;
96+prevK = ks(1);
97+i = 1;
98+while i <= I
99+ % accumulate all impacts landing on this sample
100+ s = s + A(:, i);
101+ while i < I && ks(i + 1) == ks(i)
102+ i = i + 1;
103+ s = s + A(:, i);
104+ end
105+ if i < I
106+ nextK = ks(i + 1);
107+ else
108+ nextK = T + 1;
109+ end
110+ len = nextK - prevK;
111+ if len > 0
112+ seg = real(Gt * ((s * ones(1, len)) .* exp(lam * ((0:len - 1) * dt))));
113+ p(prevK:nextK - 1) = p(prevK:nextK - 1) + seg;
114+ s = s .* exp(lam * (len * dt));
115+ end
116+ prevK = nextK;
117+ i = i + 1;
118+end
119+
120+% ---- direct impact clicks (point-driven plate approximation) ----
121+mpp = P.density * P.shell_thickness_m; % shell mass per unit area
122+dirScale = P.direct_gain * rho0 / (2 * pi * mpp * rmag);
123+directPeak = 0;
124+for i = 1:I
125+ nP = max(2, round(tauimp(i) * fsHi));
126+ k0 = ks(i) + round(-R * cosg(i) / c * fsHi); % near side arrives earlier
127+ if k0 < 1 || k0 + nP > T
128+ continue;
129+ end
130+ amp = dirScale * fimp(i) * 0.5 * (1 + cosg(i));
131+ pulse = amp * sin(pi * (0:nP) / nP);
132+ p(k0:k0 + nP) = p(k0:k0 + nP) + pulse;
133+ if amp > directPeak
134+ directPeak = amp;
135+ end
136+end
137+
138+radinfo.gainAbs = abs(G);
139+radinfo.directPeakPa = directPeak;
140+radinfo.listenerInside = listenerInside;
141+radinfo.listenerDistM = rmag;
142+end
matlab/maraca_render.madded+25−0View file
@@ -0,0 +1,25 @@
1+function [y, basis, dinfo] = maraca_render(P)
2+% MARACA_RENDER Full render: basis + impacts -> radiate -> output.
3+% Convenience composition of the stage functions; the web app calls the
4+% stages separately so it can cache intermediates.
5+
6+basis = maraca_basis(P);
7+fsHi = P.sample_rate * P.oversample;
8+[timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P);
9+[p, radinfo] = maraca_radiate(P, basis, timp, jimp, tauimp, fimp, ximp, fsHi);
10+[y, peakPa] = maraca_output(P, p, P.oversample);
11+
12+dinfo.impactCount = ginfo.count;
13+dinfo.meanV = ginfo.meanV;
14+dinfo.maxV = ginfo.maxV;
15+dinfo.meanTauUs = ginfo.meanTauUs;
16+dinfo.durationS = ginfo.durationS;
17+dinfo.breathingDryHz = basis.breathingDryHz;
18+dinfo.breathingHz = basis.breathingHz;
19+dinfo.peakPa = peakPa;
20+dinfo.peakSplDb = 20 * log10(peakPa / 20e-6 + 1e-30);
21+dinfo.directPeakPa = radinfo.directPeakPa;
22+dinfo.listenerInside = radinfo.listenerInside;
23+dinfo.listenerDistM = radinfo.listenerDistM;
24+dinfo.ginfo = ginfo;
25+end
matlab/maraca_rng.madded+12−0View file
@@ -0,0 +1,12 @@
1+function [r, state] = maraca_rng(state, n)
2+% MARACA_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+% impact process 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/maraca_validate.madded+107−0View file
@@ -0,0 +1,107 @@
1+% Validate the maraca model. There is no external reference implementation,
2+% so the checks are (a) closed-form identities the physics must satisfy,
3+% (b) exact reproducibility of the seeded impact process, and (c) pinned
4+% end-to-end diagnostics from the default configuration, recorded when the
5+% model was built and treated as a regression reference since.
6+% Run with: numbl run maraca_validate.m
7+
8+fails = 0;
9+
10+% ---- 1. spherical Hankel functions against closed forms ----
11+x = 3.7;
12+[h, hp] = maraca_hankel(4, x);
13+h0c = sin(x) / x - 1i * (-cos(x) / x);
14+h1c = (sin(x) / x^2 - cos(x) / x) - 1i * (-cos(x) / x^2 - sin(x) / x);
15+e1 = max(abs(h(1) - h0c), abs(h(2) - h1c));
16+d = 1e-6;
17+[ha, ~] = maraca_hankel(4, x + d);
18+[hb, ~] = maraca_hankel(4, x - d);
19+e2 = max(abs(hp - (ha - hb) / (2 * d)));
20+xs = 0.5; % small argument: Miller recursion
21+[hs, ~] = maraca_hankel(2, xs);
22+j2c = (3 / xs^3 - 1 / xs) * sin(xs) - 3 / xs^2 * cos(xs);
23+e3 = abs(real(hs(3)) - j2c);
24+ok = e1 < 1e-12 && e2 < 1e-6 && e3 < 1e-12;
25+fails = fails + (~ok);
26+fprintf('hankel closed forms: errs %.1e %.1e %.1e %s\n', e1, e2, e3, verdict(ok));
27+
28+% ---- 2. breathing mode: membrane closed form + cavity stiffness ----
29+P = maraca_defaults();
30+basis = maraca_basis(P);
31+nu = P.poisson_ratio;
32+cp = sqrt(P.youngs_modulus / (P.density * (1 - nu^2)));
33+fdry = sqrt(2 * (1 + nu)) * cp / (2 * pi * P.shell_radius_m) * P.stiffness_trim;
34+cav = P.cavity_coupling * 3 * P.air_density * P.speed_of_sound^2 ...
35+ / (P.density * P.shell_thickness_m * P.shell_radius_m);
36+fwet = sqrt((2 * pi * fdry)^2 + cav) / (2 * pi);
37+ok = abs(basis.breathingHz - fwet) < 1e-6 * fwet && all(basis.zeta > 0);
38+fails = fails + (~ok);
39+fprintf('breathing mode: %.1f Hz (closed form %.1f) %s\n', basis.breathingHz, fwet, verdict(ok));
40+
41+% ---- 3. seeded impact process is exactly reproducible ----
42+[timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P);
43+ok = ginfo.count == 3628 ...
44+ && abs(ginfo.meanV - 2.462374) < 1e-5 ...
45+ && abs(sum(jimp) - 3.586588617e-01) < 1e-9 ...
46+ && abs(timp(1) - 0.016172179) < 1e-8;
47+fails = fails + (~ok);
48+fprintf('impact process (seed %d): %d impacts, mean v %.4f m/s, sum J %.6f N s %s\n', ...
49+ P.seed, ginfo.count, ginfo.meanV, sum(jimp), verdict(ok));
50+
51+% ---- 4. Hertz contact-time scaling tau ~ v^(-1/5) ----
52+Q = P;
53+Q.grain_count = 150;
54+[~, ~, tau1, ~, ~, g1] = maraca_impacts(Q);
55+Q.shake_velocity = 2 * P.shake_velocity;
56+[~, ~, tau2, ~, ~, g2] = maraca_impacts(Q);
57+ratio = (sum(tau2) / length(tau2)) / (sum(tau1) / length(tau1));
58+ok = abs(ratio - 2^(-1 / 5)) < 0.03; % same seed, so same draws
59+fails = fails + (~ok);
60+fprintf('hertz scaling: mean tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), verdict(ok));
61+
62+% ---- 5. far-field falloff: -6 dB per doubling of distance ----
63+% RMS, not peak: the per-mode radiation gains scale exactly as 1/r, but the
64+% near-field phase terms of h_l shift the interference at the peak sample by
65+% a few tenths of a dB between distances, which is physics rather than error.
66+Q = P;
67+Q.normalize = 0;
68+Q.listener_x_m = 0; Q.listener_y_m = 0; Q.listener_z_m = 1.0;
69+[y1, ~, ~] = maraca_render(Q);
70+Q.listener_z_m = 2.0;
71+[y2, ~, ~] = maraca_render(Q);
72+drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
73+ok = abs(drop + 6.02) < 0.15;
74+fails = fails + (~ok);
75+fprintf('falloff 1 m -> 2 m: %.2f dB rms (target -6.02) %s\n', drop, verdict(ok));
76+
77+% ---- 6. pinned end-to-end diagnostics at the defaults ----
78+[y, basis, dinfo] = maraca_render(P);
79+ok = abs(dinfo.peakPa - 14.935221) < 0.02 ...
80+ && abs(basis.f(1) - 6445.5174) < 0.1 ...
81+ && basis.n == 11;
82+fails = fails + (~ok);
83+fprintf('pinned render: peak %.3f Pa (target 14.935), mode1 %.1f Hz (target 6445.5), %d modes %s\n', ...
84+ dinfo.peakPa, basis.f(1), basis.n, verdict(ok));
85+
86+% ---- 7. the tail decays to silence ----
87+fs = P.sample_rate;
88+n = length(y);
89+tail = y(n - round(0.08 * fs):n);
90+tailRms = sqrt(sum(tail.^2) / length(tail));
91+ok = tailRms < 0.02 * max(abs(y));
92+fails = fails + (~ok);
93+fprintf('tail decay: rms %.2e of peak %s\n', tailRms / max(abs(y)), verdict(ok));
94+
95+if fails > 0
96+ fprintf('%d CHECK(S) FAILED\n', fails);
97+else
98+ fprintf('all checks passed\n');
99+end
100+
101+function s = verdict(ok)
102+if ok
103+ s = 'ok';
104+else
105+ s = 'MISMATCH';
106+end
107+end
webapp/.gitignoreadded+2−0View file
@@ -0,0 +1,2 @@
1+node_modules/
2+dist/
webapp/README.mdadded+45−0View file
@@ -0,0 +1,45 @@
1+# Maraca 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 seeded stochastic process of grain-shell impacts, Hertz contacts driving the spherical-harmonic modes of a closed shell, exact sphere multipole radiation 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/maraca_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 maraca_demo.m` (or `maraca_validate.m`).
17+
18+## Why it renders fast
19+
20+Nothing is integrated. Between impacts every mode is a damped exponential, so the synthesis evaluates `exp(lambda t)` over each inter-impact segment in closed form and adds each impact's velocity jump as it lands. A few thousand impacts across a dozen modes is a few thousand short matrix products, which is why a two-second pattern renders in a fraction of a second rather than needing an ODE solver over 200 000 samples. The cost that remains is the oversampled time axis, and `oversample` trades it against alias rejection in the decimation.
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. Changing the Rayleigh speed distribution in `maraca_impacts.m` to a constant, for instance, makes the shake sound immediately mechanical.
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: shell, grain, and gesture parameters invalidate the mode basis and the impact process; the listener position and the direct-click gain only re-run radiation; pre-delay, normalization, and bit depth only re-run the output stage. Intermediate tensors stay inside the worker as numbl RuntimeValues. Moving the microphone after a full render reuses the cached impact process and returns a full-quality result.
31+
32+The New seed button redraws the whole impact process without touching any other parameter, which is the quickest way to hear how much of a shake's character is the particular draw and how much is the physics. The Loop button repeats the pattern 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/maraca_defaults.m` mirrors).
33+
34+## Layout
35+
36+- `../matlab/` — the algorithm: `maraca_basis.m`, `maraca_impacts.m`, `maraca_radiate.m`, `maraca_output.m`, `maraca_hankel.m`, `maraca_rng.m`, plus `maraca_render.m` / `maraca_demo.m` / `maraca_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/components/` — 3D scene (three.js), parameter panel, code panel and editor, plots (impact raster, waveform, spectrum, spectrogram, shell mode map)
41+- `src/selftest.ts` — Node harness for the MATLAB validation script
42+
43+## A note on the seeded generator
44+
45+The impact process draws from a 32-bit LCG in `maraca_rng.m` rather than MATLAB's `rand`, so a given seed produces the identical process in the numbl CLI, in Node, and in the browser. That is what lets the self-test pin exact impact counts and momentum sums.
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>Maraca 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": "maraca-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+272−0View file
@@ -0,0 +1,272 @@
1+import { useCallback, useEffect, useRef, useState } from 'react';
2+import { DEFAULTS, loadParams, saveParams, type MaracaParams } 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 ImpactRaster from './components/plots/ImpactRaster';
14+import ModeMap from './components/plots/ModeMap';
15+
16+export default function App() {
17+ const [params, setParams] = useState<MaracaParams>(() => loadParams());
18+ // Editable copy of the MATLAB implementation. Deliberately not persisted:
19+ // a reload restores the code the app shipped with.
20+ const [mfiles, setMfiles] = useState<MFile[]>(() => DEFAULT_MFILES.map((f) => ({ ...f })));
21+ const [codeOpen, setCodeOpen] = useState(false);
22+ const { result, error, busy, stage, progress, renderFull } = useRenderer(params, mfiles);
23+ const audioCtxRef = useRef<AudioContext | null>(null);
24+
25+ useEffect(() => saveParams(params), [params]);
26+
27+ const onChange = useCallback((key: keyof MaracaParams, value: number | boolean) => {
28+ setParams((p) => ({ ...p, [key]: value }));
29+ }, []);
30+
31+ const reset = useCallback(() => setParams({ ...DEFAULTS }), []);
32+
33+ // A new seed redraws the whole impact process without touching anything
34+ // else, which is the quickest way to hear how much of the character is the
35+ // particular draw and how much is the parameters.
36+ const reseed = useCallback(
37+ () => setParams((p) => ({ ...p, seed: 1 + Math.floor(Math.random() * 9999) })),
38+ [],
39+ );
40+
41+ const onChangeFile = useCallback((name: string, source: string) => {
42+ setMfiles((fs) => fs.map((f) => (f.name === name ? { ...f, source } : f)));
43+ }, []);
44+
45+ const onResetFile = useCallback((name: string) => {
46+ setMfiles((fs) =>
47+ fs.map((f) =>
48+ f.name === name
49+ ? { ...f, source: DEFAULT_MFILES.find((d) => d.name === name)?.source ?? f.source }
50+ : f,
51+ ),
52+ );
53+ }, []);
54+
55+ const onResetAllFiles = useCallback(() => setMfiles(DEFAULT_MFILES.map((f) => ({ ...f }))), []);
56+ const codeEdited = mfiles.some(
57+ (f) => f.source !== DEFAULT_MFILES.find((d) => d.name === f.name)?.source,
58+ );
59+
60+ const [looping, setLooping] = useState(false);
61+ const loopSrcRef = useRef<AudioBufferSourceNode | null>(null);
62+
63+ const getAudioCtx = useCallback((sampleRate: number) => {
64+ let ctx = audioCtxRef.current;
65+ if (!ctx || ctx.state === 'closed') {
66+ ctx = new AudioContext({ sampleRate });
67+ audioCtxRef.current = ctx;
68+ }
69+ void ctx.resume();
70+ return ctx;
71+ }, []);
72+
73+ const makeBuffer = useCallback((ctx: AudioContext, r: NonNullable<typeof result>) => {
74+ const buf = ctx.createBuffer(1, r.y.length, r.sampleRate);
75+ buf.copyToChannel(new Float32Array(r.y), 0);
76+ return buf;
77+ }, []);
78+
79+ const play = useCallback(() => {
80+ if (!result) return;
81+ const ctx = getAudioCtx(result.sampleRate);
82+ const src = ctx.createBufferSource();
83+ src.buffer = makeBuffer(ctx, result);
84+ src.connect(ctx.destination);
85+ src.start();
86+ }, [result, getAudioCtx, makeBuffer]);
87+
88+ // Loop mode: repeat the shake continuously, swapping in each new render as
89+ // parameters change so edits are audible without pressing Play.
90+ useEffect(() => {
91+ if (!looping || !result) return;
92+ const ctx = getAudioCtx(result.sampleRate);
93+ const src = ctx.createBufferSource();
94+ src.buffer = makeBuffer(ctx, result);
95+ src.loop = true;
96+ src.connect(ctx.destination);
97+ src.start();
98+ loopSrcRef.current = src;
99+ return () => {
100+ loopSrcRef.current = null;
101+ try {
102+ src.stop();
103+ } catch {
104+ // already stopped
105+ }
106+ src.disconnect();
107+ };
108+ }, [looping, result, getAudioCtx, makeBuffer]);
109+
110+ const download = useCallback(() => {
111+ if (!result) return;
112+ const blob = encodeWav(result.y, result.sampleRate, params.bit_depth === 16 ? 16 : 24);
113+ const url = URL.createObjectURL(blob);
114+ const a = document.createElement('a');
115+ a.href = url;
116+ a.download = result.quality === 'full' ? 'maraca.wav' : 'maraca_preview.wav';
117+ a.click();
118+ setTimeout(() => URL.revokeObjectURL(url), 5000);
119+ }, [result, params.bit_depth]);
120+
121+ return (
122+ <div className="app">
123+ <div className="topbar">
124+ <h1>Maraca Simulator</h1>
125+ <button className="primary" onClick={renderFull} disabled={busy}>
126+ Render
127+ </button>
128+ <button
129+ className={`play ${result?.quality === 'preview' ? 'preview' : ''}`}
130+ onClick={play}
131+ disabled={!result}
132+ title={
133+ result?.quality === 'preview'
134+ ? 'Playing the fast preview render; press Render for full quality'
135+ : 'Play the rendered shake'
136+ }
137+ >
138+ ▶ Play{result?.quality === 'preview' ? ' (preview)' : ''}
139+ </button>
140+ <button
141+ className={`loop ${looping ? 'active' : ''}`}
142+ onClick={() => setLooping((v) => !v)}
143+ disabled={!result && !looping}
144+ title="Repeat the shake continuously so parameter changes are audible immediately"
145+ >
146+ ⟳ Loop{looping ? ' on' : ''}
147+ </button>
148+ <button onClick={reseed} title="Draw a new impact process from a new random seed">
149+ ⚄ New seed
150+ </button>
151+ <button onClick={download} disabled={!result}>
152+ Download .wav{result?.quality === 'preview' ? ' (preview)' : ''}
153+ </button>
154+ <button onClick={reset}>Reset to defaults</button>
155+ <button
156+ className={codeOpen ? 'code active' : 'code'}
157+ onClick={() => setCodeOpen((v) => !v)}
158+ title="View and edit the MATLAB source that numbl executes"
159+ >
160+ {'{ }'} Code{codeEdited ? ' •' : ''}
161+ </button>
162+ <div className="status">
163+ {error && <span className="badge err">code error</span>}
164+ {result && (
165+ <span className={`badge ${result.quality}`}>
166+ {result.quality === 'full' ? 'full quality' : 'preview quality — press Render'}
167+ </span>
168+ )}
169+ {busy ? (
170+ <>
171+ <span>{stage}…</span>
172+ <div className="progressbar">
173+ <div style={{ width: `${Math.round(progress * 100)}%` }} />
174+ </div>
175+ </>
176+ ) : (
177+ <span>ready</span>
178+ )}
179+ </div>
180+ </div>
181+ <div className="main">
182+ <div className="sidebar">
183+ <ParamPanel params={params} onChange={onChange} />
184+ </div>
185+ <div className="content">
186+ {error && (
187+ <div className="card wide error-card">
188+ <h2>MATLAB error</h2>
189+ <div>
190+ <strong>
191+ {error.file ? `${error.file}${error.line ? `:${error.line}` : ''}` : 'error'}
192+ </strong>{' '}
193+ {error.message}
194+ </div>
195+ <div className="code-hint">
196+ Showing the last successful render. Open Code to fix the source, or revert it.
197+ </div>
198+ </div>
199+ )}
200+ <div className="card wide">
201+ <h2>Geometry — shell, grain impacts, microphone</h2>
202+ <Scene3D
203+ params={params}
204+ impactPos={result?.impactPos ?? null}
205+ impactV={result?.impactV ?? null}
206+ impactT={result?.impactT ?? null}
207+ />
208+ </div>
209+ {result && (
210+ <>
211+ <div className="card">
212+ <h2>Diagnostics</h2>
213+ <DiagnosticsPanel result={result} />
214+ </div>
215+ <div className="card">
216+ <h2>Impact process</h2>
217+ <ImpactRaster
218+ impactT={result.impactT}
219+ impactV={result.impactV}
220+ rateT={result.rateT}
221+ rateHz={result.rateHz}
222+ durationS={result.diag.durationS}
223+ impactCount={result.diag.impactCount}
224+ plotted={result.impactT.length}
225+ />
226+ </div>
227+ <div className="card">
228+ <h2>Waveform</h2>
229+ <WaveformPlot y={result.y} sampleRate={result.sampleRate} />
230+ </div>
231+ <div className="card">
232+ <h2>Spectrum</h2>
233+ <SpectrumPlot
234+ y={result.y}
235+ sampleRate={result.sampleRate}
236+ modeF={result.modeF}
237+ breathingHz={result.diag.breathingHz}
238+ />
239+ </div>
240+ <div className="card">
241+ <h2>Spectrogram</h2>
242+ <SpectrogramPlot y={result.y} sampleRate={result.sampleRate} />
243+ </div>
244+ <div className="card">
245+ <h2>Shell modes</h2>
246+ <ModeMap
247+ modeF={result.modeF}
248+ modeL={result.modeL}
249+ modeZeta={result.modeZeta}
250+ modeSigma={result.modeSigma}
251+ breathingHz={result.diag.breathingHz}
252+ />
253+ </div>
254+ </>
255+ )}
256+ </div>
257+ </div>
258+ {codeOpen && (
259+ <CodePanel
260+ mfiles={mfiles}
261+ onChangeFile={onChangeFile}
262+ onResetFile={onResetFile}
263+ onResetAll={onResetAllFiles}
264+ onClose={() => setCodeOpen(false)}
265+ onRenderFull={renderFull}
266+ error={error}
267+ busy={busy}
268+ />
269+ )}
270+ </div>
271+ );
272+}
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 maracaHighlight = 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(maracaHighlight),
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+64−0View file
@@ -0,0 +1,64 @@
1+import type { RenderResponse } from '../worker/renderWorker';
2+
3+export default function DiagnosticsPanel({ result }: { result: RenderResponse }) {
4+ const d = result.diag;
5+ const modes = Array.from(result.modeF);
6+ const shown = modes.slice(0, 10);
7+ return (
8+ <div className="diag">
9+ <table>
10+ <tbody>
11+ <tr>
12+ <th>Impacts</th>
13+ <td>
14+ {d.impactCount} · mean {d.meanV.toFixed(2)} m/s, max {d.maxV.toFixed(2)} m/s
15+ </td>
16+ </tr>
17+ <tr>
18+ <th>Contact time</th>
19+ <td>{d.meanTauUs.toFixed(0)} µs mean (Hertz)</td>
20+ </tr>
21+ <tr>
22+ <th>Shell modes</th>
23+ <td className="partials">
24+ {shown.map((f) => f.toFixed(0)).join(', ')}
25+ {modes.length > shown.length ? `, … (${modes.length} total)` : ''} Hz
26+ </td>
27+ </tr>
28+ <tr>
29+ <th>Breathing</th>
30+ <td>
31+ {d.breathingHz.toFixed(0)} Hz (dry {d.breathingDryHz.toFixed(0)} Hz, cavity{' '}
32+ {(d.breathingHz - d.breathingDryHz >= 0 ? '+' : '') +
33+ (d.breathingHz - d.breathingDryHz).toFixed(0)}{' '}
34+ Hz)
35+ </td>
36+ </tr>
37+ <tr>
38+ <th>Listener</th>
39+ <td>{(d.listenerDistM * 100).toFixed(0)} cm from the shell centre</td>
40+ </tr>
41+ <tr>
42+ <th>Peak pressure</th>
43+ <td>
44+ {d.peakPa.toFixed(3)} Pa ({d.peakSplDb.toFixed(1)} dB SPL)
45+ </td>
46+ </tr>
47+ <tr>
48+ <th>Length</th>
49+ <td>
50+ {d.durationS.toFixed(2)} s · {result.quality} · oversample {d.oversample}
51+ </td>
52+ </tr>
53+ </tbody>
54+ </table>
55+ <div className="cachenote">
56+ impacts: {d.solveCached ? 'cached' : `${(d.solveMs / 1000).toFixed(2)} s`} · radiation:{' '}
57+ {d.radiateCached ? 'cached' : `${(d.radiateMs / 1000).toFixed(2)} s`}
58+ </div>
59+ {d.warnings.map((w, i) => (
60+ <div key={i} className="warning">⚠ {w}</div>
61+ ))}
62+ </div>
63+ );
64+}
webapp/src/components/ParamPanel.tsxadded+123−0View file
@@ -0,0 +1,123 @@
1+import { PARAM_GROUPS, type ParamSpec, type MaracaParams } from '../params';
2+
3+interface Props {
4+ params: MaracaParams;
5+ onChange: (key: keyof MaracaParams, 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+271−0View file
@@ -0,0 +1,271 @@
1+// 3D view of the geometry: the shell, its handle, the grain impact points
2+// scattered over the surface, and the listener. Model coordinates put +x
3+// along the shake axis, so the two alternating impact poles sit at +/-x;
4+// three.js is y-up, so we map (x, y, z)_model -> (x, z, -y)_three, a proper
5+// rotation. Impact points are coloured by stroke so the alternation between
6+// the two poles is visible, and their size follows the impact speed.
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 { MaracaParams } from '../params';
12+
13+function m2t(x: number, y: number, z: number): [number, number, number] {
14+ return [x, z, -y];
15+}
16+
17+function makeLabel(text: string, color: string, scale = 1): THREE.Sprite {
18+ const canvas = document.createElement('canvas');
19+ canvas.width = 320;
20+ canvas.height = 64;
21+ const ctx = canvas.getContext('2d')!;
22+ ctx.font = '600 36px system-ui, sans-serif';
23+ ctx.fillStyle = color;
24+ ctx.textAlign = 'center';
25+ ctx.textBaseline = 'middle';
26+ ctx.fillText(text, 160, 32);
27+ const tex = new THREE.CanvasTexture(canvas);
28+ const mat = new THREE.SpriteMaterial({ map: tex, depthTest: false, transparent: true });
29+ const sprite = new THREE.Sprite(mat);
30+ sprite.scale.set(0.075 * scale, 0.015 * scale, 1);
31+ return sprite;
32+}
33+
34+interface Props {
35+ params: MaracaParams;
36+ impactPos: Float32Array | null;
37+ impactV: Float32Array | null;
38+ impactT: Float32Array | null;
39+}
40+
41+export default function Scene3D({ params, impactPos, impactV, impactT }: Props) {
42+ const mountRef = useRef<HTMLDivElement | null>(null);
43+ const sceneRef = useRef<{
44+ scene: THREE.Scene;
45+ group: THREE.Group | null;
46+ renderer: THREE.WebGLRenderer;
47+ camera: THREE.PerspectiveCamera;
48+ } | null>(null);
49+
50+ // one-time scene setup
51+ useEffect(() => {
52+ const mount = mountRef.current!;
53+ const renderer = new THREE.WebGLRenderer({ antialias: true });
54+ renderer.setPixelRatio(window.devicePixelRatio);
55+ renderer.setClearColor(0x14161a);
56+ mount.appendChild(renderer.domElement);
57+
58+ const scene = new THREE.Scene();
59+ const camera = new THREE.PerspectiveCamera(45, 1, 0.005, 20);
60+ camera.position.set(0.13, 0.1, 0.19);
61+
62+ const controls = new OrbitControls(camera, renderer.domElement);
63+ controls.target.set(0, 0.0, 0);
64+ controls.enableDamping = true;
65+
66+ scene.add(new THREE.AmbientLight(0xffffff, 0.55));
67+ const dir = new THREE.DirectionalLight(0xffffff, 1.1);
68+ dir.position.set(0.5, 1.0, 0.6);
69+ scene.add(dir);
70+ const dir2 = new THREE.DirectionalLight(0x8899ff, 0.3);
71+ dir2.position.set(-0.6, 0.3, -0.5);
72+ scene.add(dir2);
73+
74+ const grid = new THREE.GridHelper(0.6, 24, 0x3a4048, 0x262b33);
75+ grid.position.y = -0.09;
76+ scene.add(grid);
77+
78+ // model axes: x red (shake axis), y green, z blue
79+ const axisLen = 0.07;
80+ const axes: [THREE.Vector3, number, string][] = [
81+ [new THREE.Vector3(...m2t(axisLen, 0, 0)), 0xe57373, '+x'],
82+ [new THREE.Vector3(...m2t(0, axisLen, 0)), 0x81c784, '+y'],
83+ [new THREE.Vector3(...m2t(0, 0, axisLen)), 0x64b5f6, '+z'],
84+ ];
85+ for (const [v, color, label] of axes) {
86+ const geo = new THREE.BufferGeometry().setFromPoints([new THREE.Vector3(0, 0, 0), v]);
87+ scene.add(new THREE.Line(geo, new THREE.LineBasicMaterial({ color })));
88+ const s = makeLabel(label, '#' + color.toString(16).padStart(6, '0'));
89+ s.position.copy(v).multiplyScalar(1.15);
90+ scene.add(s);
91+ }
92+
93+ sceneRef.current = { scene, group: null, renderer, camera };
94+
95+ let raf = 0;
96+ const animate = () => {
97+ raf = requestAnimationFrame(animate);
98+ controls.update();
99+ renderer.render(scene, camera);
100+ };
101+ animate();
102+
103+ const resize = () => {
104+ const w = mount.clientWidth;
105+ const h = mount.clientHeight;
106+ renderer.setSize(w, h, false);
107+ camera.aspect = w / h;
108+ camera.updateProjectionMatrix();
109+ };
110+ resize();
111+ const ro = new ResizeObserver(resize);
112+ ro.observe(mount);
113+
114+ return () => {
115+ cancelAnimationFrame(raf);
116+ ro.disconnect();
117+ controls.dispose();
118+ renderer.dispose();
119+ mount.removeChild(renderer.domElement);
120+ sceneRef.current = null;
121+ };
122+ }, []);
123+
124+ // rebuild the geometry group whenever relevant parameters or impacts change
125+ useEffect(() => {
126+ const ctx = sceneRef.current;
127+ if (!ctx) return;
128+ if (ctx.group) {
129+ ctx.scene.remove(ctx.group);
130+ ctx.group.traverse((o) => {
131+ const mesh = o as THREE.Mesh;
132+ if (mesh.geometry) mesh.geometry.dispose();
133+ });
134+ }
135+
136+ const g = new THREE.Group();
137+ const R = params.shell_radius_m;
138+
139+ // shell: translucent so the impact points on the far side stay visible
140+ const shell = new THREE.Mesh(
141+ new THREE.SphereGeometry(R, 48, 32),
142+ new THREE.MeshStandardMaterial({
143+ color: 0xc08b52,
144+ roughness: 0.7,
145+ metalness: 0.0,
146+ transparent: true,
147+ opacity: 0.42,
148+ depthWrite: false,
149+ }),
150+ );
151+ g.add(shell);
152+ g.add(
153+ new THREE.LineSegments(
154+ new THREE.WireframeGeometry(new THREE.SphereGeometry(R, 24, 16)),
155+ new THREE.LineBasicMaterial({ color: 0x6d4a29, transparent: true, opacity: 0.35 }),
156+ ),
157+ );
158+
159+ // handle, along -z (the shell is shaken about the x axis)
160+ const handle = new THREE.Mesh(
161+ new THREE.CylinderGeometry(R * 0.16, R * 0.13, R * 2.6, 16),
162+ new THREE.MeshStandardMaterial({ color: 0x8d6e63, roughness: 0.75 }),
163+ );
164+ handle.position.set(...m2t(0, 0, -R - R * 1.2));
165+ g.add(handle);
166+
167+ // grain impact points
168+ if (impactPos && impactV && impactT && impactV.length > 0) {
169+ const n = impactV.length;
170+ let vMax = 1e-9;
171+ for (let i = 0; i < n; i++) vMax = Math.max(vMax, impactV[i]);
172+ const period = 1 / Math.max(1e-6, params.stroke_rate_hz);
173+ const pos = new Float32Array(n * 3);
174+ const col = new Float32Array(n * 3);
175+ const sizes = new Float32Array(n);
176+ const warm = new THREE.Color(0xffb74d);
177+ const cool = new THREE.Color(0x4dd0e1);
178+ for (let i = 0; i < n; i++) {
179+ // sit the marker just outside the surface so it is not z-fought
180+ const s = R * 1.012;
181+ const t = m2t(impactPos[i * 3] * s, impactPos[i * 3 + 1] * s, impactPos[i * 3 + 2] * s);
182+ pos[i * 3] = t[0];
183+ pos[i * 3 + 1] = t[1];
184+ pos[i * 3 + 2] = t[2];
185+ const stroke = Math.floor(impactT[i] / period);
186+ const c = stroke % 2 === 0 ? warm : cool;
187+ col[i * 3] = c.r;
188+ col[i * 3 + 1] = c.g;
189+ col[i * 3 + 2] = c.b;
190+ sizes[i] = impactV[i] / vMax;
191+ }
192+ const geo = new THREE.BufferGeometry();
193+ geo.setAttribute('position', new THREE.Float32BufferAttribute(pos, 3));
194+ geo.setAttribute('color', new THREE.Float32BufferAttribute(col, 3));
195+ g.add(
196+ new THREE.Points(
197+ geo,
198+ new THREE.PointsMaterial({
199+ size: R * 0.055,
200+ vertexColors: true,
201+ transparent: true,
202+ opacity: 0.75,
203+ sizeAttenuation: true,
204+ }),
205+ ),
206+ );
207+ }
208+
209+ // the two poles the grains alternate between
210+ for (const sgn of [1, -1]) {
211+ const dot = new THREE.Mesh(
212+ new THREE.SphereGeometry(R * 0.07, 14, 10),
213+ new THREE.MeshStandardMaterial({
214+ color: sgn > 0 ? 0xffb74d : 0x4dd0e1,
215+ emissive: sgn > 0 ? 0x4a2f0a : 0x0c3b42,
216+ }),
217+ );
218+ dot.position.set(...m2t(sgn * R * 1.04, 0, 0));
219+ g.add(dot);
220+ }
221+ const poleLabel = makeLabel('impact poles', '#ffb74d');
222+ poleLabel.position.set(...m2t(R * 1.1, 0, R * 0.55));
223+ g.add(poleLabel);
224+
225+ // listener: microphone marker + line from the shell centre
226+ const lis = new THREE.Vector3(
227+ ...m2t(params.listener_x_m, params.listener_y_m, params.listener_z_m),
228+ );
229+ const micHead = new THREE.Mesh(
230+ new THREE.SphereGeometry(0.016, 20, 16),
231+ new THREE.MeshStandardMaterial({ color: 0x4dd0e1, emissive: 0x0c3b42, roughness: 0.35 }),
232+ );
233+ micHead.position.copy(lis);
234+ g.add(micHead);
235+ const micBody = new THREE.Mesh(
236+ new THREE.CylinderGeometry(0.008, 0.008, 0.05, 12),
237+ new THREE.MeshStandardMaterial({ color: 0x455a64, roughness: 0.6 }),
238+ );
239+ micBody.position.copy(lis).add(new THREE.Vector3(0, -0.033, 0));
240+ g.add(micBody);
241+ const lineGeo = new THREE.BufferGeometry().setFromPoints([new THREE.Vector3(0, 0, 0), lis]);
242+ const line = new THREE.Line(
243+ lineGeo,
244+ new THREE.LineDashedMaterial({ color: 0x4dd0e1, dashSize: 0.012, gapSize: 0.009 }),
245+ );
246+ line.computeLineDistances();
247+ g.add(line);
248+ const dist = Math.hypot(params.listener_x_m, params.listener_y_m, params.listener_z_m);
249+ const micLabel = makeLabel(`mic ${(dist * 100).toFixed(0)} cm`, '#4dd0e1');
250+ micLabel.position.copy(lis).add(new THREE.Vector3(0, 0.045, 0));
251+ g.add(micLabel);
252+
253+ ctx.scene.add(g);
254+ ctx.group = g;
255+ }, [
256+ params.shell_radius_m,
257+ params.stroke_rate_hz,
258+ params.listener_x_m,
259+ params.listener_y_m,
260+ params.listener_z_m,
261+ impactPos,
262+ impactV,
263+ impactT,
264+ ]);
265+
266+ return (
267+ <div className="scene3d" ref={mountRef}>
268+ <div className="scene-hint">drag to orbit · scroll to zoom · right-drag to pan</div>
269+ </div>
270+ );
271+}
webapp/src/components/plots/ImpactRaster.tsxadded+96−0View file
@@ -0,0 +1,96 @@
1+// The stochastic impact process itself: every grain-shell collision as a dot
2+// at its time and speed, over the expected-rate envelope the strokes impose.
3+// The gap between bursts is where the maraca's rhythm lives, and a burst
4+// width comparable to the stroke period smears it into a continuous wash.
5+
6+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
7+
8+interface Props {
9+ impactT: Float32Array;
10+ impactV: Float32Array;
11+ rateT: Float32Array;
12+ rateHz: Float32Array;
13+ durationS: number;
14+ impactCount: number;
15+ plotted: number;
16+}
17+
18+export default function ImpactRaster({
19+ impactT,
20+ impactV,
21+ rateT,
22+ rateHz,
23+ durationS,
24+ impactCount,
25+ plotted,
26+}: Props) {
27+ const ref = useCanvas(
28+ (ctx, w, h) => {
29+ drawFrame(ctx, w, h);
30+ const padL = 40;
31+ const padB = 16;
32+ const pw = w - padL - 8;
33+ const ph = h - padB - 12;
34+ const top = 10;
35+
36+ let vMax = 0;
37+ for (let i = 0; i < impactV.length; i++) vMax = Math.max(vMax, impactV[i]);
38+ if (vMax <= 0) vMax = 1;
39+ let rMax = 0;
40+ for (let i = 0; i < rateHz.length; i++) rMax = Math.max(rMax, rateHz[i]);
41+
42+ const tx = (t: number) => padL + (t / durationS) * pw;
43+
44+ // time grid
45+ ctx.font = '10px system-ui';
46+ ctx.textAlign = 'center';
47+ const tStep = durationS > 4 ? 1.0 : durationS > 1.5 ? 0.5 : 0.1;
48+ for (let t = 0; t <= durationS + 1e-9; t += tStep) {
49+ const x = tx(t);
50+ ctx.strokeStyle = PLOT_COLORS.grid;
51+ ctx.beginPath();
52+ ctx.moveTo(x, top);
53+ ctx.lineTo(x, top + ph);
54+ ctx.stroke();
55+ ctx.fillStyle = PLOT_COLORS.axis;
56+ ctx.fillText(`${t.toFixed(tStep < 0.5 ? 1 : 1)}s`, x, h - 4);
57+ }
58+ ctx.textAlign = 'right';
59+ ctx.fillStyle = PLOT_COLORS.axis;
60+ ctx.fillText(`${vMax.toFixed(1)} m/s`, padL - 4, top + 4);
61+ ctx.fillText('0', padL - 4, top + ph);
62+
63+ // expected-rate envelope, filled under the curve
64+ if (rMax > 0) {
65+ ctx.fillStyle = 'rgba(111, 191, 115, 0.16)';
66+ ctx.beginPath();
67+ ctx.moveTo(padL, top + ph);
68+ for (let i = 0; i < rateT.length; i++) {
69+ ctx.lineTo(tx(rateT[i]), top + ph - (rateHz[i] / rMax) * ph * 0.9);
70+ }
71+ ctx.lineTo(padL + pw, top + ph);
72+ ctx.closePath();
73+ ctx.fill();
74+ }
75+
76+ // impacts
77+ ctx.fillStyle = 'rgba(77, 208, 225, 0.55)';
78+ for (let i = 0; i < impactT.length; i++) {
79+ const x = tx(impactT[i]);
80+ const y = top + ph - (impactV[i] / vMax) * ph;
81+ ctx.fillRect(x, y, 1.4, 1.4);
82+ }
83+
84+ ctx.fillStyle = PLOT_COLORS.axis;
85+ ctx.textAlign = 'left';
86+ const note =
87+ plotted < impactCount
88+ ? `${impactCount} impacts (${plotted} plotted) · peak rate ${rMax.toFixed(0)}/s`
89+ : `${impactCount} impacts · peak rate ${rMax.toFixed(0)}/s`;
90+ ctx.fillText(note, padL + 6, top + 4);
91+ },
92+ [impactT, impactV, rateT, rateHz, durationS, impactCount, plotted],
93+ 190,
94+ );
95+ return <canvas ref={ref} />;
96+}
webapp/src/components/plots/ModeMap.tsxadded+107−0View file
@@ -0,0 +1,107 @@
1+// The retained shell modes, plotted as frequency against spherical-harmonic
2+// degree l. The two membrane branches separate clearly: the lower one climbs
3+// slowly from the l = 2 mode and is what the grains mostly excite, the upper
4+// one starts at the breathing mode. Marker area is the radiation efficiency
5+// sigma_l, and the bar under each marker is its damping ratio, so a mode that
6+// is loud but short-lived is visible as a large marker with a long bar.
7+
8+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
9+
10+interface Props {
11+ modeF: Float32Array;
12+ modeL: Float32Array;
13+ modeZeta: Float32Array;
14+ modeSigma: Float32Array;
15+ breathingHz: number;
16+}
17+
18+export default function ModeMap({ modeF, modeL, modeZeta, modeSigma, breathingHz }: Props) {
19+ const ref = useCanvas(
20+ (ctx, w, h) => {
21+ drawFrame(ctx, w, h);
22+ const padL = 40;
23+ const padB = 20;
24+ const pw = w - padL - 10;
25+ const ph = h - padB - 14;
26+ const top = 12;
27+ const n = modeF.length;
28+ if (n === 0) return;
29+
30+ let lMax = 1;
31+ let fMax = 1;
32+ let zMax = 1e-9;
33+ let sMax = 1e-9;
34+ for (let i = 0; i < n; i++) {
35+ lMax = Math.max(lMax, modeL[i]);
36+ fMax = Math.max(fMax, modeF[i]);
37+ zMax = Math.max(zMax, modeZeta[i]);
38+ sMax = Math.max(sMax, modeSigma[i]);
39+ }
40+ fMax *= 1.08;
41+
42+ const lx = (l: number) => padL + (l / lMax) * pw;
43+ const fy = (f: number) => top + ph - (f / fMax) * ph;
44+
45+ ctx.font = '10px system-ui';
46+ ctx.textAlign = 'center';
47+ const lStep = lMax > 20 ? 5 : lMax > 8 ? 2 : 1;
48+ for (let l = 0; l <= lMax; l += lStep) {
49+ const x = lx(l);
50+ ctx.strokeStyle = PLOT_COLORS.grid;
51+ ctx.beginPath();
52+ ctx.moveTo(x, top);
53+ ctx.lineTo(x, top + ph);
54+ ctx.stroke();
55+ ctx.fillStyle = PLOT_COLORS.axis;
56+ ctx.fillText(`${l}`, x, h - 7);
57+ }
58+ ctx.fillText('degree l', padL + pw / 2, h - 0);
59+ ctx.textAlign = 'right';
60+ const fStep = fMax > 16000 ? 5000 : fMax > 6000 ? 2000 : 1000;
61+ for (let f = 0; f <= fMax; f += fStep) {
62+ const y = fy(f);
63+ ctx.strokeStyle = PLOT_COLORS.grid;
64+ ctx.beginPath();
65+ ctx.moveTo(padL, y);
66+ ctx.lineTo(padL + pw, y);
67+ ctx.stroke();
68+ ctx.fillStyle = PLOT_COLORS.axis;
69+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, padL - 4, y + 3);
70+ }
71+
72+ // damping bars, then markers on top
73+ for (let i = 0; i < n; i++) {
74+ const x = lx(modeL[i]);
75+ const y = fy(modeF[i]);
76+ const barH = (modeZeta[i] / zMax) * 22;
77+ ctx.strokeStyle = 'rgba(232, 163, 61, 0.5)';
78+ ctx.beginPath();
79+ ctx.moveTo(x, y);
80+ ctx.lineTo(x, y + barH);
81+ ctx.stroke();
82+ }
83+ for (let i = 0; i < n; i++) {
84+ const x = lx(modeL[i]);
85+ const y = fy(modeF[i]);
86+ const r = 2.5 + 5 * Math.sqrt(modeSigma[i] / sMax);
87+ const isBreathing = Math.abs(modeF[i] - breathingHz) < 1e-6;
88+ ctx.fillStyle = isBreathing ? PLOT_COLORS.good : PLOT_COLORS.trace;
89+ ctx.beginPath();
90+ ctx.arc(x, y, r, 0, 2 * Math.PI);
91+ ctx.fill();
92+ if (isBreathing) {
93+ ctx.fillStyle = PLOT_COLORS.good;
94+ ctx.textAlign = 'left';
95+ ctx.fillText(`breathing ${breathingHz.toFixed(0)} Hz`, x + r + 4, y + 3);
96+ }
97+ }
98+
99+ ctx.fillStyle = PLOT_COLORS.axis;
100+ ctx.textAlign = 'left';
101+ ctx.fillText(`${n} modes · marker = radiation efficiency · bar = damping`, padL + 6, top + 2);
102+ },
103+ [modeF, modeL, modeZeta, modeSigma, breathingHz],
104+ 210,
105+ );
106+ return <canvas ref={ref} />;
107+}
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+116−0View file
@@ -0,0 +1,116 @@
1+// Magnitude spectrum of the rendered shake, log-frequency axis, with the
2+// retained shell modes and the breathing mode overlaid. Unlike a struck bar,
3+// the spectrum here is broadband: thousands of short Hertz contacts fill the
4+// space between the modes, and the modal peaks sit on top of that floor.
5+
6+import { magnitudeSpectrum } from '../../lib/fft';
7+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
8+
9+interface Props {
10+ y: Float32Array;
11+ sampleRate: number;
12+ modeF: Float32Array;
13+ breathingHz: number;
14+}
15+
16+const F_LO = 100;
17+const F_HI = 24000;
18+const DB_RANGE = 80;
19+
20+export default function SpectrumPlot({ y, sampleRate, modeF, breathingHz }: Props) {
21+ const ref = useCanvas(
22+ (ctx, w, h) => {
23+ drawFrame(ctx, w, h);
24+ const padL = 34;
25+ const padB = 18;
26+ const pw = w - padL - 6;
27+ const ph = h - padB - 8;
28+ const fx = (f: number) =>
29+ padL + ((Math.log10(f) - Math.log10(F_LO)) / (Math.log10(F_HI) - Math.log10(F_LO))) * pw;
30+
31+ ctx.font = '10px system-ui';
32+ ctx.textAlign = 'center';
33+ for (const f of [100, 200, 500, 1000, 2000, 5000, 10000, 20000]) {
34+ const x = fx(f);
35+ ctx.strokeStyle = PLOT_COLORS.grid;
36+ ctx.beginPath();
37+ ctx.moveTo(x, 8);
38+ ctx.lineTo(x, 8 + ph);
39+ ctx.stroke();
40+ ctx.fillStyle = PLOT_COLORS.axis;
41+ ctx.fillText(f >= 1000 ? `${f / 1000}k` : `${f}`, x, h - 5);
42+ }
43+ ctx.textAlign = 'right';
44+ for (let db = 0; db >= -DB_RANGE; db -= 20) {
45+ const yy = 8 + (-db / DB_RANGE) * ph;
46+ ctx.strokeStyle = PLOT_COLORS.grid;
47+ ctx.beginPath();
48+ ctx.moveTo(padL, yy);
49+ ctx.lineTo(padL + pw, yy);
50+ ctx.stroke();
51+ ctx.fillStyle = PLOT_COLORS.axis;
52+ ctx.fillText(`${db}`, padL - 4, yy + 3);
53+ }
54+
55+ const { mag, size } = magnitudeSpectrum(y, 32768);
56+ let peak = 0;
57+ for (let i = 1; i < mag.length; i++) peak = Math.max(peak, mag[i]);
58+ if (peak <= 0) return;
59+
60+ // overlays first so the trace draws on top
61+ ctx.strokeStyle = PLOT_COLORS.accent;
62+ ctx.setLineDash([4, 3]);
63+ for (let i = 0; i < modeF.length; i++) {
64+ const f = modeF[i];
65+ if (f < F_LO || f > F_HI) continue;
66+ ctx.beginPath();
67+ ctx.moveTo(fx(f), 8);
68+ ctx.lineTo(fx(f), 8 + ph);
69+ ctx.stroke();
70+ }
71+ if (breathingHz >= F_LO && breathingHz <= F_HI) {
72+ ctx.strokeStyle = PLOT_COLORS.good;
73+ ctx.beginPath();
74+ ctx.moveTo(fx(breathingHz), 8);
75+ ctx.lineTo(fx(breathingHz), 8 + ph);
76+ ctx.stroke();
77+ ctx.fillStyle = PLOT_COLORS.good;
78+ ctx.textAlign = 'left';
79+ ctx.fillText(`breathing ${breathingHz.toFixed(0)}`, fx(breathingHz) + 3, 30);
80+ }
81+ ctx.setLineDash([]);
82+
83+ ctx.strokeStyle = PLOT_COLORS.trace;
84+ ctx.beginPath();
85+ let started = false;
86+ for (let px = 0; px <= pw; px++) {
87+ const f = 10 ** (Math.log10(F_LO) + (px / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
88+ const bin = Math.min(mag.length - 1, Math.round((f / sampleRate) * size));
89+ const fNext =
90+ 10 ** (Math.log10(F_LO) + ((px + 1) / pw) * (Math.log10(F_HI) - Math.log10(F_LO)));
91+ const binNext = Math.min(
92+ mag.length - 1,
93+ Math.max(bin + 1, Math.round((fNext / sampleRate) * size)),
94+ );
95+ let m = 0;
96+ for (let b = bin; b < binNext; b++) m = Math.max(m, mag[b]);
97+ const db = 20 * Math.log10(m / peak + 1e-12);
98+ const yy = 8 + Math.min(ph, (-db / DB_RANGE) * ph);
99+ if (!started) {
100+ ctx.moveTo(padL + px, yy);
101+ started = true;
102+ } else {
103+ ctx.lineTo(padL + px, yy);
104+ }
105+ }
106+ ctx.stroke();
107+
108+ ctx.fillStyle = PLOT_COLORS.accent;
109+ ctx.textAlign = 'left';
110+ ctx.fillText('shell modes', padL + 4, 18);
111+ },
112+ [y, sampleRate, modeF, breathingHz],
113+ 200,
114+ );
115+ return <canvas ref={ref} />;
116+}
webapp/src/components/plots/WaveformPlot.tsxadded+64−0View file
@@ -0,0 +1,64 @@
1+import { useCanvas, drawFrame, PLOT_COLORS } from './canvasUtils';
2+
3+interface Props {
4+ y: Float32Array;
5+ sampleRate: number;
6+}
7+
8+export default function WaveformPlot({ y, sampleRate }: Props) {
9+ const ref = useCanvas(
10+ (ctx, w, h) => {
11+ drawFrame(ctx, w, h);
12+ const padL = 34;
13+ const padB = 16;
14+ const pw = w - padL - 6;
15+ const ph = h - padB - 6;
16+ const midY = 6 + ph / 2;
17+
18+ ctx.strokeStyle = PLOT_COLORS.grid;
19+ ctx.beginPath();
20+ ctx.moveTo(padL, midY);
21+ ctx.lineTo(padL + pw, midY);
22+ ctx.stroke();
23+
24+ // time gridlines every 0.1 s
25+ const durS = y.length / sampleRate;
26+ ctx.fillStyle = PLOT_COLORS.axis;
27+ ctx.font = '10px system-ui';
28+ ctx.textAlign = 'center';
29+ for (let t = 0; t <= durS + 1e-9; t += 0.1) {
30+ const x = padL + (t / durS) * pw;
31+ ctx.strokeStyle = PLOT_COLORS.grid;
32+ ctx.beginPath();
33+ ctx.moveTo(x, 6);
34+ ctx.lineTo(x, 6 + ph);
35+ ctx.stroke();
36+ ctx.fillText(`${t.toFixed(1)}s`, x, h - 4);
37+ }
38+ ctx.textAlign = 'right';
39+ ctx.fillText('+1', padL - 4, 12);
40+ ctx.fillText('-1', padL - 4, 6 + ph);
41+
42+ // min/max envelope per pixel column
43+ ctx.strokeStyle = PLOT_COLORS.trace;
44+ ctx.beginPath();
45+ const perPx = y.length / pw;
46+ for (let px = 0; px < pw; px++) {
47+ const i0 = Math.floor(px * perPx);
48+ const i1 = Math.min(y.length, Math.max(i0 + 1, Math.floor((px + 1) * perPx)));
49+ let mn = Infinity, mx = -Infinity;
50+ for (let i = i0; i < i1; i++) {
51+ if (y[i] < mn) mn = y[i];
52+ if (y[i] > mx) mx = y[i];
53+ }
54+ if (!Number.isFinite(mn)) continue;
55+ ctx.moveTo(padL + px + 0.5, midY - mx * (ph / 2));
56+ ctx.lineTo(padL + px + 0.5, midY - mn * (ph / 2) + 0.5);
57+ }
58+ ctx.stroke();
59+ },
60+ [y, sampleRate],
61+ 170,
62+ );
63+ return <canvas ref={ref} />;
64+}
webapp/src/components/plots/canvasUtils.tsadded+66−0View file
@@ -0,0 +1,66 @@
1+import { useEffect, useRef } from 'react';
2+
3+// Canvas hook: resizes to element width at devicePixelRatio and calls draw.
4+export function useCanvas(
5+ draw: (ctx: CanvasRenderingContext2D, w: number, h: number) => void,
6+ deps: unknown[],
7+ height: number,
8+) {
9+ const ref = useRef<HTMLCanvasElement | null>(null);
10+
11+ useEffect(() => {
12+ const canvas = ref.current;
13+ if (!canvas) return;
14+ const render = () => {
15+ const dpr = window.devicePixelRatio || 1;
16+ const w = canvas.clientWidth;
17+ canvas.width = Math.max(1, Math.round(w * dpr));
18+ canvas.height = Math.max(1, Math.round(height * dpr));
19+ canvas.style.height = `${height}px`;
20+ const ctx = canvas.getContext('2d')!;
21+ ctx.setTransform(dpr, 0, 0, dpr, 0, 0);
22+ ctx.clearRect(0, 0, w, height);
23+ draw(ctx, w, height);
24+ };
25+ render();
26+ const ro = new ResizeObserver(render);
27+ ro.observe(canvas);
28+ return () => ro.disconnect();
29+ // eslint-disable-next-line react-hooks/exhaustive-deps
30+ }, deps);
31+
32+ return ref;
33+}
34+
35+export const PLOT_COLORS = {
36+ bg: '#181b20',
37+ grid: '#2a2f38',
38+ axis: '#8b93a1',
39+ trace: '#4dd0e1',
40+ accent: '#e8a33d',
41+ good: '#6fbf73',
42+ warn: '#e57373',
43+};
44+
45+export function drawFrame(ctx: CanvasRenderingContext2D, w: number, h: number) {
46+ ctx.fillStyle = PLOT_COLORS.bg;
47+ ctx.fillRect(0, 0, w, h);
48+}
49+
50+// Viridis-like colormap for the spectrogram, t in [0,1].
51+export function colormap(t: number): [number, number, number] {
52+ const stops: [number, number, number][] = [
53+ [13, 8, 60], [84, 39, 143], [140, 41, 129], [196, 66, 90],
54+ [237, 105, 37], [251, 155, 6], [247, 209, 61], [252, 255, 164],
55+ ];
56+ const x = Math.min(0.9999, Math.max(0, t)) * (stops.length - 1);
57+ const i = Math.floor(x);
58+ const f = x - i;
59+ const a = stops[i];
60+ const b = stops[i + 1];
61+ return [
62+ Math.round(a[0] + f * (b[0] - a[0])),
63+ Math.round(a[1] + f * (b[1] - a[1])),
64+ Math.round(a[2] + f * (b[2] - a[2])),
65+ ];
66+}
webapp/src/hooks/useRenderer.tsadded+114−0View file
@@ -0,0 +1,114 @@
1+// Hook owning the render worker: debounced auto-preview whenever the
2+// parameters or the .m sources change, explicit full-quality render, and a
3+// latest-wins request queue.
4+
5+import { useCallback, useEffect, useRef, useState } from 'react';
6+import type { MaracaParams } 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: MaracaParams;
26+ mfiles: MFile[];
27+ quality: 'preview' | 'full';
28+}
29+
30+export function useRenderer(params: MaracaParams, 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/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+43−0View file
@@ -0,0 +1,43 @@
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+ 'maraca_basis.m',
25+ 'maraca_hankel.m',
26+ 'maraca_impacts.m',
27+ 'maraca_radiate.m',
28+ 'maraca_rng.m',
29+ 'maraca_output.m',
30+]);
31+
32+// Cheap content hash (FNV-1a) so the worker can key its caches on the code.
33+export function hashMFiles(files: MFile[]): string {
34+ let h = 0x811c9dc5;
35+ for (const f of files) {
36+ const s = f.name + ' ' + f.source + '';
37+ for (let i = 0; i < s.length; i++) {
38+ h ^= s.charCodeAt(i);
39+ h = Math.imul(h, 0x01000193);
40+ }
41+ }
42+ return (h >>> 0).toString(36);
43+}
webapp/src/params.tsadded+260−0View file
@@ -0,0 +1,260 @@
1+// Parameter definitions for the maraca model, mirroring maraca_defaults.m.
2+// Each parameter is assigned a tier that determines how much recomputation a
3+// change requires:
4+// solve -> mode basis + impact process + radiation + output
5+// radiate -> radiation + output only (cached basis and impacts are reused)
6+// output -> output stage only (decimation product is reused)
7+//
8+// Note that the render duration is not a parameter: it is
9+// stroke_count / stroke_rate_hz + tail_s.
10+
11+export interface MaracaParams {
12+ // Shell geometry
13+ shell_radius_m: number;
14+ shell_thickness_m: number;
15+ // Shell material
16+ density: number;
17+ youngs_modulus: number;
18+ poisson_ratio: number;
19+ stiffness_trim: number;
20+ max_mode_l: number;
21+ max_mode_hz: number;
22+ // Damping
23+ loss_factor: number;
24+ loss_slope: number;
25+ radiation_loss: number;
26+ // Grains
27+ grain_radius_m: number;
28+ grain_density: number;
29+ grain_youngs: number;
30+ grain_size_spread: number;
31+ restitution: number;
32+ // Shake gesture / impact process
33+ stroke_rate_hz: number;
34+ stroke_count: number;
35+ accent_period: number;
36+ accent_gain: number;
37+ grain_count: number;
38+ bounces_per_grain: number;
39+ shake_velocity: number;
40+ slosh_lag_ms: number;
41+ burst_width_ms: number;
42+ impact_spread_deg: number;
43+ seed: number;
44+ // Air cavity
45+ cavity_coupling: number;
46+ // Radiation / listener
47+ direct_gain: number;
48+ listener_x_m: number;
49+ listener_y_m: number;
50+ listener_z_m: number;
51+ speed_of_sound: number;
52+ air_density: number;
53+ // Rendering
54+ sample_rate: number;
55+ oversample: number;
56+ tail_s: number;
57+ pre_delay_ms: number;
58+ normalize: boolean;
59+ peak_dbfs: number;
60+ full_scale_pa: number;
61+ bit_depth: number;
62+}
63+
64+export const DEFAULTS: MaracaParams = {
65+ shell_radius_m: 0.034,
66+ shell_thickness_m: 0.0022,
67+ density: 720.0,
68+ youngs_modulus: 2.5e9,
69+ poisson_ratio: 0.3,
70+ stiffness_trim: 1.0,
71+ max_mode_l: 16,
72+ max_mode_hz: 22000.0,
73+ loss_factor: 0.015,
74+ loss_slope: 6.0e-7,
75+ radiation_loss: 1.0,
76+ grain_radius_m: 0.0015,
77+ grain_density: 1300.0,
78+ grain_youngs: 6.0e8,
79+ grain_size_spread: 0.25,
80+ restitution: 0.55,
81+ stroke_rate_hz: 4.5,
82+ stroke_count: 8,
83+ accent_period: 2,
84+ accent_gain: 1.5,
85+ grain_count: 320,
86+ bounces_per_grain: 1.3,
87+ shake_velocity: 2.2,
88+ slosh_lag_ms: 14.0,
89+ burst_width_ms: 45.0,
90+ impact_spread_deg: 35.0,
91+ seed: 1,
92+ cavity_coupling: 1.0,
93+ direct_gain: 1.0,
94+ listener_x_m: 0.15,
95+ listener_y_m: 0.25,
96+ listener_z_m: 0.3,
97+ speed_of_sound: 343.0,
98+ air_density: 1.204,
99+ sample_rate: 48000,
100+ oversample: 2,
101+ tail_s: 0.4,
102+ pre_delay_ms: 3.0,
103+ normalize: true,
104+ peak_dbfs: -3.0,
105+ full_scale_pa: 20.0,
106+ bit_depth: 24,
107+};
108+
109+export type Tier = 'solve' | 'radiate' | 'output';
110+
111+// Keys whose change invalidates the cached mode basis or impact process. The
112+// listener position and direct_gain enter only the radiation stage, so they
113+// are deliberately absent.
114+export const SOLVE_KEYS: (keyof MaracaParams)[] = [
115+ 'shell_radius_m', 'shell_thickness_m',
116+ 'density', 'youngs_modulus', 'poisson_ratio', 'stiffness_trim',
117+ 'max_mode_l', 'max_mode_hz',
118+ 'loss_factor', 'loss_slope', 'radiation_loss',
119+ 'grain_radius_m', 'grain_density', 'grain_youngs', 'grain_size_spread',
120+ 'restitution',
121+ 'stroke_rate_hz', 'stroke_count', 'accent_period', 'accent_gain',
122+ 'grain_count', 'bounces_per_grain', 'shake_velocity', 'slosh_lag_ms', 'burst_width_ms',
123+ 'impact_spread_deg', 'seed',
124+ 'cavity_coupling',
125+ 'speed_of_sound', 'air_density', 'sample_rate', 'tail_s',
126+];
127+
128+export const RADIATE_KEYS: (keyof MaracaParams)[] = [
129+ 'direct_gain', 'listener_x_m', 'listener_y_m', 'listener_z_m',
130+];
131+
132+export const OUTPUT_KEYS: (keyof MaracaParams)[] = [
133+ 'pre_delay_ms', 'normalize', 'peak_dbfs', 'full_scale_pa',
134+];
135+
136+export interface ParamSpec {
137+ key: keyof MaracaParams;
138+ label: string;
139+ unit?: string;
140+ min: number;
141+ max: number;
142+ step?: number;
143+ log?: boolean;
144+ integer?: boolean;
145+ help?: string;
146+}
147+
148+export interface ParamGroup {
149+ name: string;
150+ specs: ParamSpec[];
151+}
152+
153+export const PARAM_GROUPS: ParamGroup[] = [
154+ {
155+ name: 'Shake gesture',
156+ specs: [
157+ { key: 'stroke_rate_hz', label: 'Stroke rate', unit: 'Hz', min: 0.5, max: 12, step: 0.1, help: 'Strokes per second. Sets the rhythm and the render length.' },
158+ { key: 'stroke_count', label: 'Strokes', min: 1, max: 32, step: 1, integer: true, help: 'Number of strokes in the pattern.' },
159+ { key: 'accent_period', label: 'Accent every', unit: 'strokes', min: 1, max: 8, step: 1, integer: true, help: '1 accents every stroke, 2 alternates, 4 gives a bar-like pulse.' },
160+ { key: 'accent_gain', label: 'Accent strength', min: 1.0, max: 3.0, step: 0.05, help: 'Velocity and count multiplier on accented strokes.' },
161+ { key: 'shake_velocity', label: 'Shake velocity', unit: 'm/s', min: 0.2, max: 6.0, step: 0.05, help: 'Mean grain impact speed. Dynamics and timbre together: faster grains give shorter Hertz contacts and a brighter hiss.' },
162+ { key: 'slosh_lag_ms', label: 'Slosh lag', unit: 'ms', min: 0, max: 80, step: 0.5, help: 'Delay from the hand reversing to the grain mass arriving at the far wall.' },
163+ { key: 'burst_width_ms', label: 'Burst width', unit: 'ms', min: 3, max: 200, step: 1, help: 'Spread of one stroke’s impacts in time. Small = a tight chick, large = a wash.' },
164+ { key: 'seed', label: 'Random seed', min: 1, max: 9999, step: 1, integer: true, help: 'Any change redraws the whole impact process; the same seed always reproduces it.' },
165+ ],
166+ },
167+ {
168+ name: 'Grains',
169+ specs: [
170+ { key: 'grain_count', label: 'Grain count', min: 5, max: 3000, step: 1, integer: true, help: 'How many seeds or beads are inside the shell. With bounces per grain, this sets the expected impacts per stroke.' },
171+ { key: 'bounces_per_grain', label: 'Bounces per grain', unit: '/stroke', min: 0.2, max: 6.0, step: 0.05, help: 'Mean wall contacts each grain makes in one stroke: one arrival plus whatever rattle follows before the hand reverses.' },
172+ { key: 'grain_radius_m', label: 'Grain radius', unit: 'm', min: 0.0003, max: 0.006, step: 0.0001, help: 'Seed or bead size. Bigger grains give heavier, slower, lower-pitched clicks.' },
173+ { key: 'grain_density', label: 'Grain density', unit: 'kg/m³', min: 300, max: 8000, log: true, help: 'Seed ~1300, glass ~2500, steel shot ~7800.' },
174+ { key: 'grain_youngs', label: 'Grain modulus', unit: 'Pa', min: 1e7, max: 2e11, log: true, help: 'Contact hardness; sets the Hertz contact time and therefore the brightness.' },
175+ { key: 'grain_size_spread', label: 'Size spread', min: 0.0, max: 0.8, step: 0.01, help: 'Lognormal spread of the grain radius. 0 = identical grains, which sounds mechanical.' },
176+ { key: 'restitution', label: 'Restitution', min: 0.05, max: 0.95, step: 0.01, help: 'Coefficient of restitution; scales the momentum each impact transfers.' },
177+ { key: 'impact_spread_deg', label: 'Impact spread', unit: '°', min: 5, max: 90, step: 1, help: 'Angular scatter of the impact points about the pole the grains are thrown at.' },
178+ ],
179+ },
180+ {
181+ name: 'Shell geometry',
182+ specs: [
183+ { key: 'shell_radius_m', label: 'Shell radius', unit: 'm', min: 0.015, max: 0.09, step: 0.001, help: 'Mode frequencies scale as 1/R.' },
184+ { key: 'shell_thickness_m', label: 'Shell thickness', unit: 'm', min: 0.0005, max: 0.008, step: 0.0001, help: 'Sets the bending term that spreads the high-degree modes, and the wall mass.' },
185+ ],
186+ },
187+ {
188+ name: 'Shell material',
189+ specs: [
190+ { key: 'density', label: 'Density', unit: 'kg/m³', min: 200, max: 2500, step: 10, help: 'Dried gourd ~700. Frequencies scale as 1/sqrt(rho).' },
191+ { key: 'youngs_modulus', label: "Young's modulus", unit: 'Pa', min: 2e8, max: 2e11, log: true, help: 'Frequencies scale as sqrt(E).' },
192+ { key: 'poisson_ratio', label: 'Poisson ratio', min: 0.0, max: 0.48, step: 0.01 },
193+ { key: 'stiffness_trim', label: 'Stiffness trim', min: 0.4, max: 2.0, step: 0.01, help: 'Overall pitch trim for the departures of a real gourd from a perfect sphere.' },
194+ { key: 'cavity_coupling', label: 'Cavity coupling', min: 0.0, max: 4.0, step: 0.05, help: 'Air stiffness inside the closed shell, which raises the breathing mode only.' },
195+ ],
196+ },
197+ {
198+ name: 'Damping',
199+ specs: [
200+ { key: 'loss_factor', label: 'Loss factor η₀', min: 0.001, max: 0.12, step: 0.001, help: 'Frequency-independent wall loss; sets how long each impact rings.' },
201+ { key: 'loss_slope', label: 'Loss slope η₁', unit: 's', min: 1e-8, max: 4e-6, log: true, help: 'Extra high-frequency damping; the main brightness control.' },
202+ { key: 'radiation_loss', label: 'Radiation loss', min: 0.0, max: 3.0, step: 0.05, help: 'Scales the acoustic radiation damping computed from the multipole impedance.' },
203+ ],
204+ },
205+ {
206+ name: 'Listener',
207+ specs: [
208+ { key: 'listener_x_m', label: 'Listener x', unit: 'm', min: -1.5, max: 1.5, step: 0.005, help: 'Along the shake axis. Radiation-only: no re-solve needed.' },
209+ { key: 'listener_y_m', label: 'Listener y', unit: 'm', min: -1.5, max: 1.5, step: 0.005 },
210+ { key: 'listener_z_m', label: 'Listener z', unit: 'm', min: -1.5, max: 1.5, step: 0.005 },
211+ { key: 'direct_gain', label: 'Direct click gain', min: 0.0, max: 3.0, step: 0.05, help: 'Broadband forced response at the contact point, standing in for the truncated high-degree modes. 0 leaves only the resolved modes.' },
212+ ],
213+ },
214+ {
215+ name: 'Output',
216+ specs: [
217+ { key: 'tail_s', label: 'Tail', unit: 's', min: 0.05, max: 2.0, step: 0.05, help: 'Silence-and-decay appended after the last stroke.' },
218+ { key: 'pre_delay_ms', label: 'Pre-delay', unit: 'ms', min: 0, max: 50, step: 0.5 },
219+ { key: 'peak_dbfs', label: 'Peak level', unit: 'dBFS', min: -24, max: 0, step: 0.5, help: 'Used when normalize is on.' },
220+ { key: 'full_scale_pa', label: 'Full scale', unit: 'Pa', min: 1, max: 200, step: 1, help: 'Used when normalize is off: this many pascals maps to 0 dBFS.' },
221+ ],
222+ },
223+ {
224+ name: 'Advanced / quality',
225+ specs: [
226+ { key: 'max_mode_l', label: 'Max degree l', min: 2, max: 40, step: 1, integer: true, help: 'Highest spherical-harmonic degree retained.' },
227+ { key: 'max_mode_hz', label: 'Max mode freq', unit: 'Hz', min: 2000, max: 24000, step: 500 },
228+ { key: 'oversample', label: 'Oversample', min: 1, max: 8, step: 1, integer: true, help: 'Solver rate multiple; 2 for render, 1 for preview.' },
229+ { key: 'speed_of_sound', label: 'Speed of sound', unit: 'm/s', min: 300, max: 400, step: 1 },
230+ { key: 'air_density', label: 'Air density', unit: 'kg/m³', min: 0.9, max: 1.5, step: 0.01 },
231+ ],
232+ },
233+];
234+
235+const STORAGE_KEY = 'maraca-webapp-params-v1';
236+
237+export function loadParams(): MaracaParams {
238+ try {
239+ const raw = localStorage.getItem(STORAGE_KEY);
240+ if (!raw) return { ...DEFAULTS };
241+ const parsed = JSON.parse(raw);
242+ const out = { ...DEFAULTS };
243+ for (const k of Object.keys(DEFAULTS) as (keyof MaracaParams)[]) {
244+ if (typeof parsed[k] === typeof DEFAULTS[k]) {
245+ (out as Record<string, unknown>)[k] = parsed[k];
246+ }
247+ }
248+ return out;
249+ } catch {
250+ return { ...DEFAULTS };
251+ }
252+}
253+
254+export function saveParams(p: MaracaParams): void {
255+ try {
256+ localStorage.setItem(STORAGE_KEY, JSON.stringify(p));
257+ } catch {
258+ // storage unavailable; persistence is a convenience only
259+ }
260+}
webapp/src/selftest.tsadded+25−0View file
@@ -0,0 +1,25 @@
1+// Runs the MATLAB validation script (../matlab/maraca_validate.m) through
2+// numbl in Node: closed-form identities, exact reproducibility of the seeded
3+// impact 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 === 'maraca_validate.m')!;
17+const t0 = performance.now();
18+const result = executeCode(main.source, { displayResults: false }, workspaceFiles, 'maraca_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+429−0View file
@@ -0,0 +1,429 @@
1+// Web Worker that renders the maraca shake 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 -> maraca_basis + maraca_impacts (cached on solve-tier params)
6+// radiate tier -> maraca_radiate (cached on the listener and direct_gain)
7+// output tier -> maraca_output (always cheap, always re-run)
8+// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).
9+
10+import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
11+import { hashMFiles, type MFile } from '../mfiles';
12+import { SOLVE_KEYS, RADIATE_KEYS, type MaracaParams } from '../params';
13+
14+export interface RenderRequest {
15+ type: 'render';
16+ id: number;
17+ params: MaracaParams;
18+ quality: 'preview' | 'full';
19+ mfiles: MFile[];
20+}
21+
22+export interface RenderResponse {
23+ type: 'result';
24+ id: number;
25+ quality: 'preview' | 'full';
26+ y: Float32Array; // final signal, [-1, 1]
27+ sampleRate: number;
28+ impactT: Float32Array; // impact times, s
29+ impactV: Float32Array; // impact speeds, m/s
30+ impactPos: Float32Array; // 3 x I unit positions on the shell
31+ rateT: Float32Array; // expected-rate envelope, s
32+ rateHz: Float32Array; // expected impacts per second
33+ modeF: Float32Array; // mode frequencies, Hz
34+ modeL: Float32Array; // spherical-harmonic degree of each mode
35+ modeZeta: Float32Array; // damping ratio of each mode
36+ modeSigma: Float32Array; // radiation efficiency of each mode
37+ diag: {
38+ modesKept: number;
39+ breathingHz: number;
40+ breathingDryHz: number;
41+ impactCount: number;
42+ meanV: number;
43+ maxV: number;
44+ meanTauUs: number;
45+ durationS: number;
46+ peakPa: number;
47+ peakSplDb: number;
48+ directPeakPa: number;
49+ listenerDistM: number;
50+ warnings: string[];
51+ solveMs: number; // 0 when the cached solve was reused
52+ radiateMs: number; // 0 when the cached radiation was reused
53+ solveCached: boolean;
54+ radiateCached: boolean;
55+ oversample: number;
56+ };
57+}
58+
59+export interface ProgressMessage {
60+ type: 'progress';
61+ id: number;
62+ stage: 'solve' | 'radiate' | 'output';
63+ frac: number;
64+}
65+
66+// Raised by a fault in the .m code (the user can edit it), so it is reported
67+// rather than left to stall the request.
68+export interface ErrorMessage {
69+ type: 'error';
70+ id: number;
71+ message: string;
72+ file?: string;
73+ line?: number;
74+}
75+
76+// Impacts drawn for the raster plot. The full process runs to thousands, and
77+// plotting every one of them is neither visible nor worth the transfer.
78+const MAX_PLOT_IMPACTS = 4000;
79+
80+interface SolveEntry {
81+ timp: RuntimeValue;
82+ jimp: RuntimeValue;
83+ tauimp: RuntimeValue;
84+ fimp: RuntimeValue;
85+ ximp: RuntimeValue;
86+ fsHi: number;
87+ impactT: Float32Array;
88+ impactV: Float32Array;
89+ impactPos: Float32Array;
90+ rateT: Float32Array;
91+ rateHz: Float32Array;
92+ modeF: Float32Array;
93+ modeL: Float32Array;
94+ modeZeta: Float32Array;
95+ modeSigma: Float32Array;
96+ breathingHz: number;
97+ breathingDryHz: number;
98+ impactCount: number;
99+ meanV: number;
100+ maxV: number;
101+ meanTauUs: number;
102+ durationS: number;
103+}
104+
105+interface RadiateEntry {
106+ p: RuntimeValue;
107+ directPeakPa: number;
108+ listenerDistM: number;
109+ listenerInside: boolean;
110+}
111+
112+// Small LRU caches so preview- and full-quality intermediates coexist. A
113+// radiate-tier change right after a full render then reuses the full solve,
114+// and the "preview" is silently upgraded to full quality.
115+const solveCache = new Map<string, SolveEntry>();
116+const radiateCache = new Map<string, RadiateEntry>();
117+const SOLVE_CACHE_MAX = 2;
118+const RADIATE_CACHE_MAX = 3;
119+
120+function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
121+ const v = cache.get(key);
122+ if (v !== undefined) {
123+ cache.delete(key);
124+ cache.set(key, v);
125+ }
126+ return v;
127+}
128+
129+function lruSet<T>(cache: Map<string, T>, key: string, value: T, max: number): void {
130+ cache.delete(key);
131+ cache.set(key, value);
132+ while (cache.size > max) {
133+ const oldest = cache.keys().next().value as string;
134+ cache.delete(oldest);
135+ }
136+}
137+
138+function keyFor(params: MaracaParams, keys: readonly (keyof MaracaParams)[]): string {
139+ return keys.map((k) => String(params[k])).join('|');
140+}
141+
142+// MATLAB source lines assigning every parameter onto a struct P.
143+function pLiteral(params: MaracaParams): string {
144+ const lines: string[] = [];
145+ for (const [k, v] of Object.entries(params)) {
146+ if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
147+ else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
148+ }
149+ return lines.join('\n');
150+}
151+
152+function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
153+ return executeCode(
154+ source,
155+ { initialVariableValues: inputs, displayResults: false },
156+ mfiles,
157+ 'stage.m',
158+ );
159+}
160+
161+function asTensor(v: RuntimeValue): RuntimeTensor {
162+ const t = v as RuntimeTensor;
163+ if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
164+ return t;
165+}
166+
167+function asScalar(v: RuntimeValue): number {
168+ if (typeof v === 'number') return v;
169+ const t = v as RuntimeTensor;
170+ if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
171+ if (typeof v === 'boolean') return v ? 1 : 0;
172+ throw new Error('expected scalar result from numbl stage');
173+}
174+
175+function toF32(v: RuntimeValue): Float32Array {
176+ const t = asTensor(v);
177+ const out = new Float32Array(t.data.length);
178+ for (let i = 0; i < t.data.length; i++) out[i] = t.data[i];
179+ return out;
180+}
181+
182+const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
183+
184+self.onmessage = (ev: MessageEvent<RenderRequest>) => {
185+ const msg = ev.data;
186+ if (msg.type !== 'render') return;
187+
188+ // numbl reports a file it could not parse through console.warn and then
189+ // skips it, so the failure surfaces at the call site as "Undefined
190+ // function". Capture the warning to report the actual syntax error instead.
191+ const warnings: string[] = [];
192+ const origWarn = console.warn;
193+ console.warn = (...args: unknown[]) => {
194+ warnings.push(args.map(String).join(' '));
195+ };
196+
197+ try {
198+ render(msg);
199+ } catch (err) {
200+ // Most likely a fault in the (user-editable) .m code. Report it so the
201+ // request completes and the message reaches the code panel.
202+ const e = err as { message?: string; file?: string; line?: number };
203+ let m: ErrorMessage = {
204+ type: 'error',
205+ id: msg.id,
206+ message: e?.message ? String(e.message) : String(err),
207+ file: typeof e?.file === 'string' ? e.file : undefined,
208+ line: typeof e?.line === 'number' ? e.line : undefined,
209+ };
210+ for (const w of warnings) {
211+ const hit = SKIP_WARNING.exec(w);
212+ if (hit) {
213+ m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
214+ break;
215+ }
216+ }
217+ (self as unknown as Worker).postMessage(m);
218+ } finally {
219+ console.warn = origWarn;
220+ }
221+};
222+
223+function render(msg: RenderRequest): void {
224+ const { id, params, mfiles } = msg;
225+ let quality = msg.quality;
226+
227+ // The .m sources are part of the cache identity: editing the code must
228+ // invalidate everything computed from the previous version.
229+ const codeKey = hashMFiles(mfiles);
230+ const solveParamsKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);
231+ const fullSolveKey = solveParamsKey + `|os=${params.oversample}`;
232+ if (quality === 'preview' && solveCache.has(fullSolveKey)) {
233+ quality = 'full'; // full-rate impact process already cached
234+ }
235+
236+ const oversample = quality === 'preview' ? Math.min(1, params.oversample) : params.oversample;
237+ const fsHi = params.sample_rate * oversample;
238+
239+ const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
240+
241+ // ---- solve stage: mode basis and the stochastic impact process ----
242+ const solveKey = solveParamsKey + `|os=${oversample}`;
243+ let solveMs = 0;
244+ let solveCached = true;
245+ let solveEntry = lruGet(solveCache, solveKey);
246+ if (!solveEntry) {
247+ solveCached = false;
248+ post({ type: 'progress', id, stage: 'solve', frac: 0.15 });
249+ const t0 = performance.now();
250+ const r = runStage(mfiles, `
251+${pLiteral(params)}
252+basis = maraca_basis(P);
253+[timp, jimp, tauimp, fimp, ximp, ginfo] = maraca_impacts(P);
254+modeF = basis.f;
255+modeL = basis.l;
256+modeZeta = basis.zeta;
257+modeSigma = basis.sigma;
258+breathingHz = basis.breathingHz;
259+breathingDryHz = basis.breathingDryHz;
260+impactCount = ginfo.count;
261+meanV = ginfo.meanV;
262+maxV = ginfo.maxV;
263+meanTauUs = ginfo.meanTauUs;
264+durationS = ginfo.durationS;
265+rateT = ginfo.rateT;
266+rateHz = ginfo.rateHz;
267+vimp = ginfo.vimp;
268+`);
269+ solveMs = performance.now() - t0;
270+
271+ const impactCount = asScalar(r.variableValues.impactCount);
272+ const allT = toF32(r.variableValues.timp);
273+ const allV = toF32(r.variableValues.vimp);
274+ const allX = toF32(r.variableValues.ximp);
275+ // thin uniformly for the raster plot, keeping the time ordering
276+ const stride = Math.max(1, Math.ceil(impactCount / MAX_PLOT_IMPACTS));
277+ const nPlot = Math.ceil(impactCount / stride);
278+ const impactT = new Float32Array(nPlot);
279+ const impactV = new Float32Array(nPlot);
280+ const impactPos = new Float32Array(nPlot * 3);
281+ for (let i = 0, j = 0; j < nPlot; i += stride, j++) {
282+ impactT[j] = allT[i];
283+ impactV[j] = allV[i];
284+ impactPos[j * 3] = allX[i * 3];
285+ impactPos[j * 3 + 1] = allX[i * 3 + 1];
286+ impactPos[j * 3 + 2] = allX[i * 3 + 2];
287+ }
288+
289+ solveEntry = {
290+ timp: r.variableValues.timp,
291+ jimp: r.variableValues.jimp,
292+ tauimp: r.variableValues.tauimp,
293+ fimp: r.variableValues.fimp,
294+ ximp: r.variableValues.ximp,
295+ fsHi,
296+ impactT,
297+ impactV,
298+ impactPos,
299+ rateT: toF32(r.variableValues.rateT),
300+ rateHz: toF32(r.variableValues.rateHz),
301+ modeF: toF32(r.variableValues.modeF),
302+ modeL: toF32(r.variableValues.modeL),
303+ modeZeta: toF32(r.variableValues.modeZeta),
304+ modeSigma: toF32(r.variableValues.modeSigma),
305+ breathingHz: asScalar(r.variableValues.breathingHz),
306+ breathingDryHz: asScalar(r.variableValues.breathingDryHz),
307+ impactCount,
308+ meanV: asScalar(r.variableValues.meanV),
309+ maxV: asScalar(r.variableValues.maxV),
310+ meanTauUs: asScalar(r.variableValues.meanTauUs),
311+ durationS: asScalar(r.variableValues.durationS),
312+ };
313+ lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
314+ }
315+
316+ // ---- radiate stage ----
317+ const radiateKey = solveKey + '||' + keyFor(params, RADIATE_KEYS);
318+ let radiateMs = 0;
319+ let radiateCached = true;
320+ let radEntry = lruGet(radiateCache, radiateKey);
321+ if (!radEntry) {
322+ radiateCached = false;
323+ post({ type: 'progress', id, stage: 'radiate', frac: 0.55 });
324+ const t0 = performance.now();
325+ const r = runStage(mfiles, `
326+${pLiteral(params)}
327+basis = maraca_basis(P);
328+[p, radinfo] = maraca_radiate(P, basis, timp, jimp, tauimp, fimp, ximp, ${solveEntry.fsHi});
329+directPeakPa = radinfo.directPeakPa;
330+listenerDistM = radinfo.listenerDistM;
331+listenerInside = radinfo.listenerInside;
332+`, {
333+ timp: solveEntry.timp,
334+ jimp: solveEntry.jimp,
335+ tauimp: solveEntry.tauimp,
336+ fimp: solveEntry.fimp,
337+ ximp: solveEntry.ximp,
338+ });
339+ radiateMs = performance.now() - t0;
340+ radEntry = {
341+ p: r.variableValues.p,
342+ directPeakPa: asScalar(r.variableValues.directPeakPa),
343+ listenerDistM: asScalar(r.variableValues.listenerDistM),
344+ listenerInside: asScalar(r.variableValues.listenerInside) !== 0,
345+ };
346+ lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
347+ }
348+
349+ // ---- output stage ----
350+ post({ type: 'progress', id, stage: 'output', frac: 0.9 });
351+ const rOut = runStage(mfiles, `
352+${pLiteral(params)}
353+[y, peakPa] = maraca_output(P, p, ${oversample});
354+`, { p: radEntry.p });
355+
356+ const y32 = toF32(rOut.variableValues.y);
357+ const peakPa = asScalar(rOut.variableValues.peakPa);
358+
359+ const warnings: string[] = [];
360+ if (radEntry.listenerInside) {
361+ warnings.push(
362+ 'the listener is at or inside the shell surface; the multipole field is evaluated at the nearest valid radius instead',
363+ );
364+ }
365+ if (!params.normalize) {
366+ let peak = 0;
367+ for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
368+ if (peak > 1.0) {
369+ warnings.push(
370+ `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,
371+ );
372+ }
373+ }
374+ const topMode = solveEntry.modeF[solveEntry.modeF.length - 1];
375+ if (topMode < 0.4 * params.max_mode_hz) {
376+ warnings.push(
377+ `the highest retained mode is only ${topMode.toFixed(0)} Hz; raise max degree l to keep the shell from sounding truncated`,
378+ );
379+ }
380+
381+ const resp: RenderResponse = {
382+ type: 'result',
383+ id,
384+ quality,
385+ y: y32,
386+ sampleRate: params.sample_rate,
387+ impactT: solveEntry.impactT.slice(),
388+ impactV: solveEntry.impactV.slice(),
389+ impactPos: solveEntry.impactPos.slice(),
390+ rateT: solveEntry.rateT.slice(),
391+ rateHz: solveEntry.rateHz.slice(),
392+ modeF: solveEntry.modeF.slice(),
393+ modeL: solveEntry.modeL.slice(),
394+ modeZeta: solveEntry.modeZeta.slice(),
395+ modeSigma: solveEntry.modeSigma.slice(),
396+ diag: {
397+ modesKept: solveEntry.modeF.length,
398+ breathingHz: solveEntry.breathingHz,
399+ breathingDryHz: solveEntry.breathingDryHz,
400+ impactCount: solveEntry.impactCount,
401+ meanV: solveEntry.meanV,
402+ maxV: solveEntry.maxV,
403+ meanTauUs: solveEntry.meanTauUs,
404+ durationS: solveEntry.durationS,
405+ peakPa,
406+ peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
407+ directPeakPa: radEntry.directPeakPa,
408+ listenerDistM: radEntry.listenerDistM,
409+ warnings,
410+ solveMs,
411+ radiateMs,
412+ solveCached,
413+ radiateCached,
414+ oversample,
415+ },
416+ };
417+ (self as unknown as Worker).postMessage(resp, [
418+ resp.y.buffer,
419+ resp.impactT.buffer,
420+ resp.impactV.buffer,
421+ resp.impactPos.buffer,
422+ resp.rateT.buffer,
423+ resp.rateHz.buffer,
424+ resp.modeF.buffer,
425+ resp.modeL.buffer,
426+ resp.modeZeta.buffer,
427+ resp.modeSigma.buffer,
428+ ]);
429+}
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+});