Building a Clean-Room Ephemeris

Why provenance matters for an astrological engine

Astrology software has accumulated decades of copying. A rounding convention introduced in one library propagates to another, and by the time someone traces a discrepancy to its origin, the original code is gone. A clean-room implementation breaks that chain. You can point at equation (3.1) in the IERS Conventions 2010 and say: that is exactly what the code computes.

For Vedic astrology this matters even more. Ayanamsha values differ by arcseconds between implementations — arcseconds that shift nakshatra boundaries and change dasha periods. When a practitioner asks why their chart shows a different nakshatra than another tool, the answer should be traceable to a specific ayanamsha definition and reference epoch, not to an undocumented behavioral quirk.

Starting with JPL DE440

Planetary positions start with the JPL Developmental Ephemeris. DE440 is the current long-arc solution, covering 1550 to 2650 CE, developed by Park, Folkner, Williams, and Boggs (2021). We use the SPK kernel format — binary files containing Chebyshev polynomial coefficients for each body over each time interval.

The SPK reader is implemented from scratch against the SPICE toolkit documentation. A record maps to a time window and a set of coefficients. Position evaluation is three nested Chebyshev recurrences — one per Cartesian component — using the standard recurrence relation T_n+1(x) = 2x·T_n(x) − T_n-1(x). Velocity follows from the derivative recurrence without finite differences.

This gives barycentric positions in the ICRS frame with sub-kilometer accuracy across the covered arc. The Sun is recovered as the difference between the Solar System Barycenter and the Earth-Moon Barycenter, adjusted by the lunar mass fraction.

fn cheby_eval(coeffs: &[f64], x: f64) -> f64 {
    let (mut t0, mut t1) = (1.0, x);
    let mut result = coeffs[0] * t0 + coeffs[1] * t1;
    for c in &coeffs[2..] {
        let t2 = 2.0 * x * t1 - t0;
        result += c * t2;
        (t0, t1) = (t1, t2);
    }
    result
}

Clenshaw recurrence for Chebyshev evaluation — the core of every planetary position computation.

The coordinate pipeline: ICRS to ecliptic

Barycentric ICRS coordinates are not what astrology needs. The pipeline from raw SPK output to an ecliptic longitude has six steps, each with its own error budget.

Delta T: the most underestimated problem

JPL ephemerides are tabulated in Barycentric Dynamical Time (TDB). Astrological inputs are in local civil time — UTC or a named timezone. The conversion chain is: local time → UTC → UT1 → TT → TDB. Each step has different uncertainty characteristics.

Delta T (ΔT = TT − UT1) is the accumulated difference between atomic time and Earth rotation, currently around 69 seconds. For modern dates it comes from IERS Bulletin A. For historical dates we use the Morrison-Stephenson polynomial model extended by Espenak-Meeus. For dates before 1620, the secular term dominates and uncertainty grows to minutes — a fact that matters enormously for historical chart rectification.

ΔT interpolates measured IERS values from 1620 to 2025 and falls back to the Espenak & Meeus (2006) polynomial expressions outside that range. Its unit tests check the interpolation against that table across roughly 1500–2030; they do not cross-validate the table itself against an independent ΔT source, so we make no accuracy claim for ΔT beyond what IERS publishes.

That honesty matters more than it sounds. ΔT past the measured record is not a solved number — it depends on how the Earth's rotation actually behaves. Our extrapolation and JPL Horizons' own prediction diverge by roughly 68 seconds at 2099, which lands as a longitude error proportional to a body's angular rate: about 45″ for the Moon, essentially nothing for Pluto. That is a disagreement about time, not about the ephemeris — and it is why every accuracy figure on this site is scoped to 1900–2025, where ΔT is measured rather than guessed.

The implementation uses Rust's f64 throughout — 64-bit double precision is sufficient for all ephemeris calculations given that the underlying SPK data carries ~14 significant digits.

Validation against JPL Horizons

The definitive validation target is JPL Horizons — an independent implementation serving DE441, while our SPK reader consumes DE440s, so the kernels are genuinely independent. The suite asserts against a committed oracle fixture of 24,350 Horizons reference positions (10 bodies × 2,435 dates, 1900–2100), regenerable via scripts/generate_horizons_oracle.py. The lunar theory (ELP/MPP02) is additionally checked against a live-fetched Horizons grid spanning −3000 to +3000 CE.

Results with the SpkReader provider (DE440s kernel): over 1900–2025 — the era where ΔT is measured rather than extrapolated — the mean residual against Horizons DE441 is 0.106″ and the max is 1.184″ (Uranus, 1948); 15,349 of 15,350 comparisons are sub-arcsecond. For the Moon specifically, the ELP/MPP02 residual against Horizons DE441 is 0.015″ at J2000, against a 0.06″ tolerance, and stays between 0.020″ and 0.053″ across 1500–2500 CE. The AnalyticalProvider — which carries no data file and runs in WASM — agrees with Horizons DE441 over 1900–2025 to 0.17″ mean (0.61″ max) for the Moon, 4.09″ mean (7.00″ max) for the Sun (VSOP87A truncation; Sun is approximated as the solar system barycentric origin), and 2.06″ mean overall, with Venus the worst case at 24.22″ max. Residuals trace cleanly to known truncation in VSOP87A and the Sun ≈ SSB approximation, not to algorithmic errors.

Beyond 2025, residuals against Horizons grow — up to 44.914″ for the Moon by 2099 — but this is ΔT prediction divergence, not ephemeris error. Our Espenak & Meeus extrapolation and Horizons' own ΔT model differ by roughly 68 seconds at 2099, which shows up as a longitude offset proportional to each body's angular rate: about 45″ for the fast-moving Moon, essentially nothing for Pluto. At a 2099 test date, the Sun, Moon, Mercury, Venus, and Mars residuals all imply the same 66–71 s ΔT offset. ΔT beyond the measured IERS record is unpredictable in principle — no ephemeris can do better without a time machine.

Primary sources cited in the implementation