SphericalHarmonicsCore

class SphericalHarmonicsCore(order: int, r_ref: float)

Bases: object

Body-fixed spherical-harmonics kernels (normalized coefficients only).

Computational core of SphericalHarmonicsGravity: the normalized Cunningham B_nm recurrence, the acceleration sum, the position Jacobian, and the partials with respect to the C_nm / S_nm coefficients. SPICE rotations, ephemerides, and estimated-parameter handling stay on the Python side, which passes in the body-fixed position and the effective mu / C_nm / S_nm for each evaluation.

With \(R\) = r_ref and fully normalized coefficients, the potential is [1]

\[U = \frac{\mu}{r} \sum_{n=0}^{N} \sum_{m=0}^{n} \left(\frac{R}{r}\right)^{n} \bar{P}_{nm}(\sin\phi) \left[\bar{C}_{nm} \cos m\lambda + \bar{S}_{nm} \sin m\lambda\right] = \frac{\mu}{R} \sum_{n=0}^{N} \sum_{m=0}^{n} \left(\bar{C}_{nm} V_{nm} + \bar{S}_{nm} W_{nm}\right),\]

in terms of the solid harmonics

\[B_{nm} = V_{nm} + i\,W_{nm} = \left(\frac{R}{r}\right)^{n+1} \bar{P}_{nm}(\sin\phi)\, e^{i m \lambda}.\]

These are built from \(B_{00} = R/r\) with Cartesian recurrences [2], so there are no angles and no singularity at the poles:

\[\begin{split}B_{mm} &= \sqrt{\frac{(1+\delta_{1m})(2m+1)}{2m}}\, \frac{R}{r}\,\frac{x + i y}{r}\, B_{m-1,m-1}, \\ B_{nm} &= \sqrt{\frac{4n^2-1}{n^2-m^2}}\, \frac{R}{r}\,\frac{z}{r}\, B_{n-1,m} - \sqrt{\frac{(2n+1)\left[(n-1)^2-m^2\right]}{(2n-3)(n^2-m^2)}} \left(\frac{R}{r}\right)^2 B_{n-2,m}.\end{split}\]

The body-fixed acceleration is [3]

\[\begin{split}\ddot{x} &= \frac{\mu}{R^2} \sum_{n,m} \left[ -g_1 \left(\bar{C}_{nm} V_{n+1,m+1} + \bar{S}_{nm} W_{n+1,m+1}\right) + g_2 \left(\bar{C}_{nm} V_{n+1,m-1} + \bar{S}_{nm} W_{n+1,m-1}\right)\right], \\ \ddot{y} &= \frac{\mu}{R^2} \sum_{n,m} \left[ g_1 \left(-\bar{C}_{nm} W_{n+1,m+1} + \bar{S}_{nm} V_{n+1,m+1}\right) + g_2 \left(-\bar{C}_{nm} W_{n+1,m-1} + \bar{S}_{nm} V_{n+1,m-1}\right)\right], \\ \ddot{z} &= \frac{\mu}{R^2} \sum_{n,m} g_3 \left(\bar{C}_{nm} V_{n+1,m} + \bar{S}_{nm} W_{n+1,m}\right),\end{split}\]

where the \(g_2\) terms are absent for \(m = 0\) and

\[\begin{split}g_1 &= \frac{1+\delta_{0m}}{2} \sqrt{\frac{(2-\delta_{0m})(2n+1)(n+m+2)(n+m+1)}{2(2n+3)}}, \\ g_2 &= \frac{1}{2} \sqrt{\frac{(2-\delta_{0m})(2n+1)(n-m+2)(n-m+1)}{(2-\delta_{1m})(2n+3)}}, \\ g_3 &= -\sqrt{\frac{(2n+1)(n+m+1)(n-m+1)}{2n+3}}.\end{split}\]

The acceleration is linear in \(\bar{C}_{nm}\) and \(\bar{S}_{nm}\), so their partials are the bracketed factors. The position Jacobian differentiates once more, which carries the recurrence to degree \(N + 2\). The _rf methods rotate the body-fixed results with the body-fixed to reference rotation \(\mathbf{R}\): \(\mathbf{a}_{rf} = \mathbf{R}\,\mathbf{a}\) and \(\mathbf{J}_{rf} = \mathbf{R}\,\mathbf{J}\,\mathbf{R}^T\).

Parameters:

orderint

Maximum degree/order of the expansion.

r_reffloat

Reference radius of the body, in km.

Methods

compute_accel_and_jacobian_bf(r_bf, mu, cnm, snm)	Fused acceleration + da/dr on a single B_nm pass, body-fixed frame.
compute_accel_and_jacobian_rf(r_bf, r_b2rf, ...)	Fused acceleration + da/dr, both rotated into the reference frame (R @ a and R @ J @ R^T).
compute_acceleration_bf(r_bf, mu, cnm, snm)	Acceleration in the body-fixed frame, in km/s^2.
compute_acceleration_rf(r_bf, r_b2rf, mu, ...)	Acceleration rotated into the reference frame (R @ a), km/s^2.
compute_bnm(r_bf)	Normalized exterior B_nm matrices for a body-fixed position [km].
compute_jacobian_ad_bf(r_bf, mu, cnm, snm)	da/dr in the body-fixed frame via forward-mode automatic differentiation (dual numbers) of the acceleration.
compute_jacobian_bf(r_bf, mu, cnm, snm)	da/dr in the body-fixed frame, in 1/s^2.
compute_jacobian_rf(r_bf, r_b2rf, mu, cnm, snm)	da/dr rotated into the reference frame (R @ J @ R^T), 1/s^2.
compute_partial_by_c_bf(r_bf, mu)	da/dC_nm in the body-fixed frame, in km/s^2.
compute_partial_by_cs_ad_bf(r_bf, mu, cnm, snm)	da/dC_nm and da/dS_nm in the body-fixed frame via forward-mode automatic differentiation (dual numbers), in km/s^2.
compute_partial_by_mu_ad_rf(r_bf, r_b2rf, ...)	da/dGM in the reference frame via forward-mode automatic differentiation (dual numbers), in km/s^2 per km^3/s^2.
compute_partial_by_pole_ad_rf(r_rel_rf, ...)	Pole partials in the reference frame via forward-mode automatic differentiation (dual numbers) of the rotated acceleration, the IAU pole rotation m_j2rf @ Rz(-(90+RA)) @ Rx(-(90-DEC)) @ Rz(-W) included, with W = w0 + pm_rate * d_days [deg].
compute_partial_by_s_bf(r_bf, mu)	da/dS_nm in the body-fixed frame, in km/s^2.

References

[1]

W. M. Kaula, Theory of Satellite Geodesy. Waltham, MA: Blaisdell, 1966.

[2]

L. E. Cunningham, “On the computation of the spherical harmonic terms needed during the numerical integration of the orbital motion of an artificial satellite,” Celestial Mechanics, vol. 2, no. 2, pp. 207-216, 1970.

[3]

O. Montenbruck and E. Gill, Satellite Orbits: Models, Methods and Applications. Berlin: Springer, 2000, Sec. 3.2.

compute_accel_and_jacobian_bf(r_bf: list, mu: float, cnm: ndarray, snm: ndarray, mu_jac: float | None = None, cnm_jac: TypeAliasForwardRef('numpy.ndarray') | None = None, snm_jac: TypeAliasForwardRef('numpy.ndarray') | None = None) → tuple[List[float], List[float], TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

Fused acceleration + da/dr on a single B_nm pass, body-fixed frame.

The jacobian uses mu_jac/cnm_jac/snm_jac when given (nominal coefficients while the acceleration uses estimated ones), otherwise the acceleration’s values.

Returns:

(acc, jac, bnm_real, bnm_imag)tuple

acc as [ax, ay, az]; jac flattened row-major.

compute_accel_and_jacobian_rf(r_bf: list, r_b2rf: ndarray, mu: float, cnm: ndarray, snm: ndarray, mu_jac: float | None = None, cnm_jac: TypeAliasForwardRef('numpy.ndarray') | None = None, snm_jac: TypeAliasForwardRef('numpy.ndarray') | None = None) → tuple[List[float], List[float], TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

Fused acceleration + da/dr, both rotated into the reference frame (R @ a and R @ J @ R^T). Returns (acc_rf, jac_rf, bnm_real, bnm_imag).

compute_acceleration_bf(r_bf: list, mu: float, cnm: ndarray, snm: ndarray) → tuple[List[float], TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

Acceleration in the body-fixed frame, in km/s^2.

Returns:

(acc, bnm_real, bnm_imag)tuple

acc as [ax, ay, az] plus the B_nm matrices.

compute_acceleration_rf(r_bf: list, r_b2rf: ndarray, mu: float, cnm: ndarray, snm: ndarray) → tuple[List[float], TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

Acceleration rotated into the reference frame (R @ a), km/s^2.

r_b2rf is the body-fixed -> reference-frame DCM (3x3). Returns (acc_rf, bnm_real, bnm_imag).

compute_bnm(r_bf: list) → tuple[TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

Normalized exterior B_nm matrices for a body-fixed position [km].

Returns:

(bnm_real, bnm_imag)tuple of np.ndarray

Each of shape (order+3, order+3).

compute_jacobian_ad_bf(r_bf: list, mu: float, cnm: ndarray, snm: ndarray) → List[float]

da/dr in the body-fixed frame via forward-mode automatic differentiation (dual numbers) of the acceleration.

Independent of compute_jacobian_bf; the two agree to machine precision, which verifies the analytical position-partial formulas. Returns the flattened row-major 3x3 Jacobian.

compute_jacobian_bf(r_bf: list, mu: float, cnm: ndarray, snm: ndarray) → tuple[List[float], TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

da/dr in the body-fixed frame, in 1/s^2.

Returns:

(jac, bnm_real, bnm_imag)tuple

jac flattened row-major as [d00, d01, d02, d10, d11, d12, d20, d21, d22].

compute_jacobian_rf(r_bf: list, r_b2rf: ndarray, mu: float, cnm: ndarray, snm: ndarray) → tuple[List[float], TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

da/dr rotated into the reference frame (R @ J @ R^T), 1/s^2. Returns (jac_rf, bnm_real, bnm_imag).

compute_partial_by_c_bf(r_bf: list, mu: float) → tuple[TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

da/dC_nm in the body-fixed frame, in km/s^2.

Returns:

(partials, bnm_real, bnm_imag)tuple of np.ndarray

partials of shape (3, order+1, order+1).

compute_partial_by_cs_ad_bf(r_bf: list, mu: float, cnm: ndarray, snm: ndarray) → tuple[TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

da/dC_nm and da/dS_nm in the body-fixed frame via forward-mode automatic differentiation (dual numbers), in km/s^2.

Independent of compute_partial_by_c_bf / compute_partial_by_s_bf.

Returns:

(d_c, d_s)tuple of np.ndarray

Each of shape (3, order+1, order+1); S_n0 left at zero.

compute_partial_by_mu_ad_rf(r_bf: list, r_b2rf: ndarray, mu: float, cnm: ndarray, snm: ndarray, r_origin_sph: list | None = None) → List[float]

da/dGM in the reference frame via forward-mode automatic differentiation (dual numbers), in km/s^2 per km^3/s^2.

With r_origin_sph (the harmonics body relative to the propagation origin, when they differ) the indirect term -GM r_o / |r_o|^3 is differentiated too. Independent of the analytical a / GM partial.

compute_partial_by_pole_ad_rf(r_rel_rf: list, m_j2rf: ndarray, ra: float, dec: float, w0: float, pm_rate: float, d_days: float, mu: float, cnm: ndarray, snm: ndarray) → tuple[List[float], List[List[float]]]

Pole partials in the reference frame via forward-mode automatic differentiation (dual numbers) of the rotated acceleration, the IAU pole rotation m_j2rf @ Rz(-(90+RA)) @ Rx(-(90-DEC)) @ Rz(-W) included, with W = w0 + pm_rate * d_days [deg].

Returns:

(a_rf, partials)tuple

a_rf in km/s^2 and [da/dRA, da/dDEC, da/dPM], in km/s^2 per rad, per rad, and per deg/day.

compute_partial_by_s_bf(r_bf: list, mu: float) → tuple[TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray'), TypeAliasForwardRef('numpy.ndarray')]

da/dS_nm in the body-fixed frame, in km/s^2.

Returns:

(partials, bnm_real, bnm_imag)tuple of np.ndarray

partials of shape (3, order+1, order+1).