Chapter 3

Finding lunar eclipses

Next, we examine lunar eclipses, which happen in the same way as solar eclipses, but at full moon. A full moon comes every month, but again, not every one of them is eclipsed.

Let us look at where the full moon travels in the sky in a bad month and in a good month.

In figure 1.1 the Sun moves against the stars, but of course neither the stars nor the very thin Moon can be seen in a daytime sky. Imagination is required.

The figure 3.1 above shows the night sky, so the motion of the stars and of the Moon can be observed. The Earth’s shadow, on the other hand, cannot be seen in the night sky either; it appears only when it causes a lunar eclipse.

3.1 Positions and distances

The findings of chapter 1 work for lunar eclipses as such. We’d only need to replace $S$, the direction of the Sun, with the direction of the antisolar point, that is of the Earth’s shadow (the opposite point on the celestial sphere). But it turns out this is not even needed.

It’s enough to just find all times of full moon and the moon’s ecliptic latitude ${\beta }_m$ at those times, and the formula for $\sigma$ will work as is. As an example, let us compute the least angular separation for the cases in figure 3.1.

from numpy import arctan, cos
from skyfield import almanac
from skyfield.api import load
from skyfield.framelib import ecliptic_frame

ts = load.timescale()
eph = load('de440s.bsp')
earth, sun, moon = eph['earth'], eph['sun'], eph['moon']

def lonlat(body, ti):
    lat, lon, _ = earth.at(ti).observe(body).apparent().frame_latlon(ecliptic_frame)
    return lon.degrees, lat.degrees

def wrap(d):
    return (d + 180) % 360 - 180

for y, m, d in [(2028, 10, 3), (2028, 12, 31)]:
    t, phase = almanac.find_discrete(ts.utc(y, m, d - 1), ts.utc(y, m, d + 1),
                                     almanac.moon_phases(eph))
    ti = t[phase == 2][0]                       # the full moon itself

    a, b = ts.tt_jd(ti.tt - 1 / 1440), ts.tt_jd(ti.tt + 1 / 1440)
    ml_a, mb_a = lonlat(moon, a)
    ml_b, mb_b = lonlat(moon, b)
    sl_a, _ = lonlat(sun, a)
    sl_b, _ = lonlat(sun, b)

    beta = lonlat(moon, ti)[1]
    lam = wrap(ml_b - ml_a) / wrap(sl_b - sl_a)
    tan_I = (mb_b - mb_a) / wrap(ml_b - ml_a)
    I_ = arctan(lam / (lam - 1) * tan_I)
    sigma = abs(beta) * cos(I_)

    print(f"{ti.utc_strftime('%Y-%m-%d %H:%M')} - {sigma:.3f}")

2028-10-03 16:25 - 4.939
2028-12-31 16:49 - 0.315

It seems there is a lunar eclipse in December. But the limits computed in chapter 2 do not apply. What is the angular limit for a lunar eclipse to be seen?

3.2 The limits for a lunar eclipse

The mathematics of the next exercise in angles is much like that for solar eclipses. There are new simplifications, but also complications. The calculation is made easier by the observer’s position having no bearing on whether an eclipse happens: the Moon is either eclipsed or it is not.

At least two things make it harder. The Earth is not a sphere, and it is surrounded by an atmosphere.

Let us start with the Earth’s penumbra. For the original derivation, see Explanatory Supplement section 11.2.3.

In figure 3.2 the point $M$ is out in space, on the edge of the penumbra at the Moon’s distance. The quantity we’re after is $f_1$, the semidiameter of the penumbra as seen from the Earth.

A little exercise with exterior angles gives

where ${\pi }_m$ and ${\pi }_s$ are the parallaxes of the Moon and the Sun and $s_s$ is the Sun’s apparent radius.The careful reader will notice that the parallaxes do not quite match the definitions of the previous chapter. ${\pi }_s$ is at the Sun’s limb rather than at its centre. The resulting error is 0.00001 arcseconds. The point $M$ is also in reality slightly further away than the Moon’s distance. The true distance of the Moon would meet the plane through $M$ at the middle of the shadow’s axis. The resulting error is about 0.8 arcseconds. We shall soon see that errors of this size are of very little consequence in lunar eclipses.

Figure 3.3 is the same situation for the umbra, whose semidiameter comes out as

Now the first nuisance. The Earth is not a sphere but a flattened ellipsoid. Computing the shadow of an ellipsoid in every orientation would be a dreadful job. The situation is remedied by using, instead of the Earth’s equatorial radius, the mean radius at latitude 45 degrees. This is done by using, in place of the parallax ${\pi }_m$,

The second nuisance. It was noticed long ago that these geometrical calculations depart somewhat from the observed times of eclipses. The Earth’s shadow does appear to be larger than the calculations would suggest. The cause is the Earth’s atmosphere. How much larger? One percent, five percent? Let us say two percent. Hence

No, really. These factors are not computed but determined by observation. Philippe de La Hire noticed as early as 1707 that the computed radius of the shadow has to be enlarged by about 1/41 for it to fit the observed contact times. Later 1/40 (Lambert) and 1/60 (Mayer) were proposed, and Beer and Mädler arrived at 1/50 from the eclipse of 1833. Chauvenet settled on it in 1891, and it became the practice of the ephemeris publications. Danjon pointed out in 1951 that the relative correction is being made in the wrong place. If the atmosphere acts like an opaque shell, it should enlarge the Earth’s radius and not the shadows. He estimated the thickness of that layer at about 75 km, or 1/85 of the Earth’s radius, which leaves the shadows slightly smaller. The difference from uncorrected predictions is observable, but it is hard to see any difference between the corrections themselves, especially as the edge of the umbra is in reality soft and ill-defined. The outer edge of the penumbra is so faint that the corrections have no practical effect at all. In long eclipse tables spanning centuries, a few eclipses may turn from total into partial. The Astronomical Almanac and the Explanatory Supplement (and this site) keep to 1/50, while Espenak's canons and the Connaissance des Temps use Danjon’s method. For a thorough treatment, see Espenak: Enlargement of Earth's Shadows.

The geocentric conditions for a lunar eclipse are now easy to state.

A penumbral eclipse happens when the angular separation of the centres of the Moon and of the shadow is $\sigma <f_1+s_m$; for a partial one the limit is $f_2+s_m$, and for a total one $f_2-s_m$.

In summary, the three conditions are

3.3 A worked example

Let us now find the lunar eclipses of 2028.

from numpy import arctan, cos, degrees
from skyfield import almanac
from skyfield.api import load
from skyfield.framelib import ecliptic_frame

R_EARTH = 6378.1366           # km, IAU equatorial radius
R_SUN = 696000.0
R_MOON = 0.2725076 * R_EARTH  # Moon/Earth radius ratio, IAU

ts = load.timescale()
eph = load('de440s.bsp')
earth, sun, moon = eph['earth'], eph['sun'], eph['moon']

t0, t1 = ts.utc(2028, 1, 1), ts.utc(2029, 1, 1)
t, phase = almanac.find_discrete(t0, t1, almanac.moon_phases(eph))

def lonlat(body, ti):
    lat, lon, _ = earth.at(ti).observe(body).apparent().frame_latlon(ecliptic_frame)
    return lon.degrees, lat.degrees

def wrap(d):
    return (d + 180) % 360 - 180

for ti in t[phase == 2]:                        # every full moon of the year
    a, b = ts.tt_jd(ti.tt - 1 / 1440), ts.tt_jd(ti.tt + 1 / 1440)
    ml_a, mb_a = lonlat(moon, a)
    ml_b, mb_b = lonlat(moon, b)
    sl_a, _ = lonlat(sun, a)
    sl_b, _ = lonlat(sun, b)

    beta = lonlat(moon, ti)[1]
    lam = wrap(ml_b - ml_a) / wrap(sl_b - sl_a)
    tan_I = (mb_b - mb_a) / wrap(ml_b - ml_a)
    I_ = arctan(lam / (lam - 1) * tan_I)
    sigma = abs(beta) * cos(I_)

    e = earth.at(ti)
    d_sun = e.observe(sun).apparent().distance().km
    d_moon = e.observe(moon).apparent().distance().km

    s_s = degrees(R_SUN / d_sun)
    s_m = degrees(R_MOON / d_moon)
    pi_s = degrees(R_EARTH / d_sun)
    pi_1 = 0.998340 * degrees(R_EARTH / d_moon)   # the 45-degree radius

    f_1 = 1.02 * (pi_1 + pi_s + s_s)
    f_2 = 1.02 * (pi_1 + pi_s - s_s)

    if sigma < f_2 - s_m:                       # tightest limit first, or a
        verdict = 'total'                       # total would report as partial
    elif sigma < f_2 + s_m:
        verdict = 'partial'
    elif sigma < f_1 + s_m:
        verdict = 'penumbral'
    else:
        verdict = '-'

    print(f"{ti.utc_strftime('%Y-%m-%d %H:%M')}  {sigma:6.3f}"
          f"  {f_1 + s_m:6.3f}  {f_2 + s_m:6.3f}  {f_2 - s_m:6.3f}  {verdict}")

Let us repeat the calculation for 2027 as well, to find penumbral eclipses.

At July’s eclipse the angular separation is 1.418° and the penumbral limit is 1.432°. The difference is 0.014°, that is 50 arcseconds. An eclipse even less thrilling than this one will probably be a long time coming. Mark it in your calendar!

3.4 A special penumbral eclipse

What if the case $M_P$ of figure 3.4 served as a limit of another kind: how often is the Moon completely within the penumbra but not at all within the umbra? This can happen when

Let us find the next few total penumbral eclipses.

from numpy import arctan, cos, degrees
from skyfield import almanac
from skyfield.api import load
from skyfield.framelib import ecliptic_frame

R_EARTH = 6378.1366           # km, IAU equatorial radius
R_SUN = 696000.0
R_MOON = 0.2725076 * R_EARTH  # Moon/Earth radius ratio, IAU

ts = load.timescale()
eph = load('de440s.bsp')
earth, sun, moon = eph['earth'], eph['sun'], eph['moon']

def lonlat(body, ti):
    lat, lon, _ = earth.at(ti).observe(body).apparent().frame_latlon(ecliptic_frame)
    return lon.degrees, lat.degrees

def wrap(d):
    return (d + 180) % 360 - 180

found = 0
year = 2027

while found < 3:                                # one year at a time, until three
    t, phase = almanac.find_discrete(ts.utc(year, 1, 1), ts.utc(year + 1, 1, 1),
                                     almanac.moon_phases(eph))
    for ti in t[phase == 2]:
        a, b = ts.tt_jd(ti.tt - 1 / 1440), ts.tt_jd(ti.tt + 1 / 1440)
        ml_a, mb_a = lonlat(moon, a)
        ml_b, mb_b = lonlat(moon, b)
        sl_a, _ = lonlat(sun, a)
        sl_b, _ = lonlat(sun, b)

        beta = lonlat(moon, ti)[1]
        lam = wrap(ml_b - ml_a) / wrap(sl_b - sl_a)
        tan_I = (mb_b - mb_a) / wrap(ml_b - ml_a)
        I_ = arctan(lam / (lam - 1) * tan_I)
        sigma = abs(beta) * cos(I_)

        e = earth.at(ti)
        d_sun = e.observe(sun).apparent().distance().km
        d_moon = e.observe(moon).apparent().distance().km

        s_s = degrees(R_SUN / d_sun)
        s_m = degrees(R_MOON / d_moon)
        pi_s = degrees(R_EARTH / d_sun)
        pi_1 = 0.998340 * degrees(R_EARTH / d_moon)

        f_1 = 1.02 * (pi_1 + pi_s + s_s)
        f_2 = 1.02 * (pi_1 + pi_s - s_s)

        if f_2 + s_m < sigma < f_1 - s_m:       # inside the penumbra, clear of the umbra
            print(ti.utc_strftime('%Y-%m-%d %H:%M'))
            found += 1
            if found == 3:
                break
    year += 1

2053-08-29 07:53
2066-12-31 14:42
2070-04-25 09:32

These take some waiting for. Better mark these in the calendar too, so they don’t slip past.

For more, see Meeus (1980). Note that the list in that source differs and that 2066 is missing. The difference is explained by the differences in the chosen percentages for enlarging the shadows, discussed above.