96aec9cTriangle strike simulator: MATLAB model and interactive web appJeremy Magland 1function [tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = triangle_strikes(P)
2% TRIANGLE_STRIKES Beater strike events: a single strike, or a seeded roll.
3%
4% strike_count = 1 gives one deterministic strike of strike_velocity at
5% strike_pos at t = 0; no random numbers are drawn. strike_count > 1 gives a
6% roll: nominally periodic hits at strike_rate_hz, alternating between
7% strike_pos and roll_pos2 (the beater rattling in a corner strikes the two
8% legs that form it in turn), with Gaussian timing jitter (timing_jitter_ms)
9% and lognormal velocity spread (velocity_spread). The first hit stays at
10% t = 0 so the onset is fixed. All randomness comes from triangle_rng seeded
11% with P.seed, so a given seed reproduces the identical roll in the CLI,
12% Node, and the browser.
13%
14% Each strike is a Hertz contact of the beater tip (a steel sphere of radius
15% beater_radius_m) on the rod:
16% k_h = (4/3) E* sqrt(R_eff), R_eff = sqrt(r_b R''),
17% 1/R'' = 1/r_b + 2/d (sphere on cylinder: R_eff is the
18% geometric mean of the two relative
19% principal curvatures)
20% tau = 3.218 (m*^2 / (k_h^2 v))^{1/5} (Hertz contact time)
21% J = (1 + restitution) m* v (momentum transferred)
22% Fpk = pi J / (2 tau) (half-sine pulse peak)
23% The reduced mass m* pairs the beater with the rod mass engaged during the
24% contact. The engaged length is estimated from the flexural group velocity
25% at the pulse's own bandwidth, l = 2 sqrt(pi kappa cL tau), which depends
26% on tau, so a short fixed-point iteration closes the loop. This is an
27% estimate, not a contact simulation: what it captures is that a heavy
28% beater lengthens the contact and darkens the strike, while a light beater
29% approaches the free-sphere limit and tau ~ v^(-1/5).
30%
31% Returns row vectors sorted by time: tstr (s), vstr (m/s), jstr (N s),
32% taustr (s), fstr (peak N), sstr (arclength fraction of the contact), and
33% sinfo with summary statistics.
35K = max(1, round(P.strike_count));
36rate = max(1e-3, P.strike_rate_hz);
37duration = (K - 1) / rate + P.tail_s;
39E = P.youngs_modulus;
40rho = P.density;
41d = P.rod_diameter_m;
42nub = 0.30; % beater Poisson ratio, fixed
43Estar = 1 / ((1 - nub^2) / P.beater_youngs + (1 - P.poisson_ratio^2) / E);
44rb = P.beater_radius_m;
45Rpp = 1 / (1 / rb + 2 / d);
46Reff = sqrt(rb * Rpp);
47kh = (4 / 3) * Estar * sqrt(Reff);
48kapcL = (d / 4) * sqrt(E / rho); % kappa * cL: flexural dispersion constant
49Aarea = pi * d^2 / 4;
51state = mod(floor(abs(P.seed)), 4294967296);
52z1 = zeros(K, 1);
53z2 = zeros(K, 1);
54if K > 1
55 [u, state] = triangle_rng(state, 4 * K);
56 U = reshape(u, K, 4);
57 z1 = sqrt(-2 * log(U(:, 1))) .* cos(2 * pi * U(:, 2));
58 z2 = sqrt(-2 * log(U(:, 3))) .* cos(2 * pi * U(:, 4));
59end
61tstr = zeros(1, K);
62vstr = zeros(1, K);
63jstr = zeros(1, K);
64taustr = zeros(1, K);
65fstr = zeros(1, K);
66sstr = zeros(1, K);
67for i = 1:K
68 t = (i - 1) / rate;
69 if i > 1
70 t = max(0, t + P.timing_jitter_ms * 1e-3 * z1(i));
71 end
72 v = P.strike_velocity;
73 if K > 1
74 v = v * exp(P.velocity_spread * z2(i));
75 end
76 v = max(0.01, v);
78 % engaged-mass fixed point for the Hertz contact time
79 tau = 30e-6;
80 mstar = P.beater_mass_kg;
81 for it = 1:4
82 leng = 2 * sqrt(pi * kapcL * tau);
83 mloc = rho * Aarea * leng;
84 mstar = P.beater_mass_kg * mloc / (P.beater_mass_kg + mloc);
85 tau = 3.218 * (mstar^2 / (kh^2 * v))^(1 / 5);
86 end
88 tstr(i) = t;
89 vstr(i) = v;
90 jstr(i) = (1 + P.restitution) * mstar * v;
91 taustr(i) = tau;
92 fstr(i) = pi * jstr(i) / (2 * tau);
93 if mod(i, 2) == 1
94 sstr(i) = P.strike_pos;
95 else
96 sstr(i) = P.roll_pos2;
97 end
98end
100[tstr, order] = sort(tstr);
101tstr = reshape(tstr, 1, K);
102vstr = reshape(vstr(order), 1, K);
103jstr = reshape(jstr(order), 1, K);
104taustr = reshape(taustr(order), 1, K);
105fstr = reshape(fstr(order), 1, K);
106sstr = reshape(sstr(order), 1, K);
108sinfo.count = K;
109sinfo.meanV = sum(vstr) / K;
110sinfo.maxV = max(vstr);
111sinfo.meanTauUs = sum(taustr) / K * 1e6;
112sinfo.peakForceN = max(fstr);
113sinfo.durationS = duration;
114end