// Per-mode listener gains: the Rayleigh front sum and the diffracted back sum over the elements,
// with each element's travel delay folded into a phase at the mode frequency. The pressure of
// mode m is then Re(G_m * complex acceleration), and the unbaffled estimate is front minus back.
//
// This is the only part of the precompute that depends on the listener, and it is cheap enough
// (one pass over elements times modes) to redo whenever the listener or the frequencies change,
// which is what makes the listener and the material live controls in the design panel. The rim
// points of the diffraction path are the mesh's own boundary nodes rather than a separate polar
// sampling, so that the same function serves the first computation and every later one.
import type { Listener } from "../config/schema.ts";
import type { ModesData } from "../modes.ts";

export interface Gains {
  front: Float32Array; // (Nm, 2) re, im
  back: Float32Array;
}

interface Geometry {
  area: Float64Array; // (E,)
  rFront: Float64Array; // (E,)
  rBack: Float64Array; // (E,)
}

/** Element areas and the two path lengths to the listener; depends on the listener position only
 * through the paths, so it is recomputed whenever the listener moves. */
export function radiationGeometry(data: ModesData, listener: Listener): Geometry {
  const { xy, tri, z0, nElements, nNodes } = data;
  const [Lx, Ly, Lz] = listener.position;
  const area = new Float64Array(nElements);
  const cx = new Float64Array(nElements), cy = new Float64Array(nElements), cz = new Float64Array(nElements);
  for (let e = 0; e < nElements; e++) {
    const a = tri[3 * e], b = tri[3 * e + 1], c = tri[3 * e + 2];
    const x0 = xy[2 * a], y0 = xy[2 * a + 1], x1 = xy[2 * b], y1 = xy[2 * b + 1], x2 = xy[2 * c], y2 = xy[2 * c + 1];
    area[e] = Math.abs(0.5 * ((x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0)));
    cx[e] = (x0 + x1 + x2) / 3;
    cy[e] = (y0 + y1 + y2) / 3;
    cz[e] = (z0[a] + z0[b] + z0[c]) / 3;
  }
  // rim points: the nodes of the outer ring, with their own heights
  let R = 0;
  for (let n = 0; n < nNodes; n++) R = Math.max(R, Math.hypot(xy[2 * n], xy[2 * n + 1]));
  const rimX: number[] = [], rimY: number[] = [], rimZ: number[] = [], rimD: number[] = [];
  for (let n = 0; n < nNodes; n++) {
    const x = xy[2 * n], y = xy[2 * n + 1];
    if (Math.hypot(x, y) < R * (1 - 1e-6)) continue;
    rimX.push(x);
    rimY.push(y);
    rimZ.push(z0[n]);
    rimD.push(Math.hypot(x - Lx, y - Ly, z0[n] - Lz));
  }
  const nRim = rimX.length;
  const rFront = new Float64Array(nElements), rBack = new Float64Array(nElements);
  for (let e = 0; e < nElements; e++) {
    rFront[e] = Math.hypot(Lx - cx[e], Ly - cy[e], Lz - cz[e]);
    let best = Infinity;
    for (let j = 0; j < nRim; j++) {
      const d = Math.hypot(cx[e] - rimX[j], cy[e] - rimY[j], cz[e] - rimZ[j]) + rimD[j];
      if (d < best) best = d;
    }
    rBack[e] = best;
  }
  return { area, rFront, rBack };
}

/** Gains for the given frequencies. ``phiFactor`` rescales the stored mode shapes, which are
 * mass-normalised at the reference density; a density change scales them by 1/sqrt(rho/rho0). */
export function computeGains(data: ModesData, listener: Listener, freq: Float64Array, phiFactor: number,
                             geom?: Geometry, onProgress?: (frac: number) => void): Gains {
  const g = geom ?? radiationGeometry(data, listener);
  const { tri, phiQ, phiScale, nElements, nModes } = data;
  const front = new Float32Array(2 * nModes), back = new Float32Array(2 * nModes);
  const fr = new Float64Array(2 * nModes), bk = new Float64Array(2 * nModes);
  const omega = new Float64Array(nModes);
  for (let m = 0; m < nModes; m++) omega[m] = 2 * Math.PI * freq[m];
  const c = listener.c_air, rhoAir = listener.rho_air, gBack = listener.back_diffraction_gain;
  const phiE = new Float64Array(nModes);
  for (let e = 0; e < nElements; e++) {
    const i0 = tri[3 * e] * nModes, i1 = tri[3 * e + 1] * nModes, i2 = tri[3 * e + 2] * nModes;
    const s = phiFactor / 3;
    for (let m = 0; m < nModes; m++) phiE[m] = s * phiScale[m] * (phiQ[i0 + m] + phiQ[i1 + m] + phiQ[i2 + m]);
    const cf = (rhoAir * g.area[e]) / (2 * Math.PI * g.rFront[e]);
    const cb = (gBack * rhoAir * g.area[e]) / (2 * Math.PI * g.rBack[e]);
    const tf = g.rFront[e] / c, tb = g.rBack[e] / c;
    for (let m = 0; m < nModes; m++) {
      const w = omega[m], p = phiE[m];
      const af = w * tf, ab = w * tb;
      const pf = cf * p, pb = cb * p;
      fr[2 * m] += pf * Math.cos(af);
      fr[2 * m + 1] -= pf * Math.sin(af);
      bk[2 * m] += pb * Math.cos(ab);
      bk[2 * m + 1] -= pb * Math.sin(ab);
    }
    if (onProgress && (e & 255) === 0) onProgress(e / nElements);
  }
  front.set(fr);
  back.set(bk);
  return { front, back };
}