Gravity Models

Earth gravity field models and computations.

Gravity models module.

This module provides functions for computing gravitational acceleration and potential using various models including WGS84 normal gravity and spherical harmonic expansions.

Examples

>>> from pytcl.gravity import normal_gravity, gravity_wgs84
>>> import numpy as np

>>> # Normal gravity at 45° latitude, sea level
>>> g = normal_gravity(np.radians(45), 0)
>>> print(f"Gravity: {g:.4f} m/s^2")

>>> # Full gravity vector
>>> result = gravity_wgs84(np.radians(45), 0, 1000)
>>> print(f"Gravity magnitude: {result.magnitude:.4f} m/s^2")

pytcl.gravity.associated_legendre(n_max, m_max, x, normalized=True)

Compute associated Legendre polynomials P_n^m(x).

Uses the recursive algorithm for numerical stability.

Parameters:

• n_max (int) – Maximum degree.

• m_max (int) – Maximum order (must be <= n_max).

• x (float) – Argument, typically cos(colatitude). Must be in [-1, 1].

• normalized (bool, optional) – If True, return fully normalized (geodetic) coefficients. Default True.

Returns:

P – Array of shape (n_max+1, m_max+1) containing P_n^m(x).

Return type:

ndarray

Notes

The fully normalized associated Legendre functions satisfy:

\[\int_{-1}^{1} [\bar{P}_n^m(x)]^2 dx = \frac{2}{2n+1}\]

Results are cached for repeated queries with the same parameters. Cache key quantizes x to 8 decimal places (~1e-8 precision).

Examples

>>> P = associated_legendre(2, 2, 0.5)
>>> P[2, 0]  # P_2^0(0.5)

pytcl.gravity.associated_legendre_derivative(n_max, m_max, x, P=None, normalized=True)

Compute derivatives of associated Legendre polynomials dP_n^m/dx.

Parameters:

• n_max (int) – Maximum degree.

• m_max (int) – Maximum order.

• x (float) – Argument in [-1, 1].

• P (ndarray, optional) – Precomputed P_n^m values. If None, computed internally.

• normalized (bool, optional) – If True, use fully normalized functions. Default True.

Returns:

dP – Array of shape (n_max+1, m_max+1) containing dP_n^m/dx.

Return type:

ndarray

Examples

>>> import numpy as np
>>> x = np.cos(np.radians(45))  # cos of 45 degrees
>>> dP = associated_legendre_derivative(2, 2, x)
>>> dP.shape
(3, 3)
>>> abs(dP[0, 0]) < 1e-10  # dP_0^0/dx = 0
True

pytcl.gravity.spherical_harmonic_sum(lat, lon, r, C, S, R, GM, n_max=None)

Evaluate spherical harmonic expansion for a scalar field.

Computes the value and gradient of a field represented by spherical harmonic coefficients C_nm and S_nm.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• r (float) – Radial distance from center of mass.

• C (ndarray) – Cosine coefficients C_nm, shape (n_max+1, n_max+1).

• S (ndarray) – Sine coefficients S_nm, shape (n_max+1, n_max+1).

• R (float) – Reference radius (e.g., Earth’s equatorial radius).

• GM (float) – Gravitational parameter (G * M).

• n_max (int, optional) – Maximum degree to use. Default uses full coefficient array.

Returns:

• V (float) – Potential value.

• dV_r (float) – Radial derivative dV/dr.

• dV_lat (float) – Latitudinal derivative (1/r) * dV/dlat.

Return type:

Tuple[float, float, float]

Notes

The spherical harmonic expansion of the gravitational potential is:

\[V = \frac{GM}{r} \sum_{n=0}^{N} \left(\frac{R}{r}\right)^n \sum_{m=0}^{n} \bar{P}_n^m(\sin\phi) (C_{nm}\cos m\lambda + S_{nm}\sin m\lambda)\]

Examples

>>> import numpy as np
>>> # Simple monopole (degree 0 only)
>>> C = np.array([[1.0]])
>>> S = np.array([[0.0]])
>>> R = 6.378e6  # meters
>>> GM = 3.986e14  # m^3/s^2
>>> V, dV_r, dV_lat = spherical_harmonic_sum(0, 0, R, C, S, R, GM, n_max=0)
>>> abs(V - GM/R) / (GM/R) < 1e-10  # V = GM/r for degree 0
True

pytcl.gravity.gravity_acceleration(lat, lon, h, C, S, R, GM, n_max=None)

Compute gravity acceleration vector from spherical harmonics.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float) – Height above reference ellipsoid in meters.

• C (ndarray) – Cosine coefficients.

• S (ndarray) – Sine coefficients.

• R (float) – Reference radius.

• GM (float) – Gravitational parameter.

• n_max (int, optional) – Maximum degree.

Returns:

• g_r (float) – Radial component of gravity (positive outward).

• g_lat (float) – Northward component of gravity.

• g_lon (float) – Eastward component of gravity.

Return type:

Tuple[float, float, float]

Examples

>>> import numpy as np
>>> # Simple monopole field
>>> C = np.array([[1.0]])
>>> S = np.array([[0.0]])
>>> R = 6.378e6
>>> GM = 3.986e14
>>> g_r, g_lat, g_lon = gravity_acceleration(0, 0, 0, C, S, R, GM, n_max=0)
>>> g_r < 0  # Gravity points inward (negative radial)
True

pytcl.gravity.legendre_scaling_factors(n_max)

Precompute scaling factors to prevent overflow in Legendre recursion.

For degrees > ~150, standard Legendre recursion can overflow in double precision. These scaling factors keep intermediate values in a representable range.

The scaling follows the approach of Holmes & Featherstone (2002), using a factor that grows with degree to counteract the natural growth of the Legendre functions.

Parameters:

n_max (int) – Maximum degree.

Returns:

scale – Scaling factors of shape (n_max+1,). The factor for degree n is 10^(-280 * n / n_max) for n_max > 150, else 1.0.

Return type:

ndarray

References

[1]

Holmes, S.A. and Featherstone, W.E. “A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions.” Journal of Geodesy 76.5 (2002): 279-299.

Examples

>>> scale = legendre_scaling_factors(100)
>>> len(scale)
101
>>> scale[0]  # No scaling for low degrees
1.0

>>> scale_high = legendre_scaling_factors(200)
>>> scale_high[200] < scale_high[0]  # Higher degrees scaled down
True

pytcl.gravity.associated_legendre_scaled(n_max, m_max, x, scale=None)

Compute scaled associated Legendre polynomials for high degrees.

For ultra-high degree computations (n > 150), standard Legendre recursion overflows. This function computes scaled values that stay in representable range.

Parameters:

• n_max (int) – Maximum degree.

• m_max (int) – Maximum order (must be <= n_max).

• x (float) – Argument in [-1, 1], typically cos(colatitude).

• scale (ndarray, optional) – Precomputed scaling factors from legendre_scaling_factors(). If None, computed internally.

Returns:

• P_scaled (ndarray) – Scaled Legendre values, shape (n_max+1, m_max+1). The actual value is P_scaled[n,m] * 10^scale_exp[n].

• scale_exp (ndarray) – Scale exponents for each degree, shape (n_max+1,). Set to 0 if no scaling needed.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

The returned values satisfy:

P_n^m(x) = P_scaled[n,m] * 10^scale_exp[n]

For normal operations (n < 150), scale_exp is all zeros and P_scaled equals the actual Legendre values.

Examples

>>> import numpy as np
>>> x = np.cos(np.radians(45))
>>> P_scaled, scale_exp = associated_legendre_scaled(10, 10, x)
>>> P_scaled.shape
(11, 11)
>>> all(scale_exp == 0)  # No scaling needed for n_max < 150
True

class pytcl.gravity.GravityConstants(GM, a, f, omega, J2)

Bases: NamedTuple

Constants for a gravity model.

GM

Gravitational parameter (m^3/s^2).

Type:

float

a

Semi-major axis / equatorial radius (m).

Type:

float

f

Flattening.

Type:

float

omega

Angular velocity (rad/s).

Type:

float

J2

Second degree zonal harmonic (unnormalized).

Type:

float

GM: float

Alias for field number 0

a: float

Alias for field number 1

f: float

Alias for field number 2

omega: float

Alias for field number 3

J2: float

Alias for field number 4

class pytcl.gravity.GravityResult(magnitude, g_down, g_north, g_east)

Bases: NamedTuple

Result of gravity computation.

magnitude

Total gravity magnitude (m/s^2).

Type:

float

g_down

Downward component (positive down) (m/s^2).

Type:

float

g_north

Northward component (m/s^2).

Type:

float

g_east

Eastward component (m/s^2).

Type:

float

magnitude: float

Alias for field number 0

g_down: float

Alias for field number 1

g_north: float

Alias for field number 2

g_east: float

Alias for field number 3

pytcl.gravity.normal_gravity_somigliana(lat, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute normal gravity on the ellipsoid using Somigliana’s formula.

Parameters:

• lat (float) – Geodetic latitude in radians.

• constants (GravityConstants, optional) – Gravity model constants. Default WGS84.

Returns:

gamma – Normal gravity on the ellipsoid surface (m/s^2).

Return type:

float

Notes

Somigliana’s closed formula gives the exact normal gravity on the reference ellipsoid without iteration.

Examples

>>> gamma = normal_gravity_somigliana(0)  # At equator
>>> abs(gamma - 9.7803) < 0.001
True

>>> gamma = normal_gravity_somigliana(np.pi/2)  # At pole
>>> abs(gamma - 9.8322) < 0.001
True

pytcl.gravity.normal_gravity(lat, h=0.0, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute normal gravity at a given latitude and height.

Uses the International Gravity Formula with free-air correction.

Parameters:

• lat (float) – Geodetic latitude in radians.

• h (float, optional) – Height above ellipsoid in meters. Default 0.

• constants (GravityConstants, optional) – Gravity model constants. Default WGS84.

Returns:

g – Normal gravity (m/s^2).

Return type:

float

Examples

>>> g = normal_gravity(0, 0)  # Sea level at equator
>>> abs(g - 9.78) < 0.01
True

>>> g = normal_gravity(np.radians(45), 1000)  # 1km altitude, 45° lat
>>> g < normal_gravity(np.radians(45), 0)  # Decreases with altitude
True

pytcl.gravity.gravity_wgs84(lat, lon, h=0.0)

Compute gravity using WGS84 model.

This computes the full gravity vector including the centrifugal acceleration due to Earth’s rotation.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float, optional) – Height above WGS84 ellipsoid in meters. Default 0.

Returns:

result – Gravity components and magnitude.

Return type:

GravityResult

Examples

>>> result = gravity_wgs84(0, 0, 0)
>>> abs(result.magnitude - 9.78) < 0.01
True

pytcl.gravity.gravity_j2(lat, lon, h=0.0, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute gravity using J2 (oblateness) model.

This simplified model includes only the J2 zonal harmonic, which accounts for Earth’s equatorial bulge.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float, optional) – Height above ellipsoid in meters. Default 0.

• constants (GravityConstants, optional) – Model constants. Default WGS84.

Returns:

result – Gravity components and magnitude.

Return type:

GravityResult

Examples

>>> result = gravity_j2(0, 0, 0)  # At equator, sea level
>>> abs(result.magnitude - 9.78) < 0.01
True

pytcl.gravity.geoid_height_j2(lat, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute approximate geoid height using J2 model.

Parameters:

• lat (float) – Geodetic latitude in radians.

• constants (GravityConstants, optional) – Model constants. Default WGS84.

Returns:

N – Geoid height (geoid - ellipsoid) in meters.

Return type:

float

Examples

>>> N = geoid_height_j2(0)  # At equator
>>> N > 0  # Equator bulges outward
True

Notes

This is a simplified model. For accurate geoid heights, use a full geoid model like EGM2008.

pytcl.gravity.gravitational_potential(lat, lon, r, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263), n_max=2)

Compute gravitational potential.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• r (float) – Radial distance from Earth’s center in meters.

• constants (GravityConstants, optional) – Model constants. Default WGS84.

• n_max (int, optional) – Maximum spherical harmonic degree. Default 2.

Returns:

U – Gravitational potential (m^2/s^2).

Return type:

float

Examples

>>> U = gravitational_potential(0, 0, 6.4e6)  # Near Earth surface
>>> U < 0  # Potential is negative
True

pytcl.gravity.free_air_anomaly(g_observed, lat, h, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute free-air gravity anomaly.

Parameters:

• g_observed (float) – Observed gravity in m/s^2.

• lat (float) – Geodetic latitude in radians.

• h (float) – Height above geoid in meters.

• constants (GravityConstants, optional) – Model constants. Default WGS84.

Returns:

delta_g – Free-air anomaly in m/s^2 (or mGal if multiplied by 1e5).

Return type:

float

Examples

>>> import numpy as np
>>> # Observed gravity slightly higher than normal
>>> delta_g = free_air_anomaly(9.81, np.radians(45), 100)
>>> isinstance(delta_g, float)
True

Notes

The free-air anomaly is the difference between observed gravity and normal gravity at the observation point.

pytcl.gravity.bouguer_anomaly(g_observed, lat, h, rho=2670.0, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute simple Bouguer gravity anomaly.

Parameters:

• g_observed (float) – Observed gravity in m/s^2.

• lat (float) – Geodetic latitude in radians.

• h (float) – Height above geoid in meters.

• rho (float, optional) – Crustal density in kg/m^3. Default 2670 (average crustal density).

• constants (GravityConstants, optional) – Model constants. Default WGS84.

Returns:

delta_g – Bouguer anomaly in m/s^2.

Return type:

float

Examples

>>> import numpy as np
>>> # Bouguer anomaly at mountain location
>>> delta_g = bouguer_anomaly(9.81, np.radians(45), 1000)
>>> isinstance(delta_g, float)
True

Notes

The Bouguer anomaly removes the gravitational effect of the topographic mass between the observation point and the geoid.

pytcl.gravity.clenshaw_sum_order(m, cos_theta, sin_theta, C, S, n_max)

Clenshaw summation for fixed order m, summing over degrees n=m to n_max.

Evaluates the partial sums:

sum_C = sum_{n=m}^{n_max} C[n,m] * P_n^m(cos_theta) sum_S = sum_{n=m}^{n_max} S[n,m] * P_n^m(cos_theta)

using backward recursion from n_max down to m.

Parameters:

• m (int) – Order (fixed for this summation).

• cos_theta (float) – Cosine of colatitude.

• sin_theta (float) – Sine of colatitude.

• C (ndarray) – Cosine coefficients array, shape (n_max+1, n_max+1).

• S (ndarray) – Sine coefficients array, shape (n_max+1, n_max+1).

• n_max (int) – Maximum degree.

Returns:

• sum_C (float) – Sum of C terms weighted by Legendre functions.

• sum_S (float) – Sum of S terms weighted by Legendre functions.

Return type:

Tuple[float, float]

Examples

>>> import numpy as np
>>> C = np.zeros((5, 5))
>>> S = np.zeros((5, 5))
>>> C[2, 0] = 1.0  # Only C20 term
>>> cos_theta, sin_theta = np.cos(np.pi/4), np.sin(np.pi/4)
>>> sum_C, sum_S = clenshaw_sum_order(0, cos_theta, sin_theta, C, S, 4)
>>> isinstance(sum_C, float)
True

pytcl.gravity.clenshaw_sum_order_derivative(m, cos_theta, sin_theta, C, S, n_max)

Clenshaw summation with derivative for fixed order m.

Evaluates both the partial sums and their derivatives with respect to colatitude.

Parameters:

• m (int) – Order.

• cos_theta (float) – Cosine of colatitude.

• sin_theta (float) – Sine of colatitude.

• C (ndarray) – Cosine coefficients.

• S (ndarray) – Sine coefficients.

• n_max (int) – Maximum degree.

Returns:

• sum_C (float) – Sum of C terms.

• sum_S (float) – Sum of S terms.

• dsum_C (float) – Derivative of sum_C with respect to theta.

• dsum_S (float) – Derivative of sum_S with respect to theta.

Return type:

Tuple[float, float, float, float]

Examples

>>> import numpy as np
>>> C = np.zeros((5, 5))
>>> S = np.zeros((5, 5))
>>> C[2, 0] = -0.0005  # J2-like term
>>> cos_theta, sin_theta = np.cos(np.pi/4), np.sin(np.pi/4)
>>> sum_C, sum_S, dsum_C, dsum_S = clenshaw_sum_order_derivative(
...     0, cos_theta, sin_theta, C, S, 4)
>>> len([sum_C, sum_S, dsum_C, dsum_S])
4

pytcl.gravity.clenshaw_potential(lat, lon, r, C, S, R, GM, n_max=None)

Compute gravitational potential using Clenshaw summation.

Evaluates the spherical harmonic expansion of the gravitational potential efficiently using Clenshaw’s algorithm.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• r (float) – Radial distance from Earth center in meters.

• C (ndarray) – Cosine coefficients (fully normalized).

• S (ndarray) – Sine coefficients (fully normalized).

• R (float) – Reference radius in meters.

• GM (float) – Gravitational parameter in m^3/s^2.

• n_max (int, optional) – Maximum degree.

Returns:

Gravitational potential in m^2/s^2.

Return type:

float

Examples

>>> import numpy as np
>>> C = np.zeros((5, 5))
>>> S = np.zeros((5, 5))
>>> C[0, 0] = 1.0  # Central term only
>>> R = 6.378e6
>>> GM = 3.986e14
>>> V = clenshaw_potential(0, 0, R, C, S, R, GM)
>>> abs(V - GM/R) / (GM/R) < 0.01  # ~GM/r for central term
True

pytcl.gravity.clenshaw_gravity(lat, lon, r, C, S, R, GM, n_max=None)

Compute gravity disturbance vector using Clenshaw summation.

Evaluates both the potential and its gradient efficiently using Clenshaw’s algorithm with derivative recursions.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• r (float) – Radial distance from Earth center in meters.

• C (ndarray) – Cosine coefficients (fully normalized).

• S (ndarray) – Sine coefficients (fully normalized).

• R (float) – Reference radius in meters.

• GM (float) – Gravitational parameter in m^3/s^2.

• n_max (int, optional) – Maximum degree.

Returns:

• g_r (float) – Radial component of gravity disturbance (positive outward) in m/s^2.

• g_lat (float) – Northward component of gravity disturbance in m/s^2.

• g_lon (float) – Eastward component of gravity disturbance in m/s^2.

Return type:

Tuple[float, float, float]

Examples

>>> import numpy as np
>>> C = np.zeros((5, 5))
>>> S = np.zeros((5, 5))
>>> C[0, 0] = 1.0
>>> R = 6.378e6
>>> GM = 3.986e14
>>> g_r, g_lat, g_lon = clenshaw_gravity(0, 0, R, C, S, R, GM)
>>> g_r < 0  # Gravity points inward
True

class pytcl.gravity.EGMCoefficients(C, S, GM, R, n_max, model_name)

Bases: NamedTuple

Earth Gravitational Model coefficients.

Parameters:

• C (ndarray) – Cosine coefficients (fully normalized), shape (n_max+1, n_max+1).

• S (ndarray) – Sine coefficients (fully normalized), shape (n_max+1, n_max+1).

• GM (float) – Gravitational parameter in m^3/s^2.

• R (float) – Reference radius in meters.

• n_max (int) – Maximum degree of the model.

• model_name (str) – Model name (“EGM96” or “EGM2008”).

C: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

S: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

GM: float

Alias for field number 2

R: float

Alias for field number 3

n_max: int

Alias for field number 4

model_name: str

Alias for field number 5

class pytcl.gravity.GeoidResult(height, lat, lon, model)

Bases: NamedTuple

Geoid height computation result.

Parameters:

• height (float) – Geoid height above reference ellipsoid in meters.

• lat (float) – Latitude in radians.

• lon (float) – Longitude in radians.

• model (str) – Model name used.

height: float

Alias for field number 0

lat: float

Alias for field number 1

lon: float

Alias for field number 2

model: str

Alias for field number 3

class pytcl.gravity.GravityDisturbance(delta_g_r, delta_g_lat, delta_g_lon, magnitude)

Bases: NamedTuple

Gravity disturbance result.

Parameters:

• delta_g_r (float) – Radial component of gravity disturbance in m/s^2.

• delta_g_lat (float) – Northward component in m/s^2.

• delta_g_lon (float) – Eastward component in m/s^2.

• magnitude (float) – Total magnitude of disturbance in m/s^2.

delta_g_r: float

Alias for field number 0

delta_g_lat: float

Alias for field number 1

delta_g_lon: float

Alias for field number 2

magnitude: float

Alias for field number 3

pytcl.gravity.get_data_dir()

Get the pytcl data directory for external data files.

The data directory is located at ~/.pytcl/data/ by default. Can be overridden by setting the PYTCL_DATA_DIR environment variable.

Returns:

Path to the data directory.

Return type:

Path

pytcl.gravity.load_egm_coefficients(model='EGM2008', n_max=None)

Load EGM coefficients from file.

Parameters:

• model (str) – Model name, either “EGM96” or “EGM2008”.

• n_max (int, optional) – Maximum degree to load. If None, loads the full model.

Returns:

Loaded coefficient structure.

Return type:

EGMCoefficients

Raises:

• FileNotFoundError – If the coefficient file is not found.

• ValueError – If an unknown model is specified.

Examples

>>> coef = load_egm_coefficients("EGM2008", n_max=360)
>>> coef.n_max
360

pytcl.gravity.geoid_height(lat, lon, model='EGM2008', n_max=None, coefficients=None)

Compute geoid height above WGS84 ellipsoid.

The geoid is the equipotential surface of the Earth’s gravity field that best fits mean sea level. This function computes the height of the geoid above the WGS84 reference ellipsoid.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• model (str) – Model name (“EGM96” or “EGM2008”).

• n_max (int, optional) – Maximum degree to use. Default uses full model.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients. If None, loads from file.

Returns:

Geoid height in meters.

Return type:

float

Examples

>>> # Geoid height at equator, prime meridian
>>> N = geoid_height(0, 0)  # Should be approximately 17 m

pytcl.gravity.geoid_heights(lats, lons, model='EGM2008', n_max=None, coefficients=None)

Compute geoid heights for multiple points.

Parameters:

• lats (ndarray) – Geodetic latitudes in radians.

• lons (ndarray) – Longitudes in radians.

• model (str) – Model name.

• n_max (int, optional) – Maximum degree.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients. If None, loads from file.

Returns:

Geoid heights in meters.

Return type:

ndarray

Examples

>>> import numpy as np
>>> coef = create_test_coefficients(n_max=10)
>>> lats = np.array([0.0, np.pi/4])  # Equator and 45°N
>>> lons = np.array([0.0, np.pi/2])  # Prime meridian and 90°E
>>> heights = geoid_heights(lats, lons, coefficients=coef)
>>> len(heights)
2

pytcl.gravity.gravity_disturbance(lat, lon, h=0.0, model='EGM2008', n_max=None, coefficients=None)

Compute gravity disturbance from EGM model.

The gravity disturbance is the difference between actual gravity and normal (reference ellipsoid) gravity at the same point.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float) – Height above ellipsoid in meters.

• model (str) – Model name.

• n_max (int, optional) – Maximum degree.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients.

Returns:

Gravity disturbance components.

Return type:

GravityDisturbance

Examples

>>> coef = create_test_coefficients(n_max=10)
>>> dist = gravity_disturbance(0, 0, h=0, coefficients=coef)
>>> isinstance(dist.magnitude, float)
True

pytcl.gravity.gravity_anomaly(lat, lon, h=0.0, model='EGM2008', n_max=None, coefficients=None)

Compute free-air gravity anomaly.

The gravity anomaly is the difference between observed gravity and normal gravity, with the normal gravity evaluated at the geoid rather than at the observation point.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float) – Height above ellipsoid in meters.

• model (str) – Model name.

• n_max (int, optional) – Maximum degree.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients. If None, loads from file.

Returns:

Gravity anomaly in m/s^2 (typically reported in mGal = 1e-5 m/s^2).

Return type:

float

Examples

>>> coef = create_test_coefficients(n_max=10)
>>> anomaly = gravity_anomaly(0, 0, h=0, coefficients=coef)
>>> isinstance(anomaly, float)
True

pytcl.gravity.deflection_of_vertical(lat, lon, model='EGM2008', n_max=None, coefficients=None)

Compute deflection of the vertical.

The deflection of the vertical is the angle between the direction of gravity (plumb line) and the normal to the reference ellipsoid.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• model (str) – Model name.

• n_max (int, optional) – Maximum degree.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients. If None, loads from file.

Returns:

• xi (float) – North-south deflection component in arcseconds.

• eta (float) – East-west deflection component in arcseconds.

Return type:

Tuple[float, float]

Notes

Positive xi means the plumb line points more north than the normal. Positive eta means the plumb line points more east than the normal.

Examples

>>> coef = create_test_coefficients(n_max=10)
>>> xi, eta = deflection_of_vertical(0, 0, coefficients=coef)
>>> isinstance(xi, float) and isinstance(eta, float)
True

pytcl.gravity.create_test_coefficients(n_max=36)

Create test coefficients for low-degree testing.

Creates a simplified set of coefficients based on published low-degree values from EGM2008 for testing purposes.

Parameters:

n_max (int) – Maximum degree (default 36).

Returns:

Test coefficient set.

Return type:

EGMCoefficients

Examples

>>> coef = create_test_coefficients(n_max=10)
>>> coef.n_max
10
>>> coef.C[0, 0]  # Central term
1.0
>>> abs(coef.C[2, 0]) > 0  # J2 term present
True

class pytcl.gravity.TidalDisplacement(radial, north, east)

Bases: NamedTuple

Result of tidal displacement computation.

radial

Radial (up) displacement in meters.

Type:

float

north

North displacement in meters.

Type:

float

east

East displacement in meters.

Type:

float

radial: float

Alias for field number 0

north: float

Alias for field number 1

east: float

Alias for field number 2

class pytcl.gravity.TidalGravity(delta_g, delta_g_north, delta_g_east)

Bases: NamedTuple

Result of tidal gravity computation.

delta_g

Gravity change in m/s^2 (positive = increased gravity).

Type:

float

delta_g_north

North component of gravity change in m/s^2.

Type:

float

delta_g_east

East component of gravity change in m/s^2.

Type:

float

delta_g: float

Alias for field number 0

delta_g_north: float

Alias for field number 1

delta_g_east: float

Alias for field number 2

class pytcl.gravity.OceanTideLoading(amplitude, phase, constituents)

Bases: NamedTuple

Ocean tide loading parameters for a station.

amplitude

Amplitude for each constituent (m) for radial, west, south.

Type:

NDArray

phase

Phase for each constituent (radians).

Type:

NDArray

constituents

Names of tidal constituents.

Type:

Tuple[str, …]

amplitude: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

phase: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

constituents: tuple[str, ...]

Alias for field number 2

pytcl.gravity.julian_centuries_j2000(mjd)

Convert Modified Julian Date to Julian centuries since J2000.0.

Parameters:

mjd (float) – Modified Julian Date.

Returns:

T – Julian centuries since J2000.0.

Return type:

float

Examples

>>> T = julian_centuries_j2000(51544.5)  # J2000.0 epoch
>>> abs(T) < 1e-10  # Should be zero at J2000
True
>>> T = julian_centuries_j2000(51544.5 + 36525)  # One century later
>>> abs(T - 1.0) < 1e-10
True

pytcl.gravity.fundamental_arguments(T)

Compute fundamental astronomical arguments (Delaunay variables).

Parameters:

T (float) – Julian centuries since J2000.0.

Returns:

• l_moon (float) – Mean anomaly of the Moon (radians).

• l_sun (float) – Mean anomaly of the Sun (radians).

• F (float) – Mean argument of latitude of the Moon (radians).

• D (float) – Mean elongation of the Moon from the Sun (radians).

• Omega (float) – Mean longitude of the ascending node of the Moon (radians).

Return type:

Tuple[float, float, float, float, float]

Notes

Based on IERS Conventions (2010) expressions.

Examples

>>> import numpy as np
>>> T = 0.0  # J2000.0
>>> l, lp, F, D, Om = fundamental_arguments(T)
>>> all(0 <= x < 2*np.pi for x in [l, lp, F, D, Om])  # All in [0, 2pi)
True

pytcl.gravity.moon_position_approximate(mjd)

Compute approximate geocentric Moon position.

Parameters:

mjd (float) – Modified Julian Date.

Returns:

• r (float) – Distance from Earth center (m).

• lat (float) – Geocentric latitude (radians).

• lon (float) – Geocentric longitude (radians).

Return type:

Tuple[float, float, float]

Notes

Low-precision formula adequate for tidal computations.

Examples

>>> r, lat, lon = moon_position_approximate(58000)
>>> 350000e3 < r < 410000e3  # Distance in km range
True
>>> import numpy as np
>>> -np.pi/2 <= lat <= np.pi/2  # Valid latitude
True

pytcl.gravity.sun_position_approximate(mjd)

Compute approximate geocentric Sun position.

Parameters:

mjd (float) – Modified Julian Date.

Returns:

• r (float) – Distance from Earth center (m).

• lat (float) – Geocentric latitude (radians).

• lon (float) – Geocentric longitude (radians).

Return type:

Tuple[float, float, float]

Notes

Low-precision formula adequate for tidal computations.

Examples

>>> r, lat, lon = sun_position_approximate(58000)
>>> 1.47e11 < r < 1.52e11  # Distance ~1 AU
True
>>> lat == 0.0  # Sun on ecliptic
True

pytcl.gravity.solid_earth_tide_displacement(lat, lon, mjd, h2=0.6078, l2=0.0847)

Compute solid Earth tide displacement at a station.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• mjd (float) – Modified Julian Date.

• h2 (float, optional) – Love number h2. Default is IERS value.

• l2 (float, optional) – Shida number l2. Default is IERS value.

Returns:

Displacement in radial, north, east directions (meters).

Return type:

TidalDisplacement

Notes

Computes the degree-2 solid Earth tide displacement caused by the Moon and Sun. The displacement follows the tide-generating potential with Love and Shida numbers.

Examples

>>> import numpy as np
>>> # Displacement at 45N, 0E on MJD 58000
>>> disp = solid_earth_tide_displacement(np.radians(45), 0, 58000)
>>> abs(disp.radial) < 0.5  # Radial displacement typically < 50cm
True

pytcl.gravity.solid_earth_tide_gravity(lat, lon, mjd, delta=1.1608)

Compute solid Earth tide effect on gravity.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• mjd (float) – Modified Julian Date.

• delta (float, optional) – Gravimetric factor (1 + h - 3k/2). Default is IERS value.

Returns:

Gravity change in vertical and horizontal components (m/s^2).

Return type:

TidalGravity

Notes

The tidal gravity change is computed using the tide-generating potential and the gravimetric factor delta = 1 + h - 3k/2.

Examples

>>> import numpy as np
>>> grav = solid_earth_tide_gravity(np.radians(45), 0, 58000)
>>> abs(grav.delta_g) < 3e-6  # Typically < 3 microGal
True

pytcl.gravity.ocean_tide_loading_displacement(mjd, amplitude, phase, constituents=('M2', 'S2', 'N2', 'K2', 'K1', 'O1', 'P1', 'Q1'))

Compute ocean tide loading displacement.

Parameters:

• mjd (float) – Modified Julian Date.

• amplitude (NDArray) – Loading amplitudes (3, n_constituents) for radial, north, east. Units: meters.

• phase (NDArray) – Loading phases (3, n_constituents) in radians.

• constituents (Tuple[str, ...], optional) – Tidal constituent names. Default is major constituents.

Returns:

Ocean loading displacement (meters).

Return type:

TidalDisplacement

Notes

Ocean tide loading parameters are typically obtained from services like IERS or computed from ocean tide models (e.g., FES2014, GOT4.10).

The displacement is computed as: sum_i amplitude_i * cos(omega_i * t + chi_i - phase_i)

where omega_i is the constituent frequency and chi_i is the astronomical argument.

Examples

>>> import numpy as np
>>> # Example loading coefficients (typically from IERS)
>>> amp = np.array([[0.01, 0.005], [0.002, 0.001], [0.002, 0.001]])
>>> phase = np.array([[0.1, 0.2], [0.3, 0.4], [0.5, 0.6]])
>>> disp = ocean_tide_loading_displacement(58000, amp, phase, ("M2", "S2"))
>>> isinstance(disp, TidalDisplacement)
True

pytcl.gravity.atmospheric_pressure_loading(lat, lon, pressure, reference_pressure=101325.0, admittance=-0.00035)

Compute atmospheric pressure loading displacement.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• pressure (float) – Surface atmospheric pressure in Pascals.

• reference_pressure (float, optional) – Reference pressure in Pascals. Default is 101325 Pa (1 atm).

• admittance (float, optional) – Pressure admittance in m/Pa. Default is -0.35 mm/hPa.

Returns:

Pressure loading displacement (meters).

Return type:

TidalDisplacement

Notes

Atmospheric pressure loading causes the Earth’s surface to deform under changing atmospheric mass. The effect is primarily radial (vertical) with magnitude approximately -0.35 mm per hPa of pressure above the mean.

For precise applications, use Green’s functions from models like ECMWF or NCEP pressure fields.

Examples

>>> import numpy as np
>>> # High pressure (1020 hPa) above reference
>>> disp = atmospheric_pressure_loading(np.radians(45), 0, 102000)
>>> disp.radial < 0  # Surface depressed under high pressure
True

pytcl.gravity.pole_tide_displacement(lat, lon, xp, yp, xp_mean=0.0, yp_mean=0.0, h2=0.6078, l2=0.0847)

Compute pole tide (polar motion) displacement.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• xp (float) – Pole x coordinate in arcseconds.

• yp (float) – Pole y coordinate in arcseconds.

• xp_mean (float, optional) – Mean pole x coordinate in arcseconds. Default 0.

• yp_mean (float, optional) – Mean pole y coordinate in arcseconds. Default 0.

• h2 (float, optional) – Love number. Default is IERS value.

• l2 (float, optional) – Shida number. Default is IERS value.

Returns:

Pole tide displacement (meters).

Return type:

TidalDisplacement

Notes

The pole tide results from the centrifugal effect of polar motion. The Earth’s rotation axis wobbles with a primary period of ~433 days (Chandler wobble) causing small displacements.

Examples

>>> import numpy as np
>>> # Pole offset of 0.1 arcsec
>>> disp = pole_tide_displacement(np.radians(45), 0, 0.1, 0.1)
>>> abs(disp.radial) < 0.03  # Typically < 3cm
True

pytcl.gravity.total_tidal_displacement(lat, lon, mjd, ocean_loading=None, pressure=None, xp=0.0, yp=0.0)

Compute total tidal displacement from all sources.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• mjd (float) – Modified Julian Date.

• ocean_loading (OceanTideLoading, optional) – Ocean loading parameters. If None, ocean loading is not included.

• pressure (float, optional) – Surface pressure in Pascals. If None, pressure loading is not included.

• xp (float, optional) – Pole x coordinate in arcseconds. Default 0.

• yp (float, optional) – Pole y coordinate in arcseconds. Default 0.

Returns:

Total tidal displacement (meters).

Return type:

TidalDisplacement

Examples

>>> import numpy as np
>>> disp = total_tidal_displacement(np.radians(45), 0, 58000)
>>> isinstance(disp, TidalDisplacement)
True

pytcl.gravity.tidal_gravity_correction(lat, lon, mjd)

Compute tidal correction to gravity observations.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• mjd (float) – Modified Julian Date.

Returns:

delta_g – Gravity correction in m/s^2 (add to observed gravity).

Return type:

float

Notes

This correction should be added to absolute gravity observations to remove the tidal signal.

Examples

>>> import numpy as np
>>> corr = tidal_gravity_correction(np.radians(45), 0, 58000)
>>> abs(corr) < 3e-6  # Typically < 300 microGal
True

Spherical Harmonics

Spherical harmonic gravity field representation and evaluation.

Spherical harmonic functions for geophysical models.

Spherical harmonics are used to represent gravitational and magnetic fields on a sphere. This module provides functions for evaluating associated Legendre polynomials and spherical harmonic expansions.

References

[1]

23. 1. Heiskanen and H. Moritz, “Physical Geodesy,” W. H. Freeman, 1967.

[2]

15. Montenbruck and E. Gill, “Satellite Orbits,” Springer, 2000.

pytcl.gravity.spherical_harmonics.associated_legendre(n_max, m_max, x, normalized=True)

Compute associated Legendre polynomials P_n^m(x).

Uses the recursive algorithm for numerical stability.

Parameters:

• n_max (int) – Maximum degree.

• m_max (int) – Maximum order (must be <= n_max).

• x (float) – Argument, typically cos(colatitude). Must be in [-1, 1].

• normalized (bool, optional) – If True, return fully normalized (geodetic) coefficients. Default True.

Returns:

P – Array of shape (n_max+1, m_max+1) containing P_n^m(x).

Return type:

ndarray

Notes

The fully normalized associated Legendre functions satisfy:

\[\int_{-1}^{1} [\bar{P}_n^m(x)]^2 dx = \frac{2}{2n+1}\]

Results are cached for repeated queries with the same parameters. Cache key quantizes x to 8 decimal places (~1e-8 precision).

Examples

>>> P = associated_legendre(2, 2, 0.5)
>>> P[2, 0]  # P_2^0(0.5)

pytcl.gravity.spherical_harmonics.associated_legendre_derivative(n_max, m_max, x, P=None, normalized=True)

Compute derivatives of associated Legendre polynomials dP_n^m/dx.

Parameters:

• n_max (int) – Maximum degree.

• m_max (int) – Maximum order.

• x (float) – Argument in [-1, 1].

• P (ndarray, optional) – Precomputed P_n^m values. If None, computed internally.

• normalized (bool, optional) – If True, use fully normalized functions. Default True.

Returns:

dP – Array of shape (n_max+1, m_max+1) containing dP_n^m/dx.

Return type:

ndarray

Examples

>>> import numpy as np
>>> x = np.cos(np.radians(45))  # cos of 45 degrees
>>> dP = associated_legendre_derivative(2, 2, x)
>>> dP.shape
(3, 3)
>>> abs(dP[0, 0]) < 1e-10  # dP_0^0/dx = 0
True

pytcl.gravity.spherical_harmonics.spherical_harmonic_sum(lat, lon, r, C, S, R, GM, n_max=None)

Evaluate spherical harmonic expansion for a scalar field.

Computes the value and gradient of a field represented by spherical harmonic coefficients C_nm and S_nm.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• r (float) – Radial distance from center of mass.

• C (ndarray) – Cosine coefficients C_nm, shape (n_max+1, n_max+1).

• S (ndarray) – Sine coefficients S_nm, shape (n_max+1, n_max+1).

• R (float) – Reference radius (e.g., Earth’s equatorial radius).

• GM (float) – Gravitational parameter (G * M).

• n_max (int, optional) – Maximum degree to use. Default uses full coefficient array.

Returns:

• V (float) – Potential value.

• dV_r (float) – Radial derivative dV/dr.

• dV_lat (float) – Latitudinal derivative (1/r) * dV/dlat.

Return type:

Tuple[float, float, float]

Notes

The spherical harmonic expansion of the gravitational potential is:

\[V = \frac{GM}{r} \sum_{n=0}^{N} \left(\frac{R}{r}\right)^n \sum_{m=0}^{n} \bar{P}_n^m(\sin\phi) (C_{nm}\cos m\lambda + S_{nm}\sin m\lambda)\]

Examples

>>> import numpy as np
>>> # Simple monopole (degree 0 only)
>>> C = np.array([[1.0]])
>>> S = np.array([[0.0]])
>>> R = 6.378e6  # meters
>>> GM = 3.986e14  # m^3/s^2
>>> V, dV_r, dV_lat = spherical_harmonic_sum(0, 0, R, C, S, R, GM, n_max=0)
>>> abs(V - GM/R) / (GM/R) < 1e-10  # V = GM/r for degree 0
True

pytcl.gravity.spherical_harmonics.gravity_acceleration(lat, lon, h, C, S, R, GM, n_max=None)

Compute gravity acceleration vector from spherical harmonics.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float) – Height above reference ellipsoid in meters.

• C (ndarray) – Cosine coefficients.

• S (ndarray) – Sine coefficients.

• R (float) – Reference radius.

• GM (float) – Gravitational parameter.

• n_max (int, optional) – Maximum degree.

Returns:

• g_r (float) – Radial component of gravity (positive outward).

• g_lat (float) – Northward component of gravity.

• g_lon (float) – Eastward component of gravity.

Return type:

Tuple[float, float, float]

Examples

>>> import numpy as np
>>> # Simple monopole field
>>> C = np.array([[1.0]])
>>> S = np.array([[0.0]])
>>> R = 6.378e6
>>> GM = 3.986e14
>>> g_r, g_lat, g_lon = gravity_acceleration(0, 0, 0, C, S, R, GM, n_max=0)
>>> g_r < 0  # Gravity points inward (negative radial)
True

pytcl.gravity.spherical_harmonics.legendre_scaling_factors(n_max)

Precompute scaling factors to prevent overflow in Legendre recursion.

For degrees > ~150, standard Legendre recursion can overflow in double precision. These scaling factors keep intermediate values in a representable range.

The scaling follows the approach of Holmes & Featherstone (2002), using a factor that grows with degree to counteract the natural growth of the Legendre functions.

Parameters:

n_max (int) – Maximum degree.

Returns:

scale – Scaling factors of shape (n_max+1,). The factor for degree n is 10^(-280 * n / n_max) for n_max > 150, else 1.0.

Return type:

ndarray

References

[1]

Holmes, S.A. and Featherstone, W.E. “A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions.” Journal of Geodesy 76.5 (2002): 279-299.

Examples

>>> scale = legendre_scaling_factors(100)
>>> len(scale)
101
>>> scale[0]  # No scaling for low degrees
1.0

>>> scale_high = legendre_scaling_factors(200)
>>> scale_high[200] < scale_high[0]  # Higher degrees scaled down
True

pytcl.gravity.spherical_harmonics.associated_legendre_scaled(n_max, m_max, x, scale=None)

Compute scaled associated Legendre polynomials for high degrees.

For ultra-high degree computations (n > 150), standard Legendre recursion overflows. This function computes scaled values that stay in representable range.

Parameters:

• n_max (int) – Maximum degree.

• m_max (int) – Maximum order (must be <= n_max).

• x (float) – Argument in [-1, 1], typically cos(colatitude).

• scale (ndarray, optional) – Precomputed scaling factors from legendre_scaling_factors(). If None, computed internally.

Returns:

• P_scaled (ndarray) – Scaled Legendre values, shape (n_max+1, m_max+1). The actual value is P_scaled[n,m] * 10^scale_exp[n].

• scale_exp (ndarray) – Scale exponents for each degree, shape (n_max+1,). Set to 0 if no scaling needed.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Notes

The returned values satisfy:

P_n^m(x) = P_scaled[n,m] * 10^scale_exp[n]

For normal operations (n < 150), scale_exp is all zeros and P_scaled equals the actual Legendre values.

Examples

>>> import numpy as np
>>> x = np.cos(np.radians(45))
>>> P_scaled, scale_exp = associated_legendre_scaled(10, 10, x)
>>> P_scaled.shape
(11, 11)
>>> all(scale_exp == 0)  # No scaling needed for n_max < 150
True

pytcl.gravity.spherical_harmonics.clear_legendre_cache()

Clear cached Legendre polynomial results.

Call this function to clear the cached associated Legendre polynomial arrays. Useful when memory is constrained or after processing a batch with different colatitude values.

Examples

>>> _ = associated_legendre(10, 10, 0.5)  # Populate cache
>>> info = get_legendre_cache_info()
>>> clear_legendre_cache()
>>> info_after = get_legendre_cache_info()
>>> info_after.currsize
0

pytcl.gravity.spherical_harmonics.get_legendre_cache_info()

Get cache statistics for Legendre polynomials.

Returns:

Named tuple with hits, misses, maxsize, currsize.

Return type:

CacheInfo

Examples

>>> clear_legendre_cache()  # Start fresh
>>> _ = associated_legendre(5, 5, 0.5)
>>> info = get_legendre_cache_info()
>>> info.currsize >= 1  # At least one entry cached
True

EGM Models

Earth Gravitational Model (EGM96/EGM2008) implementations.

Earth Gravitational Models (EGM96 and EGM2008).

This module provides support for high-degree Earth gravitational models including EGM96 (degree 360) and EGM2008 (degree 2190). These models represent the Earth’s gravitational field using spherical harmonic coefficients.

The models support: - Geoid height computation - Gravity disturbance/anomaly calculation - Deflection of the vertical

References

[1]

Lemoine, F.G., et al. “The Development of the Joint NASA GSFC and NIMA Geopotential Model EGM96.” NASA Technical Paper, 1998.

[2]

Pavlis, N.K., et al. “The development and evaluation of the Earth Gravitational Model 2008 (EGM2008).” JGR 117.B4 (2012).

[3]

National Geospatial-Intelligence Agency. “EGM2008 Model Coefficients.” https://earth-info.nga.mil/

class pytcl.gravity.egm.EGMCoefficients(C, S, GM, R, n_max, model_name)

Bases: NamedTuple

Earth Gravitational Model coefficients.

Parameters:

• C (ndarray) – Cosine coefficients (fully normalized), shape (n_max+1, n_max+1).

• S (ndarray) – Sine coefficients (fully normalized), shape (n_max+1, n_max+1).

• GM (float) – Gravitational parameter in m^3/s^2.

• R (float) – Reference radius in meters.

• n_max (int) – Maximum degree of the model.

• model_name (str) – Model name (“EGM96” or “EGM2008”).

C: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

S: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

GM: float

Alias for field number 2

R: float

Alias for field number 3

n_max: int

Alias for field number 4

model_name: str

Alias for field number 5

class pytcl.gravity.egm.GeoidResult(height, lat, lon, model)

Bases: NamedTuple

Geoid height computation result.

Parameters:

• height (float) – Geoid height above reference ellipsoid in meters.

• lat (float) – Latitude in radians.

• lon (float) – Longitude in radians.

• model (str) – Model name used.

height: float

Alias for field number 0

lat: float

Alias for field number 1

lon: float

Alias for field number 2

model: str

Alias for field number 3

class pytcl.gravity.egm.GravityDisturbance(delta_g_r, delta_g_lat, delta_g_lon, magnitude)

Bases: NamedTuple

Gravity disturbance result.

Parameters:

• delta_g_r (float) – Radial component of gravity disturbance in m/s^2.

• delta_g_lat (float) – Northward component in m/s^2.

• delta_g_lon (float) – Eastward component in m/s^2.

• magnitude (float) – Total magnitude of disturbance in m/s^2.

delta_g_r: float

Alias for field number 0

delta_g_lat: float

Alias for field number 1

delta_g_lon: float

Alias for field number 2

magnitude: float

Alias for field number 3

pytcl.gravity.egm.get_data_dir()

Get the pytcl data directory for external data files.

The data directory is located at ~/.pytcl/data/ by default. Can be overridden by setting the PYTCL_DATA_DIR environment variable.

Returns:

Path to the data directory.

Return type:

Path

pytcl.gravity.egm.parse_egm_file(filepath, n_max=None)

Parse an EGM coefficient file in NGA format.

The file format is:

n m C_nm S_nm [sigma_C sigma_S]

where n is degree, m is order, C_nm and S_nm are the normalized coefficients, and sigma values are optional uncertainties.

Parameters:

• filepath (Path) – Path to the coefficient file.

• n_max (int, optional) – Maximum degree to load. If None, loads all coefficients.

Returns:

• C (ndarray) – Cosine coefficients, shape (n_max+1, n_max+1).

• S (ndarray) – Sine coefficients, shape (n_max+1, n_max+1).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

pytcl.gravity.egm.create_test_coefficients(n_max=36)

Create test coefficients for low-degree testing.

Creates a simplified set of coefficients based on published low-degree values from EGM2008 for testing purposes.

Parameters:

n_max (int) – Maximum degree (default 36).

Returns:

Test coefficient set.

Return type:

EGMCoefficients

Examples

>>> coef = create_test_coefficients(n_max=10)
>>> coef.n_max
10
>>> coef.C[0, 0]  # Central term
1.0
>>> abs(coef.C[2, 0]) > 0  # J2 term present
True

pytcl.gravity.egm.load_egm_coefficients(model='EGM2008', n_max=None)

Load EGM coefficients from file.

Parameters:

• model (str) – Model name, either “EGM96” or “EGM2008”.

• n_max (int, optional) – Maximum degree to load. If None, loads the full model.

Returns:

Loaded coefficient structure.

Return type:

EGMCoefficients

Raises:

• FileNotFoundError – If the coefficient file is not found.

• ValueError – If an unknown model is specified.

Examples

>>> coef = load_egm_coefficients("EGM2008", n_max=360)
>>> coef.n_max
360

pytcl.gravity.egm.geoid_height(lat, lon, model='EGM2008', n_max=None, coefficients=None)

Compute geoid height above WGS84 ellipsoid.

The geoid is the equipotential surface of the Earth’s gravity field that best fits mean sea level. This function computes the height of the geoid above the WGS84 reference ellipsoid.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• model (str) – Model name (“EGM96” or “EGM2008”).

• n_max (int, optional) – Maximum degree to use. Default uses full model.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients. If None, loads from file.

Returns:

Geoid height in meters.

Return type:

float

Examples

>>> # Geoid height at equator, prime meridian
>>> N = geoid_height(0, 0)  # Should be approximately 17 m

pytcl.gravity.egm.geoid_heights(lats, lons, model='EGM2008', n_max=None, coefficients=None)

Compute geoid heights for multiple points.

Parameters:

• lats (ndarray) – Geodetic latitudes in radians.

• lons (ndarray) – Longitudes in radians.

• model (str) – Model name.

• n_max (int, optional) – Maximum degree.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients. If None, loads from file.

Returns:

Geoid heights in meters.

Return type:

ndarray

Examples

>>> import numpy as np
>>> coef = create_test_coefficients(n_max=10)
>>> lats = np.array([0.0, np.pi/4])  # Equator and 45°N
>>> lons = np.array([0.0, np.pi/2])  # Prime meridian and 90°E
>>> heights = geoid_heights(lats, lons, coefficients=coef)
>>> len(heights)
2

pytcl.gravity.egm.gravity_disturbance(lat, lon, h=0.0, model='EGM2008', n_max=None, coefficients=None)

Compute gravity disturbance from EGM model.

The gravity disturbance is the difference between actual gravity and normal (reference ellipsoid) gravity at the same point.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float) – Height above ellipsoid in meters.

• model (str) – Model name.

• n_max (int, optional) – Maximum degree.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients.

Returns:

Gravity disturbance components.

Return type:

GravityDisturbance

Examples

>>> coef = create_test_coefficients(n_max=10)
>>> dist = gravity_disturbance(0, 0, h=0, coefficients=coef)
>>> isinstance(dist.magnitude, float)
True

pytcl.gravity.egm.gravity_anomaly(lat, lon, h=0.0, model='EGM2008', n_max=None, coefficients=None)

Compute free-air gravity anomaly.

The gravity anomaly is the difference between observed gravity and normal gravity, with the normal gravity evaluated at the geoid rather than at the observation point.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float) – Height above ellipsoid in meters.

• model (str) – Model name.

• n_max (int, optional) – Maximum degree.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients. If None, loads from file.

Returns:

Gravity anomaly in m/s^2 (typically reported in mGal = 1e-5 m/s^2).

Return type:

float

Examples

>>> coef = create_test_coefficients(n_max=10)
>>> anomaly = gravity_anomaly(0, 0, h=0, coefficients=coef)
>>> isinstance(anomaly, float)
True

pytcl.gravity.egm.deflection_of_vertical(lat, lon, model='EGM2008', n_max=None, coefficients=None)

Compute deflection of the vertical.

The deflection of the vertical is the angle between the direction of gravity (plumb line) and the normal to the reference ellipsoid.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• model (str) – Model name.

• n_max (int, optional) – Maximum degree.

• coefficients (EGMCoefficients, optional) – Pre-loaded coefficients. If None, loads from file.

Returns:

• xi (float) – North-south deflection component in arcseconds.

• eta (float) – East-west deflection component in arcseconds.

Return type:

Tuple[float, float]

Notes

Positive xi means the plumb line points more north than the normal. Positive eta means the plumb line points more east than the normal.

Examples

>>> coef = create_test_coefficients(n_max=10)
>>> xi, eta = deflection_of_vertical(0, 0, coefficients=coef)
>>> isinstance(xi, float) and isinstance(eta, float)
True

Gravity Models

Standard gravity models (WGS84, J2, point mass).

Gravity field models.

This module provides implementations of standard gravity models including WGS84 normal gravity and spherical harmonic gravity models.

References

[1]

NIMA, “Department of Defense World Geodetic System 1984,” TR8350.2, 2000.

[2]

Pavlis et al., “The development and evaluation of the Earth Gravitational Model 2008 (EGM2008),” JGR, 2012.

class pytcl.gravity.models.GravityConstants(GM, a, f, omega, J2)

Bases: NamedTuple

Constants for a gravity model.

GM

Gravitational parameter (m^3/s^2).

Type:

float

a

Semi-major axis / equatorial radius (m).

Type:

float

f

Flattening.

Type:

float

omega

Angular velocity (rad/s).

Type:

float

J2

Second degree zonal harmonic (unnormalized).

Type:

float

GM: float

Alias for field number 0

a: float

Alias for field number 1

f: float

Alias for field number 2

omega: float

Alias for field number 3

J2: float

Alias for field number 4

class pytcl.gravity.models.GravityResult(magnitude, g_down, g_north, g_east)

Bases: NamedTuple

Result of gravity computation.

magnitude

Total gravity magnitude (m/s^2).

Type:

float

g_down

Downward component (positive down) (m/s^2).

Type:

float

g_north

Northward component (m/s^2).

Type:

float

g_east

Eastward component (m/s^2).

Type:

float

magnitude: float

Alias for field number 0

g_down: float

Alias for field number 1

g_north: float

Alias for field number 2

g_east: float

Alias for field number 3

pytcl.gravity.models.normal_gravity_somigliana(lat, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute normal gravity on the ellipsoid using Somigliana’s formula.

Parameters:

• lat (float) – Geodetic latitude in radians.

• constants (GravityConstants, optional) – Gravity model constants. Default WGS84.

Returns:

gamma – Normal gravity on the ellipsoid surface (m/s^2).

Return type:

float

Notes

Somigliana’s closed formula gives the exact normal gravity on the reference ellipsoid without iteration.

Examples

>>> gamma = normal_gravity_somigliana(0)  # At equator
>>> abs(gamma - 9.7803) < 0.001
True

>>> gamma = normal_gravity_somigliana(np.pi/2)  # At pole
>>> abs(gamma - 9.8322) < 0.001
True

pytcl.gravity.models.normal_gravity(lat, h=0.0, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute normal gravity at a given latitude and height.

Uses the International Gravity Formula with free-air correction.

Parameters:

• lat (float) – Geodetic latitude in radians.

• h (float, optional) – Height above ellipsoid in meters. Default 0.

• constants (GravityConstants, optional) – Gravity model constants. Default WGS84.

Returns:

g – Normal gravity (m/s^2).

Return type:

float

Examples

>>> g = normal_gravity(0, 0)  # Sea level at equator
>>> abs(g - 9.78) < 0.01
True

>>> g = normal_gravity(np.radians(45), 1000)  # 1km altitude, 45° lat
>>> g < normal_gravity(np.radians(45), 0)  # Decreases with altitude
True

pytcl.gravity.models.gravity_wgs84(lat, lon, h=0.0)

Compute gravity using WGS84 model.

This computes the full gravity vector including the centrifugal acceleration due to Earth’s rotation.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float, optional) – Height above WGS84 ellipsoid in meters. Default 0.

Returns:

result – Gravity components and magnitude.

Return type:

GravityResult

Examples

>>> result = gravity_wgs84(0, 0, 0)
>>> abs(result.magnitude - 9.78) < 0.01
True

pytcl.gravity.models.gravity_j2(lat, lon, h=0.0, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute gravity using J2 (oblateness) model.

This simplified model includes only the J2 zonal harmonic, which accounts for Earth’s equatorial bulge.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• h (float, optional) – Height above ellipsoid in meters. Default 0.

• constants (GravityConstants, optional) – Model constants. Default WGS84.

Returns:

result – Gravity components and magnitude.

Return type:

GravityResult

Examples

>>> result = gravity_j2(0, 0, 0)  # At equator, sea level
>>> abs(result.magnitude - 9.78) < 0.01
True

pytcl.gravity.models.geoid_height_j2(lat, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute approximate geoid height using J2 model.

Parameters:

• lat (float) – Geodetic latitude in radians.

• constants (GravityConstants, optional) – Model constants. Default WGS84.

Returns:

N – Geoid height (geoid - ellipsoid) in meters.

Return type:

float

Examples

>>> N = geoid_height_j2(0)  # At equator
>>> N > 0  # Equator bulges outward
True

Notes

This is a simplified model. For accurate geoid heights, use a full geoid model like EGM2008.

pytcl.gravity.models.gravitational_potential(lat, lon, r, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263), n_max=2)

Compute gravitational potential.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• r (float) – Radial distance from Earth’s center in meters.

• constants (GravityConstants, optional) – Model constants. Default WGS84.

• n_max (int, optional) – Maximum spherical harmonic degree. Default 2.

Returns:

U – Gravitational potential (m^2/s^2).

Return type:

float

Examples

>>> U = gravitational_potential(0, 0, 6.4e6)  # Near Earth surface
>>> U < 0  # Potential is negative
True

pytcl.gravity.models.free_air_anomaly(g_observed, lat, h, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute free-air gravity anomaly.

Parameters:

• g_observed (float) – Observed gravity in m/s^2.

• lat (float) – Geodetic latitude in radians.

• h (float) – Height above geoid in meters.

• constants (GravityConstants, optional) – Model constants. Default WGS84.

Returns:

delta_g – Free-air anomaly in m/s^2 (or mGal if multiplied by 1e5).

Return type:

float

Examples

>>> import numpy as np
>>> # Observed gravity slightly higher than normal
>>> delta_g = free_air_anomaly(9.81, np.radians(45), 100)
>>> isinstance(delta_g, float)
True

Notes

The free-air anomaly is the difference between observed gravity and normal gravity at the observation point.

pytcl.gravity.models.bouguer_anomaly(g_observed, lat, h, rho=2670.0, constants=(398600441800000.0, 6378137.0, 0.0033528106647474805, 7.292115e-05, 0.00108263))

Compute simple Bouguer gravity anomaly.

Parameters:

• g_observed (float) – Observed gravity in m/s^2.

• lat (float) – Geodetic latitude in radians.

• h (float) – Height above geoid in meters.

• rho (float, optional) – Crustal density in kg/m^3. Default 2670 (average crustal density).

• constants (GravityConstants, optional) – Model constants. Default WGS84.

Returns:

delta_g – Bouguer anomaly in m/s^2.

Return type:

float

Examples

>>> import numpy as np
>>> # Bouguer anomaly at mountain location
>>> delta_g = bouguer_anomaly(9.81, np.radians(45), 1000)
>>> isinstance(delta_g, float)
True

Notes

The Bouguer anomaly removes the gravitational effect of the topographic mass between the observation point and the geoid.

Clenshaw Summation

Efficient Clenshaw summation algorithm for spherical harmonics.

Clenshaw summation for efficient spherical harmonic evaluation.

The Clenshaw algorithm evaluates spherical harmonic series efficiently via backward recursion, avoiding direct computation of associated Legendre functions which can overflow at high degrees.

This implementation follows Holmes & Featherstone (2002) for numerical stability at ultra-high degrees (n > 2000).

Performance Notes

Recursion coefficients (_a_nm, _b_nm) are cached using lru_cache for 25-40% speedup on repeated evaluations with the same (n, m) pairs.

References

[1]

Holmes, S.A. and Featherstone, W.E. “A unified approach to the Clenshaw summation and the recursive computation of very high degree and order normalised associated Legendre functions.” Journal of Geodesy 76.5 (2002): 279-299.

[2]

Wittwer, T., et al. “Ultra-high degree spherical harmonic analysis and synthesis using extended-range arithmetic.” Journal of Geodesy 82.4-5 (2008): 223-229.

pytcl.gravity.clenshaw.clenshaw_sum_order(m, cos_theta, sin_theta, C, S, n_max)

Clenshaw summation for fixed order m, summing over degrees n=m to n_max.

Evaluates the partial sums:

sum_C = sum_{n=m}^{n_max} C[n,m] * P_n^m(cos_theta) sum_S = sum_{n=m}^{n_max} S[n,m] * P_n^m(cos_theta)

using backward recursion from n_max down to m.

Parameters:

• m (int) – Order (fixed for this summation).

• cos_theta (float) – Cosine of colatitude.

• sin_theta (float) – Sine of colatitude.

• C (ndarray) – Cosine coefficients array, shape (n_max+1, n_max+1).

• S (ndarray) – Sine coefficients array, shape (n_max+1, n_max+1).

• n_max (int) – Maximum degree.

Returns:

• sum_C (float) – Sum of C terms weighted by Legendre functions.

• sum_S (float) – Sum of S terms weighted by Legendre functions.

Return type:

Tuple[float, float]

Examples

>>> import numpy as np
>>> C = np.zeros((5, 5))
>>> S = np.zeros((5, 5))
>>> C[2, 0] = 1.0  # Only C20 term
>>> cos_theta, sin_theta = np.cos(np.pi/4), np.sin(np.pi/4)
>>> sum_C, sum_S = clenshaw_sum_order(0, cos_theta, sin_theta, C, S, 4)
>>> isinstance(sum_C, float)
True

pytcl.gravity.clenshaw.clenshaw_sum_order_derivative(m, cos_theta, sin_theta, C, S, n_max)

Clenshaw summation with derivative for fixed order m.

Evaluates both the partial sums and their derivatives with respect to colatitude.

Parameters:

• m (int) – Order.

• cos_theta (float) – Cosine of colatitude.

• sin_theta (float) – Sine of colatitude.

• C (ndarray) – Cosine coefficients.

• S (ndarray) – Sine coefficients.

• n_max (int) – Maximum degree.

Returns:

• sum_C (float) – Sum of C terms.

• sum_S (float) – Sum of S terms.

• dsum_C (float) – Derivative of sum_C with respect to theta.

• dsum_S (float) – Derivative of sum_S with respect to theta.

Return type:

Tuple[float, float, float, float]

Examples

>>> import numpy as np
>>> C = np.zeros((5, 5))
>>> S = np.zeros((5, 5))
>>> C[2, 0] = -0.0005  # J2-like term
>>> cos_theta, sin_theta = np.cos(np.pi/4), np.sin(np.pi/4)
>>> sum_C, sum_S, dsum_C, dsum_S = clenshaw_sum_order_derivative(
...     0, cos_theta, sin_theta, C, S, 4)
>>> len([sum_C, sum_S, dsum_C, dsum_S])
4

pytcl.gravity.clenshaw.clenshaw_geoid(lat, lon, C, S, R, GM, gamma, n_max=None)

Compute geoid height using Clenshaw summation.

The geoid height N is the height of the geoid above the reference ellipsoid, computed from the disturbing potential T:

N = T / gamma

where gamma is the normal gravity on the ellipsoid.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• C (ndarray) – Cosine coefficients (fully normalized), shape (n_max+1, n_max+1).

• S (ndarray) – Sine coefficients (fully normalized), shape (n_max+1, n_max+1).

• R (float) – Reference radius in meters.

• GM (float) – Gravitational parameter in m^3/s^2.

• gamma (float) – Normal gravity at the evaluation point in m/s^2.

• n_max (int, optional) – Maximum degree to use. Default uses full coefficient array.

Returns:

Geoid height in meters.

Return type:

float

Notes

The geoid height is computed as:

\[N = \frac{GM}{r \gamma} \sum_{n=2}^{n_{max}} \left(\frac{R}{r}\right)^n \sum_{m=0}^{n} P_n^m(\sin\phi) (C_{nm}\cos m\lambda + S_{nm}\sin m\lambda)\]

The n=0 and n=1 terms are excluded as they represent the reference field.

Examples

>>> import numpy as np
>>> C = np.zeros((5, 5))
>>> S = np.zeros((5, 5))
>>> C[0, 0] = 1.0
>>> R = 6.378e6
>>> GM = 3.986e14
>>> gamma = 9.81
>>> N = clenshaw_geoid(0, 0, C, S, R, GM, gamma)
>>> isinstance(N, float)
True

pytcl.gravity.clenshaw.clenshaw_potential(lat, lon, r, C, S, R, GM, n_max=None)

Compute gravitational potential using Clenshaw summation.

Evaluates the spherical harmonic expansion of the gravitational potential efficiently using Clenshaw’s algorithm.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• r (float) – Radial distance from Earth center in meters.

• C (ndarray) – Cosine coefficients (fully normalized).

• S (ndarray) – Sine coefficients (fully normalized).

• R (float) – Reference radius in meters.

• GM (float) – Gravitational parameter in m^3/s^2.

• n_max (int, optional) – Maximum degree.

Returns:

Gravitational potential in m^2/s^2.

Return type:

float

Examples

>>> import numpy as np
>>> C = np.zeros((5, 5))
>>> S = np.zeros((5, 5))
>>> C[0, 0] = 1.0  # Central term only
>>> R = 6.378e6
>>> GM = 3.986e14
>>> V = clenshaw_potential(0, 0, R, C, S, R, GM)
>>> abs(V - GM/R) / (GM/R) < 0.01  # ~GM/r for central term
True

pytcl.gravity.clenshaw.clenshaw_gravity(lat, lon, r, C, S, R, GM, n_max=None)

Compute gravity disturbance vector using Clenshaw summation.

Evaluates both the potential and its gradient efficiently using Clenshaw’s algorithm with derivative recursions.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• r (float) – Radial distance from Earth center in meters.

• C (ndarray) – Cosine coefficients (fully normalized).

• S (ndarray) – Sine coefficients (fully normalized).

• R (float) – Reference radius in meters.

• GM (float) – Gravitational parameter in m^3/s^2.

• n_max (int, optional) – Maximum degree.

Returns:

• g_r (float) – Radial component of gravity disturbance (positive outward) in m/s^2.

• g_lat (float) – Northward component of gravity disturbance in m/s^2.

• g_lon (float) – Eastward component of gravity disturbance in m/s^2.

Return type:

Tuple[float, float, float]

Examples

>>> import numpy as np
>>> C = np.zeros((5, 5))
>>> S = np.zeros((5, 5))
>>> C[0, 0] = 1.0
>>> R = 6.378e6
>>> GM = 3.986e14
>>> g_r, g_lat, g_lon = clenshaw_gravity(0, 0, R, C, S, R, GM)
>>> g_r < 0  # Gravity points inward
True

Tidal Effects

Solid Earth tides, ocean tide loading, and pole tide corrections.

Tidal Effects on Gravity and Position.

This module provides functions for computing tidal effects including: - Solid Earth tides (body tides) - Ocean tide loading - Atmospheric pressure loading - Pole tide effects

These effects cause periodic deformations of the Earth’s surface and variations in the gravity field.

References

[1]

Petit, G. and Luzum, B. (eds.), “IERS Conventions (2010),” IERS Technical Note No. 36, Frankfurt am Main, 2010.

[2]

McCarthy, D.D. and Petit, G. (eds.), “IERS Conventions (2003),” IERS Technical Note No. 32, Frankfurt am Main, 2004.

[3]

Mathews, P.M., Dehant, V., and Gipson, J.M., “Tidal station displacements,” JGR, 102, B9, 20469-20477, 1997.

class pytcl.gravity.tides.TidalDisplacement(radial, north, east)

Bases: NamedTuple

Result of tidal displacement computation.

radial

Radial (up) displacement in meters.

Type:

float

north

North displacement in meters.

Type:

float

east

East displacement in meters.

Type:

float

radial: float

Alias for field number 0

north: float

Alias for field number 1

east: float

Alias for field number 2

class pytcl.gravity.tides.TidalGravity(delta_g, delta_g_north, delta_g_east)

Bases: NamedTuple

Result of tidal gravity computation.

delta_g

Gravity change in m/s^2 (positive = increased gravity).

Type:

float

delta_g_north

North component of gravity change in m/s^2.

Type:

float

delta_g_east

East component of gravity change in m/s^2.

Type:

float

delta_g: float

Alias for field number 0

delta_g_north: float

Alias for field number 1

delta_g_east: float

Alias for field number 2

class pytcl.gravity.tides.OceanTideLoading(amplitude, phase, constituents)

Bases: NamedTuple

Ocean tide loading parameters for a station.

amplitude

Amplitude for each constituent (m) for radial, west, south.

Type:

NDArray

phase

Phase for each constituent (radians).

Type:

NDArray

constituents

Names of tidal constituents.

Type:

Tuple[str, …]

amplitude: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 0

phase: ndarray[tuple[Any, ...], dtype[floating]]

Alias for field number 1

constituents: tuple[str, ...]

Alias for field number 2

pytcl.gravity.tides.julian_centuries_j2000(mjd)

Convert Modified Julian Date to Julian centuries since J2000.0.

Parameters:

mjd (float) – Modified Julian Date.

Returns:

T – Julian centuries since J2000.0.

Return type:

float

Examples

>>> T = julian_centuries_j2000(51544.5)  # J2000.0 epoch
>>> abs(T) < 1e-10  # Should be zero at J2000
True
>>> T = julian_centuries_j2000(51544.5 + 36525)  # One century later
>>> abs(T - 1.0) < 1e-10
True

pytcl.gravity.tides.fundamental_arguments(T)

Compute fundamental astronomical arguments (Delaunay variables).

Parameters:

T (float) – Julian centuries since J2000.0.

Returns:

• l_moon (float) – Mean anomaly of the Moon (radians).

• l_sun (float) – Mean anomaly of the Sun (radians).

• F (float) – Mean argument of latitude of the Moon (radians).

• D (float) – Mean elongation of the Moon from the Sun (radians).

• Omega (float) – Mean longitude of the ascending node of the Moon (radians).

Return type:

Tuple[float, float, float, float, float]

Notes

Based on IERS Conventions (2010) expressions.

Examples

>>> import numpy as np
>>> T = 0.0  # J2000.0
>>> l, lp, F, D, Om = fundamental_arguments(T)
>>> all(0 <= x < 2*np.pi for x in [l, lp, F, D, Om])  # All in [0, 2pi)
True

pytcl.gravity.tides.moon_position_approximate(mjd)

Compute approximate geocentric Moon position.

Parameters:

mjd (float) – Modified Julian Date.

Returns:

• r (float) – Distance from Earth center (m).

• lat (float) – Geocentric latitude (radians).

• lon (float) – Geocentric longitude (radians).

Return type:

Tuple[float, float, float]

Notes

Low-precision formula adequate for tidal computations.

Examples

>>> r, lat, lon = moon_position_approximate(58000)
>>> 350000e3 < r < 410000e3  # Distance in km range
True
>>> import numpy as np
>>> -np.pi/2 <= lat <= np.pi/2  # Valid latitude
True

pytcl.gravity.tides.sun_position_approximate(mjd)

Compute approximate geocentric Sun position.

Parameters:

mjd (float) – Modified Julian Date.

Returns:

• r (float) – Distance from Earth center (m).

• lat (float) – Geocentric latitude (radians).

• lon (float) – Geocentric longitude (radians).

Return type:

Tuple[float, float, float]

Notes

Low-precision formula adequate for tidal computations.

Examples

>>> r, lat, lon = sun_position_approximate(58000)
>>> 1.47e11 < r < 1.52e11  # Distance ~1 AU
True
>>> lat == 0.0  # Sun on ecliptic
True

pytcl.gravity.tides.solid_earth_tide_displacement(lat, lon, mjd, h2=0.6078, l2=0.0847)

Compute solid Earth tide displacement at a station.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• mjd (float) – Modified Julian Date.

• h2 (float, optional) – Love number h2. Default is IERS value.

• l2 (float, optional) – Shida number l2. Default is IERS value.

Returns:

Displacement in radial, north, east directions (meters).

Return type:

TidalDisplacement

Notes

Computes the degree-2 solid Earth tide displacement caused by the Moon and Sun. The displacement follows the tide-generating potential with Love and Shida numbers.

Examples

>>> import numpy as np
>>> # Displacement at 45N, 0E on MJD 58000
>>> disp = solid_earth_tide_displacement(np.radians(45), 0, 58000)
>>> abs(disp.radial) < 0.5  # Radial displacement typically < 50cm
True

pytcl.gravity.tides.solid_earth_tide_gravity(lat, lon, mjd, delta=1.1608)

Compute solid Earth tide effect on gravity.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• mjd (float) – Modified Julian Date.

• delta (float, optional) – Gravimetric factor (1 + h - 3k/2). Default is IERS value.

Returns:

Gravity change in vertical and horizontal components (m/s^2).

Return type:

TidalGravity

Notes

The tidal gravity change is computed using the tide-generating potential and the gravimetric factor delta = 1 + h - 3k/2.

Examples

>>> import numpy as np
>>> grav = solid_earth_tide_gravity(np.radians(45), 0, 58000)
>>> abs(grav.delta_g) < 3e-6  # Typically < 3 microGal
True

pytcl.gravity.tides.ocean_tide_loading_displacement(mjd, amplitude, phase, constituents=('M2', 'S2', 'N2', 'K2', 'K1', 'O1', 'P1', 'Q1'))

Compute ocean tide loading displacement.

Parameters:

• mjd (float) – Modified Julian Date.

• amplitude (NDArray) – Loading amplitudes (3, n_constituents) for radial, north, east. Units: meters.

• phase (NDArray) – Loading phases (3, n_constituents) in radians.

• constituents (Tuple[str, ...], optional) – Tidal constituent names. Default is major constituents.

Returns:

Ocean loading displacement (meters).

Return type:

TidalDisplacement

Notes

Ocean tide loading parameters are typically obtained from services like IERS or computed from ocean tide models (e.g., FES2014, GOT4.10).

The displacement is computed as: sum_i amplitude_i * cos(omega_i * t + chi_i - phase_i)

where omega_i is the constituent frequency and chi_i is the astronomical argument.

Examples

>>> import numpy as np
>>> # Example loading coefficients (typically from IERS)
>>> amp = np.array([[0.01, 0.005], [0.002, 0.001], [0.002, 0.001]])
>>> phase = np.array([[0.1, 0.2], [0.3, 0.4], [0.5, 0.6]])
>>> disp = ocean_tide_loading_displacement(58000, amp, phase, ("M2", "S2"))
>>> isinstance(disp, TidalDisplacement)
True

pytcl.gravity.tides.atmospheric_pressure_loading(lat, lon, pressure, reference_pressure=101325.0, admittance=-0.00035)

Compute atmospheric pressure loading displacement.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• pressure (float) – Surface atmospheric pressure in Pascals.

• reference_pressure (float, optional) – Reference pressure in Pascals. Default is 101325 Pa (1 atm).

• admittance (float, optional) – Pressure admittance in m/Pa. Default is -0.35 mm/hPa.

Returns:

Pressure loading displacement (meters).

Return type:

TidalDisplacement

Notes

Atmospheric pressure loading causes the Earth’s surface to deform under changing atmospheric mass. The effect is primarily radial (vertical) with magnitude approximately -0.35 mm per hPa of pressure above the mean.

For precise applications, use Green’s functions from models like ECMWF or NCEP pressure fields.

Examples

>>> import numpy as np
>>> # High pressure (1020 hPa) above reference
>>> disp = atmospheric_pressure_loading(np.radians(45), 0, 102000)
>>> disp.radial < 0  # Surface depressed under high pressure
True

pytcl.gravity.tides.pole_tide_displacement(lat, lon, xp, yp, xp_mean=0.0, yp_mean=0.0, h2=0.6078, l2=0.0847)

Compute pole tide (polar motion) displacement.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• xp (float) – Pole x coordinate in arcseconds.

• yp (float) – Pole y coordinate in arcseconds.

• xp_mean (float, optional) – Mean pole x coordinate in arcseconds. Default 0.

• yp_mean (float, optional) – Mean pole y coordinate in arcseconds. Default 0.

• h2 (float, optional) – Love number. Default is IERS value.

• l2 (float, optional) – Shida number. Default is IERS value.

Returns:

Pole tide displacement (meters).

Return type:

TidalDisplacement

Notes

The pole tide results from the centrifugal effect of polar motion. The Earth’s rotation axis wobbles with a primary period of ~433 days (Chandler wobble) causing small displacements.

Examples

>>> import numpy as np
>>> # Pole offset of 0.1 arcsec
>>> disp = pole_tide_displacement(np.radians(45), 0, 0.1, 0.1)
>>> abs(disp.radial) < 0.03  # Typically < 3cm
True

pytcl.gravity.tides.total_tidal_displacement(lat, lon, mjd, ocean_loading=None, pressure=None, xp=0.0, yp=0.0)

Compute total tidal displacement from all sources.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• mjd (float) – Modified Julian Date.

• ocean_loading (OceanTideLoading, optional) – Ocean loading parameters. If None, ocean loading is not included.

• pressure (float, optional) – Surface pressure in Pascals. If None, pressure loading is not included.

• xp (float, optional) – Pole x coordinate in arcseconds. Default 0.

• yp (float, optional) – Pole y coordinate in arcseconds. Default 0.

Returns:

Total tidal displacement (meters).

Return type:

TidalDisplacement

Examples

>>> import numpy as np
>>> disp = total_tidal_displacement(np.radians(45), 0, 58000)
>>> isinstance(disp, TidalDisplacement)
True

pytcl.gravity.tides.tidal_gravity_correction(lat, lon, mjd)

Compute tidal correction to gravity observations.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• mjd (float) – Modified Julian Date.

Returns:

delta_g – Gravity correction in m/s^2 (add to observed gravity).

Return type:

float

Notes

This correction should be added to absolute gravity observations to remove the tidal signal.

Examples

>>> import numpy as np
>>> corr = tidal_gravity_correction(np.radians(45), 0, 58000)
>>> abs(corr) < 3e-6  # Typically < 300 microGal
True