/**
 * Double precision CPU reference for the tests: the six kernels, and brute
 * force near weights for one (target, patch) pair (a port of
 * tests/test_quadrature.py from the CUDA proof of concept).
 *
 * The six kernels, in the order used by fmm3dbie's helm_rpcomb_imp solver:
 *   0: S'_k   1: S_kappa   2: S'_kappa   3: D'_k - D'_kappa   4: S_k   5: D_k
 * with kappa = i|k|, S' the normal derivative at the target and D the normal
 * derivative at the source.
 */
import { leggauss } from './gauss.ts';
import { basis2d, type Surface } from './surface.ts';

const OVER4PI = 1 / (4 * Math.PI);

/** K (6, n) complex, interleaved per kernel: out[m][2q], out[m][2q+1]. */
export function kernels6(xt: ArrayLike<number>, nt: ArrayLike<number>, y: Float64Array, ns: Float64Array, k: number): Float64Array[] {
  const n = y.length / 3;
  const K = Array.from({ length: 6 }, () => new Float64Array(2 * n));
  for (let q = 0; q < n; q++) {
    const dx = xt[0] - y[q];
    const dy = xt[1] - y[n + q];
    const dz = xt[2] - y[2 * n + q];
    const r2 = dx * dx + dy * dy + dz * dz;
    const r = Math.sqrt(r2);
    const dnt = dx * nt[0] + dy * nt[1] + dz * nt[2];
    const dns = dx * ns[q] + dy * ns[n + q] + dz * ns[2 * n + q];
    const nn = ns[q] * nt[0] + ns[n + q] * nt[1] + ns[2 * n + q] * nt[2];
    const e0 = OVER4PI / r;
    const Ekr = e0 * Math.cos(k * r);
    const Eki = e0 * Math.sin(k * r);
    const Ec = e0 * Math.exp(-k * r);
    // fk = 1 - i k r, fc = 1 + k r
    const fkr = 1;
    const fki = -k * r;
    const fc = 1 + k * r;
    // gk = 3 - 3 i k r - k^2 r^2, gc = 3 + 3 k r + k^2 r^2
    const gkr = 3 - k * k * r2;
    const gki = -3 * k * r;
    const gc = 3 + 3 * k * r + k * k * r2;
    const fkEkr = fkr * Ekr - fki * Eki;
    const fkEki = fkr * Eki + fki * Ekr;
    const gkEkr = gkr * Ekr - gki * Eki;
    const gkEki = gkr * Eki + gki * Ekr;
    K[4][2 * q] = Ekr;
    K[4][2 * q + 1] = Eki;
    K[1][2 * q] = Ec;
    K[0][2 * q] = (-dnt * fkEkr) / r2;
    K[0][2 * q + 1] = (-dnt * fkEki) / r2;
    K[2][2 * q] = (-dnt * fc * Ec) / r2;
    K[5][2 * q] = (dns * fkEkr) / r2;
    K[5][2 * q + 1] = (dns * fkEki) / r2;
    const c1 = nn / r2;
    const c2 = (dns * dnt) / (r2 * r2);
    K[3][2 * q] = c1 * fkEkr - c2 * gkEkr - (c1 * fc * Ec - c2 * gc * Ec);
    K[3][2 * q + 1] = c1 * fkEki - c2 * gkEki;
  }
  return K;
}

/** Moments sum_nodes K_m P_n |J| w -> (6, npols) complex. */
export function momentsAt(surf: Surface, patch: number, uv: Float64Array, w: Float64Array, xt: ArrayLike<number>, nt: ArrayLike<number>, k: number): Float64Array[] {
  const n = w.length;
  const g = surf.evalGeometry(patch, uv);
  const y = g.subarray(0, 3 * n);
  const ns = new Float64Array(3 * n);
  const jac = new Float64Array(n);
  for (let q = 0; q < n; q++) {
    const ru = [g[3 * n + q], g[4 * n + q], g[5 * n + q]];
    const rv = [g[6 * n + q], g[7 * n + q], g[8 * n + q]];
    const cx = ru[1] * rv[2] - ru[2] * rv[1];
    const cy = ru[2] * rv[0] - ru[0] * rv[2];
    const cz = ru[0] * rv[1] - ru[1] * rv[0];
    jac[q] = Math.hypot(cx, cy, cz);
    ns[q] = cx / jac[q];
    ns[n + q] = cy / jac[q];
    ns[2 * n + q] = cz / jac[q];
  }
  const K = kernels6(xt, nt, y, ns, k);
  const P = basis2d(uv, surf.norder);
  const npols = surf.npols;
  const I = Array.from({ length: 6 }, () => new Float64Array(2 * npols));
  for (let m = 0; m < 6; m++) {
    for (let nn = 0; nn < npols; nn++) {
      let sr = 0;
      let si = 0;
      for (let q = 0; q < n; q++) {
        const f = P[nn * n + q] * jac[q] * w[q];
        sr += K[m][2 * q] * f;
        si += K[m][2 * q + 1] * f;
      }
      I[m][2 * nn] = sr;
      I[m][2 * nn + 1] = si;
    }
  }
  return I;
}

function addInto(I: Float64Array[], J: Float64Array[]) {
  for (let m = 0; m < 6; m++) for (let i = 0; i < I[m].length; i++) I[m][i] += J[m][i];
}

/** Self-patch moments: polar rules on the triangles (target, edge), split at the foot of the perpendicular, many nodes. */
export function refSelf(surf: Surface, patch: number, u0: number, v0: number, xt: ArrayLike<number>, nt: ArrayLike<number>, k: number, nx = 300, nr = 60): Float64Array[] {
  const gr = leggauss(nr);
  const gx = leggauss(nx);
  const I = Array.from({ length: 6 }, () => new Float64Array(2 * surf.npols));
  for (let e = 0; e < 4; e++) {
    const h = [1 - u0, -1 - u0, 1 - v0, -1 - v0][e];
    if (Math.abs(h) < 1e-14) continue;
    const [x1, x2] = e < 2 ? [-1 - v0, 1 - v0] : [-1 - u0, 1 - u0];
    for (const xend of [x1, x2]) {
      if (Math.abs(xend) < 1e-14) continue;
      const xi1 = Math.asinh(Math.abs(xend) / Math.abs(h));
      const n = nr * nx;
      const uv = new Float64Array(2 * n);
      const w = new Float64Array(n);
      for (let a = 0; a < nr; a++) {
        const rho = 0.5 * (gr.x[a] + 1);
        const wr = 0.5 * gr.w[a];
        for (let b = 0; b < nx; b++) {
          const xi = 0.5 * xi1 * (gx.x[b] + 1);
          const x = Math.sign(xend) * Math.abs(h) * Math.sinh(xi);
          const dx = Math.abs(h) * Math.cosh(xi) * 0.5 * xi1 * gx.w[b];
          const q = a * nx + b;
          if (e < 2) {
            uv[q] = u0 + rho * h;
            uv[n + q] = v0 + rho * x;
          } else {
            uv[q] = u0 + rho * x;
            uv[n + q] = v0 + rho * h;
          }
          w[q] = wr * dx * Math.abs(h) * rho;
        }
      }
      addInto(I, momentsAt(surf, patch, uv, w, xt, nt, k));
    }
  }
  return I;
}

/** Near-patch moments with a uniform quadtree of the given depth and an nq x nq rule per leaf. */
export function refLeaf(surf: Surface, patch: number, xt: ArrayLike<number>, nt: ArrayLike<number>, k: number, depth = 5, nq = 8): Float64Array[] {
  const g = leggauss(nq);
  const nl = 2 ** depth;
  const h = 1 / nl;
  const n = nl * nl * nq * nq;
  const uv = new Float64Array(2 * n);
  const w = new Float64Array(n);
  let q = 0;
  for (let i = 0; i < nl; i++) {
    for (let j = 0; j < nl; j++) {
      const cu = -1 + h + 2 * h * i;
      const cv = -1 + h + 2 * h * j;
      for (let a = 0; a < nq; a++) {
        for (let b = 0; b < nq; b++) {
          uv[q] = cu + h * g.x[a];
          uv[n + q] = cv + h * g.x[b];
          w[q] = h * h * g.w[a] * g.w[b];
          q++;
        }
      }
    }
  }
  return momentsAt(surf, patch, uv, w, xt, nt, k);
}

/** Node weights w_j = sum_n U[n, j] I_n from the moments: (6, npols) complex. */
export function momentsToWeights(surf: Surface, I: Float64Array[]): Float64Array[] {
  const npols = surf.npols;
  const W = Array.from({ length: 6 }, () => new Float64Array(2 * npols));
  for (let m = 0; m < 6; m++) {
    for (let j = 0; j < npols; j++) {
      let sr = 0;
      let si = 0;
      for (let n = 0; n < npols; n++) {
        const u = surf.umat[n * npols + j];
        sr += u * I[m][2 * n];
        si += u * I[m][2 * n + 1];
      }
      W[m][2 * j] = sr;
      W[m][2 * j + 1] = si;
    }
  }
  return W;
}