1/**
2 * Recurrence coefficients for the first theta-derivative of orthonormal
3 * associated Legendre functions (Condon-Shortley phase included), matching
4 * the alpha^+/alpha^- recurrence in evolving_surface/notes/algos.tex Sec 2.1:
5 *
6 * sin(theta) d/dtheta Y_l^m = alpha^+(l,m) Y_{l+1}^m + alpha^-(l,m) Y_{l-1}^m
7 *
8 * so the coefficients of sin(theta)*dtheta(u), by degree, are
9 *
10 * v_l^m = alpha^+(l-1,m) u_{l-1}^m + alpha^-(l+1,m) u_{l+1}^m
11 *
12 * dropping any term referring to a degree outside 0 <= l <= lmax. Baked to
13 * zero at each m-block's first/last element (rather than left undefined), so
14 * a consuming WGSL kernel needs only an in-bounds check, not a validity check.
15 */
16import { lmIndex, nlmCalc } from './layout.ts';
18export interface DerivCoeffs {
19 /** aPlus[lm] = alpha^+(l-1,m) when l>m, else 0 -- multiplies u_{l-1}^m. */
20 aPlus: Float64Array;
21 /** aMinus[lm] = alpha^-(l+1,m) when l<lmax, else 0 -- multiplies u_{l+1}^m. */
22 aMinus: Float64Array;
23 /** m of the coefficient at flat index lm (the phi-derivative needs only this). */
24 mOf: Uint32Array;
25}
27export function alphaPlus(l: number, m: number): number {
28 return l * Math.sqrt(((l - m + 1) * (l + m + 1)) / ((2 * l + 1) * (2 * l + 3)));
29}
31export function alphaMinus(l: number, m: number): number {
32 return -(l + 1) * Math.sqrt(((l - m) * (l + m)) / ((2 * l - 1) * (2 * l + 1)));
33}
35export function derivCoeffs(lmax: number, mmax: number): DerivCoeffs {
36 const nlm = nlmCalc(lmax, mmax);
37 const aPlus = new Float64Array(nlm);
38 const aMinus = new Float64Array(nlm);
39 const mOf = new Uint32Array(nlm);
40 for (let m = 0; m <= mmax; m++) {
41 for (let l = m; l <= lmax; l++) {
42 const lm = lmIndex(lmax, l, m);
43 mOf[lm] = m;
44 if (l - 1 >= m) aPlus[lm] = alphaPlus(l - 1, m);
45 if (l + 1 <= lmax) aMinus[lm] = alphaMinus(l + 1, m);
46 }
47 }
48 return { aPlus, aMinus, mOf };
49}