Neptune and the Historical Inverse Problem

This notebook is a narrative walkthrough that stays narrative-first: it measures the historical data window, shows the Uranus anomaly in arcseconds, and uses a compact Neptune search to explain why longitude tightens faster than distance before handing off to rebound_neptune_tutorial.ipynb.

The Historical Problem

Uranus was tracked from 1781 through 1846, long enough for a perturbation to accumulate into a visible departure from a smooth orbit. The bundled notebook data keeps that period self-contained: 66 annual observations, no live Horizons query required.

Code 
from discoverneptune.plot_style import (
    apply_style, palette, category_cycle, dual_render,
    SEQUENTIAL_CMAP, truth_line, annotate_stat, fig_size,
)
apply_style()

from pathlib import Path

import numpy as np
import pandas as pd

from discoverneptune.simulation import (
    NeptuneCandidate,
    StateVector,
    simulate_uranus_longitudes_from_vectors,
)

cwd = Path.cwd()
observation_paths = [
    cwd / "data" / "uranus_observations_1781_1846.csv",
    cwd.parent / "data" / "uranus_observations_1781_1846.csv",
]
state_vector_paths = [
    cwd / "data" / "outer_planet_state_vectors_1781_01_01.csv",
    cwd.parent / "data" / "outer_planet_state_vectors_1781_01_01.csv",
]

for observation_path in observation_paths:
    if observation_path.exists():
        break
else:
    raise FileNotFoundError(
        "Run this notebook from the repository root or the docs/ directory "
        "so the bundled Uranus observations can be found."
    )

for state_vector_path in state_vector_paths:
    if state_vector_path.exists():
        break
else:
    raise FileNotFoundError(
        "Run this notebook from the repository root or the docs/ directory "
        "so the bundled state vectors can be found."
    )

obs = pd.read_csv(observation_path)
state_vectors = pd.read_csv(state_vector_path)
state = {
    row.planet: StateVector(row.x, row.y, row.z, row.vx, row.vy, row.vz)
    for row in state_vectors.itertuples(index=False)
}

jd_times = obs["datetime_jd"].to_numpy()
truth_longitudes = obs["longitude_rad"].to_numpy()
start_jd = float(jd_times[0])
end_jd = float(jd_times[-1])
cadence_days = obs["datetime_jd"].diff().dropna()
timeline = pd.DataFrame(
    {
        "metric": [
            "first observation (JD)",
            "last observation (JD)",
            "span (years)",
            "median cadence (days)",
        ],
        "value": [
            start_jd,
            end_jd,
            (end_jd - start_jd) / 365.25,
            float(cadence_days.median()),
        ],
    }
)
timeline
metric value
0 first observation (JD) 2.371558e+06
1 last observation (JD) 2.395298e+06
2 span (years) 6.499658e+01
3 median cadence (days) 3.650000e+02

What Le Verrier Had Access To

Le Verrier worked from observation tables and hand-built theory, not from a modern ephemeris or a numerical integrator. In this reconstruction we keep the inner problem simple too: Jupiter, Saturn, and Uranus start from JPL state vectors, and only Neptune is left as an unknown.

From Residuals to an Inverse Problem

The forward model takes a small Neptune vector, propagates the outer planets, and compares the predicted Uranus longitude against the observations. The unknowns are deliberately compact: mass, semi-major axis, eccentricity, and mean longitude.

Why Longitude Is Easier Than Distance

Longitude is the cleanest observable in the bundled data set, so the notebook measures the perturbation in angular residuals. A modern grid search lets us ask a sharper question: if we keep only the best few fits, how tightly do they cluster in longitude compared with semi-major axis and mass? That comparison is the computational version of Le Verrier’s practical advantage: sky position tightens faster than orbital scale.

Code 
candidate = NeptuneCandidate(
    mass_solar=5.0e-5,
    a=30.1,
    e=0.009,
    l_rad=np.deg2rad(304.8),
)
model_longitudes = simulate_uranus_longitudes_from_vectors(
    jd_times,
    start_jd,
    state["jupiter"],
    state["saturn"],
    state["uranus"],
    candidate,
)
residual_arcsec = (
    np.degrees(np.unwrap(truth_longitudes - model_longitudes)) * 3600.0
)
residual_summary = pd.DataFrame(
    {
        "statistic": ["min", "max", "rms", "peak_to_peak"],
        "arcsec": [
            float(residual_arcsec.min()),
            float(residual_arcsec.max()),
            float(np.sqrt(np.mean(residual_arcsec**2))),
            float(np.ptp(residual_arcsec)),
        ],
    }
)
residual_summary
statistic arcsec
0 min -93.781478
1 max 39.626141
2 rms 28.659162
3 peak_to_peak 133.407619

Results: The Mass-Distance Degeneracy

The next cells keep one best longitude for each mass/a pair, then compare the top-performing solutions as a group. In this toy search the best fits keep longitude within a few degrees, even while semi-major axis and mass still move across a visibly broader valley. That is the mass-distance tradeoff the notebook is meant to surface.

Code 
mass_distance_valley = results.pivot(
    index="mass_solar", columns="a_au", values="rms_arcsec"
)
mass_distance_valley
a_au 29.3 29.8 30.3 30.8 31.3
mass_solar
0.000045 4.292753 3.292204 3.370921 4.023064 5.118561
0.000048 5.319473 3.923111 3.266632 3.540353 4.327330
0.000051 6.398116 4.942278 3.666436 3.297752 3.781378
0.000054 7.708198 5.865429 4.501412 3.499409 3.397773
Code 
best_toy_grid = competitive_solutions.iloc[0]
summary = pd.DataFrame(
    {
        "quantity": [
            "best semi-major axis (AU)",
            "best mass (solar masses)",
            "best eccentricity",
            "best longitude (deg)",
            "competitive longitude span (deg)",
            "competitive semi-major-axis span (AU)",
            "competitive mass span (solar masses)",
        ],
        "value": [
            float(best_toy_grid["a_au"]),
            float(best_toy_grid["mass_solar"]),
            float(best_toy_grid["e"]),
            float(best_toy_grid["l_deg"]),
            longitude_span_deg,
            semi_major_axis_span_au,
            mass_span_msun,
        ],
    }
)
summary
quantity value
0 best semi-major axis (AU) 30.300000
1 best mass (solar masses) 0.000048
2 best eccentricity 0.009000
3 best longitude (deg) 191.000000
4 competitive longitude span (deg) 3.000000
5 competitive semi-major-axis span (AU) 1.500000
6 competitive mass span (solar masses) 0.000009

Bridge to the REBOUND Tutorial

The hands-on companion notebook, rebound_neptune_tutorial.ipynb, takes the same bundled data and rebuilds the search step by step inside REBOUND.

The result in one figure

Everything above compresses to this: remove Neptune from the model and Uranus drifts ~130 arcseconds out of place; add the optimizer’s Neptune back and the drift collapses to the few-arcsecond model noise floor. The quick search budget used here lands at ~2.9″ RMS; the production optimizer (see The Uranus Problem) reaches ~2.5″.

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

from discoverneptune.simulation import simulate_uranus_longitudes_no_neptune
from discoverneptune.solver import solve

_baseline = simulate_uranus_longitudes_no_neptune(
    jd_times, start_jd, state["jupiter"], state["saturn"], state["uranus"]
)
_fit = solve(
    jd_times, truth_longitudes, start_jd,
    state["jupiter"], state["saturn"], state["uranus"],
    de_maxiter=30, de_popsize=5, polish=True,
)
_best = simulate_uranus_longitudes_from_vectors(
    jd_times, start_jd, state["jupiter"], state["saturn"], state["uranus"],
    _fit.neptune,
)


def _wrap(delta):
    return ((delta + np.pi) % (2 * np.pi)) - np.pi


_res_before = np.degrees(_wrap(truth_longitudes - _baseline)) * 3600.0
_res_after = np.degrees(_wrap(truth_longitudes - _best)) * 3600.0
_years = 1781.0 + (jd_times - start_jd) / 365.25
_rms_after = float(np.sqrt(np.mean(_res_after**2)))


def _build(theme="light"):
    apply_style(theme)
    pal = palette(theme)
    fig, ax = plt.subplots(figsize=(9.5, 4))
    ax.plot(_years, _res_before, "o-", lw=0.8, ms=3, color=pal.before,
            label="No Neptune — the anomaly")
    ax.plot(_years, _res_after, "o-", lw=0.8, ms=3, color=pal.after,
            label=f"With fitted Neptune — RMS {_rms_after:.1f}″")
    ax.axhline(0, color=pal.spine, lw=0.5, ls=":")
    ax.set_xlabel("Year")
    ax.set_ylabel("Uranus longitude residual (arcsec)")
    ax.set_title("The rediscovery, in one figure")
    ax.legend(loc="lower left")
    return fig


display(dual_render(_build, alt="Uranus longitude residuals before and after adding the fitted Neptune"))