function basis = woodblock_basis(P)
% WOODBLOCK_BASIS Modal basis: free-free Euler-Bernoulli bar.
% Mass-normalized bending modes and their couplings (woodblock.md section 3).
%   basis.f      modal frequencies (Hz), n x 1
%   basis.phiC   mode shape at the strike point (1/sqrt(kg))
%   basis.A      cavity volume coupling (m^2/sqrt(kg))
%   basis.shapes mass-normalized shapes sampled on basis.s, nGrid x n

L = P.length_m; W = P.width_m; H = P.wall_thickness_m;
mass = P.participating_mass_kg;
if mass == 0
  mass = P.density * W * H * L;
end
cLong = sqrt(P.youngs_modulus / P.density);
% f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)),  I = W H^3/12
prefactor = cLong * H / (2 * pi * sqrt(12) * L^2);

nGrid = 4096;
s = linspace(0, 1, nGrid).';

% Roots of cos(bL)cosh(bL) = 1; asymptotic beyond the 4th
blTab = [4.73004074 7.85320462 10.99560784 14.13716549];
f = zeros(P.n_modes, 1);
phiC = zeros(P.n_modes, 1);
A = zeros(P.n_modes, 1);
shapes = zeros(nGrid, P.n_modes);
f1 = prefactor * blTab(1)^2 * P.hollow_softening;
n = 0;
for m = 1:P.n_modes
  if m <= 4
    bl = blTab(m);
  else
    bl = (2 * m + 1) * pi / 2;
  end
  fm = prefactor * bl^2 * P.hollow_softening;
  if m > 1 && P.inharmonicity ~= 1
    fm = f1 * (fm / f1)^P.inharmonicity;
  end
  if fm > P.max_mode_hz
    break;
  end
  n = n + 1;
  sigma = (cosh(bl) - cos(bl)) / (sinh(bl) - sin(bl));
  x = bl * s;
  shape = (cosh(x) + cos(x)) - sigma * (sinh(x) + sin(x));
  rmsval = sqrt(trapz(s, shape.^2));
  shape = shape / (rmsval * sqrt(mass));
  f(n) = fm;
  phiC(n) = interp1(s, shape, P.strike_position);
  mask = (s >= P.cavity_span_lo) & (s <= P.cavity_span_hi);
  A(n) = W * L * trapz(s, shape .* mask);
  shapes(:, n) = shape;
end
if n == 0
  error('No modes below max_mode_hz -- check geometry/material.');
end
f = f(1:n);
omega = 2 * pi * f;
% assemble by assignment (struct() with array values would make a struct array)
basis.n = n;
basis.f = f;
basis.omega = omega;
basis.zeta = 0.5 * (P.loss_factor + P.loss_slope * omega);
basis.phiC = phiC(1:n);
basis.A = A(1:n) * P.cavity_coupling;
basis.mass = mass;
basis.s = s;
basis.shapes = shapes(:, 1:n);
end