96aec9cTriangle strike simulator: MATLAB model and interactive web appJeremy Magland 1% Validate the triangle model. There is no external reference implementation,
2% so the checks are (a) closed-form identities the physics must satisfy,
3% (b) structural identities of the FEM eigensolve, (c) exact reproducibility
4% of the seeded roll process, and (d) pinned end-to-end diagnostics from the
5% default configuration, recorded when the model was built and treated as a
6% regression reference since.
7% Run with: numbl run triangle_validate.m
9fails = 0;
11% ---- 1. straight free-free rod against closed forms ----
12% The same element routines that build the triangle, assembled on a straight
13% polyline, must reproduce the textbook free-free bar: bending frequencies
14% (beta L)^2 kappa cL / (2 pi L^2), the first axial mode cL / 2L, the first
15% torsional mode cT / 2L, and exactly three rigid-body modes per family.
16E = 2e11; rho = 7850; nu = 0.30; d = 0.008; L = 0.44;
17nodes = [linspace(0, L, 61); zeros(1, 61)];
18cL = sqrt(E / rho);
19kap = d / 4;
20blTab = [4.73004074 7.85320462 10.99560784];
21fbClosed = blTab.^2 * kap * cL / (2 * pi * L^2);
22faClosed = cL / (2 * L);
23ftClosed = sqrt(E / (2 * (1 + nu)) / rho) / (2 * L);
24f0 = tv_eigfreqs(nodes, 0, E, rho, nu, d);
25f1 = tv_eigfreqs(nodes, 1, E, rho, nu, d);
26% the rigid-body zeros come back as numerical noise well below the first
27% elastic mode (185 Hz here); 5 Hz separates them cleanly
28nRigid0 = sum(f0 < 5);
29nRigid1 = sum(f1 < 5);
30e0 = f0(f0 > 1);
31e1 = f1(f1 > 1);
32errB = max(abs(e0(1:3).' - fbClosed) ./ fbClosed);
33% rotary inertia lowers the FEM bending values slightly below the closed
34% form, which neglects it; 0.5% covers that plus discretization
35errA = abs(min(abs(e0 - faClosed)) / faClosed);
36errT = abs(min(abs(e1 - ftClosed)) / ftClosed);
37ok = nRigid0 == 3 && nRigid1 == 3 && errB < 5e-3 && errA < 2e-3 && errT < 2e-3;
38fails = fails + (~ok);
39fprintf('straight rod: bend err %.1e, axial %.1e, torsion %.1e, rigid %d/%d %s\n', ...
40 errB, errA, errT, nRigid0, nRigid1, tv_verdict(ok));
42% ---- 2. eigensolve identities on the triangle itself ----
43% The Cholesky reduction must hand back mass-orthonormal shapes that
44% diagonalize the stiffness: phi' M phi = I and phi' K phi = diag(om^2).
45P = triangle_defaults();
46geo = triangle_geometry(P);
47[K, M] = triangle_fem(geo.nodes, 0, P.youngs_modulus, P.density, P.poisson_ratio, P.rod_diameter_m);
48Rc = chol(M);
49Am = (Rc') \ (K / Rc);
50Am = (Am + Am') / 2;
51[Uv, Sv, ~] = svd(Am);
52om2 = flipud(diag(Sv));
53Vv = fliplr(Uv);
54Phi = Rc \ Vv;
55eOrth = max(max(abs(Phi' * M * Phi - eye(length(om2)))));
56KP = Phi' * K * Phi;
57eDiag = max(max(abs(KP - diag(om2)))) / max(om2);
58ok = eOrth < 1e-8 && eDiag < 1e-10;
59fails = fails + (~ok);
60fprintf('eigensolve: orthonormality %.1e, diagonalization %.1e %s\n', eOrth, eDiag, tv_verdict(ok));
62% ---- 3. Hertz contact-time scaling tau ~ v^(-1/5) ----
63% With a beater much lighter than the engaged rod mass the reduced mass is a
59c69d1Default to a twelve-strike rollJeremy Magland 64% constant and the free-sphere scaling must emerge. A single strike, so the
65% check does not depend on the default gesture or on the velocity spread.
96aec9cTriangle strike simulator: MATLAB model and interactive web appJeremy Magland 68Q.beater_mass_kg = 5e-4;
69[~, ~, ~, tau1, ~, ~, ~] = triangle_strikes(Q);
70Q.strike_velocity = 2 * P.strike_velocity;
71[~, ~, ~, tau2, ~, ~, ~] = triangle_strikes(Q);
72ratio = tau2(1) / tau1(1);
73ok = abs(ratio - 2^(-1 / 5)) < 0.03;
74fails = fails + (~ok);
75fprintf('hertz scaling: tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), tv_verdict(ok));
77% ---- 4. seeded roll is exactly reproducible ----
78Q = P;
79Q.strike_count = 12;
80[t1r, v1r, j1r, ~, ~, s1r, g1] = triangle_strikes(Q);
81[t2r, v2r, j2r, ~, ~, s2r, g2] = triangle_strikes(Q);
82same = max(abs(t1r - t2r)) == 0 && max(abs(v1r - v2r)) == 0 && max(abs(s1r - s2r)) == 0;
83ok = same ...
84 && g1.count == 12 ...
85 && abs(g1.meanV - 1.501211020) < 1e-6 ...
86 && abs(sum(j1r) - 0.383069278968) < 1e-9 ...
87 && abs(t1r(2) - 0.087732781124) < 1e-9;
88fails = fails + (~ok);
89fprintf('roll (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\n', ...
90 Q.seed, g1.count, g1.meanV, sum(j1r), tv_verdict(ok));
92% ---- 5. far-field falloff: -6 dB per doubling of distance ----
93% RMS over the render, with the two listeners on the same ray from the
94% origin so every element dipole keeps its directivity factor: off the ray
95% the in-plane dipoles fall faster than 1/r, which is physics rather than
96% error. Modes are capped at 8 kHz so 2 m is beyond the Rayleigh distance
97% of the frame; at 20 kHz it is not, and Fresnel-zone interference shifts
98% the ratio by a few tenths of a dB, again physics rather than error.
99Q = P;
100Q.normalize = 0;
96aec9cTriangle strike simulator: MATLAB model and interactive web appJeremy Magland 102Q.tail_s = 1.0;
103Q.max_mode_hz = 8000;
104Q.listener_x_m = 0.6; Q.listener_y_m = 0.5; Q.listener_z_m = 1.8;
105[y1, ~, ~] = triangle_render(Q);
106Q.listener_x_m = 1.2; Q.listener_y_m = 1.0; Q.listener_z_m = 3.6;
107[y2, ~, ~] = triangle_render(Q);
108drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
109ok = abs(drop + 6.02) < 0.2;
110fails = fails + (~ok);
111fprintf('falloff 2 m -> 4 m: %.2f dB rms (target -6.02) %s\n', drop, tv_verdict(ok));
113% ---- 6. pinned end-to-end diagnostics at the defaults ----
114[y, basis, dinfo] = triangle_render(P);
59c69d1Default to a twelve-strike rollJeremy Magland 115ok = abs(dinfo.peakPa - 13.769445) < 0.02 ...
96aec9cTriangle strike simulator: MATLAB model and interactive web appJeremy Magland 116 && abs(basis.f(1) - 163.9200) < 0.1 ...
117 && basis.n == 40;
118fails = fails + (~ok);
119fprintf('pinned render: peak %.4f Pa (target %.4f), mode1 %.2f Hz (target %.2f), %d modes %s\n', ...
59c69d1Default to a twelve-strike rollJeremy Magland 120 dinfo.peakPa, 13.769445, basis.f(1), 163.9200, basis.n, tv_verdict(ok));
122% ---- 7. the ring decays ----
123% A triangle rings for seconds, so the tail is not silent; but the energy in
59c69d1Default to a twelve-strike rollJeremy Magland 124% the last half second must sit well below that of the opening half second,
125% which at the defaults carries the first six hits of the roll.
96aec9cTriangle strike simulator: MATLAB model and interactive web appJeremy Magland 126fs = P.sample_rate;
127n = length(y);
128head = y(1:round(0.5 * fs));
129tail = y(n - round(0.5 * fs):n);
130dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
131ok = dropE < -25;
132fails = fails + (~ok);
133fprintf('ring decay: last/first half second %.1f dB %s\n', dropE, tv_verdict(ok));
135if fails > 0
136 fprintf('%d CHECK(S) FAILED\n', fails);
137else
138 fprintf('all checks passed\n');
139end
141function f = tv_eigfreqs(nodes, family, E, rho, nu, d)
142[K, M] = triangle_fem(nodes, family, E, rho, nu, d);
143Rc = chol(M);
144Am = (Rc') \ (K / Rc);
145Am = (Am + Am') / 2;
146[~, Sv, ~] = svd(Am);
147f = sort(sqrt(diag(Sv))) / (2 * pi);
148end
150function s = tv_verdict(ok)
151if ok
152 s = 'ok';
153else
154 s = 'MISMATCH';
155end
156end