function basis = triangle_basis(P)
% TRIANGLE_BASIS Mode set of the open triangle: FEM assembly and eigensolve.
%
% The planar geometry splits the frame into two exactly decoupled families
% (see triangle_fem.m): in-plane modes (bending + stretching) and
% out-of-plane modes (bending + torsion, converted into each other at the
% two closed corners). Each family gives a generalized eigenproblem
% K phi = om^2 M phi, reduced to symmetric standard form through the
% Cholesky factor of the consistent mass matrix and diagonalized by svd
% rather than eig: the reduced matrix is symmetric positive semidefinite,
% so its singular vectors are its eigenvectors, and svd keeps them
% orthonormal even for the near-degenerate pairs this geometry produces
% (eig's nonsymmetric path returns complex, non-orthogonal vectors there).
% The eigenvectors then come back mass-normalized: phi' M phi = I. The
% three rigid-body modes of each family sit at zero frequency and are
% dropped along with anything below min_mode_hz; modes above max_mode_hz
% are truncated (the broadband part they carry is reinstated by the direct
% term in triangle_radiate.m).
%
% For radiation and excitation the basis stores, per retained mode, the
% displacement at every node (NUx, NUy, NUz; out-of-plane modes have only
% NUz, in-plane modes only NUx, NUy) and the transverse displacement at
% every element center (EUpx, EUpy, EUpz), with the axial component removed:
% axial motion of a thin rod displaces almost no air. Element centers,
% tangents, and lengths are kept alongside so the radiation stage can treat
% each element as a point dipole.
%
% Damping: zeta = 0.5 (loss_factor + loss_slope om) plus radiation damping
% from the incoherent sum of the element dipole powers,
%   zeta_rad = rho0 om^3 sum_e (Vdip_e^2 |u_perp,e|^2) / (24 pi c^3),
% Vdip_e = 2 pi a^2 l_e (the factor 2 is the added-mass coefficient of a
% transversally moving cylinder). The incoherent sum ignores the
% cancellation between anti-phase sections, which matters only below the
% coincidence frequency (about 1.9 kHz for an 8 mm steel rod); above it the
% flexural waves are supersonic and the estimate is fair.

geo = triangle_geometry(P);
nn = geo.nn;
nel = nn - 1;
E = P.youngs_modulus;
rho = P.density;
nu = P.poisson_ratio;
d = P.rod_diameter_m;
c = P.speed_of_sound;
rho0 = P.air_density;
nodes = geo.nodes;

% ---- element geometry ----
ex = zeros(nel, 1);
ey = zeros(nel, 1);
tx = zeros(nel, 1);
ty = zeros(nel, 1);
el = zeros(nel, 1);
for e = 1:nel
  dx = nodes(1, e + 1) - nodes(1, e);
  dy = nodes(2, e + 1) - nodes(2, e);
  el(e) = sqrt(dx^2 + dy^2);
  tx(e) = dx / el(e);
  ty(e) = dy / el(e);
  ex(e) = (nodes(1, e) + nodes(1, e + 1)) / 2;
  ey(e) = (nodes(2, e) + nodes(2, e + 1)) / 2;
end

% ---- eigensolve per family, collect retained modes ----
fs = [];
fams = [];
NUx = [];
NUy = [];
NUz = [];
EUpx = [];
EUpy = [];
EUpz = [];
omLo = 2 * pi * P.min_mode_hz;
for family = 0:1
  [K, M] = triangle_fem(nodes, family, E, rho, nu, d);
  Rc = chol(M);
  Am = (Rc') \ (K / Rc);
  Am = (Am + Am') / 2;
  [Uv, Sv, ~] = svd(Am);
  om2 = flipud(diag(Sv));            % ascending
  Vv = fliplr(Uv);
  om = sqrt(max(0, om2)) * P.stiffness_trim;
  keep = find(om > omLo & om < 2 * pi * P.max_mode_hz);
  if length(keep) == 0
    continue;
  end
  Phi = Rc \ Vv(:, keep);            % mass-normalized shapes, 3 nn x nk
  nk = length(keep);
  if family == 0
    ux = Phi(1:3:3 * nn, :);
    uy = Phi(2:3:3 * nn, :);
    uz = zeros(nn, nk);
    eux = (ux(1:nn - 1, :) + ux(2:nn, :)) / 2;
    euy = (uy(1:nn - 1, :) + uy(2:nn, :)) / 2;
    ax = eux .* (tx * ones(1, nk)) + euy .* (ty * ones(1, nk));
    upx = eux - ax .* (tx * ones(1, nk));
    upy = euy - ax .* (ty * ones(1, nk));
    upz = zeros(nel, nk);
  else
    ux = zeros(nn, nk);
    uy = zeros(nn, nk);
    uz = Phi(1:3:3 * nn, :);
    upx = zeros(nel, nk);
    upy = zeros(nel, nk);
    upz = (uz(1:nn - 1, :) + uz(2:nn, :)) / 2;
  end
  fs = [fs; om(keep) / (2 * pi)];
  fams = [fams; family * ones(nk, 1)];
  NUx = [NUx, ux];
  NUy = [NUy, uy];
  NUz = [NUz, uz];
  EUpx = [EUpx, upx];
  EUpy = [EUpy, upy];
  EUpz = [EUpz, upz];
end
n = length(fs);
if n == 0
  error('No rod modes between min_mode_hz and max_mode_hz -- check geometry/material.');
end

[fs, order] = sort(fs);
fs = reshape(fs, n, 1);
fams = reshape(fams(order), n, 1);
NUx = NUx(:, order);
NUy = NUy(:, order);
NUz = NUz(:, order);
EUpx = EUpx(:, order);
EUpy = EUpy(:, order);
EUpz = EUpz(:, order);

% ---- damping: material + radiation ----
omega = 2 * pi * fs;
a = d / 2;
Vdip = 2 * pi * a^2 * el;            % nel x 1 dipole volume per element
up2 = EUpx.^2 + EUpy.^2 + EUpz.^2;   % nel x n
S2 = (Vdip.^2).' * up2;              % 1 x n
zetaRad = rho0 * omega.^3 .* reshape(S2, n, 1) / (24 * pi * c^3);
zeta = 0.5 * (P.loss_factor + P.loss_slope * omega) + P.radiation_loss * zetaRad;

basis.n = n;
basis.f = fs;
basis.omega = omega;
basis.zeta = zeta;
basis.zetaRad = zetaRad;
basis.family = fams;
basis.nodeX = nodes(1, :);
basis.nodeY = nodes(2, :);
basis.sfrac = geo.sfrac;
basis.totalLen = geo.totalLen;
basis.NUx = NUx;
basis.NUy = NUy;
basis.NUz = NUz;
basis.elemX = ex;
basis.elemY = ey;
basis.elemTanX = tx;
basis.elemTanY = ty;
basis.elemLen = el;
basis.EUpx = EUpx;
basis.EUpy = EUpy;
basis.EUpz = EUpz;
end