#!/usr/bin/env python
"""
legendre.py
Written by Tyler Sutterley (03/2023)
Computes associated Legendre functions of degree l evaluated for elements x
l must be a scalar integer and x must contain real values ranging -1 <= x <= 1
Parallels the MATLAB legendre function
Based on Fortran program by Robert L. Parker, Scripps Institution of
Oceanography, Institute for Geophysics and Planetary Physics, UCSD. 1993
INPUTS:
l: degree of Legendre polynomials
x: elements ranging from -1 to 1
typically cos(theta), where theta is the colatitude in radians
OUTPUT:
Pl: legendre polynomials of degree l for orders 0 to l
OPTIONS:
NORMALIZE: output Fully Normalized Associated Legendre Functions
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python (https://numpy.org)
REFERENCES:
M. Abramowitz and I.A. Stegun, "Handbook of Mathematical Functions",
Dover Publications, 1965, Ch. 8.
J. A. Jacobs, "Geomagnetism", Academic Press, 1987, Ch.4.
UPDATE HISTORY:
Updated 03/2023: improve typing for variables in docstrings
Updated 04/2022: updated docstrings to numpy documentation format
Updated 11/2021: modify normalization to prevent high degree overflows
Updated 05/2021: define int/float precision to prevent deprecation warning
Updated 02/2021: modify case with underflow
Updated 09/2020: verify dimensions of x variable
Updated 07/2020: added function docstrings
Updated 05/2020: added normalization option for output polynomials
Updated 03/2019: calculate twocot separately to avoid divide warning
Written 08/2016
"""
import numpy as np
[docs]
def legendre(l, x, NORMALIZE=False):
"""
Computes associated Legendre functions for a particular degree
following :cite:t:`Abramowitz:1965vw,Jacobs:1987vv`
Parameters
----------
l: int
degree of Legendre polynomials
x: np.ndarray
elements ranging from -1 to 1
Typically ``cos(theta)``, where ``theta`` is the colatitude in radians
NORMALIZE: bool, default False
Fully-normalize the Legendre Functions
Returns
-------
Pl: np.ndarray
legendre polynomials of degree ``l``
"""
# verify integer
l = np.int64(l)
# verify dimensions
x = np.atleast_1d(x).flatten()
# size of the x array
nx = len(x)
# for the l = 0 case
if l == 0:
Pl = np.ones((1, nx), dtype=np.float64)
return Pl
# for all other degrees greater than 0
rootl = np.sqrt(np.arange(0, 2 * l + 1)) # +1 to include 2*l
# s is sine of colatitude (cosine of latitude) so that 0 <= s <= 1
s = np.sqrt(1.0 - x**2) # for x=cos(th): s=sin(th)
P = np.zeros((l + 3, nx), dtype=np.float64)
# Find values of x,s for which there will be underflow
sn = (-s) ** l
tol = np.sqrt(np.finfo(np.float64).tiny)
count = np.count_nonzero((s > 0) & (np.abs(sn) <= tol))
if count > 0:
(ind,) = np.nonzero((s > 0) & (np.abs(sn) <= tol))
# Approximate solution of x*ln(x) = Pl
v = 9.2 - np.log(tol) / (l * s[ind])
w = 1.0 / np.log(v)
m1 = 1 + l * s[ind] * v * w * (1.0058 + w * (3.819 - w * 12.173))
m1 = np.where(l < np.floor(m1), l, np.floor(m1)).astype(np.int64)
# Column-by-column recursion
for k, mm1 in enumerate(m1):
col = ind[k]
# Calculate two*cotangent for underflow case
twocot = -2.0 * x[col] / s[col]
P[mm1 - 1 : l + 1, col] = 0.0
# Start recursion with proper sign
tstart = np.finfo(np.float64).eps
P[mm1 - 1, col] = np.sign(np.fmod(mm1, 2) - 0.5) * tstart
if x[col] < 0:
P[mm1 - 1, col] = np.sign(np.fmod(l + 1, 2) - 0.5) * tstart
# Recur from m1 to m = 0, accumulating normalizing factor.
sumsq = tol.copy()
for m in range(mm1 - 2, -1, -1):
P[m, col] = (
(m + 1) * twocot * P[m + 1, col]
- rootl[l + m + 2] * rootl[l - m - 1] * P[m + 2, col]
) / (rootl[l + m + 1] * rootl[l - m])
sumsq += P[m, col] ** 2
# calculate scale
scale = 1.0 / np.sqrt(2.0 * sumsq - P[0, col] ** 2)
P[0 : mm1 + 1, col] = scale * P[0 : mm1 + 1, col]
# Find the values of x,s for which there is no underflow, and (x != +/-1)
count = np.count_nonzero((x != 1) & (np.abs(sn) >= tol))
if count > 0:
(nind,) = np.nonzero((x != 1) & (np.abs(sn) >= tol))
# Calculate two*cotangent for normal case
twocot = -2.0 * x[nind] / s[nind]
# Produce normalization constant for the m = l function
d = np.arange(2, 2 * l + 2, 2)
c = np.prod(1.0 - 1.0 / d)
# Use sn = (-s)**l (written above) to write the m = l function
P[l, nind] = np.sqrt(c) * sn[nind]
P[l - 1, nind] = P[l, nind] * twocot * l / rootl[-1]
# Recur downwards to m = 0
for m in range(l - 2, -1, -1):
P[m, nind] = (
P[m + 1, nind] * twocot * (m + 1)
- P[m + 2, nind] * rootl[l + m + 2] * rootl[l - m - 1]
) / (rootl[l + m + 1] * rootl[l - m])
# calculate Pl from P
Pl = np.copy(P[0 : l + 1, :])
# Polar argument (x == +/-1)
count = np.count_nonzero(s == 0)
if count > 0:
(s0,) = np.nonzero(s == 0)
Pl[0, s0] = x[s0] ** l
# calculate Fully Normalized Associated Legendre functions
if NORMALIZE:
norm = np.zeros((l + 1))
norm[0] = np.sqrt(2.0 * l + 1)
m = np.arange(1, l + 1)
norm[1:] = (-1) ** m * np.sqrt(2.0 * (2.0 * l + 1.0))
Pl *= np.kron(np.ones((1, nx)), norm[:, np.newaxis])
else:
# Calculate the unnormalized Legendre functions by multiplying each row
# by: sqrt((l+m)!/(l-m)!) == sqrt(prod(n-m+1:n+m))
# following Abramowitz and Stegun
for m in range(1, l):
Pl[m, :] *= np.prod(rootl[l - m + 1 : l + m + 1])
# sectoral case (l = m) should be done separately to handle 0!
Pl[l, :] *= np.prod(rootl[1:])
return Pl