A Raster-oriented Method for Creating Eclipse Maps

Because the Sun is an extended light source, every sunlit object casts shadows that can be divided into areas where the Sun is completely blocked (the umbra) and where it is only partially blocked (the penumbra). The boundary between these areas is almost never apparent. The edge of the umbra shades imperceptibly into the adjacent penumbra, and we typically only notice that shadows in sunlight are somewhat blurry. We therefore lack intuition about what the boundary looks like.

For eclipse observers, however, this boundary is the threshold dividing profoundly different experiences (Aveni 2017; Russo 2017). And because of the Moon's rugged terrain, the shape of the umbra's edge is nothing like we might imagine. It is not an oval, as it would be if the Moon were a perfectly smooth sphere. In fact, it is not even approximately oval in some extreme cases. Nor is it simply a projection of the Moon's irregular silhouette, as it would be if the Sun were a point source. The true shape of the umbra is polygonal (Figure 1).

Zoom In Zoom Out Reset image size

The innovations in eclipse calculation that eventually revealed the shape of the umbra have largely been driven by software for eclipse photography and stellar occultation timing, in which errors of kilometers and seconds are increasingly intolerable. The authors of this software refine their predictions by accounting for the irregular shape of the lunar limb. Beginning in 2012, cartographer Michael Zeiler published eclipse maps featuring limb-corrected umbra shapes (Zeiler 2012; Zeiler & Bakich 2020), using eclipse software adapted for raster-oriented map calculations by Xavier Jubier (Jubier 2022) and William Kramer.

In December of 2016, we published maps and animations of the path of the 2017 August 21 total solar eclipse that for the first time accounted not only for the roughness of the lunar limb but also the terrain of the Earth (Wright 2016a). The methods used differ significantly from traditional eclipse-mapping approaches, relying instead on experience with NASA raster data sets and raster-oriented techniques more commonly used in three-dimensional animation and visual effects. The resulting improvement in accuracy sheds light on the true shape of the eclipse shadow, the nature of annular-total ("hybrid") eclipses, and the precise measurement of the solar photosphere.

Friedrich Wilhelm Bessel's formulation has been the standard method of eclipse calculation for two centuries (Bessel 1829a, 1829b, 1830; English translation, Bessel 1842; Hansen 1859; Chauvenet 1863; Newcomb 1882; Buchanan 1904; Comrie 1933; Wilkins & Springett 1961; Meeus1989; Montenbruck & Pfleger 2000; Urban & Seidelmann 2012). The essence of the method is the projection of the Moon's shadow onto a plane perpendicular to the shadow axis (the line through the centers of the Sun and Moon). Reduced to a circle moving on a plane (Figure 2), an eclipse shadow's path can be completely characterized by a small set of parameters, called Besselian elements, along with their time rate of change.

Zoom In Zoom Out Reset image size

Bessel's method is analytically rigorous and tractable without modern computing hardware. It made practical the production of eclipse maps in the annual international astronomical almanacs as well as large canons of eclipses (Oppolzer 1887; Meeus et al. 1966; Espenak & Meeus 2006; Espenak 2019). But it relies on a number of simplifying assumptions, for example that the Moon is a smooth sphere, that its center of mass (the position listed in ephemerides) and its center of figure (the center of the visible disk) are in the same place, and that all observers are on the surface of an ellipsoid, none of which are strictly true and each of which introduces small but significant errors.

Mountains and valleys in the lunar terrain make the limb of the Moon craggy and irregular (Figure 3), and this affects second and third contact times (the times at which totality begins and ends) or, equivalently, the shape of the umbra on the ground at any given time. Valleys along the limb admit sunlight at times when a smooth limb would eclipse it, and mountains block sunlight at times a smooth limb would not. The lunar limb varies over time due to libration, which rotates different parts of the lunar terrain out to the limb, changing its shape (Figure 4).

Zoom In Zoom Out Reset image size

Several attempts were made in the first half of the 20th century to catalog the Moon's limb (Kopal 1965; Weimer 1979). By far the most influential and comprehensive were the charts published by Chester Burleigh Watts (Watts 1963; Dick 2002; Herald 2006). Although flawed (Scott 1988), the Watts limb profiles remained in use for eclipse and occultation timing for half a century before finally being supplanted by global laser altimetry data from lunar orbiting spacecraft: JAXA's SELENE (also known as Kaguya) mission (Araki et al. 2009) and NASA's Lunar Reconnaissance Orbiter (LRO; Smith et al. 2017).

Herald (1983) describes a method of using the Watts limb profiles to adjust second and third contact times as well as the northern and southern limits of the path of totality. In its original form, Herald's approach is a purely graphical construction that involves plotting solar limb curves on transparencies. These are overlaid on exaggerated lunar limb plots like those of Julena Duncombe (1973) and shifted along the Sun's position angle, in effect sliding backwards and forwards in time to find the limb-corrected contacts. Eclipse software authors have automated this method and adapted it to newer lunar terrain data, but the details of their implementations remain largely undocumented, and the results are more commonly limited to refinements of contact times, not the production of more accurate maps.

The effect of the observer's elevation on the local circumstances of an eclipse has received much less attention. Elevation can shift the shadow track by several kilometers. This happens because, in general, the shadow cone is not perpendicular to the Earth's surface (Figure 5). Higher elevations intercept the slanted shadow cone at an offset in the direction of the Sun's azimuth. The magnitude of the shift is roughly $h\cot a$, where h is the elevation and a is the Sun's altitude. The effect is therefore most pronounced at high elevations and when the Sun and Moon are near the horizon. Where the Sun's azimuth is parallel to the eclipse track, the shift is along-track rather than across-track, a shift in time rather than space.

Zoom In Zoom Out Reset image size

Bessel's formulation can be adapted to account for elevation, typically by replacing the value of 1 Earth radius implicit in many of the formulas with the elevated radius, still expressed in units of Earth radii (Williams 1971). Maps of the eclipse central line at altitudes within the ionosphere were made in this way by the U.S. Naval Observatory as early as 1940 (Lewis 1940). For estimates in the field, Fiala et al. (1987) suggests interpolating between calculations tabulated at two or three heights above mean sea level. Schneider (2004) has written software for plotting aircraft flight paths that intercept the shadow cone and has used it for several airborne eclipse observations. Numerous applications and websites include elevation in the definition of place for the calculation of local circumstances.

But none of this work was used to map the effect of terrain on the shape of the umbra outline, nor did it need to confront any of the subtleties of geodesy. Mean sea level, the most common reference surface for terrestrial elevation, is not an ellipsoid, and the values in a digital elevation model (DEM) are only meaningful in the context of its datum: While most global DEMs are now based on the WGS84 ellipsoid and EGM96 geoid, many local and regional maps are not.

In digital cartography as in computer graphics more generally, spatial data and imagery can be divided into vector-based and raster-based representations. Vector-based images are space-bounding, using lines and curves, while raster-based images are space-filling, using pixels; these are logical duals (Peuquet 1984).

A purely Besselian approach to eclipse calculation produces results that are vector oriented. The solutions of equations in Besselian elements generate curves—path limits, sunrise and sunset curves, isochrons, and so on—that naturally lend themselves to vector-based maps. Existing refinements that account for the irregular shape of the Moon and Earth implicitly attempt to preserve this vector orientation, but perhaps uncomfortably: They are given short shrift in even the most comprehensive descriptions of eclipse calculations, and they are only infrequently used to make more accurate maps.

The fully raster-oriented approach described here is fundamentally based on finding the presence or absence of the Moon's shadow at every pixel of a raster map over a large number of time steps. In computer graphics terms, this is akin to rendering, the pixel-by-pixel process used by most computer animation software to generate synthetic images. All of the curves drawn on traditional eclipse maps are readily derived from these rendered shadow images. Pixel calculations are highly parallelizable, and they obviate the need to fit refinements into a vector-oriented scheme.

For a given time, each pixel of the map image is colored according to the eclipse condition at that location—black for total, gray for partial, and white for no eclipse, for example. The resulting image is a map of the shadow on the Earth's surface at that time. Repeating this calculation over a range of times produces a sequence of images that can be compiled into an animation of the shadow's motion. The image sequence can also be integrated to create a map of the shadow's path (Figure 6).

Zoom In Zoom Out Reset image size

This approach lends itself naturally to the precise location of every map pixel in three-dimensional space. Each map location's geodetic longitude and latitude are derived from its pixel coordinates in the shadow map according to the map projection. Its height above mean sea level is read from a DEM of the Earth. Height above sea level is translated to height above a reference ellipsoid by adding the geoid height, which is implicit in the DEM. Transformation from a point on the ellipsoid to other frames is then mathematically straightforward.

Map pixels found to be near the edge of the umbra are given special attention to account for the effect of the irregular lunar limb on the shape of the umbra. The traditional local circumstances calculation can supply three especially relevant values: the distance between the centers of the lunar and solar disks in units of apparent solar radii, the relative sizes of the disks, and the position angle of the Moon with respect to the Sun. In effect, this projects circles representing the Sun and Moon onto an image plane facing the observer. 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.

This is a direct comparison of the images of the Sun and Moon projected onto an image plane unique to each observer, in contrast to the traditional test of the observer's distance from the shadow axis, which relies crucially on the assumption that the cross section of the shadow is circular. The true shape of the cross section is complex, as explained later, and is not known a priori. It emerges from the calculation for all observers in the shadow.

Most aspects of a raster-oriented approach are easily adapted from the traditional methods. The following two sections provide details about the novel and less familiar components, the creation (Section 4) and use (Section 5) of a lunar limb profile. See also Wang et al. (2021) and the documentation of D. Herald's Occult software (Herald 2021).

A limb profile is a one-dimensional array of radii at equally spaced angular intervals around the lunar limb. For a given observer at a distance d from the Moon's center, a limb profile L can be constructed in the following way.

• 1.  
  Convert each pixel of the lunar DEM to a position in rectangular coordinates.
• 2.  
  Apply a rotation matrix encoding the topocentric libration at the observer.
• 3.  
  Convert the transformed position's (y, z) to polar (r, θ).
• 4.  
  Replace r with the angular radius $\alpha ={\tan }^{-1}(r/(d-x))$.
• 5.  
  Convert θ to an index i into the array L .
• 6.  
  If α is greater than L _i , set L _i = α.

As of this writing, the highest-quality lunar DEMs are SLDEM2015 (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 LDEM from LOLA. SLDEM is a gridded product in cylindrical projection covering latitudes between 60°N and 60°S at a resolution of 512 pixels per degree (60 m at the equator). Versions subsampled to 256 and 128 pixels per degree are also available. Polar LDEM maps are available at a range of resolutions.

SLDEM and LDEM complement each other. Elevations derived from TC imagery at 10 m per pixel were coregistered with LOLA data to fill between-track gaps in LOLA's coverage at lower latitudes, while the higher density of LOLA samples in polar LDEM overcomes the lack of shadow-free imagery from TC. Both products benefit from improved orbit determinations made possible by high-resolution gravity models from the GRAIL mission (Lemoine et al. 2014), and both are archived on the LRO node of the Planetary Data System (Neumann 2010). ⁴

For these DEMs, the transformation from pixel coordinates to longitude and latitude is defined in LRO archive files named dsmap.cat and dsmap_polar.cat. LRO data are tied to the Moon Mean Earth (Moon ME) frame as realized by the Jet Propulsion Laboratory Developmental Ephemeris 421 (DE421; LRO 2008; Folkner et al. 2009). Moon ME's latitude/longitude origin is the mean sub-Earth point. Because LDEM and SLDEM elevations are relative to a spherical datum, the latitude is planetocentric and the conversion to rectangular coordinates in step (1) is straightforward. Conventionally, the x-axis is in the equatorial plane through longitude 0, the z-axis points north, and y-axis points east (in the modern Mare Crisium sense), with the origin at the Moon's center of mass.

The angles encoded in the rotation matrix in step (2) are the Moon ME longitude and latitude of the sub-observer point, the place on the Moon where the observer is at the zenith. After step (2), the point at the center of the lunar disk as seen by the observer should have coordinates (x, 0, 0), where x is the datum radius (1737.4 km for LRO DEMs) plus the elevation at that point. The value r found in step (3) has the linear units chosen in step (1), typically kilometers. This is inconvenient for comparison with the apparent size of the Sun, so step (4) converts r to the subtended angle α, the angular offset from the center of the lunar disk. The arguments to the arctangent also perform perspective foreshortening: Of two points with equal r, the one closer to the observer will appear larger and therefore have the larger α. (Note that points on the far side from the observer can and do contribute to the limb, as long as their r is large enough to overcome the modest effect of their greater distance from the observer.) Finally, the index i in step (5) should be 0 for north and increase counterclockwise.

After applying this procedure to every point in the DEM, each L _i contains the largest α in the ith angle bin, and L as a whole comprises the shape of the limb (Figure 7). L should contain enough elements to represent the limb with high fidelity, but not more than can be supported by the DEM's resolution. An L that is too large may include points not on the limb, and one that is too small will miss small valleys. An L of 18,000 elements has a resolution of 0.02° (roughly 600 m) and pairs well with DEMs of 240 m resolution.

Zoom In Zoom Out Reset image size

Every observer sees a unique limb profile, but the differences for observers near each other are vanishingly small relative to the finite resolution of the limb profile, and in any case creating a limb 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°. It is also sufficient to process only those points within 75°–105° east and west longitude in SLDEM and a swath 450 km wide through the north and south polar stereographic LDEM. No points outside this band will contribute to the profile for solar eclipses (Figure 8).

Zoom In Zoom Out Reset image size

An observer found to be in or near the edge of the umbra is tested against the limb profile. The test requires the following quantities, many of which are known or easily found from the traditional calculation:

• n, the number of elements in L .
• d₀, the observer–Moon distance at which L was created.
• d, the current observer's Moon distance.
• c, the position angle of the Moon's north pole axis.
• δ, the Sun–Moon apparent distance in solar radii.
• ϕ, the position angle of the Sun with respect to the Moon.
• r_☉, the apparent radius of the Sun in radians.

Position angle ϕ is the reflection of the more commonly calculated position angle of the Moon with respect to the Sun. The two angles differ by π. Axis angle c is the apparent rotation of the Moon around the center of its disk, or the angle between the axis (its projection on the sky) and celestial north. Combined with the libration angles, it completes the description of the Moon's orientation relative to the observer. This is fairly inexpensive to calculate, but differences between observers near each other are extremely small, and in any case differences much smaller than the resolution of L will not change the result of the test. It is reasonable to calculate c once per time step for just the observer on the shadow axis.

The limb test examines each L _i to see whether it lies inside or outside the solar disk:

• 1.  
  s = (d₀/d)/r_☉.
• 2.  
  a = s L _i .
• 3.  
  θ = 2π i/n + c.
• 4.  
  $\rho ={a}^{2}+{\delta }^{2}-2a\delta \cos (\theta -\phi )$.
• 5.  
  If ρ < 1, then the eclipse is not total for the observer.

The Sun and Moon are located on the observer's image plane in polar coordinates, with celestial north up and in units of solar radii. The Moon is centered at the origin, and the Sun's center is (δ, ϕ). The limb scale s adjusts for the current observer's Moon distance and converts limb elements to units of apparent solar radius. The position of limb element L _i is found by the next two steps to be (a, θ). Step (4) is just the polar formula for the distance between points (δ, ϕ) and (a, θ), although the square root is omitted, since the comparison is with 1.0 (solar radius).

This test is relatively expensive. To limit its use to observers near the umbra edge, observers can be trivially excluded using the following conditions:

• if $s\max ({\boldsymbol{L}})-\delta \lt 1$, the observer is outside the umbra; or
• if $s\min ({\boldsymbol{L}})-\delta \geqslant 1$, the observer is inside the umbra.

This just says that if the smallest limb radius (minus the Sun–Moon distance) encompasses the Sun, the observer must be in totality, while if the largest radius is inside the solar disk, the observer cannot be. The observer is in the limb test zone (must be tested against the entire limb profile) only if neither of these conditions are true (Figure 9). The values of $\min ({\boldsymbol{L}})$ and $\max ({\boldsymbol{L}})$ are found and stored when the limb profile is created.

Zoom In Zoom Out Reset image size

To test for annularity, the sense of the comparison in step (5) is reversed. Instead of failing on ρ < 1, the test fails on ρ ≥ 1. Observers can be trivially excluded from the annularity limb test zone using the following:

• if $s\max ({\boldsymbol{L}})+\delta \lt 1$, the observer is inside the antumbra; or
• if $s\min ({\boldsymbol{L}})+\delta \geqslant 1$, the observer is outside the antumbra.

If the eclipse magnitude is very close to 1, then at any given point along the path it might be annular, total, or in an intermediate state discussed later. For those eclipses, both tests must be applied.

It is also possible to introduce antialiasing in the test without costly subpixel sampling. Consider a grayscale raster that serves as an alpha mask containing just the image of the umbra. If pixels are set to only one of two values, 1 (white) inside the umbra and 0 (black) outside, for example, the edges of the umbra will be jagged. To mitigate this, a small threshold value, , is introduced. If the sum of the failed ρ differences (that is, 1 − ρ for all ρ < 1) is greater than 0 but less than , the pixel can be set to an intermediate value, 1 − sum/.

Antialiasing may naively look like blurring, but its purpose is to remove false precision—to attenuate frequencies higher than the map's pixel spacing. Antialiasing also enhances the esthetics of the map and improves the performance of algorithms that infer contours from rasters. The choice of a value for depends primarily on the resolutions of the map and the limb profile. The goal is a set of edge pixels, 1 or 2 pixels wide, at the boundary of the umbra image. For n = 18,000, setting = A(p)/10 works well, where p is the pixel at the center of the shadow and A is the area in square kilometers.

6.1. The 2017 Eclipse Map

For NASA, the 2017 eclipse was a high-priority event for public engagement and education. It involved the communications teams for missions (especially LRO), divisions (especially NASA heliophysics), and NASA headquarters. The engagement effort was spearheaded by the Heliophysics Education Consortium, a NASA Science Mission Directorate Science Activation team. By an order of magnitude, this was the largest webcast, website, and social media event supported by NASA, ultimately reaching over 4 billion people (Young et al. 2019).

The map (Figure 10) was a key component of that effort. It served as an authoritative NASA source of information about the eclipse and increased public awareness of the NASA data sets it used. Its accuracy and method of construction attracted significant media attention and became the subject of a popular NASA video (Wright & Duberstein 2016). The NASA 2017 eclipse website (Young & Patel 2017) was the most visited U.S. government site over the period 2017 May 1 to August 30, logging 5.2 billion hits and 700 terabytes served to 50 million visitors—and the map was the most downloaded item. The map was visited an additional 4.6 million times at the NASA Scientific Visualization Studio website (Wright 2016b).

Zoom In Zoom Out Reset image size

6.2. The Polygon Explained

As a side effect of antialiasing, it is possible to identify map pixels containing the umbra edge and to verify that each edge of the umbra polygon is associated with a single or dominant lunar valley on the limb (Figure 11), either the last point still admitting sunlight before the start of totality or the first point at the end. In both cases, this is the point responsible for the diamond-ring effect for observers located on that edge. Observers near a vertex experience a double diamond ring (Zeiler 2024).

Zoom In Zoom Out Reset image size

Map pixels with fractional umbra coverage lie on or adjacent to the umbra outline. For each such pixel, the limb test code can report which L _i is at the greatest internal distance from the solar limb. This identifies the lunar valley either solely or predominantly responsible for the umbra edge at that pixel. When the reported elements of L for all such pixels are plotted, they form a minimal and coherent set of umbra polygon edges, as shown in Figure 11.

To understand how the polygonal shape arises, each of these limb points can be treated as a geometric focal point analogous to the aperture of a pinhole camera projecting an image of the Sun onto the surface of the Earth. Mapping the projected Sun images reveals a flower-like pattern with a lacuna in the center: the umbra (Figure 12).

Zoom In Zoom Out Reset image size

Consider rays fired from the Sun's limb through a point on the lunar limb and intersecting the Earth. The locus of all such rays forms a double cone with its apex at the limb point. The sunward cone is tangent to the solar photosphere, and the earthward cone delimits the projected image of the Sun. The Sun image cone angle is equal to the shadow cone angle, and the axis is very nearly parallel to the shadow axis, deviating by ∼24, the angle subtended by the Moon's radius at the Sun.

As Figure 12 suggests, the shape of the shadow, both umbra and penumbra, is delimited by the superposition of Sun images projected from every point on the lunar limb. 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. Only a small number of Sun images associated with deep valleys on the lunar limb (49 of them in this case) actually contribute edges to the umbra outline. Note that the umbra polygon's edges are not straight lines, even aside from their projection onto the curved surface of the Earth. Each polygon edge is an arc from the edge of a Sun image.

The same construction produces the shape of the antumbra during annular eclipses, but with two important differences. The focal points are peaks rather than valleys, and the antumbra shape is formed by the intersection of the projected Sun images.

6.3. Broken Annular

A consequence of true-limb eclipse calculations is greater recognition of 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 geometry of classical, smooth-limb eclipse calculations, the umbra and antumbra form a double cone, the apices meeting at a point. For hybrid eclipses, this implies that as the eclipse magnitude moves between values greater than and less than 1, the eclipse transitions instantly between total and annular phases. No intermediate state is possible.

In actuality, there is a sizable gap between the umbra and antumbra. When the intersection of the shadow axis and the Earth falls in this gap, observers at the intersection see a broken annular eclipse in which a very thin ring of sunlight around the Moon is interrupted by high points on the lunar limb. The Moon is too large to fit entirely inside the solar disk but too small to entirely cover it, an impossibility for a smoothly spherical Moon.

The broken annular phase connecting annular and total portions of hybrid eclipses can extend for thousands of kilometers. The nominally hybrid eclipse of 1894 April 6 is never actually total, spending 2 hr as a broken annular over a span of nearly 5400 km. The totality in the middle of the 2005 April 8 hybrid is flanked by more than 1300 and 730 km of broken annularity. And the hybrid of 2013 November 3 has no true annular phase; it begins and ends as a broken annular lasting less than a combined 5 minutes.

Meeus (2004) witnessed a broken annular from Greece during the 1966 May 20 eclipse and called it a "marvelous sight":

During a few seconds, all around the dark disk of the Moon beads of sunlight were visible through lunar valleys. As the eclipse proceeded, the beads at the eastern limb of the Moon were vanishing while those at the opposite limb increased [...] During a few seconds, the Moon was seen moving with respect to the Sun.

For the 1986 October 3 eclipse, a small group went to extraordinary lengths to witness what turned out to be a broken annular eclipse from a private jet at 40,000 feet (Saemundsson 2021). Accounts of this expedition by Saemundsson and by Schneider (2019) emphasize the difficulty of predicting the precise track and even the type of this marginal eclipse, given the data and technology available at the time. By contrast, several online and mobile applications are now available that predict the appearance and timing of the spectacular Baily's beads display during broken annulars from any location (Jubier 2012; McGlaun 2024; Trainor & Craig 2024).

To our knowledge, no one has attempted to map the locations from which this unique phenomenon can be seen, perhaps because no unambiguous and nonarbitrary criterion for identifying these locations has been proposed. We suggest that they are the locus of points for which

$$\begin{eqnarray*}&&\max ({\boldsymbol{L}})-\min ({\boldsymbol{L}})\gt \delta ,\end{eqnarray*}$$

a criterion inspired by the conditions defining the total and annular limb test zones. It says that broken annularity occurs when the range of limb radii is larger than the distance between the centers of the Sun and Moon. As Figure 13 shows, this produces a path that appears continuous and consistent with the total and annular segments it spans. Note that the test is only valid at locations that have already failed both the totality and annularity tests.

Zoom In Zoom Out Reset image size

6.4. The Solar Radius

The elegant lines and smooth ovals of most eclipse maps reflect the "spherical frictionless cow" assumptions inherent in the Besselian paradigm. Traditional eclipse calculations solve a simplified problem, one involving a smooth, ellipsoidal Earth, a smooth, spherical Moon, and time-varying geometry determined by a handful of low-order polynomials in a noninertial (true of date) coordinate frame. This approach was dictated by the severe limitations of hand calculation, and its errors of approximation were of an order consistent with the accuracy of 19th-century geodesy and ephemerides.

Attempts to bring the accuracy of traditional calculations up to modern standards include the use of multiple lunar radii, multiple Earth ellipsoids, a nudge factor to move the Moon's center of mass closer to its center of figure, and a procedure for sliding limb profiles back and forth in time to find a better intersection between the edges of the Moon and the Sun.

It seems likely that the modern era of eclipse calculation that began in the 1960s was spurred by the confluence of several factors, among them the publication of an English translation of Oppolzer (Oppolzer & Gingerich 1962), the publication of a detailed 20th century treatment of Bessel in the first edition of the Explanatory Supplement (Wilkins & Springett 1961), and the rapidly growing computational capacity available to researchers. This trend has reached an apotheosis of sorts with the proliferation of multi-thousand-year eclipse canons and real-time calculations running in Web browsers.

Rather than merely accelerating a 19th-century methodology, modern computing power can instead be used to move eclipse mapping into the realm of numerical simulation. The shift to computational science has already revolutionized countless fields of science and engineering by embracing computer modeling, statistical analysis, and machine learning (Denning 2017). A common feature of computational science is the ab initio calculation of certain quantities on the points or nodes of a dense mesh, after which the solution of the macroscopic problem emerges as a surface in solution space. A familiar example is the use of hydrocodes for weather prediction—and at least one research and forecast model includes a module for incorporating the atmospheric effects of solar eclipses (Montornès et al. 2016).

For the method described in this paper, the surface is the literal surface of the Earth, and the inputs include dense grids of numbers describing the shapes of the Earth and Moon. This approach has revealed not only the true shape of the umbra shadow but the reason it occurs. It shines light on a type of central eclipse that is neither total nor annular—an impossibility in the Besselian model—and provides an unambiguous criterion for mapping its occurrence. Its improved accuracy can breathe new life into efforts to use eclipses as solar radius yardsticks; in the process, it may surface distinct definitions of astrophysical and eclipse photosphere radii. The current novelty of the method and its results have proven efficacious for engaging the public in discussions of the nature of scientific prediction and computer modeling.

We do not claim that this is the first implementation of numerical eclipse simulation or the only one possible, but we believe it is the first detailed description in the eclipse calculation literature. Our hope is that it spurs further development in this direction.

This work was supported by NASA cooperative agreements NNH15ZDA004C and NNG11HP16A, the Lunar Reconnaissance Orbiter mission, and SAMDA contract NNG17HP01C. We thank the anonymous reviewer for thorough notes and cogent suggestions that helped improve and clarify our manuscript. The GSFC and USNO libraries greatly facilitated the historical research. E.W. thanks his coauthor, his SVS colleague W. Bridgman, and eclipse community members J. Pasachoff, F. Espenak, L. Quaglia, X. Jubier, J. Irwin, D. Herald, S. Habbal, M. Kentrianakis, and C. Fulco for helpful discussion, feedback, and technical guidance. C.A.Y. thanks his coauthor, his former program officer Kristen Erickson, J. Pasachoff, R. Fienberg, A. Speck, and C. Raftery for helpful discussions and support.

The main text of this article is intentionally agnostic about most of the implementation details of an eclipse calculation, but it may be helpful to briefly discuss the utility of high-precision astronomical software libraries, something not covered in traditional treatments. These libraries include SPICE (SPICE 2024), SOFA (SOFA 2024), and NOVAS (USNO 2024). Each of these can be linked directly to programs written in C and Fortran, and, with varying levels of effort and third-party support, to Python and other languages.

The code that generated the NASA maps relies primarily on SPICE. This library gives direct access to the high-precision ephemerides calculated by the Jet Propulsion Laboratory, and it is one of the primary tools used to plan NASA missions and understand the data they return. SPICE is specifically designed for solar system dynamics, while SOFA and NOVAS lean toward more general problems in astrometry.

As an example of the power and simplicity of using one of these libraries, consider the following three lines of C-language code which calculate the topocentric libration and lunar distance of an observer through calls to SPICE library functions. For discussions of the nontrivial mathematics hiding behind these three lines, see Chapront-Touzé & Chapront (1991), Weratschnig et al. (2011), and Liu (2024).

The inputs, functions, and outputs are listed in Table 1. The ICRF position of the Moon, mpos, is found using

The function spkezp_c searches a database of ephemerides, loaded previously, for objects with the identifiers ID_MOON and ID_EARTH and ephemeris coverage encompassing time t . Assuming the search succeeds, it calculates the position of the Moon with respect to the Earth, in the ICRF coordinate frame, accounting for both light travel time and stellar aberration ("lt+s"). The ephemeris data are stored in files called SPICE kernels, which are typically loaded at the start of the program. In this case, the kernel would be de440.bsp (Park et al. 2021). ⁵

The rotation matrix mmat, which transforms position vectors from ICRF to the DE421 Moon ME body-fixed frame, is found with

This call finds the orientation ephemeris and frame definition it needs in the kernels moon_pa_de440_200625.bpc and moon_de440_200625.tf. The complicated identifier for the Moon ME frame reflects the fact 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.

Finding the ICRF position of the observer, obs, is a two-step process that starts by converting geographic coordinates and height above mean sea level into an Earth body-fixed position vector:

As noted in the text, the map pixel coordinates i, j map to longitude and latitude according to the projection chosen for the map, and these are scaled to radians to form the first two arguments. The third argument, h , is the ellipsoid height, the elevation above the surface of the ellipsoid defined by the Earth equatorial radius earth_r and the polar flattening. This differs from the orthometric (mean sea level) elevations in Earth DEMs.

Mean sea level (MSL) is an undulating equipotential surface in the gravitational field of the Earth, the level at which the ocean would rest in the absence of wind, currents, and tides. MSL is irregular because the Earth is both lumpy and bumpy—its mass is distributed inhomogeneously. The difference between MSL and the surface of the ellipsoid is the geoid height N, and this must be added to the orthometric height H obtained from the DEM to find h , the ellipsoid height.

In practice, this means sampling a geoid height raster in addition to the DEM. The DEM and the call to georec_c must agree in both the geoid model (implicit in h ) and the ellipsoid. Most modern global Earth DEMs (e.g., SRTM, ASTER, ALOS, GMTED2010, MERIT, Earth2014) use the WGS84 ellipsoid and the EGM96 geoid model, but many local and regional DEMs do not.

The second step involves finding and applying the rotation matrix, here called omat, that transforms Earth body-fixed coordinates at time t to ICRF. SPICE returns an extremely precise matrix for times in the recent past, using the following call and the Earth orientation kernel earth_latest_high_prec.bpc:

But as of this writing, SPICE Earth orientation kernels are less reliable or nonexistent for times in the future or before 1962. An equivalent pair of calls in SOFA is

The first line actually returns the inverse transformation, but since this is a pure rotation, the inverse is just the transpose (the second line). With either library, omat can be calculated once per time step. The transformation in SPICE of the observer obs is then

The 10 lines of code exhibited here are sufficient to place the Moon and the observer in a common coordinate frame, find the distance between them, and calculate the Moon's topocentric libration at the observer. Astronomical libraries make this kind of parsimony possible with most of the calculations involved in eclipse prediction.