#!/usr/bin/env python
u"""
legendre_polynomials.py
Written by Tyler Sutterley (11/2024)
Computes fully normalized Legendre polynomials for an array of x values
and their first derivative
Calculates Legendre polynomials for zonal harmonics (order 0)
CALLING SEQUENCE:
Pl,dPl = legendre_polynomials(lmax, np.cos(theta))
INPUTS:
lmax: maximum 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 (geodesy normalization)
dpl: first derivative of Legendre polynomials
OPTIONS:
ASTYPE: output variable type. Default is np.float64
NOTES:
ptemp is a dummy array of length (0:lmax) storing unnormalized pl values
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python (https://numpy.org)
REFERENCES:
Hofmann-Wellenhof and Moritz, "Physical Geodesy" (2005)
http://www.springerlink.com/content/978-3-211-33544-4
UPDATE HISTORY:
Updated 11/2024: add polar argument for x == +/-1 to prevent drift
Updated 03/2023: improve typing for variables in docstrings
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 x variable
Updated 08/2020: prevent zero divisions by changing u==0 to eps of data type
Updated 07/2020: added function docstrings
Updated 07/2017: added first derivative of Legendre polynomials (dpl)
added option ASTYPE to output as different variable types e.g. np.float64
Written 03/2013
"""
import numpy as np
[docs]
def legendre_polynomials(lmax, x, ASTYPE=np.float64):
"""
Computes fully-normalized Legendre polynomials and their first derivative
following :cite:t:`HofmannWellenhof:2006hy`
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
ASTYPE: np.dtype, default np.float64
output variable data type
Returns
-------
pl: np.ndarray
fully-normalized Legendre polynomials
dpl: np.ndarray
first derivative of Legendre polynomials
"""
# verify dimensions
x = np.atleast_1d(x).flatten().astype(ASTYPE)
# size of the x array
nx = len(x)
# verify data type of spherical harmonic truncation
lmax = np.int64(lmax)
# output matrix of normalized legendre polynomials
pl = np.zeros((lmax+1,nx),dtype=ASTYPE)
# output matrix of First derivative of Legendre polynomials
dpl = np.zeros((lmax+1,nx),dtype=ASTYPE)
# dummy matrix for the recurrence relation
ptemp = np.zeros((lmax+1,nx),dtype=ASTYPE)
# 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
u0, = np.nonzero(u == 0)
u[u0] = np.finfo(u.dtype).eps
# Initialize the recurrence relation
ptemp[0,:] = 1.0
ptemp[1,:] = x
# Normalization is geodesy convention
pl[0,:] = ptemp[0,:]
pl[1,:] = np.sqrt(3.0)*ptemp[1,:]
for l in range(2,lmax+1):
ptemp[l,:] = (((2.0*l)-1.0)/l)*x*ptemp[l-1,:] - ((l-1.0)/l)*ptemp[l-2,:]
# Normalization is geodesy convention
pl[l,:] = np.sqrt((2.0*l)+1.0)*ptemp[l,:]
# Overwrite polar case (x == +/-1)
pl[l,u0] = np.sqrt((2.0*l)+1.0)*x[u0]**l
# First derivative of Legendre polynomials
for l in range(1,lmax+1):
fl = np.sqrt(((l**2.0) * (2.0*l + 1.0)) / (2.0*l - 1.0))
dpl[l,:] = (1.0/u)*(l*x*pl[l,:] - fl*pl[l-1,:])
# return the legendre polynomials and their first derivative
return (pl, dpl)