Interpolated Perigee Methodology #

LibEphemeris computes the interpolated lunar perigee using a three-layer architecture — mean position, 61-term harmonic perturbation series, and residual correction table — calibrated against JPL DE441 ephemeris via passage-interpolated harmonic fitting.

Table of Contents #

• Background
• Method
  • Three-Layer Architecture
  • Passage-Interpolated Harmonic Fitting (v2.2)
  • Perturbation Series (61 Terms)
  • Residual Correction Table
• Failed Approaches
  • Gaussian Smoothing (v0)
  • Raw Osculating Fit (v2.0)
  • Passage-Only Fit (v2.1)
• Precision and Validation
  • Key Coefficient Changes (v1 to v2.2)
  • Precision Evolution
• Comparison with Swiss Ephemeris
• Calibration Reproducibility
• References

Background #

The Moon’s orbit around Earth is continuously perturbed by the Sun. These perturbations cause:

1. The orbital axis (line of apsides) to precess with a period of ~8.85 years
2. The perigee (point of closest approach) to oscillate around this mean precession

These oscillations are substantial: approximately ±25 degrees from the mean position.

The term “interpolated perigee” refers to a smoothed version of the instantaneous (osculating) perigee position, with short-period oscillations removed. There are different valid approaches to defining and computing this smoothing.

The combined three-layer result achieves ~1.5 degree RMS agreement with Swiss Ephemeris (down from ~10–11 degrees in v1), while being grounded entirely in JPL ephemeris data rather than analytical lunar theory.

Method #

Three-Layer Architecture #

The interpolated perigee is computed as the sum of three layers:

interpolated_perigee = mean_perigee(t) + perturbation_series(t) + residual_correction(t)

Note: the first layer (mean perigee) includes its own JPL correction table (MEAN_APSE_CORRECTIONS).

Layer 1: Mean Perigee (Polynomial)

The mean perigee position uses a standard polynomial in Julian centuries T from J2000:

ω₀ = 83.3532° + 4069.0137° × T + ...  (higher-order terms)

This captures the secular apsidal precession (~40.7° per year). The mean perigee is 180° from the mean apogee (Mean Lilith).

Layer 2: Perturbation Series (61 terms)

The perturbation series _calc_elp2000_perigee_perturbations() adds periodic corrections organized by physical origin. See Perturbation Series (61 Terms) below.

Layer 3: Residual Correction Table

After applying the perturbation series, residual errors (~2 deg RMS near J2000, ~8 deg over the full range) are absorbed by a precomputed correction table (lunar_corrections.py, 15,195 entries). Linear interpolation between table entries brings the final error to < 0.1 degrees.

The correction table is generated by scripts/generate_lunar_corrections.py, which:

1. Computes JPL ground-truth perigee via ±28 day quadratic regression of osculating elements
2. Computes model perigee (mean + perturbation series)
3. Stores the difference at regular time intervals

Passage-Interpolated Harmonic Fitting (v2.2) #

The v2.2 method combines the strengths of earlier approaches while avoiding their failures (see Failed Approaches):

1. Find perigee passages — Locate all Earth-Moon distance minima over a 1000-year span [1500, 2500] CE using JPL DE441. At each passage, the Moon’s ecliptic longitude is an unambiguous measure of the perigee longitude (12,958 passages found).

2. Cubic spline interpolation — Fit a cubic spline through the (time, perigee_longitude) passage data, with angle unwrapping to handle the 0°/360° discontinuity. This creates a smooth, continuous perigee longitude function.

3. Daily resampling — Sample the spline at daily intervals to produce ~365K clean data points with full coverage of all Delaunay argument combinations (unlike passage-only data where M’ ≈ 0 always).

4. Harmonic least-squares fit — Construct a 120-term design matrix of candidate trigonometric terms from lunar theory and fit via ordinary least squares. Terms with |coefficient| < 0.001 degrees are discarded, yielding the final 61-term series.

The spline interpolation acts as a physically motivated smoothing: it passes exactly through the ground-truth passage points while providing well-conditioned data at arbitrary intermediate times. The daily samples have full M’ coverage (unlike passages) and no correlated osculating noise (unlike raw data).

Perturbation Series (61 Terms) #

The polynomial correction terms (const, T, T², T³) absorb secular drift between the mean perigee model and the JPL-calibrated series.

Residual Correction Table #

After the 61-term perturbation series, residual errors are absorbed by a precomputed correction table with the following parameters:

The residual captures secular drift of the perturbation coefficients and harmonics not included in the 61-term series. At dates outside the DE441 range, the system falls back to the perturbation series alone (gradual degradation, not failure).

Failed Approaches #

Three earlier approaches to computing the perturbation series coefficients were attempted and failed, each providing insight that informed the successful v2.2 method.

Gaussian Smoothing (v0) #

The first attempt used Gaussian smoothing of the osculating perigee with an 8-year window (2922 days). This approach failed because the apsidal line precesses ~325 degrees in 8 years (nearly a full revolution). Consequently, the samples within the Gaussian window covered nearly the entire circle, causing the circular mean computation (atan2(sum_sin, sum_cos)) to collapse near zero. The result was ~-0.16 degrees of perturbation — effectively zero — when the correct value is ~8–25 degrees.

The fundamental problem is that a Gaussian window spanning nearly a full apsidal precession cycle cannot recover the phase of an oscillation whose period is comparable to the window width.

Raw Osculating Fit (v2.0) #

The second attempt fitted a harmonic series directly to ~365K daily osculating perigee samples. This failed because the osculating noise is correlated with the design matrix terms — both are functions of the same Delaunay arguments (D, M, M’, F). The least-squares fit could not distinguish genuine perturbations from osculating artifacts, producing nonsensical coefficients with a passage RMS of 9.18 degrees.

Passage-Only Fit (v2.1) #

The third attempt fitted only at perigee passages (Earth-Moon distance minima), where the Moon’s longitude is the perigee longitude by definition. This failed because at every passage, the mean anomaly M’ ≈ 0 by definition, making the design matrix degenerate (condition number 3.7×10⁵). The coefficients of terms involving M’ were unconstrained.

The v2.2 method resolved this by using the passages to anchor a spline, then resampling the spline at daily intervals to recover full Delaunay argument coverage.

Precision and Validation #

Key Coefficient Changes (v1 to v2.2) #

The v1 coefficients were severely attenuated because v1 used quadratic regression smoothing of osculating data, which damped the harmonic amplitudes by ~50%:

The dramatic increase in the dominant sin(2D-2M’) term (evection) from -9.62° to -22.21° is the most significant change. The old value was roughly half the true amplitude due to signal attenuation by the quadratic regression window.

Precision Evolution #

The remaining ~1.5° difference vs. Swiss Ephemeris is expected: it reflects the fundamental methodological difference between LibEphemeris (JPL-grounded numerical approach) and Swiss Ephemeris (ELP2000-82B analytical theory).

Comparison with Swiss Ephemeris #

Swiss Ephemeris Approach (Moshier/ELP2000-82B) #

Swiss Ephemeris uses semi-analytical perturbation theory based on the ELP2000-82B lunar theory developed by Chapront-Touzé and Chapront. Perturbations are separated analytically based on their physical origin: “short-period” perturbations (periods less than a few months) are removed, and “long-period” perturbations that represent the true apsidal motion are retained.

LibEphemeris Approach (JPL DE441 + Passage-Interpolated Fitting) #

LibEphemeris uses an empirical harmonic series calibrated against numerical ephemeris. Perigee passages are identified from JPL DE441 Moon positions (physical ground truth), spline interpolation between passages creates a smooth perigee longitude function, and harmonic series coefficients are fitted empirically to this function. A residual correction table absorbs remaining model error.

Why the ~1.5° Residual Difference Exists #

The two approaches differ in what is smoothed and how:

The ~1.5° difference is a consequence of different smoothing philosophies applied to the same underlying phenomenon, not an error in either implementation.

Analogy: Consider audio signal processing with bass and treble components. The Swiss Ephemeris approach uses an intelligent filter that knows which frequencies represent “signal” vs. “noise” based on the source characteristics. The LibEphemeris approach identifies the signal at known clean points (passage times), interpolates between them, then fits a harmonic model to the result. Both produce a smoothed result, but the outputs differ because the smoothing methodology differs.

Date Range #

Set LIBEPHEMERIS_EPHEMERIS=de441.bsp for extended date range coverage.

Calibration Reproducibility #

The calibration can be reproduced using:

LIBEPHEMERIS_EPHEMERIS=de441.bsp python scripts/calibrate_perigee_perturbations.py \
    --start-year 1500 --end-year 2500 --output /tmp/perigee_v22_full.json

The script outputs calibrated coefficients as Python code ready to paste into lunar.py. The correction table can be regenerated with:

python scripts/generate_lunar_corrections.py

References #

1. Chapront-Touzé, M. & Chapront, J. (1988). “ELP 2000-82B: A semi-analytical lunar ephemeris.” Astronomy & Astrophysics, 190, 342-352.
2. Chapront-Touzé, M. & Chapront, J. (1991). Lunar Tables and Programs from 4000 B.C. to A.D. 8000. Willmann-Bell.
3. Park, R.S. et al. (2021). “The JPL Planetary and Lunar Ephemerides DE440 and DE441.” Astronomical Journal, 161(3), 105.
4. Moshier, S.L. (1989–2023). Ephemeris Calculation Software. Available at http://www.moshier.net — independent C implementation of ELP2000-82B used as reference for the analytical series coefficients.