1function [tstr, vstr, jstr, taustr, fstr, sstr, sinfo] = cymbal_strikes(P, basis)
2% CYMBAL_STRIKES Stick contacts: a single strike, or a seeded repeated figure.
3%
4% strike_count = 1 gives one deterministic strike of strike_velocity at
5% strike_r at t = 0; no random numbers are drawn. strike_count > 1 repeats it
6% at strike_rate_hz, alternating between strike_r and roll_r2 (a ride pattern
7% moving between the bow and the edge), with Gaussian timing jitter
8% (timing_jitter_ms) and lognormal velocity spread (velocity_spread). The
9% first hit stays at t = 0 so the onset is fixed. All randomness comes from
10% cymbal_rng seeded with P.seed, so a given seed reproduces the identical
11% pattern in the CLI, in Node, and in the browser. Every contact sits on the
12% meridian theta = 0; see cymbal_basis.m for why.
13%
14% Each contact is a Hertz impact of the stick tip on the plate:
15% E* = 1 / ((1 - nu_t^2) / E_t + (1 - nu_p^2) / E_p)
16% R_eff= 1 / sqrt((1/r_tip + 1/R_s)(1/r_tip + 1/R_theta)) (sphere on a
17% doubly curved surface, at the local principal curvatures)
18% k_h = (4/3) E* sqrt(R_eff)
19% tau = 3.218 (m*^2 / (k_h^2 v))^{1/5} (Hertz contact time)
20% J = (1 + restitution) m* v
21% Fpk = pi J / (2 tau) (half-sine peak)
22%
23% The reduced mass m* pairs the stick with the plate mass engaged during the
24% contact. Here the sibling triangle's flexural-group-velocity estimate is
25% replaced by something exact: the driving-point impedance of an infinite
26% plate is a pure resistance, Z = 8 sqrt(D rho t), independent of frequency,
27% so the mass the stick sees during a pulse of duration tau is simply Z tau. A
28% short fixed-point iteration closes the loop, since tau depends on m* which
29% depends on tau. What this captures is that a heavy stick, a soft tip or a
30% slow stroke all lengthen the contact and darken the strike; what it leaves
31% out is the millisecond detail of the contact itself, which is prescribed
32% rather than simulated.
33%
34% Returns row vectors sorted by time: tstr (s), vstr (m/s), jstr (N s),
35% taustr (s), fstr (peak N), sstr (arclength fraction of the contact), and
36% sinfo with summary statistics.
38K = max(1, round(P.strike_count));
39rate = max(1e-3, P.strike_rate_hz);
40duration = (K - 1) / rate + P.tail_s;
42E = P.youngs_modulus;
43nu = P.poisson_ratio;
44rho = P.density;
45nut = 0.35; % stick tip Poisson ratio, fixed
46Estar = 1 / ((1 - nut^2) / P.stick_youngs + (1 - nu^2) / E);
47rt = max(1e-4, P.tip_radius_m);
49state = mod(floor(abs(P.seed)), 4294967296);
50z1 = zeros(K, 1);
51z2 = zeros(K, 1);
52if K > 1
53 [u, state] = cymbal_rng(state, 4 * K);
54 U = reshape(u, K, 4);
55 z1 = sqrt(-2 * log(U(:, 1))) .* cos(2 * pi * U(:, 2));
56 z2 = sqrt(-2 * log(U(:, 3))) .* cos(2 * pi * U(:, 4));
57end
59tstr = zeros(1, K);
60vstr = zeros(1, K);
61jstr = zeros(1, K);
62taustr = zeros(1, K);
63fstr = zeros(1, K);
64sstr = zeros(1, K);
65for i = 1:K
66 t = (i - 1) / rate;
67 if i > 1
68 t = max(0, t + P.timing_jitter_ms * 1e-3 * z1(i));
69 end
70 v = P.strike_velocity;
71 if K > 1
72 v = v * exp(P.velocity_spread * z2(i));
73 end
74 v = max(0.01, v);
75 if mod(i, 2) == 1
76 sfrac = P.strike_r;
77 else
78 sfrac = P.roll_r2;
79 end
80 sfrac = min(0.999, max(0.001, sfrac));
82 % local plate properties at the contact
83 [~, ~, thk, ~, ~, cs, ct] = cymbal_at(basis, sfrac);
84 Rrel = 1 / sqrt(max(1e-9, (1 / rt + cs) * (1 / rt + ct)));
85 kh = (4 / 3) * Estar * sqrt(Rrel);
86 D = E * thk^3 / (12 * (1 - nu^2));
87 Zp = 8 * sqrt(D * rho * thk); % driving-point resistance of the plate
89 tau = 1e-4;
90 mstar = P.stick_mass_kg;
91 for it = 1:6
92 mloc = Zp * tau;
93 mstar = P.stick_mass_kg * mloc / (P.stick_mass_kg + mloc);
94 tau = 3.218 * (mstar^2 / (kh^2 * v))^(1 / 5);
95 end
97 tstr(i) = t;
98 vstr(i) = v;
99 jstr(i) = (1 + P.restitution) * mstar * v;
100 taustr(i) = tau;
101 fstr(i) = pi * jstr(i) / (2 * tau);
102 sstr(i) = sfrac;
103end
105[tstr, order] = sort(tstr);
106tstr = reshape(tstr, 1, K);
107vstr = reshape(vstr(order), 1, K);
108jstr = reshape(jstr(order), 1, K);
109taustr = reshape(taustr(order), 1, K);
110fstr = reshape(fstr(order), 1, K);
111sstr = reshape(sstr(order), 1, K);
113sinfo.count = K;
114sinfo.meanV = sum(vstr) / K;
115sinfo.maxV = max(vstr);
116sinfo.meanTauUs = sum(taustr) / K * 1e6;
117sinfo.peakForceN = max(fstr);
118sinfo.durationS = duration;
119end