Tutorial 04: EGM2008 & Mass Gravity Models

This tutorial demonstrates the high-fidelity gravitational mechanics modeled in OpenGNC, covering Two-Body, \(J_2\), and full Spherical Harmonic perturbations (EGM2008).

Theory Prerequisites

1. Two-Body Gravity

The simplest model assuming Earth is a point mass. Acceleration is along the nadir vector:

\[ \mathbf{a}_{2B} = -\frac{\mu}{r^3} \mathbf{r} \]

2. \(J_2\) Oblateness

Accounting for the equatorial bulge of the Earth via the dominant zonal harmonic \(J_2\):

\[ a_{J_2, x} = \frac{3}{2} J_2 \frac{\mu}{r^2} \left(\frac{R_E}{r}\right)^2 \frac{x}{r} \left(5 \frac{z^2}{r^2} - 1\right) \]

Adding this to \(\mathbf{a}_{2B}\) introduces Nodal Precession and Apsidal Rotation.

3. Spherical Harmonics (EGM2008)

For high-accuracy drag, decay, or orbital maintenance triggers, the Earth’s potential \(U\) is expanded into Associated Legendre Polynomials \(P_{nm}\):

\[ U = \frac{\mu}{r} \sum_{n=0}^{N} \left(\frac{R_E}{r}\right)^n \sum_{m=0}^{n} P_{nm}(\sin\phi) \left( C_{nm}\cos(m\lambda) + S_{nm}\sin(m\lambda) \right) \]

import numpy as np
import matplotlib.pyplot as plt
from opengnc.disturbances.gravity import TwoBodyGravity, J2Gravity, HarmonicsGravity
from opengnc.propagators.cowell import CowellPropagator

print("Gravity modules loaded.")

Gravity modules loaded.

1. Comparing Gravity Accelerations

Let’s evaluate the difference in acceleration magnitudes at a specific position in Low Earth Orbit (LEO) for Two-Body, \(J_2\), and a Degree 10 Spherical Harmonic Model.

# Spacecraft Position in ECI [m] (Approx 500km altitude)
r_eci = np.array([6878137.0, 0.0, 0.0]) # Equatorial plane
jd_epoch = 2451545.0 # J2000

# 1. Initialize models
gravity_2b = TwoBodyGravity()
gravity_j2 = J2Gravity()
gravity_harmonics = HarmonicsGravity(n_max=10, m_max=10) # Using N=10 for speed in validation

# 2. Calculate Accelerations
a_2b = gravity_2b.get_acceleration(r_eci)
a_j2 = gravity_j2.get_acceleration(r_eci)
a_har = gravity_harmonics.get_acceleration(r_eci, jd_epoch)

print(f"Two-Body Accel Magnitude:  {np.linalg.norm(a_2b):.6f} m/s^2")
print(f"J2 Gravity Accel Magnitude: {np.linalg.norm(a_j2):.6f} m/s^2")
print(f"Degree 10 Accel Magnitude:  {np.linalg.norm(a_har):.6f} m/s^2")

# Delta metrics
delta_j2 = np.linalg.norm(a_j2 - a_2b)
delta_har = np.linalg.norm(a_har - a_j2)
print(f"\nEffect of J2 alone: {delta_j2:.6f} m/s^2")
print(f"Effect of Higher Order (10x10) over J2: {delta_har:.6e} m/s^2")

Two-Body Accel Magnitude:  8.425509 m/s^2
J2 Gravity Accel Magnitude: 8.437274 m/s^2
Degree 10 Accel Magnitude:  0.016811 m/s^2

Effect of J2 alone: 0.011766 m/s^2
Effect of Higher Order (10x10) over J2: 8.420463e+00 m/s^2

2. Dynamic Impact: Orbit Precession

The \(J_2\) effect causes the right ascension of the ascending node (RAAN) to precess over time. Let’s propagate a short orbit under \(J_2\) and see the trajectory.

from opengnc.visualization.orbit import plot_orbit_3d

# Initial Orbit State [m], [m/s]
r_0 = np.array([6878137.0, 0.0, 0.0])
v_0 = np.array([0.0, 7612.0 * np.cos(np.radians(45)), 7612.0 * np.sin(np.radians(45))]) # 45 deg inclination

propagator = CowellPropagator()

# Setup perturbation function for J2
def j2_perturbation(t, r, v):
    # J2Gravity.get_acceleration already returns (TwoBody + J2)
    # propagator integrates (TwoBody + perturbation_acc_fn)
    # So we must return ONLY the perturbative acceleration
    return J2Gravity().get_acceleration(r) - TwoBodyGravity().get_acceleration(r)

# Propagate for 1 Orbit (approx 5400 seconds)
dt_span = 5400
states = []
steps = 540
dt_step = dt_span / steps

curr_r, curr_v = r_0, v_0
for i in range(steps):
    states.append(curr_r.copy())
    # Single-step propagation workaround for demo loops or use integrate() directly
    # For simplicity, we can just use the internal integrator if desired, 
    # but for visualization we want multiple points.
    curr_r, curr_v = propagator.propagate(curr_r, curr_v, dt_step, perturbation_acc_fn=j2_perturbation)

states = np.array(states)
print(f"Propagated over {steps} discrete steps.")

Propagated over 540 discrete steps.

3. Visualization

Using the opengnc.visualization module, we can display our 3D orbit around a spherical Earth model.

# Render the orbit in 3D
fig = plot_orbit_3d(states, title="LEO orbit under J2 Disturbances")
fig.show()