VSH Module

private function ASSOC_LEGENDRE_DP(L, K, X) result(ASSOC_LEGENDRE)

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).

Returns: $P_\ell^k(x)$.

Arguments

Return Value real(kind=dp)

private subroutine ASSOC_LEGENDRE_ALL_DP(P, LMAX, X)

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$).

PLM_INDEX(l,m).

Arguments

private subroutine ASSOC_LEGENDRE_NORM_ALL_DP(PNORM, LMAX, X)

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.

PLM_INDEX(l,m).

Arguments

private function CGCOEFF_DP(J1, M1, J2, M2, J3, M3) result(CGCOEFF)

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.

Returns: $\langle j_1 m_1 j_2 m_2|j_3 m_3\rangle$, real (0 if selection rules are violated).

Arguments

Return Value real(kind=dp)

private function DDX_ASSOC_LEGENDRE_DP(L, K, X) result(RES)

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

Returns: $dP_\ell^k/dx$ evaluated at X; 0 at the poles ($|x|=1$) or for $|k|>\ell$.

Arguments

Return Value real(kind=dp)

private subroutine DDX_ASSOC_LEGENDRE_ALL_DP(DP_OUT, P, LMAX, X)

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.

by PLM_INDEX(l,m); 0 at the poles ($|x|\ge1$). ASSOC_LEGENDRE_ALL, same size/indexing. P).

Arguments

private subroutine DDX_ASSOC_LEGENDRE_NORM_ALL_DP(DPNORM, PNORM, LMAX, X)

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.

by PLM_INDEX(l,m); 0 at the poles ($|x|\ge1$). ASSOC_LEGENDRE_NORM_ALL, same size/indexing. PNORM).

Arguments

private function DDX_LEGENDRE_DP(L, X) result(DDX_LEGENDRE)

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$.

Returns: $dP_\ell/dx$ evaluated at X.

Arguments

Return Value real(kind=dp)

private function DOT_DP(VSH1, VSH2) result(DOT)

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).

Returns: $\mathbf{VSH1}\cdot\mathbf{VSH2}$, complex scalar.

Arguments

Return Value complex(kind=dp)

private function GRAD_SSH_DP(L, K, THETA, PHI) result(GRAD_SSH)

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.

Returns: $(0,\,\partial_\theta Y_\ell^k,\,\frac{1}{\sin\theta} \partial_\phi Y_\ell^k)$, complex 3-vector.

Arguments

Return Value complex(kind=dp), DIMENSION(3)

private subroutine GRAD_SSH_ALL_DP(OUT, LMAX, THETA, PHI)

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.

1=$\hat r$ (always 0), 2=$\hat\theta$, 3=$\hat\phi$; mode axis indexed by YLM_INDEX(l,m).

Arguments

private function GWI_DP(J1, M1, J2, M2, L, M) result(GWI)

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).

Returns: $I_{j_1 m_1 j_2 m_2}^{\ell m}$, complex (real-valued in practice, but returned as COMPLEX for a uniform interface with GWJ).

Arguments

Return Value complex(kind=dp)

private function GWJ_DP(J1, M1, J2, M2, L, M) result(GWJ)

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).

Returns: $J_{j_1 m_1 j_2 m_2}^{\ell m}$, complex.

Arguments

Return Value complex(kind=dp)

private function LEGENDRE_DP(L, X) result(LEGENDRE)

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).$$

Returns: $P_\ell(x)$.

Arguments

Return Value real(kind=dp)

private function LOG_FACT_DP(N) result(LOG_FACT)

$\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$.

Returns: $\ln(n!)$.

Arguments

Return Value real(kind=dp)

private function L_SSH_DP(L, K, THETA, PHI) result(L_SSH)

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$.

Returns: $(0,\,-\frac{1}{\sin\theta}\partial_\phi Y_\ell^k,\, \partial_\theta Y_\ell^k)$, complex 3-vector.

Arguments

Return Value complex(kind=dp), DIMENSION(3)

private subroutine L_SSH_ALL_DP(OUT, LMAX, THETA, PHI)

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.

1=$\hat r$ (always 0), 2=$\hat\theta$, 3=$\hat\phi$; mode axis indexed by YLM_INDEX(l,m).

Arguments

private function PVSH_POL_DP(L, K, THETA, PHI) result(PVSH_POL)

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).

Returns: $(0,\,\frac{1}{\sqrt{\ell(\ell+1)}}\partial_\theta Y_\ell^m,\, \frac{1}{\sqrt{\ell(\ell+1)}}\frac{1}{\sin\theta}\partial_\phi Y_\ell^m)$, complex 3-vector.

Arguments

Return Value complex(kind=dp), DIMENSION(3)

private subroutine PVSH_POL_ALL_DP(OUT, LMAX, THETA, PHI)

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$.

1=$\hat r$ (always 0), 2=$\hat\theta$, 3=$\hat\phi$; mode axis indexed by YLM_INDEX(l,m).

Arguments

private function PVSH_RAD_DP(L, K, THETA, PHI) result(PVSH_RAD)

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$.

Returns: $(Y_\ell^m,\,0,\,0)$, complex 3-vector.

Arguments

Return Value complex(kind=dp), DIMENSION(3)

private subroutine PVSH_RAD_ALL_DP(OUT, LMAX, THETA, PHI)

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

1=$\hat r$=$Y_\ell^m$, 2=$\hat\theta$ (always 0), 3=$\hat\phi$ (always 0); mode axis indexed by YLM_INDEX(l,m).

Arguments

private function PVSH_TOR_DP(L, K, THETA, PHI) result(PVSH_TOR)

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).

Returns: $(0,\,\frac{i}{\sqrt{\ell(\ell+1)}}\frac{1}{\sin\theta} \partial_\phi Y_\ell^m,\,-\frac{i}{\sqrt{\ell(\ell+1)}}\partial_\theta Y_\ell^m)$, complex 3-vector.

Arguments

Return Value complex(kind=dp), DIMENSION(3)

private subroutine PVSH_TOR_ALL_DP(OUT, LMAX, THETA, PHI)

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$.

1=$\hat r$ (always 0), 2=$\hat\theta$, 3=$\hat\phi$; mode axis indexed by YLM_INDEX(l,m).

Arguments

private subroutine SHGLQ_DP(ZERO, W, LMAX)

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]$.

ascending from -1 to +1. \sum_i w_i f(x_i) ), size LMAX+1. uses LMAX+1 quadrature points.

Arguments

private function SSH_DP(L, K, THETA, PHI) result(SSH)

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}.$$

Returns: $Y_\ell^k(\theta,\phi)$, complex.

Arguments

Return Value complex(kind=dp)

private subroutine SSH_ALL_DP(YLM, LMAX, THETA, PHI)

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.

Arguments

private function SYMBOL3J_DP(J1, M1, J2, M2, J3, M3) result(SYMBOL3J)

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.

Returns: The Wigner 3-j symbol, real (0 if selection rules are violated).

Arguments

Return Value real(kind=dp)

private function VSH_POL_DN_DP(L, K, THETA, PHI) result(VSH_POL_DN)

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$).

Returns: $\mathbf{Y}_{\ell m}^{\ell-1}(\theta,\phi)$, complex 3-vector.

Arguments

Return Value complex(kind=dp), DIMENSION(3)

private subroutine VSH_POL_DN_ALL_DP(OUT, LMAX, THETA, PHI)

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$.

1=$\hat r$, 2=$\hat\theta$, 3=$\hat\phi$; mode axis indexed by YLM_INDEX(l,m).

Arguments

private function VSH_POL_UP_DP(L, K, THETA, PHI) result(VSH_POL_UP)

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$.

Returns: $\mathbf{Y}_{\ell m}^{\ell+1}(\theta,\phi)$, complex 3-vector.

Arguments

Return Value complex(kind=dp), DIMENSION(3)

private subroutine VSH_POL_UP_ALL_DP(OUT, LMAX, THETA, PHI)

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$.

1=$\hat r$, 2=$\hat\theta$, 3=$\hat\phi$; mode axis indexed by YLM_INDEX(l,m).

Arguments

private function VSH_TOR_DP(L, K, THETA, PHI) result(VSH_TOR)

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.

Returns: $\mathbf{Y}_{\ell m}^{J=\ell}(\theta,\phi)$, complex 3-vector (identical to PVSH_TOR).

Arguments

Return Value complex(kind=dp), DIMENSION(3)

private subroutine VSH_TOR_ALL_DP(OUT, LMAX, THETA, PHI)

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.

1=$\hat r$ (always 0), 2=$\hat\theta$, 3=$\hat\phi$; mode axis indexed by YLM_INDEX(l,m).

Arguments

private subroutine ASSOC_LEGENDRE_AND_DERIV_DP(L, K, X, PLK, DPLK)

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.

Arguments

private function FACTORIAL_DP(K) result(FACTORIAL)

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

Returns: $k!$ as a double-precision real (avoids integer overflow for moderate K; see LOG_FACT for large-K stability).

Arguments

Return Value real(kind=dp)

private subroutine VSH_CORE_DP(YLM_OUT, GSSH_TH, GSSH_PH, LMAX, THETA, PHI)

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.

YLM_INDEX(l,m). Y_\ell^m ), same size/indexing. Y_\ell^m ), same size/indexing.

Arguments