/**
 * The cymbal problem: a reciprocity solve for a listener near a rigid cymbal
 * (problems/cymbal.py in the CUDA proof of concept). A point source at the
 * listener's position L radiates in free space, and the only obstacle is a
 * stationary, sound-hard cymbal:
 *
 *     (Delta + k^2) u = -delta(x - L)   outside the cymbal,
 *     du/dn = 0                          on the cymbal,
 *
 * with the radiation condition. By reciprocity u(x) is also the pressure at the
 * listener due to a unit source at x. We solve the exterior mode of the same
 * impedance equation with lambda = 0, and show the field on three planes through
 * a 1 m box around the cymbal. Lengths in meters, sound speed 343 m/s.
 */
import { CymbalShape, insideTest } from './cymbal.ts';
import type { CVec, Slice } from './problem.ts';

export const C_SOUND = 343;
export const RADIUS = 0.2;
export const LISTENER: [number, number, number] = [0.25, 0, 0.45];
/** the evaluation box: [lo, hi] in x, y, z */
export const BOX: [number, number][] = [
  [-0.5, 0.5],
  [-0.5, 0.5],
  [-0.35, 0.65],
];
/** height of the horizontal slice, below the cymbal */
export const SLICE_BELOW = -0.1;

/**
 * Point sources inside the cymbal whose field is the exact solution outside, for
 * the self-test: positions (3, n) and strengths. A single source would be at most
 * T/2 = 1 cm from the surface, and its data would vary on that scale, which the
 * main problem's data does not. Instead the sources are spread 3 mm apart over
 * the midsurface out to r = rMax, with strengths tapering smoothly to zero.
 */
export function selftestSources(shape: CymbalShape, rMax = 0.12, spacing = 0.003): { x0s: Float64Array; ws: Float64Array } {
  const g: number[] = [];
  for (let t = -rMax; t <= rMax + spacing / 2; t += spacing) g.push(t);
  const pts: number[][] = [];
  const ws: number[] = [];
  for (const x of g) {
    for (const y of g) {
      const rho = x * x + y * y;
      if (rho >= rMax * rMax) continue;
      pts.push([x, y, shape.F(rho).v]);
      ws.push(((1 - rho / rMax ** 2) ** 3 * spacing ** 2) / ((Math.PI * rMax ** 2) / 4));
    }
  }
  const n = pts.length;
  const x0s = new Float64Array(3 * n);
  for (let j = 0; j < n; j++) for (let d = 0; d < 3; d++) x0s[d * n + j] = pts[j][d];
  return { x0s, ws: Float64Array.from(ws) };
}

/** sum_j ws[j] G(x, x0s_j) at points x (3, n), or with nrm given its normal derivative. */
export function sourceSum(x: Float64Array, x0s: Float64Array, ws: Float64Array, k: number, nrm?: Float64Array): CVec {
  const n = x.length / 3;
  const m = ws.length;
  const out = new Float64Array(2 * n);
  for (let i = 0; i < n; i++) {
    let sr = 0;
    let si = 0;
    for (let j = 0; j < m; j++) {
      const dx = x[i] - x0s[j];
      const dy = x[n + i] - x0s[m + j];
      const dz = x[2 * n + i] - x0s[2 * m + j];
      const r = Math.hypot(dx, dy, dz);
      const s = ws[j] / (4 * Math.PI * r);
      let gr = s * Math.cos(k * r);
      let gi = s * Math.sin(k * r);
      if (nrm) {
        const dn = (dx * nrm[i] + dy * nrm[n + i] + dz * nrm[2 * n + i]) / r;
        // g (ik - 1/r) dn
        const tr = (-gr / r - gi * k) * dn;
        const ti = (gr * k - gi / r) * dn;
        gr = tr;
        gi = ti;
      }
      sr += gr;
      si += gi;
    }
    out[2 * i] = sr;
    out[2 * i + 1] = si;
  }
  return out;
}

/** n points around the cymbal, at least 5 mm from its surface (deterministic). */
export function selftestPoints(shape: CymbalShape, n: number): Float64Array {
  const near = insideTest(shape, 0.005);
  const pts: number[][] = [];
  let s = 12345;
  const rnd = () => {
    s = (Math.imul(s, 1664525) + 1013904223) >>> 0;
    return s / 4294967296;
  };
  while (pts.length < n) {
    const x = [-0.3 + 0.6 * rnd(), -0.3 + 0.6 * rnd(), -0.15 + 0.35 * rnd()];
    if (!near(x)) pts.push(x);
  }
  const out = new Float64Array(3 * n);
  for (let j = 0; j < n; j++) for (let d = 0; d < 3; d++) out[d * n + j] = pts[j][d];
  return out;
}

/**
 * Power radiated to infinity, k int |u_inf|^2 dOmega, for u = G(., x0) + S_k sigma +
 * i alpha D_k phi1, with u ~ u_inf exp(ikr)/r. The far field pattern is computed
 * with the smooth rule at the nodes, and the sphere integral with Gauss-Legendre
 * in cos(theta) and the trapezoid rule in phi, exact for the band limit. For a
 * rigid cymbal it equals the power emitted by the source, k/(4 pi) + Im v(x0).
 */
export function farFieldPower(
  k: number,
  x0: ArrayLike<number>,
  xs: Float64Array,
  nrm: Float64Array,
  wts: Float64Array,
  sigma: CVec,
  phi1: CVec,
  alpha: number,
  leggauss: (n: number) => { x: Float64Array; w: Float64Array },
): number {
  const n = wts.length;
  let rmax = Math.hypot(x0[0], x0[1], x0[2]);
  for (let i = 0; i < n; i++) rmax = Math.max(rmax, Math.hypot(xs[i], xs[n + i], xs[2 * n + i]));
  const L = Math.ceil(k * rmax) + 20;
  const { x: ct, w: wct } = leggauss(L + 1);
  const nphi = 2 * L + 2;
  // qs = w sigma, qd = i alpha w phi1
  const qs = new Float64Array(2 * n);
  const qd = new Float64Array(2 * n);
  for (let i = 0; i < n; i++) {
    qs[2 * i] = wts[i] * sigma[2 * i];
    qs[2 * i + 1] = wts[i] * sigma[2 * i + 1];
    qd[2 * i] = -alpha * wts[i] * phi1[2 * i + 1];
    qd[2 * i + 1] = alpha * wts[i] * phi1[2 * i];
  }
  let total = 0;
  for (let a = 0; a <= L; a++) {
    const st = Math.sqrt(1 - ct[a] ** 2);
    for (let b = 0; b < nphi; b++) {
      const ph = (2 * Math.PI * b) / nphi;
      const d = [st * Math.cos(ph), st * Math.sin(ph), ct[a]];
      const p0 = -k * (d[0] * x0[0] + d[1] * x0[1] + d[2] * x0[2]);
      let ur = Math.cos(p0);
      let ui = Math.sin(p0);
      for (let i = 0; i < n; i++) {
        const ph1 = -k * (d[0] * xs[i] + d[1] * xs[n + i] + d[2] * xs[2 * n + i]);
        const er = Math.cos(ph1);
        const ei = Math.sin(ph1);
        const dn = d[0] * nrm[i] + d[1] * nrm[n + i] + d[2] * nrm[2 * n + i];
        // e (qs + (-i k dn) qd)
        const cr = qs[2 * i] + k * dn * qd[2 * i + 1];
        const ci = qs[2 * i + 1] - k * dn * qd[2 * i];
        ur += er * cr - ei * ci;
        ui += er * ci + ei * cr;
      }
      total += ((wct[a] * 2 * Math.PI) / nphi) * (ur * ur + ui * ui);
    }
  }
  return (k * total) / (16 * Math.PI * Math.PI);
}

/**
 * Three planes through the box: xz through the listener, yz through the cymbal's
 * axis, and xy below the cymbal. The mask is 1 outside the cymbal and away from
 * the source.
 */
export function boxSlices(ngrid: number, x0: ArrayLike<number>, shape: CymbalShape): Slice[] {
  const inside = insideTest(shape, 1e-3);
  const defs: [Slice['name'], [number, number], number, number][] = [
    ['xz', [0, 2], 1, x0[1]],
    ['yz', [1, 2], 0, 0],
    ['xy', [0, 1], 2, SLICE_BELOW],
  ];
  const out: Slice[] = [];
  for (const [name, axes, fixed, value] of defs) {
    const N = ngrid * ngrid;
    const pts = new Float64Array(3 * N);
    const mask = new Uint8Array(N);
    const [ea, eb] = [BOX[axes[0]], BOX[axes[1]]];
    const h = (ea[1] - ea[0]) / (ngrid - 1);
    for (let row = 0; row < ngrid; row++) {
      for (let col = 0; col < ngrid; col++) {
        const q = row * ngrid + col;
        const x = [0, 0, 0];
        x[axes[0]] = ea[0] + ((ea[1] - ea[0]) * col) / (ngrid - 1);
        x[axes[1]] = eb[0] + ((eb[1] - eb[0]) * row) / (ngrid - 1);
        x[fixed] = value;
        pts[q] = x[0];
        pts[N + q] = x[1];
        pts[2 * N + q] = x[2];
        const atSource = Math.hypot(x[0] - x0[0], x[1] - x0[1], x[2] - x0[2]) < 1e-3 * h;
        mask[q] = inside(x) || atSource ? 0 : 1;
      }
    }
    out.push({ name, axes, fixed, value, extent: [[ea[0], ea[1]], [eb[0], eb[1]]], pts, mask });
  }
  return out;
}