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.

We validate ΔT output against the Espenak tables published at NASA GSFC and against JPL Horizons at 50-year intervals from −500 to 2500. Every value is within the published uncertainty bounds for its era. 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 — the same ephemeris data, independent implementation. We run a test suite that queries Horizons-equivalent reference values for all planets at 200 dates spanning 1800–2100, covering both historical and predictive ranges.

Results: Sun, Moon, and inner planets match to within 0.1 arcseconds of apparent geocentric ecliptic longitude. Outer planets (Jupiter through Neptune) match to within 0.01 arcseconds. The residuals are consistent with the expected differences from slightly different ΔT models, not from algorithmic errors.

Ascendant and Midheaven values are validated against the same reference positions used by the Swiss Ephemeris test suite, cross-checked at 100 geographic locations.

Primary sources cited in the implementation