function S = cymbal2_system(P)
% CYMBAL2_SYSTEM Precompute everything the time loop needs. The linear
% updates block-diagonalize over azimuthal wavenumber, so the implicit
% shell solve and the Airy (stress function) solve are stored as one small
% radial matrix per transform column: the time step becomes
%   Uhat+ = (A^-1 B) Uhat + (A^-1 C) Uhat- + A^-1 (ht^2/rho h) NLhat
% with A, B, C the matrices of Samejima (2021) Eq. (20) restricted to one
% wavenumber, the boundary-condition rows of Eqs. (22)-(30) built in, and
% NL the explicit nonlinear + forcing terms evaluated in physical space.
% Rows are equilibrated before inversion; the stored inverses absorb the
% scaling so the loop applies them to raw right-hand sides.

G = cymbal2_grid(P);
nrt = G.nrt;
Nr = G.Nr;
nth = G.nth;

h = P.thickness_m;
rhoh = P.density * h;
Eh = P.youngs_modulus * h;
Rs = P.radius_m^2 / (2 * P.rise_m);
ht = 1 / (P.sample_rate * P.oversample);
g = P.g_param;
phiw = P.phi_param;
kap = 0;   % set below from cymbal2_perm (needs Dflex)
s1 = P.sigma1 * P.radius_m^2;
s0 = P.sigma0;

MW = P.washer_mass;
RW = P.washer_damping;
KW = P.washer_stiffness;
% Spring term split as phiw/2 on the (k+1, k-1) average and (1 - phiw) on
% step k. Note the paper prints the complementary split, which is unstable
% for a massless washer (the roots of (1-phi)/2 (z^2+1) + phi z include
% |z| > 1 whenever phi > 1/2); with the implicit weight on the average,
% phi = 0.7 is unconditionally stable.
AW = MW + (ht / 2) * RW + (phiw / 2) * ht^2 * KW;
BW = 2 * MW - (1 - phiw) * ht^2 * KW;
CW = -MW + (ht / 2) * RW - (phiw / 2) * ht^2 * KW;

AiB = zeros(nrt, nrt, nth);
AiC = zeros(nrt, nrt, nth);
Ai = zeros(nrt, nrt, nth);
Fi = zeros(nrt, nrt, nth);
nDropped = 0;

e1 = zeros(1, nrt); e1(1) = 1;
e3 = zeros(1, nrt); e3(3) = 1;
I = eye(nrt);

for m = 0:nth / 2
  lam2 = (2 * cos(m * G.hth) - 2) / G.hth^2;
  op = cymbal2_perm(P, G, lam2);
  kap = op.Dflex / rhoh;

  Am = zeros(nrt, nrt);
  Bm = zeros(nrt, nrt);
  Cm = zeros(nrt, nrt);
  if P.support_clamped ~= 0
    l0 = 4;
  else
    l0 = 3;
  end
  for l = l0:Nr
    Am(l, :) = (1 + s0 * ht) * I(l, :) + g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);
    Bm(l, :) = 2 * I(l, :) - (1 - 2 * g) * kap * ht^2 * op.Bih(l, :);
    Cm(l, :) = (s0 * ht - 1) * I(l, :) - g * kap * ht^2 * op.Bih(l, :) - s1 * ht * op.Lap(l, :);
  end
  if P.support_clamped ~= 0
    Am(1, :) = e1;
    Am(2, :) = e3;
    Am(3, :) = G.D1(3, :);
  else
    % The shear side of the washer balance gets the same implicit split as
    % the spring; the paper evaluates it at step k only, which feeds a
    % spurious mode on the virtual lines when the washer is massless.
    Am(1, :) = op.MomW;
    Am(2, :) = AW * e3 - (phiw / 2) * ht^2 * op.ShW;
    Bm(2, :) = BW * e3 + (1 - phiw) * ht^2 * op.ShW;
    Cm(2, :) = CW * e3 + (phiw / 2) * ht^2 * op.ShW;
  end
  Am(Nr + 1, :) = op.MomE;
  Am(Nr + 2, :) = op.ShE;

  % Row equilibration, absorbed into the stored inverses
  sc = 1.0 ./ max(abs(Am), [], 2);
  Asc = sc .* Am;
  Aisc = inv(Asc);
  AiBm = Aisc * (sc .* Bm);
  AiCm = Aisc * (sc .* Cm);
  Aim = Aisc .* sc';   % inv(Am) applied to raw rhs

  % On grids too coarse to resolve a wavenumber's inner boundary layer,
  % the discretized boundary conditions can hand that wavenumber a
  % negative stiffness; such columns are truncated rather than integrated.
  % At the default grid none are.
  sr = max(abs(eig([AiBm, AiCm; I, zeros(nrt)])));
  if sr > 1 + 1e-7
    AiBm = zeros(nrt, nrt);
    AiCm = zeros(nrt, nrt);
    Aim = zeros(nrt, nrt);
    nDropped = nDropped + 1 + (m > 0 && m < nth / 2);
  end

  % Airy solve: biharmonic rows inside, phi = phi_r = 0 at both circles
  Fm = zeros(nrt, nrt);
  Fm(1, 1) = 1;
  Fm(2, :) = G.D1(3, :);
  Fm(3, 3) = 1;
  for l = 4:Nr - 1
    Fm(l, :) = op.Bih(l, :);
  end
  Fm(Nr, Nr) = 1;
  Fm(Nr + 1, :) = G.D1(Nr, :);
  Fm(Nr + 2, Nr + 2) = 1;
  scf = 1.0 ./ max(abs(Fm), [], 2);
  Fim = inv(scf .* Fm) .* scf';

  if m == 0
    cols = 1;
  elseif m == nth / 2
    cols = nth;
  else
    cols = [2 * m, 2 * m + 1];
  end
  for cc = 1:length(cols)
    AiB(:, :, cols(cc)) = AiBm;
    AiC(:, :, cols(cc)) = AiCm;
    Ai(:, :, cols(cc)) = Aim;
    Fi(:, :, cols(cc)) = Fim;
  end
end

if nDropped > nth / 2
  error(['cymbal2: the radial grid cannot resolve the inner boundary layer ' ...
    '(need roughly washer_radius > 3 (a - c) / (n_radial - 3)); raise n_radial']);
end

% Strike and pickup nodes (radial fraction 0 at the washer, 1 at the rim)
ls = 3 + round(min(1, max(0, P.strike_r)) * (Nr - 3));
lp = 3 + round(min(1, max(0, P.pickup_r)) * (Nr - 3));
js = 1 + mod(round(P.strike_theta_deg / 360 * nth), nth);
jp = 1 + mod(round(P.pickup_theta_deg / 360 * nth), nth);

S.G = G;
S.AiB = AiB; S.AiC = AiC; S.Ai = Ai; S.Fi = Fi;
S.rhoh = rhoh; S.Eh = Eh; S.Rs = Rs; S.ht = ht;
S.ls = ls; S.lp = lp; S.js = js; S.jp = jp;
S.eomRow0 = 3 + (P.support_clamped ~= 0);   % first equation-of-motion row
S.nDropped = nDropped;                      % truncated theta columns
S.strikeRadiusM = G.r(ls);
S.pickupRadiusM = G.r(lp);
S.forceScale = 1 / (G.r(ls) * G.hr * G.hth);   % point force -> pressure
end