Historical Comparison: Le Verrier’s Observations vs Modern Simulation

This notebook compares the historical residuals reported by Le Verrier in 1846 with our modernized numerical simulation. It validates that the ‘signal’ Le Verrier saw—the systematic drift of Uranus—is reproduced by a modern gravitational model when Neptune is omitted.

Code

import numpy as np
import pandas as pd
from astropy.time import Time
from IPython.display import display

from discoverneptune.data import find_bundled_data_dir
from discoverneptune.simulation import (
    simulate_uranus_longitudes_no_neptune,
    StateVector,
)
from discoverneptune.plot_style import (
    apply_style,
    dual_render,
    fig_size,
    palette,
    truth_line,
)

apply_style()

1. Load Historical Data

We load the 1846 residuals extracted from Le Verrier’s memoirs. We use the project’s data-loading convention to locate the files.

Code

data_dir = find_bundled_data_dir()
hist_df = pd.read_csv(data_dir / 'leverrier_historical_observations.csv')

# Convert dates to Julian Dates for the integrator
# Note: Astropy Time requires a list or array of strings
hist_df['jd'] = Time(hist_df['date'].tolist()).jd
hist_df.head()

	date	year	observer	residual_arcsec	type	jd
0	1690-12-23	1690	Flamsteed	67.1	ancient	2338676.5
1	1712-03-22	1712	Flamsteed	92.7	ancient	2346435.5
2	1715-03-10	1715	Flamsteed	110.0	ancient	2347518.5
3	1750-09-25	1750	Lemonnier	19.0	ancient	2360501.5
4	1750-10-14	1750	Lemonnier	12.4	ancient	2360520.5

2. Run Modern Simulation

We integrate the solar system without Neptune, using JPL state vectors from 1781-01-01 as our initial conditions. This ‘forward model’ allows us to see how Uranus would have moved according to Newton if Neptune didn’t exist.

Code

svs = pd.read_csv(data_dir / 'outer_planet_state_vectors_1781_01_01.csv')
state = {
    row.planet: StateVector(row.x, row.y, row.z, row.vx, row.vy, row.vz)
    for row in svs.itertuples(index=False)
}

# Modern no-Neptune residuals over the span we have proxy observations for
# (1781–1846). Residual = observed − modeled, angle-wrapped to (−π, π].
#
# NOTE: do NOT reuse the notebook's existing `start_jd = 2371628.5` — that is
# the 13 Mar 1781 discovery-date row. The state vectors and the proxy
# observations share the 1781-01-01 epoch (JD 2371557.5), which is the first
# proxy row.
obs = pd.read_csv(data_dir / "uranus_observations_1781_1846.csv")
proxy_jd = obs["datetime_jd"].to_numpy(dtype=float)
proxy_lon = obs["longitude_rad"].to_numpy(dtype=float)
proxy_start_jd = float(proxy_jd[0])  # 2371557.5 = 1781-01-01

model_lons = simulate_uranus_longitudes_no_neptune(
    proxy_jd, proxy_start_jd, state["jupiter"], state["saturn"], state["uranus"]
)

wrapped = ((proxy_lon - model_lons + np.pi) % (2 * np.pi)) - np.pi
sim_residuals_arcsec = np.degrees(wrapped) * 3600.0
proxy_years = 1781.0 + (proxy_jd - proxy_start_jd) / 365.25

3. Visualize the Signal

The plot below shows Le Verrier’s reduced residuals (points) alongside the residual of a modern no-Neptune simulation against JPL-derived longitudes over 1781–1846 (dashed line). The two series are measured against different reference orbits, so their signs and shapes are not directly comparable: Le Verrier’s residuals are quoted against Bouvard’s fitted tables, while the modern curve compares to an exact 1781 state propagated without Neptune. What they share is what matters — both anomalies grow to order 100–130 arcseconds by 1846, surging after the 1822 Uranus–Neptune conjunction, far beyond what refitting any Keplerian orbit can absorb. The ancient points push the same conflict back another ninety years.

Code

def _build_historical(hist_df, proxy_years, sim_residuals_arcsec, *, theme="light"):
    import matplotlib.pyplot as plt

    pal = palette(theme)
    fig, ax = plt.subplots(figsize=fig_size("wide"))

    ancient = hist_df[hist_df["type"] == "ancient"]
    modern = hist_df[hist_df["type"].isin(["modern_discovery", "illustrative"])]

    ax.plot(
        ancient["year"], ancient["residual_arcsec"], "o",
        label="Ancient observations (Le Verrier 1846)",
        color=pal.before, alpha=0.7,
    )
    ax.plot(
        modern["year"], modern["residual_arcsec"], "s",
        label="Modern / illustrative trend",
        color=pal.after,
    )
    ax.plot(
        proxy_years, sim_residuals_arcsec, "--",
        label="Observed − no-Neptune model (1781–1846)",
        color=pal.spine, lw=1.5,
    )

    truth_line(ax, 0.0, axis="y", theme=theme)
    ax.set_title("Historical residuals vs. modern simulation")
    ax.set_xlabel("Year")
    ax.set_ylabel("Heliocentric longitude residual (arcseconds)")
    ax.legend()
    return fig


display(
    dual_render(
        _build_historical,
        hist_df,
        proxy_years,
        sim_residuals_arcsec,
        alt="Le Verrier's historical Uranus residuals overlaid with a modern no-Neptune simulation",
    )
)