1function MB = cymbal2_mallet(P)
2% CYMBAL2_MALLET Precompute the striking mechanism: the mallet stick as an
3% Euler-Bernoulli beam (Hermite cubic elements, consistent mass,
4% mass-proportional damping), held at one end by the hand, with the mallet
5% head attached at the other end as a separate SDOF body (Samejima 2021,
6% Sec. 3). The hand grips either rigidly (grip_clamped = 1: transverse
7% displacement and rotation both prescribed) or through a rotational
8% spring (grip_clamped = 0: displacement prescribed, rotation free against
9% grip_spring). The head couples to the beam tip by the rigid connection
10% u_M = u_S. The paper integrates the head as its own SDOF (Eq. 34) and
11% solves the connection force from the constraint; because that couples an
12% explicit head update to the implicit beam, the discrete pair is
13% unstable, so the head mass is lumped onto the tip node instead, which is
14% the same physics under the rigid connection (substituting u_M = u_S and
15% f_S = -f_M2 into Eq. (31) puts M_M on the tip and f_M1 on the beam).
17ne = P.n_beam;
18nn = ne + 1;
19ndof = 2 * nn;
20le = P.stick_length_m / ne;
21EI = P.stick_youngs * P.stick_inertia_m4;
22rA = P.stick_density * P.stick_area_m2;
24Ke = EI / le^3 * [12, 6 * le, -12, 6 * le; ...
25 6 * le, 4 * le^2, -6 * le, 2 * le^2; ...
26 -12, -6 * le, 12, -6 * le; ...
27 6 * le, 2 * le^2, -6 * le, 4 * le^2];
28Me = rA * le / 420 * [156, 22 * le, 54, -13 * le; ...
29 22 * le, 4 * le^2, 13 * le, -3 * le^2; ...
30 54, 13 * le, 156, -22 * le; ...
31 -13 * le, -3 * le^2, -22 * le, 4 * le^2];
33K = zeros(ndof, ndof);
34M = zeros(ndof, ndof);
35for e = 1:ne
36 ix = (2 * e - 1):(2 * e + 2);
37 K(ix, ix) = K(ix, ix) + Ke;
38 M(ix, ix) = M(ix, ix) + Me;
39end
40R = P.stick_damping_rate * M;
41M(ndof - 1, ndof - 1) = M(ndof - 1, ndof - 1) + P.head_mass_kg;
43% Node 1 is the grip, node nn the tip. Dof 1 (grip displacement) is always
44% prescribed to the hand motion; dof 2 (grip rotation) is prescribed to 0
45% for the clamped hand, or free with the rotational spring added.
46if P.grip_clamped ~= 0
47 fi = 3:ndof;
48else
49 K(2, 2) = K(2, 2) + P.grip_spring;
50 fi = 2:ndof;
51end
52Kf = K(fi, fi); Mf = M(fi, fi); Rf = R(fi, fi);
53Kc = K(fi, 1); Mc = M(fi, 1); Rc = R(fi, 1);
55ht = 1 / (P.sample_rate * P.oversample);
56ps = P.phi_param;
57% Stiffness split with the implicit weight ps/2 on the (k+1, k-1) average,
58% as in the shell's washer row; the complementary split printed in the
59% paper is unstable for the beam's highest elements at audio rates.
60AS = Mf + (ht / 2) * Rf + (ps / 2) * ht^2 * Kf;
61BS = 2 * Mf - (1 - ps) * ht^2 * Kf;
62CS = -Mf + (ht / 2) * Rf - (ps / 2) * ht^2 * Kf;
63ASi = inv(AS);
64MB.ASiB = ASi * BS;
65MB.ASiC = ASi * CS;
67% Hand coupling, folded so the per-step load is
68% vM * hA(k) + vR * hB(k) + vK * hC(k)
69% with hA = -(uH+ - 2 uH + uH-), hB = -(ht/2)(uH+ - uH-),
70% hC = -ht^2 (ps uH + ((1-ps)/2)(uH+ + uH-)).
71MB.vM = ASi * Mc;
72MB.vR = ASi * Rc;
73MB.vK = ASi * Kc;
75tipi = find(fi == ndof - 1);
76MB.tipi = tipi;
77MB.tipcol = ASi(:, tipi);
78MB.gtip = ASi(tipi, tipi);
79MB.nfree = length(fi);
80MB.ht = ht;
81MB.ps = ps;
82end