// Rayleigh-integral listener with the rim-diffraction back term, as in cymsim/radiation.py, but
// with the retarded-time sum evaluated at the audio rate: the normal acceleration of every node
// is decimated to the audio rate inside the time loop (the same Kaiser-windowed sinc kernel as
// the Python package, applied through DECIMATION_SPAN accumulator slots per node), the element
// means of the finished samples are kept in a short history, and each element's history is read
// at its fractional delay with a windowed-sinc interpolator of DELAY_TAPS taps. The per-element
// coefficients, the delay taps and the decimation taps are computed here. See THEORY.md.
import type { Listener } from "../config/schema.ts";
import type { CymbalProfile } from "./geometry.ts";
import type { Mesh } from "./mesh.ts";

/** Audio samples spanned by the decimation kernel (accumulator slots per node); as in cymsim/decimation.py. */
export const DECIMATION_SPAN = 32;
export const KAISER_BETA = 8.0;
export const CUTOFF = 0.46;
/** Taps of the fractional-delay interpolator, centred: lags i0 - DELAY_TAPS/2 + 1 ... i0 + DELAY_TAPS/2. */
export const DELAY_TAPS = 16;

/** Modified Bessel function I0 by its power series (converges quickly for x < 30). */
export function besselI0(x: number): number {
  let sum = 1, term = 1;
  const y = (x * x) / 4;
  for (let k = 1; k < 60; k++) { term *= y / (k * k); sum += term; if (term < 1e-17 * sum) break; }
  return sum;
}

/** h(d) for d = -W m / 2 ... W m / 2 - 1, unit DC gain; identical to cymsim.decimation.decimation_kernel. */
export function decimationKernel(m: number, W = DECIMATION_SPAN, beta = KAISER_BETA, cutoff = CUTOFF): Float64Array {
  const L = W * m, half = L / 2, fc = cutoff / m, h = new Float64Array(L), i0 = besselI0(beta);
  let sum = 0;
  for (let i = 0; i < L; i++) {
    const d = i - half, x = 2 * Math.PI * fc * d, sinc = d === 0 ? 1 : Math.sin(x) / x;
    const u = d / half, win = besselI0(beta * Math.sqrt(Math.max(1 - u * u, 0))) / i0;
    h[i] = sinc * win; sum += h[i];
  }
  for (let i = 0; i < L; i++) h[i] /= sum;
  return h;
}

/**
 * Decimation taps in polyphase layout: taps[w * m + kk] is the weight with which the sub-step of
 * phase kk (0..m-1) within audio sample j0 contributes to audio sample j0 + w - W/2 + 1, for slot
 * w = 0..W-1. This is the (W, m) layout stored in a .cymsim file.
 */
export function decimationTaps(m: number, W = DECIMATION_SPAN): Float64Array {
  const h = decimationKernel(m, W), taps = new Float64Array(W * m);
  for (let w = 0; w < W; w++) for (let kk = 0; kk < m; kk++) taps[w * m + kk] = h[kk + (W / 2 - 1 - w) * m + (W * m) / 2];
  return taps;
}

/** Windowed-sinc taps that read x(j - i0 - fr) from x(j - i0 - t + T/2 - 1), t = 0..T-1 (Kaiser, beta 6, unit sum). */
export function fractionalDelayTaps(fr: number, T = DELAY_TAPS): Float64Array {
  const taps = new Float64Array(T), half = T / 2, i0 = besselI0(6);
  let sum = 0;
  for (let t = 0; t < T; t++) {
    const k = t - half + 1, x = k - fr;              // lag offset relative to i0 (k = -T/2 + 1 ... T/2)
    const sinc = x === 0 ? 1 : Math.sin(Math.PI * x) / (Math.PI * x);
    const u = x / half, win = besselI0(6 * Math.sqrt(Math.max(1 - u * u, 0))) / i0;
    taps[t] = sinc * win; sum += taps[t];
  }
  for (let t = 0; t < T; t++) taps[t] /= sum;
  return taps;
}

export interface Radiation {
  coefF: Float64Array; coefB: Float64Array;     // rho A / (2 pi r), front and back
  i0F: Uint32Array; i0B: Uint32Array;           // integer part of the delay in audio samples
  frF: Float64Array; frB: Float64Array;         // fractional part
  tapsF: Float64Array; tapsB: Float64Array;     // (E, DELAY_TAPS) interpolator weights
  ringLength: number;                           // element history ring length in audio samples
  taps: Float64Array;                           // (DECIMATION_SPAN, substeps) decimation weights
  rFront: Float64Array; rBack: Float64Array;
}

export function buildRadiation(mesh: Mesh, profile: CymbalProfile, z0Node: Float64Array, listener: Listener, fs: number, substeps: number, nRim = 720): Radiation {
  const E = mesh.e, [Lx, Ly, Lz] = listener.position;
  const cx = new Float64Array(E), cy = new Float64Array(E), cz = new Float64Array(E);
  for (let e = 0; e < E; e++) {
    cx[e] = mesh.centroid[2 * e]; cy[e] = mesh.centroid[2 * e + 1];
    cz[e] = (z0Node[mesh.tri[3 * e]] + z0Node[mesh.tri[3 * e + 1]] + z0Node[mesh.tri[3 * e + 2]]) / 3;
  }
  let R = 0;
  for (let n = 0; n < mesh.n; n++) R = Math.max(R, Math.hypot(mesh.xy[2 * n], mesh.xy[2 * n + 1]));
  const rimX = new Float64Array(nRim), rimY = new Float64Array(nRim), rimZ = new Float64Array(nRim), d2 = new Float64Array(nRim);
  for (let j = 0; j < nRim; j++) {
    const phi = (2 * Math.PI * j) / nRim;
    rimX[j] = R * Math.cos(phi); rimY[j] = R * Math.sin(phi); rimZ[j] = profile.z0(rimX[j], rimY[j]);
    d2[j] = Math.hypot(rimX[j] - Lx, rimY[j] - Ly, rimZ[j] - Lz);
  }
  const rFront = new Float64Array(E), rBack = new Float64Array(E);
  const coefF = new Float64Array(E), coefB = new Float64Array(E), i0F = new Uint32Array(E), frF = new Float64Array(E);
  const i0B = new Uint32Array(E), frB = new Float64Array(E), tapsF = new Float64Array(E * DELAY_TAPS), tapsB = new Float64Array(E * DELAY_TAPS);
  const { c_air: c, rho_air: rho, back_diffraction_gain: g } = listener;
  let maxI0 = 0;
  for (let e = 0; e < E; e++) {
    rFront[e] = Math.hypot(Lx - cx[e], Ly - cy[e], Lz - cz[e]);
    let best = Infinity;
    for (let j = 0; j < nRim; j++) best = Math.min(best, Math.hypot(cx[e] - rimX[j], cy[e] - rimY[j], cz[e] - rimZ[j]) + d2[j]);
    rBack[e] = best;
    coefF[e] = (rho * mesh.area[e]) / (2 * Math.PI * rFront[e]);
    coefB[e] = (g * rho * mesh.area[e]) / (2 * Math.PI * rBack[e]);
    const dF = (rFront[e] / c) * fs, dB = (rBack[e] / c) * fs;
    i0F[e] = Math.floor(dF); frF[e] = dF - i0F[e]; i0B[e] = Math.floor(dB); frB[e] = dB - i0B[e];
    tapsF.set(fractionalDelayTaps(frF[e]), e * DELAY_TAPS); tapsB.set(fractionalDelayTaps(frB[e]), e * DELAY_TAPS);
    maxI0 = Math.max(maxI0, i0F[e], i0B[e]);
  }
  return { coefF, coefB, i0F, i0B, frF, frB, tapsF, tapsB, ringLength: maxI0 + DELAY_TAPS + 4, taps: decimationTaps(substeps), rFront, rBack };
}

/** The delayed, interpolated value of an element history ring at sample j: sum_t taps[t] x(j - i0 - t + T/2 - 1). */
export function delayedValue(hist: ArrayLike<number>, stride: number, e: number, D: number, j: number, i0: number, taps: ArrayLike<number>, tapBase: number): number {
  let s = 0;
  for (let t = 0; t < DELAY_TAPS; t++) {
    const lag = i0 + t - DELAY_TAPS / 2 + 1, jj = j - lag;
    if (jj >= 0 && lag >= 0) s += taps[tapBase + t] * hist[(jj % D) * stride + e];
  }
  return s;
}