function [K, M] = cymbal_shell(G, n, E, nu, rho, Lref)
% CYMBAL_SHELL Free-free thin-shell system for one circumferential order n.
%
% The cymbal is a shell of revolution, so its modes separate exactly in the
% azimuth: every mode has n nodal diameters and is
%   u(s, th) = U(s) cos(n th)    (meridional, along the surface toward the rim)
%   v(s, th) = V(s) sin(n th)    (circumferential)
%   w(s, th) = W(s) cos(n th)    (along the outward normal)
% with a partner rotated by pi / 2n. A strike on the meridian th = 0 excites
% only the cos member of each pair, which is why the model solves one member
% per (n, k) and says so in the README. The 1-D energy functional below is the
% th integral of the Love-Kirchhoff shell energy; every quadratic term carries
% cos^2 or sin^2, so the integral contributes the same constant factor (pi, or
% 2 pi at n = 0) to K and M and drops out of the eigenvalues. At n = 0 the V
% field decouples and gives the axisymmetric torsional modes, which have no
% normal motion and therefore radiate nothing.
%
% Strains (Kraus, with the curvatures signed so that a sphere of radius R has
% 1 / R_s = 1 / R_theta = +1 / R for the outward normal):
%   eps_s   = U' + W cs
%   eps_th  = (n V + U r') / r + W ct
%   gamma   = -n U / r + V' - V r' / r            (engineering shear)
%   beta_s  = U cs - W'                           (rotation of the normal)
%   beta_th = n W / r + V ct
%   kap_s   = beta_s'  = U' cs + U cs' - W''
%   kap_th  = (n beta_th + beta_s r') / r
%   kap_sth = -n beta_s / r + beta_th' - beta_th r' / r     (the full twist,
%                                                            factor 2 included)
% and the energy densities are the standard C = E t / (1 - nu^2) membrane and
% D = E t^3 / (12 (1 - nu^2)) bending forms, integrated with weight r ds.
%
% Elements have two nodes and four degrees of freedom each: U and V are linear
% (C0), W is cubic Hermite (C1) with dW/ds as its second freedom, so beta_s and
% its derivative are both well defined. Four-point Gauss quadrature per
% element evaluates the tabulated geometry at the quadrature points rather
% than freezing it per element, so a curved meridian is integrated rather than
% approximated by conical frusta.
%
% Rotary inertia rho t^3 / 12 (beta_s^2 + beta_th^2) is kept. It shifts the
% retained frequencies by well under a percent, but it is what gives the
% dW/ds freedoms mass and so makes M positive definite, which is what lets the
% eigensolve reduce through chol(M). The dW/ds freedoms are then scaled by
% Lref (a congruence, so the eigenvalues are untouched) to bring them onto the
% same numerical scale as the displacements.
%
% Returns K and M, each 4 nn x 4 nn, symmetric, M positive definite. The
% returned system is in the SCALED freedoms; multiply an eigenvector by
% diag(1, 1, 1, Lref) per node to recover physical dW/ds.

nn = G.nn;
ndof = 4 * nn;
K = zeros(ndof);
M = zeros(ndof);

% four-point Gauss on [0, 1]
gx = [0.5 - 0.5 * 0.8611363115940526, 0.5 - 0.5 * 0.3399810435848563, ...
      0.5 + 0.5 * 0.3399810435848563, 0.5 + 0.5 * 0.8611363115940526];
gw = [0.3478548451374538, 0.6521451548625461, ...
      0.6521451548625461, 0.3478548451374538] / 2;

for e = 1:nn - 1
  s0 = G.snode(e);
  s1 = G.snode(e + 1);
  L = s1 - s0;
  Ke = zeros(8);
  Me = zeros(8);
  for g = 1:4
    xi = gx(g);
    sg = s0 + xi * L;
    rr = interp1(G.s, G.r, sg);
    tt = interp1(G.s, G.thk, sg);
    dr = interp1(G.s, G.drds, sg);
    cs = interp1(G.s, G.cs, sg);
    ct = interp1(G.s, G.ct, sg);
    csp = interp1(G.s, G.csp, sg);
    ctp = interp1(G.s, G.ctp, sg);

    % shape functions, DOF order [U1 V1 W1 W1' U2 V2 W2 W2']
    NU = [1 - xi, 0, 0, 0, xi, 0, 0, 0];
    NUp = [-1 / L, 0, 0, 0, 1 / L, 0, 0, 0];
    NV = [0, 1 - xi, 0, 0, 0, xi, 0, 0];
    NVp = [0, -1 / L, 0, 0, 0, 1 / L, 0, 0];
    h1 = 1 - 3 * xi^2 + 2 * xi^3;
    h2 = xi - 2 * xi^2 + xi^3;
    h3 = 3 * xi^2 - 2 * xi^3;
    h4 = -xi^2 + xi^3;
    d1 = -6 * xi + 6 * xi^2;
    d2 = 1 - 4 * xi + 3 * xi^2;
    d3 = 6 * xi - 6 * xi^2;
    d4 = -2 * xi + 3 * xi^2;
    e1 = -6 + 12 * xi;
    e2 = -4 + 6 * xi;
    e3 = 6 - 12 * xi;
    e4 = -2 + 6 * xi;
    NW = [0, 0, h1, L * h2, 0, 0, h3, L * h4];
    NWp = [0, 0, d1 / L, d2, 0, 0, d3 / L, d4];
    NWpp = [0, 0, e1 / L^2, e2 / L, 0, 0, e3 / L^2, e4 / L];

    Bes = NUp + cs * NW;
    Bet = (n * NV + dr * NU) / rr + ct * NW;
    Bg = -n * NU / rr + NVp - dr * NV / rr;
    Bbs = cs * NU - NWp;
    Bbt = n * NW / rr + ct * NV;
    Bks = csp * NU + cs * NUp - NWpp;
    Bkt = (n * Bbt + dr * Bbs) / rr;
    Bbtp = n * (NWp / rr - dr * NW / rr^2) + ctp * NV + ct * NVp;
    Bkst = -n * Bbs / rr + Bbtp - dr * Bbt / rr;

    Cm = E * tt / (1 - nu^2);
    Db = E * tt^3 / (12 * (1 - nu^2));
    wq = gw(g) * L * rr;

    Ke = Ke + wq * (Cm * (Bes' * Bes + Bet' * Bet ...
                          + nu * (Bes' * Bet + Bet' * Bes) ...
                          + (1 - nu) / 2 * (Bg' * Bg)) ...
                    + Db * (Bks' * Bks + Bkt' * Bkt ...
                            + nu * (Bks' * Bkt + Bkt' * Bks) ...
                            + (1 - nu) / 2 * (Bkst' * Bkst)));
    Me = Me + wq * (rho * tt * (NU' * NU + NV' * NV + NW' * NW) ...
                    + rho * tt^3 / 12 * (Bbs' * Bbs + Bbt' * Bbt));
  end
  d = [4 * e - 3, 4 * e - 2, 4 * e - 1, 4 * e, 4 * e + 1, 4 * e + 2, 4 * e + 3, 4 * e + 4];
  K(d, d) = K(d, d) + Ke;
  M(d, d) = M(d, d) + Me;
end

% scale the dW/ds freedoms so every freedom is numerically comparable
sc = ones(ndof, 1);
sc(4:4:ndof) = Lref;
Sd = sc * sc';
K = (K .* Sd + (K .* Sd)') / 2;
M = (M .* Sd + (M .* Sd)') / 2;
end