Eclipse seasons

In the age of computing on paper especially, other conditions were used as eclipse criteria: faster ones, and sometimes less exact.

In the situation of figure 1.5, the criteria could be set at the moment of new moon both for the Moon’s height in the vertical direction (the ecliptic latitude) ${\beta }_m$ and for the Sun’s distance from the node in the horizontal direction (the segment SN).

Because the inclination of the orbit $I$ and the semidiameters $s_s$ and $s_m$ vary a little over time, limits of that kind are not quite exact, unlike the $\sigma$ limits above.

These quantities can still be given limits that say when an eclipse is certain and when it is possible. The older edition of the Explanatory Supplement (1961) gives the limits in the vertical direction: a solar eclipse is possible when ${\beta }_m<1.58^{\circ }$ and certain when ${\beta }_m<1.41^{\circ }$.

Espenak gives limits for the Sun’s horizontal distance for solar eclipses ($18.59^{\circ }$ and $15.39^{\circ }$) and for lunar eclipses ($16.8^{\circ }$ and $15.2^{\circ }$).

These horizontal distances are convenient if the eclipse seasons introduced below are to be marked on a calendar.

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(2039, 1, 1), ts.utc(2040, 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, ph in zip(t, phase):
    if ph not in (0, 2):                        # new and full moons only
        continue

    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_m = degrees(R_EARTH / d_moon)
    pi_s = degrees(R_EARTH / d_sun)

    if ph == 0:                                 # new moon: chapter 2's limits
        partial = s_s + s_m + pi_m - pi_s
        central = abs(s_s - s_m + pi_m - pi_s)
        if sigma > partial:
            verdict = '-'
        elif sigma > central:
            verdict = 'partial solar'
        else:
            verdict = 'total solar' if s_m > s_s else 'annular solar'
    else:                                       # full moon: chapter 3's limits
        pi_1 = 0.998340 * pi_m
        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:
            verdict = 'total lunar'
        elif sigma < f_2 + s_m:
            verdict = 'partial lunar'
        elif sigma < f_1 + s_m:
            verdict = 'penumbral lunar'
        else:
            verdict = '-'

    print(f"{ti.utc_strftime('%Y-%m-%d')}  {verdict}")

from numpy import cross
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 sun_south_of_orbit(t):
    """Is the Sun on the south side of the plane of the Moon's orbit?"""
    r, v = (moon - earth).at(t).frame_xyz_and_velocity(ecliptic_frame)
    h = cross(r.km, v.km_per_s, axis=0)   # normal to the plane of the orbit
    s = earth.at(t).observe(sun).apparent().frame_xyz(ecliptic_frame).km
    return (s * h).sum(axis=0) < 0

sun_south_of_orbit.step_days = 5     # the Sun needs half a year between nodes

t, south = almanac.find_discrete(ts.utc(2039, 1, 1), ts.utc(2040, 1, 1),
                                 sun_south_of_orbit)

for ti, s in zip(t, south):
    # the Sun passes to the south side at the ascending node
    print(f"{ti.utc_strftime('%Y-%m-%d %H:%M')}  "
          f"{'ascending' if s else 'descending'}")