Wright and Young 2024, the raster method

The question. What does Ernie Wright's published method actually specify, with which datasets and constants, and how much of an implementation can a developer reconstruct from the paper alone?

The publication

The published method is Wright, E. and Young, C. A. 2024, "A Raster-oriented Method for Creating Eclipse Maps", The Astronomical Journal 168:163 (13 pages), published 2024 September 19, DOI 10.3847/1538-3881/ad6b23 ¹. Wright's affiliation is Universities Space Research Association and the Scientific Visualization Studio, Goddard Space Flight Center, Code 606.4. Young is Associate Director for Science, Heliophysics Science Division, Code 670. Wright's ORCID is 0000-0002-5970-5019. The article is open access under the Creative Commons Attribution 4.0 licence and carries an online animation as supplementary material ¹.

Crossref confirms the DOI, the authorship and the issue date ², and the paper is indexed as gold open access under a CC BY licence ¹. NASA's September 2024 announcement ³ and the USRA newsroom release ⁴ both give the title and journal. No preprint exists: there is no arXiv version, and the ADS EPRINT link for bibcode 2024AJ....168..163W does not resolve.

Quotations below are from the publisher's PDF as captured by the Internet Archive on 2024 November 19; the article page itself is not retrievable without a browser session ¹.

What the paper claims

The abstract states the thesis in one sentence: "Ignoring the terrain of both bodies introduces errors on the order of kilometers in the ground track of the umbra and seconds in the duration and contact times of totality" ¹. The paper's own summary of its contribution is that it is "the first detailed description in the eclipse calculation literature" of numerical eclipse simulation, while explicitly not claiming to be the first implementation ¹.

The background section credits prior limb-corrected maps to Michael Zeiler from 2012, "using eclipse software adapted for raster-oriented map calculations by Xavier Jubier and William Kramer", and dates the SVS work to December 2016, when maps and animations of the 2017 August 21 eclipse "for the first time accounted not only for the roughness of the lunar limb but also the terrain of the Earth" ¹.

The raster approach

The method is defined against the Besselian elementsBesselian elementsThe time-dependent quantities $x$, $y$, $d$, $\mu$, $l_1$, $l_2$ and the constants $f_1$, $f_2$ that describe the Moon's shadow relative to the fundamental plane, from which any eclipse circumstance can be computed. tradition. Bessel's method projects the shadow onto the fundamental planefundamental planeThe plane through the Earth's centre perpendicular to the axis of the Moon's shadow. Its x axis lies in the equator pointing east, its y axis points north, and the shadow's cross-section on it is an exact circle. and characterises the eclipse by a few elements and their rates. The paper lists its simplifying assumptions: the Moon is a smooth sphere, its centre of mass coincides with its centre of figurecentre of figureThe centre of the Moon's mean limb figure. The ephemeris tabulates the centre of mass, which lies about 1.9 km away, mostly along the Earth-Moon line, so the sky-plane part is about 0.5 arcseconds, historically applied as +0.50 arcseconds in longitude and −0.25 arcseconds in latitude. A limb profile referenced to the centre of mass, as the LOLA and Kaguya grids are, supersedes the correction., and all observers lie on an ellipsoid ¹.

The raster method instead finds "the presence or absence of the Moon's shadow at every pixel of a raster map over a large number of time steps". Each pixel is coloured "black for total, gray for partial, and white for no eclipse". Repeating over time gives an animation, and integrating the image sequence gives the path of totality ¹. This is how the path limits and the central line are derived: they are not solved for as curves. They emerge as the boundary of the union of the umbra rasters, and the central line is the track of the shadow-axis intersection.

Each pixel is placed in three-dimensional space. Its geodetic longitude and latitude follow from the map projection. Its height above mean sea level is read from an Earth DEM. Height above the ellipsoid is obtained "by adding the geoid height". The paper states that most global DEMs "are now based on the WGS84 ellipsoid and EGM96 geoid" and that the DEM and the coordinate conversion must agree on both ¹. The Earth DEM used for the NASA maps is not named in the paper. The SVS pages name SRTM ⁵.

The terrain effect is quantified as a shift of roughly $h \cot a$, where $h$ is the elevation and $a$ is the Sun's altitude, in the direction of the Sun's azimuth ¹. Figure 1 shows the umbra west of Idaho Falls on 2017 August 21 "perturbed and shifted roughly 1.5 km southeast by the elevation of the Earth's terrain". The 2016 SVS page puts the western-states shift at "as much as 3 kilometers" ⁵.

For pixels near the umbra edge the traditional local-circumstances calculation supplies three values: the distance between the centres of the lunar and solar discs in apparent solar radii, the relative sizes of the discs, and the position angleposition angleAn angle measured on the sky at the centre of a disk, from celestial north through east, used to locate a point on the limb of the Sun or Moon. Contact points and limb heights are indexed by it. of the Moon with respect to the Sun. "To test for totality, the Moon circle is replaced by a limb profile, an array of radii representing the actual shape of the limb. The eclipse is total for a given observer if the limb profile completely encompasses the Sun circle" ¹. The paper stresses that this is a direct comparison on an image plane unique to each observer, whereas the Besselian test of distance from the shadow axis "relies crucially on the assumption that the cross section of the shadow is circular".

Creating the limb profile

Section 4 defines a lunar limb profilelunar limb profileThe height of the Moon's silhouette edge above a reference sphere, tabulated as a function of position angle around the disk for a given libration. It is what turns a smooth-disk eclipse prediction into one that knows where the valleys are. $L$ as "a one-dimensional array of radii at equally spaced angular intervals around the lunar limb". For an observer at distance $d$ from the Moon's centre the construction is ¹:

1. Convert each pixel of the lunar DEM to rectangular coordinates.
2. Apply a rotation matrix encoding the topocentric librationlibrationThe slow apparent rocking of the Moon that changes which part of its surface faces an observer, expressed as a longitude and a latitude of the sub-observer point in selenographic coordinates. It ranges about plus or minus 8 degrees in each coordinate and sets which terrain lies on the limb. at the observer.
3. Convert the transformed position's $(y, z)$ to polar $(r, \theta)$.
4. Replace $r$ with the angular radius $\alpha = \tan^{-1}\!\big(r / (d - x)\big)$.
5. Convert $\theta$ to an index $i$ into the array $L$.
6. If $\alpha$ is greater than $L_i$, set $L_i = \alpha$.

The rotation in step 2 encodes the Moon mean-Earth/polar-axis (ME) frameMoon mean-Earth/polar-axis (ME) frameThe lunar body-fixed frame whose z axis is the mean rotation pole and whose x axis points to the mean sub-Earth point, as realised by a JPL ephemeris. LRO LOLA, SLDEM2015 and Kaguya LALT grids use the ME frame of DE421; the SPICE frame MOON_ME_DE440_ME421 reproduces it from DE440. longitude and latitude of the sub-observer point. After the rotation the centre of the disc as seen by the observer has coordinates $(x, 0, 0)$, "where $x$ is the datum radius (1737.4 km for LRO DEMs) plus the elevation at that point". The arctangent in step 4 performs perspective foreshortening, and the paper notes that far-side points "can and do contribute to the limb". Index 0 is north and $i$ increases counterclockwise ¹.

The datasets are named precisely. "As of this writing, the highest-quality lunar DEMs are SLDEM2015SLDEM2015A merged DEM of LOLA altimetry co-registered with SELENE Terrain Camera stereo, covering 60 S to 60 N at 512 pixels per degree, about 60 m at the equator, with a typical vertical accuracy of 3 to 4 m. (Barker et al. 2016), produced by combining data from the terrain camera (TC) on SELENE and the LRO laser altimeter (LOLA) on LRO, and the polar stereographic LDEMLDEMThe LOLA Lunar Digital Elevation Model series in the LRO GDR archive, gridded in equirectangular projection at 4, 16, 64, 128, 256, 512 and 1024 pixels per degree with heights in metres above a 1737.4 km sphere. from LOLA." SLDEM2015 covers 60 N to 60 S at 512 pixels per degree (60 m at the equator), with subsampled versions at 256 and 128 pixels per degree. The polar LDEMs fill the poles. Both are archived on the LRO node of the Planetary Data System, and the paper gives the link https://pds-geosciences.wustl.edu/missions/lro/lola.htm ^{1 6}. Barker et al. describe SLDEM2015 as built from about 4.5 billion LOLA heights co-registered with 43,200 SELENE TC stereo DEMs, with typical vertical accuracy of 3 to 4 m ⁷.

The coordinate frame is stated: "LRO data are tied to the Moon Mean Earth (Moon ME) frame as realized by the Jet Propulsion Laboratory Developmental Ephemeris 421 (DE421)". The origin of longitude is the mean sub-Earth point, latitudes are planetocentric latitudeplanetocentric latitudeLatitude measured as the angle from the equatorial plane at the body's centre, as used for a spherical datum. Lunar DEMs referenced to the 1737.4 km sphere use it, so no ellipsoidal correction is needed when converting pixels to rectangular coordinates. because the DEM elevations are relative to a sphere, and the axes follow the modern convention with $x$ through longitude 0, $z$ north, $y$ east ¹. The LRO project's coordinate document confirms the ME system, the 1737.4 km reference radius and the selection of DE421 ⁸.

Resolution of the profile: "An L of 18,000 elements has a resolution of 0.02 degrees (roughly 600 m) and pairs well with DEMs of 240 m resolution." Too large an array "may include points not on the limb", too small "will miss small valleys" ¹. The 240 m figure corresponds to the 128 pixel per degree subsample of SLDEM2015, not the full 512 pixel per degree product.

Cost control: creating a profile for every observer "is prohibitively expensive on modest computing hardware. A reasonable compromise is to create a profile for the observer on the shadow axis at time steps of a few minutes, or whenever a libration angle has changed by some small threshold, say 0.01 degrees." Only DEM points within 75 to 105 degrees east and west longitude in SLDEM and a 450 km wide swath through each polar LDEM are processed. Libration in latitude is always small during solar eclipses because they occur near the nodes ¹.

The limb test

Section 5 lists the required quantities ¹:

• $n$, the number of elements in $L$.
• $d_0$, the observer to Moon distance at which $L$ was created.
• $d$, the current observer's Moon distance.
• $c$, the position angle of the lunar axisposition angle of the lunar axisThe angle, measured eastward from celestial north, of the Moon's rotation axis projected on the sky. Limb profiles are indexed from this axis, so it must be computed for the same time and observer as the libration. Written c in Wright and Young's limb test and stored as MoonC in the SVS umbra shapefiles., the position angle of the Moon's north pole axis.
• $\delta$, the Sun to Moon apparent distance in solar radii.
• $\phi$, the position angle of the Sun with respect to the Moon (the more usual angle of the Moon with respect to the Sun differs by $\pi$).
• $r$, the apparent radius of the Sun in radians.

The test examines each $L_i$:

1. $s = (d_0/d)/r$
2. $a = s L_i$
3. $\theta = 2\pi i/n + c$
4. $\rho = a^2 + \delta^2 - 2a\delta\cos(\theta - \phi)$
5. If $\rho < 1$, the eclipse is not total for the observer.

The Sun and Moon are placed in polar coordinates on the observer's image plane, celestial north up, units of apparent solar radius, Moon at the origin and Sun at $(\delta, \phi)$. The scale $s$ corrects for the observer's distance and converts limb radii to solar radii. Step 4 is the polar distance formula with the square root omitted because the comparison is against 1 ¹. The axis angle $c$ and the profile are computed once per time step for the observer on the shadow axis, since differences between nearby observers are far smaller than the profile resolution.

Trivial exclusion defines the limb test zonelimb test zoneIn Wright and Young's raster method, the band of observers near the umbra edge for which totality cannot be decided from the mean lunar radius and every limb element must be tested against the solar disk in units of the apparent solar radius. ¹:

• if $s\max(L) - \delta < 1$, the observer is outside the umbra;
• if $s\min(L) - \delta \ge 1$, the observer is inside the umbra.

Only observers satisfying neither condition need the full test. For annularity the sense reverses, failing on $\rho \ge 1$, with exclusions $s\max(L) + \delta < 1$ (inside the antumbraantumbraThe extension of the Moon's shadow cone beyond its vertex, within which the Moon appears smaller than the Sun. Where the antumbra reaches the ground the eclipse is annular.) and $s\min(L) + \delta \ge 1$ (outside). For eclipses with magnitude near 1 both tests are applied ¹.

antialiasingantialiasingAssigning intermediate values to raster pixels that straddle a boundary, so that the edge is represented at sub-pixel precision rather than as a jagged step. In the SVS method the value comes from the sum of limb-test failures against a threshold set by the pixel area. is folded into the test. For failed elements the sum of $1 - \rho$ is accumulated. If the sum is greater than 0 but less than a threshold $\epsilon$, the pixel gets the intermediate value $1 - \text{sum}/\epsilon$. "For $n$ = 18,000, setting $\epsilon = A(p)/10$ works well, where $p$ is the pixel at the center of the shadow and $A$ is the area in square kilometers" ¹. The stated purpose is to remove false precision and to help contour-tracing algorithms, which is how the vector shapefiles are extracted from the rasters.

Contact times at a site

The paper does not describe a separate local-circumstances routine. Second and third contact at a site follow from the same test applied to that site across time steps: totality begins at the first step where the limb profile encompasses the Sun and ends at the last. The cities JSON released with the 2024 map gives contact times to the second for towns in totality and appears to be produced this way ⁹. The time resolution of the released umbra polygons, 1 s for the high-resolution set, bounds the timing precision of anything derived from them ¹⁰.

The polygon explained

Section 6.2 uses the antialiasing pass to identify, for each edge pixel, "which $L_i$ is at the greatest internal distance from the solar limb". Plotting those elements "form[s] a minimal and coherent set of umbra polygon edges". Each limb point is then treated as a pinhole projecting an image of the Sun onto the Earth. "The locus of all such rays forms a double cone with its apex at the limb point." The Sun-image cone angle equals the shadow-cone angle and its axis deviates from the shadow axis by 2.4 arcseconds, the angle subtended by the Moon's radius at the Sun ¹.

The shape follows: "The outer edge of the penumbra is defined by the outer edges of the Sun images. The umbra is the set of interior points not reached by any of these Sun images." At 18:00 UTC on 2017 August 21, 49 valleys contributed edges, and Figure 12 shows the 49 projected Sun images with the umbra as the hole in the middle. The edges are arcs, not straight lines. For the antumbra the focal points are peaks and the shape is the intersection of the Sun images ¹. The 2024 SVS explainer page reproduces this account and the 49-valley count ¹¹. Observers on an edge see the diamond ring from that valley, and observers at a vertex see a double diamond ring, a point the paper credits to Zeiler ¹².

Broken annular

Section 6.3 introduces a "third type of central eclipse in which the shadow axis intersects the Earth's surface but the eclipse at that location is neither total nor annular". In the smooth-limb model the umbra and antumbra cones meet at a point, so a hybrid switches instantly. In reality there is a gap in which "a very thin ring of sunlight around the Moon is interrupted by high points on the lunar limb" ¹.

The proposed criterion is that broken annularbroken annularA central eclipse phase, occurring in the gap between the umbra and the antumbra of a hybrid eclipse, in which the Moon is too large to fit inside the Sun and too small to cover it, so a thin ring of sunlight is interrupted by lunar peaks. Wright and Young map it where the range of limb radii exceeds the Sun to Moon separation. locations are the locus of points for which

valid only at locations that have already failed both the totality and annularity tests ¹. Three examples are given. The 1894 April 6 eclipse "is never actually total, spending 2 hr as a broken annular over a span of nearly 5400 km". The 2005 April 8 hybrid has totality "flanked by more than 1300 and 730 km of broken annularity". The 2013 November 3 hybrid "has no true annular phase". Figure 13 maps the 2005 path assuming $r$ = 696,000 km and simulates Baily's beads at 1 s intervals from 21:55:26 to 21:55:32 UTC as seen from 94.02587 W, 6.45677 N. The SVS release of that simulation runs in real time from 21:55:20.5 to 21:55:35.5 UTC ¹³.

The solar radius

Section 6.4 is the paper's statement on the constant that matters most near the limits. "Nearly all eclipse calculations use a de facto standard value of the solar radius first published more than a century ago and expressed as either 959.63 arcseconds at 1 au (Auwers 1891) or 696,000 km (Yatskiv 1978)." The IAU 2015 nominal unit is 695,700 km, equivalent to 959.22 arcseconds, and the paper quotes Prsa et al. that nominal values "should be understood as conversion factors only" ¹.

The paper argues the nominal solar radiusnominal solar radiusThe IAU 2015 conversion constant of exactly 695,700 km (959.22 arcseconds at 1 au), chosen to match the helioseismic photosphere. It is a unit for expressing stellar radii, not a measured limb, and is too small for eclipse geometry. "does not provide the solar radius required for accurate eclipse prediction, which defines totality as the complete extinction of the photosphere". It lists recent eclipse solar radiuseclipse solar radiusThe solar radius that, used in an accurate eclipse model with a lunar limb profile, reproduces observed second and third contacts. Flash-spectrum and photometer measurements put it at 959.95 to 960.0 arcseconds at 1 au, about 0.3 arcseconds larger than the canonical value. determinations: 959.99 plus or minus 0.06 arcseconds (Lamy et al. 2015), 959.95 plus or minus 0.05 (Quaglia et al. 2021), 959.98 (Jubier et al. 2021) and 960.01 (Guhl 2023) ¹. Quaglia et al. obtained their value from flash-spectrum video at the southern limit of the 2017 path and report no significant wavelength dependence ¹⁴. Guhl's value comes from Baily's beads photometry at the northern limit of the 2023 April 20 hybrid, a correction of plus 0.38 arcseconds to 959.63, giving 960.01 plus or minus 0.12 arcseconds ¹⁵.

The sensitivity statement is quotable: "a discrepancy of 1 s in duration, for example, corresponds to an error as small as 20 km (0.03 arcseconds) in solar radius" for points very close to the path limits ¹. The paper recounts the March 2024 episode in which John Irwin's map, "based on a solar radius of 959.95 arcseconds", shifted the northern limit "to the southeast by several city blocks compared to maps calculated with the conventional radius", and dozens of journalists asked whether NASA would revise its map ¹. It did not. See svs-vs-other-predictors.md.

The paper does not say which solar radius the 2024 NASA map used. The SVS product pages for 2017 state 696,000 km ¹⁶, and the paper's Figure 13 uses the same value, so the working assumption is that all SVS maps to date use 696,000 km, which is 959.63 to 959.645 arcseconds depending on the au conversion used.

Appendix A, the software stack

The appendix is the only public description of the implementation. "The code that generated the NASA maps relies primarily on SPICE." It prints C-language calls ¹:

vsub_c( obs, mpos, pos );          /* translate to Moon-centered */
mxv_c( mmat, pos, pos );           /* rotate from ICRF to Moon ME */
reclat_c( pos, dist, lam, bet );   /* to spherical coords */

which return the observer's Moon distance and the Moon ME sub-observer longitude and latitude, that is, the topocentric libration. The Moon's ICRF position comes from spkezp_c(ID_MOON, t, "ICRF", "lt+s", ID_EARTH, mpos, &lt) with light time and stellar aberration, from the kernel de440.bsp. The rotation to the body-fixed frame is pxform_c("ICRF", "MOON_ME_DE440_ME421", t, mmat) using moon_pa_de440_200625.bpc and moon_de440_200625.tf. That frame kernel is the 2020 release, and NAIF's current file for the same frame is moon_de440_250416.tf. The frame name records that "the origin of selenographic coordinates for LRO data products is the mean sub-Earth point found by the DE421 calculation, not that of DE440" ¹.

The observer is placed with georec_c(lon, lat, h, earth_r, flattening, obs), where $h$ is ellipsoid height, and rotated to the ICRF with pxform_c("ITRF93", "ICRF", t, omat) using earth_latest_high_prec.bpc. The paper warns that "SPICE Earth orientation kernels are less reliable or nonexistent for times in the future or before 1962" and gives the SOFA equivalent, iauC2t06a(tta, ttb, uta, utb, xp, yp, omat_inv) followed by a transpose ¹. This is where ΔTΔTThe difference TT − UT1 between uniform Terrestrial Time, in which the elements are computed, and Universal Time, which follows the Earth's irregular rotation. About 69 s in 2024. It converts the hour angle μ to a geographic longitude, is known only by prediction for future eclipses, and is the largest error source for historical ones. enters: it is implicit in the Earth orientation parameters (EOP)Earth orientation parameters (EOP)The measured quantities UT1 − UTC, polar motion (x, y), length of day and celestial pole offsets, published by the IERS, that relate the Earth-fixed frame to the celestial frame. They cannot be known in advance and are the only prediction inputs that must be re-fetched before an eclipse. kernel rather than a stated number. The SVS 2017 pages state the value that resulted, $\Delta T$ = 68.917 s, and describe the kernel as "earth_070425_370426_predict.bpc (ΔT corrected)" ¹⁶.

The ephemeris is therefore DE440 at the time of writing the paper, while every SVS product page from 2015 through 2025 lists DE421 as its dataset ^{10 17}. Both cannot be true for the same run. The likely reading is that the maps were made with DE421 and the paper's appendix was written against the current kernel set. The difference between DE421 and DE440 in the Moon's geocentric position in 2024 is far below the 600 m limb-profile resolution, so the choice has no visible effect on the products.

Nothing in the paper or on any SVS page mentions Maya, RenderMan or any commercial renderer for the eclipse calculation itself. The 2017 shapefile page says only that "the map was rendered in animation software, but maps are more typically created using GIS tools" ¹⁸. The 2015 page says the Besselian method "was adapted to the routines available in NAIF's SPICE software library" ¹⁹.

What the paper does not say

• No $\Delta T$ value for 2024. It is embedded in the Earth orientation kernel.
• No time step for the 2024 map. The released umbra_hi shapefile is at 1 s and umbra_lo at 10 s ¹⁰.
• No refraction model. The March 2025 map page says the sunrise and sunset lines "are based on an idealized model" and that "the apparent positions of the Sun and Moon relative to the horizon are affected by local terrain and atmospheric refraction" ¹⁷.
• No Earth DEM named. SVS pages credit SRTM ⁵.
• No numerical comparison with Espenak or Jubier. The only comparison is the qualitative "several city blocks" against Irwin's map.
• No stated uncertainty in kilometres for its own path edges. NASA's public statement during the 2024 dispute acknowledged "a tiny but real uncertainty about the size of the Sun" and that "uncertainty in the Earth's rotation can also affect eclipse predictions on this level" ²⁰.
• No code release. See the searches described in svs-products-and-data.md.

Prior art the paper builds on

The reference list is the best map of the field's genealogy. On the limb: Kopal 1965 and Weimer 1979 for early catalogues, then the Watts chartsWatts chartsThe lunar limb profiles published by Chester Burleigh Watts in 1963 from machine tracings of about 700 photographs covering the range of librations, used for eclipse and occultation timing until laser altimetry replaced them. of 1963, which "remained in use for eclipse and occultation timing for half a century" despite flaws documented by Scott 1988, before being supplanted by SELENE and LRO altimetry ^{1 21 22 23 24}. On applying limb profiles: Herald 1983 in the Journal of the BAA, "a purely graphical construction that involves plotting solar limb curves on transparencies" over exaggerated limb plots like Duncombe's 1973 USNO circular, sliding along the position angle to find the limb-corrected contacts ^{25 26}. The paper says software authors "have automated this method and adapted it to newer lunar terrain data, but the details of their implementations remain largely undocumented". It points to Herald's Occult documentation and to Wang et al. 2021 for details of limb-profile use ^{27 28}.

On elevation: Williams 1971 for adapting Bessel by replacing 1 Earth radius with the elevated radius, Lewis 1940 for USNO central-line maps at ionospheric altitudes, Fiala et al. 1987 for interpolating between tabulated heights, and Schneider 2004 for aircraft intercepts ^{29 30}. On the canonical formulation: Bessel 1829 and 1830, Chauvenet 1863, both editions of the Explanatory Supplement, Meeus 1989 and Montenbruck and Pfleger 2000 ³¹. On ephemerides and frames: Folkner et al. 2009 for DE421, Park et al. 2021 for DE440, the LRO coordinate-system memo of 2008 ^{32 33 8}.

Two names a reader of the eclipse-radius literature would expect are absent. Sigismondi is not cited. Dunham is not cited by name, though IOTA appears through Herald and through Guhl's paper in the Journal for Occultation Astronomy. The Kaguya (SELENE) work enters only as Araki et al. 2009 and through SLDEM2015's use of the SELENE terrain camera.

Sources compared

What a developer should do

Read the paper's sections 4 and 5 first. They are complete enough to implement. Then:

1. Download SLDEM2015 at 128 or 256 pixels per degree and the polar LDEMs from the PDS LOLA node ⁶. Read dsmap.cat and dsmap_polar.cat for the pixel-to-longitude mapping. Use the 1737.4 km datum and planetocentric latitude.
2. Use SPICE with de440.bsp, moon_pa_de440_200625.bpc, the frame kernel moon_de440_250416.tf, which is NAIF's current release of the moon_de440_200625.tf the paper prints, and the frame MOON_ME_DE440_ME421, exactly as printed in Appendix A, so that the DEM and the frame agree ¹.
3. Build $L$ with 18,000 elements for the shadow-axis observer, rebuilt when libration moves by 0.01 degrees. Store $\min(L)$ and $\max(L)$.
4. Implement the limb test as written, with the trivial exclusions, and the antialiasing with $\epsilon = A(p)/10$.
5. Convert SRTM or a successor DEM's orthometric heights to ellipsoid heights with EGM96 before calling georec_c.
6. Choose and state the solar radius. The paper gives you the sensitivity: 0.03 arcseconds is 1 s of duration at the limits.
7. Validate against the released umbra_hi polygons for 2024 ³⁴, which are the paper's own output.

What this changes

For the pipeline design this replaces the "limb correction as a post-processing step on Besselian contacts" model with a per-pixel test. The path limits, central line, duration contours and obscuration contours are all derived products of a raster stack, which means one code path produces every map layer. It also establishes that the limb profile, not the ephemeris, sets the resolution floor at about 600 m, and that the solar radius, not the limb, is the dominant unresolved uncertainty at the limits.

Open questions

• The 2024 run's $\Delta T$ and ephemeris. Obtain the SVS API JSON for page 5123 or an author statement to settle whether DE421 or DE440 and which EOP kernel produced the 2024 shapefiles.
• The Earth DEM resolution used. SRTM exists at 1 and 3 arcseconds. The paper's 240 m lunar pairing suggests a coarse Earth grid, but the map pixel size for 2024 is unstated. The insolation dataset is at 360/8192 degrees per pixel ³⁵, which may be the shadow-raster grid.
• A test of the broken-annular criterion against observation. The 1986 October 3 and 2023 April 20 eclipses have observers. Compare the criterion's locus with the reported bead sequences.
• The unpublished code. Ask the author whether the C code can be released. Until then the shapefiles are the only ground truth.