Reimplement the algorithm in MATLAB syntax, run it with numbl
The model is now written in MATLAB (matlab/*.m) and executed by numbl in the
browser; the TypeScript physics port is removed. The worker runs each stage as
a separate numbl execution, passing tensors between them as RuntimeValues and
keeping the existing solve/radiate/output cache tiers.
The solver exploits the fact that the system is nonlinear only during contact:
RK4 runs until the mallet has provably separated, then the linear body+cavity
system is evolved in closed form from one eigendecomposition, evaluated
vectorized over the time axis. Validation against woodblock.py is unchanged
(contact 0.301 ms at 794 N, 6.132 Pa, and the four listener targets from the
spec), and the slowest eigenvalue comes out at the documented -34 s^-1.
Adds a code panel: the .m sources are viewable and editable in the app with
MATLAB syntax highlighting, edits re-render as you type, and they are held in
memory only so a reload restores the shipped code. The worker now reports
faults in the .m code instead of stalling, preferring numbl's parse warning
over the misleading "Undefined function" it raises at the call site.
28 changed files+3081−739
matlab/woodblock_basis.madded+69−0View file
@@ -0,0 +1,69 @@
1+function basis = woodblock_basis(P)
2+% WOODBLOCK_BASIS Modal basis: free-free Euler-Bernoulli bar.
3+% Mass-normalized bending modes and their couplings (woodblock.md section 3).
4+% basis.f modal frequencies (Hz), n x 1
5+% basis.phiC mode shape at the strike point (1/sqrt(kg))
6+% basis.A cavity volume coupling (m^2/sqrt(kg))
7+% basis.shapes mass-normalized shapes sampled on basis.s, nGrid x n
8+
9+L = P.length_m; W = P.width_m; H = P.wall_thickness_m;
10+mass = P.participating_mass_kg;
11+if mass == 0
12+ mass = P.density * W * H * L;
13+end
14+cLong = sqrt(P.youngs_modulus / P.density);
15+% f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)), I = W H^3/12
16+prefactor = cLong * H / (2 * pi * sqrt(12) * L^2);
17+
18+nGrid = 4096;
19+s = linspace(0, 1, nGrid).';
20+
21+% Roots of cos(bL)cosh(bL) = 1; asymptotic beyond the 4th
22+blTab = [4.73004074 7.85320462 10.99560784 14.13716549];
23+f = zeros(P.n_modes, 1);
24+phiC = zeros(P.n_modes, 1);
25+A = zeros(P.n_modes, 1);
26+shapes = zeros(nGrid, P.n_modes);
27+f1 = prefactor * blTab(1)^2 * P.hollow_softening;
28+n = 0;
29+for m = 1:P.n_modes
30+ if m <= 4
31+ bl = blTab(m);
32+ else
33+ bl = (2 * m + 1) * pi / 2;
34+ end
35+ fm = prefactor * bl^2 * P.hollow_softening;
36+ if m > 1 && P.inharmonicity ~= 1
37+ fm = f1 * (fm / f1)^P.inharmonicity;
38+ end
39+ if fm > P.max_mode_hz
40+ break;
41+ end
42+ n = n + 1;
43+ sigma = (cosh(bl) - cos(bl)) / (sinh(bl) - sin(bl));
44+ x = bl * s;
45+ shape = (cosh(x) + cos(x)) - sigma * (sinh(x) + sin(x));
46+ rmsval = sqrt(trapz(s, shape.^2));
47+ shape = shape / (rmsval * sqrt(mass));
48+ f(n) = fm;
49+ phiC(n) = interp1(s, shape, P.strike_position);
50+ mask = (s >= P.cavity_span_lo) & (s <= P.cavity_span_hi);
51+ A(n) = W * L * trapz(s, shape .* mask);
52+ shapes(:, n) = shape;
53+end
54+if n == 0
55+ error('No modes below max_mode_hz -- check geometry/material.');
56+end
57+f = f(1:n);
58+omega = 2 * pi * f;
59+% assemble by assignment (struct() with array values would make a struct array)
60+basis.n = n;
61+basis.f = f;
62+basis.omega = omega;
63+basis.zeta = 0.5 * (P.loss_factor + P.loss_slope * omega);
64+basis.phiC = phiC(1:n);
65+basis.A = A(1:n) * P.cavity_coupling;
66+basis.mass = mass;
67+basis.s = s;
68+basis.shapes = shapes(:, 1:n);
69+end
matlab/woodblock_defaults.madded+54−0View file
@@ -0,0 +1,54 @@
1+function P = woodblock_defaults()
2+% WOODBLOCK_DEFAULTS Default parameters, mirroring the web app's params.ts.
3+% Note woodblock.py's reference default for mallet_stiffness is 3.0e8; the
4+% app (and this file) default to a softer 1.0e8.
5+
6+P.length_m = 0.200;
7+P.width_m = 0.055;
8+P.wall_thickness_m = 0.010;
9+P.block_height_m = 0.050;
10+
11+P.density = 850.0;
12+P.youngs_modulus = 1.6e10;
13+P.hollow_softening = 0.90;
14+P.participating_mass_kg = 0.35;
15+P.inharmonicity = 1.00;
16+P.n_modes = 26;
17+P.max_mode_hz = 20000.0;
18+
19+P.loss_factor = 0.010;
20+P.loss_slope = 1.2e-7;
21+
22+P.strike_position = 0.42;
23+P.mallet_mass = 0.025;
24+P.mallet_velocity = 2.5;
25+P.mallet_stiffness = 1.0e8;
26+P.contact_exponent = 1.5;
27+P.contact_damping = 0.8;
28+
29+P.cavity_volume_m3 = 2.6e-4;
30+P.slit_area_m2 = 8.0e-4;
31+P.slit_length_m = 0.012;
32+P.cavity_damping = 0.06;
33+P.cavity_coupling = 1.0;
34+P.cavity_span_lo = 0.30;
35+P.cavity_span_hi = 0.70;
36+
37+P.listener_x_m = 0.15;
38+P.listener_y_m = 0.25;
39+P.listener_z_m = 0.30;
40+P.slit_x_m = 0.0;
41+P.slit_y_m = 0.0275;
42+P.slit_z_m = 0.0;
43+P.radiating_elements = 48;
44+P.speed_of_sound = 343.0;
45+P.air_density = 1.204;
46+
47+P.sample_rate = 48000;
48+P.oversample = 16;
49+P.duration_s = 0.45;
50+P.pre_delay_ms = 3.0;
51+P.normalize = 1;
52+P.peak_dbfs = -3.0;
53+P.full_scale_pa = 20.0;
54+end
matlab/woodblock_demo.madded+24−0View file
@@ -0,0 +1,24 @@
1+% Render a woodblock strike with default parameters and print diagnostics.
2+% Run with: numbl run woodblock_demo.m
3+
4+P = woodblock_defaults();
5+tic;
6+[y, basis, dinfo] = woodblock_render(P);
7+elapsed = toc;
8+
9+fprintf('modes kept : %d\n', basis.n);
10+fprintf('partials Hz: ');
11+fprintf('%.0f ', basis.f);
12+fprintf('\nratios : ');
13+fprintf('%.2f ', basis.f / basis.f(1));
14+fprintf('\nhelmholtz : %.0f Hz\n', dinfo.helmholtzHz);
15+fprintf('contact : %.3f ms, peak %.0f N\n', dinfo.contactMs, dinfo.peakForceN);
16+fprintf('peak press.: %.3f Pa (%.1f dB SPL)\n', dinfo.peakPa, dinfo.peakSplDb);
17+fprintf('length : %.3f s @ %d Hz\n', length(y) / P.sample_rate, P.sample_rate);
18+fprintf('render time: %.2f s\n', elapsed);
19+if dinfo.spacingNeed > 0
20+ fprintf('WARNING : element spacing is coarse; consider %d radiating elements\n', dinfo.spacingNeed);
21+end
22+if dinfo.listenerInside
23+ fprintf('WARNING : listener is inside the source region\n');
24+end
matlab/woodblock_output.madded+62−0View file
@@ -0,0 +1,62 @@
1+function [y, peakPa] = woodblock_output(P, p, oversample)
2+% WOODBLOCK_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 (woodblock.md section 8).
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 .* (wb_i0(beta * sqrt(max(0, 1 - r.^2))) / wb_i0(beta));
21+ h = h / sum(h); % unity DC gain
22+ pf = conv(p, h); % pf(i) = sum_k p(k) h(i-k+1)
23+ outLen = floor(T / R);
24+ y = pf((0:outLen - 1) * R + M2 + 1); % centered polyphase 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+total = floor(P.duration_s * P.sample_rate);
32+if padN > 0
33+ y = [zeros(1, padN) y];
34+end
35+y = y(1:min(length(y), total));
36+
37+% Squared fade on the tail
38+fade = min(256, floor(length(y) / 8));
39+if fade > 1
40+ y(end - fade + 1:end) = y(end - fade + 1:end) .* linspace(1, 0, fade).^2;
41+end
42+
43+if P.normalize ~= 0
44+ pk = max(abs(y));
45+ if pk > 0
46+ y = y / pk * 10^(P.peak_dbfs / 20);
47+ end
48+else
49+ y = y / P.full_scale_pa;
50+end
51+end
52+
53+function v = wb_i0(x)
54+% Modified Bessel function of the first kind, order zero (series expansion).
55+v = ones(size(x));
56+term = ones(size(x));
57+hx = x / 2;
58+for k = 1:40
59+ term = term .* (hx / k).^2;
60+ v = v + term;
61+end
62+end
matlab/woodblock_radiate.madded+57−0View file
@@ -0,0 +1,57 @@
1+function [p, spacingNeed, listenerInside] = woodblock_radiate(P, basis, Qdd, Xdd, fsHi, M)
2+% WOODBLOCK_RADIATE Retarded free-field pressure at the listener point.
3+% The top and bottom faces move together, so each strip is a local dipole;
4+% the strip surface accelerations are computed in one matrix product and
5+% reused for both faces with opposite sign (woodblock.md section 7).
6+% spacingNeed is 0 when the element spacing is fine enough, otherwise the
7+% recommended element count; listenerInside flags |r| < L/2.
8+
9+N = basis.n;
10+T = size(Qdd, 2);
11+L = P.length_m; W = P.width_m;
12+halfH = 0.5 * P.block_height_m;
13+patchA = W * L / M;
14+c = P.speed_of_sound;
15+rho0 = P.air_density;
16+
17+sElem = ((1:M).' - 0.5) / M;
18+phiElem = zeros(M, N);
19+for n = 1:N
20+ phiElem(:, n) = interp1(basis.s, basis.shapes(:, n), sElem);
21+end
22+xElem = (sElem - 0.5) * L;
23+
24+lx = P.listener_x_m; ly = P.listener_y_m; lz = P.listener_z_m;
25+dTop = sqrt((xElem - lx).^2 + ly^2 + (halfH - lz)^2);
26+dBot = sqrt((xElem - lx).^2 + ly^2 + (-halfH - lz)^2);
27+slitD = sqrt((P.slit_x_m - lx)^2 + (P.slit_y_m - ly)^2 + (P.slit_z_m - lz)^2);
28+
29+maxDelay = ceil(max([dTop; dBot; slitD]) / c * fsHi) + 2;
30+out = zeros(1, T + maxDelay);
31+
32+wddAll = phiElem * Qdd; % M x T strip surface accelerations
33+for j = 1:M
34+ aT = patchA * rho0 / (4 * pi * dTop(j));
35+ aB = -patchA * rho0 / (4 * pi * dBot(j));
36+ kT = round(dTop(j) / c * fsHi) + 1;
37+ kB = round(dBot(j) / c * fsHi) + 1;
38+ out(kT:kT + T - 1) = out(kT:kT + T - 1) + aT * wddAll(j, :);
39+ out(kB:kB + T - 1) = out(kB:kB + T - 1) + aB * wddAll(j, :);
40+end
41+
42+% Slit monopole
43+aS = P.slit_area_m2 * rho0 / (4 * pi * slitD);
44+kS = round(slitD / c * fsHi) + 1;
45+out(kS:kS + T - 1) = out(kS:kS + T - 1) + aS * Xdd;
46+
47+p = out(1:T);
48+
49+% Setup checks (check_setup in the reference)
50+fMax = basis.f(N);
51+if L / M > 0.25 * c / fMax
52+ spacingNeed = ceil(L / (0.25 * c / fMax));
53+else
54+ spacingNeed = 0;
55+end
56+listenerInside = sqrt(lx^2 + ly^2 + lz^2) < 0.5 * L;
57+end
matlab/woodblock_render.madded+23−0View file
@@ -0,0 +1,23 @@
1+function [y, basis, dinfo] = woodblock_render(P)
2+% WOODBLOCK_RENDER Full render: solve -> 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 = woodblock_basis(P);
7+fsHi = P.sample_rate * P.oversample;
8+[Qdd, Xdd, force, helmholtzHz, contactMs, peakForceN, maxRealLam] = ...
9+ woodblock_solve(P, basis, P.oversample);
10+[p, spacingNeed, listenerInside] = ...
11+ woodblock_radiate(P, basis, Qdd, Xdd, fsHi, P.radiating_elements);
12+[y, peakPa] = woodblock_output(P, p, P.oversample);
13+
14+dinfo.helmholtzHz = helmholtzHz;
15+dinfo.contactMs = contactMs;
16+dinfo.peakForceN = peakForceN;
17+dinfo.maxRealLam = maxRealLam;
18+dinfo.peakPa = peakPa;
19+dinfo.peakSplDb = 20 * log10(peakPa / 20e-6 + 1e-30);
20+dinfo.spacingNeed = spacingNeed;
21+dinfo.listenerInside = listenerInside;
22+dinfo.force = force;
23+end
matlab/woodblock_solve.madded+145−0View file
@@ -0,0 +1,145 @@
1+function [Qdd, Xdd, force, helmholtzHz, contactMs, peakForceN, maxRealLam] = woodblock_solve(P, basis, oversample)
2+% WOODBLOCK_SOLVE Integrate the coupled mallet / body / cavity system.
3+%
4+% The system is nonlinear only while the mallet is in contact (a fraction of
5+% a millisecond), so the solve is split in two:
6+% contact phase: RK4 on the full nonlinear system, until the mallet has
7+% rebounded far enough that no future contact is possible;
8+% free phase: the remaining coupled body+cavity system is linear, so it is
9+% evolved in closed form from one eigendecomposition of the (2N+2)
10+% system matrix, evaluated vectorized over the whole time axis.
11+% This departs from the reference woodblock.py, which runs RK4 throughout;
12+% the free-phase eigen solution is exact in time, so the two agree to the
13+% reference's own integration tolerance.
14+%
15+% Returns modal accelerations Qdd (N x T), slit-plug acceleration Xdd (1 x T),
16+% contact force (1 x T), and scalar diagnostics. maxRealLam should be
17+% strictly negative; anything else means the coupling was assembled wrongly.
18+
19+N = basis.n;
20+fsHi = P.sample_rate * oversample;
21+dt = 1 / fsHi;
22+nSteps = floor(P.duration_s * fsHi);
23+
24+om = basis.omega;
25+om2 = om.^2;
26+tzo = 2 * basis.zeta .* om;
27+phiC = basis.phiC;
28+A = basis.A;
29+
30+rho0 = P.air_density; c = P.speed_of_sound;
31+V = P.cavity_volume_m3; S = P.slit_area_m2;
32+% End correction on the slit: ~0.85 * equivalent radius per open end.
33+Leff = P.slit_length_m + 2 * 0.85 * sqrt(S / pi);
34+omH = c * sqrt(S / (V * Leff));
35+tzoH = 2 * P.cavity_damping * omH;
36+pScale = rho0 * c^2 / V;
37+cavDrive = c^2 / (V * Leff);
38+K = P.mallet_stiffness;
39+alpha = P.contact_exponent;
40+mu = P.contact_damping;
41+invMh = 1 / P.mallet_mass;
42+
43+% ---- contact phase: RK4 on [q; qd; xi; xid; y; yd] ----
44+maxContact = min(nSteps, ceil(0.02 * fsHi)); % 20 ms hard cap
45+z = zeros(2 * N + 4, 1);
46+z(2 * N + 4) = P.mallet_velocity;
47+Qc = zeros(N, maxContact);
48+Xc = zeros(1, maxContact);
49+forceC = zeros(1, maxContact);
50+Tc = 0;
51+for k = 1:maxContact
52+ [d1, qdd1, xidd1, fc1] = wb_deriv(z, N, om2, tzo, phiC, A, S, pScale, omH, tzoH, cavDrive, K, alpha, mu, invMh);
53+ Qc(:, k) = qdd1;
54+ Xc(k) = xidd1;
55+ forceC(k) = fc1;
56+ Tc = k;
57+ [d2, ~, ~, ~] = wb_deriv(z + 0.5 * dt * d1, N, om2, tzo, phiC, A, S, pScale, omH, tzoH, cavDrive, K, alpha, mu, invMh);
58+ [d3, ~, ~, ~] = wb_deriv(z + 0.5 * dt * d2, N, om2, tzo, phiC, A, S, pScale, omH, tzoH, cavDrive, K, alpha, mu, invMh);
59+ [d4, ~, ~, ~] = wb_deriv(z + dt * d3, N, om2, tzo, phiC, A, S, pScale, omH, tzoH, cavDrive, K, alpha, mu, invMh);
60+ z = z + (dt / 6) * (d1 + 2 * d2 + 2 * d3 + d4);
61+
62+ % No-recontact exit: after separation the mallet coasts at constant
63+ % velocity while the surface stays within B (a conservative bound from
64+ % the current modal energy), so once the mallet is receding and further
65+ % than 1.5*B below the surface rest plane, contact can never resume.
66+ y = z(2 * N + 3);
67+ yd = z(2 * N + 4);
68+ if fc1 == 0 && yd <= 0
69+ q = z(1:N);
70+ qd = z(N + 1:2 * N);
71+ B = sum(abs(phiC) .* sqrt(q.^2 + (qd ./ om).^2));
72+ if y <= -1.5 * B
73+ break;
74+ end
75+ end
76+end
77+
78+% ---- free phase: eigendecomposition of the linear body+cavity system ----
79+% state x = [q; qd; xi; xid], xdot = MM x
80+T2 = nSteps - Tc;
81+MM = zeros(2 * N + 2);
82+MM(1:N, N + 1:2 * N) = eye(N);
83+MM(N + 1:2 * N, 1:N) = -(diag(om2) + pScale * (A * A.'));
84+MM(N + 1:2 * N, N + 1:2 * N) = -diag(tzo);
85+MM(N + 1:2 * N, 2 * N + 1) = -pScale * S * A;
86+MM(2 * N + 1, 2 * N + 2) = 1;
87+MM(2 * N + 2, 1:N) = -cavDrive * A.';
88+MM(2 * N + 2, 2 * N + 1) = -omH^2;
89+MM(2 * N + 2, 2 * N + 2) = -tzoH;
90+
91+[VV, DD] = eig(MM);
92+lam = diag(DD);
93+maxRealLam = max(real(lam));
94+cvec = VV \ z(1:2 * N + 2);
95+w = lam .* cvec; % acceleration weights per eigenmode
96+t2 = (0:T2 - 1) * dt; % t2 = 0 is the state right after contact
97+E = exp(lam * t2); % (2N+2) x T2, complex
98+Qf = real((VV(N + 1:2 * N, :) * diag(w)) * E);
99+Xf = real((VV(2 * N + 2, :) * diag(w)) * E);
100+
101+Qdd = [Qc(:, 1:Tc) Qf];
102+Xdd = [Xc(1:Tc) Xf];
103+force = [forceC(1:Tc) zeros(1, T2)];
104+
105+helmholtzHz = omH / (2 * pi);
106+contactMs = sum(force > 0) * dt * 1e3;
107+peakForceN = max(force);
108+end
109+
110+function [d, qdd, xidd, fc] = wb_deriv(z, N, om2, tzo, phiC, A, S, pScale, omH, tzoH, cavDrive, K, alpha, mu, invMh)
111+% Derivative of the coupled system with Hunt-Crossley mallet contact.
112+q = z(1:N);
113+qd = z(N + 1:2 * N);
114+xi = z(2 * N + 1);
115+xid = z(2 * N + 2);
116+y = z(2 * N + 3);
117+yd = z(2 * N + 4);
118+
119+wsurf = phiC.' * q;
120+wdsurf = phiC.' * qd;
121+delta = y - wsurf;
122+if delta > 0
123+ % f = K d^a (1 + lambda * ddot); lambda in s/m keeps the correction O(1)
124+ fc = K * delta^alpha * (1 + mu * (yd - wdsurf));
125+ if fc < 0
126+ fc = 0;
127+ end
128+else
129+ fc = 0;
130+end
131+
132+Aq = A.' * q;
133+pc = pScale * (S * xi + Aq);
134+
135+qdd = -tzo .* qd - om2 .* q + fc * phiC - A * pc;
136+xidd = -tzoH * xid - omH^2 * xi - cavDrive * Aq;
137+
138+d = zeros(2 * N + 4, 1);
139+d(1:N) = qd;
140+d(N + 1:2 * N) = qdd;
141+d(2 * N + 1) = xid;
142+d(2 * N + 2) = xidd;
143+d(2 * N + 3) = yd;
144+d(2 * N + 4) = -fc * invMh;
145+end
matlab/woodblock_validate.madded+34−0View file
@@ -0,0 +1,34 @@
1+% Validate against woodblock.py reference diagnostics (its defaults use
2+% mallet_stiffness = 3.0e8) and the listener targets in woodblock.md §7.
3+% Run with: numbl run woodblock_validate.m
4+
5+P = woodblock_defaults();
6+P.mallet_stiffness = 3.0e8;
7+
8+[y, basis, dinfo] = woodblock_render(P);
9+fprintf('reference config:\n');
10+fprintf(' contact %.3f ms (target 0.301), peak %.0f N (target 794)\n', ...
11+ dinfo.contactMs, dinfo.peakForceN);
12+fprintf(' peak %.3f Pa (target 6.132), %.1f dB SPL (target 109.7)\n', ...
13+ dinfo.peakPa, dinfo.peakSplDb);
14+fprintf(' max real eigenvalue: %.1f (must be < 0)\n', dinfo.maxRealLam);
15+
16+targets = [0.15 0.25 0.30 109.7
17+ 0.00 0.00 0.40 100.9
18+ 0.00 0.40 0.00 90.9
19+ 0.00 0.05 0.02 125.4];
20+for i = 1:size(targets, 1)
21+ Q = P;
22+ Q.listener_x_m = targets(i, 1);
23+ Q.listener_y_m = targets(i, 2);
24+ Q.listener_z_m = targets(i, 3);
25+ [~, ~, di] = woodblock_render(Q);
26+ err = abs(di.peakSplDb - targets(i, 4));
27+ if err < 0.3
28+ verdict = 'ok';
29+ else
30+ verdict = 'MISMATCH';
31+ end
32+ fprintf('listener (%.2f, %.2f, %.2f): %.1f dB SPL (target %.1f) %s\n', ...
33+ targets(i, 1), targets(i, 2), targets(i, 3), di.peakSplDb, targets(i, 4), verdict);
34+end
webapp/README.mdmodified+22−8View file
@@ -1,6 +1,8 @@
11 # Woodblock Simulator — web app
22
3-Interactive browser port of the physical model in `../woodblock.py` (see `../woodblock.md` for the physics). The full chain runs client-side in a Web Worker: Hunt-Crossley mallet contact, free-free bar modes two-way coupled to a Helmholtz cavity, retarded-monopole radiation to the listener point, and Kaiser-windowed anti-aliased decimation to 48 kHz.
3+Interactive browser implementation of the physical model specified in `../woodblock.md` (reference implementation: `../woodblock.py`). 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: Hunt-Crossley mallet contact, free-free bar modes two-way coupled to a Helmholtz cavity, retarded-monopole radiation to the listener point, and Kaiser-windowed anti-aliased decimation to 48 kHz.
46
57 ## Usage
68
@@ -8,20 +10,32 @@ Interactive browser port of the physical model in `../woodblock.py` (see `../woo
810 npm install
911 npm run dev # development server
1012 npm run build # production build in dist/
11-npm run selftest # validate the physics port against the Python reference values
13+npm run selftest # run ../matlab/woodblock_validate.m through numbl in Node
1214 ```
1315
14-The self-test reproduces the reference diagnostics (partials 1003/2766/5423/8964/13390/18702 Hz, Helmholtz 484 Hz, contact 0.301 ms at 794 N, 109.7 dB SPL) and the four validation listener positions from `woodblock.md` §7.
16+The self-test reproduces the woodblock.py reference diagnostics (contact 0.301 ms at 794 N, 6.132 Pa / 109.7 dB SPL at the default listener) and the four validation listener positions from `woodblock.md` §7. The MATLAB code can also be run directly: `cd ../matlab && numbl run woodblock_demo.m` (or `woodblock_validate.m`).
17+
18+## The solver
19+
20+The reference implementation integrates the whole strike with RK4. The MATLAB implementation instead exploits that the system is nonlinear only while the mallet is in contact (a fraction of a millisecond): it runs RK4 until the mallet has provably separated (receding, and further from the surface than a bound derived from the current modal energy), then evolves the remaining coupled linear body+cavity system in closed form from one eigendecomposition of the (2N+2) system matrix, evaluated vectorized over the time axis. The free-phase solution is exact in time, and the results match the reference to display precision on every published diagnostic.
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. Editing `sqrt(12)` to `sqrt(48)` in the pitch prefactor, for instance, halves every partial.
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.
1527
1628 ## How rendering is tiered
1729
18-Changing any parameter triggers a debounced fast preview (oversample 8, 24 radiating elements); the Render button runs the full-quality pass (oversample 16, 48 elements). While the current result is a preview, the Play and Download buttons and the status badge say so. The worker keys its caches by tier: geometry, material, damping, mallet, and cavity parameters invalidate the modal solve; listener and slit positions and block height only re-run radiation; pre-delay, normalization, and bit depth only re-run the output stage. Small LRU caches mean that moving the microphone after a full render reuses the cached solve and returns a full-quality result in a fraction of a second.
30+Changing any parameter triggers a debounced fast preview (oversample 8, 24 radiating elements); the Render button runs the full-quality pass (oversample 16, 48 elements). 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: geometry, material, damping, mallet, and cavity parameters invalidate the solve; listener and slit positions and block height 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 solve and returns a full-quality result.
1931
20-The Loop button repeats the strike continuously, swapping in each new render as parameters change, so edits are audible without pressing Play. Parameters persist in localStorage; "Reset to defaults" restores the values from `woodblock.py`.
32+The Loop button repeats the strike continuously, swapping in each new render as parameters change, so edits are audible without pressing Play. Parameters persist in localStorage; "Reset to defaults" restores the values in `src/params.ts` (which `../matlab/woodblock_defaults.m` mirrors).
2133
2234 ## Layout
2335
36+- `../matlab/` — the algorithm: `woodblock_basis.m`, `woodblock_solve.m`, `woodblock_radiate.m`, `woodblock_output.m`, plus `woodblock_render.m` / `woodblock_demo.m` / `woodblock_validate.m` for CLI use
2437 - `src/params.ts` — parameter schema, defaults, tier assignment, persistence
25-- `src/physics/` — basis, RK4 solver, radiation, decimation (straight ports of the Python)
26-- `src/worker/renderWorker.ts` — off-thread rendering with tiered caching
27-- `src/components/` — 3D scene (three.js), parameter panel, plots (waveform, spectrum, spectrogram, contact force, animated mode shapes)
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 (waveform, spectrum, spectrogram, contact force, animated mode shapes)
41+- `src/selftest.ts` — Node harness for the MATLAB validation script
webapp/package-lock.jsonmodified+1885−16View file
This diff is 2,127 lines long and is not shown.
webapp/package.jsonmodified+9−1View file
@@ -7,9 +7,17 @@
77 "dev": "vite",
88 "build": "tsc --noEmit && vite build",
99 "preview": "vite preview",
10- "selftest": "tsx src/physics/selftest.ts"
10+ "selftest": "tsx src/selftest.ts"
1111 },
1212 "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",
1321 "react": "^18.3.1",
1422 "react-dom": "^18.3.1",
1523 "three": "^0.168.0"
webapp/src/App.tsxmodified+60−1View file
@@ -1,8 +1,10 @@
11 import { useCallback, useEffect, useRef, useState } from 'react';
22 import { DEFAULTS, loadParams, saveParams, type WoodblockParams } from './params';
3+import { DEFAULT_MFILES, type MFile } from './mfiles';
34 import { useRenderer } from './hooks/useRenderer';
45 import { encodeWav } from './lib/wav';
56 import ParamPanel from './components/ParamPanel';
7+import CodePanel from './components/CodePanel';
68 import Scene3D from './components/Scene3D';
79 import DiagnosticsPanel from './components/DiagnosticsPanel';
810 import WaveformPlot from './components/plots/WaveformPlot';
@@ -13,7 +15,11 @@ import ModeShapes from './components/plots/ModeShapes';
1315
1416 export default function App() {
1517 const [params, setParams] = useState<WoodblockParams>(() => loadParams());
16- const { result, busy, stage, progress, renderFull } = useRenderer(params);
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);
1723 const audioCtxRef = useRef<AudioContext | null>(null);
1824
1925 useEffect(() => saveParams(params), [params]);
@@ -24,6 +30,25 @@ export default function App() {
2430
2531 const reset = useCallback(() => setParams({ ...DEFAULTS }), []);
2632
33+ const onChangeFile = useCallback((name: string, source: string) => {
34+ setMfiles((fs) => fs.map((f) => (f.name === name ? { ...f, source } : f)));
35+ }, []);
36+
37+ const onResetFile = useCallback((name: string) => {
38+ setMfiles((fs) =>
39+ fs.map((f) =>
40+ f.name === name
41+ ? { ...f, source: DEFAULT_MFILES.find((d) => d.name === name)?.source ?? f.source }
42+ : f,
43+ ),
44+ );
45+ }, []);
46+
47+ const onResetAllFiles = useCallback(() => setMfiles(DEFAULT_MFILES.map((f) => ({ ...f }))), []);
48+ const codeEdited = mfiles.some(
49+ (f) => f.source !== DEFAULT_MFILES.find((d) => d.name === f.name)?.source,
50+ );
51+
2752 const [looping, setLooping] = useState(false);
2853 const loopSrcRef = useRef<AudioBufferSourceNode | null>(null);
2954
@@ -116,7 +141,15 @@ export default function App() {
116141 Download .wav{result?.quality === 'preview' ? ' (preview)' : ''}
117142 </button>
118143 <button onClick={reset}>Reset to defaults</button>
144+ <button
145+ className={codeOpen ? 'code active' : 'code'}
146+ onClick={() => setCodeOpen((v) => !v)}
147+ title="View and edit the MATLAB source that numbl executes"
148+ >
149+ {'{ }'} Code{codeEdited ? ' •' : ''}
150+ </button>
119151 <div className="status">
152+ {error && <span className="badge err">code error</span>}
120153 {result && (
121154 <span className={`badge ${result.quality}`}>
122155 {result.quality === 'full' ? 'full quality' : 'preview quality — press Render'}
@@ -139,6 +172,20 @@ export default function App() {
139172 <ParamPanel params={params} onChange={onChange} />
140173 </div>
141174 <div className="content">
175+ {error && (
176+ <div className="card wide error-card">
177+ <h2>MATLAB error</h2>
178+ <div>
179+ <strong>
180+ {error.file ? `${error.file}${error.line ? `:${error.line}` : ''}` : 'error'}
181+ </strong>{' '}
182+ {error.message}
183+ </div>
184+ <div className="code-hint">
185+ Showing the last successful render. Open Code to fix the source, or revert it.
186+ </div>
187+ </div>
188+ )}
142189 <div className="card wide">
143190 <h2>Geometry — block, strike, slit, microphone</h2>
144191 <Scene3D params={params} />
@@ -188,6 +235,18 @@ export default function App() {
188235 )}
189236 </div>
190237 </div>
238+ {codeOpen && (
239+ <CodePanel
240+ mfiles={mfiles}
241+ onChangeFile={onChangeFile}
242+ onResetFile={onResetFile}
243+ onResetAll={onResetAllFiles}
244+ onClose={() => setCodeOpen(false)}
245+ onRenderFull={renderFull}
246+ error={error}
247+ busy={busy}
248+ />
249+ )}
191250 </div>
192251 );
193252 }
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 woodblockHighlight = 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(woodblockHighlight),
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/hooks/useRenderer.tsmodified+69−33View file
@@ -1,78 +1,114 @@
1-// Hook owning the render worker: debounced auto-preview on parameter change,
2-// explicit full-quality render, and a latest-wins request queue.
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.
34
45 import { useCallback, useEffect, useRef, useState } from 'react';
56 import type { WoodblockParams } from '../params';
6-import type { ProgressMessage, RenderResponse } from '../worker/renderWorker';
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+}
715
816 export interface RenderState {
917 result: RenderResponse | null;
18+ error: RenderError | null;
1019 busy: boolean;
1120 stage: string;
1221 progress: number; // 0..1 of the current stage
1322 }
1423
15-export function useRenderer(params: WoodblockParams) {
24+interface Request {
25+ params: WoodblockParams;
26+ mfiles: MFile[];
27+ quality: 'preview' | 'full';
28+}
29+
30+export function useRenderer(params: WoodblockParams, mfiles: MFile[]) {
1631 const workerRef = useRef<Worker | null>(null);
1732 const nextId = useRef(1);
1833 const inFlight = useRef<number | null>(null);
19- const queued = useRef<{ params: WoodblockParams; quality: 'preview' | 'full' } | null>(null);
34+ const queued = useRef<Request | null>(null);
35+ const sendRef = useRef<(r: Request) => void>(() => {});
2036 const [state, setState] = useState<RenderState>({
2137 result: null,
38+ error: null,
2239 busy: false,
2340 stage: '',
2441 progress: 0,
2542 });
2643
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+
2758 useEffect(() => {
2859 const worker = new Worker(new URL('../worker/renderWorker.ts', import.meta.url), {
2960 type: 'module',
3061 });
3162 workerRef.current = worker;
32- worker.onmessage = (ev: MessageEvent<RenderResponse | ProgressMessage>) => {
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>) => {
3376 const msg = ev.data;
3477 if (msg.type === 'progress') {
3578 setState((s) => ({ ...s, stage: msg.stage, progress: msg.frac }));
3679 return;
3780 }
38- inFlight.current = null;
39- const q = queued.current;
40- if (q) {
41- queued.current = null;
42- // show the result we have, then immediately start the queued request
43- setState({ result: msg, busy: true, stage: 'queued', progress: 0 });
44- send(q.params, q.quality);
45- } else {
46- setState({ result: msg, busy: false, stage: '', progress: 0 });
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;
4791 }
92+ finish(() => ({ result: msg, error: null, busy: false, stage: '', progress: 0 }));
4893 };
4994 return () => {
5095 worker.terminate();
5196 workerRef.current = null;
97+ inFlight.current = null;
98+ queued.current = null;
5299 };
53- // eslint-disable-next-line react-hooks/exhaustive-deps
54- }, []);
55-
56- const send = useCallback((p: WoodblockParams, quality: 'preview' | 'full') => {
57- const worker = workerRef.current;
58- if (!worker) return;
59- if (inFlight.current !== null) {
60- queued.current = { params: p, quality };
61- return;
62- }
63- const id = nextId.current++;
64- inFlight.current = id;
65- setState((s) => ({ ...s, busy: true, stage: 'starting', progress: 0 }));
66- worker.postMessage({ type: 'render', id, params: p, quality });
67100 }, []);
68101
69- const renderFull = useCallback(() => send(params, 'full'), [params, send]);
102+ const renderFull = useCallback(
103+ () => send({ params, mfiles, quality: 'full' }),
104+ [params, mfiles, send],
105+ );
70106
71- // Debounced auto-preview whenever parameters change (including on mount).
107+ // Debounced auto-preview on any change to the parameters or the code.
72108 useEffect(() => {
73- const t = setTimeout(() => send(params, 'preview'), 250);
109+ const t = setTimeout(() => send({ params, mfiles, quality: 'preview' }), 250);
74110 return () => clearTimeout(t);
75- }, [params, send]);
111+ }, [params, mfiles, send]);
76112
77113 return { ...state, renderFull };
78114 }
webapp/src/mfiles.tsadded+0−0View file
Binary file not shown.
webapp/src/physics/basis.tsdeleted+0−145View file
@@ -1,145 +0,0 @@
1-// Modal basis: free-free Euler-Bernoulli bar, mass-normalized bending modes.
2-// Port of ModalBasis in woodblock.py.
3-
4-import type { WoodblockParams } from '../params';
5-
6-// Roots of cos(bL)cosh(bL) = 1 for the free-free beam.
7-const BETA_L = [4.73004074, 7.85320462, 10.99560784, 14.13716549];
8-
9-export function betaL(n: number): number {
10- if (n <= BETA_L.length) return BETA_L[n - 1];
11- return ((2 * n + 1) * Math.PI) / 2.0;
12-}
13-
14-// Free-free mode shape on s in [0,1], unnormalized.
15-export function modeShape(bl: number, s: Float64Array): Float64Array {
16- const sigma = (Math.cosh(bl) - Math.cos(bl)) / (Math.sinh(bl) - Math.sin(bl));
17- const out = new Float64Array(s.length);
18- for (let i = 0; i < s.length; i++) {
19- const x = bl * s[i];
20- out[i] = Math.cosh(x) + Math.cos(x) - sigma * (Math.sinh(x) + Math.sin(x));
21- }
22- return out;
23-}
24-
25-function trapz(y: Float64Array, ds: number): number {
26- let sum = 0;
27- for (let i = 0; i < y.length; i++) sum += y[i];
28- sum -= 0.5 * (y[0] + y[y.length - 1]);
29- return sum * ds;
30-}
31-
32-function interp1(x: number, xs: Float64Array, ys: Float64Array): number {
33- if (x <= xs[0]) return ys[0];
34- const n = xs.length;
35- if (x >= xs[n - 1]) return ys[n - 1];
36- // uniform grid
37- const dx = xs[1] - xs[0];
38- const j = Math.min(n - 2, Math.floor((x - xs[0]) / dx));
39- const t = (x - xs[j]) / dx;
40- return ys[j] * (1 - t) + ys[j + 1] * t;
41-}
42-
43-export interface ModalBasisData {
44- n: number;
45- f: Float64Array; // modal frequencies (Hz)
46- omega: Float64Array;
47- zeta: Float64Array;
48- phiC: Float64Array; // mode shape at strike point (1/sqrt(kg))
49- A: Float64Array; // cavity volume coupling (m^2/sqrt(kg))
50- mass: number;
51- // Mass-normalized shapes sampled on the dense s grid, for later
52- // interpolation at arbitrary radiating-element positions.
53- sGrid: Float64Array;
54- shapes: Float64Array[]; // length n, each sampled on sGrid
55-}
56-
57-const NS = 4096;
58-
59-export function buildBasis(cfg: WoodblockParams): ModalBasisData {
60- const L = cfg.length_m;
61- const W = cfg.width_m;
62- const H = cfg.wall_thickness_m;
63- const area = W * H;
64- const mass = cfg.participating_mass_kg || cfg.density * area * L;
65- const cLong = Math.sqrt(cfg.youngs_modulus / cfg.density);
66-
67- // f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)), I = W H^3/12
68- const prefactor = (cLong * H) / (2.0 * Math.PI * Math.sqrt(12.0) * L * L);
69-
70- const s = new Float64Array(NS);
71- for (let i = 0; i < NS; i++) s[i] = i / (NS - 1);
72- const ds = 1.0 / (NS - 1);
73-
74- const freqs: number[] = [];
75- const phiC: number[] = [];
76- const A: number[] = [];
77- const shapes: Float64Array[] = [];
78- const lo = cfg.cavity_span_lo;
79- const hi = cfg.cavity_span_hi;
80-
81- for (let n = 1; n <= cfg.n_modes; n++) {
82- const bl = betaL(n);
83- let f = prefactor * bl * bl * cfg.hollow_softening;
84- if (n > 1 && cfg.inharmonicity !== 1.0) {
85- const f0 = prefactor * betaL(1) ** 2 * cfg.hollow_softening;
86- f = f0 * (f / f0) ** cfg.inharmonicity;
87- }
88- if (f > cfg.max_mode_hz) break;
89-
90- const raw = modeShape(bl, s);
91- const sq = new Float64Array(NS);
92- for (let i = 0; i < NS; i++) sq[i] = raw[i] * raw[i];
93- const rms = Math.sqrt(trapz(sq, ds));
94- const norm = 1.0 / (rms * Math.sqrt(mass));
95- const shape = new Float64Array(NS);
96- for (let i = 0; i < NS; i++) shape[i] = raw[i] * norm;
97-
98- freqs.push(f);
99- phiC.push(interp1(cfg.strike_position, s, shape));
100- const masked = new Float64Array(NS);
101- for (let i = 0; i < NS; i++) masked[i] = s[i] >= lo && s[i] <= hi ? shape[i] : 0.0;
102- A.push(W * L * trapz(masked, ds));
103- shapes.push(shape);
104- }
105-
106- const nKept = freqs.length;
107- const fArr = Float64Array.from(freqs);
108- const omega = new Float64Array(nKept);
109- const zeta = new Float64Array(nKept);
110- const AArr = new Float64Array(nKept);
111- for (let i = 0; i < nKept; i++) {
112- omega[i] = 2 * Math.PI * fArr[i];
113- zeta[i] = 0.5 * (cfg.loss_factor + cfg.loss_slope * omega[i]);
114- AArr[i] = A[i] * cfg.cavity_coupling;
115- }
116-
117- return {
118- n: nKept,
119- f: fArr,
120- omega,
121- zeta,
122- phiC: Float64Array.from(phiC),
123- A: AArr,
124- mass,
125- sGrid: s,
126- shapes,
127- };
128-}
129-
130-// Interpolate the stored mode shapes at M radiating-element centres.
131-// Returns phiElem as an array of n Float64Arrays of length M, plus sElem.
132-export function elementShapes(basis: ModalBasisData, M: number): {
133- phiElem: Float64Array[];
134- sElem: Float64Array;
135-} {
136- const sElem = new Float64Array(M);
137- for (let j = 0; j < M; j++) sElem[j] = (j + 0.5) / M;
138- const phiElem: Float64Array[] = [];
139- for (let n = 0; n < basis.n; n++) {
140- const row = new Float64Array(M);
141- for (let j = 0; j < M; j++) row[j] = interp1(sElem[j], basis.sGrid, basis.shapes[n]);
142- phiElem.push(row);
143- }
144- return { phiElem, sElem };
145-}
webapp/src/physics/radiate.tsdeleted+0−105View file
@@ -1,105 +0,0 @@
1-// Radiation: retarded free-field pressure at the listener point.
2-// Port of radiate() and check_setup() in woodblock.py. The top and bottom
3-// faces move together, so each strip is a local dipole; the surface velocity
4-// sum is computed once per strip and reused for both faces with opposite sign.
5-
6-import type { WoodblockParams } from '../params';
7-import type { ModalBasisData } from './basis';
8-import { elementShapes } from './basis';
9-
10-export interface RadiateResult {
11- p: Float64Array; // pascals at the listener, at fsHi
12- warnings: string[];
13-}
14-
15-export function radiate(
16- cfg: WoodblockParams,
17- Qdd: Float64Array[],
18- Xdd: Float64Array,
19- basis: ModalBasisData,
20- fsHi: number,
21- M: number,
22- onProgress?: (frac: number) => void,
23-): RadiateResult {
24- const rho0 = cfg.air_density;
25- const c = cfg.speed_of_sound;
26- const N = basis.n;
27- const T = Xdd.length;
28- const L = cfg.length_m;
29- const W = cfg.width_m;
30- const halfH = 0.5 * cfg.block_height_m;
31- const patch = (W * L) / M;
32-
33- const { phiElem, sElem } = elementShapes(basis, M);
34-
35- const lx = cfg.listener_x_m;
36- const ly = cfg.listener_y_m;
37- const lz = cfg.listener_z_m;
38-
39- const distTop = new Float64Array(M);
40- const distBot = new Float64Array(M);
41- let maxDist = 0;
42- for (let j = 0; j < M; j++) {
43- const x = (sElem[j] - 0.5) * L;
44- const dx = x - lx;
45- const dy = 0.0 - ly;
46- const dzT = halfH - lz;
47- const dzB = -halfH - lz;
48- distTop[j] = Math.sqrt(dx * dx + dy * dy + dzT * dzT);
49- distBot[j] = Math.sqrt(dx * dx + dy * dy + dzB * dzB);
50- maxDist = Math.max(maxDist, distTop[j], distBot[j]);
51- }
52- const sx = cfg.slit_x_m - lx;
53- const sy = cfg.slit_y_m - ly;
54- const sz = cfg.slit_z_m - lz;
55- const slitDist = Math.sqrt(sx * sx + sy * sy + sz * sz);
56- maxDist = Math.max(maxDist, slitDist);
57-
58- const maxDelay = Math.ceil((maxDist / c) * fsHi) + 2;
59- const out = new Float64Array(T + maxDelay);
60- const wdd = new Float64Array(T);
61-
62- for (let j = 0; j < M; j++) {
63- // Surface acceleration of strip j, shared by the +z and -z sources
64- wdd.fill(0);
65- for (let n = 0; n < N; n++) {
66- const phi = phiElem[n][j];
67- if (phi === 0) continue;
68- const q = Qdd[n];
69- for (let t = 0; t < T; t++) wdd[t] += phi * q[t];
70- }
71- const ampTop = (patch * rho0) / (4.0 * Math.PI * distTop[j]);
72- const ampBot = -(patch * rho0) / (4.0 * Math.PI * distBot[j]);
73- const kTop = Math.round((distTop[j] / c) * fsHi);
74- const kBot = Math.round((distBot[j] / c) * fsHi);
75- for (let t = 0; t < T; t++) {
76- out[kTop + t] += ampTop * wdd[t];
77- out[kBot + t] += ampBot * wdd[t];
78- }
79- if (onProgress) onProgress((j + 1) / M);
80- }
81-
82- // Slit monopole
83- const ampSlit = (cfg.slit_area_m2 * rho0) / (4.0 * Math.PI * slitDist);
84- const kSlit = Math.round((slitDist / c) * fsHi);
85- for (let t = 0; t < T; t++) out[kSlit + t] += ampSlit * Xdd[t];
86-
87- // Setup warnings (check_setup in the reference)
88- const warnings: string[] = [];
89- const spacing = L / M;
90- const fMax = basis.f[basis.n - 1];
91- if (spacing > (0.25 * c) / fMax) {
92- const need = Math.ceil(L / ((0.25 * c) / fMax));
93- warnings.push(
94- `element spacing is coarse at ${fMax.toFixed(0)} Hz; consider ${need} radiating elements`,
95- );
96- }
97- const d = Math.sqrt(lx * lx + ly * ly + lz * lz);
98- if (d < 0.5 * L) {
99- warnings.push(
100- `listener is ${(d * 100).toFixed(0)} cm from the block centre, inside the source region; point-source elements are crude here`,
101- );
102- }
103-
104- return { p: out.subarray(0, T).slice(), warnings };
105-}
webapp/src/physics/render.tsdeleted+0−114View file
@@ -1,114 +0,0 @@
1-// Rendering orchestration: solve -> radiate -> output stage. Mirrors render()
2-// in woodblock.py, split into stages so the worker can cache intermediates.
3-
4-import type { WoodblockParams } from '../params';
5-import { buildBasis, type ModalBasisData } from './basis';
6-import { simulate, type SolveResult } from './solver';
7-import { radiate } from './radiate';
8-import { decimate } from './resample';
9-
10-export interface Diagnostics {
11- modesKept: number;
12- partialsHz: number[];
13- ratios: number[];
14- helmholtzHz: number;
15- contactMs: number;
16- peakForceN: number;
17- peakPa: number;
18- peakSplDb: number;
19- listenerDistM: number;
20- warnings: string[];
21-}
22-
23-export interface OutputStageResult {
24- y: Float64Array;
25- peakPa: number;
26- peakSplDb: number;
27- clipWarning: string | null;
28-}
29-
30-// Steps 3-5 of §8: decimate, pre-roll, fade, normalize or absolute scale.
31-export function outputStage(cfg: WoodblockParams, p: Float64Array, oversample: number): OutputStageResult {
32- let peakPa = 0;
33- for (let i = 0; i < p.length; i++) {
34- const a = Math.abs(p[i]);
35- if (a > peakPa) peakPa = a;
36- }
37- const peakSplDb = 20.0 * Math.log10(peakPa / 20e-6 + 1e-30);
38-
39- let y = decimate(p, oversample);
40-
41- const pad = Math.floor((cfg.pre_delay_ms * 1e-3) * cfg.sample_rate);
42- const total = Math.floor(cfg.duration_s * cfg.sample_rate);
43- if (pad > 0) {
44- const padded = new Float64Array(Math.min(pad + y.length, total));
45- padded.set(y.subarray(0, padded.length - pad), pad);
46- y = padded;
47- } else if (y.length > total) {
48- y = y.subarray(0, total).slice();
49- }
50-
51- const fade = Math.min(256, Math.floor(y.length / 8));
52- if (fade > 1) {
53- for (let i = 0; i < fade; i++) {
54- const g = 1.0 - i / (fade - 1);
55- y[y.length - fade + i] *= g * g;
56- }
57- }
58-
59- let clipWarning: string | null = null;
60- if (cfg.normalize) {
61- let peak = 0;
62- for (let i = 0; i < y.length; i++) peak = Math.max(peak, Math.abs(y[i]));
63- if (peak > 0) {
64- const g = (10 ** (cfg.peak_dbfs / 20.0)) / peak;
65- for (let i = 0; i < y.length; i++) y[i] *= g;
66- }
67- } else {
68- let peak = 0;
69- for (let i = 0; i < y.length; i++) {
70- y[i] /= cfg.full_scale_pa;
71- peak = Math.max(peak, Math.abs(y[i]));
72- }
73- if (peak > 1.0) {
74- clipWarning = `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${cfg.full_scale_pa.toFixed(2)}`;
75- }
76- }
77-
78- return { y, peakPa, peakSplDb, clipWarning };
79-}
80-
81-export interface FullRenderResult {
82- y: Float64Array;
83- basis: ModalBasisData;
84- solve: SolveResult;
85- diag: Diagnostics;
86-}
87-
88-// One-shot render, used by the self-test; the worker calls the stages
89-// individually so it can cache them.
90-export function renderAll(cfg: WoodblockParams, oversample?: number, elements?: number): FullRenderResult {
91- const os = oversample ?? cfg.oversample;
92- const M = elements ?? cfg.radiating_elements;
93- const basis = buildBasis(cfg);
94- const solve = simulate(cfg, basis, os);
95- const rad = radiate(cfg, solve.Qdd, solve.Xdd, basis, solve.fsHi, M);
96- const out = outputStage(cfg, rad.p, os);
97-
98- const lx = cfg.listener_x_m, ly = cfg.listener_y_m, lz = cfg.listener_z_m;
99- const warnings = [...rad.warnings];
100- if (out.clipWarning) warnings.push(out.clipWarning);
101- const diag: Diagnostics = {
102- modesKept: basis.n,
103- partialsHz: Array.from(basis.f),
104- ratios: Array.from(basis.f).map((f) => f / basis.f[0]),
105- helmholtzHz: solve.helmholtzHz,
106- contactMs: solve.contactMs,
107- peakForceN: solve.peakForceN,
108- peakPa: out.peakPa,
109- peakSplDb: out.peakSplDb,
110- listenerDistM: Math.sqrt(lx * lx + ly * ly + lz * lz),
111- warnings,
112- };
113- return { y: out.y, basis, solve, diag };
114-}
webapp/src/physics/resample.tsdeleted+0−56View file
@@ -1,56 +0,0 @@
1-// Anti-aliased decimation: Kaiser-windowed-sinc FIR lowpass at ~0.45 * fs_lo,
2-// evaluated only at the retained output samples (polyphase, so the cost is
3-// taps * T / R rather than taps * T). Replaces scipy.resample_poly(p, 1, R).
4-
5-function besselI0(x: number): number {
6- // Series expansion; converges quickly for the argument range used here.
7- let sum = 1.0;
8- let term = 1.0;
9- const half = x / 2.0;
10- for (let k = 1; k < 64; k++) {
11- term *= (half / k) * (half / k);
12- sum += term;
13- if (term < 1e-16 * sum) break;
14- }
15- return sum;
16-}
17-
18-export function designDecimator(R: number, beta = 8.0): Float64Array {
19- // Cutoff at 0.45 of the output Nyquist-setting rate: fc/fs_hi = 0.45/R.
20- const taps = 12 * R + 1; // odd length, linear phase
21- const h = new Float64Array(taps);
22- const M2 = (taps - 1) / 2;
23- const fcNorm = 0.45 / R; // cycles per input sample
24- const i0beta = besselI0(beta);
25- let sum = 0;
26- for (let i = 0; i < taps; i++) {
27- const m = i - M2;
28- const sinc = m === 0 ? 2 * fcNorm : Math.sin(2 * Math.PI * fcNorm * m) / (Math.PI * m);
29- const r = m / M2;
30- const w = besselI0(beta * Math.sqrt(Math.max(0, 1 - r * r))) / i0beta;
31- h[i] = sinc * w;
32- sum += h[i];
33- }
34- for (let i = 0; i < taps; i++) h[i] /= sum; // unity DC gain
35- return h;
36-}
37-
38-export function decimate(p: Float64Array, R: number): Float64Array {
39- if (R === 1) return p.slice();
40- const h = designDecimator(R);
41- const taps = h.length;
42- const M2 = (taps - 1) / 2;
43- const T = p.length;
44- const outLen = Math.floor(T / R);
45- const y = new Float64Array(outLen);
46- for (let m = 0; m < outLen; m++) {
47- const centre = m * R;
48- let acc = 0;
49- // y[m] = sum_i h[i] * p[centre + M2 - i], zero outside the buffer
50- const iLo = Math.max(0, centre + M2 - (T - 1));
51- const iHi = Math.min(taps - 1, centre + M2);
52- for (let i = iLo; i <= iHi; i++) acc += h[i] * p[centre + M2 - i];
53- y[m] = acc;
54- }
55- return y;
56-}
webapp/src/physics/selftest.tsdeleted+0−42View file
@@ -1,42 +0,0 @@
1-// Validation against the Python reference (woodblock.py at defaults):
2-// modes kept : 6
3-// partials Hz: 1003, 2766, 5423, 8964, 13390, 18702
4-// helmholtz : 484 Hz
5-// contact : 0.301 ms, peak 794 N
6-// peak press.: 6.132 Pa (109.7 dB SPL)
7-// Run with: npm run selftest
8-
9-import { DEFAULTS } from '../params';
10-import { renderAll } from './render';
11-
12-// The app's default mallet stiffness is softer than woodblock.py's; pin the
13-// reference value here so the comparison against the Python output holds.
14-const REFERENCE = { ...DEFAULTS, mallet_stiffness: 3.0e8 };
15-
16-const t0 = performance.now();
17-const { y, basis, diag } = renderAll(REFERENCE);
18-const t1 = performance.now();
19-
20-console.log(`modes kept : ${basis.n}`);
21-console.log('partials Hz: ' + Array.from(basis.f).map((f) => f.toFixed(0)).join(', '));
22-console.log('ratios : ' + Array.from(basis.f).map((f) => (f / basis.f[0]).toFixed(2)).join(', '));
23-console.log(`helmholtz : ${diag.helmholtzHz.toFixed(0)} Hz`);
24-console.log(`contact : ${diag.contactMs.toFixed(3)} ms, peak ${diag.peakForceN.toFixed(0)} N`);
25-console.log(`peak press.: ${diag.peakPa.toFixed(3)} Pa (${diag.peakSplDb.toFixed(1)} dB SPL)`);
26-console.log(`length : ${(y.length / REFERENCE.sample_rate).toFixed(3)} s @ ${REFERENCE.sample_rate} Hz`);
27-for (const w of diag.warnings) console.log(`WARNING : ${w}`);
28-console.log(`render time: ${((t1 - t0) / 1000).toFixed(2)} s`);
29-
30-// Validation listener positions from the spec (§7), peak SPL targets.
31-const targets: [number, number, number, number][] = [
32- [0.15, 0.25, 0.30, 109.7],
33- [0.00, 0.00, 0.40, 100.9],
34- [0.00, 0.40, 0.00, 90.9],
35- [0.00, 0.05, 0.02, 125.4],
36-];
37-for (const [x, yy, z, target] of targets) {
38- const cfg = { ...REFERENCE, listener_x_m: x, listener_y_m: yy, listener_z_m: z };
39- const r = renderAll(cfg);
40- const ok = Math.abs(r.diag.peakSplDb - target) < 0.5 ? 'ok' : 'MISMATCH';
41- console.log(`listener (${x}, ${yy}, ${z}): ${r.diag.peakSplDb.toFixed(1)} dB SPL (target ${target}) ${ok}`);
42-}
webapp/src/physics/solver.tsdeleted+0−157View file
@@ -1,157 +0,0 @@
1-// Time-domain RK4 integration of the coupled mallet / body / cavity system.
2-// Port of simulate() in woodblock.py. State: [q(N), qdot(N), xi, xidot, y, ydot].
3-
4-import type { WoodblockParams } from '../params';
5-import type { ModalBasisData } from './basis';
6-
7-export interface SolveResult {
8- Qdd: Float64Array[]; // N arrays of length T: modal accelerations
9- Xdd: Float64Array; // slit-plug acceleration
10- force: Float64Array; // contact force f_c(t)
11- fsHi: number;
12- helmholtzHz: number;
13- peakForceN: number;
14- contactMs: number;
15-}
16-
17-export function simulate(
18- cfg: WoodblockParams,
19- basis: ModalBasisData,
20- oversample: number,
21- onProgress?: (frac: number) => void,
22-): SolveResult {
23- const N = basis.n;
24- if (N === 0) throw new Error('No modes below max_mode_hz -- check geometry/material.');
25-
26- const fsHi = cfg.sample_rate * oversample;
27- const dt = 1.0 / fsHi;
28- const nSteps = Math.floor(cfg.duration_s * fsHi);
29-
30- const om2 = new Float64Array(N);
31- const twoZetaOm = new Float64Array(N);
32- for (let i = 0; i < N; i++) {
33- om2[i] = basis.omega[i] * basis.omega[i];
34- twoZetaOm[i] = 2.0 * basis.zeta[i] * basis.omega[i];
35- }
36- const phiC = basis.phiC;
37- const A = basis.A;
38-
39- const rho0 = cfg.air_density;
40- const c = cfg.speed_of_sound;
41- const V = cfg.cavity_volume_m3;
42- const S = cfg.slit_area_m2;
43- // End correction on the slit: ~0.85 * equivalent radius per open end.
44- const Leff = cfg.slit_length_m + 2.0 * 0.85 * Math.sqrt(S / Math.PI);
45- const omegaH = c * Math.sqrt(S / (V * Leff));
46- const twoZetaHOm = 2.0 * cfg.cavity_damping * omegaH;
47- const omegaH2 = omegaH * omegaH;
48- const pScale = (rho0 * c * c) / V;
49- const cavDrive = (c * c) / (V * Leff);
50-
51- const K = cfg.mallet_stiffness;
52- const alpha = cfg.contact_exponent;
53- const mu = cfg.contact_damping;
54- const invMh = 1.0 / cfg.mallet_mass;
55-
56- const dim = 2 * N + 4;
57- const IX = 2 * N;
58- const IXD = 2 * N + 1;
59- const IY = 2 * N + 2;
60- const IYD = 2 * N + 3;
61-
62- const state = new Float64Array(dim);
63- state[IYD] = cfg.mallet_velocity;
64-
65- let lastFc = 0.0;
66- let lastXidd = 0.0;
67-
68- const qddScratch = new Float64Array(N);
69-
70- function deriv(s: Float64Array, d: Float64Array): void {
71- // Mallet contact
72- let w = 0.0;
73- let wd = 0.0;
74- for (let i = 0; i < N; i++) {
75- w += phiC[i] * s[i];
76- wd += phiC[i] * s[N + i];
77- }
78- const y = s[IY];
79- const yd = s[IYD];
80- const delta = y - w;
81- let fc = 0.0;
82- if (delta > 0.0) {
83- // Hunt-Crossley: f = K d^a (1 + lambda * ddot); lambda is in s/m so the
84- // correction stays dimensionless and O(1) over a sub-ms collision.
85- fc = K * Math.pow(delta, alpha) * (1.0 + mu * (yd - wd));
86- if (fc < 0.0) fc = 0.0;
87- }
88- lastFc = fc;
89-
90- // Cavity pressure from combined volume change of walls and slit plug
91- let Aq = 0.0;
92- for (let i = 0; i < N; i++) Aq += A[i] * s[i];
93- const pc = pScale * (S * s[IX] + Aq);
94-
95- for (let i = 0; i < N; i++) {
96- const qdd = -twoZetaOm[i] * s[N + i] - om2[i] * s[i] + fc * phiC[i] - A[i] * pc;
97- d[i] = s[N + i];
98- d[N + i] = qdd;
99- qddScratch[i] = qdd;
100- }
101- const xidd = -twoZetaHOm * s[IXD] - omegaH2 * s[IX] - cavDrive * Aq;
102- d[IX] = s[IXD];
103- d[IXD] = xidd;
104- d[IY] = yd;
105- d[IYD] = -fc * invMh;
106- lastXidd = xidd;
107- }
108-
109- const Qdd: Float64Array[] = [];
110- for (let i = 0; i < N; i++) Qdd.push(new Float64Array(nSteps));
111- const Xdd = new Float64Array(nSteps);
112- const force = new Float64Array(nSteps);
113-
114- const k1 = new Float64Array(dim);
115- const k2 = new Float64Array(dim);
116- const k3 = new Float64Array(dim);
117- const k4 = new Float64Array(dim);
118- const tmp = new Float64Array(dim);
119-
120- const progressEvery = Math.max(1, Math.floor(nSteps / 50));
121-
122- for (let k = 0; k < nSteps; k++) {
123- deriv(state, k1);
124- for (let i = 0; i < N; i++) Qdd[i][k] = qddScratch[i];
125- Xdd[k] = lastXidd;
126- force[k] = lastFc;
127-
128- for (let i = 0; i < dim; i++) tmp[i] = state[i] + 0.5 * dt * k1[i];
129- deriv(tmp, k2);
130- for (let i = 0; i < dim; i++) tmp[i] = state[i] + 0.5 * dt * k2[i];
131- deriv(tmp, k3);
132- for (let i = 0; i < dim; i++) tmp[i] = state[i] + dt * k3[i];
133- deriv(tmp, k4);
134- for (let i = 0; i < dim; i++) {
135- state[i] += (dt / 6.0) * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);
136- }
137-
138- if (onProgress && k % progressEvery === 0) onProgress(k / nSteps);
139- }
140-
141- let peakForce = 0.0;
142- let contactSamples = 0;
143- for (let k = 0; k < nSteps; k++) {
144- if (force[k] > peakForce) peakForce = force[k];
145- if (force[k] > 0.0) contactSamples++;
146- }
147-
148- return {
149- Qdd,
150- Xdd,
151- force,
152- fsHi,
153- helmholtzHz: omegaH / (2 * Math.PI),
154- peakForceN: peakForce,
155- contactMs: (contactSamples * dt) * 1e3,
156- };
157-}
webapp/src/selftest.tsadded+24−0View file
@@ -0,0 +1,24 @@
1+// Runs the MATLAB validation script (../matlab/woodblock_validate.m) through
2+// numbl in Node, checking the implementation against the woodblock.py
3+// reference diagnostics and the listener targets in woodblock.md §7.
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 === 'woodblock_validate.m')!;
17+const t0 = performance.now();
18+const result = executeCode(main.source, { displayResults: false }, workspaceFiles, 'woodblock_validate.m');
19+process.stdout.write(result.output.join(''));
20+console.log(`total time : ${((performance.now() - t0) / 1000).toFixed(2)} s`);
21+
22+if (result.output.join('').includes('MISMATCH')) {
23+ process.exit(1);
24+}
webapp/src/style.cssmodified+123−0View file
@@ -79,6 +79,8 @@ button.loop.active {
7979 border-color: var(--accent2);
8080 font-weight: 600;
8181 }
82+button.code { font-family: ui-monospace, SFMono-Regular, Menlo, monospace; }
83+button.code.active { border-color: var(--accent2); color: var(--accent2); }
8284
8385 .status {
8486 margin-left: auto;
@@ -99,6 +101,7 @@ button.loop.active {
99101 }
100102 .badge.preview { color: var(--accent); border-color: var(--accent); }
101103 .badge.full { color: var(--good); border-color: var(--good); }
104+.badge.err { color: var(--warn); border-color: var(--warn); }
102105
103106 .progressbar {
104107 width: 90px;
@@ -230,6 +233,126 @@ details.group[open] summary { border-bottom: 1px solid var(--border); }
230233 .mode-list button { padding: 3px 10px; font-size: 12px; }
231234 .mode-list button.active { border-color: var(--accent2); color: var(--accent2); }
232235
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+
233356 @media (max-width: 1100px) {
234357 .content { grid-template-columns: 1fr; }
235358 }
webapp/src/vite-env.d.tsadded+1−0View file
@@ -0,0 +1 @@
1+/// <reference types="vite/client" />
webapp/src/worker/renderWorker.tsmodified+217−60View file
@@ -1,20 +1,22 @@
1-// Web Worker that runs the simulation off the main thread and caches
2-// intermediates so cheap parameter changes skip the expensive stages:
3-// solve tier -> ODE solve (cached on solve-tier params + oversample)
4-// radiate tier -> retarded-source sum (cached on radiate-tier params + M)
5-// output tier -> decimation + normalization (always cheap, always re-run)
1+// Web Worker that renders the woodblock strike by executing the MATLAB
2+// implementation in ../matlab through numbl. The stages run as separate
3+// executeCode calls so intermediates can be cached and cheap parameter
4+// changes skip the expensive ones:
5+// solve tier -> woodblock_solve (cached on solve-tier params + oversample)
6+// radiate tier -> woodblock_radiate (cached on radiate-tier params + M)
7+// output tier -> woodblock_output (always cheap, always re-run)
8+// Tensors cross the boundary as numbl RuntimeValues (column-major Float64).
69
10+import { executeCode, type RuntimeValue, type RuntimeTensor } from 'numbl';
11+import { hashMFiles, type MFile } from '../mfiles';
712 import { SOLVE_KEYS, RADIATE_KEYS, type WoodblockParams } from '../params';
8-import { buildBasis, type ModalBasisData } from '../physics/basis';
9-import { simulate, type SolveResult } from '../physics/solver';
10-import { radiate } from '../physics/radiate';
11-import { outputStage } from '../physics/render';
1213
1314 export interface RenderRequest {
1415 type: 'render';
1516 id: number;
1617 params: WoodblockParams;
1718 quality: 'preview' | 'full';
19+ mfiles: MFile[];
1820 }
1921
2022 export interface RenderResponse {
@@ -25,7 +27,7 @@ export interface RenderResponse {
2527 sampleRate: number;
2628 force: Float32Array; // contact force, first forceMs at fsHi
2729 forceFsHz: number;
28- modeShapes: Float32Array; // n x SHAPE_PTS, concatenated
30+ modeShapes: Float32Array; // n x shapePts, concatenated
2931 shapePts: number;
3032 diag: {
3133 modesKept: number;
@@ -54,14 +56,41 @@ export interface ProgressMessage {
5456 frac: number;
5557 }
5658
57-const SHAPE_PTS = 129;
59+// Raised by a fault in the .m code (the user can edit it), so it is reported
60+// rather than left to stall the request.
61+export interface ErrorMessage {
62+ type: 'error';
63+ id: number;
64+ message: string;
65+ file?: string;
66+ line?: number;
67+}
68+
69+const SHAPE_PTS = 128;
5870 const FORCE_MS = 8;
5971
72+interface SolveEntry {
73+ Qdd: RuntimeValue;
74+ Xdd: RuntimeValue;
75+ fsHi: number;
76+ f: number[];
77+ modeShapes: Float32Array;
78+ force: Float32Array;
79+ helmholtzHz: number;
80+ contactMs: number;
81+ peakForceN: number;
82+}
83+
84+interface RadiateEntry {
85+ p: RuntimeValue;
86+ warnings: string[];
87+}
88+
6089 // Small LRU caches so preview- and full-quality intermediates coexist. A
6190 // radiate-tier change right after a full render then reuses the full solve,
6291 // and the "preview" is silently upgraded to full quality.
63-const solveCache = new Map<string, { basis: ModalBasisData; solve: SolveResult }>();
64-const radiateCache = new Map<string, { p: Float64Array; warnings: string[] }>();
92+const solveCache = new Map<string, SolveEntry>();
93+const radiateCache = new Map<string, RadiateEntry>();
6594 const SOLVE_CACHE_MAX = 2;
6695 const RADIATE_CACHE_MAX = 3;
6796
@@ -69,7 +98,7 @@ function lruGet<T>(cache: Map<string, T>, key: string): T | undefined {
6998 const v = cache.get(key);
7099 if (v !== undefined) {
71100 cache.delete(key);
72- cache.set(key, v); // refresh recency
101+ cache.set(key, v);
73102 }
74103 return v;
75104 }
@@ -87,13 +116,88 @@ function keyFor(params: WoodblockParams, keys: readonly (keyof WoodblockParams)[
87116 return keys.map((k) => String(params[k])).join('|');
88117 }
89118
119+// MATLAB source lines assigning every parameter onto a struct P.
120+function pLiteral(params: WoodblockParams): string {
121+ const lines: string[] = [];
122+ for (const [k, v] of Object.entries(params)) {
123+ if (typeof v === 'number') lines.push(`P.${k} = ${v};`);
124+ else if (typeof v === 'boolean') lines.push(`P.${k} = ${v ? 1 : 0};`);
125+ }
126+ return lines.join('\n');
127+}
128+
129+function runStage(mfiles: MFile[], source: string, inputs?: Record<string, RuntimeValue>) {
130+ return executeCode(
131+ source,
132+ { initialVariableValues: inputs, displayResults: false },
133+ mfiles,
134+ 'stage.m',
135+ );
136+}
137+
138+function asTensor(v: RuntimeValue): RuntimeTensor {
139+ const t = v as RuntimeTensor;
140+ if (!t || t.kind !== 'tensor') throw new Error('expected tensor result from numbl stage');
141+ return t;
142+}
143+
144+function asScalar(v: RuntimeValue): number {
145+ if (typeof v === 'number') return v;
146+ const t = v as RuntimeTensor;
147+ if (t && t.kind === 'tensor' && t.data.length === 1) return t.data[0];
148+ if (typeof v === 'boolean') return v ? 1 : 0;
149+ throw new Error('expected scalar result from numbl stage');
150+}
151+
152+const SKIP_WARNING = /skipping (\S+\.m) \(syntax error at line (\d+)\)/;
153+
90154 self.onmessage = (ev: MessageEvent<RenderRequest>) => {
91155 const msg = ev.data;
92156 if (msg.type !== 'render') return;
93- const { id, params } = msg;
157+
158+ // numbl reports a file it could not parse through console.warn and then
159+ // skips it, so the failure surfaces at the call site as "Undefined
160+ // function". Capture the warning to report the actual syntax error instead.
161+ const warnings: string[] = [];
162+ const origWarn = console.warn;
163+ console.warn = (...args: unknown[]) => {
164+ warnings.push(args.map(String).join(' '));
165+ };
166+
167+ try {
168+ render(msg);
169+ } catch (err) {
170+ // Most likely a fault in the (user-editable) .m code. Report it so the
171+ // request completes and the message reaches the code panel.
172+ const e = err as { message?: string; file?: string; line?: number };
173+ let m: ErrorMessage = {
174+ type: 'error',
175+ id: msg.id,
176+ message: e?.message ? String(e.message) : String(err),
177+ file: typeof e?.file === 'string' ? e.file : undefined,
178+ line: typeof e?.line === 'number' ? e.line : undefined,
179+ };
180+ for (const w of warnings) {
181+ const hit = SKIP_WARNING.exec(w);
182+ if (hit) {
183+ m = { type: 'error', id: msg.id, message: 'syntax error', file: hit[1], line: Number(hit[2]) };
184+ break;
185+ }
186+ }
187+ (self as unknown as Worker).postMessage(m);
188+ } finally {
189+ console.warn = origWarn;
190+ }
191+};
192+
193+function render(msg: RenderRequest): void {
194+ const { id, params, mfiles } = msg;
94195 let quality = msg.quality;
95196
96- const solveParamsKey = keyFor(params, SOLVE_KEYS);
197+ // The .m sources are part of the cache identity: editing the code must
198+ // invalidate everything computed from the previous version.
199+ const codeKey = hashMFiles(mfiles);
200+ const solveParamsKey = codeKey + '|' + keyFor(params, SOLVE_KEYS);
97201 const fullSolveKey = solveParamsKey + `|os=${params.oversample}`;
98202 if (quality === 'preview' && solveCache.has(fullSolveKey)) {
99203 quality = 'full'; // full solve already cached; only cheap stages remain
@@ -103,67 +207,124 @@ self.onmessage = (ev: MessageEvent<RenderRequest>) => {
103207 const elements = quality === 'preview'
104208 ? Math.min(24, params.radiating_elements)
105209 : params.radiating_elements;
210+ const fsHi = params.sample_rate * oversample;
106211
107212 const post = (m: ProgressMessage) => (self as unknown as Worker).postMessage(m);
108213
214+ // ---- solve stage ----
109215 const solveKey = solveParamsKey + `|os=${oversample}`;
110216 let solveMs = 0;
111217 let solveCached = true;
112218 let solveEntry = lruGet(solveCache, solveKey);
113219 if (!solveEntry) {
114220 solveCached = false;
221+ post({ type: 'progress', id, stage: 'solve', frac: 0.15 });
115222 const t0 = performance.now();
116- const basis = buildBasis(params);
117- const solve = simulate(params, basis, oversample, (frac) =>
118- post({ type: 'progress', id, stage: 'solve', frac }),
119- );
223+ const r = runStage(mfiles, `
224+${pLiteral(params)}
225+basis = woodblock_basis(P);
226+[Qdd, Xdd, force, helmholtzHz, contactMs, peakForceN, maxRealLam] = woodblock_solve(P, basis, ${oversample});
227+f = basis.f;
228+shapes = basis.shapes(1:32:end, :);
229+`);
120230 solveMs = performance.now() - t0;
121- solveEntry = { basis, solve };
231+
232+ const fT = asTensor(r.variableValues.f);
233+ const f = Array.from(fT.data);
234+ const n = f.length;
235+
236+ // decimated mode shapes, normalized per mode for the viewer
237+ const shapesT = asTensor(r.variableValues.shapes);
238+ const rows = shapesT.shape[0];
239+ const modeShapes = new Float32Array(n * SHAPE_PTS);
240+ for (let m = 0; m < n; m++) {
241+ let peak = 0;
242+ for (let i = 0; i < rows; i++) peak = Math.max(peak, Math.abs(shapesT.data[m * rows + i]));
243+ for (let i = 0; i < SHAPE_PTS; i++) {
244+ const j = Math.round((i / (SHAPE_PTS - 1)) * (rows - 1));
245+ modeShapes[m * SHAPE_PTS + i] = peak > 0 ? shapesT.data[m * rows + j] / peak : 0;
246+ }
247+ }
248+
249+ const forceT = asTensor(r.variableValues.force);
250+ const forceLen = Math.min(forceT.data.length, Math.round((FORCE_MS / 1000) * fsHi));
251+ const force = new Float32Array(forceLen);
252+ for (let i = 0; i < forceLen; i++) force[i] = forceT.data[i];
253+
254+ solveEntry = {
255+ Qdd: r.variableValues.Qdd,
256+ Xdd: r.variableValues.Xdd,
257+ fsHi,
258+ f,
259+ modeShapes,
260+ force,
261+ helmholtzHz: asScalar(r.variableValues.helmholtzHz),
262+ contactMs: asScalar(r.variableValues.contactMs),
263+ peakForceN: asScalar(r.variableValues.peakForceN),
264+ };
122265 lruSet(solveCache, solveKey, solveEntry, SOLVE_CACHE_MAX);
123266 }
124- const { basis, solve } = solveEntry;
125267
268+ // ---- radiate stage ----
126269 const radiateKey = solveKey + '||' + keyFor(params, RADIATE_KEYS) + `|M=${elements}`;
127270 let radiateMs = 0;
128271 let radiateCached = true;
129272 let radEntry = lruGet(radiateCache, radiateKey);
130273 if (!radEntry) {
131274 radiateCached = false;
275+ post({ type: 'progress', id, stage: 'radiate', frac: 0.55 });
132276 const t0 = performance.now();
133- const rad = radiate(params, solve.Qdd, solve.Xdd, basis, solve.fsHi, elements, (frac) =>
134- post({ type: 'progress', id, stage: 'radiate', frac }),
135- );
277+ const r = runStage(mfiles, `
278+${pLiteral(params)}
279+basis = woodblock_basis(P);
280+[p, spacingNeed, listenerInside] = woodblock_radiate(P, basis, Qdd, Xdd, ${solveEntry.fsHi}, ${elements});
281+`, { Qdd: solveEntry.Qdd, Xdd: solveEntry.Xdd });
136282 radiateMs = performance.now() - t0;
137- radEntry = { p: rad.p, warnings: rad.warnings };
283+
284+ const warnings: string[] = [];
285+ const spacingNeed = asScalar(r.variableValues.spacingNeed);
286+ if (spacingNeed > 0) {
287+ const fMax = solveEntry.f[solveEntry.f.length - 1];
288+ warnings.push(
289+ `element spacing is coarse at ${fMax.toFixed(0)} Hz; consider ${spacingNeed} radiating elements`,
290+ );
291+ }
292+ if (asScalar(r.variableValues.listenerInside) !== 0) {
293+ const d = Math.hypot(params.listener_x_m, params.listener_y_m, params.listener_z_m);
294+ warnings.push(
295+ `listener is ${(d * 100).toFixed(0)} cm from the block centre, inside the source region; point-source elements are crude here`,
296+ );
297+ }
298+ radEntry = { p: r.variableValues.p, warnings };
138299 lruSet(radiateCache, radiateKey, radEntry, RADIATE_CACHE_MAX);
139300 }
140301
141- const out = outputStage(params, radEntry.p, oversample);
302+ // ---- output stage ----
303+ post({ type: 'progress', id, stage: 'output', frac: 0.9 });
304+ const rOut = runStage(mfiles, `
305+${pLiteral(params)}
306+[y, peakPa] = woodblock_output(P, p, ${oversample});
307+`, { p: radEntry.p });
142308
143- // Contact force excerpt for plotting
144- const forceLen = Math.min(solve.force.length, Math.round((FORCE_MS / 1000) * solve.fsHi));
145- const force = new Float32Array(forceLen);
146- for (let i = 0; i < forceLen; i++) force[i] = solve.force[i];
309+ const yT = asTensor(rOut.variableValues.y);
310+ const y32 = new Float32Array(yT.data.length);
311+ for (let i = 0; i < yT.data.length; i++) y32[i] = yT.data[i];
312+ const peakPa = asScalar(rOut.variableValues.peakPa);
147313
148- // Mode shapes decimated for the mode viewer (unnormalized visual scale)
149- const modeShapes = new Float32Array(basis.n * SHAPE_PTS);
150- for (let n = 0; n < basis.n; n++) {
151- const src = basis.shapes[n];
314+ const warnings = [...radEntry.warnings];
315+ if (!params.normalize) {
152316 let peak = 0;
153- for (let i = 0; i < src.length; i++) peak = Math.max(peak, Math.abs(src[i]));
154- for (let i = 0; i < SHAPE_PTS; i++) {
155- const j = Math.round((i / (SHAPE_PTS - 1)) * (src.length - 1));
156- modeShapes[n * SHAPE_PTS + i] = peak > 0 ? src[j] / peak : 0;
317+ for (let i = 0; i < y32.length; i++) peak = Math.max(peak, Math.abs(y32[i]));
318+ if (peak > 1.0) {
319+ warnings.push(
320+ `clipping: peak is ${peakPa.toFixed(2)} Pa but full_scale_pa is ${params.full_scale_pa.toFixed(2)}`,
321+ );
157322 }
158323 }
159-
160- const y32 = new Float32Array(out.y.length);
161- for (let i = 0; i < out.y.length; i++) y32[i] = out.y[i];
162-
163- const warnings = [...radEntry.warnings];
164- if (out.clipWarning) warnings.push(out.clipWarning);
165324 const lx = params.listener_x_m, ly = params.listener_y_m, lz = params.listener_z_m;
166325
326+ const force = solveEntry.force.slice();
327+ const modeShapes = solveEntry.modeShapes.slice();
167328 const resp: RenderResponse = {
168329 type: 'result',
169330 id,
@@ -171,18 +332,18 @@ self.onmessage = (ev: MessageEvent<RenderRequest>) => {
171332 y: y32,
172333 sampleRate: params.sample_rate,
173334 force,
174- forceFsHz: solve.fsHi,
335+ forceFsHz: solveEntry.fsHi,
175336 modeShapes,
176337 shapePts: SHAPE_PTS,
177338 diag: {
178- modesKept: basis.n,
179- partialsHz: Array.from(basis.f),
180- ratios: Array.from(basis.f).map((f) => f / basis.f[0]),
181- helmholtzHz: solve.helmholtzHz,
182- contactMs: solve.contactMs,
183- peakForceN: solve.peakForceN,
184- peakPa: out.peakPa,
185- peakSplDb: out.peakSplDb,
339+ modesKept: solveEntry.f.length,
340+ partialsHz: solveEntry.f,
341+ ratios: solveEntry.f.map((f) => f / solveEntry.f[0]),
342+ helmholtzHz: solveEntry.helmholtzHz,
343+ contactMs: solveEntry.contactMs,
344+ peakForceN: solveEntry.peakForceN,
345+ peakPa,
346+ peakSplDb: 20 * Math.log10(peakPa / 20e-6 + 1e-30),
186347 listenerDistM: Math.sqrt(lx * lx + ly * ly + lz * lz),
187348 warnings,
188349 solveMs,
@@ -193,9 +354,5 @@ self.onmessage = (ev: MessageEvent<RenderRequest>) => {
193354 elements,
194355 },
195356 };
196- (self as unknown as Worker).postMessage(resp, [
197- y32.buffer,
198- force.buffer,
199- modeShapes.buffer,
200- ]);
201-};
357+ (self as unknown as Worker).postMessage(resp, [y32.buffer, force.buffer, modeShapes.buffer]);
358+}
webapp/tsconfig.jsonmodified+2−1View file
@@ -15,5 +15,6 @@
1515 "noUnusedParameters": true,
1616 "noFallthroughCasesInSwitch": true
1717 },
18- "include": ["src"]
18+ "include": ["src"],
19+ "exclude": ["src/selftest.ts"]
1920 }
webapp/vite.config.tsmodified+1−0View file
@@ -5,4 +5,5 @@ export default defineConfig({
55 plugins: [react()],
66 base: './',
77 worker: { format: 'es' },
8+ server: { fs: { allow: ['..'] } }, // the .m sources live in ../matlab
89 });