Custom gravity models#
Added in version 7.12.0.
In addition to the builtin EGM 2008 geopotential model, heyoka.py also provides the ability to define custom gravity models.
Theoretical background#
We assume that the gravitational potential generated by a body can be expressed via a spherical harmonics expansion:
Here \(\mu\) is the gravitational parameter of the body, \(a\) is a reference radius, \(\bar{P}_{nm}\) are normalised associated Legendre functions, \(\bar{C}_{nm}\) and \(\bar{S}_{nm}\) are the normalised spherical harmonics coefficients (these are non-dimensional quantities depending on the distribution of mass, \(a\), and the chosen reference frame), and \(\left( r, \phi, \lambda \right)\) are spherical coordinates.
Note
The definitions and conventions adopted here - in particular, the normalisation of the Legendre functions and of the harmonic coefficients - follow §3.2 of [MG00].
Note that while the gravitational potential is most commonly introduced in terms of spherical coordinates, heyoka.py always computes gravitational potentials and accelerations directly in Cartesian coordinates. This avoids singularities near the poles and ensures good numerical conditioning for all orbits.
Example 1: simple central potential#
As a first example, let us see the computation of the acceleration due to a simple central potential via the sh_gravity_acc() function:
import heyoka as hy
# Introduce the Cartesian variables.
x, y, z = hy.make_vars("x", "y", "z")
# Compute the acceleration in (x, y, z) due to a point-like body situated in the origin.
acc = hy.model.sh_gravity_acc(
xyz=[x, y, z], sh_coefficients=[[1.0, 0.0]], mu=1.0, a=1.0
)
In addition to the Cartesian coordinates, \(\mu\) and \(a\), sh_gravity_acc() takes in input the list of normalised spherical harmonics coefficients sh_coefficients. These are expected to be provided as a list of \(\left(\bar{C}, \bar{S} \right)\) pairs, ordered first by degree and then by order. That is:
In a central potential, and with the point-like body assumed to be in the origin of the reference frame, \(\bar{C}_{0,0}=1\) and \(\bar{S}_{0,0}=0\), while all coefficients of degree \(n \geq 1\) are zero. Hence, we pass in sh_coefficients=[[1.0, 0.0]].
We can now proceed to a numerical evaluation of the computed acceleration against the analytical point-mass acceleration:
# Compile the acceleration into a function of the Cartesian coordinates.
acc_cf = hy.cfunc(acc, [x, y, z])
# Evaluate at a test point and compare against the
# analytical point-mass acceleration -mu*r/|r|^3 (here mu=1).
import numpy as np
r = np.array([1.5, -0.7, 0.3])
print("sh_gravity_acc :", acc_cf(r))
print("analytical :", -r / np.linalg.norm(r) ** 3)
sh_gravity_acc : [-0.3150733 0.14703421 -0.06301466]
analytical : [-0.3150733 0.14703421 -0.06301466]
Example 2: the EGM96 model#
In this second example, we show how we can implement the EGM96 geopotential model in heyoka.py via the custom gravity model machinery.
As a first step, let us introduce the list of EGM96 coefficients up to degree \(n=5\) (in order to keep things short):
# EGM96 normalised [C, S] harmonic coefficients up to degree/order 5, in the
# degree-then-order layout expected by sh_gravity_acc()/sh_gravity_pot()
# (contiguous, starting at (n, m) = (0, 0)).
#
# NOTE: the degree-1 terms (1,0) and (1,1) are not present in the EGM96 file
# (they are zero when the origin coincides with the geocentre). They are
# filled here with zeros to keep the layout gap-free.
egm96_CS = np.array(
[
[1.000000000000e00, 0.000000000000e00], # (0, 0)
[0.000000000000e00, 0.000000000000e00], # (1, 0) (filled)
[0.000000000000e00, 0.000000000000e00], # (1, 1) (filled)
[-0.484165371736e-03, 0.000000000000e00], # (2, 0)
[-0.186987635955e-09, 0.119528012031e-08], # (2, 1)
[0.243914352398e-05, -0.140016683654e-05], # (2, 2)
[0.957254173792e-06, 0.000000000000e00], # (3, 0)
[0.202998882184e-05, 0.248513158716e-06], # (3, 1)
[0.904627768605e-06, -0.619025944205e-06], # (3, 2)
[0.721072657057e-06, 0.141435626958e-05], # (3, 3)
[0.539873863789e-06, 0.000000000000e00], # (4, 0)
[-0.536321616971e-06, -0.473440265853e-06], # (4, 1)
[0.350694105785e-06, 0.662671572540e-06], # (4, 2)
[0.990771803829e-06, -0.200928369177e-06], # (4, 3)
[-0.188560802735e-06, 0.308853169333e-06], # (4, 4)
[0.685323475630e-07, 0.000000000000e00], # (5, 0)
[-0.621012128528e-07, -0.944226127525e-07], # (5, 1)
[0.652438297612e-06, -0.323349612668e-06], # (5, 2)
[-0.451955406071e-06, -0.214847190624e-06], # (5, 3)
[-0.295301647654e-06, 0.496658876769e-07], # (5, 4)
[0.174971983203e-06, -0.669384278219e-06], # (5, 5)
]
)
Then, correspondingly, we introduce the \(\mu\) and \(a\) parameters from the EGM96 model (SI units):
mu_egm96 = 0.3986004415e15
a_egm96 = 0.6378136300e07
We can then proceed to formulate the gravitational acceleration according to EGM96:
acc_egm96 = hy.model.sh_gravity_acc(
xyz=[x, y, z], sh_coefficients=egm96_CS, mu=mu_egm96, a=a_egm96
)
Let us evaluate numerically how much the accelerations computed by EGM2008 and EGM96 differ at a sample point in LEO:
# Construct the acceleration according to EGM2008.
#
# NOTE: truncated to n=5 in order to match acc_egm96.
acc_egm2008 = np.array(hy.model.egm2008_acc([x, y, z], n=5, m=5))
# Create a compiled function evaluating the acceleration difference
# between EGM96 and EGM2008.
cf_acc_diff = hy.cfunc(acc_egm2008 - acc_egm96, [x, y, z])
# Evaluate the difference numerically.
print(f"Acceleration difference in LEO: {cf_acc_diff([6700e3, 0.0, 0.0])} m/s**2")
Acceleration difference in LEO: [6.17662972e-08 6.55763275e-09 2.96841722e-08] m/s**2
As expected, the difference between the two models is ~1 part in \(10^8\).
Example 3: symbolic coefficients#
In this last example, we show how the custom gravity models machinery supports providing the spherical harmonics coefficients as symbolic, rather than numerical, quantities.
As a first step, we take the original EGM96 coefficients array and we replace the \(\bar{C}_{2,0}\) coefficient with a runtime parameter:
# Convert the dtype to object.
egm96_CS_sym = egm96_CS.astype(object)
# Replace C20 with par[0].
egm96_CS_sym[3, 0] = hy.par[0]
Then, we can compute the gravitational acceleration according to this new partially-symbolic model:
acc_egm96_sym = hy.model.sh_gravity_acc(
xyz=[x, y, z], sh_coefficients=egm96_CS_sym, mu=mu_egm96, a=a_egm96
)
Let us now build compiled functions for numerical evaluation, and let us confirm that the purely numerical and the partially-symbolic models agree:
# Build compiled functions for the evaluation of the two EGM96 models.
cf_egm96_sym = hy.cfunc(acc_egm96_sym, [x, y, z])
cf_egm96 = hy.cfunc(acc_egm96, [x, y, z])
# Fetch the original C20 numerical value.
C20_egm96 = egm96_CS[3, 0]
print(
f"Acceleration difference: {cf_egm96_sym([6700e3, 0.0, 0.0], pars=[C20_egm96]) - cf_egm96([6700e3, 0.0, 0.0])}"
)
Acceleration difference: [0. 0. 0.]
As we can see, the two models are in agreement. With symbolic coefficients, we can build more flexible models that avoid hard-coding numerical constants at creation time.
A second advantage of symbolic coefficients is that we can compute the derivatives of potential and accelerations with respect to the coefficients themselves, which is important in applications such as gravity field recovery.
As a simple example, let us compute \(\partial \mathbf{a}/\partial \bar C_{2,0}\) for the EGM96 model via the diff_tensors() API:
# Compute the Jacobian of the acceleration wrt C20.
da_dC20 = hy.diff_tensors(acc_egm96_sym, [hy.par[0]]).jacobian
# Create a compiled function for the evaluation of da_dC20.
da_dC20_cf = hy.cfunc(da_dC20.flatten(), [x, y, z])
Now let us evaluate the derivative at the usual point in LEO:
da_dC20_cf([6700e3, 0.0, 0.0])
array([26.98997097, -0. , 0. ])
Since our sample point lies on the equator, the derivative of the acceleration is, as expected, purely radial.