// Modal basis: free-free Euler-Bernoulli bar, mass-normalized bending modes.
// Port of ModalBasis in woodblock.py.

import type { WoodblockParams } from '../params';

// Roots of cos(bL)cosh(bL) = 1 for the free-free beam.
const BETA_L = [4.73004074, 7.85320462, 10.99560784, 14.13716549];

export function betaL(n: number): number {
  if (n <= BETA_L.length) return BETA_L[n - 1];
  return ((2 * n + 1) * Math.PI) / 2.0;
}

// Free-free mode shape on s in [0,1], unnormalized.
export function modeShape(bl: number, s: Float64Array): Float64Array {
  const sigma = (Math.cosh(bl) - Math.cos(bl)) / (Math.sinh(bl) - Math.sin(bl));
  const out = new Float64Array(s.length);
  for (let i = 0; i < s.length; i++) {
    const x = bl * s[i];
    out[i] = Math.cosh(x) + Math.cos(x) - sigma * (Math.sinh(x) + Math.sin(x));
  }
  return out;
}

function trapz(y: Float64Array, ds: number): number {
  let sum = 0;
  for (let i = 0; i < y.length; i++) sum += y[i];
  sum -= 0.5 * (y[0] + y[y.length - 1]);
  return sum * ds;
}

function interp1(x: number, xs: Float64Array, ys: Float64Array): number {
  if (x <= xs[0]) return ys[0];
  const n = xs.length;
  if (x >= xs[n - 1]) return ys[n - 1];
  // uniform grid
  const dx = xs[1] - xs[0];
  const j = Math.min(n - 2, Math.floor((x - xs[0]) / dx));
  const t = (x - xs[j]) / dx;
  return ys[j] * (1 - t) + ys[j + 1] * t;
}

export interface ModalBasisData {
  n: number;
  f: Float64Array;            // modal frequencies (Hz)
  omega: Float64Array;
  zeta: Float64Array;
  phiC: Float64Array;         // mode shape at strike point (1/sqrt(kg))
  A: Float64Array;            // cavity volume coupling (m^2/sqrt(kg))
  mass: number;
  // Mass-normalized shapes sampled on the dense s grid, for later
  // interpolation at arbitrary radiating-element positions.
  sGrid: Float64Array;
  shapes: Float64Array[];     // length n, each sampled on sGrid
}

const NS = 4096;

export function buildBasis(cfg: WoodblockParams): ModalBasisData {
  const L = cfg.length_m;
  const W = cfg.width_m;
  const H = cfg.wall_thickness_m;
  const area = W * H;
  const mass = cfg.participating_mass_kg || cfg.density * area * L;
  const cLong = Math.sqrt(cfg.youngs_modulus / cfg.density);

  // f_n = (beta_n L)^2 / (2 pi L^2) * sqrt(E I / (rho A)),  I = W H^3/12
  const prefactor = (cLong * H) / (2.0 * Math.PI * Math.sqrt(12.0) * L * L);

  const s = new Float64Array(NS);
  for (let i = 0; i < NS; i++) s[i] = i / (NS - 1);
  const ds = 1.0 / (NS - 1);

  const freqs: number[] = [];
  const phiC: number[] = [];
  const A: number[] = [];
  const shapes: Float64Array[] = [];
  const lo = cfg.cavity_span_lo;
  const hi = cfg.cavity_span_hi;

  for (let n = 1; n <= cfg.n_modes; n++) {
    const bl = betaL(n);
    let f = prefactor * bl * bl * cfg.hollow_softening;
    if (n > 1 && cfg.inharmonicity !== 1.0) {
      const f0 = prefactor * betaL(1) ** 2 * cfg.hollow_softening;
      f = f0 * (f / f0) ** cfg.inharmonicity;
    }
    if (f > cfg.max_mode_hz) break;

    const raw = modeShape(bl, s);
    const sq = new Float64Array(NS);
    for (let i = 0; i < NS; i++) sq[i] = raw[i] * raw[i];
    const rms = Math.sqrt(trapz(sq, ds));
    const norm = 1.0 / (rms * Math.sqrt(mass));
    const shape = new Float64Array(NS);
    for (let i = 0; i < NS; i++) shape[i] = raw[i] * norm;

    freqs.push(f);
    phiC.push(interp1(cfg.strike_position, s, shape));
    const masked = new Float64Array(NS);
    for (let i = 0; i < NS; i++) masked[i] = s[i] >= lo && s[i] <= hi ? shape[i] : 0.0;
    A.push(W * L * trapz(masked, ds));
    shapes.push(shape);
  }

  const nKept = freqs.length;
  const fArr = Float64Array.from(freqs);
  const omega = new Float64Array(nKept);
  const zeta = new Float64Array(nKept);
  const AArr = new Float64Array(nKept);
  for (let i = 0; i < nKept; i++) {
    omega[i] = 2 * Math.PI * fArr[i];
    zeta[i] = 0.5 * (cfg.loss_factor + cfg.loss_slope * omega[i]);
    AArr[i] = A[i] * cfg.cavity_coupling;
  }

  return {
    n: nKept,
    f: fArr,
    omega,
    zeta,
    phiC: Float64Array.from(phiC),
    A: AArr,
    mass,
    sGrid: s,
    shapes,
  };
}

// Interpolate the stored mode shapes at M radiating-element centres.
// Returns phiElem as an array of n Float64Arrays of length M, plus sElem.
export function elementShapes(basis: ModalBasisData, M: number): {
  phiElem: Float64Array[];
  sElem: Float64Array;
} {
  const sElem = new Float64Array(M);
  for (let j = 0; j < M; j++) sElem[j] = (j + 0.5) / M;
  const phiElem: Float64Array[] = [];
  for (let n = 0; n < basis.n; n++) {
    const row = new Float64Array(M);
    for (let j = 0; j < M; j++) row[j] = interp1(sElem[j], basis.sGrid, basis.shapes[n]);
    phiElem.push(row);
  }
  return { phiElem, sElem };
}