Earth-Moon Free-Return Trajectory¶

In this example we'll design and fly an Earth-Moon free-return trajectory: a path that departs a low Earth parking orbit, coasts out to the Moon, swings around its far side, and lets lunar gravity bend it back onto an Earth-return leg without any dedicated return burn. This is the trajectory class that gave the early Apollo missions their abort safety margin - if the service propulsion system had failed on the way out, the spacecraft would still have looped around the Moon and returned to a survivable re-entry. Apollo 13's abort after its oxygen tank ruptured depended on exactly this property. Artemis I did not fly a strict free return; Artemis II, the first crewed Artemis flight, does.

We'll use the NumericalOrbitPropagator integrated about the Earth-Moon barycenter (the EMBI frame) with a 5x5 Earth spherical-harmonic field plus Moon and Sun third-body perturbations, target the translunar injection (TLI) delta-v with a bisection search, and apply the burn with a TimeEvent callback. A terminal ValueEvent on geodetic altitude stops the flight at the atmospheric entry interface on the way home. Finally we'll visualize the result in the Earth-Moon Rotating (EMR) frame, where the trajectory traces its characteristic figure-8.

What "Free Return" Means¶

A free-return trajectory is a solution of the Earth-Moon-spacecraft three-body problem whose outbound leg is aimed so that the lunar flyby rotates the velocity vector back toward Earth. No propulsion is needed after the initial injection: the Moon's gravity does the work of turning the spacecraft around. The defining property is passive safety - once on the trajectory, the spacecraft returns to Earth's vicinity on its own.

Geometry¶

The departure geometry fixes everything about the transfer except its energy. We depart a 400 km (ISS-like) circular parking orbit from a point near the antipode of the Moon's position at the expected arrival time, burning prograde in the Moon's instantaneous orbital plane so the transfer apogee reaches toward the Moon roughly half an orbit later.

A real mission targets a two-dimensional B-plane at the Moon (a miss distance and an approach angle). Here AIM_OFFSET_DEG is a simplified stand-in for that second dimension: rotating the departure point ahead of the pure antipode about the orbit normal sets up a flyby that swings around the far side of the Moon rather than passing in front of it. We use spk_state to query the Moon's Earth-relative state from the DE440s ephemeris and build an orthonormal departure frame from it. Because the mission is integrated about the Earth-Moon barycenter, the departure state is translated from ECI into the EMBI frame with state_eci_to_emb.

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 departure geometry.

# Free-return geometry: depart a 400 km (ISS-like) parking orbit from a point
# near the Moon's antipode at the expected arrival time, in the Moon's
# instantaneous orbital plane, with a prograde TLI burn. Rotating the departure
# point ahead of the pure antipode by AIM_OFFSET_DEG sets up a flyby that swings
# around the far side of the Moon and bends the trajectory back onto an
# Earth-return leg. A real mission aims a two-dimensional B-plane target (miss
# distance and approach angle); AIM_OFFSET_DEG is a simplified stand-in for that
# second dimension.
epoch = bh.Epoch.from_datetime(2024, 3, 1, 0, 0, 0.0, 0.0, bh.TimeSystem.UTC)
TRANSFER_TIME = 3.1 * 86400.0
AIM_OFFSET_DEG = 10.0
DEPART_ALT = 400e3


def _rodrigues(vec, axis, angle):
    """Rotate ``vec`` about unit ``axis`` by ``angle`` (Rodrigues' formula)."""
    return (
        vec * np.cos(angle)
        + np.cross(axis, vec) * np.sin(angle)
        + axis * np.dot(axis, vec) * (1.0 - np.cos(angle))
    )


x_moon = bh.spk_state(bh.NAIFId.MOON, bh.NAIFId.EARTH, epoch + TRANSFER_TIME)
r_moon, v_moon = x_moon[:3], x_moon[3:]
h_hat = np.cross(r_moon, v_moon)
h_hat /= np.linalg.norm(h_hat)

r0 = bh.R_EARTH + DEPART_ALT
u_antipode = -r_moon / np.linalg.norm(r_moon)  # opposite the arrival point
u_hat = _rodrigues(u_antipode, h_hat, np.radians(AIM_OFFSET_DEG))
t_hat = np.cross(h_hat, u_hat)  # prograde, same sense as the Moon
v_circ = np.sqrt(bh.GM_EARTH / r0)

# Back the parking-orbit state up by a short coast so the TLI burn fires exactly
# at the designed departure point and epoch. The mission is integrated about the
# Earth-Moon barycenter (EMBI frame), so the departure state is translated from
# ECI into EMBI with ``state_eci_to_emb``.
T_PARK = 1800.0  # 30 min of parking-orbit coast before TLI
n_park = np.sqrt(bh.GM_EARTH / r0**3)
start_epoch = epoch - T_PARK
u0 = _rodrigues(u_hat, h_hat, -n_park * T_PARK)
t0 = np.cross(h_hat, u0)
state_park = bh.state_eci_to_emb(start_epoch, np.concatenate([r0 * u0, v_circ * t0]))

Force Model¶

A free return is shaped by three bodies. We integrate about the Earth-Moon barycenter, so the central gravity term is zero - the barycenter has no mass of its own. Earth carries a 5x5 spherical-harmonic field as an attributed third body (evaluated at the spacecraft's Earth-relative position), and the Moon and Sun are point-mass perturbers from the DE440s ephemeris. This is the EMB-centered force-model pattern used for cislunar propagation: the integration state stays barycentric while Earth-fidelity gravity still acts near perigee. The lunar term is what bends the trajectory home; the Sun is a smaller but non-negligible perturbation accumulated over the multi-day flight.

# Integrate about the Earth-Moon barycenter (EMBI): the barycenter has no mass
# of its own, so the central gravity term is zero. Earth carries a 5x5
# spherical-harmonic field as an attributed third body (evaluated at the
# object's Earth-relative position), and the Moon and Sun are point-mass
# perturbers from DE440s. The lunar term is what bends the path home; the Sun is
# a smaller but non-negligible perturbation over the multi-day flight.
force_config = bh.ForceModelConfig.for_body(
    bh.CentralBody.EMB,
    bh.GravityConfiguration.zero(),
    third_body=[
        bh.ThirdBodyConfiguration(
            bh.ThirdBody.EARTH,
            gravity=bh.GravityConfiguration.spherical_harmonic(degree=5, order=5),
        ),
        bh.ThirdBody.MOON,
        bh.ThirdBody.SUN,
    ],
)
force_config.validate()

Targeting the TLI Delta-V¶

This example's purpose is showing the Earth-Moon rotating frame, so we take some shortcuts to generate the trajectory - a fixed departure geometry and a one-dimensional bisection on the TLI delta-v - which in practice should not be done. Real mission design uses Lambert solvers and B-plane targeting.

With the departure geometry fixed, the only free parameter is the TLI delta-v, which sets the transfer's energy. The miss distance at the Moon is a V-shaped function of that delta-v: too little energy and the transfer apogee never reaches lunar distance, too much and the spacecraft races out ahead of the Moon. The free-return solutions live on the ascending branch of the V, where the perilune radius grows with delta-v. We run a coarse scan to reveal the V and locate its minimum, then bracket the target perilune on the ascending branch and refine the delta-v with a bisection.

The same builder that scores a candidate during targeting flies the final mission, so the trajectory the search converges on is exactly the one flown.

# Geodetic altitude above Earth from an EMBI-centered state. The integration
# state is barycentric, so it is translated to ECI before the altitude is
# computed - a plain AltitudeEvent would measure altitude above the barycenter,
# not the Earth. This scalar drives the terminal re-entry event.
def geodetic_altitude(event_epoch, event_state):
    x_eci = bh.state_emb_to_eci(event_epoch, event_state)
    x_ecef = bh.position_eci_to_ecef(event_epoch, x_eci[:3])
    return bh.position_ecef_to_geodetic(x_ecef, bh.AngleFormat.DEGREES)[2]


# Fly a candidate mission: coast the parking orbit up to the design epoch, apply
# the TLI impulsively through a TimeEvent callback, then integrate. The same
# builder is used to score candidates during targeting and to fly the final
# mission, so the trajectory the targeter searches is exactly the one flown.
def fly(dv, duration, terminal=False):
    """Propagate a candidate mission with a ``dv`` TLI at ``epoch``."""

    def tli_callback(event_epoch, event_state):
        # Burn prograde relative to Earth (the parking-orbit velocity), not
        # relative to the EMBI integration frame: translate to ECI, add the
        # delta-v along the Earth-relative velocity, translate back.
        x_eci = bh.state_emb_to_eci(event_epoch, event_state)
        x_eci[3:6] += dv * x_eci[3:6] / np.linalg.norm(x_eci[3:6])
        return (bh.state_eci_to_emb(event_epoch, x_eci), bh.EventAction.CONTINUE)

    prop = bh.NumericalOrbitPropagator.builder(
        start_epoch, state_park, force_config
    ).build()
    prop.add_event_detector(bh.TimeEvent(epoch, "TLI").with_callback(tli_callback))
    if terminal:
        prop.add_event_detector(
            bh.ValueEvent(
                "Re-entry interface",
                geodetic_altitude,
                120e3,
                bh.EventDirection.DECREASING,
            ).set_terminal()
        )
    prop.propagate_to(epoch + duration)
    return prop


# The miss distance at the Moon is a V-shaped function of the TLI delta-v: too
# little energy and the transfer apogee never reaches lunar distance, too much
# and the spacecraft races past ahead of the Moon. The free-return branch is the
# ascending side of the V, where the perilune radius grows with delta-v. A
# coarse scan locates that branch; a bisection then refines the delta-v to a
# target perilune. This scalar search stands in for the Lambert solvers and
# differential-correction targeters a real mission uses.
def min_moon_distance(dv):
    """Propagate a candidate TLI and return the closest lunar approach [m]."""
    prop = fly(dv, 6.0 * 86400.0)
    return min(
        np.linalg.norm(
            prop.state_eci(epoch + t)[:3]
            - bh.spk_state(bh.NAIFId.MOON, bh.NAIFId.EARTH, epoch + t)[:3]
        )
        for t in np.arange(60.0, 6.0 * 86400.0, 600.0)
    )


TARGET_PERILUNE = bh.R_MOON + 2000e3

# Coarse scan over the near-escape delta-v range to reveal the V and locate its
# minimum (the closest reachable approach for this geometry).
dv_grid = np.arange(3.06e3, 3.13e3, 5.0)
perilunes = np.array([min_moon_distance(dv) for dv in dv_grid])
i_min = int(np.argmin(perilunes))

# Bracket the target on the ascending (free-return) branch, then bisect. The
# V-minimum sits below the target; walk up the ascending side until the next
# grid point crosses the target and refine within that interval.
if TARGET_PERILUNE <= perilunes[i_min]:
    raise ValueError("Target perilune is below the closest reachable approach")
i = i_min
while i + 1 < len(dv_grid) and perilunes[i + 1] < TARGET_PERILUNE:
    i += 1
if i + 1 >= len(dv_grid):
    raise ValueError("Target perilune not reached within the scanned delta-v range")
dv_lo, dv_hi = dv_grid[i], dv_grid[i + 1]
for _ in range(40):
    dv_mid = 0.5 * (dv_lo + dv_hi)
    if min_moon_distance(dv_mid) < TARGET_PERILUNE:
        dv_lo = dv_mid
    else:
        dv_hi = dv_mid
dv_tli = 0.5 * (dv_lo + dv_hi)
print(f"TLI delta-v: {dv_tli / 1e3:.4f} km/s")

Why a coarse scan first?

A naive bisection on delta-v cannot converge here: the miss distance is not monotonic, so a bracket chosen blindly may straddle the bottom of the V, and the descending branch reaches the same perilune values without ever returning to Earth. The coarse scan finds the V's minimum so the bisection can be restricted to the ascending branch, where perilune increases monotonically with delta-v and the free-return solution for this geometry lives. The final propagation, terminated by the entry event, is what actually confirms the Earth return.

Flying the Mission¶

We fly the tuned design as it would be flown. The propagator starts in the parking orbit, coasts one short arc, and applies the TLI impulsively through a TimeEvent callback. Because the integration state is barycentric, the burn is applied along the spacecraft's Earth-relative velocity - the state is translated to ECI, the delta-v is added, and it is translated back. The flight is stopped by a terminal ValueEvent on geodetic altitude, triggered on decreasing altitude at 120 km - the atmospheric entry interface where a real capsule would begin re-entry. A plain AltitudeEvent would measure altitude above the barycenter rather than the Earth, so the custom value function translates the barycentric state to ECI before computing altitude. The flight time to that point falls out of the propagation rather than being prescribed.

# Fly the tuned design to completion: the terminal 120 km altitude event stops
# the propagation at the atmospheric entry interface on the return leg. The
# flight time to that point falls out of the propagation rather than being
# prescribed.
MISSION_TIME = 12.0 * 86400.0
prop = fly(dv_tli, MISSION_TIME, terminal=True)

# The terminal re-entry event ends the flight before MISSION_TIME; measure the
# flown time from TLI and sample only the arc the propagator actually flew.
flight_time = prop.current_epoch() - epoch
print(f"Flight time to re-entry: {flight_time / 86400.0:.2f} days")

Event geometry follows the integration center, not Earth

Because the state is integrated about the Earth-Moon barycenter, any detector or maneuver that assumes an Earth-centered state is wrong by the Earth-barycenter offset (thousands of kilometers). A plain AltitudeEvent measures altitude above the barycenter, and an impulsive burn added to the raw integration velocity would be along the barycentric velocity - so both must translate to ECI first. This example's terminal entry cutoff uses a custom ValueEvent that re-centers on Earth before computing geodetic altitude.

Distance History¶

Sampling the recorded trajectory gives the distance from Earth and from the Moon over the whole flight. The Earth distance climbs to nearly lunar distance and returns to the entry interface; the Moon distance dips sharply at the flyby. The final epoch is appended to the sample times so the re-entry point itself is captured.

# Distance from Earth and from the Moon over the whole flight, sampled from the
# trajectory the propagator recorded. A 1 s offset keeps every sample off the
# TLI event epoch, where the state is discontinuous (pre- vs. post-burn); the
# final epoch is appended so the re-entry point itself is captured.
dt = 600.0
sample_times = np.append(np.arange(1.0, flight_time, dt), flight_time)
sample_epochs = [epoch + t for t in sample_times]
times_days = sample_times / 86400.0
earth_dists = np.array([np.linalg.norm(prop.state_eci(e)[:3]) for e in sample_epochs])
moon_dists = np.array(
    [
        np.linalg.norm(
            prop.state_eci(e)[:3] - bh.spk_state(bh.NAIFId.MOON, bh.NAIFId.EARTH, e)[:3]
        )
        for e in sample_epochs
    ]
)

fig_distance = go.Figure()

fig_distance.add_trace(
    go.Scatter(
        x=times_days.tolist(),
        y=(earth_dists / 1e3).tolist(),
        mode="lines",
        line={"color": "steelblue", "width": 2},
        name="Distance from Earth",
    )
)
fig_distance.add_trace(
    go.Scatter(
        x=times_days.tolist(),
        y=(moon_dists / 1e3).tolist(),
        mode="lines",
        line={"color": "gray", "width": 2, "dash": "dash"},
        name="Distance from Moon",
    )
)

fig_distance.update_layout(
    title="Free-Return Trajectory: Distance from Earth and Moon",
    xaxis_title="Time (days)",
    yaxis_title="Distance (km)",
    height=500,
    margin={"l": 60, "r": 40, "t": 60, "b": 60},
)

The Moon-distance curve reaches its minimum - the perilune - a little over three days after departure, and the Earth-distance curve turns over shortly after, marking the moment the lunar flyby has bent the trajectory back toward home. The spacecraft reaches the entry interface after about six and a half days.

Figure-8 in the Rotating Frame¶

In an inertial frame the free-return path is an unremarkable elongated loop. Its structure only becomes visible in the Earth-Moon Rotating (EMR) frame, which co-rotates with the Earth-Moon line so the Moon sits fixed on one axis. In that frame the outbound leg, the far-side lunar swing-by, and the return leg trace the characteristic figure-8 that is the signature of a free-return trajectory - the outbound and return legs cross near Earth. We build two views: a 3D view around the textured Earth and Moon, and a top-down (X-Y) view with the bodies drawn to scale. The top-down view carries direction-of-travel arrows along the path, and the fixed body spheres are placed at the perilune epoch so the swing-by aligns with the Moon.

# In the Earth-Moon Rotating (EMR) frame the Moon is held fixed on the axis, so
# the free-return path traces the characteristic figure-8 that is invisible in
# an inertial frame. Place the fixed body spheres at the perilune epoch so the
# lunar swing-by aligns with the Moon, and sample the trajectory in the EMR
# frame for both a 3D and a top-down (X-Y) view.
i_perilune = int(np.argmin(moon_dists))
reference_epoch = sample_epochs[i_perilune]

emr_states = np.array(
    [prop.state_in_frame(bh.CelestialFrame.EMR, e) for e in sample_epochs]
)
emr_xyz_km = emr_states[:, :3] / 1e3
emr_vel = emr_states[:, 3:6]

# Direction-of-travel arrows at a handful of points evenly spaced along the arc.
arrow_idx = np.linspace(len(sample_epochs) * 0.05, len(sample_epochs) * 0.96, 8)
arrow_idx = arrow_idx.astype(int)


def _emr_body_xy_km(naif_id):
    """Body position in the EMR frame at the perilune epoch [km]."""
    return (
        bh.position_frame_to_frame(
            bh.CelestialFrame.BodyCenteredICRF(naif_id),
            bh.CelestialFrame.EMR,
            reference_epoch,
            np.zeros(3),
        )
        / 1e3
    )


earth_xy = _emr_body_xy_km(bh.NAIFId.EARTH)
moon_xy = _emr_body_xy_km(bh.NAIFId.MOON)

# 3D view: textured Earth and Moon with the trajectory.
fig_emr = bh.plot_earth_moon_rotating_3d(
    [{"trajectory": prop.trajectory, "color": "#fc3d21", "label": "Free return"}],
    backend="plotly",
    reference_epoch=reference_epoch,
    view_elevation=50.0,
    view_azimuth=-120.0,
    view_distance=2.4,
)

# 2D top-down (X-Y) view: the figure-8 with Earth and Moon drawn to scale and
# arrow markers rotated to the local direction of travel. The self-crossing near
# Earth is the signature of the circumlunar free return. Plotly's arrow marker
# points along +Y at angle 0 and rotates clockwise, so the heading is measured
# clockwise from +Y.
heading_deg = np.degrees(np.arctan2(emr_vel[arrow_idx, 0], emr_vel[arrow_idx, 1]))
fig_emr_2d = go.Figure()
fig_emr_2d.add_trace(
    go.Scatter(
        x=emr_xyz_km[:, 0].tolist(),
        y=emr_xyz_km[:, 1].tolist(),
        mode="lines",
        line={"color": "#fc3d21", "width": 2},
        name="Free return",
    )
)
fig_emr_2d.add_trace(
    go.Scatter(
        x=emr_xyz_km[arrow_idx, 0].tolist(),
        y=emr_xyz_km[arrow_idx, 1].tolist(),
        mode="markers",
        marker={
            "symbol": "arrow",
            "size": 13,
            "angle": heading_deg.tolist(),
            "color": "#fc3d21",
        },
        showlegend=False,
    )
)
theta = np.linspace(0.0, 2.0 * np.pi, 120)
for center, radius, color, name in [
    (earth_xy, bh.R_EARTH / 1e3, "#3f7bc2", "Earth"),
    (moon_xy, bh.R_MOON / 1e3, "#9aa0a6", "Moon"),
]:
    fig_emr_2d.add_trace(
        go.Scatter(
            x=(center[0] + radius * np.cos(theta)).tolist(),
            y=(center[1] + radius * np.sin(theta)).tolist(),
            mode="lines",
            fill="toself",
            fillcolor=color,
            line={"color": color, "width": 1},
            name=name,
        )
    )
fig_emr_2d.update_layout(
    title="Free Return in the Earth-Moon Rotating Frame (top-down)",
    xaxis_title="X (km)",
    yaxis_title="Y (km)",
    yaxis={"scaleanchor": "x", "scaleratio": 1},
    height=600,
    margin={"l": 60, "r": 40, "t": 60, "b": 60},
    legend={"x": 0.99, "y": 0.99, "xanchor": "right", "yanchor": "top"},
)

The top-down view makes the figure-8 unmistakable: the outbound leg swings around the far side of the Moon and the return leg crosses it near Earth, closing the loop.

Body textures: Solar System Scope, CC BY 4.0.

Full Code Example¶

Full Code

earth_moon_free_return.py
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)

# Free-return geometry: depart a 400 km (ISS-like) parking orbit from a point
# near the Moon's antipode at the expected arrival time, in the Moon's
# instantaneous orbital plane, with a prograde TLI burn. Rotating the departure
# point ahead of the pure antipode by AIM_OFFSET_DEG sets up a flyby that swings
# around the far side of the Moon and bends the trajectory back onto an
# Earth-return leg. A real mission aims a two-dimensional B-plane target (miss
# distance and approach angle); AIM_OFFSET_DEG is a simplified stand-in for that
# second dimension.
epoch = bh.Epoch.from_datetime(2024, 3, 1, 0, 0, 0.0, 0.0, bh.TimeSystem.UTC)
TRANSFER_TIME = 3.1 * 86400.0
AIM_OFFSET_DEG = 10.0
DEPART_ALT = 400e3


def _rodrigues(vec, axis, angle):
    """Rotate ``vec`` about unit ``axis`` by ``angle`` (Rodrigues' formula)."""
    return (
        vec * np.cos(angle)
        + np.cross(axis, vec) * np.sin(angle)
        + axis * np.dot(axis, vec) * (1.0 - np.cos(angle))
    )


x_moon = bh.spk_state(bh.NAIFId.MOON, bh.NAIFId.EARTH, epoch + TRANSFER_TIME)
r_moon, v_moon = x_moon[:3], x_moon[3:]
h_hat = np.cross(r_moon, v_moon)
h_hat /= np.linalg.norm(h_hat)

r0 = bh.R_EARTH + DEPART_ALT
u_antipode = -r_moon / np.linalg.norm(r_moon)  # opposite the arrival point
u_hat = _rodrigues(u_antipode, h_hat, np.radians(AIM_OFFSET_DEG))
t_hat = np.cross(h_hat, u_hat)  # prograde, same sense as the Moon
v_circ = np.sqrt(bh.GM_EARTH / r0)

# Back the parking-orbit state up by a short coast so the TLI burn fires exactly
# at the designed departure point and epoch. The mission is integrated about the
# Earth-Moon barycenter (EMBI frame), so the departure state is translated from
# ECI into EMBI with ``state_eci_to_emb``.
T_PARK = 1800.0  # 30 min of parking-orbit coast before TLI
n_park = np.sqrt(bh.GM_EARTH / r0**3)
start_epoch = epoch - T_PARK
u0 = _rodrigues(u_hat, h_hat, -n_park * T_PARK)
t0 = np.cross(h_hat, u0)
state_park = bh.state_eci_to_emb(start_epoch, np.concatenate([r0 * u0, v_circ * t0]))

# Integrate about the Earth-Moon barycenter (EMBI): the barycenter has no mass
# of its own, so the central gravity term is zero. Earth carries a 5x5
# spherical-harmonic field as an attributed third body (evaluated at the
# object's Earth-relative position), and the Moon and Sun are point-mass
# perturbers from DE440s. The lunar term is what bends the path home; the Sun is
# a smaller but non-negligible perturbation over the multi-day flight.
force_config = bh.ForceModelConfig.for_body(
    bh.CentralBody.EMB,
    bh.GravityConfiguration.zero(),
    third_body=[
        bh.ThirdBodyConfiguration(
            bh.ThirdBody.EARTH,
            gravity=bh.GravityConfiguration.spherical_harmonic(degree=5, order=5),
        ),
        bh.ThirdBody.MOON,
        bh.ThirdBody.SUN,
    ],
)
force_config.validate()


# Geodetic altitude above Earth from an EMBI-centered state. The integration
# state is barycentric, so it is translated to ECI before the altitude is
# computed - a plain AltitudeEvent would measure altitude above the barycenter,
# not the Earth. This scalar drives the terminal re-entry event.
def geodetic_altitude(event_epoch, event_state):
    x_eci = bh.state_emb_to_eci(event_epoch, event_state)
    x_ecef = bh.position_eci_to_ecef(event_epoch, x_eci[:3])
    return bh.position_ecef_to_geodetic(x_ecef, bh.AngleFormat.DEGREES)[2]


# Fly a candidate mission: coast the parking orbit up to the design epoch, apply
# the TLI impulsively through a TimeEvent callback, then integrate. The same
# builder is used to score candidates during targeting and to fly the final
# mission, so the trajectory the targeter searches is exactly the one flown.
def fly(dv, duration, terminal=False):
    """Propagate a candidate mission with a ``dv`` TLI at ``epoch``."""

    def tli_callback(event_epoch, event_state):
        # Burn prograde relative to Earth (the parking-orbit velocity), not
        # relative to the EMBI integration frame: translate to ECI, add the
        # delta-v along the Earth-relative velocity, translate back.
        x_eci = bh.state_emb_to_eci(event_epoch, event_state)
        x_eci[3:6] += dv * x_eci[3:6] / np.linalg.norm(x_eci[3:6])
        return (bh.state_eci_to_emb(event_epoch, x_eci), bh.EventAction.CONTINUE)

    prop = bh.NumericalOrbitPropagator.builder(
        start_epoch, state_park, force_config
    ).build()
    prop.add_event_detector(bh.TimeEvent(epoch, "TLI").with_callback(tli_callback))
    if terminal:
        prop.add_event_detector(
            bh.ValueEvent(
                "Re-entry interface",
                geodetic_altitude,
                120e3,
                bh.EventDirection.DECREASING,
            ).set_terminal()
        )
    prop.propagate_to(epoch + duration)
    return prop


# The miss distance at the Moon is a V-shaped function of the TLI delta-v: too
# little energy and the transfer apogee never reaches lunar distance, too much
# and the spacecraft races past ahead of the Moon. The free-return branch is the
# ascending side of the V, where the perilune radius grows with delta-v. A
# coarse scan locates that branch; a bisection then refines the delta-v to a
# target perilune. This scalar search stands in for the Lambert solvers and
# differential-correction targeters a real mission uses.
def min_moon_distance(dv):
    """Propagate a candidate TLI and return the closest lunar approach [m]."""
    prop = fly(dv, 6.0 * 86400.0)
    return min(
        np.linalg.norm(
            prop.state_eci(epoch + t)[:3]
            - bh.spk_state(bh.NAIFId.MOON, bh.NAIFId.EARTH, epoch + t)[:3]
        )
        for t in np.arange(60.0, 6.0 * 86400.0, 600.0)
    )


TARGET_PERILUNE = bh.R_MOON + 2000e3

# Coarse scan over the near-escape delta-v range to reveal the V and locate its
# minimum (the closest reachable approach for this geometry).
dv_grid = np.arange(3.06e3, 3.13e3, 5.0)
perilunes = np.array([min_moon_distance(dv) for dv in dv_grid])
i_min = int(np.argmin(perilunes))

# Bracket the target on the ascending (free-return) branch, then bisect. The
# V-minimum sits below the target; walk up the ascending side until the next
# grid point crosses the target and refine within that interval.
if TARGET_PERILUNE <= perilunes[i_min]:
    raise ValueError("Target perilune is below the closest reachable approach")
i = i_min
while i + 1 < len(dv_grid) and perilunes[i + 1] < TARGET_PERILUNE:
    i += 1
if i + 1 >= len(dv_grid):
    raise ValueError("Target perilune not reached within the scanned delta-v range")
dv_lo, dv_hi = dv_grid[i], dv_grid[i + 1]
for _ in range(40):
    dv_mid = 0.5 * (dv_lo + dv_hi)
    if min_moon_distance(dv_mid) < TARGET_PERILUNE:
        dv_lo = dv_mid
    else:
        dv_hi = dv_mid
dv_tli = 0.5 * (dv_lo + dv_hi)
print(f"TLI delta-v: {dv_tli / 1e3:.4f} km/s")

# Fly the tuned design to completion: the terminal 120 km altitude event stops
# the propagation at the atmospheric entry interface on the return leg. The
# flight time to that point falls out of the propagation rather than being
# prescribed.
MISSION_TIME = 12.0 * 86400.0
prop = fly(dv_tli, MISSION_TIME, terminal=True)

# The terminal re-entry event ends the flight before MISSION_TIME; measure the
# flown time from TLI and sample only the arc the propagator actually flew.
flight_time = prop.current_epoch() - epoch
print(f"Flight time to re-entry: {flight_time / 86400.0:.2f} days")

# Distance from Earth and from the Moon over the whole flight, sampled from the
# trajectory the propagator recorded. A 1 s offset keeps every sample off the
# TLI event epoch, where the state is discontinuous (pre- vs. post-burn); the
# final epoch is appended so the re-entry point itself is captured.
dt = 600.0
sample_times = np.append(np.arange(1.0, flight_time, dt), flight_time)
sample_epochs = [epoch + t for t in sample_times]
times_days = sample_times / 86400.0
earth_dists = np.array([np.linalg.norm(prop.state_eci(e)[:3]) for e in sample_epochs])
moon_dists = np.array(
    [
        np.linalg.norm(
            prop.state_eci(e)[:3] - bh.spk_state(bh.NAIFId.MOON, bh.NAIFId.EARTH, e)[:3]
        )
        for e in sample_epochs
    ]
)

fig_distance = go.Figure()

fig_distance.add_trace(
    go.Scatter(
        x=times_days.tolist(),
        y=(earth_dists / 1e3).tolist(),
        mode="lines",
        line={"color": "steelblue", "width": 2},
        name="Distance from Earth",
    )
)
fig_distance.add_trace(
    go.Scatter(
        x=times_days.tolist(),
        y=(moon_dists / 1e3).tolist(),
        mode="lines",
        line={"color": "gray", "width": 2, "dash": "dash"},
        name="Distance from Moon",
    )
)

fig_distance.update_layout(
    title="Free-Return Trajectory: Distance from Earth and Moon",
    xaxis_title="Time (days)",
    yaxis_title="Distance (km)",
    height=500,
    margin={"l": 60, "r": 40, "t": 60, "b": 60},
)

# In the Earth-Moon Rotating (EMR) frame the Moon is held fixed on the axis, so
# the free-return path traces the characteristic figure-8 that is invisible in
# an inertial frame. Place the fixed body spheres at the perilune epoch so the
# lunar swing-by aligns with the Moon, and sample the trajectory in the EMR
# frame for both a 3D and a top-down (X-Y) view.
i_perilune = int(np.argmin(moon_dists))
reference_epoch = sample_epochs[i_perilune]

emr_states = np.array(
    [prop.state_in_frame(bh.CelestialFrame.EMR, e) for e in sample_epochs]
)
emr_xyz_km = emr_states[:, :3] / 1e3
emr_vel = emr_states[:, 3:6]

# Direction-of-travel arrows at a handful of points evenly spaced along the arc.
arrow_idx = np.linspace(len(sample_epochs) * 0.05, len(sample_epochs) * 0.96, 8)
arrow_idx = arrow_idx.astype(int)


def _emr_body_xy_km(naif_id):
    """Body position in the EMR frame at the perilune epoch [km]."""
    return (
        bh.position_frame_to_frame(
            bh.CelestialFrame.BodyCenteredICRF(naif_id),
            bh.CelestialFrame.EMR,
            reference_epoch,
            np.zeros(3),
        )
        / 1e3
    )


earth_xy = _emr_body_xy_km(bh.NAIFId.EARTH)
moon_xy = _emr_body_xy_km(bh.NAIFId.MOON)

# 3D view: textured Earth and Moon with the trajectory.
fig_emr = bh.plot_earth_moon_rotating_3d(
    [{"trajectory": prop.trajectory, "color": "#fc3d21", "label": "Free return"}],
    backend="plotly",
    reference_epoch=reference_epoch,
    view_elevation=50.0,
    view_azimuth=-120.0,
    view_distance=2.4,
)

# 2D top-down (X-Y) view: the figure-8 with Earth and Moon drawn to scale and
# arrow markers rotated to the local direction of travel. The self-crossing near
# Earth is the signature of the circumlunar free return. Plotly's arrow marker
# points along +Y at angle 0 and rotates clockwise, so the heading is measured
# clockwise from +Y.
heading_deg = np.degrees(np.arctan2(emr_vel[arrow_idx, 0], emr_vel[arrow_idx, 1]))
fig_emr_2d = go.Figure()
fig_emr_2d.add_trace(
    go.Scatter(
        x=emr_xyz_km[:, 0].tolist(),
        y=emr_xyz_km[:, 1].tolist(),
        mode="lines",
        line={"color": "#fc3d21", "width": 2},
        name="Free return",
    )
)
fig_emr_2d.add_trace(
    go.Scatter(
        x=emr_xyz_km[arrow_idx, 0].tolist(),
        y=emr_xyz_km[arrow_idx, 1].tolist(),
        mode="markers",
        marker={
            "symbol": "arrow",
            "size": 13,
            "angle": heading_deg.tolist(),
            "color": "#fc3d21",
        },
        showlegend=False,
    )
)
theta = np.linspace(0.0, 2.0 * np.pi, 120)
for center, radius, color, name in [
    (earth_xy, bh.R_EARTH / 1e3, "#3f7bc2", "Earth"),
    (moon_xy, bh.R_MOON / 1e3, "#9aa0a6", "Moon"),
]:
    fig_emr_2d.add_trace(
        go.Scatter(
            x=(center[0] + radius * np.cos(theta)).tolist(),
            y=(center[1] + radius * np.sin(theta)).tolist(),
            mode="lines",
            fill="toself",
            fillcolor=color,
            line={"color": color, "width": 1},
            name=name,
        )
    )
fig_emr_2d.update_layout(
    title="Free Return in the Earth-Moon Rotating Frame (top-down)",
    xaxis_title="X (km)",
    yaxis_title="Y (km)",
    yaxis={"scaleanchor": "x", "scaleratio": 1},
    height=600,
    margin={"l": 60, "r": 40, "t": 60, "b": 60},
    legend={"x": 0.99, "y": 0.99, "xanchor": "right", "yanchor": "top"},
)