// Radiation: retarded free-field pressure at the listener point.
// Port of radiate() and check_setup() in woodblock.py. The top and bottom
// faces move together, so each strip is a local dipole; the surface velocity
// sum is computed once per strip and reused for both faces with opposite sign.

import type { WoodblockParams } from '../params';
import type { ModalBasisData } from './basis';
import { elementShapes } from './basis';

export interface RadiateResult {
  p: Float64Array;      // pascals at the listener, at fsHi
  warnings: string[];
}

export function radiate(
  cfg: WoodblockParams,
  Qdd: Float64Array[],
  Xdd: Float64Array,
  basis: ModalBasisData,
  fsHi: number,
  M: number,
  onProgress?: (frac: number) => void,
): RadiateResult {
  const rho0 = cfg.air_density;
  const c = cfg.speed_of_sound;
  const N = basis.n;
  const T = Xdd.length;
  const L = cfg.length_m;
  const W = cfg.width_m;
  const halfH = 0.5 * cfg.block_height_m;
  const patch = (W * L) / M;

  const { phiElem, sElem } = elementShapes(basis, M);

  const lx = cfg.listener_x_m;
  const ly = cfg.listener_y_m;
  const lz = cfg.listener_z_m;

  const distTop = new Float64Array(M);
  const distBot = new Float64Array(M);
  let maxDist = 0;
  for (let j = 0; j < M; j++) {
    const x = (sElem[j] - 0.5) * L;
    const dx = x - lx;
    const dy = 0.0 - ly;
    const dzT = halfH - lz;
    const dzB = -halfH - lz;
    distTop[j] = Math.sqrt(dx * dx + dy * dy + dzT * dzT);
    distBot[j] = Math.sqrt(dx * dx + dy * dy + dzB * dzB);
    maxDist = Math.max(maxDist, distTop[j], distBot[j]);
  }
  const sx = cfg.slit_x_m - lx;
  const sy = cfg.slit_y_m - ly;
  const sz = cfg.slit_z_m - lz;
  const slitDist = Math.sqrt(sx * sx + sy * sy + sz * sz);
  maxDist = Math.max(maxDist, slitDist);

  const maxDelay = Math.ceil((maxDist / c) * fsHi) + 2;
  const out = new Float64Array(T + maxDelay);
  const wdd = new Float64Array(T);

  for (let j = 0; j < M; j++) {
    // Surface acceleration of strip j, shared by the +z and -z sources
    wdd.fill(0);
    for (let n = 0; n < N; n++) {
      const phi = phiElem[n][j];
      if (phi === 0) continue;
      const q = Qdd[n];
      for (let t = 0; t < T; t++) wdd[t] += phi * q[t];
    }
    const ampTop = (patch * rho0) / (4.0 * Math.PI * distTop[j]);
    const ampBot = -(patch * rho0) / (4.0 * Math.PI * distBot[j]);
    const kTop = Math.round((distTop[j] / c) * fsHi);
    const kBot = Math.round((distBot[j] / c) * fsHi);
    for (let t = 0; t < T; t++) {
      out[kTop + t] += ampTop * wdd[t];
      out[kBot + t] += ampBot * wdd[t];
    }
    if (onProgress) onProgress((j + 1) / M);
  }

  // Slit monopole
  const ampSlit = (cfg.slit_area_m2 * rho0) / (4.0 * Math.PI * slitDist);
  const kSlit = Math.round((slitDist / c) * fsHi);
  for (let t = 0; t < T; t++) out[kSlit + t] += ampSlit * Xdd[t];

  // Setup warnings (check_setup in the reference)
  const warnings: string[] = [];
  const spacing = L / M;
  const fMax = basis.f[basis.n - 1];
  if (spacing > (0.25 * c) / fMax) {
    const need = Math.ceil(L / ((0.25 * c) / fMax));
    warnings.push(
      `element spacing is coarse at ${fMax.toFixed(0)} Hz; consider ${need} radiating elements`,
    );
  }
  const d = Math.sqrt(lx * lx + ly * ly + lz * lz);
  if (d < 0.5 * L) {
    warnings.push(
      `listener is ${(d * 100).toFixed(0)} cm from the block centre, inside the source region; point-source elements are crude here`,
    );
  }

  return { p: out.subarray(0, T).slice(), warnings };
}