/**
 * Node-side tests for the double-precision planner and reference transform.
 * Run: node scripts/test-node.ts
 */
import { gaussNodesWeights } from '../src/gauss.ts';
import { legendreCoeffs, legendreRow } from '../src/coeffs.ts';
import { lmIndex } from '../src/layout.ts';
import { ShtReference, randomSpectrum } from '../src/reference.ts';

let failures = 0;
function check(name: string, ok: boolean, detail = '') {
  console.log(`${ok ? 'PASS' : 'FAIL'}  ${name}${detail ? '  (' + detail + ')' : ''}`);
  if (!ok) failures++;
}

// --- Gauss quadrature ---
{
  const { x, w } = gaussNodesWeights(64);
  let sw = 0, sx2 = 0;
  for (let i = 0; i < 64; i++) {
    sw += w[i];
    sx2 += w[i] * x[i] * x[i];
  }
  check('gauss: sum(w) == 2', Math.abs(sw - 2) < 1e-13, `err=${Math.abs(sw - 2).toExponential(2)}`);
  check('gauss: int x^2 == 2/3', Math.abs(sx2 - 2 / 3) < 1e-13, `err=${Math.abs(sx2 - 2 / 3).toExponential(2)}`);
  check('gauss: nodes decreasing', x[0] > x[1] && x[1] > x[63]);
}

// --- Analytic Legendre values (orthonormal + Condon-Shortley) ---
{
  const lmax = 8;
  const coeffs = legendreCoeffs(lmax, lmax);
  const theta = 0.7;
  const ct = Math.cos(theta), st = Math.sin(theta);
  const row = new Float64Array(lmax + 1);

  legendreRow(coeffs, lmax, 0, ct, st, row);
  const y00 = Math.sqrt(1 / (4 * Math.PI));
  const y10 = Math.sqrt(3 / (4 * Math.PI)) * ct;
  const y20 = Math.sqrt(5 / (16 * Math.PI)) * (3 * ct * ct - 1);
  check('Y_0^0', Math.abs(row[0] - y00) < 1e-14);
  check('Y_1^0', Math.abs(row[1] - y10) < 1e-14);
  check('Y_2^0', Math.abs(row[2] - y20) < 1e-14);

  legendreRow(coeffs, lmax, 1, ct, st, row);
  const y11 = -Math.sqrt(3 / (8 * Math.PI)) * st; // CS phase => negative
  const y21 = -Math.sqrt(15 / (8 * Math.PI)) * st * ct;
  check('Y_1^1 (CS phase)', Math.abs(row[0] - y11) < 1e-14, `got ${row[0]}, want ${y11}`);
  check('Y_2^1', Math.abs(row[1] - y21) < 1e-14);

  legendreRow(coeffs, lmax, 2, ct, st, row);
  const y22 = Math.sqrt(15 / (32 * Math.PI)) * st * st;
  check('Y_2^2', Math.abs(row[0] - y22) < 1e-14);
}

// --- Orthonormality under Gauss quadrature ---
{
  const lmax = 42, nlat = 48;
  const coeffs = legendreCoeffs(lmax, lmax);
  const { x, w } = gaussNodesWeights(nlat);
  let worst = 0;
  for (const m of [0, 1, 7, 25]) {
    const rowsI = new Float64Array(lmax + 1);
    const rowsJ = new Float64Array(lmax + 1);
    for (const [la, lb] of [[m, m], [m, m + 3], [lmax, lmax], [m + 1, lmax]] as const) {
      if (la > lmax || lb > lmax) continue;
      let s = 0;
      for (let i = 0; i < nlat; i++) {
        const st = Math.sqrt(1 - x[i] * x[i]);
        legendreRow(coeffs, lmax, m, x[i], st, rowsI);
        legendreRow(coeffs, lmax, m, x[i], st, rowsJ);
        s += w[i] * rowsI[la - m] * rowsJ[lb - m];
      }
      const want = la === lb ? 1 / (2 * Math.PI) : 0;
      worst = Math.max(worst, Math.abs(s - want));
    }
  }
  check('orthonormality: max err < 1e-12', worst < 1e-12, `worst=${worst.toExponential(2)}`);
}

// --- Reference round trip ---
{
  const cfg = { lmax: 31, mmax: 31, nlat: 34, nphi: 64 };
  const ref = new ShtReference(cfg);
  const q0 = randomSpectrum(cfg, 999);
  const spat = ref.synth(q0);
  const q1 = ref.analys(spat);
  let num = 0, den = 0;
  for (let k = 0; k < q0.length; k++) {
    num += (q1[k] - q0[k]) ** 2;
    den += q0[k] ** 2;
  }
  const rel = Math.sqrt(num / den);
  check('reference round trip rel L2 < 1e-12', rel < 1e-12, `rel=${rel.toExponential(2)}`);
}

// --- Mean value: Y_00 coefficient of a constant field ---
{
  const cfg = { lmax: 15, mmax: 15, nlat: 16, nphi: 32 };
  const ref = new ShtReference(cfg);
  const spat = new Float64Array(cfg.nlat * cfg.nphi).fill(1.0);
  const q = ref.analys(spat);
  const want = Math.sqrt(4 * Math.PI); // <1, Y00> = sqrt(4pi)
  check('constant field -> Q_00 = sqrt(4pi)', Math.abs(q[0] - want) < 1e-12, `got ${q[0]}`);
  let rest = 0;
  for (let k = 2; k < q.length; k++) rest = Math.max(rest, Math.abs(q[k]));
  check('constant field -> other coeffs ~ 0', rest < 1e-12, `max=${rest.toExponential(2)}`);
}

console.log(failures === 0 ? '\nALL NODE TESTS PASSED' : `\n${failures} TEST(S) FAILED`);
process.exit(failures === 0 ? 0 : 1);