#!/usr/bin/env python
u"""
legendre.py
Written by Tyler Sutterley (05/2022)
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 Legrendre 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)
scipy: Scientific Tools for Python (https://docs.scipy.org/doc/)
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 05/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 of degree ``l``
following [Abramowitz1965]_ and [Jacobs1987]_
Parameters
----------
l: int
degree of Legrendre polynomials
x: float
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: legendre polynomials of degree ``l``
References
----------
.. [Abramowitz1965] M. Abramowitz and I. A. Stegun,
*Handbook of Mathematical Functions*, 1046 pp., (1965).
.. [Jacobs1987] J. A. Jacobs, *Geomagnetism*,
Volume 1, 1st Edition, 832 pp., (1987).
"""
#-- 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 twocot 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 twocot 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