Solar Eclipse Predictions¶
On 8 April 2024 the Moon's shadow crossed North America from Mexico's Pacific coast to Newfoundland. This tutorial computes where that eclipse was total: the central line and the path of totality between its northern and southern limits, and the duration of totality. It then checks the results against NASA's published predictions.
We use the classical method of Besselian elements, the method used by NASA and the national almanacs (Urban & Seidelmann 2013; Meeus 1989). The tutorial also walks through the algorithm. Each quantity is derived before it is used, and each piece of code sits next to the formulas it implements. You need vectors and a little trigonometry, and nothing else about eclipses.
- The idea, in plain language.
- Inputs and conventions: the Sun and Moon, the Earth's rotation, time scales, and constants.
- The fundamental plane, where the shadow becomes simple.
- The Besselian elements, derived one at a time.
- Fitting the elements with polynomials, and comparing them with NASA's.
- From the plane to the ground: observer coordinates, the central line, the limits and the duration of totality.
- The 2024 path: a map, local circumstances for eleven cities, a comparison with NASA's path table, and an independent cross-check.
- The same code applied to 12 August 2045.
- Limitations and references.
The idea in plain language¶
The Moon, lit by the Sun, casts two cones of shadow.
The umbra is the full shadow, where the Moon hides the whole Sun. It is shaped like the sharpened tip of a pencil. At the Moon it is as wide as the Moon, 3,474 km, and it narrows steadily behind the Moon to a point about 374,000 km away. The Moon is between 356,000 and 407,000 km from the Earth, so the tip sometimes reaches past the Earth's surface and sometimes falls short. Totality happens wherever the umbral cone meets the ground. When the tip falls short, observers beyond it, in the antumbra, see a ring of Sun around the Moon: an annular eclipse.
The penumbra widens away from the Moon. Inside it the Sun is only partly hidden: a partial eclipse.
As the Moon moves along its orbit and the Earth turns, the umbra's footprint slides across the ground. The strip it sweeps out is the path of totality. The two edges of that strip are the northern and southern limits. The track of the cone's axis is the central line, where totality lasts longest. How long a place stays inside the footprint is its duration of totality.
Why a plane? Seen from the rotating Earth, the geometry is awkward. The cone is tilted, and the ground is a curved, flattened, turning surface. The footprint is an ellipse that stretches as the Sun gets lower. The method of Besselian elements splits the problem in two:
- The shadow. Take the plane through the Earth's centre that is perpendicular to the shadow axis: the fundamental plane. Each cone cuts it in an exact circle. The circle's centre moves along an almost straight line at an almost constant speed, and its radius barely changes. A handful of smooth functions of time, the Besselian elements, describe it completely. They depend only on the Sun and the Moon, and they are the same for every observer.
- The Earth. Where a given place is, relative to that plane, at a given time. This is pure geometry: rotate the place into the plane's axes and compare its distance from the shadow axis with the circle's radius.
The umbral cone is very narrow: its half-angle is about 0.27°. Its radius therefore changes by only about 30 km between the fundamental plane and the ground, a correction handled with one term.
The sketch below shows how the Sun and the Moon define the two cones. The edges of the umbra (solid lines) are the external tangents to the two limbs: they meet at the umbral vertex, and past it they open out again into the antumbra. The edges of the penumbra (dashed) are the internal tangents, which cross the axis between the Sun and the Moon. $f_1$ and $f_2$ are the cones' half-angles, the angles between their edges and the shadow axis. Where the cones reach the ground, the umbra makes a small spot of total eclipse inside the penumbra's larger patch of partial eclipse. The enlarged panel follows the cones through the Earth to the fundamental plane, where their radii are $l_1$ and $l_2$. It shows a total eclipse: the umbral vertex lies beyond the plane, so the plane cuts the converging umbra, and in the sign convention derived below $l_2 < 0$, with radius $|l_2|$. The sketch is not to scale: the Moon is really about 56 Earth radii from the Earth, the Sun 400 times farther still, and the cones are far thinner.
import datetime
import numpy as np
import satkit as sk
import pandas as pd
import ssl, certifi
ssl._create_default_https_context = lambda: ssl.create_default_context(cafile=certifi.where())
import matplotlib.pyplot as plt
from matplotlib.patches import Arc
import scienceplots # noqa: F401
plt.style.use(["science", "no-latex", "../satkit.mplstyle"])
%config InlineBackend.figure_formats = ['svg']
import warnings
warnings.filterwarnings("ignore", "Downloading")
import cartopy.crs as ccrs
import cartopy.feature as cfeature
# Schematic of the shadow cones and the fundamental plane (not to scale)
from matplotlib.patches import Circle, Polygon, Rectangle, Wedge
from matplotlib.lines import Line2D
from matplotlib.patches import Ellipse
from matplotlib.colors import to_rgb
TEXTURED = False # True: shaded, lightly textured bodies instead of flat fills
# Lengths in Earth radii. The radii of the Earth, the Moon and the cones on the plane are in rough
# proportion (l1 about 0.7, |l2| about 0.1); the distances are far too short, so the Sun is too
# small and the angles f1 and f2 are exaggerated.
sun_r, moon_r, earth_r = 1.4, 0.33, 1.0
moon_x, a = 8.8, 1.8 # Sun at x = 0; a is the Moon's height z above the plane
earth_x = moon_x + a
INK, PLANE, SHADE, ANGLE = "#333333", "#009988", "#4a5878", "#CC3311"
SUN, MOON, EARTH, EARTH_EDGE, LIGHT = "#f3d9a4", "#c9c9c9", "#dbe7f1", "#0077BB", "#fbf3dd"
# Umbra: the external tangents to the Sun and the Moon, meeting at the umbral vertex
sin_f2 = (sun_r - moon_r) / moon_x
f2 = np.arcsin(sin_f2)
vertex2 = moon_x + moon_r / sin_f2
# Penumbra: the internal tangents, crossing the axis between the Sun and the Moon
sin_f1 = (sun_r + moon_r) / moon_x
f1 = np.arcsin(sin_f1)
vertex1 = moon_x - moon_r / sin_f1
# The umbral vertex lies beyond the plane, so the plane cuts the converging umbra: total, l2 < 0
l1 = (earth_x - vertex1) * np.tan(f1)
l2 = (earth_x - vertex2) * np.tan(f2)
x_end = earth_x + 2.15
def umbra_y(x): # upper edge of the umbra (and, past the vertex, of the antumbra)
return (vertex2 - x) * np.tan(f2)
def penumbra_y(x): # upper edge of the penumbra
return (x - vertex1) * np.tan(f1)
def hits_earth(edge):
# x where an upper cone edge meets the Earth's near side
xs = np.linspace(earth_x - earth_r, earth_x, 20001)
return xs[np.argmin(np.abs((xs - earth_x) ** 2 + edge(xs) ** 2 - earth_r**2))]
xu, xp = hits_earth(umbra_y), hits_earth(penumbra_y)
ang_u = np.degrees(np.arctan2(umbra_y(xu), earth_x - xu)) # half-angle of the spots, seen
ang_p = np.degrees(np.arctan2(penumbra_y(xp), earth_x - xp)) # from the Earth's centre
def disc(ax, x, r, kind, zorder, edge=None):
# A body: a flat fill, or (TEXTURED) a shaded disc lit from the Sun's side, on the left
base = {"sun": SUN, "moon": MOON, "earth": EARTH}[kind]
if not TEXTURED:
ax.add_patch(Circle((x, 0), r, facecolor=base, edgecolor=edge or "none", lw=0.8, zorder=zorder))
return
X, Y = np.meshgrid(np.linspace(-1, 1, 301), np.linspace(-1, 1, 301))
mu = np.sqrt(np.clip(1 - X**2 - Y**2, 0, 1))
if kind == "sun": # warm core, orange limb, limb darkening
core, limb = np.array(to_rgb("#f8dc8c")), np.array(to_rgb("#eaa45a"))
img = core + (limb - core) * (1 - mu[..., None]) ** 0.6
for k, al in ((1.25, 0.10), (1.12, 0.14)): # faint glow
ax.add_patch(Circle((x, 0), k * r, facecolor="#f6c56a", alpha=al, lw=0, zorder=zorder - 0.1))
else:
shade = 0.82 + 0.18 * (0.5 * mu - 0.5 * X) # brighter towards the Sun
img = np.array(to_rgb(base)) * shade[..., None]
im = ax.imshow(np.clip(img, 0, 1), extent=(x - r, x + r, -r, r), origin="lower", zorder=zorder,
interpolation="bilinear")
clip = Circle((x, 0), r, transform=ax.transData)
im.set_clip_path(clip)
spots = {"moon": [(-0.3, 0.35, 0.18), (0.25, -0.2, 0.24), (0.35, 0.45, 0.12), (-0.4, -0.45, 0.14)],
"earth": [(-0.35, 0.3, 0.5), (0.3, -0.35, 0.6), (0.45, 0.55, 0.35)]}.get(kind, [])
for cx, cy, cr in spots: # low-contrast craters, or continent-like blobs
colour = "#7d7d7d" if kind == "moon" else "#8fb58a"
p = Ellipse((x + cx * r, cy * r), 2 * cr * r, 1.4 * cr * r, angle=30 * cx, facecolor=colour,
alpha=0.18 if kind == "moon" else 0.3, lw=0, zorder=zorder + 0.01)
ax.add_patch(p)
p.set_clip_path(clip)
if edge:
ax.add_patch(Circle((x, 0), r, facecolor="none", edgecolor=edge, lw=0.8, zorder=zorder + 0.02))
def draw_shadow(ax, inside):
# Cones from the Moon to the Earth's surface; with inside=True they continue through the Earth
# to the fundamental plane (penumbra) and past the umbral vertex (antumbra)
xm = moon_x - moon_r * sin_f1
x_pen = x_end if inside else earth_x
ax.add_patch(Polygon([(xm, penumbra_y(xm)), (x_pen, penumbra_y(x_pen)), (x_pen, -penumbra_y(x_pen)),
(xm, -penumbra_y(xm))], facecolor=SHADE, alpha=0.2, lw=0, zorder=1))
xm = moon_x + moon_r * sin_f2
ax.add_patch(Polygon([(xm, umbra_y(xm)), (vertex2, 0), (xm, -umbra_y(xm))],
facecolor=SHADE, alpha=0.9, lw=0, zorder=1))
ax.add_patch(Polygon([(vertex2, 0), (x_end, -umbra_y(x_end)), (x_end, umbra_y(x_end))],
facecolor=SHADE, alpha=0.12, lw=0, zorder=1))
for s in (1, -1):
ax.plot([vertex2, x_end], [0, -s * umbra_y(x_end)], color=INK, lw=0.6, alpha=0.6, zorder=4)
fig = plt.figure(figsize=(10, 5.8))
ax = fig.add_axes([0.0, 0.45, 1.0, 0.5])
axd = fig.add_axes([0.5, 0.02, 0.46, 0.36])
# ---- Overview: Sun, Moon and Earth ----
# Sunlight grazing the Moon, between the Sun's limb and the Moon's limb
ax.add_patch(Polygon([(sun_r * sin_f2, sun_r * np.cos(f2)), (moon_x + moon_r * sin_f2, moon_r * np.cos(f2)),
(moon_x + moon_r * sin_f2, -moon_r * np.cos(f2)), (sun_r * sin_f2, -sun_r * np.cos(f2))],
facecolor=LIGHT, lw=0, zorder=0))
draw_shadow(ax, inside=False)
disc(ax, 0, sun_r, "sun", 3)
disc(ax, moon_x, moon_r, "moon", 3)
disc(ax, earth_x, earth_r, "earth", 2, EARTH_EDGE)
# Where the shadow meets the ground: the umbra's small spot inside the penumbra's larger patch
ax.add_patch(Arc((earth_x, 0), 2 * earth_r, 2 * earth_r, theta1=180 - ang_p, theta2=180 + ang_p,
color="#9aa6bd", lw=5, zorder=3))
ax.add_patch(Arc((earth_x, 0), 2 * earth_r, 2 * earth_r, theta1=180 - ang_u, theta2=180 + ang_u,
color=SHADE, lw=5.5, zorder=3))
# Tangent lines from the Sun's limb to the Moon's limb, as far as the Earth
for s in (1, -1):
ax.plot([sun_r * sin_f2, xu], [s * sun_r * np.cos(f2), s * umbra_y(xu)], color=INK, lw=0.9, zorder=4)
ax.plot([sun_r * sin_f1, xp], [-s * sun_r * np.cos(f1), s * penumbra_y(xp)], color=INK, lw=0.9, ls="--",
zorder=4)
ax.plot([-sun_r - 0.3, x_end + 0.1], [0, 0], color="0.45", lw=0.6, ls=":", zorder=5) # shadow axis
ax.plot([earth_x, earth_x], [-1.25, 1.25], color=PLANE, lw=1.4, zorder=5)
ax.text(earth_x, 1.3, "fundamental plane", color=PLANE, fontsize=9, ha="center", va="bottom")
# f1, where the penumbral edges cross the axis
r = 2.6
ax.add_patch(Arc((vertex1, 0), 2 * r, 2 * r, theta1=180 - np.degrees(f1), theta2=180, color=ANGLE, lw=1.4,
zorder=5))
ax.text(vertex1 - r - 0.08, 0.5 * r * np.sin(f1), r"$f_1$", fontsize=11, color=ANGLE, ha="right", va="center")
# Labels in columns, with thin leaders ending on the feature
def leader(text, xy, xytext, ha="right", target=ax):
target.annotate(text, xy=xy, xytext=xytext, fontsize=9, ha=ha, va="center", zorder=6,
arrowprops=dict(arrowstyle="-", lw=0.5, color="0.35", shrinkA=2, shrinkB=0))
target.plot(*xy, "o", color="0.2", ms=2, zorder=6)
xg = moon_x + moon_r + 0.12
leader("penumbra — partial", (xg + 0.05, 0.5 * (umbra_y(xg + 0.05) + penumbra_y(xg + 0.05))), (7.6, 1.75))
leader("umbra — total", (xg, 0.0), (7.6, 1.3))
leader("total eclipse", (earth_x - earth_r * np.cos(np.radians(0.5 * ang_u)), -earth_r * np.sin(np.radians(0.5 * ang_u))),
(8.3, -1.45))
leader("partial eclipse", (earth_x - earth_r * np.cos(np.radians(0.5 * (ang_u + ang_p))),
-earth_r * np.sin(np.radians(0.5 * (ang_u + ang_p)))), (8.3, -1.85))
leader("antumbra — annular", (earth_x + 1.55, 0.0), (x_end - 0.1, 1.3), ha="center")
ax.text(0, -sun_r - 0.1, "Sun", fontsize=10, ha="center", va="top")
ax.text(moon_x, -0.5, "Moon", fontsize=10, ha="center", va="top")
ax.text(earth_x + 0.75, -1.0, "Earth", fontsize=10, ha="left", va="top")
ax.set_xlim(-sun_r - 0.3, x_end + 0.4)
ax.set_ylim(-2.05, 2.0)
ax.set_aspect("equal")
ax.axis("off")
ax.set_title("Shadow cones and the fundamental plane (schematic, not to scale)", fontsize=10)
# Key for the line styles, in the space under the Sun
fig.legend(handles=[Line2D([], [], color=INK, lw=0.9, label="external tangents: edges of the umbra"),
Line2D([], [], color=INK, lw=0.9, ls="--", label="internal tangents: edges of the penumbra")],
loc="upper left", bbox_to_anchor=(0.03, 0.36), fontsize=9, frameon=False, handlelength=2.2)
# ---- Detail at the Earth: the cones continued to the fundamental plane ----
box = (earth_x - 0.95, earth_x + 2.05, -0.9, 0.9)
disc(axd, earth_x, earth_r, "earth", 0.5, EARTH_EDGE)
draw_shadow(axd, inside=True)
for s in (1, -1):
axd.plot([box[0], vertex2], [s * umbra_y(box[0]), 0], color=INK, lw=0.9, zorder=4)
axd.plot([box[0], box[1]], [s * penumbra_y(box[0]), s * penumbra_y(box[1])], color=INK, lw=0.9, ls="--", zorder=4)
axd.plot([box[0], box[1]], [0, 0], color="0.45", lw=0.6, ls=":", zorder=4)
axd.plot([earth_x, earth_x], [box[2] + 0.1, box[3] - 0.1], color=PLANE, lw=1.4, zorder=3)
def dimension(x, y, colour):
# dimension line at x from the axis to height y, with short end ticks and a leader to the plane
axd.annotate("", xy=(x, y), xytext=(x, 0), zorder=5,
arrowprops=dict(arrowstyle="<|-|>", lw=0.8, color=colour, shrinkA=0, shrinkB=0,
mutation_scale=7))
for yy in (0, y):
axd.plot([earth_x, x + 0.03], [yy, yy], color=colour, lw=0.6, zorder=5)
dimension(earth_x + 0.12, l1, INK)
axd.text(earth_x + 0.17, 0.5 * l1, r"$l_1$", fontsize=11, ha="left", va="center")
dimension(earth_x + 0.12, l2, "white")
axd.text(earth_x + 0.17, l2 - 0.04, r"$|l_2|$", fontsize=11, ha="left", va="top")
# zeta-hat from the Earth's centre towards the Moon and the Sun
axd.annotate("", xy=(earth_x - 0.7, 0), xytext=(earth_x, 0), zorder=6,
arrowprops=dict(arrowstyle="-|>", lw=1.2, color="white", mutation_scale=10, shrinkA=0, shrinkB=0))
axd.text(box[0] + 0.1, 0.37, r"$\hat{\zeta}$ (towards" + "\nMoon & Sun)", fontsize=9, ha="left", va="center",
linespacing=1.3)
# The umbral vertex and f2, on the antumbra side where there is room
axd.plot(vertex2, 0, "o", color=INK, ms=3.5, zorder=6)
axd.text(vertex2 - 0.02, 0.1, "umbral vertex", fontsize=9, ha="right", va="bottom")
r = 0.95
axd.add_patch(Wedge((vertex2, 0), r, -np.degrees(f2), 0, facecolor=ANGLE, alpha=0.25, lw=0, zorder=5))
axd.add_patch(Arc((vertex2, 0), 2 * r, 2 * r, theta1=-np.degrees(f2), theta2=0, color=ANGLE, lw=1.6, zorder=6))
axd.text(vertex2 + r + 0.03, -0.5 * r * np.sin(f2), r"$f_2$", fontsize=11, color=ANGLE, ha="left", va="center")
axd.text(box[1] - 0.04, umbra_y(box[1]) - 0.02, "antumbra", fontsize=9, ha="right", va="top")
axd.set_xlim(box[0], box[1])
axd.set_ylim(box[2], box[3])
axd.set_aspect("equal")
axd.set_xticks([])
axd.set_yticks([])
for sp in axd.spines.values():
sp.set_color("0.6")
sp.set_linewidth(0.6)
axd.set_title("Boxed region, enlarged: the cones continued through the Earth to the fundamental plane\n"
r"(through the Earth's centre, $\perp$ shadow axis)", fontsize=9, linespacing=1.4)
# Mark the detail's region on the overview and join the two
ax.add_patch(Rectangle((box[0], box[2]), box[1] - box[0], box[3] - box[2], fill=False, ec="0.6", lw=0.6,
zorder=6))
plt.show()
Inputs and conventions¶
The calculation needs four ingredients: where the Sun and the Moon are, how the Earth is oriented, which time scale each of those uses, and a few physical constants.
The Sun and the Moon¶
satkit.jplephem.geocentric_pos gives the geocentric positions of the Sun and the Moon in the GCRF, the inertial frame centred on the Earth. The positions come from the JPL DE440 ephemeris (Park et al. 2021). These are geometric positions: where the bodies are at the instant $t$. The shadow falling on the Earth at $t$, however, is made of light that left the Sun about 499 s earlier and grazed the Moon about 1.2 s earlier. We therefore need the positions at those two earlier instants.
We work in the geocentric frame, which moves with the Earth. Over a few minutes this frame is inertial to far better than we need, so light travels through it in straight lines at speed $c$. The ray that reaches the Earth at time $t$ left the geocentric Sun at $t - \tau_\odot$ and passed the geocentric Moon at $t - \tau_{\rm M}$, where
$$ \tau_\odot = \frac{|\mathbf r_\odot|}{c} \approx 499\ {\rm s}, \qquad \tau_{\rm M} = \frac{|\mathbf r_{\rm M}|}{c} \approx 1.2\ {\rm s}. $$
One evaluation of each is enough: in 499 s the Sun's distance changes by at most a few hundred kilometres, which changes $\tau_\odot$ by less than a millisecond.
- Sun: light time and aberration. The Sun itself hardly moves in 499 s. The geocentric Sun, however, moves because the Earth does: $\mathbf v_\oplus \tau_\odot \approx 30\ {\rm km\,s^{-1}} \times 499\ {\rm s} \approx 15{,}000$ km, or 20.5″ as seen from the Earth. This is the annual aberration, and $\mathbf r_\odot(t-\tau_\odot)$ is the apparent Sun of the almanacs. The shadow axis passes through the Moon, so turning the Sun's direction by 20.5″ swings the axis about the Moon by the same angle. At the Earth, about 360,000 km beyond the Moon, that moves the shadow by about 36 km.
- Moon: light time. In 1.2 s the Moon moves about 1.2 km relative to the Earth, which moves the shadow by about 1.3 km.
Both effects are computed below, once the elements are defined.
From the inertial frame to the rotating Earth¶
satkit.frametransform.qgcrf2itrf rotates both vectors into the Earth-fixed ITRF. The rotation has three parts: precession and nutation (IAU 2006/2000A), the Earth rotation angle, and polar motion (Petit & Luzum 2010). The Earth rotation angle is a linear function of UT1, the time scale that measures how far the Earth has actually turned. The rotation uses UT1 and polar motion from the IERS Earth orientation parameters.
In the ITRF the Earth stands still. The direction of the shadow axis then gives two of the elements directly, like the latitude and longitude of a point on the sky:
- $d$ is the declination of the shadow axis: its angle above the ITRF equator.
- $\mu$ is its Greenwich hour angle: the angle from the Greenwich meridian to the meridian of the axis, measured westward. It is minus the ITRF longitude of the axis direction. $\mu$ grows by about $15.004^\circ$ per hour. That is the Earth's rotation rate relative to the stars, $15.041^\circ$ per hour, less the Sun's eastward drift along the ecliptic, about $0.04^\circ$ per hour.
The classical calculation gets $\mu$ as the Greenwich sidereal time minus the axis's right ascension. Rotating into the ITRF gives the same angle without sidereal time.
Time scales and $\Delta T$¶
Three time scales appear:
- TT (Terrestrial Time) is uniform atomic time, TAI + 32.184 s. The JPL ephemeris uses TDB, which differs from TT by less than 2 ms. The positions of the Sun and Moon, and so the elements $x$, $y$, $d$, $l_1$ and $l_2$, are functions of TT.
- UT1 is the Earth's rotation angle expressed as a time. It is measured (by VLBI), not predicted. $\mu$ depends on UT1.
- UTC ticks in atomic seconds and is kept within 0.9 s of UT1 by leap seconds. Clock times in this tutorial are UTC. Since 2017, TT − UTC = 69.184 s.
The difference $\Delta T = {\rm TT} - {\rm UT1} = 69.184\ {\rm s} - ({\rm UT1} - {\rm UTC})$ is what an eclipse prediction has to guess. Tidal friction slows the Earth's rotation, and the core, oceans and atmosphere add irregular changes. As a result $\Delta T$ drifts by fractions of a second per year and cannot be predicted to better than seconds, years ahead. Predictions made in advance therefore compute the elements in TT, which the ephemeris fixes, and assume a $\Delta T$ to place them on the rotating Earth. NASA assumed $\Delta T$ = 70.6 s for 2024. One second of $\Delta T$ is $15.04''$ of Earth rotation, which shifts the whole path $15''$ in longitude: about 0.46 km × cos(latitude).
We compute the eclipse after the fact, with the measured UT1, so our $\mu$ corresponds to the actual $\Delta T$ of 69.2 s. When we compare our elements with NASA's, we first convert NASA's $\mu$ to the same $\Delta T$.
In the code, the reference epoch $t_0$ is 18:00:00 TT (17:58:50.816 UTC), and times are handled as hours from $t_0$ (th). These are elapsed SI hours, which are the same in TT and UTC when no leap second falls in between.
Constants¶
| Quantity | Code | Value | Source |
|---|---|---|---|
| Earth's equatorial radius, the unit of length ($a_e$) | AE |
6378.137 km | WGS84 (NGA 2014) |
| Earth's flattening $f$ | FLAT |
1/298.257223563 | WGS84 |
| Lunar radius for the penumbra, $k_1$ | K1 |
0.272488 $a_e$ = 1737.97 km | NASA (Espenak) |
| Lunar radius for the umbra, $k_2$ | K2 |
0.272281 $a_e$ = 1736.65 km | NASA |
| Solar radius $R_\odot$ | R_SUN |
959.63″ at 1 au = 696,000 km = 109.12 $a_e$ | NASA |
Why two lunar radii? The Moon's limb is not smooth. Mountains and valleys make its profile ragged by a few kilometres. A partial eclipse begins when the first mountain touches the Sun's disk, so the penumbra uses a larger radius. Totality needs every valley closed. Sunlight shining through the lowest valleys (Baily's beads) ends totality early, so the umbra uses a radius close to the valley floors. The value $k = 0.2725076$ (1738.1 km), adopted by the IAU in 1982 as the Moon's mean radius, would make the path about 3 km wider than NASA's. The current IAU mean radius from lunar altimetry is 1737.4 km. We use NASA's two values so that our elements can be compared with NASA's.
The Sun's radius. 959.63″ is Auwers's 19th-century measurement of the solar semi-diameter at 1 au. It has been the almanac standard ever since, and NASA uses it. satkit.consts.sun_radius is the IAU 2015 nominal photospheric radius, 695,700 km. It would make the umbra about 1.5 km wider.
AE = sk.consts.wgs84_a # unit of length: Earth equatorial radius, m
FLAT = sk.consts.wgs84_f
K1, K2 = 0.272488, 0.272281 # lunar radius for the penumbra and umbra, Earth radii
ARCSEC = np.pi / 180 / 3600
R_SUN = sk.consts.au * np.tan(959.63 * ARCSEC) / AE # solar radius, Earth radii
def sun_moon_itrf(times):
# Apparent Sun and light-time-corrected Moon in the ITRF, in Earth radii
c = sk.consts.c
sun0 = sk.jplephem.geocentric_pos(sk.solarsystem.Sun, times)
moon0 = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, times)
t_sun = [t - sk.duration.from_seconds(tau) for t, tau in zip(times, np.linalg.norm(sun0, axis=1) / c)]
t_moon = [t - sk.duration.from_seconds(tau) for t, tau in zip(times, np.linalg.norm(moon0, axis=1) / c)]
sun = sk.jplephem.geocentric_pos(sk.solarsystem.Sun, t_sun)
moon = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, t_moon)
q = sk.frametransform.qgcrf2itrf(times)
sun = np.array([qi * v for qi, v in zip(q, sun)]) / AE
moon = np.array([qi * v for qi, v in zip(q, moon)]) / AE
return sun, moon
The fundamental plane¶
Let $\mathbf r_\odot$ and $\mathbf r_{\rm M}$ be the Sun and Moon from sun_moon_itrf, in the ITRF and in Earth radii. The shadow axis is the line through the centres of the Sun and the Moon. Its direction, pointing from the Moon towards the Sun, is
$$ \mathbf G = \mathbf r_\odot - \mathbf r_{\rm M}, \qquad \mathbf g = \frac{\mathbf G}{|\mathbf G|}. $$
The fundamental plane passes through the Earth's centre, perpendicular to $\mathbf g$. We give it three unit vectors:
- $\hat{\boldsymbol\zeta} = \mathbf g$, along the axis, towards the Moon and the Sun. Writing $\mathbf g$ in terms of its declination $d$ and Greenwich hour angle $\mu$ (hour angles are measured westward, hence the minus sign in the $y$ component) defines these two elements:
$$ \hat{\boldsymbol\zeta} = (\cos d\cos\mu,\ -\cos d\sin\mu,\ \sin d) \quad\Longrightarrow\quad d = \arcsin g_z, \qquad \mu = \operatorname{atan2}(-g_y,\ g_x). $$
- $\hat{\boldsymbol\xi}$, east, in the plane. It is perpendicular both to the axis and to the Earth's rotation axis $\hat{\mathbf k} = (0, 0, 1)$, so it lies in the equatorial plane:
$$ \hat{\boldsymbol\xi} = \frac{\hat{\mathbf k}\times\hat{\boldsymbol\zeta}}{|\hat{\mathbf k}\times\hat{\boldsymbol\zeta}|} = (\sin\mu,\ \cos\mu,\ 0). $$
- $\hat{\boldsymbol\eta}$, north, in the plane, completing the right-handed set:
$$ \hat{\boldsymbol\eta} = \hat{\boldsymbol\zeta}\times\hat{\boldsymbol\xi} = (-\sin d\cos\mu,\ \sin d\sin\mu,\ \cos d). $$
Since $\hat{\mathbf k} = \cos d\,\hat{\boldsymbol\eta} + \sin d\,\hat{\boldsymbol\zeta}$, $\hat{\boldsymbol\eta}$ is the direction of the Earth's north pole projected onto the plane.
Any vector $\mathbf r$ (ITRF, in Earth radii) then has fundamental-plane coordinates $(\mathbf r\cdot\hat{\boldsymbol\xi},\ \mathbf r\cdot\hat{\boldsymbol\eta},\ \mathbf r\cdot\hat{\boldsymbol\zeta})$. Seen from the Sun, the fundamental plane is a map of the Earth's disk: $\xi$ to the right (east), $\eta$ up (north), and $\zeta$ towards the viewer, measuring height above the plane.
The plane goes through the Earth's centre, not the ground, so that it is the same for every observer. An observer's height above it, $\zeta$, lies between 0 and 1 on the day side, and it enters through the narrow opening angle of the cones.
In the code the three unit vectors are xhat, yhat and zhat (in elements_from and plane_axes). The sketch below shows the plane in three dimensions. It is a schematic, not to scale: the tilt $d$ and the umbra are exaggerated.
# The fundamental plane in 3-D (schematic, not to scale). Vectors are (xi, eta, zeta) components.
e1 = np.array([-0.8, 0.0, 1.0])
e1 /= np.linalg.norm(e1)
e2 = np.array([-0.35, 1.0, -0.28])
e2 /= np.linalg.norm(e2)
tilt = np.arctan2(e2[2], e1[2]) # turn the page so that zeta points to the right
page = np.array([[np.cos(tilt), np.sin(tilt)], [-np.sin(tilt), np.cos(tilt)]])
def proj(v):
# orthographic projection onto the page
v = np.asarray(v, dtype=float)
return np.stack([v @ e1, v @ e2], axis=-1) @ page.T
q = np.linspace(0, 2 * np.pi, 200)
def ring(r, centre=(0, 0, 0)):
# circle of radius r parallel to the fundamental plane
return proj(np.stack([centre[0] + r * np.cos(q), centre[1] + r * np.sin(q), np.full_like(q, centre[2])], axis=-1))
fig, ax = plt.subplots(figsize=(9, 5.2))
ax.add_patch(plt.Circle((0, 0), 1.0, color="#0077BB", alpha=0.15, lw=0))
ax.fill(*ring(1.75).T, color="#009988", alpha=0.12, lw=0)
ax.plot(*ring(1.75).T, color="#009988", lw=1.0)
ax.text(*proj([0.9, -1.75, 0]), "fundamental plane", color="#009988", fontsize=10, va="top")
ax.text(0.45, 0.55, "Earth", color="#0077BB", fontsize=10)
# Shadow axis, the Moon and the umbral cone
A = np.array([0.55, 0.45, 0.0]) # where the axis crosses the plane: (x, y)
zeta_hat = np.array([0.0, 0.0, 1.0])
moon_c, moon_r, l2_draw = A + 3.4 * zeta_hat, 0.35, 0.12
ax.plot(*np.array([proj(A - 1.3 * zeta_hat), proj(A + 4.6 * zeta_hat)]).T, color="0.3", lw=0.8, ls=":")
ax.text(*(proj(A + 2.0 * zeta_hat) + [0, -0.45]), "shadow axis (dotted),\ntowards the Moon and the Sun",
fontsize=9, ha="center", va="top", color="0.3")
ax.fill(*ring(moon_r, moon_c).T, color="#888888", alpha=0.6, lw=0)
ax.text(*(proj(moon_c) + [0, -0.5]), "Moon", fontsize=10, ha="center", va="top")
for s in (1, -1):
ax.plot(*np.array([proj(moon_c + [0, s * moon_r, 0]), proj(A + [0, s * l2_draw, 0])]).T, color="k", lw=0.9)
ax.fill(*ring(l2_draw, A).T, color="k", alpha=0.7, lw=0)
ax.plot(*proj(A), "o", color="w", ms=2)
ax.annotate("umbra on the plane: a circle\nof radius $|l_2|$ about $(x, y)$", xy=proj(A + [0, -l2_draw, 0]),
xytext=(0.35, -1.35), fontsize=9, arrowprops=dict(arrowstyle="->", lw=0.6))
# Unit vectors and the Earth's north pole
for v, label, off in [((1.35, 0, 0), r"$\xi$ (east)", (-0.1, -0.25)), ((0, 1.35, 0), r"$\eta$ (north)", (0.06, 0.02)),
((0, 0, 1.6), r"$\zeta$", (0.02, 0.08))]:
ax.annotate("", xy=proj(v), xytext=(0, 0), arrowprops=dict(arrowstyle="-|>", lw=1.3, color="#CC3311"))
ax.text(*(proj(v) + off), label, fontsize=11, color="#CC3311")
d_draw = np.radians(25)
pole = 1.35 * np.array([0, np.cos(d_draw), np.sin(d_draw)])
ax.plot(*np.array([[0, 0], proj(pole)]).T, color="#0077BB", lw=1.0, ls="--")
ax.text(*(proj(pole) + [0.05, 0.02]), "N pole", fontsize=9, color="#0077BB")
ax.plot(0, 0, "k.", ms=4)
ax.set_xlim(-2.0, 4.6)
ax.set_ylim(-1.9, 1.9)
ax.set_aspect("equal")
ax.axis("off")
plt.show()
The Besselian elements, one at a time¶
| Element | Meaning | Code |
|---|---|---|
| $x$, $y$ | where the shadow axis crosses the fundamental plane | x, y |
| $d$ | declination of the shadow axis | d (degrees) |
| $\mu$ | Greenwich hour angle of the shadow axis | mu (degrees) |
| $f_1$, $f_2$ | half-angles of the penumbral and umbral cones | tan_f1, tan_f2 |
| $l_1$, $l_2$ | radii of the penumbral and umbral cones on the fundamental plane | l1, l2 |
Lengths are in Earth radii. $d$ and $\mu$ were defined above; the rest follow.
$x$ and $y$: where the axis crosses the plane¶
The shadow axis is parallel to $\hat{\boldsymbol\zeta}$ and passes through the Moon's centre. Every point on it therefore has the same $\xi$ and $\eta$ as the Moon, and it crosses the plane at the Moon's projection:
$$ x = \hat{\boldsymbol\xi}\cdot\mathbf r_{\rm M}, \qquad y = \hat{\boldsymbol\eta}\cdot\mathbf r_{\rm M}, \qquad z = \hat{\boldsymbol\zeta}\cdot\mathbf r_{\rm M}. $$
$z$, the Moon's height above the plane (about 56 Earth radii), is not a published element, but the cone radii need it. $x$ increases by about 0.51 Earth radii per hour: that is the Moon's orbital motion relative to the direction of the Sun, about 3,300 km/h. When $\sqrt{x^2 + y^2} < 1$, roughly, the axis hits the Earth and the eclipse is central. The minimum of $\sqrt{x^2+y^2}$ is the eclipse's $\gamma$.
$f_1$ and $f_2$: the cone half-angles¶
Each cone is bounded by the lines that touch both the Sun and the Moon. The line and the two radii drawn to its points of contact form two right triangles with a common angle at the cone's vertex. Let $G = |\mathbf G|$ be the Sun–Moon distance, and $s_1$, $s_2$ the distances from the Moon's centre to the two vertices.
- Umbra. The vertex lies beyond the Moon, so the Moon's triangle is nested inside the Sun's: $$ \sin f_2 = \frac{k_2}{s_2} = \frac{R_\odot}{G + s_2} \quad\Longrightarrow\quad \sin f_2 = \frac{R_\odot - k_2}{G}. $$
- Penumbra. The bounding lines cross between the Sun and the Moon: $$ \sin f_1 = \frac{k_1}{s_1} = \frac{R_\odot}{G - s_1} \quad\Longrightarrow\quad \sin f_1 = \frac{R_\odot + k_1}{G}. $$
(From the first equality, $s_2 = k_2/\sin f_2$. Put it into the second and solve for $\sin f_2$; likewise for $f_1$.) Both angles are close to 16′, the Sun's apparent radius seen from the Moon: the edge of each cone is a ray from the Sun's limb that grazes the Moon's limb.
The diagram shows both constructions. It is not to scale.
# Similar triangles for the cone half-angles (schematic, not to scale)
R, k, G = 2.0, 0.8, 6.0 # Sun radius, Moon radius, Sun-Moon distance (Sun at 0, Moon at G)
fig, axes_ = plt.subplots(2, 1, figsize=(9, 6.4))
def right_angle(ax, T, a, b, size=0.18):
# small square at T between unit directions a and b
pts = np.array([T + size * a, T + size * (a + b), T + size * b])
ax.plot(pts[:, 0], pts[:, 1], color="k", lw=0.6)
for ax, kind in zip(axes_, ("umbra", "penumbra")):
ax.add_patch(plt.Circle((0, 0), R, color="#EE7733", alpha=0.3, lw=0))
ax.add_patch(plt.Circle((G, 0), k, color="#888888", alpha=0.45, lw=0))
ax.plot([-2.4, 12.2], [0, 0], color="0.5", lw=0.6, ls=":")
if kind == "umbra":
sin_f = (R - k) / G
s = k / sin_f # Moon centre to vertex, beyond the Moon
V = G + s
f = np.arcsin(sin_f)
n = np.array([np.sin(f), np.cos(f)])
T_sun, T_moon = np.array([0, 0]) + R * n, np.array([G, 0]) + k * n
line = np.array([T_sun - 0.8 * np.array([np.cos(f), -np.sin(f)]), [V + 1.2, -1.2 * np.tan(f)]])
name, sub, ks = "umbra", "2", "k_2"
title = (r"umbra: $\sin f_2 = k_2/s_2 = R_\odot/(G+s_2)$"
r"$\ \Rightarrow\ \sin f_2 = (R_\odot-k_2)/G$")
else:
sin_f = (R + k) / G
s = k / sin_f # Moon centre to vertex, towards the Sun
V = G - s
f = np.arcsin(sin_f)
n = np.array([np.sin(f), np.cos(f)])
T_sun, T_moon = R * n, np.array([G, 0]) - k * n
line = np.array([T_sun - 0.8 * np.array([np.cos(f), -np.sin(f)]), [12.0, -(12.0 - V) * np.tan(f)]])
name, sub, ks = "penumbra", "1", "k_1"
title = (r"penumbra: $\sin f_1 = k_1/s_1 = R_\odot/(G-s_1)$"
r"$\ \Rightarrow\ \sin f_1 = (R_\odot+k_1)/G$")
ax.plot(line[:, 0], line[:, 1], color="k", lw=1.1)
# radii to the points of contact, perpendicular to the line
ax.plot([0, T_sun[0]], [0, T_sun[1]], color="#CC3311", lw=1.0)
ax.plot([G, T_moon[0]], [0, T_moon[1]], color="#CC3311", lw=1.0)
along = np.array([np.cos(f), -np.sin(f)])
right_angle(ax, T_sun, -n, along)
right_angle(ax, T_moon, -np.sign(T_moon[1]) * n, along)
ax.text(*(0.5 * T_sun + [-0.45, 0.0]), r"$R_\odot$", fontsize=11, color="#CC3311")
ax.text(*(0.5 * (T_moon + [G, 0]) + [-0.42, 0.0]), f"${ks}$", fontsize=11, color="#CC3311")
# the angle f at the vertex
ang = np.degrees(f)
if kind == "umbra":
ax.add_patch(Arc((V, 0), 2.6, 2.6, theta1=180 - ang, theta2=180, lw=0.8))
ax.text(V - 1.6, 0.12, f"$f_{sub}$", fontsize=11)
else:
ax.add_patch(Arc((V, 0), 1.6, 1.6, theta1=180 - ang, theta2=180, lw=0.8))
ax.text(V - 1.15, 0.12, f"$f_{sub}$", fontsize=11)
ax.plot(V, 0, "k.", ms=5)
ax.text(V, -0.3, "vertex", ha="center", va="top", fontsize=9)
# distances G and s along the axis
yb = -2.35
ax.annotate("", xy=(0, yb), xytext=(G, yb), arrowprops=dict(arrowstyle="<->", lw=0.7))
ax.text(G / 2, yb - 0.1, "$G$", ha="center", va="top", fontsize=11)
ax.annotate("", xy=(G, yb + 0.6), xytext=(V, yb + 0.6), arrowprops=dict(arrowstyle="<->", lw=0.7))
ax.text((G + V) / 2, yb + 0.5, f"$s_{sub}$", ha="center", va="top", fontsize=11)
ax.text(0, 2.25, "Sun", ha="center", fontsize=10)
ax.text(G, k + 0.2, "Moon", ha="center", fontsize=10)
ax.set_title(title, fontsize=10)
ax.set_xlim(-2.5, 12.3)
ax.set_ylim(-3.0, 2.8)
ax.set_aspect("equal")
ax.axis("off")
plt.show()
Vertices, $l_1$ and $l_2$: the radii on the plane¶
The vertices are at $s_1 = k_1/\sin f_1$ towards the Sun and $s_2 = k_2/\sin f_2$ away from it, measured from the Moon's centre. Their heights above the fundamental plane are
$$ c_1 = z + \frac{k_1}{\sin f_1}, \qquad c_2 = z - \frac{k_2}{\sin f_2}. $$
A cone's radius at height $\zeta$ is the distance from the vertex times $\tan f$:
$$ L_1(\zeta) = (c_1 - \zeta)\tan f_1 = l_1 - \zeta\tan f_1, \qquad L_2(\zeta) = (c_2 - \zeta)\tan f_2 = l_2 - \zeta\tan f_2, $$
and on the fundamental plane itself ($\zeta = 0$):
$$ l_1 = c_1\tan f_1 = \left(z + \frac{k_1}{\sin f_1}\right)\tan f_1, \qquad l_2 = c_2\tan f_2 = \left(z - \frac{k_2}{\sin f_2}\right)\tan f_2 . $$
Sign convention. $L_1$ is always positive. $L_2$ is a signed radius:
- $L_2 < 0$: the observer is between the Moon and the umbral vertex ($\zeta > c_2$), inside the converging umbra. The eclipse is total, and the umbra's radius there is $|L_2|$.
- $L_2 > 0$: the observer is beyond the vertex, in the diverging antumbra. The eclipse is annular.
In 2024, $c_2$ was about −2.2: the vertex lay 14,100 km beyond the Earth's centre. Every point on the day side ($0 \le \zeta \le 1$) had $L_2 < 0$, so the eclipse was total along its whole path. When $c_2$ falls between 0 and 1 the sign can change along the path, and the eclipse is hybrid. NASA publishes $l_2$ with this sign, so a negative $l_2$ in a table of elements usually means a total eclipse.
In the code. elements_from (below) follows these formulas line by line. G is $\mathbf G$ and Gn is $G$. zhat, xhat and yhat are the unit vectors. x, y and z are the Moon's coordinates. Then come sin_f1, sin_f2, tan_f1 and tan_f2, and finally l1 and l2. raw_elements evaluates the elements at a list of times.
def elements_from(sun, moon):
# Besselian elements from ITRF Sun and Moon positions in Earth radii (d and mu in degrees)
G = sun - moon
Gn = np.linalg.norm(G, axis=1)
zhat = G / Gn[:, None]
d = np.arcsin(zhat[:, 2])
mu = np.unwrap(np.arctan2(-zhat[:, 1], zhat[:, 0]))
xhat = np.stack([np.sin(mu), np.cos(mu), np.zeros_like(mu)], axis=1)
yhat = np.stack([-np.sin(d) * np.cos(mu), np.sin(d) * np.sin(mu), np.cos(d)], axis=1)
x, y, z = (np.sum(v * moon, axis=1) for v in (xhat, yhat, zhat))
sin_f1 = (R_SUN + K1) / Gn
sin_f2 = (R_SUN - K2) / Gn
tan_f1 = sin_f1 / np.sqrt(1 - sin_f1**2)
tan_f2 = sin_f2 / np.sqrt(1 - sin_f2**2)
return {
"x": x, "y": y, "d": np.degrees(d), "mu": np.degrees(mu),
"l1": (z + K1 / sin_f1) * tan_f1, "l2": (z - K2 / sin_f2) * tan_f2,
"tan_f1": tan_f1, "tan_f2": tan_f2,
}
def raw_elements(times):
# Besselian elements at each time, from the JPL ephemerides
return elements_from(*sun_moon_itrf(times))
# The intermediate quantities at t0 = 2024-04-08 18:00 TT, NASA's reference epoch for this eclipse
t0 = sk.time.from_jd(2460409.25, sk.timescale.TT)
print(f"t0 = {t0} UTC (18:00:00 TT)")
sun, moon = sun_moon_itrf([t0])
Gn = np.linalg.norm(sun - moon)
zhat = (sun - moon)[0] / Gn
z = zhat @ moon[0]
sin_f1, sin_f2 = (R_SUN + K1) / Gn, (R_SUN - K2) / Gn
e0 = raw_elements([t0])
km = AE / 1e3
pd.DataFrame(
[("G, Sun-Moon distance", Gn), ("R_sun", R_SUN), ("z, Moon's height above the plane", z),
("s1 = k1 / sin f1", K1 / sin_f1), ("s2 = k2 / sin f2", K2 / sin_f2),
("c1, penumbral vertex", z + K1 / sin_f1), ("c2, umbral vertex", z - K2 / sin_f2),
("x", e0["x"][0]), ("y", e0["y"][0]), ("l1", e0["l1"][0]), ("l2", e0["l2"][0])],
columns=["quantity at t0", "Earth radii"],
).assign(km=lambda df: df["Earth radii"] * km).style.format(
{"Earth radii": "{:.6f}", "km": "{:,.1f}"}).hide(axis="index")
t0 = 2024-04-08T17:58:50.816000Z UTC (18:00:00 TT)
| quantity at t0 | Earth radii | km |
|---|---|---|
| G, Sun-Moon distance | 23433.606649 | 149,462,753.6 |
| R_sun | 109.122268 | 695,996.8 |
| z, Moon's height above the plane | 56.406806 | 359,770.3 |
| s1 = k1 / sin f1 | 58.370043 | 372,292.1 |
| s2 = k2 / sin f2 | 58.617608 | 373,871.1 |
| c1, penumbral vertex | 114.776849 | 732,062.5 |
| c2, umbral vertex | -2.210802 | -14,100.8 |
| x | -0.318249 | -2,029.8 |
| y | 0.219765 | 1,401.7 |
| l1 | 0.535817 | 3,417.5 |
| l2 | -0.010269 | -65.5 |
The umbral cone is $s_2 \approx 58.6$ Earth radii long (374,000 km), and the Moon was $z \approx 56.4$ Earth radii above the plane (360,000 km). The vertex therefore lay at $c_2 \approx -2.2$, beyond the Earth. On the fundamental plane the umbra's radius was $|l_2| \approx 65$ km. On the ground under the axis, about 0.9 Earth radii closer to the Moon, it was 93 km: $\zeta\tan f_2$ adds 27 km. The half-angles were $f_1 = 16.05'$ and $f_2 = 15.97'$.
Here are the two light-time corrections described earlier, measured by what they do to $x$ and $y$:
# Size of the light-time corrections: elements from geometric (uncorrected) positions at t0
q0 = sk.frametransform.qgcrf2itrf(t0)
geo_sun = (q0 * sk.jplephem.geocentric_pos(sk.solarsystem.Sun, t0) / AE)[None, :]
geo_moon = (q0 * sk.jplephem.geocentric_pos(sk.solarsystem.Moon, t0) / AE)[None, :]
full = elements_from(sun, moon)
for label, e in [("Sun: light time and aberration", elements_from(geo_sun, moon)),
("Moon: light time", elements_from(sun, geo_moon))]:
dx, dy = (full["x"] - e["x"])[0] * km, (full["y"] - e["y"])[0] * km
print(f"{label:32s} moves the shadow axis by ({dx:+.1f}, {dy:+.1f}) km, {np.hypot(dx, dy):.1f} km")
Sun: light time and aberration moves the shadow axis by (+33.1, +13.5) km, 35.8 km Moon: light time moves the shadow axis by (-1.2, -0.6) km, 1.3 km
Fitting the elements with polynomials¶
We could evaluate raw_elements every time the elements are needed. Instead, like the almanacs, we evaluate them every 10 minutes over six hours centred on $t_0$ (37 samples). We then fit each of $x$, $y$, $d$, $l_1$, $l_2$ and $\mu$ with a least-squares cubic polynomial in $t - t_0$, in hours. There are four reasons:
- The elements are smooth. Over six hours the Moon covers a short arc of its orbit, and the Sun barely moves. $x$ and $y$ are nearly linear, and $\mu$ is nearly linear. The cubic term of $x$ contributes only about 1.5 km at ±3 h, and the fit residuals printed below show that everything beyond the cubic adds less than 0.5 m.
- Speed. The mapping to the ground evaluates the elements a very large number of times. The duration map alone runs six Newton iterations at each of 300,000 places, each at its own time. That is 1.8 million evaluations. Each call to
raw_elementscosts about 20 µs per instant (two ephemeris look-ups and an Earth-orientation rotation), which adds up to half a minute. The polynomials take milliseconds. - Derivatives for free. The limits and the duration need the rates $\dot x$, $\dot y$, $\dot d$ and $\dot\mu$. The polynomials differentiate exactly.
- The published form. NASA and the almanacs publish their elements as polynomials in the same form. The same code therefore runs on NASA's coefficients, which we use below to test our mapping independently of our ephemeris.
$\tan f_1$ and $\tan f_2$ change by less than one part in 10,000 over the six hours, as the Sun–Moon distance changes. As NASA does, we treat them as constants, the means over the window. The largest error this causes in $L_2$ is about a metre.
Besselian.fit does the sampling and fitting. Calling a Besselian object at hours th returns the elements and their rates per hour (x_dot, mu_dot, ...).
class Besselian:
# Cubic polynomials in t - t0 (hours, TT) for x, y, d, l1, l2 and mu
names = ["x", "y", "d", "l1", "l2", "mu"]
P = np.polynomial.polynomial
def __init__(self, t0, coef, tan_f1, tan_f2):
self.t0, self.coef, self.tan_f1, self.tan_f2 = t0, coef, tan_f1, tan_f2
@classmethod
def fit(cls, t0, span_hours=3.0, step_hours=1 / 6):
th = np.arange(-span_hours, span_hours + 1e-9, step_hours)
e = raw_elements([t0 + sk.duration.from_hours(h) for h in th])
coef = {k: cls.P.polyfit(th, e[k], 3) for k in cls.names}
coef["mu"][0] %= 360
return cls(t0, coef, e["tan_f1"].mean(), e["tan_f2"].mean())
def __call__(self, th):
# Elements and their rates (per hour) at hours th from t0
out = {k: self.P.polyval(th, c) for k, c in self.coef.items()}
out.update({k + "_dot": self.P.polyval(th, self.P.polyder(c)) for k, c in self.coef.items()})
return out
def hours(self, t):
# Hours from t0 of a satkit time
return (t - self.t0).hours
def time(self, th):
# satkit time at hours th from t0
return self.t0 + sk.duration.from_hours(float(th))
bes = Besselian.fit(t0)
# Fit residuals, halfway between the samples
th_mid = np.arange(-3 + 1 / 12, 3, 1 / 6)
exact, fitted = raw_elements([bes.time(h) for h in th_mid]), bes(th_mid)
print("largest fit residual:")
for k in ("x", "y", "l1", "l2"):
print(f" {k:3s} {np.max(np.abs(exact[k] - fitted[k])) * AE:.3f} m")
for k in ("d", "mu"):
r = (exact[k] - fitted[k] + 180) % 360 - 180
print(f" {k:3s} {np.max(np.abs(r)) * 3600:.1e} arcsec")
largest fit residual: x 0.079 m y 0.353 m l1 0.010 m l2 0.010 m d 7.5e-04 arcsec mu 1.7e-05 arcsec
What the elements look like¶
The figure shows the fundamental plane face on, as seen from the Sun. It uses the fitted elements. The Earth's disk and graticule are drawn as they were at $t_0$, and only the sunward hemisphere is shown. The outline is an ellipse with semi-axes 1 and $\rho_1 = \sqrt{1 - e^2\cos^2 d}$, the Explanatory Supplement's name for it, where $e^2 = f(2-f)$. The shadow axis runs along the straight red line during the six hours, and the umbra is the dot at its greatest-eclipse position. The large circle is the penumbra at the same instant. Meanwhile the graticule turns about the projected pole by 15° per hour. Mapping the straight track onto that turning globe is the job of the next section.
plane_axes builds the three unit vectors from $d$ and $\mu$. The mapping code uses it from here on.
def plane_axes(e):
# ITRF unit vectors of the fundamental-plane axes, shape (..., 3)
d, mu = np.radians(e["d"]), np.radians(e["mu"])
xhat = np.stack([np.sin(mu), np.cos(mu), np.zeros_like(mu)], axis=-1)
yhat = np.stack([-np.sin(d) * np.cos(mu), np.sin(d) * np.sin(mu), np.cos(d)], axis=-1)
zhat = np.stack([np.cos(d) * np.cos(mu), -np.cos(d) * np.sin(mu), np.sin(d)], axis=-1)
return xhat, yhat, zhat
def geodetic_grid(lon_deg, lat_deg):
# ITRF positions (m) of a lat/lon grid on the ellipsoid
lon, lat = np.meshgrid(np.radians(lon_deg), np.radians(lat_deg))
e2 = FLAT * (2 - FLAT)
N = AE / np.sqrt(1 - e2 * np.sin(lat) ** 2)
return np.stack([N * np.cos(lat) * np.cos(lon), N * np.cos(lat) * np.sin(lon),
N * (1 - e2) * np.sin(lat)], axis=-1)
th_ge = 0.29 # about greatest eclipse, 18:17 TT
e_t0, e_ge = bes(np.array([0.0])), bes(np.array([th_ge]))
axes_t0 = [v[0] for v in plane_axes(e_t0)]
def on_plane(p):
# (xi, eta) of ITRF points (m) at t0, NaN on the far side of the Earth
xi, eta, zeta = (p / AE @ v for v in axes_t0)
return np.where(zeta >= 0, xi, np.nan), np.where(zeta >= 0, eta, np.nan)
fig, ax = plt.subplots(figsize=(8.5, 6.4))
q = np.linspace(0, 2 * np.pi, 361)
rho1 = np.sqrt(1 - FLAT * (2 - FLAT) * np.cos(np.radians(e_t0["d"][0])) ** 2)
ax.fill(np.cos(q), rho1 * np.sin(q), color="#0077BB", alpha=0.1, lw=0)
ax.plot(np.cos(q), rho1 * np.sin(q), color="#0077BB", lw=0.8)
fine = np.linspace(-180, 180, 721)
for lat in range(-60, 61, 30):
ax.plot(*on_plane(geodetic_grid(fine, [lat])[0]), color="#0077BB", lw=1.0 if lat == 0 else 0.4)
for lon in range(-180, 180, 30):
ax.plot(*on_plane(geodetic_grid([lon], np.linspace(-90, 90, 361))[:, 0]), color="#0077BB", lw=0.4)
xi, eta = on_plane(geodetic_grid([lon], [3])[0, 0])
if np.isfinite(xi) and abs(xi) < 0.9:
ax.text(xi, eta, f"{abs(lon)}°{'W' if lon < 0 else 'E'}", fontsize=7, color="#0077BB", ha="center")
ax.plot(*on_plane(geodetic_grid([0], [90])[0, 0]), "o", color="#0077BB", ms=3)
ax.text(*np.array(on_plane(geodetic_grid([0], [90])[0, 0])) + [0.04, 0.03], "N pole", fontsize=8, color="#0077BB")
th_track = np.linspace(-3, 3, 121)
track = bes(th_track)
ax.plot(track["x"], track["y"], color="#CC3311", lw=1.3, label="shadow axis $(x, y)$, ±3 h")
for h in range(-3, 4):
eh = bes(np.array([float(h)]))
ax.plot(eh["x"], eh["y"], "o", color="#CC3311", ms=3)
ax.text(eh["x"][0] + 0.05, eh["y"][0] - 0.1, bes.time(h).strftime("%H:%M"), fontsize=8, color="#CC3311")
ax.plot(e_ge["x"][0] + e_ge["l1"][0] * np.cos(q), e_ge["y"][0] + e_ge["l1"][0] * np.sin(q), color="0.35",
lw=0.8, ls="--", label="penumbra, radius $l_1$")
ax.fill(e_ge["x"][0] + abs(e_ge["l2"][0]) * np.cos(q), e_ge["y"][0] + abs(e_ge["l2"][0]) * np.sin(q),
color="k", label=r"umbra, radius $|l_2|$ (to scale)")
for v, label in [((0.35, 0), r"$\xi$"), ((0, 0.35), r"$\eta$")]:
ax.annotate("", xy=v, xytext=(0, 0), arrowprops=dict(arrowstyle="-|>", lw=1.0, color="k"))
ax.text(v[0] + 0.03, v[1] + 0.03, label, fontsize=11)
ax.set_xlim(-2.0, 1.7)
ax.set_ylim(-1.2, 1.4)
ax.set_aspect("equal")
ax.set_xlabel(r"$\xi$ (Earth radii, east)")
ax.set_ylabel(r"$\eta$ (Earth radii, north)")
ax.legend(loc="lower right", fontsize=8.5, framealpha=0.9)
ax.set_title("The fundamental plane seen from the Sun, 8 April 2024 (times UTC)", fontsize=10)
plt.show()
Comparison with NASA's elements¶
NASA publishes polynomial Besselian elements for each eclipse (Fred Espenak, eclipse.gsfc.nasa.gov, SE2024Apr08T Besselian elements). The coefficients are copied below, with $\tan f_1 = 0.0046683$ and $\tan f_2 = 0.0046450$. NASA computed them from the VSOP87 and ELP2000-85 analytical theories of the Sun and Moon (Bretagnon & Francou 1988; Chapront-Touzé & Chapront 1988), not from the JPL numerical ephemerides.
Time scales. NASA's polynomials use Terrestrial Dynamical Time, $t_0$ = 18:00 TT, which is also our $t_0$. As explained above, the Earth's rotation follows UT1, and NASA had to assume $\Delta T$ = 70.6 s. Following the almanac convention, NASA's $\mu$ is the ephemeris hour angle, computed as if UT1 were TT. The Greenwich hour angle is
$$ \mu_{\rm Greenwich} = \mu - 1.002738 \times 15''\,{\rm s}^{-1} \times \Delta T = \mu - 0.00417807^\circ \times \Delta T[{\rm s}] .$$
The factor 1.002738 converts seconds of time into seconds of Earth rotation: the Earth turns 360° in a sidereal day, about 4 minutes shorter than a solar day.
In April 2024 the actual $\Delta T$ was 69.2 s. NASA's value was 1.4 s too large, which moves its path about 0.5 km east. The function delta_t gets $\Delta T$ from satkit: TT − UTC from the leap-second table, minus UT1 − UTC from the IERS. Below we convert NASA's $\mu$ with this measured $\Delta T$, so that the comparison shows only the ephemeris differences. The last column then shows the effect of NASA's own $\Delta T$. Lengths are converted to kilometres ($\times\,a_e$).
# NASA / Espenak polynomial Besselian elements, 2024 Apr 08, t0 = 18:00 TDT, Delta T = 70.6 s
# "Eclipse Predictions by Fred Espenak, NASA's GSFC"
# https://eclipse.gsfc.nasa.gov/SEbeselm/SEbeselm2001/SE2024Apr08Tbeselm.html
NASA_COEF = {
"x": [-0.318157, 0.5117105, 0.0000326, -0.0000085],
"y": [0.219747, 0.2709586, -0.0000594, -0.0000047],
"d": [7.58620, 0.014844, -0.000002, 0.0],
"l1": [0.535813, 0.0000618, -0.0000128, 0.0],
"l2": [-0.010274, 0.0000615, -0.0000127, 0.0],
"mu": [89.59122, 15.004084, 0.0, 0.0],
}
NASA_TAN_F1, NASA_TAN_F2 = 0.0046683, 0.0046450
NASA_DELTA_T = 70.6
def delta_t(t):
# TT - UT1 in seconds, from satkit's leap seconds and the IERS UT1 - UTC
tt_minus_utc = (t.to_mjd(sk.timescale.TT) - t.to_mjd(sk.timescale.UTC)) * 86400
return tt_minus_utc - sk.frametransform.earth_orientation_params(t)[0]
def nasa_besselian(delta_t_seconds):
# NASA's elements with mu converted to the Greenwich hour angle for a given Delta T
coef = {k: np.array(v, dtype=float) for k, v in NASA_COEF.items()}
coef["mu"][0] -= 0.00417807 * delta_t_seconds
return Besselian(t0, coef, NASA_TAN_F1, NASA_TAN_F2)
dT = delta_t(t0)
print(f"Delta T from satkit: {dT:.2f} s (NASA assumed {NASA_DELTA_T} s)")
nasa_measured_dt = nasa_besselian(dT)
nasa_own_dt = nasa_besselian(NASA_DELTA_T)
th = np.arange(-3.0, 3.01, 1.0)
ours, theirs, theirs_own = bes(th), nasa_measured_dt(th), nasa_own_dt(th)
df = pd.DataFrame({
"TT": [f"{18 + h:02.0f}:00" for h in th],
"x (km)": (ours["x"] - theirs["x"]) * km,
"y (km)": (ours["y"] - theirs["y"]) * km,
"d (arcsec)": (ours["d"] - theirs["d"]) * 3600,
"mu (arcsec)": (ours["mu"] - theirs["mu"]) * 3600,
"l1 (km)": (ours["l1"] - theirs["l1"]) * km,
"l2 (km)": (ours["l2"] - theirs["l2"]) * km,
"mu, NASA's Delta T (arcsec)": (ours["mu"] - theirs_own["mu"]) * 3600,
})
print(f"tan f1: satkit {bes.tan_f1:.7f}, NASA {NASA_TAN_F1:.7f}")
print(f"tan f2: satkit {bes.tan_f2:.7f}, NASA {NASA_TAN_F2:.7f}")
print("satkit minus NASA:")
df.style.format(precision=3).hide(axis="index")
Delta T from satkit: 69.20 s (NASA assumed 70.6 s) tan f1: satkit 0.0046683, NASA 0.0046683 tan f2: satkit 0.0046451, NASA 0.0046450 satkit minus NASA:
| TT | x (km) | y (km) | d (arcsec) | mu (arcsec) | l1 (km) | l2 (km) | mu, NASA's Delta T (arcsec) |
|---|---|---|---|---|---|---|---|
| 15:00 | -0.621 | 0.080 | -0.295 | 0.018 | 0.027 | 0.030 | 21.067 |
| 16:00 | -0.606 | 0.097 | -0.360 | 0.015 | 0.026 | 0.030 | 21.064 |
| 17:00 | -0.596 | 0.108 | -0.400 | 0.005 | 0.025 | 0.030 | 21.054 |
| 18:00 | -0.588 | 0.114 | -0.416 | -0.012 | 0.025 | 0.030 | 21.037 |
| 19:00 | -0.578 | 0.118 | -0.405 | -0.034 | 0.025 | 0.029 | 21.015 |
| 20:00 | -0.563 | 0.119 | -0.369 | -0.063 | 0.026 | 0.029 | 20.986 |
| 21:00 | -0.540 | 0.120 | -0.306 | -0.096 | 0.027 | 0.030 | 20.953 |
Reading the table. Each row is a whole hour of TT across the six-hour window. Each column is satkit's element minus NASA's. For $x$, $y$, $l_1$ and $l_2$ the difference is a length on the fundamental plane, in km. When the Sun is high, a shift on the plane is roughly the same shift on the ground; when the Sun is low, the ground shift is larger, by $1/\sin(\text{altitude})$ in the direction away from the Sun. For $d$ and $\mu$ the difference is an angle, in arcseconds. One arcsecond of $\mu$ moves the path one arcsecond of longitude, 31 m × cos(latitude).
- $x$ and $y$, the position of the shadow axis, differ by about 0.6 km, almost constant over the six hours. That is 0.3″ as seen from the Earth, the expected level of disagreement between the ELP2000-85 lunar theory and the JPL DE440 ephemeris for 2024.
- $d$ differs by about 0.4″. Our frame is the ITRF, whose pole differs from the celestial pole NASA uses by the polar motion, about 0.36″ on that day. This is a difference of reference, not of accuracy: 0.4″ of $d$ is about 10 m on the ground.
- $\mu$ agrees to a few hundredths of an arcsecond once both use the measured $\Delta T$. With NASA's $\Delta T$ (last column) the difference is 21″: the 1.4 s of Earth rotation, or 0.5 km of longitude at the path's latitudes.
- $l_1$, $l_2$, $\tan f_1$ and $\tan f_2$ agree to 30 m and to the last digit NASA gives. The constants match, and the shadow cones are the same.
From the plane to the ground¶
So far we have described the shadow on the fundamental plane. Now we bring in the Earth: where a given place lies relative to that plane, and how it moves as the Earth turns.
Observer coordinates¶
Take a place at geodetic latitude $\varphi$, longitude $\lambda$ (east positive) and height $h$ on the WGS84 ellipsoid. Its geocentric position, in Earth radii, is conventionally written
$$ \mathbf p = (\rho\cos\varphi'\cos\lambda,\ \rho\cos\varphi'\sin\lambda,\ \rho\sin\varphi'), $$
where $\rho$ is the distance from the Earth's centre and $\varphi'$ is the geocentric latitude. The two products follow from the ellipsoid, using the reduced latitude $u$, with $\tan u = (1-f)\tan\varphi$:
$$ \rho\sin\varphi' = (1-f)\sin u + \frac{h}{a_e}\sin\varphi, \qquad \rho\cos\varphi' = \cos u + \frac{h}{a_e}\cos\varphi . $$
Projecting $\mathbf p$ onto the three unit vectors gives its fundamental-plane coordinates. With the hour angle of the shadow axis at the place,
$$ H = \mu + \lambda, $$
they are
$$ \begin{aligned} \xi &= \rho\cos\varphi'\,\sin H, \\ \eta &= \rho\sin\varphi'\,\cos d - \rho\cos\varphi'\,\cos H\,\sin d, \\ \zeta &= \rho\sin\varphi'\,\sin d + \rho\cos\varphi'\,\cos H\,\cos d . \end{aligned} $$
Sign convention. $\mu$ is measured westward from Greenwich and $\lambda$ eastward, so $H = \mu + \lambda$ is the local hour angle, measured westward. $H = 0$ means the shadow axis is on the place's meridian. Meeus and the Explanatory Supplement count longitude positive west, and so write $H = \mu - \lambda$; it is the same angle.
The code does not use these formulas directly. It computes $\mathbf p$ once (geodetic_grid, or satkit.itrfcoord(...).vector) and takes dot products with the unit vectors from plane_axes. That is the same thing, as the cell below checks for Dallas.
Why $\zeta$ matters. $\zeta$ is the place's height above the fundamental plane, towards the Moon. At that height the umbra's radius is not $l_2$ but
$$ L_2 = l_2 - \zeta\tan f_2 . $$
$\tan f_2 \approx 0.0046$ Earth radii per Earth radius, or 29.6 km per 6378 km. So a place with the Sun overhead ($\zeta \approx 1$) sits in an umbra about 30 km larger in radius than the one on the plane. Since $\zeta \approx \sin(\text{Sun's altitude})$, it also tells us whether the Sun is up: $\zeta > 0$ on the day side.
The observer moves. The Earth turns, so a fixed place moves on the fundamental plane. Differentiating $\xi$ and $\eta$ with respect to time, with $\mu$ and $d$ changing, gives
$$ \dot\xi = \dot\mu\,\rho\cos\varphi'\cos H = \dot\mu\,(\zeta\cos d - \eta\sin d), \qquad \dot\eta = \dot\mu\,\xi\sin d - \dot d\,\zeta, $$
with $\dot\mu$ and $\dot d$ in radians per hour. On the day side ($\zeta > 0$) the place moves east ($\dot\xi > 0$), at up to $\dot\mu \approx 0.26$ Earth radii per hour: 1,670 km/h at the equator. The shadow moves east too, at about 0.58 Earth radii per hour, so the Earth's rotation slows the shadow down relative to the ground. This is observer_rates.
The central line¶
The central line is where the shadow axis meets the ground. We need the point with $\xi = x$ and $\eta = y$ whose $\zeta$ puts it on the ellipsoid, on the side facing the Sun. The point is
$$ \mathbf p(\zeta) = x\,\hat{\boldsymbol\xi} + y\,\hat{\boldsymbol\eta} + \zeta\,\hat{\boldsymbol\zeta}, $$
and the ellipsoid is $X^2 + Y^2 + Z^2/(1-f)^2 = 1$ (Earth radii, ITRF). Stretching the $Z$ axis by $1/(1-f)$, with $S = {\rm diag}(1, 1, 1/(1-f))$, turns the ellipsoid into the unit sphere. The condition $|S\,\mathbf p(\zeta)|^2 = 1$ is then a quadratic in $\zeta$. Write $\mathbf A = S(x\hat{\boldsymbol\xi} + y\hat{\boldsymbol\eta})$ and $\mathbf B = S\hat{\boldsymbol\zeta}$:
$$ (\mathbf B\cdot\mathbf B)\,\zeta^2 + 2(\mathbf A\cdot\mathbf B)\,\zeta + (\mathbf A\cdot\mathbf A - 1) = 0, \qquad \zeta = \frac{-\mathbf A\cdot\mathbf B + \sqrt{(\mathbf A\cdot\mathbf B)^2 - (\mathbf B\cdot\mathbf B)(\mathbf A\cdot\mathbf A - 1)}}{\mathbf B\cdot\mathbf B} . $$
The larger root is the intersection facing the Sun; the smaller one is on the night side. A negative discriminant means the line misses the Earth. The function returns NaN there, and those instants drop out of the path. For a sphere ($f = 0$) this reduces to the textbook $\zeta = \sqrt{1 - \xi^2 - \eta^2}$. The Explanatory Supplement handles the flattening with its auxiliary quantities $\rho_1$, $\rho_2$, $d_1$ and $d_2$; the stretched quadratic is the same calculation in vector form. In the code, surface_zeta solves it (qa, qb and qc are $\mathbf B\cdot\mathbf B$, $\mathbf A\cdot\mathbf B$ and $\mathbf A\cdot\mathbf A - 1$), and central_line applies it to $(x, y)$.
Back to latitude and longitude. to_itrf turns $(\xi, \eta, \zeta)$ back into an ITRF vector in metres. The longitude is $\lambda = \operatorname{atan2}(p_Y, p_X)$. The latitude needs care. The angle $\arctan\big(p_Z/\sqrt{p_X^2+p_Y^2}\big)$ is the geocentric latitude $\varphi'$, while maps use the geodetic latitude $\varphi$: the angle between the ellipsoid's normal and the equator. They differ by up to 11.5′ at 45°, or 21 km on the ground, far more than any other error here. On the ellipsoid ($h = 0$) the conversion is exact and closed-form,
$$ \tan\varphi = \frac{\tan\varphi'}{(1-f)^2}, $$
and satkit.itrfcoord, which latlon uses, handles any height with Bowring's iteration (Bowring 1976).
STRETCH = np.array([1.0, 1.0, 1.0 / (1.0 - FLAT)]) # maps the WGS84 ellipsoid to the unit sphere
def surface_zeta(xi, eta, axes):
# zeta at which (xi, eta) meets the ellipsoid on the sunward side; NaN off the limb
xhat, yhat, zhat = axes
a = (xi[..., None] * xhat + eta[..., None] * yhat) * STRETCH
b = zhat * STRETCH
qa, qb, qc = np.sum(b * b, -1), np.sum(a * b, -1), np.sum(a * a, -1) - 1
with np.errstate(invalid="ignore"):
return (-qb + np.sqrt(qb**2 - qa * qc)) / qa
def observer_rates(e, xi, eta, zeta):
# Velocity of a point fixed on the Earth, in fundamental-plane coordinates per hour
d = np.radians(e["d"])
mu_dot, d_dot = np.radians(e["mu_dot"]), np.radians(e["d_dot"])
return mu_dot * (zeta * np.cos(d) - eta * np.sin(d)), mu_dot * xi * np.sin(d) - d_dot * zeta
def to_itrf(xi, eta, zeta, axes):
xhat, yhat, zhat = axes
return (xi[..., None] * xhat + eta[..., None] * yhat + zeta[..., None] * zhat) * AE
def central_line(bes, th):
e = bes(th)
axes = plane_axes(e)
return to_itrf(e["x"], e["y"], surface_zeta(e["x"], e["y"], axes), axes)
def latlon(p):
# Geodetic latitude and longitude (deg) of ITRF points; NaN rows stay NaN
flat = p.reshape(-1, 3)
out = np.full((len(flat), 2), np.nan)
ok = np.all(np.isfinite(flat), axis=1)
if ok.any():
out[ok] = [(c.latitude_deg, c.longitude_deg) for c in map(sk.itrfcoord, flat[ok])]
return out.reshape(*p.shape[:-1], 2)
# Observer coordinates two ways, for Dallas at 18:17 TT
lat, lon = np.radians(32.7767), np.radians(-96.7970)
u = np.arctan((1 - FLAT) * np.tan(lat))
rho_sin, rho_cos = (1 - FLAT) * np.sin(u), np.cos(u)
d, H = np.radians(e_ge["d"][0]), np.radians(e_ge["mu"][0]) + lon
textbook = (rho_cos * np.sin(H), rho_sin * np.cos(d) - rho_cos * np.cos(H) * np.sin(d),
rho_sin * np.sin(d) + rho_cos * np.cos(H) * np.cos(d))
p_dallas = sk.itrfcoord(latitude_deg=32.7767, longitude_deg=-96.7970).vector / AE
dots = [p_dallas @ v[0] for v in plane_axes(e_ge)]
print(f"hour angle H = {np.degrees(H):.3f} deg")
print("(xi, eta, zeta) from the formulas: ", np.round(textbook, 9))
print("(xi, eta, zeta) from the dot products: ", np.round(dots, 9))
print(f"shadow axis (x, y) = ({e_ge['x'][0]:.6f}, {e_ge['y'][0]:.6f})")
hour angle H = -3.144 deg (xi, eta, zeta) from the formulas: [-0.04615475 0.42255259 0.90408317] (xi, eta, zeta) from the dot products: [-0.04615475 0.42255259 0.90408317] shadow axis (x, y) = (-0.169850, 0.298338)
The two agree to every digit. At 18:17 TT Dallas was 0.18 Earth radii (1,100 km) north-east of the shadow axis on the plane, and $\zeta = 0.90$ put the Sun about 65° above its horizon. The shadow reached Dallas about 25 minutes later: totality there began at 18:40:43 UTC (the table of contacts below).
The northern and southern limits¶
The tangency condition. Follow one place through the eclipse, from its own point of view. It is fixed on the ground, and the umbra, a circle of radius $|L_2|$ about $(x, y)$, drifts past it with the relative velocity
$$ (a, b) = (\dot x - \dot\xi,\ \dot y - \dot\eta), \qquad n = \sqrt{a^2 + b^2}. $$
Let $(u, v) = (x - \xi,\ y - \eta)$ be the offset of the shadow axis from the place. The squared distance $u^2 + v^2$ is smallest when its time derivative $2(ua + vb)$ vanishes, that is, when $(u, v)$ is perpendicular to $(a, b)$. A place is on a limit when, at that moment of closest approach, it is exactly on the edge of the umbra: $|(u, v)| = |L_2|$. A place a little farther out is never inside the umbra; a place a little farther in is. The limit points at time $t$ are therefore the two points of the umbral circle whose radius is perpendicular to the relative motion:
$$ \xi = x \mp |L_2|\,\frac{b}{n}, \qquad \eta = y \pm |L_2|\,\frac{a}{n} \qquad \text{(upper signs: northern limit).} $$
Geometrically, the limits are the envelope of the moving circles. A disk sliding in a direction sweeps a strip whose edges are offset by its radius perpendicular to that direction. That is why the limit points are not simply $(x, y \pm |L_2|)$, the top and bottom of the circle. The track is inclined, so the top and bottom of the circle lie inside the strip, at distance $|L_2|\cos\alpha$ from its centre line, where $\alpha$ is the track's inclination. It is also why the motion must be taken relative to the observer. The strip that matters is the one swept over the turning Earth, and the Earth's rotation changes the direction of motion: $(\dot\xi, \dot\eta)$ points mostly east, so subtracting it tilts the track.
The figure uses the actual numbers at $t_0$, in kilometres on the plane and relative to the shadow axis. The black circle is the umbra at $t_0$, and the grey circles are the umbra at earlier and later instants, as seen from the ground. The umbra's radius on the ground is $|L_2| \approx 93$ km. On the plane the shadow moves 28° north of east, but relative to the turning Earth it moves 45° north of east. The naive top and bottom points would put the limits only 65 km from the central line instead of 93 km. Using the shadow's own motion on the plane, without the Earth's rotation, would put them at 89 km, and in the wrong places along the path.
# The tangency condition at t0, drawn with the actual numbers (km, relative to the shadow axis)
e = bes(np.array([0.0]))
axes = plane_axes(e)
zeta = surface_zeta(e["x"], e["y"], axes)
xi_dot, eta_dot = observer_rates(e, e["x"], e["y"], zeta)
a, b = (e["x_dot"] - xi_dot)[0], (e["y_dot"] - eta_dot)[0]
L = abs(e["l2"] - zeta * bes.tan_f2)[0] * km
alpha = np.arctan2(b, a) # direction of the motion relative to the ground
beta = np.arctan2(e["y_dot"][0], e["x_dot"][0]) # direction of the shadow's own motion on the plane
along, perp = np.array([np.cos(alpha), np.sin(alpha)]), np.array([-np.sin(alpha), np.cos(alpha)])
print(f"|L2| on the ground {L:.1f} km; relative motion {np.degrees(alpha):.1f} deg, "
f"shadow's own motion {np.degrees(beta):.1f} deg north of east")
print(f"distance from the central line: tangent points {L:.1f} km, top/bottom of the circle "
f"{L * np.cos(alpha):.1f} km, tangent points for the shadow's own motion {L * np.cos(alpha - beta):.1f} km")
fig, ax = plt.subplots(figsize=(8, 6))
q = np.linspace(0, 2 * np.pi, 200)
for s in np.arange(-240, 241, 80):
c = s * along
ax.plot(c[0] + L * np.cos(q), c[1] + L * np.sin(q), color="k" if s == 0 else "0.7", lw=1.1 if s == 0 else 0.6)
ends = np.array([-340, 340])[:, None] * along
ax.plot(*ends.T, color="#CC3311", lw=1.2, label="central line (track of the axis)")
for sign in (1, -1):
ax.plot(*(ends + sign * L * perp).T, color="#4a5878", lw=1.2, label="limits" if sign == 1 else None)
ax.plot(*(sign * L * perp), "o", color="#4a5878", ms=6, label="tangent points" if sign == 1 else None)
ax.plot(0, sign * L, "x", color="#EE7733", ms=8, mew=2, label="top and bottom of the circle" if sign == 1 else None)
ax.annotate("", xy=L * perp, xytext=(0, 0), arrowprops=dict(arrowstyle="->", lw=0.9))
ax.text(*(0.5 * L * perp + [-6, -6]), "$|L_2|$", fontsize=11, ha="right", va="top")
ax.annotate("", xy=150 * along, xytext=(0, 0), arrowprops=dict(arrowstyle="-|>", lw=1.6, color="#CC3311"))
ax.text(*(150 * along + [8, -18]), "motion relative\nto the ground $(a, b)$", fontsize=9, color="#CC3311")
bdir = np.array([np.cos(beta), np.sin(beta)])
ax.annotate("", xy=150 * bdir, xytext=(0, 0), arrowprops=dict(arrowstyle="-|>", lw=1.0, ls="--", color="0.4"))
ax.text(*(150 * bdir + [8, -4]), "shadow's own motion\non the plane", fontsize=9, color="0.4", va="top")
ax.set_xlim(-260, 330)
ax.set_ylim(-230, 250)
ax.set_aspect("equal")
ax.set_xlabel(r"$\xi - x$ (km, east)")
ax.set_ylabel(r"$\eta - y$ (km, north)")
ax.legend(loc="lower right", fontsize=8.5, framealpha=0.95)
plt.show()
|L2| on the ground 92.8 km; relative motion 45.2 deg, shadow's own motion 27.9 deg north of east distance from the central line: tangent points 92.8 km, top/bottom of the circle 65.4 km, tangent points for the shadow's own motion 88.6 km
The iteration. $L_2$ and the observer's velocity $(\dot\xi, \dot\eta)$ depend on the limit point's own $\zeta$, which is not known until the point is found. limit_lines therefore iterates. It starts from the central line's $\zeta$, computes the limit point, finds that point's $\zeta$ on the ellipsoid with surface_zeta, and repeats. Each limit point belongs to an instant: the moment at which that place sees the umbra's edge pass over it, which is also its moment of maximum eclipse. NASA's tables are organised the same way. This is the construction of the Explanatory Supplement (Urban & Seidelmann 2013) and of Meeus (1989), written with vectors.
(Strictly, $L_2$ also changes slowly with time as the place's $\zeta$ changes, which tilts the tangency condition slightly. The effect is to slide the point about 0.2 km along the path, which changes its time label by a fraction of a second but leaves the line where it is.)
Convergence, and sunrise and sunset. The cell below compares 2 to 8 iterations with a 30-iteration reference, on points every 10 s along both limits. Over most of the path two iterations are already good to a metre. The slow cases are near sunrise and sunset. There the Sun is on the horizon and the ground is seen edge-on from the Sun, so a small step in $(\xi, \eta)$ changes $\zeta$, and with it $L_2$ and the observer's velocity, a great deal. Even there each iteration reduces the error about sevenfold, and eight iterations are good to about 0.2 m.
At the very ends of the path the limit point can fall beyond the Earth's limb. There the quadratic has no real root, surface_zeta returns NaN, and the point drops out. That is why the northern and southern limits begin and end at slightly different times from the central line. The short curves that close the path at its two ends, where the umbra's edge meets the sunrise and sunset terminator, are not drawn here. Near those curves the Sun is on the horizon, where atmospheric refraction, which we ignore, also matters.
umbra_outline uses the same iteration to draw the whole umbral circle on the ground at one instant: the dashed outlines on the map.
def limit_lines(bes, th, niter=8):
# Northern and southern limits of the umbra at hours th (ITRF, m; NaN where off the Earth)
e = bes(th)
axes = plane_axes(e)
out = {}
for side, sign in (("north", 1), ("south", -1)):
xi, eta = e["x"], e["y"]
zeta = np.nan_to_num(surface_zeta(xi, eta, axes))
for _ in range(niter):
xi_dot, eta_dot = observer_rates(e, xi, eta, zeta)
a, b = e["x_dot"] - xi_dot, e["y_dot"] - eta_dot
n = np.hypot(a, b)
L2 = np.abs(e["l2"] - zeta * bes.tan_f2)
xi, eta = e["x"] - sign * L2 * b / n, e["y"] + sign * L2 * a / n
zeta = surface_zeta(xi, eta, axes)
out[side] = to_itrf(xi, eta, zeta, axes)
return out
def umbra_outline(bes, th, npts=181, niter=6):
# Outline of the umbra on the ground at one instant (ITRF, m)
e = bes(np.array([th]))
axes = plane_axes(e)
q = np.linspace(0, 2 * np.pi, npts)
zeta = np.zeros_like(q)
for _ in range(niter):
L2 = np.abs(e["l2"] - np.nan_to_num(zeta) * bes.tan_f2)
xi, eta = e["x"] + L2 * np.cos(q), e["y"] + L2 * np.sin(q)
zeta = surface_zeta(xi, eta, axes)
return to_itrf(xi, eta, zeta, axes)
# Convergence of the limit iteration, against a 30-iteration reference
th_check = np.arange(-3.0, 3.0, 10 / 3600)
reference = limit_lines(bes, th_check, niter=30)
rows = []
for n in (2, 4, 6, 8):
lim = limit_lines(bes, th_check, niter=n)
err = np.stack([np.linalg.norm(lim[s] - reference[s], axis=1) for s in ("north", "south")])
worst = np.nanargmax(np.fmax(err[0], err[1]))
rows.append({"iterations": n, "median error (m)": np.nanmedian(err), "largest error (m)": np.nanmax(err),
"largest error at (UTC)": bes.time(th_check[worst]).strftime("%H:%M:%S")})
pd.DataFrame(rows).style.format({"median error (m)": "{:.1e}", "largest error (m)": "{:.2g}"}).hide(axis="index")
| iterations | median error (m) | largest error (m) | largest error at (UTC) |
|---|---|---|---|
| 2 | 9.0e-01 | 1.3e+04 | 16:40:30 |
| 4 | 1.5e-05 | 3e+02 | 16:40:30 |
| 6 | 2.3e-10 | 6.8 | 16:40:30 |
| 8 | 0.0e+00 | 0.15 | 16:40:30 |
The largest errors all come from the first minutes of the path, at sunrise over the Pacific (the central line reaches the Earth at 16:40 UTC).
Duration of totality¶
For a fixed place, totality begins at second contact and ends at third contact: the two instants when the place is exactly on the edge of the umbra. They are the two roots in $t$ of
$$ u(t)^2 + v(t)^2 = L_2(t)^2, \qquad (u, v) = (x - \xi,\ y - \eta). $$
Around the moment of greatest eclipse the relative motion is very nearly uniform and straight, and the solution uses that.
- Closest approach, by Newton's method. The distance is smallest where $g(t) = ua + vb$ vanishes, since $g$ is half the time derivative of $u^2+v^2$. With straight-line motion $\dot g \approx a^2 + b^2 = n^2$, so Newton's step is $$ \Delta t = -\frac{ua + vb}{n^2}. $$ Every place starts from the same guess, $t_0$, even though totality on this map runs from 20 minutes before to nearly two hours after it. Because the motion is so nearly straight, $g$ is nearly linear in $t$ and Newton's method converges quickly from there: three steps are good to 0.1 s of duration and four to 0.1 ms. The code uses six.
- The miss distance. At closest approach the offset is perpendicular to the motion, and its length is the cross product $$ m = \frac{|ub - va|}{n}. $$
- The contacts. With the motion straight at speed $n$, the place is within $|L_2|$ of the axis for $\pm\sqrt{L_2^2 - m^2}/n$ about closest approach. The duration is $$ \tau = \frac{2\sqrt{L_2^2 - m^2}}{n}, $$ converted from hours to seconds. It is defined when $m < |L_2|$ (the place is inside the path) and $\zeta > 0$ (the Sun is above the horizon). Otherwise it is NaN.
This is totality_duration. Its argument is an array of ITRF positions of any shape, and each place gets its own Newton iterate. Because it takes ITRF positions, it also works for places above the ellipsoid: pass sk.itrfcoord(latitude_deg=..., longitude_deg=..., altitude=h).vector.
Duration contours. On the map below, geodetic_grid generates a 0.1° grid of 300,000 places. totality_duration evaluates all of them in one vectorised call, which takes about a tenth of a second, and matplotlib's contour draws the lines of equal duration at whole minutes.
The straight-line approximation in step 3 can be checked. Following Meeus (1989), the cell below refines each contact by re-evaluating the elements and the observer's position at the contact and solving the same quadratic again. For a few cities the refined duration differs from the one-step value by less than 0.1 s.
def totality_duration(bes, p_itrf, niter=6):
# Duration of totality (s) and hours from t0 of maximum eclipse at ITRF points (m)
p = np.asarray(p_itrf, dtype=float) / AE
th = np.zeros(p.shape[:-1])
with np.errstate(invalid="ignore", divide="ignore"):
return _duration_newton(bes, p, th, niter)
def _duration_newton(bes, p, th, niter):
for _ in range(niter):
e = bes(th)
xi, eta, zeta = (np.sum(p * v, -1) for v in plane_axes(e))
xi_dot, eta_dot = observer_rates(e, xi, eta, zeta)
u, v = e["x"] - xi, e["y"] - eta
a, b = e["x_dot"] - xi_dot, e["y_dot"] - eta_dot
n2 = a * a + b * b
th = th - (u * a + v * b) / n2
m = np.abs(u * b - v * a) / np.sqrt(n2)
L2 = np.abs(e["l2"] - zeta * bes.tan_f2)
tau = np.where((m < L2) & (zeta > 0), 2 * np.sqrt(np.abs(L2**2 - m**2) / n2) * 3600, np.nan)
return tau, th
def refined_contacts(bes, p_itrf, th_max, niter=4, cone=2):
# Contacts (hours from t0), each re-linearised at the contact itself: with the umbra
# (cone=2, second and third contacts) or the penumbra (cone=1, first and fourth)
l, tan_f = ("l1", bes.tan_f1) if cone == 1 else ("l2", bes.tan_f2)
p = np.asarray(p_itrf, dtype=float) / AE
contacts = []
for sign in (-1, 1):
th = th_max
for _ in range(niter):
e = bes(np.array([th]))
xi, eta, zeta = (np.sum(p * v, -1) for v in plane_axes(e))
xi_dot, eta_dot = observer_rates(e, xi, eta, zeta)
u, v = e["x"] - xi, e["y"] - eta
a, b = e["x_dot"] - xi_dot, e["y_dot"] - eta_dot
L = np.abs(e[l] - zeta * tan_f)
n2, g = a * a + b * b, u * a + v * b
th = th + ((-g + sign * np.sqrt(g * g - n2 * (u * u + v * v - L**2))) / n2)[0]
contacts.append(th)
return contacts
rows = []
for name, lat, lon in [("Mazatlán", 23.2494, -106.4111), ("Dallas", 32.7767, -96.7970),
("Cleveland", 41.4993, -81.6944), ("Burlington", 44.4759, -73.2121)]:
p = sk.itrfcoord(latitude_deg=lat, longitude_deg=lon).vector
tau, th_max = totality_duration(bes, p)
c2, c3 = refined_contacts(bes, p, float(th_max))
rows.append({"place": name, "second contact (UTC)": bes.time(c2).strftime("%H:%M:%S"),
"third contact (UTC)": bes.time(c3).strftime("%H:%M:%S"), "duration (s)": float(tau),
"refined (s)": (c3 - c2) * 3600, "difference (s)": float(tau) - (c3 - c2) * 3600})
pd.DataFrame(rows).style.format(precision=2).hide(axis="index")
| place | second contact (UTC) | third contact (UTC) | duration (s) | refined (s) | difference (s) |
|---|---|---|---|---|---|
| Mazatlán | 18:07:31 | 18:11:48 | 257.48 | 257.49 | -0.01 |
| Dallas | 18:40:43 | 18:44:34 | 231.49 | 231.45 | 0.04 |
| Cleveland | 19:13:45 | 19:17:34 | 228.83 | 228.86 | -0.03 |
| Burlington | 19:26:07 | 19:29:22 | 194.69 | 194.75 | -0.06 |
The 2024 path¶
Evaluate the elements every 2 seconds over the six-hour window; the shadow is on the Earth for about 3 h 17 min of it. Points off the Earth come back as NaN and are dropped. This handles the ends of the path, where only part of the umbra touches the Earth at sunrise and sunset. The code uses $|L_2|$, so for an annular eclipse ($L_2 > 0$, the vertex above the ground) the same functions give the path of annularity. The sign of $L_2$ tells which one it is, and it changes along the path of a hybrid eclipse. We also evaluate the duration of totality on a 0.1° grid for the contours.
def greatest_eclipse(bes):
# Minimum distance of the shadow axis from the Earth's centre (gamma) and its time
th = np.linspace(-3, 3, 21601)
e = bes(th)
gamma = np.hypot(e["x"], e["y"])
i = np.argmin(gamma)
return gamma[i], bes.time(th[i])
def compute_path(bes, step_seconds=2.0, grid=None):
th = np.arange(-3.0, 3.0, step_seconds / 3600)
path = {"bes": bes, "th": th, "central": central_line(bes, th)}
path.update(limit_lines(bes, th))
path["ll"] = {k: latlon(path[k]) for k in ("central", "north", "south")}
if grid is not None:
lon, lat = grid
tau, _ = totality_duration(bes, geodetic_grid(lon, lat))
path["grid"] = (lon, lat, tau)
return path
path2024 = compute_path(bes, grid=(np.arange(-125, -50, 0.1), np.arange(14, 54, 0.1)))
ok = np.isfinite(path2024["ll"]["central"][:, 0])
t_first, t_last = (bes.time(path2024["th"][i]) for i in np.flatnonzero(ok)[[0, -1]])
print(f"central line on the Earth from {t_first.strftime('%H:%M:%S')} to {t_last.strftime('%H:%M:%S')} UTC")
gamma, t_ge = greatest_eclipse(bes)
print(f"greatest eclipse at {t_ge.strftime('%H:%M:%S.%f')[:10]} UTC, gamma = {gamma:.4f} "
"(NASA: 18:17:18.3 UT with its Delta T, gamma = 0.3431)")
central line on the Earth from 16:40:00 to 19:54:28 UTC greatest eclipse at 18:17:19.8 UTC, gamma = 0.3431 (NASA: 18:17:18.3 UT with its Delta T, gamma = 0.3431)
Map¶
The map below shows the path of totality over North America:
- the band between the limits;
- the central line, with the time (UTC) every 10 minutes;
- the outline of the umbra every 20 minutes;
- contours of equal duration of totality, one colour per minute. The path is only about 40 pixels wide at this scale, too narrow for labels on the contours.
The cities are those of the table of local circumstances below, labelled with the duration of totality where they saw it. The inset shows the whole track.
PC = ccrs.PlateCarree()
LAMBERT = ccrs.LambertConformal(central_longitude=-96, central_latitude=39, standard_parallels=(33, 45))
LAND, OCEAN = "#f2efe6", "#dfe9f1"
BAND, CENTRE, OUTLINE = "#4a5878", "#CC3311", "#222222"
CONTOUR_COLORS = ["#c6dbef", "#9ecae1", "#6baed6", "#3182bd", "#08519c", "#08306b"] # 1..6 min
LABEL_BOX = dict(boxstyle="round,pad=0.12", fc="white", ec="none", alpha=0.75)
def fmt_duration(s):
return f"{int(s // 60)}m{int(round(s % 60)):02d}s"
def base_map(ax, extent):
# Lambert conformal map of North America with Natural Earth 1:50m features
ax.set_extent(extent, crs=PC)
ax.set_facecolor(OCEAN)
ax.add_feature(cfeature.LAND.with_scale("50m"), facecolor=LAND, edgecolor="none").set_rasterized(True)
ax.add_feature(cfeature.LAKES.with_scale("50m"), facecolor=OCEAN, edgecolor="#9aa8b3", linewidth=0.3)
ax.add_feature(cfeature.STATES.with_scale("50m"), edgecolor="#a39d8f", linewidth=0.35, facecolor="none")
ax.add_feature(cfeature.BORDERS.with_scale("50m"), edgecolor="#555555", linewidth=0.7)
ax.add_feature(cfeature.COASTLINE.with_scale("50m"), edgecolor="#6b7a86", linewidth=0.5)
gl = ax.gridlines(draw_labels=True, x_inline=False, y_inline=False, rotate_labels=False,
linewidth=0.3, color="#888888", alpha=0.5, linestyle="--")
gl.top_labels = gl.right_labels = False
gl.xlabel_style = gl.ylabel_style = {"size": 9}
def utc_marks(bes, th, minutes):
# (hours from t0, datetime) at each whole multiple of `minutes` UTC within th
first, last = bes.time(th[0]).to_datetime(), bes.time(th[-1]).to_datetime()
mark = first.replace(second=0, microsecond=0)
mark += datetime.timedelta(minutes=(-mark.minute) % minutes)
marks = []
while mark <= last:
marks.append((bes.hours(sk.time.from_datetime(mark)), mark))
mark += datetime.timedelta(minutes=minutes)
return marks
def draw_path(ax, path, cities=(), tick_minutes=10, outline_minutes=20, tick_side="south"):
bes, th = path["bes"], path["th"]
ll = {k: v[::5] for k, v in path["ll"].items()} # every 10 s is plenty for drawing
both = np.isfinite(ll["north"][:, 0]) & np.isfinite(ll["south"][:, 0])
band = np.concatenate([ll["north"][both], ll["south"][both][::-1]])
ax.fill(band[:, 1], band[:, 0], transform=PC, color=BAND, alpha=0.25, lw=0, label="path of totality")
for side in ("north", "south"):
ax.plot(ll[side][:, 1], ll[side][:, 0], transform=PC, color=BAND, lw=0.8)
ax.plot(ll["central"][:, 1], ll["central"][:, 0], transform=PC, color=CENTRE, lw=1.4,
label="central line")
# Contours of equal duration of totality, one colour per minute
lon, lat, tau = path["grid"]
levels = np.arange(1, int(np.nanmax(tau) // 60) + 1)
ax.contour(lon, lat, tau / 60, levels=levels, transform=PC, linewidths=0.7,
colors=CONTOUR_COLORS[: len(levels)])
for lev in levels:
ax.plot([], [], color=CONTOUR_COLORS[lev - 1], lw=1.2, label=f"{lev} min of totality")
# Ticks across the central line every `tick_minutes`, labelled beyond one limit
proj = ax.projection
x0, x1, y0, y1 = ax.get_extent(crs=proj)
for h, when in utc_marks(bes, th, tick_minutes):
lim = limit_lines(bes, np.array([h]))
pts = np.array([central_line(bes, np.array([h]))[0], lim["north"][0], lim["south"][0]])
if not np.all(np.isfinite(pts)):
continue
c, n, s = proj.transform_points(PC, *latlon(pts)[:, ::-1].T)[:, :2]
if not (x0 < c[0] < x1 and y0 < c[1] < y1):
continue
edge = s if tick_side == "south" else n
ax.plot(*np.array([c + 0.35 * (n - c), c + 0.35 * (s - c)]).T, transform=proj, color=CENTRE, lw=1.0)
ax.text(*(c + 1.55 * (edge - c)), when.strftime("%H:%M"), transform=proj, fontsize=8,
ha="center", va="center", color=CENTRE, bbox=LABEL_BOX)
if outline_minutes and when.minute % outline_minutes == 0:
o = latlon(umbra_outline(bes, h))
ax.plot(o[:, 1], o[:, 0], transform=PC, color=OUTLINE, lw=0.6, ls=(0, (3, 2)))
if outline_minutes:
ax.plot([], [], color=OUTLINE, lw=0.6, ls=(0, (3, 2)), label=f"umbra every {outline_minutes} min")
# Cities: (name, lat, lon, (label dx, dy in degrees, horizontal alignment))
for name, lat, lon, (dx, dy, ha) in cities:
tau_city, _ = totality_duration(bes, sk.itrfcoord(latitude_deg=lat, longitude_deg=lon).vector)
label = name if np.isnan(tau_city) else f"{name} {fmt_duration(tau_city)}"
ax.plot(lon, lat, "o", transform=PC, ms=3.5, color="k", mec="white", mew=0.6, zorder=5)
ax.text(lon + dx, lat + dy, label, transform=PC, fontsize=8, ha=ha, va="center",
bbox=LABEL_BOX, zorder=5)
def inset_globe(fig, path, rect, centre):
# Orthographic view of the whole track
ax = fig.add_axes(rect, projection=ccrs.Orthographic(*centre))
ax.set_global()
ax.set_facecolor(OCEAN)
ax.add_feature(cfeature.LAND, facecolor=LAND, edgecolor="none")
ax.add_feature(cfeature.COASTLINE, linewidth=0.3, edgecolor="#6b7a86")
ll = path["ll"]
both = np.isfinite(ll["north"][:, 0]) & np.isfinite(ll["south"][:, 0])
band = np.concatenate([ll["north"][both], ll["south"][both][::-1]])
ax.fill(band[:, 1], band[:, 0], transform=PC, color=BAND, alpha=0.6, lw=0)
ax.plot(ll["central"][:, 1], ll["central"][:, 0], transform=PC, color=CENTRE, lw=0.8)
ax.spines["geo"].set_linewidth(0.5)
return ax
# The cities of the local-circumstances table below, with label offsets
CITIES_2024 = [
("Mexico City", 19.4326, -99.1332, (0.5, 0, "left")),
("Austin", 30.2672, -97.7431, (0.5, -0.3, "left")),
("Dallas", 32.7767, -96.7970, (-0.5, 0.2, "right")),
("St. Louis", 38.6270, -90.1994, (-0.5, 0.2, "right")),
("Cleveland", 41.4993, -81.6944, (-0.5, 0.4, "right")),
("Burlington", 44.4759, -73.2121, (-0.5, -0.3, "right")),
("Montreal", 45.5017, -73.5673, (-0.5, 0.3, "right")),
("Quebec City", 46.8139, -71.2080, (0.4, 0.6, "left")),
("New York", 40.7128, -74.0060, (0.5, -0.2, "left")),
("Boston", 42.3601, -71.0589, (0.5, -0.2, "left")),
("Halifax", 44.6488, -63.5752, (0.5, -0.2, "left")),
]
fig, ax = plt.subplots(figsize=(10, 7.4), subplot_kw={"projection": LAMBERT})
base_map(ax, [-116, -62, 17, 50])
draw_path(ax, path2024, CITIES_2024)
ax.legend(loc="upper left", fontsize=8.5, framealpha=0.9)
ax.set_title("Total solar eclipse of 8 April 2024: path of totality")
inset_globe(fig, path2024, [0.67, 0.12, 0.21, 0.21], (-90, 25))
plt.show()
/opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/shapely/creation.py:225: RuntimeWarning: invalid value encountered in linestrings return lib.linestrings(coords, np.intc(handle_nan), out=out, **kwargs) /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/shapely/creation.py:225: RuntimeWarning: invalid value encountered in linestrings return lib.linestrings(coords, np.intc(handle_nan), out=out, **kwargs) /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/shapely/creation.py:225: RuntimeWarning: invalid value encountered in linestrings return lib.linestrings(coords, np.intc(handle_nan), out=out, **kwargs)
/opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/shapely/creation.py:225: RuntimeWarning: invalid value encountered in linestrings return lib.linestrings(coords, np.intc(handle_nan), out=out, **kwargs)
Local circumstances for cities¶
Local circumstances are what an observer at one place sees: when the partial eclipse begins and ends, when totality begins and ends if the place is inside the path, and how much of the Sun is covered at maximum. They follow from the geometry of the duration of totality, with the penumbra added.
- Contacts. First and fourth contact, the beginning and end of the partial eclipse, are the instants when the place is on the edge of the penumbra: $u^2 + v^2 = L_1^2$, with $L_1 = l_1 - \zeta\tan f_1$. Second and third contact use $L_2$, as before.
refined_contactssolves either one. The penumbral contacts are more than an hour from maximum, so the relative motion is no longer straight over that interval, and the re-linearisation at each contact matters: it converges in a few steps. - Maximum eclipse is the closest approach found by the Newton iteration in
totality_duration, at distance $m$ from the shadow axis. - Magnitude is the fraction of the Sun's diameter covered at maximum (Meeus 1989): $$ M = \frac{L_1 - m}{L_1 + L_2}, $$ with $L_2$ signed. $M = 0$ on the edge of the penumbra ($m = L_1$) and $M = 1$ on the edge of the umbra ($m = -L_2$); inside a total eclipse's umbra $L_2 < 0$ and $M > 1$. The ratio of the apparent diameters of the Moon and the Sun is $A = (L_1 - L_2)/(L_1 + L_2)$.
- Obscuration is the fraction of the Sun's area covered: the overlap of a disc of radius 1 (the Sun) and a disc of radius $A$ (the Moon) whose centres are $c = 1 + A - 2M$ apart, divided by $\pi$. It is 1 during totality.
The times are for sea level and ignore the lunar limb profile, like the rest of the tutorial. All these contacts happen with the Sun above the horizon ($\zeta > 0$); for places near sunrise or sunset, check $\zeta$ at each contact.
def obscuration(mag, A):
# Fraction of the Sun's disc (radius 1) covered by the Moon's (radius A), centres 1 + A - 2 mag apart
c = 1 + A - 2 * mag
if c <= A - 1:
return 1.0
if c >= 1 + A:
return 0.0
a1 = np.arccos((c**2 + 1 - A**2) / (2 * c))
a2 = np.arccos((c**2 + A**2 - 1) / (2 * c * A))
lens = a1 + A**2 * a2 - 0.5 * np.sqrt((1 + A - c) * (c + 1 - A) * (c - 1 + A) * (c + 1 + A))
return lens / np.pi
def local_circumstances(bes, p_itrf):
# Contacts (hours from t0), magnitude and obscuration at maximum, for one ITRF position (m)
_, th_max = totality_duration(bes, p_itrf)
th_max = float(th_max)
e = bes(np.array([th_max]))
xi, eta, zeta = (np.asarray(p_itrf) / AE @ v[0] for v in plane_axes(e))
m = np.hypot(e["x"][0] - xi, e["y"][0] - eta)
L1 = e["l1"][0] - zeta * bes.tan_f1
L2 = e["l2"][0] - zeta * bes.tan_f2 # signed: negative inside a total eclipse's umbra
mag = (L1 - m) / (L1 + L2)
A = (L1 - L2) / (L1 + L2)
c1, c4 = refined_contacts(bes, p_itrf, th_max, niter=6, cone=1)
c2, c3 = refined_contacts(bes, p_itrf, th_max) if m < -L2 else (np.nan, np.nan)
return {"c1": c1, "c2": c2, "max": th_max, "c3": c3, "c4": c4, "magnitude": mag,
"obscuration": obscuration(mag, A)}
def clock(th):
return None if np.isnan(th) else bes.time(th).strftime("%H:%M:%S")
rows = []
for name, lat, lon, _ in CITIES_2024:
lc = local_circumstances(bes, sk.itrfcoord(latitude_deg=lat, longitude_deg=lon).vector)
rows.append({"place": name, "partial begins": clock(lc["c1"]), "totality begins": clock(lc["c2"]),
"maximum": clock(lc["max"]), "totality ends": clock(lc["c3"]), "partial ends": clock(lc["c4"]),
"totality (s)": (lc["c3"] - lc["c2"]) * 3600, "magnitude": lc["magnitude"],
"obscuration": lc["obscuration"]})
pd.DataFrame(rows).style.format({"totality (s)": "{:.1f}", "magnitude": "{:.3f}", "obscuration": "{:.3f}"},
na_rep="–").hide(axis="index")
| place | partial begins | totality begins | maximum | totality ends | partial ends | totality (s) | magnitude | obscuration |
|---|---|---|---|---|---|---|---|---|
| Mexico City | 16:55:23 | – | 18:14:17 | – | 19:36:27 | – | 0.790 | 0.747 |
| Austin | 17:17:13 | 18:36:12 | 18:37:02 | 18:37:52 | 19:58:07 | 100.2 | 1.002 | 1.000 |
| Dallas | 17:23:18 | 18:40:43 | 18:42:38 | 18:44:34 | 20:02:41 | 231.4 | 1.015 | 1.000 |
| St. Louis | 17:42:58 | – | 19:00:49 | – | 20:17:22 | – | 0.988 | 0.991 |
| Cleveland | 17:59:22 | 19:13:45 | 19:15:40 | 19:17:34 | 20:28:59 | 228.9 | 1.022 | 1.000 |
| Burlington | 18:14:15 | 19:26:07 | 19:27:45 | 19:29:22 | 20:37:20 | 194.8 | 1.015 | 1.000 |
| Montreal | 18:14:28 | 19:26:50 | 19:27:33 | 19:28:16 | 20:36:52 | 85.8 | 1.002 | 1.000 |
| Quebec City | 18:18:21 | – | 19:30:17 | – | 20:38:23 | – | 0.983 | 0.987 |
| New York | 18:10:36 | – | 19:25:35 | – | 20:36:24 | – | 0.910 | 0.899 |
| Boston | 18:16:07 | – | 19:29:48 | – | 20:39:08 | – | 0.932 | 0.925 |
| Halifax | 18:27:20 | – | 19:38:04 | – | 20:44:08 | – | 0.946 | 0.942 |
Times are UTC. Five of the eleven cities saw totality. Austin, Montreal and Dallas lie inside the path, Austin and Montreal close to its southern and northern limits, so their totality is much shorter than Dallas's 3 min 51 s. St. Louis and Quebec City, just outside, saw 99% of the Sun's diameter covered. A direct calculation, which scans the topocentric angular separation of the Sun and the Moon every half second with their instantaneous apparent radii, reproduces every contact to within 1 s.
Comparison with NASA's path table¶
NASA's path table (Espenak, SE2024Apr08T path) lists the northern limit, southern limit and central line every two minutes, with the path width and the duration of totality on the central line. The eight rows below, between Mexico and Maine, are copied from it.
Each of NASA's points belongs to an instant in UT computed with $\Delta T$ = 70.6 s, so a timing difference moves a point along the path. What matters for an observer is the distance across it. For each of NASA's points we therefore give its distance from our corresponding line, which does not depend on timing.
Two checks separate the sources of error:
- Mapping only. Our mapping from elements to the ground, fed NASA's own elements and $\Delta T$, compared with NASA's points. This tests the formulas of the previous section and nothing else.
- satkit. The whole calculation, with the JPL ephemerides and the measured $\Delta T$.
# NASA / Espenak path table, 2024 Apr 08 ("Eclipse Predictions by Fred Espenak, NASA's GSFC")
# https://eclipse.gsfc.nasa.gov/SEpath/SEpath2001/SE2024Apr08Tpath.html
# UT, northern limit, southern limit, central line (deg, arcmin, N/S, deg, arcmin, E/W),
# path width (km), central duration (s)
NASA_PATH = [
("18:10", (23, 47.7, "N", 106, 50.7, "W"), (22, 36.2, "N", 105, 23.3, "W"), (23, 11.9, "N", 106, 6.8, "W"), 199, 267.2),
("18:30", (29, 32.2, "N", 101, 11.9, "W"), (28, 15.1, "N", 99, 49.8, "W"), (28, 53.6, "N", 100, 30.7, "W"), 195, 267.0),
("18:40", (32, 23.5, "N", 98, 1.3, "W"), (31, 3.1, "N", 96, 43.3, "W"), (31, 43.2, "N", 97, 22.2, "W"), 192, 263.7),
("18:50", (35, 14.1, "N", 94, 26.9, "W"), (33, 50.1, "N", 93, 14.6, "W"), (34, 32.0, "N", 93, 50.7, "W"), 189, 258.0),
("19:00", (38, 3.8, "N", 90, 18.7, "W"), (36, 36.0, "N", 89, 14.5, "W"), (37, 19.7, "N", 89, 46.6, "W"), 186, 250.0),
("19:10", (40, 51.6, "N", 85, 22.6, "W"), (39, 19.9, "N", 84, 29.8, "W"), (40, 5.6, "N", 84, 56.3, "W"), 183, 239.4),
("19:20", (43, 35.6, "N", 79, 16.9, "W"), (42, 0.4, "N", 78, 40.9, "W"), (42, 47.8, "N", 78, 59.2, "W"), 179, 225.9),
("19:30", (46, 11.6, "N", 71, 24.2, "W"), (44, 33.8, "N", 71, 13.8, "W"), (45, 22.4, "N", 71, 19.7, "W"), 173, 209.1),
]
def nasa_point(p):
lat = (p[0] + p[1] / 60) * (1 if p[2] == "N" else -1)
lon = (p[3] + p[4] / 60) * (1 if p[5] == "E" else -1)
return sk.itrfcoord(latitude_deg=lat, longitude_deg=lon)
def offset_km(point, line):
# Signed distance (km) from an itrfcoord to a polyline of ITRF points; + if the line is north
L = line[np.all(np.isfinite(line), axis=1)]
p = point.vector
i = np.argmin(np.linalg.norm(L - p, axis=1))
candidates = []
for j in (max(i - 1, 0), min(i, len(L) - 2)):
a, b = L[j], L[j + 1]
s = np.clip(np.dot(p - a, b - a) / np.dot(b - a, b - a), 0, 1)
candidates.append(a + s * (b - a))
q = min(candidates, key=lambda c: np.linalg.norm(c - p))
north = point.qenu2itrf * np.array([0.0, 1.0, 0.0])
return np.sign(np.dot(q - p, north)) * np.linalg.norm(q - p) / 1e3
# NASA's elements with NASA's Delta T, mapped with our formulas, as lines
nasa_path = compute_path(nasa_own_dt)
rows = []
for ut, pn, ps, pc, width, dur in NASA_PATH:
hh, mm = map(int, ut.split(":"))
t_ut = sk.time(2024, 4, 8, hh, mm, 0)
N, S, C = nasa_point(pn), nasa_point(ps), nasa_point(pc)
# our width, across the path through our central point at this UT
c = sk.itrfcoord(central_line(bes, np.array([bes.hours(t_ut)]))[0])
our_width = abs(offset_km(c, path2024["north"])) + abs(offset_km(c, path2024["south"]))
tau, _ = totality_duration(bes, C.vector[None, :])
rows.append({
"UT": ut,
"mapping only: max |offset| (km)": max(abs(offset_km(P, nasa_path[k]))
for P, k in ((N, "north"), (C, "central"), (S, "south"))),
"north limit (km)": offset_km(N, path2024["north"]),
"central line (km)": offset_km(C, path2024["central"]),
"south limit (km)": offset_km(S, path2024["south"]),
"width NASA (km)": width,
"width satkit (km)": our_width,
"duration NASA (s)": dur,
"duration satkit (s)": tau[0],
})
pd.DataFrame(rows).style.format(precision=2).format(precision=1, subset=["width satkit (km)", "duration NASA (s)", "duration satkit (s)"]).hide(axis="index")
| UT | mapping only: max |offset| (km) | north limit (km) | central line (km) | south limit (km) | width NASA (km) | width satkit (km) | duration NASA (s) | duration satkit (s) |
|---|---|---|---|---|---|---|---|---|
| 18:10 | 0.08 | 0.86 | 1.03 | 1.04 | 199 | 198.9 | 267.2 | 267.0 |
| 18:30 | 0.06 | 0.88 | 0.93 | 1.01 | 195 | 194.6 | 267.0 | 266.8 |
| 18:40 | 0.08 | 0.78 | 0.84 | 0.92 | 192 | 192.1 | 263.7 | 263.5 |
| 18:50 | 0.06 | 0.77 | 0.79 | 0.93 | 189 | 189.3 | 258.0 | 257.9 |
| 19:00 | 0.06 | 0.69 | 0.79 | 0.77 | 186 | 186.2 | 250.0 | 249.8 |
| 19:10 | 0.09 | 0.66 | 0.72 | 0.86 | 183 | 182.6 | 239.4 | 239.2 |
| 19:20 | 0.05 | 0.59 | 0.64 | 0.68 | 179 | 178.3 | 225.9 | 225.8 |
| 19:30 | 0.10 | 0.48 | 0.70 | 0.62 | 173 | 173.2 | 209.1 | 209.0 |
Reading the table. Each row is one of NASA's two-minute rows, from Mexico (18:10 UT) to Maine (19:30 UT).
- Mapping only: max |offset| is the largest distance of NASA's three points (north, centre, south) from the lines our code draws from NASA's own elements and $\Delta T$. NASA's coordinates are rounded to 0.1′, up to 185 m, so values below about 0.1 km mean that the two mappings agree to the rounding of the table.
- North limit, central line, south limit are the signed distances from each of NASA's points to our full calculation's line, measured across the path. They are positive when our line lies north of NASA's point.
- Width is NASA's tabulated width and ours, measured across the path through our central point at the same UT.
- Duration is the duration of totality at NASA's central-line point: NASA's value, and ours from
totality_duration.
What the numbers say. With NASA's elements, our mapping reproduces NASA's points to the rounding of the table, so the formulas agree. The full satkit calculation puts all three lines about 0.5 to 1 km north (more precisely, north-west, across the path) of NASA's. The width agrees to within a kilometre (NASA rounds it to whole kilometres), and the central duration to 0.2 s. The shift has two parts of similar size:
- $\Delta T$. The Earth turned 1.4 s further than NASA assumed. That moves the shadow track 21″ of longitude, about 0.5 km, west. Across a path heading north-east, most of that appears as a shift to the north-west.
- Ephemeris. The shadow axis ($x$) is 0.6 km further west on the fundamental plane with DE440 than with ELP2000-85, which also moves the path north-west.
Our central durations are 0.1 to 0.2 s shorter than NASA's. On the central line $\tau = 2|L_2|/n$, so the duration changes by about 3 s per kilometre of umbral radius. Our $|l_2|$ is 30 m smaller than NASA's, which alone accounts for about 0.1 s; NASA's durations are rounded to 0.1 s.
Both differences are smaller than the effects neither calculation includes. The real lunar limb is not a circle: its mountains and valleys move the limits by 1 to 3 km (NASA's page gives this range). The centre of the Moon's figure is also offset from its centre of mass by about 0.5 km. Observers who want the limits to better than a kilometre use limb-profile corrections from lunar laser altimetry. For choosing where to stand, both calculations give the same path.
Cross-check without the fundamental plane¶
The Besselian mapping can be checked with a direct 3-D calculation that uses neither the fundamental plane nor the polynomial fit. Intersect the Sun–Moon line with the WGS84 ellipsoid, using the same apparent Sun and light-time-corrected Moon.
def direct_central_line(times):
sun, moon = sun_moon_itrf(times)
u = (moon - sun) / np.linalg.norm(moon - sun, axis=1)[:, None] # direction of the light
a, b = moon * STRETCH, u * STRETCH
qa, qb, qc = np.sum(b * b, 1), np.sum(a * b, 1), np.sum(a * a, 1) - 1
s = (-qb - np.sqrt(qb**2 - qa * qc)) / qa # first intersection along the light
return (moon + s[:, None] * u) * AE
times = [sk.time(2024, 4, 8, 17, 0, 0) + sk.duration.from_minutes(m) for m in range(0, 150, 5)]
direct = direct_central_line(times)
besselian = central_line(bes, np.array([bes.hours(t) for t in times]))
print(f"max difference, direct vs Besselian central line: {np.max(np.linalg.norm(direct - besselian, axis=1)):.1f} m")
max difference, direct vs Besselian central line: 0.2 m
The two agree to metres: the polynomial fit and the mapping lose nothing at the scale of the map.
Another date: 12 August 2045¶
The next total eclipses visible from the contiguous United States are on 23 August 2044 and 12 August 2045. The 2044 eclipse crosses Canada and ends at sunset in Montana and North Dakota. The 2045 eclipse crosses the country from northern California to Florida, with up to 6 minutes of totality.
The code above only needs a reference epoch $t_0$ near greatest eclipse. eclipse_t0 finds one from the elements themselves. It scans the day for the minimum distance of the shadow axis from the Earth's centre, $\sqrt{x^2+y^2}$ (the eclipse's $\gamma$), and rounds to the nearest hour of TT.
Earth orientation for 2045 is not known yet. satkit holds the last UT1−UTC value of its table (and warns; the warning is turned off below), and future leap seconds are unknown too. See the limitations below for what that means for the path.
def eclipse_t0(year, month, day):
# Whole hour of TT nearest to the minimum of sqrt(x^2 + y^2) on the given date
start = sk.time(year, month, day)
times = [start + sk.duration.from_minutes(m) for m in range(0, 1440, 10)]
e = raw_elements(times)
t_min = times[int(np.argmin(np.hypot(e["x"], e["y"])))]
jd = t_min.to_jd(sk.timescale.TT)
return sk.time.from_jd(np.round((jd - 0.5) * 24) / 24 + 0.5, sk.timescale.TT)
sk.frametransform.disable_eop_time_warning()
bes2045 = Besselian.fit(eclipse_t0(2045, 8, 12))
gamma, t_ge = greatest_eclipse(bes2045)
print(f"t0 = {bes2045.t0} UTC; greatest eclipse at {t_ge.strftime('%H:%M:%S')} UTC, gamma = {gamma:.4f}")
path2045 = compute_path(bes2045, grid=(np.arange(-136, -62, 0.1), np.arange(17, 52, 0.1)))
print(f"longest totality on the map: {fmt_duration(np.nanmax(path2045['grid'][2]))}")
CITIES_2045 = [
("Reno", 39.5296, -119.8138, (0.4, 0.8, "left")),
("Salt Lake City", 40.7608, -111.8910, (0.4, 0.8, "left")),
("Denver", 39.7392, -104.9903, (0.5, 0.6, "left")),
("Tulsa", 36.1540, -95.9928, (0.5, 0.6, "left")),
("Little Rock", 34.7465, -92.2896, (0.5, 0.6, "left")),
("Atlanta", 33.7490, -84.3880, (0.5, 0.3, "left")),
("Orlando", 28.5384, -81.3789, (0.5, 0.5, "left")),
("Miami", 25.7617, -80.1918, (0.5, 0.2, "left")),
]
fig, ax = plt.subplots(figsize=(10, 6.8), subplot_kw={"projection": LAMBERT})
base_map(ax, [-125, -70, 21, 49])
draw_path(ax, path2045, CITIES_2045, outline_minutes=0)
ax.legend(loc="lower left", fontsize=8.5, framealpha=0.9)
ax.set_title("Total solar eclipse of 12 August 2045: path of totality")
plt.show()
t0 = 2045-08-12T17:58:50.816000Z UTC; greatest eclipse at 17:41:29 UTC, gamma = 0.2116
longest totality on the map: 6m06s
/opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/shapely/creation.py:225: RuntimeWarning: invalid value encountered in linestrings return lib.linestrings(coords, np.intc(handle_nan), out=out, **kwargs) /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/shapely/creation.py:225: RuntimeWarning: invalid value encountered in linestrings return lib.linestrings(coords, np.intc(handle_nan), out=out, **kwargs) /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/shapely/creation.py:225: RuntimeWarning: invalid value encountered in linestrings return lib.linestrings(coords, np.intc(handle_nan), out=out, **kwargs)
Limitations¶
The calculation is the classical one, and it shares the classical simplifications.
- A smooth lunar limb. The Moon is modelled as a sphere of radius $k_2$. Its real limb has mountains and valleys a few kilometres high, which move the true limits by 1 to 3 km and change the duration by a second or so on the central line, and by much more near the limits. At the start and end of totality, sunlight shines through the valleys as Baily's beads. Precise limit predictions use limb profiles from lunar laser altimetry, which this tutorial does not attempt. The single radius $k_2$, fitted to the valleys, is a good average.
- A sea-level ellipsoid. The limits, the central line and the duration map are for places on the WGS84 ellipsoid, $h = 0$. For places at height $h$ the path shifts by roughly $h/\tan(\text{altitude})$ towards the Sun's azimuth. That is 0.6 km for $h$ = 1 km with the Sun at 60°, and more when the Sun is low.
totality_durationaccepts any ITRF position, so local circumstances at altitude are one argument away. The drawn lines, however, assume $h = 0$. The geoid, which is up to about 100 m from the ellipsoid, is ignored as well. - $\Delta T$ for the future. For 2024 we use the measured UT1. For 2045
satkitholds the last known UT1−UTC, and no future leap seconds are known, so it effectively assumes $\Delta T$ stays near today's 69 s. The real value is uncertain by several seconds. Each second of $\Delta T$ error shifts the whole path $15''$ in longitude, about 0.35 km at 40° latitude, and shifts every contact time by about a second. The path's shape, width and durations are not affected. - Refraction. The geometry ignores the atmosphere. Near sunrise and sunset refraction lifts the Sun by about half a degree, which slightly moves the ends of the path and the limits where the Sun is low.
- Cubic fits over ±3 h. The polynomials are valid only inside the six-hour window they were fitted to. That window covers the whole of these two eclipses, but a longer event would need a longer window or a higher degree.
References¶
- Explanatory Supplement to the Astronomical Almanac, 3rd ed. (Urban & Seidelmann 2013): the eclipse chapter defines the Besselian elements, the fundamental plane and the formulas for the central line, limits and local circumstances used here.
- Meeus, Elements of Solar Eclipses 1951–2200 (Meeus 1989): the same method in a compact form, with the contact-time iteration used for the refined contacts.
- Espenak, NASA eclipse web site (Espenak): the Besselian elements and path table used for validation, and the lunar radii $k_1$, $k_2$ and solar semi-diameter.
- JPL DE440 (Park et al. 2021): the Sun and Moon positions.
- VSOP87 (Bretagnon & Francou 1988) and ELP2000-85 (Chapront-Touzé & Chapront 1988): the analytical theories behind NASA's elements.
- IERS Conventions (2010) (Petit & Luzum 2010): the GCRF→ITRF rotation, UT1 and polar motion.
- WGS84 (NGA 2014): the Earth ellipsoid, $a_e$ and $f$.
- Bowring (1976) (Bowring 1976): the geodetic-latitude conversion in
satkit.itrfcoord.