#!/usr/bin/env python
"""
gen_point_load.py
Written by Tyler Sutterley (07/2026)
Calculates gravitational spherical harmonic coefficients for point masses
CALLING SEQUENCE:
Ylms = gen_point_load(data, lon, lat, LMAX=LMAX)
INPUTS:
data: data magnitude
lon: longitude of points
lat: latitude of points
OUTPUTS:
clm: cosine spherical harmonic coefficients (geodesy normalization)
slm: sine spherical harmonic coefficients (geodesy normalization)
l: spherical harmonic degree to LMAX
m: spherical harmonic order to MMAX
OPTIONS:
LMAX: Upper bound of Spherical Harmonic Degrees
MMAX: Upper bound of Spherical Harmonic Orders
UNITS: input data units
1: grams of mass (default)
2: gigatonnes of mass
list: custom degree-dependent unit conversion factor
LOVE: input load Love numbers up to degree LMAX (hl,kl,ll)
PYTHON DEPENDENCIES:
numpy: Scientific Computing Tools For Python (https://numpy.org)
scipy: Scientific Tools for Python (https://docs.scipy.org/doc/)
PROGRAM DEPENDENCIES:
legendre.py: Computes associated Legendre polynomials for degree l
units.py: class for converting spherical harmonic data to specific units
harmonics.py: spherical harmonic data class for processing GRACE/GRACE-FO
destripe_harmonics.py: calculates the decorrelation (destriping) filter
and filters the GRACE/GRACE-FO coefficients for striping errors
REFERENCES:
I. M. Longman, Journal of Geophysical Research, 67(2), 1962
https://doi.org/10.1029/JZ067i002p00845
W. E. Farrell, Reviews of Geophysics and Space Physics, 10(3), 1972
https://doi.org/10.1029/RG010i003p00761
H. N. Pollack, Journal of Geophysical Research, 78(11), 1973
https://doi.org/10.1029/JB078i011p01760
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 from custom units
Updated 01/2021: use harmonics class for spherical harmonic operations
Updated 07/2020: added function docstrings
Written 05/2020
"""
import numpy as np
import gravity_toolkit.units
import gravity_toolkit.harmonics
from gravity_toolkit.legendre import legendre
[docs]
def gen_point_load(data, lon, lat, LMAX=60, MMAX=None, UNITS=1, LOVE=None):
"""
Calculates spherical harmonic coefficients for point masses
:cite:p:`Longman:1962ev,Farrell:1972cm,Pollack:1973gi`
Parameters
----------
data: np.ndarray
data magnitude
lon: np.ndarray
longitude of points
lat: np.ndarray
latitude of points
LMAX: int, default 60
Upper bound of Spherical Harmonic Degrees
MMAX: int or NoneType, default None
Upper bound of Spherical Harmonic Orders
UNITS: int, default 1
Input data units
- ``1``: grams of mass (g)
- ``2``: gigatonnes of mass (Gt)
- list: custom degree-dependent unit conversion factor
LOVE: tuple or NoneType, default None
Input load Love numbers up to degree LMAX (``hl``, ``kl``, ``ll``)
Returns
-------
clm: np.ndarray
cosine spherical harmonic coefficients
slm: np.ndarray
sine spherical harmonic coefficients
l: np.ndarray
spherical harmonic degree to LMAX
m: np.ndarray
spherical harmonic order to MMAX
"""
# upper bound of spherical harmonic orders (default == LMAX)
if MMAX is None:
MMAX = np.copy(LMAX)
# number of input data points
npts = len(data.flatten())
# convert output longitude and latitude into radians
phi = np.radians(lon.flatten())
theta = np.radians(90.0 - lat.flatten())
# extract degree dependent factor for specific units
factors = gravity_toolkit.units(lmax=LMAX)
int_fact = np.zeros((npts))
if isinstance(UNITS, (list, np.ndarray)):
# custom units
dfactor = np.copy(UNITS)
int_fact[:] = 1.0
elif UNITS == 1:
# Default Parameter: Input in grams (g)
dfactor = factors.spatial(*LOVE).cmwe / (factors.rad_e**2)
int_fact[:] = 1.0
elif UNITS == 2:
# Input in gigatonnes (Gt)
dfactor = factors.spatial(*LOVE).cmwe / (factors.rad_e**2)
int_fact[:] = 1e15
else:
raise ValueError(f'Unknown units {UNITS}')
# flattened form of data converted to units
D = int_fact * data.flatten()
# Initializing output spherical harmonic matrices
Ylms = gravity_toolkit.harmonics(lmax=LMAX, mmax=MMAX)
Ylms.clm = np.zeros((LMAX + 1, MMAX + 1))
Ylms.slm = np.zeros((LMAX + 1, MMAX + 1))
# for each degree l
for l in range(LMAX + 1):
m1 = np.min([l, MMAX]) + 1
SPH = _complex_harmonics(l, D, phi, theta, dfactor[l])
# truncate to spherical harmonic order and save to output
Ylms.clm[l, :m1] = SPH.real[:m1]
Ylms.slm[l, :m1] = SPH.imag[:m1]
# return the output spherical harmonics object
return Ylms
# calculate spherical harmonics of degree l evaluated at (theta,phi)
[docs]
def _complex_harmonics(l, data, phi, theta, coeff):
"""
Calculates the spherical harmonics for a particular degree
evaluated from data at coordinates
Parameters
----------
l: int
spherical harmonic degree
data: np.ndarray
data magnitude in grams
phi: np.ndarray
longitude of points in radians
theta: np.ndarray
colatitude of points in radians
coeff: np.ndarray
degree-dependent factor for converting units
Returns
-------
Ylms: np.ndarray
spherical harmonic coefficients in Eulerian form
"""
# calculate normalized legendre polynomials (order, points)
Pl = legendre(l, np.cos(theta), NORMALIZE=True)
# spherical harmonic orders up to degree l
m = np.arange(0, l + 1)
# calculate Euler's of order m multiplied by azimuth phi
m_phi = np.exp(1j * np.einsum('m...,p...->mp...', m, phi))
# reshape data to (order, points)
D = np.kron(np.ones((l + 1, 1)), data[np.newaxis, :])
# calculate spherical harmonics summing over all points
Yl = np.einsum('mp...,mp...,mp...->m...', D, Pl, m_phi)
# return harmonics for degree l multiplied by coefficients
return coeff * Yl