Skip to content

LRO Lunar Orbit

In this example we'll set up an LRO-like low lunar science orbit and propagate it with brahe's lunar force model. The Lunar Reconnaissance Orbiter (LRO) has flown a polar, near-frozen low lunar science orbit around the Moon since 2009, mapping the surface at high resolution; this example uses an LRO-like ~30 x 180 km class orbit rather than reproducing its exact orbital history. We'll compare a full-fidelity propagation against a point-mass Moon to quantify how the Moon's lumpy gravity field ("mascons" - dense, gravitationally anomalous regions beneath several lunar maria) perturbs a low lunar orbit, then visualize the trajectory in 3D around a textured Moon.

Why Lunar Orbits Are Different

Unlike Earth, the Moon's gravity field is dominated by large, irregular mass concentrations rather than a smooth oblateness term. Orbit designers counter the resulting perturbations by choosing a "frozen orbit" - an inclination and argument of perilune where the long-period perturbation from the gravity field's dominant harmonics averages to zero, so eccentricity and perilune altitude stay bounded over many orbits instead of drifting monotonically. Even so, these mascon-driven gravity anomalies are strong enough that most low lunar orbits are unstable on timescales of weeks to months without station-keeping. Proper modeling of the high-order gravity dynamics enable design of station-keeping maneuvers that minimize the growth of these instabilities, conserving fuel and extending mission-lifetime.

Orbit Setup

The propagator integrates in the Moon-Centered Inertial (LCI) frame, whose axes are ICRF-aligned: the LCI z-axis is the ICRF pole, which sits about 22 degrees from the Moon's spin pole. Passing the elements straight to state_koe_to_eci would measure the 85.2 degree inclination and the south-pole perilune against the ICRF pole, not the lunar equator. state_koe_to_inertial_for_body instead references the elements to the body's mean equator at J2000 - the plane normal to the body's IAU pole, with the x-axis on the ascending node of that equator on the ICRF equator - and returns the state directly in the body-centered inertial frame (LCI for the Moon). It takes a CentralBody (which supplies both the Moon's gravitational parameter and its pole), so the orbit is placed against the lunar equator with no manual basis construction. We also evaluate the lunar mean pole at J2000 (the third row of the ICRF-to-lunar-body-fixed rotation) to confirm the geometry below:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
import os
import pathlib
import sys

import numpy as np
import plotly.graph_objects as go

import brahe as bh

bh.initialize_eop()
bh.load_common_spice_kernels()

With the standard preamble in place, the next step sets up the frozen orbit geometry.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# LRO-like frozen science orbit: ~30 x 180 km polar orbit. Perilune is
# kept over the southern hemisphere (argument of perilune 270 deg).
epoch = bh.Epoch.from_datetime(2024, 3, 1, 0, 0, 0.0, 0.0, bh.TimeSystem.UTC)
params = np.array([1000.0, 0.0, 0.0, 10.0, 1.3])  # mass, -, -, srp_area, Cr

r_p = bh.R_MOON + 30e3
r_a = bh.R_MOON + 180e3
a = (r_p + r_a) / 2
e = (r_a - r_p) / (r_a + r_p)
oe = np.array([a, e, 85.2, 0.0, 270.0, 0.0])  # [m, -, deg, deg, deg, deg]

# state_koe_to_inertial_for_body references the elements to the Moon's mean
# equator at J2000 (the plane normal to the lunar IAU pole) and returns the
# state directly in the LCI frame the propagator integrates in, so the 85.2 deg
# inclination and south-pole perilune are measured against the Moon's equator
# rather than the ICRF pole.
state0 = bh.state_koe_to_inertial_for_body(
    oe, bh.CentralBody.Moon, bh.AngleFormat.DEGREES
)

# Lunar mean pole (ICRF) at J2000, used below to confirm the orbit geometry:
# the third row of the ICRF -> lunar body-fixed (IAU) rotation is the spin pole.
j2000 = bh.Epoch.from_datetime(2000, 1, 1, 12, 0, 0.0, 0.0, bh.TimeSystem.TDB)
moon_pole = np.asarray(bh.rotation_icrf_to_body_fixed_iau(301, j2000))[2, :]

print(
    f"Perilune radius: {r_p / 1e3:.1f} km (altitude: {(r_p - bh.R_MOON) / 1e3:.0f} km)"
)
print(
    f"Apolune radius:  {r_a / 1e3:.1f} km (altitude: {(r_a - bh.R_MOON) / 1e3:.0f} km)"
)
print(f"Semi-major axis: {a / 1e3:.1f} km, eccentricity: {e:.4f}")

Propagation

We propagate the same initial state under two force models: lunar_default() (50x50 GRGM660PRIM gravity, SRP occulted by the Moon and Earth, and Earth/Sun third-body perturbations) and a point-mass Moon via ForceModelConfig.for_body with GravityConfiguration.point_mass(). Both propagators integrate in the Moon-Centered Inertial (LCI) frame:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Full lunar force model: 50x50 GRGM660PRIM gravity, SRP, Earth+Sun third
# bodies (downloads the gravity model and kernels on first run).
prop_full = bh.NumericalOrbitPropagator(
    epoch,
    state0,
    bh.NumericalPropagationConfig.default(),
    bh.ForceModelConfig.lunar_default(),
    params,
)
# Point-mass comparison: same orbit, Moon treated as a point mass.
prop_pm = bh.NumericalOrbitPropagator(
    epoch,
    state0,
    bh.NumericalPropagationConfig.default(),
    bh.ForceModelConfig.for_body(
        bh.CentralBody.Moon, bh.GravityConfiguration.point_mass()
    ),
    params,
)
duration = 7 * 86400.0
print(f"\nPropagating {duration / 86400.0:.0f} days (full model vs. point mass)...")
prop_full.propagate_to(epoch + duration)
prop_pm.propagate_to(epoch + duration)
print("  Complete!")

state_bci vs. state_eci

state_bci (API) returns the propagator's native state in the central body's body-centered inertial frame (LCI for a Moon-centered propagator). state_eci (API) instead always returns an Earth-centered state - for a Moon-centered propagator it adds the Moon's Earth-relative position, which would report altitudes near the Earth-Moon distance rather than above the lunar surface. Use state_bci whenever you need the distance from the body being orbited.

Altitude Comparison

Sampling both trajectories over the 7-day propagation and computing altitude above R_MOON from the Moon-centered radius, then plotting the difference between the full-gravity and point-mass solutions, shows how steadily the mascons perturb the orbit:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
# state_bci returns the state in the propagator's native Moon-centered
# inertial (LCI) frame. state_eci would instead add the Moon's Earth-relative
# position, which is not what we want when measuring altitude above the
# lunar surface.
dt = 120.0
epochs = [epoch + t for t in np.arange(0.0, duration, dt)]
r_full = np.array([np.linalg.norm(prop_full.state_bci(e)[:3]) for e in epochs])
r_pm = np.array([np.linalg.norm(prop_pm.state_bci(e)[:3]) for e in epochs])
alt_full_km = (r_full - bh.R_MOON) / 1e3
alt_pm_km = (r_pm - bh.R_MOON) / 1e3
times_days = np.arange(0.0, duration, dt) / 86400.0

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# Full minus point-mass altitude: overlaying the two raw traces over 7 days
# (~90 orbits) is an unreadable smear, so plot their difference instead.
diff_km = alt_full_km - alt_pm_km

fig_altitude = go.Figure()

fig_altitude.add_trace(
    go.Scatter(
        x=times_days.tolist(),
        y=diff_km.tolist(),
        mode="lines",
        line={"color": "red", "width": 2},
        name="Full − point mass",
    )
)

fig_altitude.update_layout(
    title="Modeled Altitude Difference: 50x50 Gravity Field minus Point Mass",
    xaxis_title="Time (days)",
    yaxis_title="Altitude difference (km)",
    height=500,
    margin={"l": 60, "r": 40, "t": 60, "b": 60},
)

Over 7 days, the difference between the full and point-mass solutions grows to about 17 km at its peak - a substantial fraction of the orbit's own 150 km altitude range between perilune and apolune - while the full-model orbit's perilune altitude stays above 26 km, remaining bound to the Moon.

3D Visualization

plot_trajectory_3d accepts central_body="moon" to render an interactive 3D view of the trajectory around a textured Moon. Non-Earth central bodies plot the trajectory in that body's centered-inertial frame, converting through the reference frame router with to_frame() when the trajectory is declared in another frame; a Moon-centered NumericalOrbitPropagator's .trajectory is already in CelestialFrame.LCI, so no conversion is needed. We plot the final 12 hours of the full-gravity trajectory:

1
2
3
4
5
6
fig_3d = bh.plot_trajectory_3d(
    [{"trajectory": prop_full.trajectory, "color": "red", "label": "LRO"}],
    time_range=(epoch + (duration - 12 * 3600.0), epoch + duration),
    central_body="moon",
    backend="plotly",
)

Body textures: Solar System Scope, CC BY 4.0.

Full Code Example

Full Code
lro_lunar_orbit.py
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
import os
import pathlib
import sys

import numpy as np
import plotly.graph_objects as go

import brahe as bh

bh.initialize_eop()
bh.load_common_spice_kernels()

# Configuration for output files
SCRIPT_NAME = pathlib.Path(__file__).stem
OUTDIR = pathlib.Path(os.getenv("BRAHE_FIGURE_OUTPUT_DIR", "./docs/figures/"))
os.makedirs(OUTDIR, exist_ok=True)

# LRO-like frozen science orbit: ~30 x 180 km polar orbit. Perilune is
# kept over the southern hemisphere (argument of perilune 270 deg).
epoch = bh.Epoch.from_datetime(2024, 3, 1, 0, 0, 0.0, 0.0, bh.TimeSystem.UTC)
params = np.array([1000.0, 0.0, 0.0, 10.0, 1.3])  # mass, -, -, srp_area, Cr

r_p = bh.R_MOON + 30e3
r_a = bh.R_MOON + 180e3
a = (r_p + r_a) / 2
e = (r_a - r_p) / (r_a + r_p)
oe = np.array([a, e, 85.2, 0.0, 270.0, 0.0])  # [m, -, deg, deg, deg, deg]

# state_koe_to_inertial_for_body references the elements to the Moon's mean
# equator at J2000 (the plane normal to the lunar IAU pole) and returns the
# state directly in the LCI frame the propagator integrates in, so the 85.2 deg
# inclination and south-pole perilune are measured against the Moon's equator
# rather than the ICRF pole.
state0 = bh.state_koe_to_inertial_for_body(
    oe, bh.CentralBody.Moon, bh.AngleFormat.DEGREES
)

# Lunar mean pole (ICRF) at J2000, used below to confirm the orbit geometry:
# the third row of the ICRF -> lunar body-fixed (IAU) rotation is the spin pole.
j2000 = bh.Epoch.from_datetime(2000, 1, 1, 12, 0, 0.0, 0.0, bh.TimeSystem.TDB)
moon_pole = np.asarray(bh.rotation_icrf_to_body_fixed_iau(301, j2000))[2, :]

print(
    f"Perilune radius: {r_p / 1e3:.1f} km (altitude: {(r_p - bh.R_MOON) / 1e3:.0f} km)"
)
print(
    f"Apolune radius:  {r_a / 1e3:.1f} km (altitude: {(r_a - bh.R_MOON) / 1e3:.0f} km)"
)
print(f"Semi-major axis: {a / 1e3:.1f} km, eccentricity: {e:.4f}")

# Full lunar force model: 50x50 GRGM660PRIM gravity, SRP, Earth+Sun third
# bodies (downloads the gravity model and kernels on first run).
prop_full = bh.NumericalOrbitPropagator(
    epoch,
    state0,
    bh.NumericalPropagationConfig.default(),
    bh.ForceModelConfig.lunar_default(),
    params,
)
# Point-mass comparison: same orbit, Moon treated as a point mass.
prop_pm = bh.NumericalOrbitPropagator(
    epoch,
    state0,
    bh.NumericalPropagationConfig.default(),
    bh.ForceModelConfig.for_body(
        bh.CentralBody.Moon, bh.GravityConfiguration.point_mass()
    ),
    params,
)
duration = 7 * 86400.0
print(f"\nPropagating {duration / 86400.0:.0f} days (full model vs. point mass)...")
prop_full.propagate_to(epoch + duration)
prop_pm.propagate_to(epoch + duration)
print("  Complete!")

# state_bci returns the state in the propagator's native Moon-centered
# inertial (LCI) frame. state_eci would instead add the Moon's Earth-relative
# position, which is not what we want when measuring altitude above the
# lunar surface.
dt = 120.0
epochs = [epoch + t for t in np.arange(0.0, duration, dt)]
r_full = np.array([np.linalg.norm(prop_full.state_bci(e)[:3]) for e in epochs])
r_pm = np.array([np.linalg.norm(prop_pm.state_bci(e)[:3]) for e in epochs])
alt_full_km = (r_full - bh.R_MOON) / 1e3
alt_pm_km = (r_pm - bh.R_MOON) / 1e3
times_days = np.arange(0.0, duration, dt) / 86400.0

# Full minus point-mass altitude: overlaying the two raw traces over 7 days
# (~90 orbits) is an unreadable smear, so plot their difference instead.
diff_km = alt_full_km - alt_pm_km

fig_altitude = go.Figure()

fig_altitude.add_trace(
    go.Scatter(
        x=times_days.tolist(),
        y=diff_km.tolist(),
        mode="lines",
        line={"color": "red", "width": 2},
        name="Full − point mass",
    )
)

fig_altitude.update_layout(
    title="Modeled Altitude Difference: 50x50 Gravity Field minus Point Mass",
    xaxis_title="Time (days)",
    yaxis_title="Altitude difference (km)",
    height=500,
    margin={"l": 60, "r": 40, "t": 60, "b": 60},
)

fig_3d = bh.plot_trajectory_3d(
    [{"trajectory": prop_full.trajectory, "color": "red", "label": "LRO"}],
    time_range=(epoch + (duration - 12 * 3600.0), epoch + duration),
    central_body="moon",
    backend="plotly",
)

# Validation
divergence_km = np.abs(alt_full_km - alt_pm_km).max()

# Inclination of the initial state relative to the Moon's spin pole confirms
# the orbit plane was built about the lunar equator, not the ICRF pole.
h0 = np.cross(state0[:3], state0[3:])
h0 /= np.linalg.norm(h0)
inc0_deg = np.degrees(np.arccos(np.clip(np.dot(h0, moon_pole), -1.0, 1.0)))
print(f"\nInitial inclination (rel. Moon pole): {inc0_deg:.4f} deg (target: 85.2 deg)")
print(f"Min perilune altitude (full model): {alt_full_km.min():.2f} km")
print(f"Max altitude divergence (full vs. point mass): {divergence_km:.3f} km")

assert abs(inc0_deg - 85.2) < 0.1, (
    f"Initial inclination not at design value rel. Moon pole: {inc0_deg:.4f} deg"
)
assert alt_full_km.min() > 0, "Orbit impacted the lunar surface"
assert divergence_km > 12.0, (
    f"Full and point-mass solutions did not diverge as expected: {divergence_km:.3f} km"
)

print("\nExample validated successfully!")

See Also