Ephemeris Engines: Data Sources, Interpolation, and Numerical Methods #

When an astrology engine needs the position of Mars at a specific instant, it doesn’t solve Newton’s equations of motion in real time. Instead, it looks up the answer in pre-computed data — a set of high-precision planetary positions calculated by astronomers using months of supercomputer time. The software that reads this data and returns a position for any requested time is the ephemeris engine.

This article covers the three main data sources (JPL DE, VSOP87, and Moshier), the mathematics of Chebyshev polynomial interpolation used to read JPL files, and the architecture of a typical engine.

Data Sources #

JPL Development Ephemerides (DE440/DE441) #

The NASA Jet Propulsion Laboratory produces the definitive solar system ephemerides. DE440 (covering 1550–2650 CE) and DE441 (covering −13200 to +17191 CE) are the current standards.

These are generated by numerically integrating the equations of motion for all major solar system bodies simultaneously, including:

• Newtonian gravity between all bodies
• General relativistic corrections (PPN formalism)
• Solar oblateness ($J_2$)
• Lunar libration and solid-body tides
• Perturbations from 343 asteroids

The result is a set of barycentric position and velocity vectors for the Sun, Moon, and planets, sampled at fine time intervals and then compressed into Chebyshev polynomial coefficients for compact storage and fast interpolation.

File format: the time span is divided into fixed-length blocks (typically 32 days for planets, 4–8 days for the Moon). Each block stores a set of Chebyshev coefficients for each coordinate $(X, Y, Z)$ of each body. The file size is ~100 MB for DE440.

Precision: sub-milliarcsecond for the modern era, degrading to ~arcsecond level for dates thousands of years from present.

VSOP87 #

VSOP87 (Variations Séculaires des Orbites Planétaires, Bretagnon & Francou 1988) is an analytical planetary theory. Instead of numerically integrated data files, it expresses planetary coordinates as sums of trigonometric series:

$$X(T) = \sum_{i=0}^{5} T^i \sum_{j} A_{i,j} \cos(B_{i,j} + C_{i,j},T)$$

where $T$ is Julian centuries from J2000.0, and the coefficients $A, B, C$ are tabulated. The series for Earth alone has over 2,000 terms.

Advantages: no external data files needed — the coefficients can be compiled directly into the program. This makes VSOP87 ideal for embedded systems or environments where a 100 MB file is impractical.

Disadvantages: evaluating thousands of trigonometric terms is slower than a Chebyshev lookup, and precision degrades outside the ±4000 year range from J2000 (~0.01" at J2000, increasing to ~1" at ±2000 years).

VSOP87 has several versions (VSOP87A through VSOP87E) providing coordinates in different reference frames. VSOP87C gives heliocentric ecliptic of-date, while VSOP87A gives heliocentric ecliptic J2000 — make sure you know which version your engine uses.

ELP/MPP02 #

ELP2000-82 and its successor MPP02 are the analytical theories for the Moon. The Moon’s motion is too complex for VSOP87 to handle (due to strong solar perturbations and the Moon’s close proximity to Earth). ELP/MPP02 has over 37,000 terms.

Moshier #

Steve Moshier produced a truncated, standalone C implementation of VSOP87 and ELP2000 that requires no external data files and fits in a few hundred KB of compiled code. It trades ~0.1–2" of precision for zero file dependencies.

Many ephemeris engines use Moshier as the fallback when JPL data files are not found:

Query: position of Mars at JD(TT)
    │
    ├─ JPL DE440 file present? ──YES──→ Read Chebyshev block, evaluate
    │
    └─ NO ──→ Evaluate Moshier analytical series

Chebyshev Polynomial Interpolation #

The JPL ephemeris files encode positions as Chebyshev polynomial coefficients. Understanding how to evaluate them is essential for anyone implementing or debugging an ephemeris engine.

Chebyshev Polynomials #

The Chebyshev polynomials of the first kind $T_n(x)$ are defined on $[-1, 1]$ by the recurrence:

$$T_0(x) = 1, \quad T_1(x) = x, \quad T_n(x) = 2x,T_{n-1}(x) - T_{n-2}(x)$$

A function $f(x)$ on $[-1, 1]$ can be approximated as:

$$f(x) \approx \frac{c_0}{2} + \sum_{i=1}^{N-1} c_i,T_i(x)$$

Chebyshev approximation has the minimax property: among all polynomial approximations of degree $N$, the Chebyshev expansion minimizes the maximum error over the interval. This is why JPL uses them — a degree-14 Chebyshev polynomial can represent 32 days of planetary motion to sub-milliarcsecond accuracy.

Reading a JPL Block #

To evaluate the $X$-coordinate of Mars at time $t$:

1. Find the block covering time $t$. Each block spans $[t_{\text{start}}, t_{\text{end}}]$.
2. Normalize $t$ to $[-1, 1]$: $$x = \frac{2(t - t_{\text{start}})}{t_{\text{end}} - t_{\text{start}}} - 1$$
3. Read the $N$ Chebyshev coefficients $c_0, c_1, \ldots, c_{N-1}$ for this block.
4. Evaluate the sum using the Clenshaw algorithm (below).

Clenshaw Algorithm #

Direct evaluation of $\sum c_i,T_i(x)$ by computing each $T_i(x)$ explicitly is wasteful. The Clenshaw algorithm evaluates the sum using only the recurrence relation, requiring $O(N)$ multiplications and no storage of intermediate polynomials:

Set $b_{N+1} = 0$, $b_N = 0$. Then for $i = N-1, N-2, \ldots, 1$:

$$b_i = c_i + 2x,b_{i+1} - b_{i+2}$$

The result is:

$$f(x) = c_0 + x,b_1 - b_2$$

(Note: some implementations use the convention $f(x) = \frac{c_0}{2} + x,b_1 - b_2$. Verify which convention your data source uses.)

Velocity from Chebyshev Coefficients #

The derivative of the Chebyshev expansion gives the velocity. The derivative of $T_n(x)$ follows a similar recurrence:

$$T’0(x) = 0, \quad T’1(x) = 1, \quad T’n(x) = 2T{n-1}(x) + 2x,T’{n-1}(x) - T’{n-2}(x)$$

Or, more efficiently, using the identity $T’n(x) = n,U{n-1}(x)$ where $U$ is the Chebyshev polynomial of the second kind, and evaluating via a modified Clenshaw. Most implementations compute position and velocity in a single pass.

Engine Architecture #

A complete ephemeris engine follows this pipeline:

Input: JD(TT), body ID
    │
    ▼
[Time normalization: JD(TT) → T → block index → x ∈ [-1,1]]
    │
    ▼
[Chebyshev evaluation: c[i] → X, Y, Z (barycentric)]
    │
    ▼
[Frame shift: barycentric → geocentric (subtract Earth vector)]
    │
    ▼
[Coordinate conversion: cartesian → spherical (λ, β, r)]
    │
    ▼
Output: ecliptic longitude, latitude, distance, velocities

The geocentric step is critical and easy to get wrong. JPL DE files store positions relative to the solar system barycenter. To get geocentric positions (what astrology needs), you must:

1. Evaluate the Earth-Moon barycenter (EMB) position from the Chebyshev data
2. Evaluate the Moon’s position relative to EMB
3. Compute Earth = EMB − Moon × (Moon mass ratio)
4. Subtract Earth from the target body

Precision Comparison #

For astrological purposes, even Moshier’s ~0.1" precision is vastly more than needed — the Ascendant changes by about 1° per 4 minutes of clock time, so a 0.1" planetary position error is negligible. The choice of data source matters more for validation and for historically distant dates.

References #

• Standish, E. M. (1998). “JPL Planetary and Lunar Ephemerides, DE405/LE405.” JPL IOM 312.F-98-048.
• Park, R. S. et al. (2021). “The JPL Planetary and Lunar Ephemerides DE440 and DE441.” The Astronomical Journal, 161, 105.
• Bretagnon, P., & Francou, G. (1988). “Planetary theories in rectangular and spherical variables. VSOP87 solutions.” Astronomy and Astrophysics, 202, 309–315.
• Meeus, J. (1998). Astronomical Algorithms, 2nd ed. Willmann-Bell. Chapters 33, 47 (VSOP87 implementation).
• Press, W. H. et al. (2007). Numerical Recipes, 3rd ed. Cambridge University Press. Section 5.8 (Chebyshev approximation).