function MB = cymbal2_mallet(P)
% CYMBAL2_MALLET Precompute the striking mechanism: the mallet stick as an
% Euler-Bernoulli beam (Hermite cubic elements, consistent mass,
% mass-proportional damping), held at one end by the hand, with the mallet
% head attached at the other end as a separate SDOF body (Samejima 2021,
% Sec. 3). The hand grips either rigidly (grip_clamped = 1: transverse
% displacement and rotation both prescribed) or through a rotational
% spring (grip_clamped = 0: displacement prescribed, rotation free against
% grip_spring). The head couples to the beam tip by the rigid connection
% u_M = u_S. The paper integrates the head as its own SDOF (Eq. 34) and
% solves the connection force from the constraint; because that couples an
% explicit head update to the implicit beam, the discrete pair is
% unstable, so the head mass is lumped onto the tip node instead, which is
% the same physics under the rigid connection (substituting u_M = u_S and
% f_S = -f_M2 into Eq. (31) puts M_M on the tip and f_M1 on the beam).

ne = P.n_beam;
nn = ne + 1;
ndof = 2 * nn;
le = P.stick_length_m / ne;
EI = P.stick_youngs * P.stick_inertia_m4;
rA = P.stick_density * P.stick_area_m2;

Ke = EI / le^3 * [12, 6 * le, -12, 6 * le; ...
  6 * le, 4 * le^2, -6 * le, 2 * le^2; ...
  -12, -6 * le, 12, -6 * le; ...
  6 * le, 2 * le^2, -6 * le, 4 * le^2];
Me = rA * le / 420 * [156, 22 * le, 54, -13 * le; ...
  22 * le, 4 * le^2, 13 * le, -3 * le^2; ...
  54, 13 * le, 156, -22 * le; ...
  -13 * le, -3 * le^2, -22 * le, 4 * le^2];

K = zeros(ndof, ndof);
M = zeros(ndof, ndof);
for e = 1:ne
  ix = (2 * e - 1):(2 * e + 2);
  K(ix, ix) = K(ix, ix) + Ke;
  M(ix, ix) = M(ix, ix) + Me;
end
R = P.stick_damping_rate * M;
M(ndof - 1, ndof - 1) = M(ndof - 1, ndof - 1) + P.head_mass_kg;

% Node 1 is the grip, node nn the tip. Dof 1 (grip displacement) is always
% prescribed to the hand motion; dof 2 (grip rotation) is prescribed to 0
% for the clamped hand, or free with the rotational spring added.
if P.grip_clamped ~= 0
  fi = 3:ndof;
else
  K(2, 2) = K(2, 2) + P.grip_spring;
  fi = 2:ndof;
end
Kf = K(fi, fi); Mf = M(fi, fi); Rf = R(fi, fi);
Kc = K(fi, 1); Mc = M(fi, 1); Rc = R(fi, 1);

ht = 1 / (P.sample_rate * P.oversample);
ps = P.phi_param;
% Stiffness split with the implicit weight ps/2 on the (k+1, k-1) average,
% as in the shell's washer row; the complementary split printed in the
% paper is unstable for the beam's highest elements at audio rates.
AS = Mf + (ht / 2) * Rf + (ps / 2) * ht^2 * Kf;
BS = 2 * Mf - (1 - ps) * ht^2 * Kf;
CS = -Mf + (ht / 2) * Rf - (ps / 2) * ht^2 * Kf;
ASi = inv(AS);
MB.ASiB = ASi * BS;
MB.ASiC = ASi * CS;

% Hand coupling, folded so the per-step load is
%   vM * hA(k) + vR * hB(k) + vK * hC(k)
% with hA = -(uH+ - 2 uH + uH-), hB = -(ht/2)(uH+ - uH-),
% hC = -ht^2 (ps uH + ((1-ps)/2)(uH+ + uH-)).
MB.vM = ASi * Mc;
MB.vR = ASi * Rc;
MB.vK = ASi * Kc;

tipi = find(fi == ndof - 1);
MB.tipi = tipi;
MB.tipcol = ASi(:, tipi);
MB.gtip = ASi(tipi, tipi);
MB.nfree = length(fi);
MB.ht = ht;
MB.ps = ps;
end