// Losses given as decay times. The design script gives T60, the time in
// which a partial decays by 60 dB, as a number, a table [f T60] or a
// function of frequency; the design layer samples it on LOSS_FREQS, and this
// module interpolates it, converts it to decay rates for the modes, and fits
// the two-parameter model of the full simulation to it.

import { type Loss, type Material, bendingConstant } from "./params.ts";

/** Frequencies at which T60 is sampled: 128 points from 10 Hz to 40 kHz, log-spaced. */
export const LOSS_FREQS: Float64Array = (() => {
  const n = 128, f = new Float64Array(n);
  for (let i = 0; i < n; i++) f[i] = 10 * Math.pow(4000, i / (n - 1));
  return f;
})();

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

/** The default losses: 4 s at low frequencies, falling to about 1 s at 5 kHz. */
export function defaultT60(f: number): number {
  return 4 / (1 + f / 2000);
}

export function lossFromFunction(fn: (f: number) => number): Loss {
  const t60 = new Float64Array(LOSS_FREQS.length);
  for (let i = 0; i < t60.length; i++) t60[i] = fn(LOSS_FREQS[i]);
  return { freqs: LOSS_FREQS, t60 };
}

export function defaultLoss(): Loss {
  return lossFromFunction(defaultT60);
}

export function constantLoss(t60: number): Loss {
  return lossFromFunction(() => t60);
}

/** Log-log interpolation of a table of (frequency, T60) rows, clamped at the ends. */
export function lossFromTable(freqs: ArrayLike<number>, t60s: ArrayLike<number>): Loss {
  const n = freqs.length;
  const idx = Array.from({ length: n }, (_, i) => i).sort((a, b) => freqs[a] - freqs[b]);
  const fs = idx.map((i) => freqs[i]), ts = idx.map((i) => t60s[i]);
  return lossFromFunction((f) => interpLogLog(fs, ts, f));
}

function interpLogLog(fs: number[], ts: number[], f: number): number {
  const n = fs.length;
  if (n === 1 || f <= fs[0]) return ts[0];
  if (f >= fs[n - 1]) return ts[n - 1];
  let i = 0;
  while (i < n - 2 && fs[i + 1] < f) i++;
  const t = (Math.log(f) - Math.log(fs[i])) / (Math.log(fs[i + 1]) - Math.log(fs[i]));
  return Math.exp(Math.log(ts[i]) * (1 - t) + Math.log(ts[i + 1]) * t);
}

/** T60 at a frequency, by log-log interpolation of the sampled curve, clamped at its ends. */
export function t60At(loss: Loss, f: number): number {
  return interpLogLog(Array.from(loss.freqs), Array.from(loss.t60), f);
}

/** Decay rate sigma [1/s] of a partial at frequency f: its amplitude decays as exp(-sigma t). */
export function decayRate(loss: Loss, f: number): number {
  return LN1000 / t60At(loss, f);
}

/**
 * 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 design's decay 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(loss: Loss, mat: Material, hMean: number, fLo = 50, fHi = 10000, n = 40): SigmaFit {
  const kappa = bendingConstant(mat, hMean);
  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(decayRate(loss, f));
  }
  // Unconstrained least squares, then clamp: with two parameters the
  // non-negative solution is the unconstrained one, or one parameter at
  // zero with the other refit.
  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;
  let s1 = (n * sxy - sx * sy) / det, s0 = (sy - s1 * sx) / n;
  if (s1 < 0) { s1 = 0; s0 = sy / n; }
  if (s0 < 0) { s0 = 0; s1 = sxy / sxx; }
  return { sigma0: s0, sigma1: s1, 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;
}