Tutorial 06: Environmental Perturbations

This tutorial demonstrates how to model and simulate non-gravitational and third-body environmental perturbations in OpenGNC.

Theory Prerequisites

1. Atmospheric Drag

The drag acceleration opposes the spacecraft’s velocity relative to the rotating atmosphere:

\[ \mathbf{a}_{drag} = -\frac{1}{2} \rho v_{rel}^2 \frac{C_D A}{m} \hat{\mathbf{v}}_{rel} \]

where \(\rho\) is local density, \(C_D\) is drag coefficient, \(A\) is area, and \(m\) is mass.

2. Solar Radiation Pressure (SRP)

The force exerted by solar photons. It pushes away from the Sun:

\[ \mathbf{a}_{srp} = - \nu P_{sun} C_R \frac{A}{m} \hat{\mathbf{u}}_{sun} \]

where \(\nu\) is the shadow function (1 in sunlight, 0 in eclipse), \(P_{sun}\) is radiation pressure at distance, and \(C_R\) is reflectivity coefficient.

3. Third-Body Gravity

Gravitational pull from other celestial bodies (Sun, Moon):

\[ \mathbf{a}_{3B} = \sum \mu_k \left( \frac{\mathbf{s}_k - \mathbf{r}}{|\mathbf{s}_k - \mathbf{r}|^3} - \frac{\mathbf{s}_k}{|\mathbf{s}_k|^3} \right) \]

import numpy as np
import matplotlib.pyplot as plt
from opengnc.disturbances.drag import LumpedDrag
from opengnc.environment.density import HarrisPriester
from opengnc.disturbances.srp import Canonball
from opengnc.disturbances.gravity import ThirdBodyGravity, TwoBodyGravity
from opengnc.propagators.cowell import CowellPropagator

print("Environment modules loaded.")

Environment modules loaded.

1. Atmospheric Drag & Density

Let’s evaluate the atmospheric density and resulting drag acceleration at different altitudes and compare them.

density_model = HarrisPriester()
drag_model = LumpedDrag(density_model)

# Spacecraft Properties
mass = 4.0 # 4 kg (3U CubeSat)
area = 0.03 # 0.03 m^2 cross section
cd = 2.2 # Drag coefficient

# State: 400km Altitude
r_400 = np.array([6778137.0, 0.0, 0.0])
v_400 = np.array([0.0, 7670.0, 0.0])
jd = 2451545.0

a_drag = drag_model.get_acceleration(r_400, v_400, jd, mass, area, cd)
print(f"Drag at 400km: Accordance Vector: {a_drag} m/s^2")
print(f"Magnitude: {np.linalg.norm(a_drag):.6e} m/s^2")

Drag at 400km: Accordance Vector: [-0.00000000e+00 -2.74518523e-06 -0.00000000e+00] m/s^2
Magnitude: 2.745185e-06 m/s^2

2. Solar Radiation Pressure (SRP)

Evaluating are drag acceleration vs SRP at equivalent conditions.

srp_model = Canonball()
cr = 1.3 # Reflectivity

a_srp = srp_model.get_acceleration(r_400, jd, mass, area, cr)
print(f"SRP Accel Vector: {a_srp} m/s^2")
print(f"Magnitude: {np.linalg.norm(a_srp):.6e} m/s^2")

print(f"\nRatio (Drag/SRP): {np.linalg.norm(a_drag)/np.linalg.norm(a_srp):.2f}x")
print("At 400km LEO, Atmospheric Drag dominates heavily. SRP becomes relevant at higher altitudes (MEO/GEO).")

SRP Accel Vector: [-8.00476320e-09  4.01254285e-08  1.73962601e-08] m/s^2
Magnitude: 4.446073e-08 m/s^2

Ratio (Drag/SRP): 61.74x
At 400km LEO, Atmospheric Drag dominates heavily. SRP becomes relevant at higher altitudes (MEO/GEO).

3. Third-Body Gravity

Lastly, let’s see the acceleration from the Moon and Sun.

third_body = ThirdBodyGravity()
a_3b = third_body.get_acceleration(r_400, jd)
print(f"Third-Body Accel Vector: {a_3b} m/s^2")
print(f"Magnitude: {np.linalg.norm(a_3b):.6e} m/s^2")

Third-Body Accel Vector: [5.58127159e-08 5.88545936e-07 1.51971346e-07] m/s^2
Magnitude: 6.104070e-07 m/s^2

4. Total Disturbance Propagation

We can aggregate these and feed into a propagator.

propagator = CowellPropagator()

# Combined Perturbations Function
def environmental_perturbations(t, r, v):
    # Note: t is local inside integrator loop relative to integration start
    curr_jd = jd + (t / 86400.0) 
    
    acc_drag = drag_model.get_acceleration(r, v, curr_jd, mass, area, cd)
    acc_srp  = srp_model.get_acceleration(r, curr_jd, mass, area, cr)
    acc_3b   = third_body.get_acceleration(r, curr_jd)
    
    return acc_drag + acc_srp + acc_3b

# Propagate short window
r_f, v_f = propagator.propagate(r_400, v_400, 3600, perturbation_acc_fn=environmental_perturbations, dt_step=60)
print(f"\nPropagated position after 1 hour: {r_f} m")

Propagated position after 1 hour: [-4.06313714e+06 -5.43038765e+06  1.11069242e-01] m