function [M, Kmat, omS, zS, meta] = bass_assemble(P, s, st)
% BASS_ASSEMBLE Mass and stiffness matrices of one fingering's coupled
% string-body system, plus the uncoupled string modes and damping.
%
% The displacement ansatz and the resulting matrices are documented in
% bass_basis.m, which calls this per fingering; bass_validate.m calls it
% directly to check the closed-form string limit and the eigensolve
% identities. Coordinates are [q_1..q_ns, xi_1..xi_nb]: mass-normalized
% pinned-pinned stiff-string modes, then the mass-normalized body modes of
% bass_body.m. M and Kmat are symmetric and M is positive definite (the
% coupling is weak: the off-diagonal mass terms are of order
% sqrt(mu Lv / m_body) / k pi, a few percent here).

f0open = [97.9989 73.4162 55.0000 41.2034] * P.tuning_trim;
tension = [P.tension_g_n P.tension_d_n P.tension_a_n P.tension_e_n];
L0 = P.scale_length_m;

body = bass_body(P);
nb = body.n;
omRes = 2 * pi * 1.25 * max(body.f);   % residual admittance high-pass corner

T = tension(s);
cw = 2 * L0 * f0open(s);               % transverse wave speed of the string
mu = T / cw^2;
Lv = L0 * 2^(-st / 12);
EI = P.inharmonicity * T * L0^2 / pi^2;

% string modes up to max_mode_hz (at least one, at most 160)
kEst = min(160, max(1, floor(2 * Lv * P.max_mode_hz / cw) + 2));
kk = (1:kEst)' * pi / Lv;
om2 = kk.^2 * (T / mu) + kk.^4 * (EI / mu);
omAll = sqrt(om2);
ns = max(1, sum(omAll <= 2 * pi * P.max_mode_hz));
ns = min(ns, kEst);
omS = omAll(1:ns);
kv = (1:ns)';

% uncoupled string damping: internal loss, the finger termination on
% stopped notes, and the resistive residual bridge admittance (high-passed
% above the retained body modes so the modal region is not double-counted;
% a resistive end gives a constant decay RATE, zeta ~ Y T / (Lv om))
hp = omS.^2 ./ (omS.^2 + omRes^2);
zS = 0.5 * (P.loss_factor + P.loss_slope * omS) ...
  + 0.5 * P.residual_y * (T / Lv) * hp ./ omS;
if st > 0
  zS = zS + 0.5 * P.finger_loss;
end

% ---- assemble ----
n = ns + nb;
gk = sqrt(2 * mu * Lv) * ((-1).^(kv + 1)) ./ (kv * pi);
M = eye(n);
Kmat = zeros(n, n);
for k = 1:ns
  Kmat(k, k) = omS(k)^2;
end
for j = 1:nb
  for k = 1:ns
    M(k, ns + j) = gk(k) * body.a(j);
    M(ns + j, k) = M(k, ns + j);
  end
  for i = 1:nb
    M(ns + i, ns + j) = M(ns + i, ns + j) + (mu * Lv / 3) * body.a(i) * body.a(j);
    Kmat(ns + i, ns + j) = (T / Lv) * body.a(i) * body.a(j);
  end
  Kmat(ns + j, ns + j) = Kmat(ns + j, ns + j) + body.omega(j)^2;
end

meta.ns = ns;
meta.nb = nb;
meta.n = n;
meta.Lv = Lv;
meta.mu = mu;
meta.T = T;
meta.f0 = f0open(s) * 2^(st / 12);
meta.body = body;
end