1function [M, Kmat, omS, zS, meta] = bass_assemble(P, s, st)
2% BASS_ASSEMBLE Mass and stiffness matrices of one fingering's coupled
3% string-body system, plus the uncoupled string modes and damping.
4%
5% The displacement ansatz and the resulting matrices are documented in
6% bass_basis.m, which calls this per fingering; bass_validate.m calls it
7% directly to check the closed-form string limit and the eigensolve
8% identities. Coordinates are [q_1..q_ns, xi_1..xi_nb]: mass-normalized
9% pinned-pinned stiff-string modes, then the mass-normalized body modes of
10% bass_body.m. M and Kmat are symmetric and M is positive definite (the
11% coupling is weak: the off-diagonal mass terms are of order
12% sqrt(mu Lv / m_body) / k pi, a few percent here).
14f0open = [97.9989 73.4162 55.0000 41.2034] * P.tuning_trim;
15tension = [P.tension_g_n P.tension_d_n P.tension_a_n P.tension_e_n];
16L0 = P.scale_length_m;
18body = bass_body(P);
19nb = body.n;
20omRes = 2 * pi * 1.25 * max(body.f); % residual admittance high-pass corner
22T = tension(s);
23cw = 2 * L0 * f0open(s); % transverse wave speed of the string
24mu = T / cw^2;
25Lv = L0 * 2^(-st / 12);
26EI = P.inharmonicity * T * L0^2 / pi^2;
28% string modes up to max_mode_hz (at least one, at most 160)
29kEst = min(160, max(1, floor(2 * Lv * P.max_mode_hz / cw) + 2));
30kk = (1:kEst)' * pi / Lv;
31om2 = kk.^2 * (T / mu) + kk.^4 * (EI / mu);
32omAll = sqrt(om2);
33ns = max(1, sum(omAll <= 2 * pi * P.max_mode_hz));
34ns = min(ns, kEst);
35omS = omAll(1:ns);
36kv = (1:ns)';
38% uncoupled string damping: internal loss, the finger termination on
39% stopped notes, and the resistive residual bridge admittance (high-passed
40% above the retained body modes so the modal region is not double-counted;
41% a resistive end gives a constant decay RATE, zeta ~ Y T / (Lv om))
42hp = omS.^2 ./ (omS.^2 + omRes^2);
43zS = 0.5 * (P.loss_factor + P.loss_slope * omS) ...
44 + 0.5 * P.residual_y * (T / Lv) * hp ./ omS;
45if st > 0
46 zS = zS + 0.5 * P.finger_loss;
47end
49% ---- assemble ----
50n = ns + nb;
51gk = sqrt(2 * mu * Lv) * ((-1).^(kv + 1)) ./ (kv * pi);
52M = eye(n);
53Kmat = zeros(n, n);
54for k = 1:ns
55 Kmat(k, k) = omS(k)^2;
56end
57for j = 1:nb
58 for k = 1:ns
59 M(k, ns + j) = gk(k) * body.a(j);
60 M(ns + j, k) = M(k, ns + j);
61 end
62 for i = 1:nb
63 M(ns + i, ns + j) = M(ns + i, ns + j) + (mu * Lv / 3) * body.a(i) * body.a(j);
64 Kmat(ns + i, ns + j) = (T / Lv) * body.a(i) * body.a(j);
65 end
66 Kmat(ns + j, ns + j) = Kmat(ns + j, ns + j) + body.omega(j)^2;
67end
69meta.ns = ns;
70meta.nb = nb;
71meta.n = n;
72meta.Lv = Lv;
73meta.mu = mu;
74meta.T = T;
75meta.f0 = f0open(s) * 2^(st / 12);
76meta.body = body;
77end