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
64% constant and the free-sphere scaling must emerge.
65Q = P;
66Q.beater_mass_kg = 5e-4;
67[~, ~, ~, tau1, ~, ~, ~] = triangle_strikes(Q);
68Q.strike_velocity = 2 * P.strike_velocity;
69[~, ~, ~, tau2, ~, ~, ~] = triangle_strikes(Q);
70ratio = tau2(1) / tau1(1);
71ok = abs(ratio - 2^(-1 / 5)) < 0.03;
72fails = fails + (~ok);
73fprintf('hertz scaling: tau ratio %.4f (target %.4f) %s\n', ratio, 2^(-1 / 5), tv_verdict(ok));
75% ---- 4. seeded roll is exactly reproducible ----
76Q = P;
77Q.strike_count = 12;
78[t1r, v1r, j1r, ~, ~, s1r, g1] = triangle_strikes(Q);
79[t2r, v2r, j2r, ~, ~, s2r, g2] = triangle_strikes(Q);
80same = max(abs(t1r - t2r)) == 0 && max(abs(v1r - v2r)) == 0 && max(abs(s1r - s2r)) == 0;
81ok = same ...
82 && g1.count == 12 ...
83 && abs(g1.meanV - 1.501211020) < 1e-6 ...
84 && abs(sum(j1r) - 0.383069278968) < 1e-9 ...
85 && abs(t1r(2) - 0.087732781124) < 1e-9;
86fails = fails + (~ok);
87fprintf('roll (seed %d): %d strikes, mean v %.6f m/s, sum J %.9f N s %s\n', ...
88 Q.seed, g1.count, g1.meanV, sum(j1r), tv_verdict(ok));
90% ---- 5. far-field falloff: -6 dB per doubling of distance ----
91% RMS over the render, with the two listeners on the same ray from the
92% origin so every element dipole keeps its directivity factor: off the ray
93% the in-plane dipoles fall faster than 1/r, which is physics rather than
94% error. Modes are capped at 8 kHz so 2 m is beyond the Rayleigh distance
95% of the frame; at 20 kHz it is not, and Fresnel-zone interference shifts
96% the ratio by a few tenths of a dB, again physics rather than error.
97Q = P;
98Q.normalize = 0;
99Q.tail_s = 1.0;
100Q.max_mode_hz = 8000;
101Q.listener_x_m = 0.6; Q.listener_y_m = 0.5; Q.listener_z_m = 1.8;
102[y1, ~, ~] = triangle_render(Q);
103Q.listener_x_m = 1.2; Q.listener_y_m = 1.0; Q.listener_z_m = 3.6;
104[y2, ~, ~] = triangle_render(Q);
105drop = 10 * log10(sum(y2.^2) / sum(y1.^2));
106ok = abs(drop + 6.02) < 0.2;
107fails = fails + (~ok);
108fprintf('falloff 2 m -> 4 m: %.2f dB rms (target -6.02) %s\n', drop, tv_verdict(ok));
110% ---- 6. pinned end-to-end diagnostics at the defaults ----
111[y, basis, dinfo] = triangle_render(P);
112ok = abs(dinfo.peakPa - 7.217688) < 0.02 ...
113 && abs(basis.f(1) - 163.9200) < 0.1 ...
114 && basis.n == 40;
115fails = fails + (~ok);
116fprintf('pinned render: peak %.4f Pa (target %.4f), mode1 %.2f Hz (target %.2f), %d modes %s\n', ...
117 dinfo.peakPa, 7.217688, basis.f(1), 163.9200, basis.n, tv_verdict(ok));
119% ---- 7. the ring decays ----
120% A triangle rings for seconds, so the tail is not silent; but the energy in
121% the last half second must sit well below the first half second's.
122fs = P.sample_rate;
123n = length(y);
124head = y(1:round(0.5 * fs));
125tail = y(n - round(0.5 * fs):n);
126dropE = 10 * log10(sum(tail.^2) / sum(head.^2));
127ok = dropE < -25;
128fails = fails + (~ok);
129fprintf('ring decay: last/first half second %.1f dB %s\n', dropE, tv_verdict(ok));
131if fails > 0
132 fprintf('%d CHECK(S) FAILED\n', fails);
133else
134 fprintf('all checks passed\n');
135end
137function f = tv_eigfreqs(nodes, family, E, rho, nu, d)
138[K, M] = triangle_fem(nodes, family, E, rho, nu, d);
139Rc = chol(M);
140Am = (Rc') \ (K / Rc);
141Am = (Am + Am') / 2;
142[~, Sv, ~] = svd(Am);
143f = sort(sqrt(diag(Sv))) / (2 * pi);
144end
146function s = tv_verdict(ok)
147if ok
148 s = 'ok';
149else
150 s = 'MISMATCH';
151end
152end