1/**
2 * Recurrence coefficients for orthonormal associated Legendre functions
3 * ytilde_l^m(theta) (spherical-harmonic normalized, Condon-Shortley phase
4 * included), matching SHTNS legendre_precomp() with norm=sht_orthonormal:
5 *
6 * ytilde_m^m(theta) = amm * sin(theta)^m
7 * ytilde_{m+1}^m = a_{m+1}^m * cos(theta) * ytilde_m^m
8 * ytilde_l^m = a_l^m * cos(theta) * ytilde_{l-1}^m + b_l^m * ytilde_{l-2}^m
9 *
10 * with (cf. sht_legendre.c lines 442-447):
11 * a_{m+1}^m = sqrt(2m+3)
12 * a_l^m = sqrt( (2l+1)(2l-1) / ((l+m)(l-m)) )
13 * b_l^m = -sqrt( (2l+1)/(2l-3) * ((l-1+m)(l-1-m)) / ((l+m)(l-m)) )
14 * amm = cs^m * sqrt( 1/(4pi) * prod_{k=1..m} (2k+1)/(2k) )
15 *
16 * With this normalization, Y_lm(theta,phi) = ytilde_l^m(theta) e^{i m phi}
17 * and integral |Y_lm|^2 dOmega = 1.
18 */
19import { lmIndex, nlmCalc } from './layout.ts';
21export interface LegendreCoeffs {
22 /** amm[m]: seed value (includes Condon-Shortley phase (-1)^m). */
23 amm: Float64Array;
24 /** ab[2*lm], ab[2*lm+1] = (a_l^m, b_l^m); entries at l=m unused (0), b at l=m+1 unused (0). */
25 ab: Float64Array;
26}
28export function legendreCoeffs(lmax: number, mmax: number): LegendreCoeffs {
29 const nlm = nlmCalc(lmax, mmax);
30 const amm = new Float64Array(mmax + 1);
31 const ab = new Float64Array(2 * nlm);
33 let t = 1.0 / (4.0 * Math.PI);
34 amm[0] = Math.sqrt(t);
35 for (let m = 1; m <= mmax; m++) {
36 t *= (2 * m + 1) / (2 * m);
37 amm[m] = -Math.sqrt(t); // (-1)^m accumulates: Condon-Shortley phase
38 if (m % 2 === 0) amm[m] = -amm[m];
39 }
41 for (let m = 0; m <= mmax; m++) {
42 if (m + 1 <= lmax) {
43 const lm = lmIndex(lmax, m + 1, m);
44 ab[2 * lm] = Math.sqrt(2 * m + 3); // a_{m+1}^m
45 ab[2 * lm + 1] = 0;
46 }
47 for (let l = m + 2; l <= lmax; l++) {
48 const lm = lmIndex(lmax, l, m);
49 const t1 = (l + m) * (l - m);
50 const t2 = (l - 1 + m) * (l - 1 - m);
51 ab[2 * lm] = Math.sqrt(((2 * l + 1) * (2 * l - 1)) / t1);
52 ab[2 * lm + 1] = -Math.sqrt(((2 * l + 1) / (2 * l - 3)) * (t2 / t1));
53 }
54 }
55 return { amm, ab };
56}
58/**
59 * Evaluate ytilde_l^m(theta) for l = m..lmax at one point, in f64.
60 * ct = cos(theta), st = sin(theta). Plain (unscaled) recurrence: fine in
61 * f64 for the moderate lmax this library targets (underflow of st^m only
62 * matters for m of several hundred very close to the poles).
63 */
64export function legendreRow(
65 coeffs: LegendreCoeffs,
66 lmax: number,
67 m: number,
68 ct: number,
69 st: number,
70 out: Float64Array, // length lmax - m + 1
71): void {
72 let y0 = coeffs.amm[m] * Math.pow(st, m);
73 out[0] = y0;
74 if (m === lmax) return;
75 const base = lmIndex(lmax, m, m);
76 let y1 = coeffs.ab[2 * (base + 1)] * ct * y0;
77 out[1] = y1;
78 for (let l = m + 2; l <= lmax; l++) {
79 const lm = base + (l - m);
80 const y2 = coeffs.ab[2 * lm] * ct * y1 + coeffs.ab[2 * lm + 1] * y0;
81 y0 = y1;
82 y1 = y2;
83 out[l - m] = y2;
84 }
85}