// The losses derived from physical constants: Zener's thermoelastic loss
// factor has its peak where it should; the outline perimeter of a disc and
// a rectangle; the limits of the average radiation efficiency; the fit of
// the two-parameter operator; the per-mode radiation rate against the
// closed form for a baffled piston and against the direct double sum; the
// per-mode thermoelastic rate of a uniform plate against eta_te omega / 2.

import { test } from "node:test";
import assert from "node:assert/strict";
import { makeGrid, withMask, nodeWeight, RHO_AIR, SPEED_OF_SOUND, ROOM_TEMPERATURE, type Material } from "../src/params.ts";
import { uniformPlate } from "../src/model/reference.ts";
import { assemble } from "../src/model/assemble.ts";
import { computeModes, localRunner } from "../src/modal/modes.ts";
import { modalModel, modalRates, sliceModel } from "../src/modal/synth.ts";
import { dampingPlate, radiationRates, radiationRateDirect, thermoelasticRates, computeDamping, bendingEnergyPerNode } from "../src/modal/damping.ts";
import { thermoelasticRelaxation, thermoelasticLossFactor, outlinePerimeter, averageRadiationEfficiency, coincidenceFrequency, fitLine, fitSigma, fittedRate, plateGeometry, estimatedRate, materialRate, LN1000 } from "../src/loss.ts";

const steel: Material = { E: 200e9, nu: 0.3, rho: 7850, eta: 3e-4, thermal: { alpha: 12e-6, cp: 460, k: 50 } };
const bronze: Material = { E: 110e9, nu: 0.34, rho: 8800, eta: 2e-4, thermal: { alpha: 18e-6, cp: 370, k: 40 } };

test("Zener's thermoelastic loss peaks at omega tau = 1 with the value delta / 2", () => {
  const h = 1e-3;
  const { delta, tau } = thermoelasticRelaxation(bronze, h, ROOM_TEMPERATURE);
  const th = bronze.thermal!;
  assert.ok(Math.abs(tau - (h * h * bronze.rho * th.cp) / (Math.PI * Math.PI * th.k)) < 1e-15);
  assert.ok(Math.abs(delta - ((bronze.E * th.alpha ** 2 * ROOM_TEMPERATURE) / (bronze.rho * th.cp)) * ((1 + bronze.nu) / (1 - bronze.nu))) < 1e-15);
  const peak = thermoelasticLossFactor(bronze, h, 1 / tau, ROOM_TEMPERATURE);
  assert.ok(Math.abs(peak - delta / 2) < 1e-15);
  assert.ok(thermoelasticLossFactor(bronze, h, 0.1 / tau, ROOM_TEMPERATURE) < peak);
  assert.ok(thermoelasticLossFactor(bronze, h, 10 / tau, ROOM_TEMPERATURE) < peak);
  // For a millimeter of bronze the peak lies below 100 Hz and the factor is of order 1e-3.
  assert.ok(1 / (2 * Math.PI * tau) < 100 && delta / 2 > 1e-4 && delta / 2 < 1e-2);
  // Twice the thickness quarters the peak frequency; without thermal constants there is no loss.
  assert.ok(Math.abs(thermoelasticRelaxation(bronze, 2 * h, ROOM_TEMPERATURE).tau / tau - 4) < 1e-12);
  assert.equal(thermoelasticLossFactor({ ...bronze, thermal: null }, h, 1000, ROOM_TEMPERATURE), 0);
});

test("the outline perimeter of a rectangle and of a disc", () => {
  const rect = makeGrid(0.3, 0.2, 0.01);
  // The contour cuts each corner diagonally, which costs about 0.3 cells per corner.
  assert.ok(Math.abs(outlinePerimeter(rect) / (2 * (0.3 + 0.2)) - 1) < 0.02, `rectangle perimeter ${outlinePerimeter(rect)}`);
  const R = 0.2;
  const disc = withMask(makeGrid(2 * R, 2 * R, 0.008), (x, y) => x * x + y * y <= R * R);
  // The contour's segments run in eight directions only, which overstates the length of a curve by a few percent (a staircase would overstate it by 27 percent).
  assert.ok(Math.abs(outlinePerimeter(disc) / (2 * Math.PI * R) - 1) < 0.08, `disc perimeter ${outlinePerimeter(disc)}`);
  assert.ok(outlinePerimeter(disc) > 2 * Math.PI * R, "an overestimate, never an underestimate");
  const ring = withMask(makeGrid(2 * R, 2 * R, 0.008), (x, y) => x * x + y * y <= R * R && x * x + y * y >= 0.05 ** 2);
  assert.ok(Math.abs(outlinePerimeter(ring) / (2 * Math.PI * (R + 0.05)) - 1) < 0.1, "a hole adds its own perimeter");
});

test("the average radiation efficiency has the textbook limits", () => {
  const fc = coincidenceFrequency(bronze, 1e-3);
  assert.ok(fc > 15000 && fc < 20000, `coincidence ${fc} Hz for 1 mm bronze`);
  const A = 0.13, U = 1.28, a = 0.4, b = 0.4;
  // Simply supported: about linear in frequency far below coincidence.
  const s1 = averageRadiationEfficiency(100, fc, A, U, a, b, false), s2 = averageRadiationEfficiency(200, fc, A, U, a, b, false);
  assert.ok(s1 > 0 && s1 < 0.1);
  assert.ok(Math.abs(s2 / s1 - 2) < 0.1, "grows about linearly with frequency far below coincidence");
  assert.ok(averageRadiationEfficiency(fc / 2, fc, A, U, a, b, false) > s2);
  assert.ok(Math.abs(averageRadiationEfficiency(4 * fc, fc, A, U, a, b) - 1 / Math.sqrt(0.75)) < 1e-9, "above coincidence");
  assert.ok(averageRadiationEfficiency(fc * 1.0001, fc, A, U, a, b) <= (0.5 - 0.15) * Math.sqrt(a / (SPEED_OF_SOUND / fc)) + 1e-9, "capped at the coincidence value");
  // Free edges: the same times (f / fc)^2 below coincidence, continuous at coincidence.
  const f1 = averageRadiationEfficiency(100, fc, A, U, a, b), f2 = averageRadiationEfficiency(200, fc, A, U, a, b);
  assert.ok(Math.abs(f1 / s1 - (100 / fc) ** 2) < 1e-12 && Math.abs(f2 / f1 - 8) < 0.4, "free edges: about cubic in frequency far below coincidence");
  assert.ok(Math.abs(averageRadiationEfficiency(fc * 0.9999, fc, A, U, a, b) / averageRadiationEfficiency(fc * 0.9999, fc, A, U, a, b, false) - 1) < 1e-3);
});

test("the two-parameter fit reproduces a line and bounds the estimate", () => {
  const kappa = 1.1;
  const x = [1, 2, 3, 5, 8], y = x.map((v) => 0.7 + 2e-4 * v);
  const fit = fitLine(x, y, kappa);
  assert.ok(Math.abs(fit.sigma0 - 0.7) < 1e-9 && Math.abs(fit.sigma1 - 2e-4) < 1e-12);
  assert.ok(fitLine(x, x.map((v) => 3 - v), kappa).sigma1 === 0, "a falling curve gives sigma1 = 0");
  const g = makeGrid(0.4, 0.4, 0.02);
  const pl = uniformPlate(g, bronze, 1e-3);
  const geom = plateGeometry(pl);
  const f2 = fitSigma(pl, geom);
  assert.ok(f2.sigma0 >= 0 && f2.sigma1 >= 0);
  const e100 = estimatedRate(pl, geom, 100), e5k = estimatedRate(pl, geom, 5000);
  assert.ok(e5k.total > e100.total, "the estimate rises with frequency");
  assert.ok(Math.abs(e100.material - materialRate(bronze.eta, 2 * Math.PI * 100)) < 1e-12);
  assert.ok(e100.radiation > 0 && e100.thermoelastic > 0);
  assert.ok(fittedRate(f2, 1000) > 0.3 * e100.total && fittedRate(f2, 1000) < 3 * e5k.total);
});

/** A small free steel plate with its modes and per-mode losses. */
async function plateWithModes(mat: Material, spacing = 0.0125) {
  const g = makeGrid(0.2, 0.15, spacing);
  const pl = uniformPlate(g, mat, 1e-3);
  const sys = assemble(pl);
  const modes = await computeModes(sys, 4000, localRunner(sys), { sliceSize: 20 });
  const model = modalModel(pl, sys, modes);
  return { pl, model, dp: dampingPlate(pl) };
}

test("the radiation rate of a uniform piston at low frequency equals the closed form", async () => {
  const { model, dp } = await plateWithModes(steel);
  // One fake mode of uniform shape 1 / sqrt(M) at 20 Hz, where the plate is small compared with the wavelength.
  const one = sliceModel(model, 0, 1);
  let M = 0, S = 0;
  for (let k = 0; k < one.nGrid; k++) { M += one.masses[k]; S += dp.area[k]; }
  for (let r = 0; r < one.nRows; r++) one.wShapes[r] = 1 / Math.sqrt(M);
  one.freqs = new Float64Array([20]);
  const omega = 2 * Math.PI * 20;
  const expected = (RHO_AIR * omega * omega * S * S) / (2 * Math.PI * SPEED_OF_SOUND * M); // both faces of a baffled piston
  const fftRate = radiationRates(one, dp, 0, 1)[0];
  assert.ok(Math.abs(fftRate / expected - 1) < 2e-3, `fft ${fftRate} vs ${expected}`);
  assert.ok(Math.abs(radiationRateDirect(one, dp, 0) / expected - 1) < 2e-3);
});

test("the FFT convolution equals the direct double sum for every mode", async () => {
  const { model, dp } = await plateWithModes(steel);
  const rates = radiationRates(model, dp, 0, model.count);
  for (let m = 0; m < model.count; m += Math.max(1, Math.floor(model.count / 12))) {
    const direct = radiationRateDirect(model, dp, m);
    assert.ok(Math.abs(rates[m] - direct) < 1e-9 * Math.max(1, direct), `mode ${m}: ${rates[m]} vs ${direct}`);
  }
  // The high modes radiate faster than the low ones in absolute terms, and the rates are positive.
  for (let m = 0; m < model.count; m++) assert.ok(rates[m] >= 0);
  assert.ok(rates[model.count - 1] > rates[0]);
});

test("simply supported shapes reproduce Wallace's low-frequency radiation efficiencies", async () => {
  // sin(m pi x / a) sin(n pi y / b) on a baffled rectangle at k a << 1 (Wallace 1972):
  // odd-odd 32 k^2 ab / (pi^5 m^2 n^2), odd-even 8 k^4 a b^3 / (3 pi^5 m^2 n^2), even-even 2 k^6 a^3 b^3 / (15 pi^5 m^2 n^2).
  const a = 0.3, b = 0.45, h = 1e-3, f = 30, c = SPEED_OF_SOUND, k = (2 * Math.PI * f) / c;
  const g = makeGrid(a, b, 0.01), pl = uniformPlate(g, steel, h), dp = dampingPlate(pl);
  const rhoH = steel.rho * h;
  for (const [m, n] of [[1, 1], [1, 2], [2, 1], [2, 2], [3, 3]]) {
    const rows: number[] = [], vals: number[] = [];
    let norm = 0;
    for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) {
      const kk = j * g.nx + i, w = nodeWeight(g, i, j);
      if (!w) continue;
      const v = Math.sin((m * Math.PI * i * g.hx) / a) * Math.sin((n * Math.PI * j * g.hy) / b);
      rows.push(kk); vals.push(v); norm += rhoH * g.hx * g.hy * w * v * v;
    }
    const one = { nGrid: g.n, nRows: rows.length, count: 1, freqs: new Float64Array([f]), wShapes: Float32Array.from(vals.map((v) => v / Math.sqrt(norm))), wRow: new Int32Array(0), rowNode: Int32Array.from(rows), sigmaMat: new Float64Array(1), sigmaTe: new Float64Array(1), sigmaRad: new Float64Array(1), masses: new Float64Array(0) };
    const efficiency = (radiationRateDirect(one, dp, 0) * rhoH) / (RHO_AIR * c); // sigma = rho_air c sigma_bar / (rho h)
    const odd = (q: number) => q % 2 === 1;
    const wallace = odd(m) && odd(n) ? (32 * k * k * a * b) / (Math.PI ** 5 * m * m * n * n)
      : !odd(m) && !odd(n) ? (2 * k ** 6 * a ** 3 * b ** 3) / (15 * Math.PI ** 5 * m * m * n * n)
      : odd(m) ? (8 * k ** 4 * a * b ** 3) / (3 * Math.PI ** 5 * m * m * n * n) : (8 * k ** 4 * a ** 3 * b) / (3 * Math.PI ** 5 * m * m * n * n);
    // The 1 cm grid samples the mean of the (3,3) shape about 3 percent low; the 8 mm grid is within 2 percent.
    assert.ok(Math.abs(efficiency / wallace - 1) < 0.04, `(${m},${n}): ${efficiency} vs ${wallace}`);
  }
});

test("the thermoelastic rate of a flat uniform plate is eta_te omega / 2, and the total rate is the sum", async () => {
  const { pl, model, dp } = await plateWithModes(steel);
  const rates = thermoelasticRates(model, dp, 0, model.count);
  for (let m = 0; m < model.count; m++) {
    const omega = 2 * Math.PI * model.freqs[m];
    const expected = (thermoelasticLossFactor(steel, 1e-3, omega, ROOM_TEMPERATURE) * omega) / 2;
    assert.ok(Math.abs(rates[m] / expected - 1) < 1e-4, `mode ${m}: ${rates[m]} vs ${expected}`);
  }
  // The bending energy of a mass-normalized mode of a flat free plate is omega^2 / 2.
  const w = new Float64Array(dp.nx * dp.ny);
  for (let r = 0; r < model.nRows; r++) w[model.rowNode[r]] = model.wShapes[r * model.count + 3];
  const U = bendingEnergyPerNode(dp, w);
  let total = 0;
  for (const u of U) total += u;
  const omega3 = 2 * Math.PI * model.freqs[3];
  assert.ok(Math.abs(total / (omega3 * omega3 / 2) - 1) < 1e-4, `bending energy ${total} vs ${omega3 * omega3 / 2}`);
  computeDamping(model, dp);
  const sigma = modalRates(model);
  for (let m = 0; m < model.count; m++) {
    const omega = 2 * Math.PI * model.freqs[m];
    assert.ok(Math.abs(sigma[m] - (materialRate(steel.eta, omega) + model.sigmaTe[m] + model.sigmaRad[m])) < 1e-12);
    assert.ok(model.sigmaMat[m] > 0 && model.sigmaRad[m] > 0 && model.sigmaTe[m] > 0);
    assert.ok(LN1000 / sigma[m] > 0.1 && LN1000 / sigma[m] < 1000, `T60 ${LN1000 / sigma[m]} s`);
  }
  void pl;
});