#!/usr/bin/env python
"""
clenshaw_summation.py
Written by Tyler Sutterley (07/2026)
Calculates the spatial field for a series of spherical harmonics for a
sequence of ungridded points
CALLING SEQUENCE:
spatial = clenshaw_summation(clm, slm, lon, lat, UNITS=1,
LMAX=60, LOVE=(hl,kl,ll))
INPUTS:
clm: cosine spherical harmonic coefficients
slm: sine spherical harmonic coefficients
lon: longitude of points
lat: latitude of points
OPTIONS:
RAD: Gaussian smoothing radius (km)
UNITS: output data units
1: cm of water thickness
2: mm of geoid height
3: mm of elastic crustal deformation [Davis 2004]
4: microGal gravitational perturbation
5: mbar equivalent surface pressure
6: cm of viscoelastic crustal uplift (GIA) [See Wahr 1995 or Wahr 2000]
list: custom degree-dependent unit conversion factor
LMAX: Upper bound of Spherical Harmonic Degrees
LOVE: input load Love numbers up to degree LMAX (hl,kl,ll)
ASTYPE: floating point precision for calculating Clenshaw summation
SCALE: scaling factor to prevent underflow in Clenshaw summation
OUTPUTS:
spatial: spatial field for lon/lat
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python (https://numpy.org)
PROGRAM DEPENDENCIES:
gauss_weights.py: Computes the Gaussian weights as a function of degree
units.py: class for converting spherical harmonic data to specific units
REFERENCES:
Holmes and Featherstone, "A Unified Approach to the Clenshaw Summation and
the Recursive Computation of Very High Degree and Order Normalised
Associated Legendre Functions", Journal of Geodesy (2002)
https://doi.org/10.1007/s00190-002-0216-2
Tscherning and Poder, "Some Geodetic Applications of Clenshaw Summation",
Bollettino di Geodesia e Scienze (1982)
UPDATE HISTORY:
Updated 07/2026: use np.einsum for spherical harmonic summations
use np.radians to convert from degrees to radians
Updated 04/2023: allow love numbers to be None for custom units case
Updated 03/2023: improve typing for variables in docstrings
Updated 02/2023: set custom units as top option in if/else statements
Updated 11/2022: use f-strings for formatting verbose or ascii output
Updated 04/2022: updated docstrings to numpy documentation format
Updated 11/2021: added UNITS list option for converting to custom units
Updated 09/2021: fix passing SCALE keyword argument to clenshaw_s_m
Updated 06/2021: output equivalent pressure in pascals
Updated 08/2020: parameterize float precision to improve computational time
Updated 07/2020: added function docstrings
Updated 04/2020: reading load love numbers outside of this function
using the units class for converting normalized spherical harmonics
Updated 03/2018: added option for output in equivalent pressure (UNITS=5)
simplified love number extrapolation if LMAX is greater than 696
Written 08/2017
"""
import numpy as np
from gravity_toolkit.gauss_weights import gauss_weights
from gravity_toolkit.units import units
[docs]
def clenshaw_summation(
clm,
slm,
lon,
lat,
RAD=0,
UNITS=0,
LMAX=0,
LOVE=None,
ASTYPE=np.longdouble,
SCALE=1e-280,
):
r"""
Calculates the spatial field for a series of spherical harmonics for a
sequence of ungridded points :cite:p:`Holmes:2002ff,Tscherning:1982tu`
Parameters
----------
clm: np.ndarray
cosine spherical harmonic coefficients
slm: np.ndarray
sine spherical harmonic coefficients
lon: np.ndarray
longitude of points
lat: np.ndarray
latitude of points
RAD: int or float, default 0
Gaussian smoothing radius (km)
UNITS: int, str, list or np.ndarray, default 0
Output data units
- ``1``: cm water equivalent thickness (cm w.e., g/cm\ :sup:`2`)
- ``2``: mm geoid height
- ``3``: mm elastic crustal deformation :cite:p:`Davis:2004il`
- ``4``: microGal gravitational perturbation
- ``5``: mbar equivalent surface pressure
- ``6``: cm viscoelastic crustal uplift (GIA) :cite:p:`Wahr:2000ek`
- list: custom degree-dependent unit conversion factor
LMAX: int, default 0
Upper bound of Spherical Harmonic Degrees
LOVE: tuple or NoneType, default None
Load Love numbers up to degree LMAX (``hl``, ``kl``, ``ll``)
ASTYPE: np.dtype, default np.longdouble
floating point precision for calculating Clenshaw summation
SCALE: float, default 1e-280
scaling factor to prevent underflow in Clenshaw summation
Returns
-------
spatial: np.ndarray
calculated spatial field for latitude and longitude
"""
# check if lat and lon are the same size
if len(lat) != len(lon):
raise ValueError('Incompatible vector dimensions (lon, lat)')
# calculate colatitude and longitude in radians
th = np.radians(90.0 - lat)
phi = np.squeeze(np.radians(lon))
# calculate cos and sin of colatitudes
t = np.cos(th)
u = np.sin(th)
# dimensions of theta and phi
npts = len(th)
# Gaussian Smoothing
if RAD != 0:
wl = 2.0 * np.pi * gauss_weights(RAD, LMAX)
else:
# else = 1
wl = np.ones((LMAX + 1))
# Setting units factor for output
# dfactor is the degree dependent coefficients
factors = units(lmax=LMAX)
if isinstance(UNITS, (list, np.ndarray)):
# custom units
dfactor = np.copy(UNITS)
elif isinstance(UNITS, str):
# named units
dfactor = factors.harmonic(*LOVE).get(UNITS)
elif isinstance(UNITS, int):
# use named unit codes
dfactor = factors.harmonic(*LOVE).get(units.bycode(UNITS))
else:
raise ValueError(f'Unknown units {UNITS}')
# complex spherical harmonics
ylm = clm - 1j * slm
# smooth degree dependent factors
f = dfactor * wl
# calculating cos(m*phi) and sin(m*phi) using Euler's formula
mm = np.arange(0, LMAX + 1)
m_phi = np.exp(1j * np.einsum('m...,p...->pm...', mm, phi))
# initiate summation
s_m = 0.0
# iterate to calculate complete summation
for m in range(LMAX, 0, -1):
# calculate summation for order m
a_m = np.sqrt((2.0 * m + 3.0) / (2.0 * m + 2.0))
cs_m = _clenshaw(t, f, m, ylm, LMAX, SCALE=SCALE)
# update summation and discard imaginary component
s_m = a_m * u * s_m + (cs_m * m_phi[:, m]).real
# add the final terms to calculate spatial field
cs_m = _clenshaw(t, f, 0, ylm, LMAX, SCALE=SCALE)
spatial = np.sqrt(3.0) * u * s_m + cs_m.real
# return the calculated spatial field
return spatial
# PURPOSE: compute Clenshaw summation of the fully normalized associated
# Legendre's function for constant order m
def _clenshaw(t, f, m, Ylm1, lmax, SCALE=1e-280):
"""
Compute conditioned arrays for Clenshaw summation from the fully-normalized
associated Legendre's function for an order m
Parameters
----------
t: np.ndarray
elements ranging from -1 to 1, typically cos(th)
f: np.ndarray
degree dependent factors
m: int
spherical harmonic order
Ylm1: np.ndarray
complex form of spherical harmonics
lmax: int
maximum spherical harmonic degree
SCALE: float, default 1e-280
scaling factor to prevent underflow in Clenshaw summation
Returns
-------
cs_m: np.ndarray
conditioned array for clenshaw summation
"""
# allocate for output matrix
N = len(t)
cs_m = np.zeros((N), dtype=np.clongdouble)
# scaling to prevent overflow
ylm = SCALE * Ylm1.astype(np.clongdouble)
# convert lmax and m to float
lm = np.float64(lmax)
mm = np.float64(m)
if m == lmax:
cs_m[:] = f[lmax] * ylm[lmax, lmax]
elif m == (lmax - 1):
a_lm = (
np.sqrt(
((2.0 * lm - 1.0) * (2.0 * lm + 1.0)) / ((lm - mm) * (lm + mm))
)
* t
)
cs_m[:] = (
a_lm * f[lmax] * ylm[lmax, lmax - 1]
+ f[lmax - 1] * ylm[lmax - 1, lmax - 1]
)
elif (m <= (lmax - 2)) and (m >= 1):
s_mm_minus_2 = f[lmax] * ylm[lmax, m]
a_lm = (
np.sqrt(
((2.0 * lm - 1.0) * (2.0 * lm + 1.0)) / ((lm - mm) * (lm + mm))
)
* t
)
s_mm_minus_1 = a_lm * s_mm_minus_2 + f[lmax - 1] * ylm[lmax - 1, m]
for l in range(lmax - 2, m - 1, -1):
ll = np.float64(l)
a_lm = (
np.sqrt(
((2.0 * ll + 1.0) * (2.0 * ll + 3.0))
/ ((ll + 1.0 - mm) * (ll + 1.0 + mm))
)
* t
)
b_lm = np.sqrt(
((2.0 * ll + 5.0) * (ll + mm + 1.0) * (ll - mm + 1.0))
/ ((ll + 2.0 - mm) * (ll + 2.0 + mm) * (2.0 * ll + 1.0))
)
s_mm_l = (
a_lm * s_mm_minus_1 - b_lm * s_mm_minus_2 + f[l] * ylm[l, m]
)
s_mm_minus_2 = np.copy(s_mm_minus_1)
s_mm_minus_1 = np.copy(s_mm_l)
cs_m[:] = np.copy(s_mm_l)
elif m == 0:
s_mm_minus_2 = f[lmax] * ylm[lmax, 0]
a_lm = np.sqrt(((2.0 * lm - 1.0) * (2.0 * lm + 1.0)) / (lm * lm)) * t
s_mm_minus_1 = a_lm * s_mm_minus_2 + f[lmax - 1] * ylm[lmax - 1, 0]
for l in range(lmax - 2, m - 1, -1):
ll = np.float64(l)
a_lm = (
np.sqrt(
((2.0 * ll + 1.0) * (2.0 * ll + 3.0))
/ ((ll + 1.0) * (ll + 1.0))
)
* t
)
b_lm = np.sqrt(
((2.0 * ll + 5.0) * (ll + 1.0) * (ll + 1.0))
/ ((ll + 2.0) * (ll + 2.0) * (2.0 * ll + 1.0))
)
s_mm_l = (
a_lm * s_mm_minus_1 - b_lm * s_mm_minus_2 + f[l] * ylm[l, 0]
)
s_mm_minus_2 = np.copy(s_mm_minus_1)
s_mm_minus_1 = np.copy(s_mm_l)
cs_m[:] = np.copy(s_mm_l)
# return rescaled cs_m
return cs_m / SCALE