function [f, W] = cymbal2_modes(P, m)
% CYMBAL2_MODES Eigenfrequencies (Hz) of the linearized shell for one
% azimuthal wavenumber m, with the membrane (curvature) stiffening and the
% support boundary condition included; the washer is taken as its static
% stiffness (its mass and damping are ignored here, as in the paper's
% linear FEM reference). Used by the demo and the validation, not by the
% renderer. W returns the mode shapes on the real radial nodes l = 3..Nr.

G = cymbal2_grid(P);
nrt = G.nrt;
Nr = G.Nr;
hth = G.hth;
lam2 = (2 * cos(m * hth) - 2) / hth^2;
op = cymbal2_perm(P, G, lam2);

E = P.youngs_modulus;
h = P.thickness_m;
rhoh = P.density * h;
Rs = P.radius_m^2 / (2 * P.rise_m);

% Airy solve operator: Fm * phi = Z * (E h / Rs) * Lap * u, with
% phi = phi_r = 0 at both circles and identity rows on the outermost
% virtual lines.
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;
Zl = zeros(nrt, nrt);
for l = 4:Nr - 1
  Zl(l, :) = op.Lap(l, :);
end
Phimap = Fm \ ((E * h / Rs) * Zl);

% Full stiffness on the equation-of-motion rows l = 3..Nr
Kfull = zeros(nrt, nrt);
Kterm = op.Dflex * op.Bih + (1 / Rs) * op.Lap * Phimap;
for l = 3:Nr
  Kfull(l, :) = Kterm(l, :);
end

% Boundary rows, static. With the washer the shell's equation of motion
% holds on l = 3..Nr and the four virtual lines are eliminated; with the
% rigid clamp u and u_r vanish at l = 3, so u(3) itself joins the
% eliminated set and the equations of motion start at l = 4 (the innermost
% virtual line is then never engaged).
e3 = zeros(1, nrt); e3(3) = 1;
if P.support_clamped ~= 0
  Cmat = [e3; G.D1(3, :); op.MomE; op.ShE];
  keep = 4:Nr;
  elim = [2, 3, Nr + 1, Nr + 2];
else
  Cmat = [op.MomW; op.ShW - P.washer_stiffness * e3; op.MomE; op.ShE];
  keep = 3:Nr;
  elim = [1, 2, Nr + 1, Nr + 2];
end
X = -Cmat(:, elim) \ Cmat(:, keep);
Kred = Kfull(keep, keep) + Kfull(keep, elim) * X;

[V, Dm] = eig(Kred / rhoh);
om2 = real(diag(Dm));
[om2, idx] = sort(om2);
om2 = max(om2, 0);
f = sqrt(om2) / (2 * pi);
W = real(V(:, idx));
end