REBOUND and the Neptune Search

This tutorial is the hands-on companion to the historical Neptune notebook. It is designed for readers who want to see the numerical pieces directly: how the bundled data is loaded, how a small REBOUND simulation is assembled, and how the Uranus anomaly becomes a search problem for Neptune.

The notebook stays offline-friendly and only uses the data files bundled under data/.

REBOUND Basics

We start with a minimal integrator setup and a few helper paths so this notebook works from either the repository root or the docs/ directory.

from pathlib import Path

import numpy as np
import pandas as pd
import rebound
from scipy.optimize import differential_evolution

cwd = Path.cwd()
# Works from the repository root or the `docs/` directory.
candidate_data_dirs = [cwd / "data", cwd.parent / "data"]
data_dir = next((path for path in candidate_data_dirs if (path / "outer_planet_state_vectors_1781_01_01.csv").exists()), None)

if data_dir is None:
    raise FileNotFoundError(
        "Could not find bundled notebook data in `data/`."
    )

outer_planets = pd.read_csv(data_dir / "outer_planet_state_vectors_1781_01_01.csv")
observations = pd.read_csv(data_dir / "uranus_observations_1781_1846.csv")
observation_times_jd = observations["datetime_jd"].to_numpy(dtype=float)
simulation_epoch_jd = float(observation_times_jd[0])

Build the Outer Solar System

The next step is to put the Sun and the known outer planets into a small simulation. Neptune is still missing on purpose: this baseline gives us the reference system that the search cells will perturb and compare against the observations.

planet_masses_msun = {
    "jupiter": 9.547919e-4,
    "saturn": 2.858859806e-4,
    "uranus": 4.366244e-5,
}
days_per_year = 365.25


def build_simulation(
    neptune_parameters: dict[str, float] | None = None,
) -> rebound.Simulation:
    sim = rebound.Simulation()
    sim.units = ("yr", "AU", "Msun")
    sim.add(m=1.0)

    for row in outer_planets.itertuples(index=False):
        body_name = row.planet.title()
        vx_yr = float(row.vx) * days_per_year
        vy_yr = float(row.vy) * days_per_year
        vz_yr = float(row.vz) * days_per_year

        sim.add(
            m=planet_masses_msun[row.planet],
            x=float(row.x),
            y=float(row.y),
            z=float(row.z),
            vx=vx_yr,
            vy=vy_yr,
            vz=vz_yr,
            hash=body_name,
        )

    if neptune_parameters is not None:
        sim.add(
            m=neptune_parameters["mass_msun"],
            a=neptune_parameters["semi_major_axis_au"],
            e=neptune_parameters["eccentricity"],
            inc=0.0,
            Omega=0.0,
            omega=0.0,
            M=np.deg2rad(neptune_parameters["mean_longitude_deg"]),
            hash="Neptune",
        )

    sim.move_to_com()
    sim.integrator = "whfast"
    sim.dt = 0.02
    return sim


sim = build_simulation()

Compare Uranus Longitudes

With observations loaded, we can compare modeled and reference Uranus longitudes directly. The same residuals feed the objective below and the final fit summary.

observed_longitudes_rad = observations["longitude_rad"].to_numpy(dtype=float)
observed_longitudes_unwrapped_rad = np.unwrap(observed_longitudes_rad)
observed_offsets_arcsec = np.rad2deg(
    observed_longitudes_unwrapped_rad - observed_longitudes_unwrapped_rad[0]
) * 3600.0


def simulate_uranus_longitudes(
    neptune_parameters: dict[str, float] | None,
    times_jd: np.ndarray = observation_times_jd,
) -> np.ndarray:
    times_jd = np.asarray(times_jd, dtype=float)
    if times_jd.ndim != 1 or times_jd.size == 0:
        raise ValueError(
            "times_jd must be a one-dimensional array with at least one epoch"
        )

    sim = build_simulation(neptune_parameters)
    longitudes_rad = np.empty(times_jd.size, dtype=float)

    for index, jd in enumerate(times_jd):
        years = (float(jd) - simulation_epoch_jd) / days_per_year
        sim.integrate(years, exact_finish_time=1)
        sun = sim.particles[0]
        uranus = sim.particles["Uranus"]
        longitudes_rad[index] = np.arctan2(
            uranus.y - sun.y, uranus.x - sun.x
        )

    longitudes_rad = np.unwrap(longitudes_rad)
    return np.rad2deg(longitudes_rad - longitudes_rad[0]) * 3600.0


observed_offsets_arcsec[:5]

array([    0.        , 15739.33475885, 31583.91811035, 47537.10199437,
       63646.38967435])

Search for Neptune

The optimization below uses a reduced budget so the notebook runs quickly, but it is still the same Neptune search workflow. We report the subset objective and the full-data check afterward.

neptune_parameter_names = (
    "mass_msun",
    "semi_major_axis_au",
    "eccentricity",
    "mean_longitude_deg",
)

neptune_guess = {
    "mass_msun": 5.15e-5,
    "semi_major_axis_au": 30.1,
    "eccentricity": 0.01,
    "mean_longitude_deg": 300.0,
}

search_bounds = {
    "mass_msun": (4.5e-5, 6.0e-5),
    "semi_major_axis_au": (28.0, 32.0),
    "eccentricity": (0.0, 0.1),
    "mean_longitude_deg": (0.0, 360.0),
}
parameter_bounds = [search_bounds[name] for name in neptune_parameter_names]
fit_indices = np.linspace(0, observation_times_jd.size - 1, 12, dtype=int)
fit_times_jd = observation_times_jd[fit_indices]
fit_observed_offsets_arcsec = observed_offsets_arcsec[fit_indices]


def vector_to_neptune_parameters(vector: np.ndarray) -> dict[str, float]:
    values = np.asarray(vector, dtype=float)
    if values.shape != (4,):
        raise ValueError(
            "Neptune search vectors must contain exactly four values"
        )

    return {
        name: float(value)
        for name, value in zip(neptune_parameter_names, values, strict=True)
    }


def objective(neptune_vector: np.ndarray) -> float:
    modeled_offsets_arcsec = simulate_uranus_longitudes(
        vector_to_neptune_parameters(neptune_vector),
        times_jd=fit_times_jd,
    )
    residuals_arcsec = modeled_offsets_arcsec - fit_observed_offsets_arcsec
    return float(np.sqrt(np.mean(residuals_arcsec**2)))


nominal_modeled_offsets_arcsec = simulate_uranus_longitudes(neptune_guess)
comparison = pd.DataFrame(
    {
        "datetime_jd": observation_times_jd,
        "reference_offset_arcsec": observed_offsets_arcsec,
        "modeled_offset_arcsec": nominal_modeled_offsets_arcsec,
    }
)
comparison["residual_arcsec"] = (
    comparison["modeled_offset_arcsec"] - comparison["reference_offset_arcsec"]
)
comparison.head()

	datetime_jd	reference_offset_arcsec	modeled_offset_arcsec	residual_arcsec
0	2371557.5	0.000000	0.000000	0.000000
1	2371922.5	15739.334759	15739.487756	0.152997
2	2372287.5	31583.918110	31584.226640	0.308529
3	2372652.5	47537.101994	47537.637604	0.535609
4	2373018.5	63646.389674	63647.236812	0.847138

optimizer_result = differential_evolution(
    objective,
    bounds=parameter_bounds,
    maxiter=2,
    popsize=4,
    seed=7,
    polish=False,
    updating="deferred",
    workers=1,
)

best_neptune_parameters = vector_to_neptune_parameters(optimizer_result.x)
best_fit_offsets_arcsec = simulate_uranus_longitudes(best_neptune_parameters)
fit_subset_rms_arcsec = float(optimizer_result.fun)
best_fit_comparison = pd.DataFrame(
    {
        "datetime_jd": observation_times_jd,
        "reference_offset_arcsec": observed_offsets_arcsec,
        "modeled_offset_arcsec": best_fit_offsets_arcsec,
    }
)
best_fit_comparison["residual_arcsec"] = (
    best_fit_comparison["modeled_offset_arcsec"]
    - best_fit_comparison["reference_offset_arcsec"]
)
best_fit_full_data_rms_arcsec = float(
    np.sqrt(np.mean(best_fit_comparison["residual_arcsec"]**2))
)
best_fit_summary = pd.Series(
    {
        "fit_subset_rms_arcsec": fit_subset_rms_arcsec,
        "best_fit_full_data_rms_arcsec": best_fit_full_data_rms_arcsec,
        "mass_msun": best_neptune_parameters["mass_msun"],
        "semi_major_axis_au": best_neptune_parameters["semi_major_axis_au"],
        "eccentricity": best_neptune_parameters["eccentricity"],
        "mean_longitude_deg": best_neptune_parameters["mean_longitude_deg"],
    }
)
best_fit_summary

fit_subset_rms_arcsec              6.988950
best_fit_full_data_rms_arcsec      6.593544
mass_msun                          0.000056
semi_major_axis_au                31.920728
eccentricity                       0.090777
mean_longitude_deg               208.711104
dtype: float64

Optional: Horizons Fetch

The tutorial stays offline by default, but the next code cell shows how to rebuild the same bundled tables directly from JPL Horizons in a networked environment. If live fetch fails or stays disabled, the notebook continues to use the bundled CSV files.

ENABLE_HORIZONS_FETCH = False

HORIZONS_PLANET_IDS = {
    "jupiter": "599",
    "saturn": "699",
    "uranus": "799",
}
HORIZONS_QUERY = {
    "start": "1781-01-01",
    "stop": "1846-09-23",
    "step": "1y",
}

def query_horizons_vectors(
    horizons_id: str,
    *,
    start: str,
    stop: str,
    step: str,
) -> pd.DataFrame:
    from astroquery.jplhorizons import Horizons

    table = Horizons(
        id=horizons_id,
        location="@sun",
        epochs={"start": start, "stop": stop, "step": step},
    ).vectors()
    if len(table) == 0:
        raise RuntimeError(f"Horizons returned no rows for id={horizons_id}")

    frame = pd.DataFrame(
        {
            "datetime_jd": np.asarray(table["datetime_jd"], dtype=float),
            "x": np.asarray(table["x"], dtype=float),
            "y": np.asarray(table["y"], dtype=float),
            "z": np.asarray(table["z"], dtype=float),
            "vx": np.asarray(table["vx"], dtype=float),
            "vy": np.asarray(table["vy"], dtype=float),
            "vz": np.asarray(table["vz"], dtype=float),
        }
    )
    frame["longitude_rad"] = np.arctan2(frame["y"], frame["x"])
    return frame

def fetch_horizons_observations() -> pd.DataFrame:
    uranus_vectors = query_horizons_vectors(
        HORIZONS_PLANET_IDS["uranus"],
        **HORIZONS_QUERY,
    )
    return uranus_vectors.loc[:, ["datetime_jd", "longitude_rad"]].copy()

def fetch_horizons_state_vectors() -> pd.DataFrame:
    rows: list[dict[str, float | str]] = []
    for planet_name in ("jupiter", "saturn", "uranus"):
        vectors = query_horizons_vectors(
            HORIZONS_PLANET_IDS[planet_name],
            **HORIZONS_QUERY,
        )
        first_row = vectors.iloc[0]
        rows.append(
            {
                "planet": planet_name,
                "x": float(first_row["x"]),
                "y": float(first_row["y"]),
                "z": float(first_row["z"]),
                "vx": float(first_row["vx"]),
                "vy": float(first_row["vy"]),
                "vz": float(first_row["vz"]),
            }
        )
    return pd.DataFrame(rows)

live_observations = observations.copy()
live_outer_planets = outer_planets.copy()
live_fetch_status = (
    "Using bundled data from data/. Set ENABLE_HORIZONS_FETCH = "
    "True to refresh the tutorial inputs from JPL Horizons."
)

if ENABLE_HORIZONS_FETCH:
    try:
        live_observations = fetch_horizons_observations()
        live_outer_planets = fetch_horizons_state_vectors()
        live_fetch_status = (
            "Fetched live Horizons data with the same observation and state-"
            "vector columns as the bundled CSV files."
        )
    except Exception as exc:
        live_observations = observations.copy()
        live_outer_planets = outer_planets.copy()
        live_fetch_status = (
            "Using bundled data after a live Horizons fetch failure: "
            f"{exc}"
        )

print(live_fetch_status)
pd.Series(
    {
        "observation_rows": int(live_observations.shape[0]),
        "state_vector_rows": int(live_outer_planets.shape[0]),
        "observation_columns": ", ".join(live_observations.columns),
        "state_vector_columns": ", ".join(live_outer_planets.columns),
    }
)

Using bundled data from data/. Set ENABLE_HORIZONS_FETCH = True to refresh the tutorial inputs from JPL Horizons.

observation_rows                                 66
state_vector_rows                                 3
observation_columns      datetime_jd, longitude_rad
state_vector_columns    planet, x, y, z, vx, vy, vz
dtype: object

Seeing the system

A top-down view ties the pieces together: integrate the system for one Uranus period (84 years) and trace each orbit. This is the geometry every longitude in this project is measured on.

import matplotlib.pyplot as plt
from IPython.display import display

from discoverneptune.plot_style import apply_style, category_cycle, dual_render

orbit_sim = build_simulation(best_neptune_parameters)

trace_names = ["Jupiter", "Saturn", "Uranus", "Neptune"]
trace_times_yr = np.linspace(0.0, 84.0, 600)
orbit_tracks = {name: np.empty((trace_times_yr.size, 2)) for name in trace_names}
for index, t_yr in enumerate(trace_times_yr):
    orbit_sim.integrate(t_yr)
    sun = orbit_sim.particles[0]
    for name in trace_names:
        body = orbit_sim.particles[name]
        orbit_tracks[name][index] = (body.x - sun.x, body.y - sun.y)


def build_orbit_figure(theme="light"):
    apply_style(theme)
    cycle = category_cycle(theme)
    fig, ax = plt.subplots(figsize=(7, 7))
    # Gruvbox gold for the Sun — distinct from every orbit color in both themes.
    ax.plot(0, 0, "o", ms=8, color="#d79921", zorder=5)
    ax.annotate("Sun", (0, 0), xytext=(6, 6),
                textcoords="offset points", fontsize=10, color="#d79921")
    for color, name in zip(cycle, trace_names, strict=False):
        ax.plot(orbit_tracks[name][:, 0], orbit_tracks[name][:, 1], lw=1.0, color=color)
        ax.plot(*orbit_tracks[name][-1], "o", ms=5, color=color)
        ax.annotate(name, orbit_tracks[name][-1], xytext=(6, 6),
                    textcoords="offset points", fontsize=10, color=color)
    ax.set_aspect("equal")
    ax.set_xlabel("x (AU)")
    ax.set_ylabel("y (AU)")
    ax.set_title("One Uranus period, top-down (1781–1865)")
    return fig


display(dual_render(build_orbit_figure, alt="Top-down traces of Jupiter, Saturn, Uranus, and Neptune orbits over one 84-year Uranus period"))