1/**
2 * Node-side tests for the double-precision planner and reference transform.
3 * Run: node scripts/test-node.ts
4 */
5import { gaussNodesWeights } from '../src/gauss.ts';
6import { legendreCoeffs, legendreRow } from '../src/coeffs.ts';
7import { lmIndex } from '../src/layout.ts';
8import { ShtReference, randomSpectrum } from '../src/reference.ts';
10let failures = 0;
11function check(name: string, ok: boolean, detail = '') {
12 console.log(`${ok ? 'PASS' : 'FAIL'} ${name}${detail ? ' (' + detail + ')' : ''}`);
13 if (!ok) failures++;
14}
16// --- Gauss quadrature ---
17{
18 const { x, w } = gaussNodesWeights(64);
19 let sw = 0, sx2 = 0;
20 for (let i = 0; i < 64; i++) {
21 sw += w[i];
22 sx2 += w[i] * x[i] * x[i];
23 }
24 check('gauss: sum(w) == 2', Math.abs(sw - 2) < 1e-13, `err=${Math.abs(sw - 2).toExponential(2)}`);
25 check('gauss: int x^2 == 2/3', Math.abs(sx2 - 2 / 3) < 1e-13, `err=${Math.abs(sx2 - 2 / 3).toExponential(2)}`);
26 check('gauss: nodes decreasing', x[0] > x[1] && x[1] > x[63]);
27}
29// --- Analytic Legendre values (orthonormal + Condon-Shortley) ---
30{
31 const lmax = 8;
32 const coeffs = legendreCoeffs(lmax, lmax);
33 const theta = 0.7;
34 const ct = Math.cos(theta), st = Math.sin(theta);
35 const row = new Float64Array(lmax + 1);
37 legendreRow(coeffs, lmax, 0, ct, st, row);
38 const y00 = Math.sqrt(1 / (4 * Math.PI));
39 const y10 = Math.sqrt(3 / (4 * Math.PI)) * ct;
40 const y20 = Math.sqrt(5 / (16 * Math.PI)) * (3 * ct * ct - 1);
41 check('Y_0^0', Math.abs(row[0] - y00) < 1e-14);
42 check('Y_1^0', Math.abs(row[1] - y10) < 1e-14);
43 check('Y_2^0', Math.abs(row[2] - y20) < 1e-14);
45 legendreRow(coeffs, lmax, 1, ct, st, row);
46 const y11 = -Math.sqrt(3 / (8 * Math.PI)) * st; // CS phase => negative
47 const y21 = -Math.sqrt(15 / (8 * Math.PI)) * st * ct;
48 check('Y_1^1 (CS phase)', Math.abs(row[0] - y11) < 1e-14, `got ${row[0]}, want ${y11}`);
49 check('Y_2^1', Math.abs(row[1] - y21) < 1e-14);
51 legendreRow(coeffs, lmax, 2, ct, st, row);
52 const y22 = Math.sqrt(15 / (32 * Math.PI)) * st * st;
53 check('Y_2^2', Math.abs(row[0] - y22) < 1e-14);
54}
56// --- Orthonormality under Gauss quadrature ---
57{
58 const lmax = 42, nlat = 48;
59 const coeffs = legendreCoeffs(lmax, lmax);
60 const { x, w } = gaussNodesWeights(nlat);
61 let worst = 0;
62 for (const m of [0, 1, 7, 25]) {
63 const rowsI = new Float64Array(lmax + 1);
64 const rowsJ = new Float64Array(lmax + 1);
65 for (const [la, lb] of [[m, m], [m, m + 3], [lmax, lmax], [m + 1, lmax]] as const) {
66 if (la > lmax || lb > lmax) continue;
67 let s = 0;
68 for (let i = 0; i < nlat; i++) {
69 const st = Math.sqrt(1 - x[i] * x[i]);
70 legendreRow(coeffs, lmax, m, x[i], st, rowsI);
71 legendreRow(coeffs, lmax, m, x[i], st, rowsJ);
72 s += w[i] * rowsI[la - m] * rowsJ[lb - m];
73 }
74 const want = la === lb ? 1 / (2 * Math.PI) : 0;
75 worst = Math.max(worst, Math.abs(s - want));
76 }
77 }
78 check('orthonormality: max err < 1e-12', worst < 1e-12, `worst=${worst.toExponential(2)}`);
79}
81// --- Reference round trip ---
82{
83 const cfg = { lmax: 31, mmax: 31, nlat: 34, nphi: 64 };
84 const ref = new ShtReference(cfg);
85 const q0 = randomSpectrum(cfg, 999);
86 const spat = ref.synth(q0);
87 const q1 = ref.analys(spat);
88 let num = 0, den = 0;
89 for (let k = 0; k < q0.length; k++) {
90 num += (q1[k] - q0[k]) ** 2;
91 den += q0[k] ** 2;
92 }
93 const rel = Math.sqrt(num / den);
94 check('reference round trip rel L2 < 1e-12', rel < 1e-12, `rel=${rel.toExponential(2)}`);
95}
97// --- Mean value: Y_00 coefficient of a constant field ---
98{
99 const cfg = { lmax: 15, mmax: 15, nlat: 16, nphi: 32 };
100 const ref = new ShtReference(cfg);
101 const spat = new Float64Array(cfg.nlat * cfg.nphi).fill(1.0);
102 const q = ref.analys(spat);
103 const want = Math.sqrt(4 * Math.PI); // <1, Y00> = sqrt(4pi)
104 check('constant field -> Q_00 = sqrt(4pi)', Math.abs(q[0] - want) < 1e-12, `got ${q[0]}`);
105 let rest = 0;
106 for (let k = 2; k < q.length; k++) rest = Math.max(rest, Math.abs(q[k]));
107 check('constant field -> other coeffs ~ 0', rest < 1e-12, `max=${rest.toExponential(2)}`);
108}
110console.log(failures === 0 ? '\nALL NODE TESTS PASSED' : `\n${failures} TEST(S) FAILED`);
111process.exit(failures === 0 ? 0 : 1);