Procedures

All unnormalized associated Legendre functions $P_\ell^m(x)$ for $0\le\ell\le\ell_{max},\,0\le m\le\ell$, via the Bonnet recurrence, in one pass (amortizing shared recurrence terms across every $\ell$ at fixed $m$ -- far cheaper than $(\ell_{max}{+}1)^2/2$ independent ASSOC_LEGENDRE calls). Condon-Shortley phase; same recurrence and same stability limit as ASSOC_LEGENDRE (see ASSOC_LEGENDRE_NORM_ALL for the far-more-stable normalized alternative, itself safe to $\ell\le2000$).

Private helper: $P_\ell^k(x)$ and $dP_\ell^k/dx$ together, from a single upward recurrence pass. Avoids the redundant separate ASSOC_LEGENDRE calls that DDX_ASSOC_LEGENDRE would otherwise need, and the repeated SSH + DDX_ASSOC_LEGENDRE evaluation that GRAD_SSH and L_SSH both build on.

Associated Legendre function $P_\ell^k(x)$, Condon-Shortley phase convention, computed by the classic Bonnet-type upward recurrence starting from $P_m^m(x) = (-1)^m(2m{-}1)!!\,(1{-}x^2)^{m/2}$. For negative order, $P_\ell^{-k}(x) = (-1)^k\,\frac{(\ell-k)!}{(\ell+k)!} P_\ell^k(x)$. Unnormalized (unlike ASSOC_LEGENDRE_NORM_ALL); loses numerical stability well before the intermediate $P_m^m$ values literally overflow, since the later upward-in-$\ell$ recurrence involves a subtraction of two comparably huge terms (catastrophic cancellation).

All $4\pi$-normalized associated Legendre functions $N_\ell^m P_\ell^m(\cos\theta)$, $N_\ell^m=\sqrt{\frac{2\ell+1} {4\pi}\frac{(\ell-m)!}{(\ell+m)!}}$ (i.e. the real-valued radial part of SSH, for $0\le\ell\le\ell_{max},\,0\le m\le\ell$), using the Holmes & Featherstone (2002) modified forward-column recurrence. The normalization is folded directly into the recurrence coefficients, keeping every intermediate value $O(1/\sqrt{4\pi})$ rather than letting $P_m^m$ itself grow/shrink combinatorially like the unnormalized ASSOC_LEGENDRE_ALL/Bonnet recurrence does. This is the recurrence SSH_ALL and every batch VSH routine build on.

Batch-consistency check: ASSOC_LEGENDRE_ALL vs. looped ASSOC_LEGENDRE calls, for $0\le\ell\le\ell_{max},\,0\le m\le\ell$, swept over NX points $x\in[-0.9,0.9]$. Writes every (l,m,x) triple's batch value, single-mode value, and absolute difference to OUTUNIT; fails (STATUS=1) if any relative error exceeds $10^{-12}$.

Batch-consistency check: DDX_ASSOC_LEGENDRE_ALL vs. looped DDX_ASSOC_LEGENDRE calls, for $0\le\ell\le\ell_{max},\,0\le m\le\ell$, swept over NX points $x\in[-0.9,0.9]$. Same tolerance/output convention as BATCH_ALM_CONS.

Batch-consistency check: GRAD_SSH_ALL vs. looped GRAD_SSH calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NTH colatitude points at fixed $\phi=\pi/4$. Fails (STATUS=1) if any component's absolute difference exceeds $10^{-12}$.

Batch-consistency check: L_SSH_ALL vs. looped L_SSH calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NTH colatitude points at fixed $\phi=\pi/4$. Fails (STATUS=1) if any component's absolute difference exceeds $10^{-12}$.

Batch-consistency check: PVSH_POL_ALL vs. looped PVSH_POL calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NTH colatitude points at fixed $\phi=\pi/4$. Fails (STATUS=1) if any component's absolute difference exceeds $10^{-12}$.

Batch-consistency check: PVSH_RAD_ALL vs. looped PVSH_RAD calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NTH colatitude points at fixed $\phi=\pi/4$. Fails (STATUS=1) if any component's absolute difference exceeds $10^{-12}$.

Batch-consistency check: PVSH_TOR_ALL vs. looped PVSH_TOR calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NTH colatitude points at fixed $\phi=\pi/4$ (both routines agree at $\ell=0$ too -- PVSH_TOR explicitly branches to return zero there rather than evaluating a 0/0 form). Fails (STATUS=1) if any component's absolute difference exceeds $10^{-12}$.

Batch-consistency check: SSH_ALL vs. looped SSH calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NTH colatitude points at fixed $\phi=\pi/4$. Fails (STATUS=1) if any $|Y_{batch}-Y_{single}|$ exceeds $10^{-13}$.

Batch-consistency check: VSH_POL_DN_ALL vs. looped VSH_POL_DN calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NTH colatitude points at fixed $\phi=\pi/4$. Fails (STATUS=1) if any component's absolute difference exceeds $10^{-12}$.

Batch-consistency check: VSH_POL_UP_ALL vs. looped VSH_POL_UP calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NTH colatitude points at fixed $\phi=\pi/4$. Fails (STATUS=1) if any component's absolute difference exceeds $10^{-12}$.

Batch-consistency check: VSH_TOR_ALL vs. looped VSH_TOR calls, for $1\le\ell\le\ell_{max},\,-\ell\le m\le\ell$ (starts at $\ell=1$; VSH_TOR is identical to PVSH_TOR, already covered at $\ell=0$ by BATCH_PVSH_TOR_CONS), swept over NTH colatitude points at fixed $\phi=\pi/4$. Fails (STATUS=1) if any component's absolute difference exceeds $10^{-12}$.

Benchmarks ASSOC_LEGENDRE_NORM_ALL (batch) against looped ASSOC_LEGENDRE (single-mode) calls, for $0\le\ell\le\ell_{max},\, 0\le m\le\ell$, swept over NX points $x\in[-0.9,0.9]$ (one grid sweep per timing, repeated/calibrated per WALL_TIME).

Benchmarks SSH_ALL (batch) against looped SSH (single-mode) calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NX colatitude points in $(0,\pi)$ at a fixed longitude.

Benchmarks VSH_POL_DN_ALL (batch) against looped VSH_POL_DN (single-mode) calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NX colatitude points.

Benchmarks VSH_POL_UP_ALL (batch) against looped VSH_POL_UP (single-mode) calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NX colatitude points.

Benchmarks VSH_TOR_ALL (batch) against looped VSH_TOR (single-mode) calls, for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, swept over NX colatitude points.

Clebsch-Gordan coefficient $\langle j_1 m_1 j_2 m_2 | j_3 m_3\rangle$ for coupling two angular momenta $j_1,j_2$ to a total $j_3$, via the closed-form Racah sum (evaluated in log-factorial form for numerical stability at large arguments -- see LOG_FACT): $$\langle j_1 m_1 j_2 m_2|j_3 m_3\rangle = \delta_{m_3,m_1+m_2} \sqrt{2j_3+1}\,\Delta(j_1j_2j_3)\times$$ $$\sqrt{(j_1{+}m_1)!(j_1{-}m_1)!(j_2{+}m_2)!(j_2{-}m_2)!(j_3{+}m_3)!(j_3{-}m_3)!} \sum_k \frac{(-1)^k}{k!\,(j_1{+}j_2{-}j_3{-}k)!(j_3{-}j_1{-}m_2{+}k)! (j_3{-}j_2{+}m_1{+}k)!(j_1{-}m_1{-}k)!(j_2{+}m_2{-}k)!},$$ where $\Delta(j_1j_2j_3) = \sqrt{\frac{(j_1+j_2-j_3)!(j_2+j_3-j_1)! (j_3+j_1-j_2)!}{(j_1+j_2+j_3+1)!}}$. Adapted from David Simpson (NASA GSFC). Automatically returns 0 for any invalid combination -- violated triangle inequality $|j_1-j_2|\le j_3\le j_1+j_2$, $|m_i|>j_i$, or $m_1+m_2\ne m_3$ -- so it is always safe to call with an arbitrary sextuple of integers, e.g. when sweeping a parameter range.

All derivatives $dP_\ell^m/dx$ for $0\le\ell\le\ell_{max},\, 0\le m\le\ell$, reusing a precomputed ASSOC_LEGENDRE_ALL table rather than recomputing $P_\ell^m$ from scratch.

Derivative $dP_\ell^k/dx$ of the associated Legendre function, via ASSOC_LEGENDRE_AND_DERIV (thin wrapper that discards $P_\ell^k(x)$ itself).

All derivatives $d(N_\ell^m P_\ell^m)/dx$ of the normalized associated Legendre functions, for $0\le\ell\le\ell_{max},\,0\le m\le\ell$, reusing a precomputed ASSOC_LEGENDRE_NORM_ALL table: $$\frac{d(N_\ell^mP_\ell^m)}{dx} = \frac{\sqrt{\frac{(2\ell+1) (\ell^2-m^2)}{2\ell-1}}\,N_{\ell-1}^mP_{\ell-1}^m - \ell\,x\, N_\ell^mP_\ell^m}{1-x^2}.$$ This is the recurrence VSH_CORE (and therefore every batch VSH routine) uses for the $\hat\theta$ component.

Derivative $dP_\ell/dx$ of the Legendre polynomial, via $(x^2{-}1)\,P_\ell'(x) = \ell\,[x P_\ell(x) - P_{\ell-1}(x)]$ away from the poles, and the closed-form pole values $P_\ell'(\pm1) = (\pm1)^{\ell+1}\,\ell(\ell{+}1)/2$ at $x=\pm1$.

Decomposes FIELD_FN into $(a_{rad},a_{pol},a_{tor})$ for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, via 2D Gauss-Legendre($\theta$) times uniform-trapezoidal($\phi$) quadrature: $$a_{rad}(\ell,m) = \sum_{quad} \mathbf{F}(\theta,\phi)\cdot \mathrm{[[PVSH_RAD]]}(\ell,m,\theta,\phi)^{*}\,w_\theta\,w_\phi$$ (and similarly for $a_{pol},a_{tor}$ against PVSH_POL and PVSH_TOR), using the batch _ALL routines so every mode at a quadrature point comes from one call.

Decomposes FIELD_FN into $(a_{rad},a_{pol},a_{tor})$ for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, via 2D Gauss-Legendre($\theta$) times uniform-trapezoidal($\phi$) quadrature. Identical to vsh_decomposition.f90's DECOMPOSE -- duplicated rather than shared, matching this repository's convention of fully independent worked-example programs (see CMakeLists.txt).

Dot product of two 3-vectors (typically two VSH values evaluated at the same $(\theta,\phi)$): $\mathbf{V}_1\cdot\mathbf{V}_2 = \sum_{c=1}^3 V_{1,c}V_{2,c}$. Note this is a plain component-wise product, not a Hermitian inner product -- callers wanting $\mathbf{V}_1\cdot\mathbf{V}_2^{*}$ must conjugate VSH2 themselves first (as done throughout src/tests.f90 and the examples/ drivers).

Factorial $k!$ for $k \ge 0$, computed by direct multiplication.

Target field for Part B: a smooth angular Gaussian bump with slight $\phi$ modulation (carries $m=0$ and $m=\pm1$ content, on top of effectively unbounded $\ell$ content from the Gaussian profile in $\theta$). Purely radial, for simplicity.

Target field for Part A: equal unit weight on all three basis families at each of the eight known modes; zero elsewhere. A genuine finite sum over eight independent modes, not a shortcut.

Target field for Part A: equal unit weight on PVSH_TOR alone at each of the eight known modes; zero elsewhere. Purely toroidal by construction (each term is a direct PVSH_TOR call).

Target field for Part B: an exact axisymmetric $(\ell,m)=(2,0)$ "quadrupole" toroidal background plus a broadband toroidal field built from $\mathbf{B}=\hat r\times\nabla_\perp T$ applied to an equator-centered Gaussian bump generating function (see module docstring for the closed-form derivatives and why $\sin^2\theta$, not $\sin\theta$, is needed to force $B_\theta$/$B_\phi$ to vanish exactly at both poles) -- purely toroidal by construction in both pieces.

All angular gradients $\nabla_\perp Y_\ell^m$ (GRAD_SSH) for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, from one VSH_CORE pass.

Angular (surface) gradient of the scalar spherical harmonic, $\nabla_\perp Y_\ell^k = \hat\theta\,\partial_\theta Y_\ell^k + \hat\phi\,\frac{1}{\sin\theta}\partial_\phi Y_\ell^k$. Returned as a 3-vector $(\hat r,\hat\theta,\hat\phi)$ with a zero radial component, since the gradient of a purely angular function is purely tangential.

Geppert-Wiebicke axisymmetric coupling coefficient $I_{j_1 m_1 j_2 m_2}^{\ell m}$ (Geppert & Wiebicke, 1991), a compact shorthand for the generation/amplification amplitude of a two-mode interaction in the Hall-MHD induction equation: $$I_{j_1 m_1 j_2 m_2}^{\ell m} = \sqrt{\frac{(2j_1+1)(2j_2+1)} {4\pi(2\ell+1)}}\,\langle j_1 0\,j_2 0|\ell 0\rangle\, \langle j_1 m_1\,j_2 m_2|\ell m\rangle.$$ This is algebraically the standard Gaunt coefficient -- the closed form for $\int Y_{j_1}^{m_1}Y_{j_2}^{m_2}(Y_\ell^m)^{*}\,d\Omega$ -- written in terms of two CGCOEFF evaluations. Built entirely from CGCOEFF, which guards all selection rules, so GWI is safe to call with any integer sextuple and returns 0 for invalid combinations (unlike GWJ, see its own selection-rule note).

Geppert-Wiebicke non-axisymmetric coupling coefficient $J_{j_1 m_1 j_2 m_2}^{\ell m}$ (Geppert & Wiebicke, 1991), the companion to GWI encoding the amplitude of a two-mode Hall-MHD interaction that generates or amplifies a non-axisymmetric component: $$J_{j_1 m_1 j_2 m_2}^{\ell m} = -\frac{i}{2} \sqrt{\frac{(2j_1+1)(2j_2+1)}{4\pi(2\ell+1)}}\, \sqrt{(j_1{+}j_2{+}\ell{+}2)(j_2{+}\ell{-}j_1)(j_1{+}j_2{-}\ell{+}1) (j_1{-}j_2{+}\ell{+}1)}\; \times$$ $$\langle j_1{+}1,0\,j_2 0|\ell 0\rangle\, \langle j_1 m_1\,j_2 m_2|\ell m\rangle.$$ Cross-validated against an independent VSH-dot-product expansion in TEST2_GW (src/tests.f90).

All angular-momentum-operator fields $\hat r\times\nabla_\perp Y_\ell^m$ (L_SSH) for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, from one VSH_CORE pass.

The angular-momentum-operator field $\hat r \times \nabla_\perp Y_\ell^k$ (proportional to the quantum-mechanical orbital angular-momentum operator $\mathbf{L}\,Y_\ell^k$ acting on the scalar harmonic). Simply GRAD_SSH rotated 90 degrees within the tangent plane: $\hat r\times(\hat\theta\,G_\theta+\hat\phi\,G_\phi) = \hat\theta\,(-G_\phi) + \hat\phi\,G_\theta$.

Legendre polynomial $P_\ell(x)$, degree $\ell \ge 0$, via the standard three-term recurrence $$(\ell{+}1)\,P_{\ell+1}(x) = (2\ell{+}1)\,x\,P_\ell(x) - \ell\,P_{\ell-1}(x).$$

$\ln(n!)$, via $\ln(n!) = \ln\Gamma(n{+}1)$ (the intrinsic LOG_GAMMA). Used throughout CGCOEFF, SSH, and ASSOC_LEGENDRE to evaluate factorial ratios like $\sqrt{(\ell-k)!/(\ell+k)!}$ as $\exp(\tfrac12[\texttt{LOG\_FACT}(\ell{-}k)-\texttt{LOG\_FACT}(\ell{+}k)])$ without overflowing the individual factorials for large $\ell$.

1D packed-triangular index for $(\ell,m)$ into the output arrays of ASSOC_LEGENDRE_ALL, DDX_ASSOC_LEGENDRE_ALL, ASSOC_LEGENDRE_NORM_ALL, and DDX_ASSOC_LEGENDRE_NORM_ALL (all of which only need $0\le m\le\ell$, by conjugate/parity symmetry). Matches the SHTOOLS PlmIndex convention.

All poloidal (horizontal) polar VSH members PVSH_POL for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, from one VSH_CORE pass. Zero for $\ell=0$.

Poloidal (horizontal) member of the polar vector spherical harmonic basis, $\mathbf{Y}_{\ell m}^{(+1)} = \frac{1}{\sqrt{\ell(\ell+1)}}\, \nabla_\perp Y_\ell^m$ -- purely horizontal (zero radial component), i.e. GRAD_SSH normalized by $1/\sqrt{\ell(\ell+1)}$. Combined with a radial stream function $f(r)$, this is the horizontal-field piece of a poloidal field, $\mathbf{B}_{\theta,\phi} = \frac{1}{r}\frac{d(rf)}{dr}\, \sqrt{\ell(\ell+1)}\,\mathbf{Y}_{\ell m}^{(+1)}$. Zero for $\ell=0$ (a constant has no horizontal gradient).

Verifies pointwise orthogonality of the poloidal and toroidal polar VSH members, $\mathbf{Y}_{\ell m}^{(+1)}\cdot\mathbf{Y}_{\ell m}^{(0)} = 0$ for every $\ell\ge1,\,m,\,\theta,\,\phi$ -- an exact analytic identity (poloidal is a pure gradient, toroidal a pure $\hat r\times$gradient, and $\nabla_\perp f\cdot\hat r\times \nabla_\perp f=0$ always), so any nonzero DOT value here is purely numerical error.

All radial polar VSH members PVSH_RAD for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, from one VSH_CORE pass.

Radial member of the polar vector spherical harmonic basis, $\mathbf{Y}_{\ell m}^{(-1)} = \hat r\,Y_\ell^m(\theta,\phi)$. Purely radial by construction (zero horizontal components) -- combined with a radial stream function $f(r)$, this is the $B_r$-generating piece of a poloidal magnetic (or any divergence-free) field, $B_r = \frac{\ell(\ell+1)}{r^2}\,f(r)\,Y_\ell^m$.

All toroidal polar VSH members PVSH_TOR for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, from one VSH_CORE pass. Zero for $\ell=0$.

Toroidal member of the polar vector spherical harmonic basis, $\mathbf{Y}_{\ell m}^{(0)} = \frac{1}{\sqrt{\ell(\ell+1)}}\, \hat r\times\nabla_\perp Y_\ell^m$ -- purely horizontal (zero radial component), identical in form to L_SSH but normalized by $1/\sqrt{\ell(\ell+1)}$. This is the horizontal-field piece of a toroidal magnetic field, $\mathbf{B}_{tor} = T(r)\, \mathbf{Y}_{\ell m}^{(0)}$ for a toroidal stream function $T(r)$. Zero for $\ell=0$ (no toroidal monopole).

Reconstructs the field at $(\theta,\phi)$ from coefficients up to LMAX: $\mathbf{F}\approx\sum_{\ell,m} a_{rad}\,\mathrm{[[PVSH_RAD]]}+a_{pol}\,\mathrm{[[PVSH_POL]]}+ a_{tor}\,\mathrm{[[PVSH_TOR]]}$, the inverse of DECOMPOSE.

Reconstructs the field at $(\theta,\phi)$ from coefficients up to LMAX. Identical to vsh_decomposition.f90's RECONSTRUCT.

GWI vs. the brute-force Gaunt-coefficient integral (via SHGLQ theta-quadrature times a uniform phi sum), swept over $j_1,j_2=0\ldots$ JMAX, all valid $m_1,m_2,\ell,m$. Writes examples/gwi_gwj_sweep/gwi_check.dat.

GWJ vs. the $\mathrm{[[PVSH_POL]]}\cdot\mathrm{[[PVSH_TOR]]}$ pointwise identity, swept over $j_1,j_2=0\ldots$ JMAX, all valid $m_1,m_2$, at N_POINTS fixed $(\theta,\phi)$. Writes examples/gwi_gwj_sweep/gwj_check.dat.

Part A: DECOMPOSE FIELD_KNOWN, confirm exact recovery of the eight known unit coefficients and near-zero everywhere else. Writes examples/vsh_decomposition/known_modes.dat.

Part A: DECOMPOSE FIELD_KNOWN_TOR, confirm exact recovery of the eight known unit a_tor coefficients and near-zero everywhere else, and confirm a_rad/a_pol are near-zero at every mode (the direct purity check). Writes examples/vsh_decomposition_tor/known_modes.dat.

Part B: DECOMPOSE+RECONSTRUCT FIELD_GAUSSIAN at increasing $\ell_{max}$ (quadrature order scaled with each: LMAX_QUAD = 3*LMAX_TEST+10, NPHI = 4*LMAX_TEST+16), tracking max reconstruction error against the true closed-form field on an independent evaluation grid. Self-check: error at the largest $\ell_{max}$ must be much smaller than at the smallest (genuine convergence, not a flat/noisy curve -- which would indicate a quadrature-order-vs-$\ell_{max}$ scaling bug). Writes examples/vsh_decomposition/convergence.dat.

Part B: DECOMPOSE+RECONSTRUCT FIELD_QUAD_BUMP at increasing $\ell_{max}$ (quadrature order scaled with each: LMAX_QUAD = 3*LMAX_TEST+30, NPHI = 4*LMAX_TEST+32 -- a higher floor than example 2d needs, since the equator-centered bump's sharper angular structure otherwise leaves a small quadrature-truncation leak in a_rad/a_pol at the smallest LMAX_TEST, shrinking with quadrature order rather than indicating an impure construction), tracking max reconstruction error against the true closed-form field on an independent evaluation grid, plus max a_rad/a_pol leakage across the sweep. Self-check: error at the largest $\ell_{max}$ must be much smaller than at the smallest (genuine convergence). Writes examples/vsh_decomposition_tor/convergence.dat.

Gauss-Legendre quadrature nodes and weights, exact for integrating any polynomial of degree $\le 2\ell_{max}{+}1$ using only $\ell_{max}{+}1$ points -- the natural quadrature for spherical harmonic transforms, since $P_\ell^m(\cos\theta)$ is a polynomial in $x=\cos\theta$ times $(1{-}x^2)^{|m|/2}$, and the transform integral separates exactly into this $x$-quadrature times a uniform sum over $\phi$ (itself exact for finite Fourier content). Nodes are the zeros of $P_{\ell_{max}+1}(x)$, found by Newton iteration (converges to machine precision in under 10 steps); weights are $w_i = 2/[(1-x_i^2)\,P_{\ell_{max}+1}'(x_i)^2]$.

All scalar spherical harmonics $Y_\ell^m(\theta,\phi)$ for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, in one call. Computes ASSOC_LEGENDRE_NORM_ALL once, then fills $m\ge0$ directly and $m<0$ via the conjugate symmetry $Y_\ell^{-m} = (-1)^m\,(Y_\ell^m)^{*}$ rather than a second recurrence pass.

Scalar spherical harmonic $Y_\ell^k(\theta,\phi)$, fully normalized with the Condon-Shortley phase: $$Y_\ell^k(\theta,\phi) = \sqrt{\frac{2\ell+1}{4\pi}\frac{(\ell-k)!}{(\ell+k)!}}\, P_\ell^k(\cos\theta)\,e^{ik\phi}.$$

Verifies scalar spherical harmonic orthonormality, $\int Y_{\ell_1}^m(Y_{\ell_2}^m)^{*}\,d\Omega = \delta_{\ell_1\ell_2}$, for every $0\le m\le\ell_1\le\ell_2\le\ell_{max}$. The $\phi$ integral is done analytically ($2\pi$ for matching $m$, by construction 0 otherwise, so only same-$m$ pairs are checked); the $\theta$ integral (in $u=\cos\theta$) uses a plain NQUAD-point midpoint rule via ASSOC_LEGENDRE_ALL, not SHGLQ -- sufficient here since only a scalar pass/fail diagonal check is needed, not a high-order transform.

Wigner 3-j symbol $\begin{pmatrix}j_1&j_2&j_3\\m_1&m_2&m_3\end{pmatrix}$, via the standard conversion from CGCOEFF: $$\begin{pmatrix}j_1&j_2&j_3\\m_1&m_2&m_3\end{pmatrix} = \frac{(-1)^{j_1-j_2-m_3}}{\sqrt{2j_3+1}}\, \langle j_1,-m_1,j_2,-m_2|j_3,-m_3\rangle.$$ Inherits CGCOEFF's selection-rule guards -- safe to call with any integer sextuple, returning 0 for invalid combinations.

Hand-derived closed form for Test 1's target quantity (a specific product of VSH inner products), evaluated over a colatitude grid at $\phi=0$ for comparison against TEST1_VSH and TEST1_GW.

Test 1's target quantity computed via the GWI-coefficient expansion in scalar spherical harmonics, rather than direct VSH DOT products -- must agree with TEST1_ANALYTIC and TEST1_VSH to machine precision. This is the check that specifically exercises GWI (see also TEST2_GW for the GWJ counterpart).

Test 1's target quantity computed directly from DOT products of VSH_TOR/VSH_POL_DN/VSH_POL_UP values -- must agree with TEST1_ANALYTIC to machine precision.

Hand-derived closed form for Test 2's target quantity (a different product of VSH inner products than Test 1, with nonzero $m$ content), evaluated over a colatitude grid at a fixed $\phi$ for comparison against TEST2_VSH and TEST2_GW.

Test 2's target quantity computed via the GWJ-coefficient expansion (combined with a GWI term) in scalar spherical harmonics -- must agree with TEST2_ANALYTIC and TEST2_VSH to machine precision. This is the check that specifically exercises GWJ.

Test 2's target quantity computed directly from DOT products of PVSH_POL/PVSH_TOR/PVSH_RAD values -- must agree with TEST2_ANALYTIC to machine precision.

Private helper shared by every batch VSH routine (GRAD_SSH_ALL, L_SSH_ALL, PVSH_RAD_ALL, PVSH_POL_ALL, PVSH_TOR_ALL, VSH_TOR_ALL, VSH_POL_UP_ALL, VSH_POL_DN_ALL): computes $Y_\ell^m$ and its two angular gradient components (GRAD_SSH's $\hat\theta,\hat\phi$ parts) for every $(\ell,m)$ in one pass, so each caller only has to apply its own linear combination/normalization on top.

All $J=\ell{-}1$ standard-VSH members VSH_POL_DN for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, from one VSH_CORE pass: $\sqrt{\ell/(2\ell{+}1)}\,\mathbf{Y}_{\ell m}^{(-1)} + \sqrt{(\ell{+}1)/(2\ell{+}1)}\,\mathbf{Y}_{\ell m}^{(+1)}$. Zero for $\ell=0$ -- see VSH_POL_DN for the dipole-ratio caveat at $\ell=1$.

The $J=\ell{-}1$ member of the standard (total-angular-momentum coupled) vector spherical harmonic basis, built from the polar basis by $$\mathbf{Y}_{\ell m}^{\ell-1} = \sqrt{\tfrac{\ell}{2\ell+1}}\, \mathbf{Y}_{\ell m}^{(-1)} + \sqrt{\tfrac{\ell+1}{2\ell+1}}\, \mathbf{Y}_{\ell m}^{(+1)},$$ i.e. PVSH_RAD and PVSH_POL rotated together by an $\ell$-dependent angle (see VSH_POL_UP for the companion $J=\ell{+}1$ rotation). Zero for $\ell=0$ (no $J=-1$).

Verifies the orthogonal-rotation identity linking the standard and polar VSH bases can be inverted exactly: $$\sqrt{\tfrac{\ell}{2\ell+1}}\,\mathbf{Y}_{\ell m}^{\ell+1} + \sqrt{\tfrac{\ell+1}{2\ell+1}}\,\mathbf{Y}_{\ell m}^{\ell-1} = \mathbf{Y}_{\ell m}^{(+1)},$$ $$-\sqrt{\tfrac{\ell+1}{2\ell+1}}\,\mathbf{Y}_{\ell m}^{\ell+1} + \sqrt{\tfrac{\ell}{2\ell+1}}\,\mathbf{Y}_{\ell m}^{\ell-1} = \mathbf{Y}_{\ell m}^{(-1)},$$ i.e. rotating VSH_POL_UP/VSH_POL_DN back by the same angle used to build them from PVSH_POL/PVSH_RAD recovers the originals exactly.

All $J=\ell{+}1$ standard-VSH members VSH_POL_UP for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$, from one VSH_CORE pass: $-\sqrt{(\ell{+}1)/(2\ell{+}1)}\,\mathbf{Y}_{\ell m}^{(-1)} + \sqrt{\ell/(2\ell{+}1)}\,\mathbf{Y}_{\ell m}^{(+1)}$. Well-defined (nonzero) at $\ell=0$ -- see VSH_POL_UP for the dipole-ratio note at $\ell=1$.

The $J=\ell{+}1$ member of the standard (total-angular-momentum coupled) vector spherical harmonic basis, built from the polar basis by $$\mathbf{Y}_{\ell m}^{\ell+1} = \sqrt{\tfrac{\ell+1}{2\ell+1}}\, \mathbf{Y}_{\ell m}^{(-1)} - \sqrt{\tfrac{\ell}{2\ell+1}}\, \mathbf{Y}_{\ell m}^{(+1)},$$ the companion rotation to VSH_POL_DN's $J=\ell{-}1$. Unlike VSH_POL_DN this is well-defined (nonzero radial part) at $\ell=0$, since $J=1$ still exists when $\ell=0$.

All toroidal standard-VSH members VSH_TOR for $0\le\ell\le\ell_{max},\,-\ell\le m\le\ell$. Identical to PVSH_TOR_ALL (thin wrapper) -- the toroidal member is common to both bases.

Toroidal member of the "standard" (total-angular-momentum-coupled) vector spherical harmonic basis, $\mathbf{Y}_{\ell m}^{J=\ell}$. Identical to PVSH_TOR -- the toroidal member is common to both the polar and standard bases, since $\hat r\times\nabla_\perp Y_\ell^m$ already carries total angular momentum $J=\ell$ with no $J=\ell\pm1$ admixture.

Wall-clock seconds elapsed, via SYSTEM_CLOCK with an INTEGER(i8) count so the calibration loops in every BENCH_* routine below cannot roll the counter over during a long-running sweep.

1D index for $(\ell,m)$ into the output arrays of every batch _ALL routine that needs the full $-\ell\le m\le\ell$ range (SSH_ALL, GRAD_SSH_ALL, L_SSH_ALL, and every batch VSH routine).