Dawn at Ceres¶

import os
import pathlib
import sys

import numpy as np

import brahe as bh

# No `bh.initialize_eop()` call is needed here: this example never converts
# to or from Earth-relative frames, so it has no dependency on Earth
# orientation data. SPICE kernels are loaded below (de440s for Sun/Jupiter
# positions plus a Ceres SPK from Horizons) so that the third-body and solar
# radiation pressure perturbations resolve around Ceres.

# 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)

# Resolve Ceres in the JPL Small-Body Database to get its NAIF/SPK ID and SI
# physical parameters, then generate and load a targeted SPK from Horizons so
# third-body and SRP perturbations resolve around Ceres. Both the SBDB query
# and the SPK are cached under the brahe cache directory and reused on repeat
# runs, so this hits the network only once per machine.
ceres_obj = bh.datasets.sbdb.SBDBClient().lookup("Ceres")
CERES_NAIF_ID = ceres_obj.naif_id()  # 20000001 (SBDB SPK-ID scheme)
GM_CERES = ceres_obj.gm  # m^3/s^2
R_CERES = ceres_obj.radius  # m
print(f"Ceres: {ceres_obj.full_name}, NAIF ID {CERES_NAIF_ID}")
print(f"  GM = {GM_CERES:.6e} m^3/s^2, radius = {R_CERES / 1e3:.1f} km")

# Load the common SPICE kernels (includes de440s) so Sun/Jupiter positions are
# available, then fetch and load a Ceres SPK covering the propagation span.
bh.load_common_spice_kernels()
spk_t0 = bh.Epoch.from_datetime(2015, 12, 1, 0, 0, 0.0, 0.0, bh.TimeSystem.TDB)
spk_t1 = bh.Epoch.from_datetime(2016, 3, 1, 0, 0, 0.0, 0.0, bh.TimeSystem.TDB)
spk = bh.datasets.horizons.HorizonsClient().get_spk(
    bh.datasets.horizons.HorizonsSPKRequest.for_spkid(CERES_NAIF_ID, spk_t0, spk_t1)
)
spk.load()
print(f"Loaded Ceres SPK: {spk.path}")

# Ceres orientation model (IAU WGCCRE 2015 pole/rotation). SBDB does not
# supply spin-pole or prime-meridian constants, so they are set here:
CERES_POLE_RA = 291.418  # deg, ICRF right ascension of the spin pole
CERES_POLE_DEC = 66.764  # deg, ICRF declination of the spin pole
CERES_W0 = 170.650  # deg, prime meridian angle at J2000 TT
CERES_W_RATE = 952.1532635  # deg/day

# A body-fixed frame from the IAU-style pole/prime-meridian model, registered
# as a custom frame: rotation(epoch) -> ICRF-to-body DCM.
_J2000_TT_MJD = 51544.5


def ceres_rotation(epc):
    d = epc.mjd_as_time_system(bh.TimeSystem.TT) - _J2000_TT_MJD
    alpha = np.radians(CERES_POLE_RA)
    delta = np.radians(CERES_POLE_DEC)
    w = np.radians((CERES_W0 + CERES_W_RATE * d) % 360.0)
    return (
        bh.Rz(w, bh.AngleFormat.RADIANS)
        @ bh.Rx(np.pi / 2 - delta, bh.AngleFormat.RADIANS)
        @ bh.Rz(np.pi / 2 + alpha, bh.AngleFormat.RADIANS)
    )


def ceres_omega(epc=None):
    return np.array([0.0, 0.0, np.radians(CERES_W_RATE) / 86400.0])


# The `1` is an arbitrarily chosen registry key for this custom frame; it
# just needs to be unique among registered custom frames.
bh.register_custom_frame(1, ceres_rotation, ceres_omega)
ceres_fixed = bh.CelestialFrame.BodyFixedCustom(CERES_NAIF_ID, 1)

ceres = bh.CentralBody.Custom(
    "Ceres",
    CERES_NAIF_ID,
    GM_CERES,
    radius=R_CERES,
    omega=ceres_omega(),
    fixed_frame=ceres_fixed,
)

# With the Ceres SPK loaded, enable solar and Jovian third-body perturbations
# and solar radiation pressure. The Sun dominates the third-body signal at
# Ceres; Jupiter is added via the built-in `ThirdBody.JUPITER_BARYCENTER`. The
# occulting body for eclipse geometry is Ceres itself (an SRP shadow cast by
# Earth is meaningless 2.8 AU away), supplied as a Custom occulting body.
SPACECRAFT_MASS = 750.0  # kg (Dawn dry mass, approximate)
SRP_AREA = 20.0  # m^2 (solar-array-dominated cross section, approximate)
CR = 1.3  # reflectivity coefficient

srp = bh.SolarRadiationPressureConfiguration(
    area=bh.ParameterSource.value(SRP_AREA),
    cr=bh.ParameterSource.value(CR),
    eclipse_model=bh.EclipseModel.CONICAL,
    occulting_bodies=[
        bh.OccultingBody.Custom(name="Ceres", naif_id=CERES_NAIF_ID, radius=R_CERES)
    ],
)

force_config = bh.ForceModelConfig.for_body(
    ceres,
    bh.GravityConfiguration.point_mass(),
    srp=srp,
    third_body=[
        bh.ThirdBody.SUN,
        bh.ThirdBody.JUPITER_BARYCENTER,
    ],
    mass=bh.ParameterSource.value(SPACECRAFT_MASS),
)

# Dawn LAMO (Low Altitude Mapping Orbit): ~375 km circular polar orbit.
epoch = bh.Epoch.from_datetime(2016, 1, 1, 0, 0, 0.0, 0.0, bh.TimeSystem.UTC)

oe = np.array(
    [R_CERES + 375e3, 0.001, 90.0, 0.0, 0.0, 0.0]
)  # [a [m], e [-], i [deg], RAAN [deg], argp [deg], M [deg]]

state0 = bh.state_koe_to_inertial_for_body(oe, ceres, bh.AngleFormat.DEGREES)

# Ceres spin pole (ICRF) from the IAU pole constants, used below to confirm the
# orbit geometry (declination / right ascension -> unit pole vector). This is
# the same pole the registered frame carries, evaluated directly here.
_pole_ra = np.radians(CERES_POLE_RA)
_pole_dec = np.radians(CERES_POLE_DEC)
ceres_pole = np.array(
    [
        np.cos(_pole_dec) * np.cos(_pole_ra),
        np.cos(_pole_dec) * np.sin(_pole_ra),
        np.sin(_pole_dec),
    ]
)

print(f"LAMO altitude: {(oe[0] - R_CERES) / 1e3:.1f} km, inclination: {oe[2]:.1f} deg")
print(f"Ceres spin pole (ICRF): {np.array2string(ceres_pole, precision=4)}")

prop = bh.NumericalOrbitPropagator.builder(epoch, state0, force_config).build()
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 (here, Ceres-centered) -- the frame the orbit
# was actually integrated in. state_eci would instead re-center the state onto
# Earth via SPK ephemeris (now possible since the Ceres SPK is loaded), which
# is not the frame we want here.
dt = 120.0
epochs = [epoch + t for t in np.arange(0.0, duration, dt)]
radii_km = np.array([np.linalg.norm(prop.state_bci(epc)[:3]) for epc in epochs]) / 1e3
print(f"\nMin radius: {radii_km.min():.2f} km, max radius: {radii_km.max():.2f} km")
print(f"(Ceres radius: {R_CERES / 1e3:.1f} km)")

# Body-fixed position through the registered frame:
x_fixed = prop.state_in_frame(ceres_fixed, epoch + duration)
print(
    f"\nBody-fixed state at t+{duration / 86400.0:.0f} d: {np.array2string(x_fixed, precision=1)}"
)

fig_3d = bh.plot_trajectory_3d(
    [{"trajectory": prop.trajectory, "color": "cyan", "label": "Dawn (LAMO)"}],
    central_body="ceres",
    backend="plotly",
)

# Quantify the perturbations by propagating a two-body (point-mass, no
# third-body, no SRP) baseline from the same initial state and comparing. The
# solar/Jovian third-body and SRP accelerations are small at a ~375 km, ~5.4 h
# LAMO orbit, so over 2 days (~9 revolutions) they perturb the trajectory by
# tens of meters rather than reshaping it -- but the divergence is measurable
# and grows with time, which a pure two-body run cannot reproduce.
baseline_config = bh.ForceModelConfig.for_body(
    ceres, bh.GravityConfiguration.point_mass()
)
baseline = bh.NumericalOrbitPropagator.builder(epoch, state0, baseline_config).build()
baseline.propagate_to(epoch + duration)

pos_divergence_m = np.array(
    [
        np.linalg.norm(prop.state_bci(epc)[:3] - baseline.state_bci(epc)[:3])
        for epc in epochs
    ]
)
vel_divergence_mps = np.array(
    [
        np.linalg.norm(prop.state_bci(epc)[3:] - baseline.state_bci(epc)[3:])
        for epc in epochs
    ]
)
max_divergence_m = pos_divergence_m.max()
print(
    f"\nMax perturbed-vs-two-body position divergence: {max_divergence_m:.2f} m "
    f"(0 for a two-body-only run)"
)

# Perturbation signature: the perturbed run's position and velocity divergence
# from the two-body baseline, growing over the propagation. A pure two-body run
# would sit at zero; the growth is the solar/Jovian third-body and SRP signal.
import plotly.graph_objects as go

hours = np.arange(0.0, duration, dt) / 3600.0
fig_div = go.Figure()
fig_div.add_trace(
    go.Scatter(x=hours, y=pos_divergence_m, name="Position (m)", mode="lines")
)
fig_div.add_trace(
    go.Scatter(
        x=hours,
        y=vel_divergence_mps * 1e3,
        name="Velocity (mm/s)",
        mode="lines",
        yaxis="y2",
    )
)
fig_div.update_layout(
    title="Perturbed vs two-body divergence",
    xaxis_title="Time since epoch (hours)",
    yaxis={"title": "Position divergence (m)"},
    yaxis2={"title": "Velocity divergence (mm/s)", "overlaying": "y", "side": "right"},
    legend={"x": 0.02, "y": 0.98},
)

# Validation
# Inclination relative to the Ceres spin pole confirms the orbit plane was
# built about the Ceres equator, not the ICRF pole. The perturbations barely
# tilt the plane over 2 days (a few 1e-5 deg), so the orbit stays polar; the
# measurable signature of the perturbations is the position divergence above.
incs_deg = []
for epc in epochs:
    x = prop.state_bci(epc)
    h = np.cross(x[:3], x[3:])
    h /= np.linalg.norm(h)
    incs_deg.append(np.degrees(np.arccos(np.clip(np.dot(h, ceres_pole), -1.0, 1.0))))
incs_deg = np.array(incs_deg)
max_excursion = np.max(np.abs(incs_deg - 90.0))
print(f"Inclination rel. Ceres pole: {incs_deg.min():.6f} to {incs_deg.max():.6f} deg")
print(f"Max excursion from 90 deg polar design: {max_excursion:.6f} deg")

# The perturbations move the trajectory measurably off the two-body baseline
# (they are active) while staying small over 2 days: the orbit remains polar to
# well under a hundredth of a degree and inside the LAMO altitude band, and the
# body-fixed state is finite.
print(
    f"\nPerturbations active: {max_divergence_m:.1f} m divergence from the "
    f"two-body baseline; orbit stays polar (excursion {max_excursion:.1e} deg) "
    f"and bounded ({radii_km.min():.1f}-{radii_km.max():.1f} km radius)."
)