// The losses of the plate, derived from physical constants rather than
// prescribed. Three mechanisms drain the energy of a vibrating plate in this
// model. Internal friction of the metal is a hysteretic loss with the
// material's loss factor eta: a mode of angular frequency omega decays at the
// rate eta omega / 2. Thermoelastic damping (Zener) comes from the heat that
// flows across the thickness between the compressed and the stretched face
// of a bent plate; it follows from the thermal constants of the metal and
// the local thickness, with a peak at the frequency set by the thermal
// diffusion time across the thickness. Radiation into the air is the energy
// the plate loses as sound; the linear model computes it exactly for each
// mode from the mode's shape (modal/damping.ts), and this module gives the
// analytic average of Leppington for a plate of given area and perimeter,
// which the full model's loss operator is fitted to. The felt of a stand is
// a fourth mechanism, handled with the stand (params.ts, synth.ts).

import { type Grid, type Material, type Plate, RHO_AIR, SPEED_OF_SOUND, bendingConstant, insideCell, massProperties, meanThickness, nodeWeight } from "./params.ts";

/** ln(1000): a rate sigma [1/s] and a decay time relate by sigma = LN1000 / T60. */
export const LN1000 = Math.log(1000);

/** A decay time in seconds for display: an infinity sign for no decay, two decimals below 10 s, whole seconds below 1000 s, a power of ten beyond (a loss that is effectively absent). */
export function t60Text(seconds: number): string {
  if (!isFinite(seconds) || seconds <= 0) return "∞";
  if (seconds < 10) return seconds.toFixed(2);
  if (seconds < 1000) return seconds.toFixed(0);
  return seconds.toExponential(0).replace("e+", "e");
}

/** The decay time of a rate for display (see t60Text). */
export function rateText(sigma: number): string {
  return t60Text(sigma > 0 ? LN1000 / sigma : Infinity);
}

// ---------- Internal friction ----------

/** Decay rate of a mode of angular frequency omega from a hysteretic loss factor eta: the energy lost per cycle is 2 pi eta times the stored energy. */
export function materialRate(eta: number, omega: number): number {
  return (eta * omega) / 2;
}

// ---------- Thermoelastic damping ----------

/**
 * Zener's thermoelastic relaxation of a plate of thickness h: the relaxation
 * strength delta and the relaxation time tau. The strength carries the
 * factor (1 + nu) / (1 - nu) of the biaxial stress state of a bent plate
 * (Zener's beam value is E alpha^2 T / (rho cp)); the time is the thermal
 * diffusion time across the thickness, h^2 / (pi^2 chi) with the diffusivity
 * chi = k / (rho cp).
 */
export function thermoelasticRelaxation(mat: Material, h: number, T: number): { delta: number; tau: number } {
  const th = mat.thermal;
  if (!th) return { delta: 0, tau: 1 };
  const delta = ((mat.E * th.alpha * th.alpha * T) / (mat.rho * th.cp)) * ((1 + mat.nu) / (1 - mat.nu));
  const tau = (h * h * mat.rho * th.cp) / (Math.PI * Math.PI * th.k);
  return { delta, tau };
}

/** Thermoelastic loss factor of bending at angular frequency omega in a plate of thickness h: delta omega tau / (1 + (omega tau)^2), largest (delta / 2) at omega tau = 1. */
export function thermoelasticLossFactor(mat: Material, h: number, omega: number, T: number): number {
  const { delta, tau } = thermoelasticRelaxation(mat, h, T);
  const x = omega * tau;
  return (delta * x) / (1 + x * x);
}

// ---------- Geometry of the outline ----------

/**
 * The length of the outline of the plate, from the marching-squares contour
 * of the cell mask at level one half (cells outside the grid count as
 * outside). A staircase count would overstate the perimeter of a disc by
 * 4 / pi; the contour, whose segments run in eight directions, overstates
 * it by a few percent, which is well within the accuracy of the average
 * radiation efficiency it feeds.
 */
export function outlinePerimeter(g: Grid): number {
  const cx = g.nx - 1, cy = g.ny - 1;
  const inside = (ci: number, cj: number) => ci >= 0 && cj >= 0 && ci < cx && cj < cy && insideCell(g, ci, cj);
  const hx = g.hx, hy = g.hy;
  let length = 0;
  // Squares between the centers of cells (ci, cj), (ci+1, cj), (ci+1, cj+1), (ci, cj+1), one ring beyond the grid.
  for (let cj = -1; cj < cy; cj++) {
    for (let ci = -1; ci < cx; ci++) {
      const a = inside(ci, cj) ? 1 : 0, b = inside(ci + 1, cj) ? 1 : 0, c = inside(ci + 1, cj + 1) ? 1 : 0, d = inside(ci, cj + 1) ? 1 : 0;
      const code = a | (b << 1) | (c << 2) | (d << 3);
      if (code === 0 || code === 15) continue;
      // Midpoints of the square's edges: bottom (a-b), right (b-c), top (c-d), left (d-a), as offsets in cells.
      const bottom = [0.5, 0], right = [1, 0.5], top = [0.5, 1], left = [0, 0.5];
      const seg = (p: number[], q: number[]) => { length += Math.hypot((q[0] - p[0]) * hx, (q[1] - p[1]) * hy); };
      switch (code) {
        case 1: case 14: seg(left, bottom); break;
        case 2: case 13: seg(bottom, right); break;
        case 3: case 12: seg(left, right); break;
        case 4: case 11: seg(right, top); break;
        case 6: case 9: seg(bottom, top); break;
        case 7: case 8: seg(left, top); break;
        case 5: seg(left, bottom); seg(right, top); break; // saddle: two corners
        case 10: seg(bottom, right); seg(left, top); break;
      }
    }
  }
  return length;
}

// ---------- Radiation ----------

/** Coincidence frequency of bending waves in a plate of thickness h: where the bending wavelength equals the wavelength of sound, f_c = c^2 / (2 pi kappa). */
export function coincidenceFrequency(mat: Material, h: number): number {
  return (SPEED_OF_SOUND * SPEED_OF_SOUND) / (2 * Math.PI * bendingConstant(mat, h));
}

/**
 * The radiation efficiency of a baffled plate averaged over its modes near
 * the frequency f. The simply supported case is that of Leppington,
 * Broadbent and Heron (1982), as given by Fahy and Gardonio: below
 * coincidence the edges radiate and the efficiency grows about linearly with
 * frequency, at coincidence it peaks at a value set by the size of the plate
 * in coincidence wavelengths, and above it tends to one. A plate with free
 * edges radiates much less below coincidence: its modes are orthogonal to
 * the rigid translation and rotations, so they have no monopole or dipole
 * moment, and at a free edge the bending moment and the shear vanish, so
 * the boundary terms that make a simply supported edge radiate are smaller
 * by the factor (k / k_b)^2 in amplitude, k the acoustic and k_b the bending
 * wavenumber, that is by (f / f_c)^2 in power. The per-mode computation of
 * modal/damping.ts on the presets follows the simply supported curve times
 * this factor within a factor of two from 100 Hz to 8 kHz, and that is what
 * `freeEdges` applies. `a` and `b` are the shorter and the longer side.
 */
export function averageRadiationEfficiency(f: number, fc: number, area: number, perimeter: number, a: number, b: number, freeEdges = true): number {
  const c = SPEED_OF_SOUND;
  const lambdaC = c / fc;
  const cap = Math.max(1, (0.5 - 0.15 * (a / b)) * Math.sqrt(a / lambdaC));
  if (f >= fc) return Math.min(cap, 1 / Math.sqrt(Math.max(1e-12, 1 - fc / f)));
  const mu = Math.sqrt(fc / f);
  const s = ((perimeter * c) / (4 * Math.PI * Math.PI * area * fc)) * (1 / Math.sqrt(mu * mu - 1)) * (Math.log((mu + 1) / (mu - 1)) + (2 * mu) / (mu * mu - 1));
  return Math.min(cap, s) * (freeEdges ? (f / fc) ** 2 : 1);
}

/** Decay rate from a radiation efficiency: the two faces of a baffled plate of surface mass rho h radiate the power 2 rho_air c sigma_rad <v^2> A, so sigma = rho_air c sigma_rad / (rho h). */
export function radiationRateFromEfficiency(efficiency: number, mat: Material, h: number): number {
  return (RHO_AIR * SPEED_OF_SOUND * efficiency) / (mat.rho * h);
}

// ---------- The estimate for a whole plate ----------

/** Geometric quantities of a plate that the loss estimate needs. */
export interface PlateGeometry {
  area: number;
  perimeter: number;
  hMean: number;
  kappa: number;
  fc: number;
  /** Shorter and longer side of the bounding box [m]. */
  a: number;
  b: number;
  /** Node areas over the grid [m^2] (zero outside the plate), for the area-weighted thermoelastic loss. */
  nodeArea: Float64Array;
}

export function plateGeometry(pl: Plate): PlateGeometry {
  const g = pl.grid;
  const { area } = massProperties(pl);
  const hMean = meanThickness(pl);
  const nodeArea = new Float64Array(g.n);
  for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) nodeArea[j * g.nx + i] = g.hx * g.hy * nodeWeight(g, i, j);
  return { area, perimeter: outlinePerimeter(g), hMean, kappa: bendingConstant(pl.mat, hMean), fc: coincidenceFrequency(pl.mat, hMean), a: Math.min(pl.Lx, pl.Ly), b: Math.max(pl.Lx, pl.Ly), nodeArea };
}

/** The three contributions to the decay rate of a bending partial at frequency f [1/s], and their sum. */
export interface LossEstimate {
  material: number;
  thermoelastic: number;
  radiation: number;
  total: number;
}

/**
 * The decay rate of a bending partial at frequency f estimated without the
 * modes: the material rate, the thermoelastic rate averaged over the plate
 * area with the local thickness, and the radiation rate from the average
 * efficiency at the mean thickness. The modal computation replaces the last
 * two by exact per-mode values; this estimate serves the full model's loss
 * operator and the summary of a design before its modes exist.
 */
export function estimatedRate(pl: Plate, geom: PlateGeometry, f: number): LossEstimate {
  const omega = 2 * Math.PI * f;
  const material = materialRate(pl.mat.eta, omega);
  let te = 0;
  if (pl.mat.thermal) {
    let num = 0, den = 0;
    for (let k = 0; k < geom.nodeArea.length; k++) {
      const a = geom.nodeArea[k];
      if (a === 0) continue;
      num += a * thermoelasticLossFactor(pl.mat, pl.h[k], omega, pl.temperature);
      den += a;
    }
    te = den > 0 ? ((num / den) * omega) / 2 : 0;
  }
  const radiation = radiationRateFromEfficiency(averageRadiationEfficiency(f, geom.fc, geom.area, geom.perimeter, geom.a, geom.b), pl.mat, geom.hMean);
  return { material, thermoelastic: te, radiation, total: material + te + radiation };
}

/** The in-plane motion of the full model has no radiation and no thermoelastic loss in this model; it decays by internal friction at the rate of the lowest in-plane resonance, c_p / (2 L). */
export function inPlaneRate(pl: Plate): number {
  const cp = Math.sqrt(pl.mat.E / (pl.mat.rho * (1 - pl.mat.nu * pl.mat.nu)));
  const fIp = cp / (2 * Math.max(pl.Lx, pl.Ly));
  return materialRate(pl.mat.eta, 2 * Math.PI * fIp);
}

// ---------- The two-parameter operator of the full model ----------

/**
 * The two-parameter loss model of the full simulation. A bending mode of
 * angular frequency omega decays at the rate sigma0 + sigma1 k^2 with
 * k^2 = omega / kappa; the fit is a non-negative least-squares fit of that
 * line to the estimated rates on log-spaced frequencies from 50 Hz to
 * 10 kHz, at the mean thickness of the plate.
 */
export interface SigmaFit {
  sigma0: number; // [1/s]
  sigma1: number; // [m^2/s]
  kappa: number; // [m^2/s] bending constant used for the fit
}

export function fitSigma(pl: Plate, geom: PlateGeometry = plateGeometry(pl), fLo = 50, fHi = 10000, n = 40): SigmaFit {
  const kappa = geom.kappa;
  const x: number[] = [], y: number[] = [];
  for (let i = 0; i < n; i++) {
    const f = fLo * Math.pow(fHi / fLo, i / (n - 1));
    x.push((2 * Math.PI * f) / kappa); // k^2
    y.push(estimatedRate(pl, geom, f).total);
  }
  return fitLine(x, y, kappa);
}

/** Non-negative least squares of y = s0 + s1 x: with two parameters the solution is the unconstrained one when it is feasible, else the best of the fits with one parameter at zero, else zero. */
export function fitLine(x: number[], y: number[], kappa: number): SigmaFit {
  const n = x.length;
  let sx = 0, sy = 0, sxx = 0, sxy = 0;
  for (let i = 0; i < n; i++) { sx += x[i]; sy += y[i]; sxx += x[i] * x[i]; sxy += x[i] * y[i]; }
  const det = n * sxx - sx * sx;
  const candidates: [number, number][] = [[0, 0]];
  if (det > 0) { const s1 = (n * sxy - sx * sy) / det; candidates.push([(sy - s1 * sx) / n, s1]); }
  candidates.push([sy / n, 0]);
  if (sxx > 0) candidates.push([0, sxy / sxx]);
  let best: [number, number] = [0, 0], bestRes = Infinity;
  for (const [s0, s1] of candidates) {
    if (!(s0 >= 0 && s1 >= 0)) continue;
    let res = 0;
    for (let i = 0; i < n; i++) res += (y[i] - s0 - s1 * x[i]) ** 2;
    if (res < bestRes) { bestRes = res; best = [s0, s1]; }
  }
  return { sigma0: best[0], sigma1: best[1], kappa };
}

/** Decay rate of the fitted model at a frequency. */
export function fittedRate(fit: SigmaFit, f: number): number {
  return fit.sigma0 + (fit.sigma1 * 2 * Math.PI * f) / fit.kappa;
}