Skip to content

MRO Mars Orbit

In this example we'll set up an MRO-like sun-synchronous science orbit and propagate it with brahe's Mars force model. The Mars Reconnaissance Orbiter (MRO) has flown a sun-synchronous, near-polar science orbit since 2006, imaging the surface at a consistent local solar time on every pass (Zurek & Smrekar, 2007); this example uses an MRO-like ~255 x 320 km, 92.6 degree inclination orbit. We'll track how the osculating orbital elements evolve under the full force model over a two-day propagation, then visualize the trajectory in 3D around a textured Mars.

Why Sun-Synchronous Orbits Work at Mars

A sun-synchronous orbit keeps its orbital plane at a fixed orientation relative to the Sun, so the spacecraft crosses each latitude at the same local solar time on every pass - valuable for consistent lighting in surface imagery. This only works because a planet's oblateness (the \(J_2\) zonal harmonic) causes the orbital plane to precess at a rate

\[ \dot{\Omega} = -\frac{3}{2} n J_2 \left(\frac{R}{p}\right)^2 \cos i \]

where \(n\) is the orbit's mean motion, \(R\) is the planet's radius, \(p = a(1-e^2)\) is the semi-latus rectum, and \(i\) is the inclination. For a near-polar, slightly retrograde inclination, this precession rate can be tuned to exactly match the planet's mean motion around the Sun. Setting \(\dot{\Omega}\) equal to Mars's heliocentric mean motion (about 0.524 deg/day, from its 687-day year) and solving for \(i\) at this example's semi-major axis and eccentricity, with Mars's \(J_2 = 1.96045 \times 10^{-3}\) and \(R = R_{\text{Mars}}\), gives \(i \approx 92.6\) degrees - close to polar, with the small excess over 90 degrees giving just enough retrograde precession to track Mars's slower year. This example's 255 x 320 km, 92.6 degree orbit is designed around this resonance.

Orbit Setup

The propagator integrates in the Mars-Centered Inertial (MCI) frame, whose axes are ICRF-aligned: the MCI z-axis is the ICRF pole, which sits about 37 degrees from Mars's spin pole. Passing 92.6 degrees straight to state_koe_to_eci would measure the inclination against the ICRF pole, not Mars's equator. state_koe_to_inertial_for_body instead references the elements to Mars's mean equator at J2000 - the plane normal to the Mars IAU pole, with the x-axis on the ascending node of that equator on the ICRF equator - and returns the state directly in MCI. It takes a CentralBody, which supplies both Mars's gravitational parameter and its pole, so the orbit is placed against the Mars equator with no manual basis construction.

The ascending node is still derived from the Sun's direction rather than hardcoded: the right ascension of the ascending node is set to the Sun's right ascension in the Mars-equatorial plane plus 45 degrees, placing the node near a 15:00 (mid-afternoon) local solar time. The Sun direction is projected onto the equatorial basis vectors (the ascending node and its quadrature axis) evaluated at J2000, the same reference plane the elements are referenced to:

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

import numpy as np

import brahe as bh

bh.initialize_eop()
bh.load_common_spice_kernels()

With the standard preamble in place, the next step sets up the sun-synchronous 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
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
# MRO-like sun-synchronous science orbit: ~255 x 320 km, i = 92.6 deg,
# node placed for a mid-afternoon local solar time.
epoch = bh.Epoch.from_datetime(2024, 3, 1, 0, 0, 0.0, 0.0, bh.TimeSystem.UTC)
params = np.array([2180.0, 20.0, 2.2, 20.0, 1.3])  # mass, drag_area, Cd, srp_area, Cr

# Mars mean pole (ICRF) at J2000: the reference plane state_koe_to_inertial_for_body
# uses is the Mars equator at J2000, whose pole is the third row of the ICRF ->
# Mars body-fixed (IAU) rotation. x_eq is the ascending node of that equator on
# the ICRF equator (ICRF pole x spin pole) and y_eq completes the triad; they
# span the equatorial plane the RAAN below is measured in.
j2000 = bh.Epoch.from_datetime(2000, 1, 1, 12, 0, 0.0, 0.0, bh.TimeSystem.TDB)
mars_pole = np.asarray(bh.rotation_icrf_to_body_fixed_iau(499, j2000))[2, :]
x_eq = np.cross([0.0, 0.0, 1.0], mars_pole)
x_eq /= np.linalg.norm(x_eq)
y_eq = np.cross(mars_pole, x_eq)

# Sun-synchronous node placement. At the ascending node the local solar time
# is 12:00 + (RAAN - sun_ra) / (15 deg/hr), where sun_ra is the Sun's right
# ascension in the Mars-equatorial basis, so a 15:00 (mid-afternoon) node
# needs RAAN = sun_ra + 45 deg. The Sun direction is taken from the Mars
# barycenter (NAIF 4); the barycenter-to-center offset is negligible here.
sun_mci = bh.spk_state(10, 4, epoch)[:3]
sun_ra = np.degrees(np.arctan2(sun_mci @ y_eq, sun_mci @ x_eq)) % 360.0
raan = (sun_ra + 45.0) % 360.0

r_p = bh.R_MARS + 255e3
r_a = bh.R_MARS + 320e3
a = (r_p + r_a) / 2
e = (r_a - r_p) / (r_a + r_p)
oe = np.array([a, e, 92.6, raan, 270.0, 0.0])  # [m, -, deg, deg, deg, deg]

# state_koe_to_inertial_for_body references the elements to Mars's mean equator
# at J2000 and returns the state directly in the MCI frame the propagator
# integrates in, so the 92.6 deg inclination is measured against Mars's equator.
state0 = bh.state_koe_to_inertial_for_body(
    oe, bh.CentralBody.Mars, bh.AngleFormat.DEGREES
)

print(
    f"Periapsis radius: {r_p / 1e3:.1f} km (altitude: {(r_p - bh.R_MARS) / 1e3:.0f} km)"
)
print(
    f"Apoapsis radius:  {r_a / 1e3:.1f} km (altitude: {(r_a - bh.R_MARS) / 1e3:.0f} km)"
)
print(f"Semi-major axis: {a / 1e3:.1f} km, eccentricity: {e:.4f}")
print(f"Mars spin pole (ICRF): {np.array2string(mars_pole, precision=4)}")
print(f"RAAN for 15:00 LTAN: {raan:.2f} deg (Sun RA {sun_ra:.2f} deg)")

Propagation

We propagate under ForceModelConfig.mars_default(): 50x50 GMM-2B gravity, exponential atmospheric drag, SRP occulted by Mars, and Sun third-body perturbations. The propagator integrates in the Mars-Centered Inertial (MCI) frame:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
# mars_default(): 50x50 GMM-2B gravity, exponential atmospheric drag, SRP
# occulted by Mars, and Sun third-body perturbations (downloads the gravity
# model on first run).
force_config = bh.ForceModelConfig.mars_default()
prop = bh.NumericalOrbitPropagator(
    epoch, state0, bh.NumericalPropagationConfig.default(), force_config, params
)
duration = 2 * 86400.0
print(f"\nPropagating {duration / 86400.0:.0f} days...")
prop.propagate_to(epoch + duration)
print("  Complete!")

state_bci vs. state_eci for Mars orbits

A Mars-centered propagator's state_bci returns the state in the Mars body-centered inertial frame (MCI), while state_eci always returns an Earth-centered state - for a Mars-centered propagator it adds Mars's Earth-relative position, reporting distances near the Earth-Mars range rather than above the Martian surface. Use state_bci whenever you need the distance from the body being orbited.

Element Evolution

Sampling the trajectory over the 2-day propagation and converting each Cartesian state back to Keplerian elements with state_inertial_to_koe_for_body, then plotting the full six-element set with plot_keplerian_trajectory, shows how the orbit evolves under the full force model:

 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
# state_bci returns the propagator's native state in the central body's
# body-centered inertial frame (MCI for a Mars-centered propagator).
# state_eci would instead always return an Earth-centered state, adding
# Mars's Earth-relative position, which is not what we want here.
dt = 120.0
epochs = [epoch + t for t in np.arange(0.0, duration, dt)]
sma_km, ecc, inc_deg, raan_deg, argp_deg, anom_deg = [], [], [], [], [], []
for epc in epochs:
    x = prop.state_bci(epc)
    # state_inertial_to_koe_for_body references the elements to Mars's mean
    # equator at J2000, so koe[2] is already the Mars-pole-relative inclination
    # (no manual re-measurement against the ICRF pole needed).
    koe = bh.state_inertial_to_koe_for_body(
        x, bh.CentralBody.Mars, bh.AngleFormat.DEGREES
    )
    sma_km.append(koe[0] / 1e3)
    ecc.append(koe[1])
    inc_deg.append(koe[2])
    raan_deg.append(koe[3])
    argp_deg.append(koe[4])
    anom_deg.append(koe[5])
sma_km = np.array(sma_km)
ecc = np.array(ecc)
inc_deg = np.array(inc_deg)
raan_deg = np.array(raan_deg)
argp_deg = np.array(argp_deg)
anom_deg = np.array(anom_deg)
times_sec = np.arange(0.0, duration, dt)
alt_p_km = sma_km * (1 - ecc) - bh.R_MARS / 1e3

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
# plot_keplerian_trajectory's raw-array input assumes SI/radians (sma in
# meters, angles in radians) regardless of the display units requested via
# sma_units/angle_units. The time column must be elapsed seconds so the
# plotter's auto hours/minutes axis labeling applies correctly. The
# inclination here is already Mars-pole-relative, since
# state_inertial_to_koe_for_body references elements to the Mars equator.
koe_history = np.column_stack(
    (
        times_sec,
        sma_km * 1e3,
        ecc,
        np.radians(inc_deg),
        np.radians(raan_deg),
        np.radians(argp_deg),
        np.radians(anom_deg),
    )
)
fig_elements = bh.plot_keplerian_trajectory(
    [{"trajectory": koe_history, "label": "MRO"}],
    angle_units="deg",
    sma_units="km",
    backend="plotly",
)

plot_keplerian_trajectory renders all six elements in a 2x3 grid: semi-major axis, eccentricity, and inclination on the top row, RAAN, argument of periapsis, and mean anomaly on the bottom row. Because state_inertial_to_koe_for_body references all six elements to Mars's mean equator at J2000 - the same plane the orbit was designed in - the inclination, RAAN, and argument of periapsis are directly comparable to their design values (92.6 degrees, the sun-synchronous node, and 270 degrees). The inclination stays within about 0.1 degree of 92.6 degrees over this 2-day window: \(J_2\) drives the nodal precession that makes the orbit sun-synchronous but produces no secular change in inclination, so the residual motion is a bounded short-period oscillation rather than decay. Semi-major axis and eccentricity likewise show short-period oscillation without net secular decay over this timespan, and periapsis altitude stays comfortably above the Mars surface throughout.

3D Visualization

plot_trajectory_3d accepts central_body="mars" to render an interactive 3D view of the trajectory around a textured Mars. 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 Mars-centered NumericalOrbitPropagator's .trajectory is already in CelestialFrame.MCI, so no conversion is needed:

1
2
3
4
5
fig_3d = bh.plot_trajectory_3d(
    [{"trajectory": prop.trajectory, "color": "orange", "label": "MRO"}],
    central_body="mars",
    backend="plotly",
)

Body textures: Solar System Scope, CC BY 4.0.

Full Code Example

Full Code
mro_mars_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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
import os
import pathlib
import sys

import numpy as np

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)

# MRO-like sun-synchronous science orbit: ~255 x 320 km, i = 92.6 deg,
# node placed for a mid-afternoon local solar time.
epoch = bh.Epoch.from_datetime(2024, 3, 1, 0, 0, 0.0, 0.0, bh.TimeSystem.UTC)
params = np.array([2180.0, 20.0, 2.2, 20.0, 1.3])  # mass, drag_area, Cd, srp_area, Cr

# Mars mean pole (ICRF) at J2000: the reference plane state_koe_to_inertial_for_body
# uses is the Mars equator at J2000, whose pole is the third row of the ICRF ->
# Mars body-fixed (IAU) rotation. x_eq is the ascending node of that equator on
# the ICRF equator (ICRF pole x spin pole) and y_eq completes the triad; they
# span the equatorial plane the RAAN below is measured in.
j2000 = bh.Epoch.from_datetime(2000, 1, 1, 12, 0, 0.0, 0.0, bh.TimeSystem.TDB)
mars_pole = np.asarray(bh.rotation_icrf_to_body_fixed_iau(499, j2000))[2, :]
x_eq = np.cross([0.0, 0.0, 1.0], mars_pole)
x_eq /= np.linalg.norm(x_eq)
y_eq = np.cross(mars_pole, x_eq)

# Sun-synchronous node placement. At the ascending node the local solar time
# is 12:00 + (RAAN - sun_ra) / (15 deg/hr), where sun_ra is the Sun's right
# ascension in the Mars-equatorial basis, so a 15:00 (mid-afternoon) node
# needs RAAN = sun_ra + 45 deg. The Sun direction is taken from the Mars
# barycenter (NAIF 4); the barycenter-to-center offset is negligible here.
sun_mci = bh.spk_state(10, 4, epoch)[:3]
sun_ra = np.degrees(np.arctan2(sun_mci @ y_eq, sun_mci @ x_eq)) % 360.0
raan = (sun_ra + 45.0) % 360.0

r_p = bh.R_MARS + 255e3
r_a = bh.R_MARS + 320e3
a = (r_p + r_a) / 2
e = (r_a - r_p) / (r_a + r_p)
oe = np.array([a, e, 92.6, raan, 270.0, 0.0])  # [m, -, deg, deg, deg, deg]

# state_koe_to_inertial_for_body references the elements to Mars's mean equator
# at J2000 and returns the state directly in the MCI frame the propagator
# integrates in, so the 92.6 deg inclination is measured against Mars's equator.
state0 = bh.state_koe_to_inertial_for_body(
    oe, bh.CentralBody.Mars, bh.AngleFormat.DEGREES
)

print(
    f"Periapsis radius: {r_p / 1e3:.1f} km (altitude: {(r_p - bh.R_MARS) / 1e3:.0f} km)"
)
print(
    f"Apoapsis radius:  {r_a / 1e3:.1f} km (altitude: {(r_a - bh.R_MARS) / 1e3:.0f} km)"
)
print(f"Semi-major axis: {a / 1e3:.1f} km, eccentricity: {e:.4f}")
print(f"Mars spin pole (ICRF): {np.array2string(mars_pole, precision=4)}")
print(f"RAAN for 15:00 LTAN: {raan:.2f} deg (Sun RA {sun_ra:.2f} deg)")

# mars_default(): 50x50 GMM-2B gravity, exponential atmospheric drag, SRP
# occulted by Mars, and Sun third-body perturbations (downloads the gravity
# model on first run).
force_config = bh.ForceModelConfig.mars_default()
prop = bh.NumericalOrbitPropagator(
    epoch, state0, bh.NumericalPropagationConfig.default(), force_config, params
)
duration = 2 * 86400.0
print(f"\nPropagating {duration / 86400.0:.0f} days...")
prop.propagate_to(epoch + duration)
print("  Complete!")

# state_bci returns the propagator's native state in the central body's
# body-centered inertial frame (MCI for a Mars-centered propagator).
# state_eci would instead always return an Earth-centered state, adding
# Mars's Earth-relative position, which is not what we want here.
dt = 120.0
epochs = [epoch + t for t in np.arange(0.0, duration, dt)]
sma_km, ecc, inc_deg, raan_deg, argp_deg, anom_deg = [], [], [], [], [], []
for epc in epochs:
    x = prop.state_bci(epc)
    # state_inertial_to_koe_for_body references the elements to Mars's mean
    # equator at J2000, so koe[2] is already the Mars-pole-relative inclination
    # (no manual re-measurement against the ICRF pole needed).
    koe = bh.state_inertial_to_koe_for_body(
        x, bh.CentralBody.Mars, bh.AngleFormat.DEGREES
    )
    sma_km.append(koe[0] / 1e3)
    ecc.append(koe[1])
    inc_deg.append(koe[2])
    raan_deg.append(koe[3])
    argp_deg.append(koe[4])
    anom_deg.append(koe[5])
sma_km = np.array(sma_km)
ecc = np.array(ecc)
inc_deg = np.array(inc_deg)
raan_deg = np.array(raan_deg)
argp_deg = np.array(argp_deg)
anom_deg = np.array(anom_deg)
times_sec = np.arange(0.0, duration, dt)
alt_p_km = sma_km * (1 - ecc) - bh.R_MARS / 1e3

# plot_keplerian_trajectory's raw-array input assumes SI/radians (sma in
# meters, angles in radians) regardless of the display units requested via
# sma_units/angle_units. The time column must be elapsed seconds so the
# plotter's auto hours/minutes axis labeling applies correctly. The
# inclination here is already Mars-pole-relative, since
# state_inertial_to_koe_for_body references elements to the Mars equator.
koe_history = np.column_stack(
    (
        times_sec,
        sma_km * 1e3,
        ecc,
        np.radians(inc_deg),
        np.radians(raan_deg),
        np.radians(argp_deg),
        np.radians(anom_deg),
    )
)
fig_elements = bh.plot_keplerian_trajectory(
    [{"trajectory": koe_history, "label": "MRO"}],
    angle_units="deg",
    sma_units="km",
    backend="plotly",
)

fig_3d = bh.plot_trajectory_3d(
    [{"trajectory": prop.trajectory, "color": "orange", "label": "MRO"}],
    central_body="mars",
    backend="plotly",
)

# Validation
print(f"\nMin periapsis altitude: {alt_p_km.min():.2f} km")
print(f"Mean inclination (rel. Mars pole): {inc_deg.mean():.4f} deg (target: 92.6 deg)")

assert alt_p_km.min() > 0, "Orbit descended below the Mars surface"
assert abs(inc_deg[0] - 92.6) < 0.01, (
    f"Initial inclination not at design value rel. Mars pole: {inc_deg[0]:.4f} deg"
)
# Measured against Mars's spin pole, the osculating inclination stays within
# ~0.1 deg of the 92.6 deg design value: J2 drives the nodal precession that
# makes the orbit sun-synchronous but no secular change in inclination, so
# the remaining variation is a bounded short-period oscillation.
assert abs(inc_deg.mean() - 92.6) < 0.1, (
    f"Mean inclination drifted too far from 92.6 deg: {inc_deg.mean():.4f} deg"
)

print("\nExample validated successfully!")

See Also