#!/usr/bin/env python
"""
associated_legendre.py
Written by Tyler Sutterley (03/2023)
Computes fully-normalized associated Legendre Polynomials
UPDATE HISTORY:
Updated 03/2023: improve typing for variables in docstrings
Updated 01/2023: refactored associated legendre polynomials
Updated 04/2022: updated docstrings to numpy documentation format
Updated 05/2021: define int/float precision to prevent deprecation warning
Updated 09/2020: verify dimensions of input x variable
Updated 08/2020: prevent zero divisions by changing u==0 to eps of data type
Updated 07/2020: added function docstrings
Updated 10/2018: using future division for python3 Compatibility
Updated 07/2017: output first differential of legendre polynomials
Updated 05/2015: added parameter MMAX for MMAX != LMAX
Updated 09/2013: new format for file headers
Written 03/2013
"""
from __future__ import division
import numpy as np
[docs]
def associated_legendre(LMAX, x, method='holmes', MMAX=None, astype=np.float64):
"""
Computes fully-normalized associated Legendre Polynomials and their
first derivative
Parameters
----------
LMAX: int
maximum degree of Legendre polynomials
x: np.ndarray
elements ranging from -1 to 1
Typically ``cos(theta)``, where ``theta`` is the colatitude in radians
method: str, default 'holmes'
Method for computing the associated Legendre polynomials
- ``'columbo'``
- ``'holmes'``
- ``'mohlenkamp'``
MMAX: int or NoneType, default None
maximum order of Associated Legendre polynomials
astype: np.dtype, default np.float64
output variable data type
Returns
-------
plms: np.ndarray
fully-normalized Legendre polynomials
dplms: np.ndarray
first derivative of Legendre polynomials
"""
if method.lower() == 'colombo':
return plm_colombo(LMAX, x, MMAX=MMAX, astype=astype)
elif method.lower() == 'holmes':
return plm_holmes(LMAX, x, MMAX=MMAX, astype=astype)
elif method.lower() == 'mohlenkamp':
return plm_mohlenkamp(LMAX, x, MMAX=MMAX, astype=astype)
raise ValueError(f'Unknown method {method}')
[docs]
def plm_colombo(LMAX, x, MMAX=None, astype=np.float64):
"""
Computes fully-normalized associated Legendre Polynomials and their
first derivative using a Standard forward column method :cite:p:`Colombo:1981vh`
Parameters
----------
LMAX: int
maximum degree of Legendre polynomials
x: np.ndarray
elements ranging from -1 to 1
Typically ``cos(theta)``, where ``theta`` is the colatitude in radians
MMAX: int or NoneType, default None
maximum order of Associated Legendre polynomials
astype: np.dtype, default np.float64
output variable data type
Returns
-------
plm: np.ndarray
fully-normalized Legendre polynomials
dplm: np.ndarray
first derivative of Legendre polynomials
"""
# removing singleton dimensions of x
x = np.atleast_1d(x).flatten().astype(astype)
# length of the x array
jm = len(x)
# verify data type of spherical harmonic truncation
LMAX = np.int64(LMAX)
# upper bound of spherical harmonic orders (default = LMAX)
if MMAX is None:
MMAX = np.copy(LMAX)
# allocating for the plm matrix and differentials
plm = np.zeros((LMAX + 1, LMAX + 1, jm))
dplm = np.zeros((LMAX + 1, LMAX + 1, jm))
# u is sine of colatitude (cosine of latitude) so that 0 <= s <= 1
# for x=cos(th): u=sin(th)
u = np.sqrt(1.0 - x**2)
# update where u==0 to eps of data type to prevent invalid divisions
u[u == 0] = np.finfo(u.dtype).eps
# Calculating the initial polynomials for the recursion
plm[0, 0, :] = 1.0
plm[1, 0, :] = np.sqrt(3.0) * x
plm[1, 1, :] = np.sqrt(3.0) * u
# calculating first derivatives for harmonics of degree 1
dplm[1, 0, :] = (1.0 / u) * (x * plm[1, 0, :] - np.sqrt(3) * plm[0, 0, :])
dplm[1, 1, :] = (x / u) * plm[1, 1, :]
for l in range(2, LMAX + 1):
for m in range(0, l): # Zonal and Tesseral harmonics (non-sectorial)
# Computes the non-sectorial terms from previously computed
# sectorial terms.
alm = np.sqrt(
((2.0 * l - 1.0) * (2.0 * l + 1.0)) / ((l - m) * (l + m))
)
blm = np.sqrt(
((2.0 * l + 1.0) * (l + m - 1.0) * (l - m - 1.0))
/ ((l - m) * (l + m) * (2.0 * l - 3.0))
)
# if (m == l-1): plm[l-2,m,:] will be 0
plm[l, m, :] = alm * x * plm[l - 1, m, :] - blm * plm[l - 2, m, :]
# calculate first derivatives
flm = np.sqrt(
((l**2.0 - m**2.0) * (2.0 * l + 1.0)) / (2.0 * l - 1.0)
)
dplm[l, m, :] = (1.0 / u) * (
l * x * plm[l, m, :] - flm * plm[l - 1, m, :]
)
# Sectorial harmonics
# The sectorial harmonics serve as seed values for the recursion
# starting with P00 and P11 (outside the loop)
plm[l, l, :] = (
u
* np.sqrt((2.0 * l + 1.0) / (2.0 * l))
* np.squeeze(plm[l - 1, l - 1, :])
)
# calculate first derivatives for sectorial harmonics
dplm[l, l, :] = np.longdouble(l) * (x / u) * plm[l, l, :]
# return the legendre polynomials and their first derivative
# truncating orders to MMAX
return plm[:, : MMAX + 1, :], dplm[:, : MMAX + 1, :]
[docs]
def plm_holmes(LMAX, x, MMAX=None, astype=np.float64):
"""
Computes fully-normalized associated Legendre Polynomials and their
first derivative using the recursion relation from :cite:p:`Holmes:2002ff`
Parameters
----------
LMAX: int
maximum degree of Legendre polynomials
x: np.ndarray
elements ranging from -1 to 1
Typically ``cos(theta)``, where ``theta`` is the colatitude in radians
MMAX: int or NoneType, default None
maximum order of Associated Legendre polynomials
astype: np.dtype, default np.float64
output variable data type
Returns
-------
plm: np.ndarray
fully-normalized Legendre polynomials
dplm: np.ndarray
first derivative of Legendre polynomials
"""
# removing singleton dimensions of x
x = np.atleast_1d(x).flatten().astype(astype)
# length of the x array
jm = len(x)
# verify data type of spherical harmonic truncation
LMAX = np.int64(LMAX)
# upper bound of spherical harmonic orders (default = LMAX)
if MMAX is None:
MMAX = np.copy(LMAX)
# scaling factor
scalef = 1.0e-280
# allocate for multiplicative factors, and plms
f1 = np.zeros(((LMAX + 1) * (LMAX + 2) // 2), dtype=astype)
f2 = np.zeros(((LMAX + 1) * (LMAX + 2) // 2), dtype=astype)
p = np.zeros(((LMAX + 1) * (LMAX + 2) // 2, jm), dtype=astype)
plm = np.zeros((LMAX + 1, LMAX + 1, jm), dtype=astype)
dplm = np.zeros((LMAX + 1, LMAX + 1, jm), dtype=astype)
# Precompute multiplicative factors used in recursion relationships
# Note that prefactors are not used for the case when m=l and m=l-1,
# as a different recursion is used for these two values.
k = 2 # k = l*(l+1)/2 + m
for l in range(2, LMAX + 1):
k += 1
f1[k] = (
np.sqrt(2.0 * l - 1.0) * np.sqrt(2.0 * l + 1.0) / np.longdouble(l)
)
f2[k] = (
np.longdouble(l - 1.0)
* np.sqrt(2.0 * l + 1.0)
/ (np.sqrt(2.0 * l - 3.0) * np.longdouble(l))
)
for m in range(1, l - 1):
k += 1
f1[k] = (
np.sqrt(2.0 * l + 1.0)
* np.sqrt(2.0 * l - 1.0)
/ (np.sqrt(l + m) * np.sqrt(l - m))
)
f2[k] = (
np.sqrt(2.0 * l + 1.0)
* np.sqrt(l - m - 1.0)
* np.sqrt(l + m - 1.0)
/ (np.sqrt(2.0 * l - 3.0) * np.sqrt(l + m) * np.sqrt(l - m))
)
k += 2
# u is sine of colatitude (cosine of latitude) so that 0 <= s <= 1
# for x=cos(th): u=sin(th)
u = np.sqrt(1.0 - x**2)
# update where u==0 to eps of data type to prevent invalid divisions
u[u == 0] = np.finfo(u.dtype).eps
# Calculate P(l,0). These are not scaled.
p[0, :] = 1.0
p[1, :] = np.sqrt(3.0) * x
k = 1
for l in range(2, LMAX + 1):
k += l
p[k, :] = f1[k] * x * p[k - l, :] - f2[k] * p[k - 2 * l + 1, :]
# Calculate P(m,m), P(m+1,m), and P(l,m)
pmm = np.sqrt(2.0) * scalef
rescalem = 1.0 / scalef
kstart = 0
for m in range(1, LMAX):
rescalem = rescalem * u
# Calculate P(m,m)
kstart += m + 1
pmm = pmm * np.sqrt(2 * m + 1) / np.sqrt(2 * m)
p[kstart, :] = pmm
# Calculate P(m+1,m)
k = kstart + m + 1
p[k, :] = x * np.sqrt(2 * m + 3) * pmm
# Calculate P(l,m)
for l in range(m + 2, LMAX + 1):
k += l
p[k, :] = x * f1[k] * p[k - l, :] - f2[k] * p[k - 2 * l + 1, :]
p[k - 2 * l + 1, :] = p[k - 2 * l + 1, :] * rescalem
# rescale
p[k, :] = p[k, :] * rescalem
p[k - LMAX, :] = p[k - LMAX, :] * rescalem
# Calculate P(LMAX,LMAX)
rescalem = rescalem * u
kstart += m + 2
p[kstart, :] = pmm * np.sqrt(2 * LMAX + 1) / np.sqrt(2 * LMAX) * rescalem
# reshape Legendre polynomials to output dimensions
for m in range(LMAX + 1):
for l in range(m, LMAX + 1):
lm = (l * (l + 1)) // 2 + m
plm[l, m, :] = p[lm, :]
# calculate first derivatives
if l == m:
dplm[l, m, :] = np.longdouble(m) * (x / u) * plm[l, m, :]
else:
flm = np.sqrt(
((l**2.0 - m**2.0) * (2.0 * l + 1.0)) / (2.0 * l - 1.0)
)
dplm[l, m, :] = (1.0 / u) * (
l * x * plm[l, m, :] - flm * plm[l - 1, m, :]
)
# return the legendre polynomials and their first derivative
# truncating orders to MMAX
return plm[:, : MMAX + 1, :], dplm[:, : MMAX + 1, :]
[docs]
def plm_mohlenkamp(LMAX, x, MMAX=None, astype=np.float64):
"""
Computes fully-normalized associated Legendre Polynomials and their
first derivative using the recursion relation from :cite:p:`Mohlenkamp:2016vv`
Derived from :cite:p:`Szego:1939tn` recurrence formula for Jacobi Polynomials
Parameters
----------
LMAX: int
maximum degree of Legendre polynomials
x: np.ndarray
elements ranging from -1 to 1
Typically ``cos(theta)``, where ``theta`` is the colatitude in radians
MMAX: int or NoneType, default None
maximum order of Associated Legendre polynomials
astype: np.dtype, default np.float64
output variable data type
Returns
-------
plm: np.ndarray
fully-normalized Legendre polynomials
dplm: np.ndarray
first derivative of Legendre polynomials
"""
# Verify LMAX as integer
LMAX = np.int64(LMAX)
# upper bound of spherical harmonic orders (default = LMAX)
if MMAX is None:
MMAX = np.copy(LMAX)
# removing singleton dimensions of x
x = np.atleast_1d(x).flatten()
# length of the x array
sx = len(x)
# Initialize the output Legendre polynomials
plm = np.zeros((LMAX + 1, MMAX + 1, sx), dtype=astype)
dplm = np.zeros((LMAX + 1, LMAX + 1, sx), dtype=astype)
# Jacobi polynomial for the recurrence relation
jlmm = np.zeros((LMAX + 1, MMAX + 1, sx))
# for x=cos(th): u= sin(th)
u = np.sqrt(1.0 - x**2)
# update where u==0 to eps of data type to prevent invalid divisions
u[u == 0] = np.finfo(u.dtype).eps
# for all spherical harmonic orders of interest
for mm in range(0, MMAX + 1): # equivalent to 0:MMAX
# Initialize the recurrence relation
# J-1,m,m Term == 0
# J0,m,m Term
if mm > 0:
# j ranges from 1 to mm for the product
j = np.arange(0, mm) + 1.0
jlmm[0, mm, :] = np.prod(np.sqrt(1.0 + 1.0 / (2.0 * j))) / np.sqrt(
2.0
)
else: # if mm == 0: jlmm = 1/sqrt(2)
jlmm[0, mm, :] = 1.0 / np.sqrt(2.0)
# Jk,m,m Terms
for k in range(1, LMAX + 1): # computation for SH degrees
# Initialization begins at -1
# this is to make the formula parallel the function written in
# Martin Mohlenkamp's Guide to Spherical Harmonics
# Jacobi General Terms
if k == 1: # for degree 1 terms
jlmm[k, mm, :] = (
2.0
* x
* jlmm[k - 1, mm, :]
* np.sqrt(1.0 + (mm - 0.5) / k)
* np.sqrt(1.0 - (mm - 0.5) / (k + 2.0 * mm))
)
else: # for all other spherical harmonic degrees
jlmm[k, mm, :] = 2.0 * x * jlmm[k - 1, mm, :] * np.sqrt(
1.0 + (mm - 0.5) / k
) * np.sqrt(1.0 - (mm - 0.5) / (k + 2.0 * mm)) - jlmm[
k - 2, mm, :
] * np.sqrt(1.0 + 4.0 / (2.0 * k + 2.0 * mm - 3.0)) * np.sqrt(
1.0 - (1.0 / k)
) * np.sqrt(1.0 - 1.0 / (k + 2.0 * mm))
# Normalization is geodesy convention
for l in range(mm, LMAX + 1): # equivalent to mm:LMAX
if mm == 0: # Geodesy normalization (m=0) == sqrt(2)*sin(th)^0
# u^mm term is dropped as u^0 = 1
plm[l, mm, :] = np.sqrt(2.0) * jlmm[l - mm, mm, :]
else: # Geodesy normalization all others == 2*sin(th)^mm
plm[l, mm, :] = 2.0 * (u**mm) * jlmm[l - mm, mm, :]
# calculate first derivatives
if l == mm:
dplm[l, mm, :] = np.longdouble(mm) * (x / u) * plm[l, mm, :]
else:
flm = np.sqrt(
((l**2.0 - mm**2.0) * (2.0 * l + 1.0)) / (2.0 * l - 1.0)
)
dplm[l, mm, :] = (1.0 / u) * (
l * x * plm[l, mm, :] - flm * plm[l - 1, mm, :]
)
# return the legendre polynomials and their first derivative
return plm, dplm