From ephemeris to elements

The question. Given a modern ephemeris such as JPL DE440, what sequence of equations produces a table of Besselian elements, which reference frame the inputs must be in, which of light-time, aberration, precession and nutation they must carry, how the published polynomial is fitted, and when a developer should skip the elements and compute topocentric circumstances directly. The definitions and constants are in Besselian elements and the fundamental plane and are not repeated.

What the inputs are

The 1961 Supplement is explicit: "The basic quantities from which the calculations of solar and lunar eclipses are derived are: the apparent right ascension $\alpha_\odot$, declination $\delta_\odot$, and radius vector $R$ of the Sun, and the right ascension $\alpha_\mathrm{m}$, declination $\delta_\mathrm{m}$, and horizontal parallax of the Moon, for every hour of E.T. during the eclipse; and the ephemeris sidereal time" ¹. An apparent placeapparent placeA body's direction as seen from the Earth's centre at an instant, including light-time, aberration, precession and nutation, referred to the true equator and equinox of date. is "the proper place of an object expressed with respect to the true (intermediate) equator and equinox of date", as opposed to a geometric positiongeometric positionThe straight-line direction from the observer to a body at an instant with no correction for light-time or aberration., which has no "corrections for light-time, aberration, etc." ¹¹. The same Supplement notes that "Besselian elements, however, are rigorously independent of refraction" ¹.

Three points follow for anyone starting from a JPL binary ephemeris, which delivers barycentric ICRF positions in TDB.

1. Frame. The fundamental system is defined through the equator and through $\mu$, an hour angle from a meridian. The x axis lies in the equatorial plane. Both need the true equator and equinox of date. A developer must apply precession and nutation to rotate GCRS vectors to the equator of date before forming $a$ and $d$. Current practice is IAU 2006/2000A ¹⁰. The Swiss Ephemeris does this by requesting equatorial coordinates without the SEFLG_J2000 flag ². The solareclipses.com generator says "Use Equatorial Apparent of the Date coordinates" ³. The Swiss Ephemeris comment "nutation need not be in lunar and solar positions, if mean sidereal time will be used" records the one consistent shortcut: nutation can be omitted from both the positions and the sidereal time, but not from one alone ².
2. Light-time and aberration. The almanac inputs are apparent places, so they carry light-time and annual aberration. The Swiss Ephemeris eclipse_where() calls swe_calc() with its default flags, which return apparent positions with light-time and aberration ^{2 12}. The solareclipses.com page states "Light-time already embedded in ephemeris data" ³. Quaglia et al., who work topocentrically, apply "(a) light time to the body in question; (b) gravitational bending of the light path by the mass of the Sun; and (c) planetary aberration" ¹⁰. Annual aberration displaces the Sun and the Moon by nearly the same angle, up to 20.5 arcseconds, so it moves the point Z but changes $x, y$ by only the differential part, of order the parallax factor times 20 arcseconds. Whether to use apparent or astrometric places is therefore a question of consistency with the sidereal time and the observer model rather than of shadow physics. None of the sources read for this note quantifies the difference. That is an open question below.
3. Time. Evaluate the ephemeris at instants of TT. The difference between TT and TDB, under 2 milliseconds, is negligible for the Moon's motion of 0.5 arcseconds per second. Do not evaluate at UT: ΔT belongs in the conversion of $\mu$ to longitude, not in the argument of the ephemeris ^{1 13}.

Espenak's VSOP87/ELP2000-82 series are "referred to the mean equinox of the date" for the Sun ¹⁴, and his later series use JPL DE405 with the Moon's centre of mass ¹⁵. The lunar theory was truncated at 0.0005 arcseconds and 1 metre, giving mean errors of "0.0006 second of time in right ascension and 0.006 arcsecond in declination" for the Moon and contact-time errors "of the order of 1/40 second" ¹⁴.

The equations, in order

Symbols used below.

• $\alpha_\odot, \delta_\odot, R$: right ascension, declination and distance of the Sun, distance in au.
• $\alpha_\mathrm{m}, \delta_\mathrm{m}, \pi_\mathrm{m}$: right ascension, declination and equatorial horizontal parallax of the Moon.
• $\pi_0 = 8.794$ arcseconds: the solar parallax at 1 au.
• $s_0 = 959.63$ arcseconds: the solar semidiameter at 1 au.
• $k$: the lunar radius ratio, one value per cone.
• $\theta$: the Greenwich apparent sidereal time at the TT instant less ΔT.

All lengths are in Earth equatorial radii unless stated. Sources: the 1961 Supplement section 9B and the 1992 Supplement sections 8.322 and 8.323 ^{1 5}.

Step 0. Positions. From the ephemeris at TT instant $t$, obtain apparent geocentric $(\alpha_\odot, \delta_\odot, R)$ and $(\alpha_\mathrm{m}, \delta_\mathrm{m}, r_\mathrm{m})$ on the true equator and equinox of date. If the Moon's distance is given in km, form $\sin\pi_\mathrm{m} = a_e / r_\mathrm{m}$ with $a_e = 6378.137$ km, so that $r_\mathrm{m} = 1/\sin\pi_\mathrm{m}$ in Earth radii ⁵.

Step 1. Scale. The 1961 Supplement works with $g = G/R$ and $b = r_\mathrm{m}/R$, evaluated as

Step 2. Direction of the axis (point Z).

Step 3. Moon in the fundamental system.

Step 4. Hour angle.

Step 5. Cone half-angles.

Step 6. Vertex heights.

Step 7. Shadow radii on the fundamental plane.

Step 8. Repeat steps 0 to 7 at each sample instant. The 1992 Supplement used a 10-minute grid "but for some instances additional calculations at an interval of 30 seconds or less are necessary", and took derivatives numerically ⁵.

Step 9. Fit. For publication, sample at $t = -2, -1, 0, +1, +2$ hours about $t_0$ and fit $a_0 + a_1 t + a_2 t^2 + a_3 t^3$ by least squares to each of $x, y, d, \mu, l_1, l_2$. NASA describes "a least-squares fit to elements calculated at five uniformly spaced times over a six hour period centered at $t_0$" ^{6 17}. Miller's page shows the normal equations solved by Gauss-Jordan elimination and reports that $\mu$ needs only first degree and $d, l_1, l_2$ only second ⁴. The NASA CSV confirms: $\mu$ has a zero quadratic coefficient in every row inspected, $d, l_1, l_2$ have three coefficients, $x, y$ have four ¹⁸. Five samples and a cubic leave one degree of freedom, so the fit is nearly an interpolation. Publish $\tan f_1, \tan f_2$ at $t_0$.

Step 10. Record metadata. Ephemeris and version, TT reference $t_0$, ΔT, $k_1$, $k_2$, $s_0$, $\pi_0$, ellipsoid, and whether the lunar position is centre of mass ^{6 15}.

A check for the developer: with the NASA 2024 April 8 elements, $x^2 + y^2$ at $t = 0$ gives $\sqrt{0.318244^2 + 0.219764^2} = 0.3868$, and the catalogue's gamma of 0.34314 is the minimum of that distance at 18:18:29 TDT, the catalogued time of greatest eclipse ¹⁸.

Vector shortcut used by the Swiss Ephemeris

The Swiss Ephemeris skips $a, d, \mu$ and works in au with the unit vector $\mathbf{e}$ from Sun to Moon. Its code variable s0 is the almanac's $z$, not the almanac's $s_0$, and it forms the shadow diameters on the fundamental plane directly:

/* distance of moon from fundamental plane */
s0 = -dot_prod(rm, e);
/* distance of shadow axis from geocenter */
r0 = sqrt(dm * dm - s0 * s0);
/* diameter of core shadow on fundamental plane */
d0 = (s0 / dsm * (drad * 2 - dmoon) - dmoon) / cosf1;
/* diameter of half-shadow on fundamental plane */
D0 = (s0 / dsm * (drad * 2 + dmoon) + dmoon) / cosf2;

Here drad is the solar radius, dmoon the lunar diameter, dsm the Sun-Moon distance, and sinf1 = (drad - rmoon)/dsm, sinf2 = (drad + rmoon)/dsm, so f1 is the umbral angle and f2 the penumbral, the reverse of the almanac subscripts ². Algebraically d0/2 equals $c_2\tan f_2$ of the almanac and D0/2 equals $c_1\tan f_1$. The flattening is applied before this by rm[2] /= earthobl; rs[2] /= earthobl; with earthobl = 1 - EARTH_OBLATENESS and EARTH_OBLATENESS = 1/298.25642 ^{2 19}. That trick maps the ellipsoid to a sphere of equatorial radius and is exact for deciding whether the axis meets the surface, but it is not how the almanacs treat the observer, and a developer mixing the two must not apply the flattening twice. The constants are DSUN = 1392000000.0 / AUNIT, DMOON = 3476300.0 / AUNIT, DEARTH = 6378140.0 * 2 / AUNIT, AUNIT = 1.49597870700e+11 ^{2 19}. The library documentation says the northern and southern limits of the umbra and penumbra "are not implemented yet" ¹².

Open-source generators found

The full repository inventory and its comparison matrix are in GitHub repositories. This section keeps only what each repository does about generating elements. Searches were run for "besselian elements" with python, javascript, fortran, rust and C, for Skyfield-based eclipse code, and for the Swiss Ephemeris, Occult, USNO and IMCCE. Findings:

Negative findings: no Fortran or C code that generates elements from a JPL ephemeris with documented validation against the almanac was found. Occult's public tutorial mentions solar eclipses only in its feature list ²³. No IMCCE page describing its own Besselian-element computation was located in the searches run. The USNO eclipse computer describes a topocentric method, not elements ⁹. Yuk Tung Liu's derivation page returned an access error and a JavaScript-only stub, and the Observable notebook by Yamahata returned rate-limit errors on two attempts, so neither was read.

Elements versus direct topocentric computation

The Besselian method answers the global questions cheaply: where the axis meets the ellipsoid, the outline of the umbra, the limits of the penumbra, the duration on the central line, all from six polynomials ²⁴. Espenak's reason for the method is that "using a full ephemeris is more computationally involved", and the elements "are basically an ephemeris valid only over a short period of time (usually about five hours), reduced to a few polynomials" ²⁴. Local circumstances from elements need the observer's $\xi, \eta, \zeta$ and an iteration on the equations of condition, as in sections 8.33 to 8.36 of the 1992 Supplement ⁵.

Direct topocentric computation skips the plane. At each trial time one computes the apparent topocentric positions of the Sun and Moon for the observer, their apparent semidiameters, and the angular separation of their centres, then finds the times when the separation equals the sum or difference of the semidiameters. The Swiss Ephemeris eclipse_how() is a compact example:

rmoon = asin(RMOON / lm[2]) * RADTODEG;
rsun = asin(drad / ls[2]) * RADTODEG;
rsplusrm = rsun + rmoon;
rsminusrm = rsun - rmoon;
...
dctr = acos(swi_dot_prod_unit(x1, x2)) * RADTODEG;
if (dctr < rsminusrm) retc = SE_ECL_ANNULAR;
else if (dctr < fabs(rsminusrm)) retc = SE_ECL_TOTAL;
else if (dctr < rsplusrm) retc = SE_ECL_PARTIAL;
...
attr[0] = lsunleft / rsun / 2;   /* magnitude = (rsun + rmoon - dctr) / (2 rsun) */

with the positions requested after swe_set_topo() for the observer and at TT $= $ UT $+ \Delta T$ ². The USNO eclipse computer does the same: "The computation of Eclipse Local Circumstances is started by iteratively computing topocentric positions of the Sun and Moon to find the time of Maximum Eclipse", using IAU radii of 696,000 km and 1737.4 km and no limb profile ⁹. Quaglia et al. go furthest: "Our computational model does not depend, like the traditional method, on Besselian elements ... to compute nominal predictions and then on applying limb corrections ... Instead, it performs its computations directly with the solar limb and the lunar limb profile to find times of contact and other related quantities. No approximating adjustments are involved" ¹⁰.

When each is used:

• Use elements for catalogues, maps, path limits, and any product that must be reproducible from a published table, and for matching NASA or the Almanac, because their constants and fits are documented ⁶.
• Use direct topocentric computation when the ephemeris is at hand and the question is local: contact times at one site, Baily's beads, flash-spectrum timing, eclipse solar radius. It removes the polynomial truncation, the constant $\tan f$ approximation, the ellipsoid auxiliaries, and the split of $k$ into two values, and it lets a limb profile replace $k$ entirely ¹⁰.
• The two agree to well under a second when the same ephemeris, radii and ΔT are used. The observed differences between predictors come from those inputs, not from the method. Quaglia et al. report Occult's Baily's beads tool giving durations 2 s or more longer than their model at the same site, and Solar Eclipse Maestro 4 s or more, attributed to limb handling ¹⁰.

Sources compared

What a developer should do

1. Build a function elements_at(t_tt) that returns $(x, y, d, \mu, l_1, l_2, \tan f_1, \tan f_2, z)$ from apparent geocentric places of date, following steps 0 to 7 above, with $k_1, k_2, s_0, \pi_0, a_e, \Delta T$ as parameters.
2. Validate it against the 1961 Supplement worked example (1961 February 15 08h ET) ¹ and against the NASA 2024 April 8 polynomial at $t = 0$ ¹⁸, expecting agreement at the $10^{-5}$ Earth-radius level once the same ephemeris, $k$ and ΔT are used. Differences larger than $10^{-4}$ Earth radii in $x, y$ (about 640 m) point to a frame or aberration inconsistency.
3. Only then add the polynomial fit as an export format, and keep the rigorous per-instant function as the computational path for local circumstances.
4. Implement the direct topocentric method as a second, independent path and use the agreement of the two as a regression test.

What this changes

The pipeline gets two computational paths from one ephemeris layer: an elements path for global products and a topocentric path for local products. The ephemeris layer must expose apparent geocentric places of date, apparent topocentric places, and a sidereal-time function that is consistent in its treatment of nutation. Nothing else in the design changes.

Open questions

• Quantify the effect of using astrometric versus apparent places for the Sun and Moon on $x, y$: obtain or write a test comparing both against the NASA 2024 elements with the same DE ephemeris and ΔT.
• Obtain Meeus's Elements of Solar Eclipses 1951-2200 to confirm its sampling scheme and whether it applies nutation to $\mu$.
• Obtain the Python chain behind Eclipse-Engine's data/eclipses.json, or the author's description of it, to see how DE440s was reduced to elements and what $k$ and ΔT were used ⁸.
• Read libephemeris's eclipse module source to see whether it is a port of swecl.c. The provenance document names no file ²⁰.
• Obtain the 1992 Supplement's section 8.35 (general solar eclipse phenomena) in a clean scan to check the equations for the axis-ellipsoid intersection, which the optical character recognition of the available scan garbles.