Determining lunar eclipse in skyfield

Question

I am given a list of dates in UTC, all hours cast to 00:00.

I'd like to determine if a (lunar) eclipse occurred in a given day (ie past 24 hours)

Considering the python snippet

from sykfield.api import load
eph = load('de421.bsp')
def eclipticangle(t):

    moon, earth = eph['moon'], eph['earth']
    e = earth.at(t)
    x, y, _ = e.observe(moon).apparent().ecliptic_latlon()

    return x.degrees

I am assuming one is able to determine if an eclipse occurred within 24hrs of a time t by

1. Checking that the first angle is close enough to 180 (easy)
2. Checking if the second degree is close enough to 0 (not so essy?)

Now as far as the answer in the comment suggests it is not so trivial to solve the second problem simply by testing if the angle is close to 0.

Therefore, my question is

Can someone provide a function to determine if a lunar eclipse occurred on a given day t?

Edit. This question was edited to reflect the feedback from Brandon Rhodes left in the comments below.

Answer 1

I just went through section 11.2.3 of the Explanatory Supplement to the Astronomical Almanac and tried turning it into Skyfield Python code. Here is what I came up with:

import numpy as np

from skyfield.api import load
from skyfield.constants import ERAD
from skyfield.functions import angle_between, length_of
from skyfield.searchlib import find_maxima

eph = load('de421.bsp')
earth = eph['earth']
moon = eph['moon']
sun = eph['sun']

def f(t):
    e = earth.at(t).position.au
    s = sun.at(t).position.au
    m = moon.at(t).position.au
    return angle_between(s - e, m - e)

f.step_days = 5.0

ts = load.timescale()
start_time = ts.utc(2019, 1, 1)
end_time = ts.utc(2020, 1, 1)

t, y = find_maxima(start_time, end_time, f)

e = earth.at(t).position.m
m = moon.at(t).position.m
s = sun.at(t).position.m

solar_radius_m = 696340e3
moon_radius_m = 1.7371e6

pi_m = np.arcsin(ERAD / length_of(m - e))
pi_s = np.arcsin(ERAD / length_of(s - e))
s_s = np.arcsin(solar_radius_m / length_of(s - e))

pi_1 = 0.998340 * pi_m

sigma = angle_between(s - e, e - m)
s_m = np.arcsin(moon_radius_m / length_of(e - m))

penumbral = sigma < 1.02 * (pi_1 + pi_s + s_s) + s_m
partial = sigma < 1.02 * (pi_1 + pi_s - s_s) + s_m
total = sigma < 1.02 * (pi_1 + pi_s - s_s) - s_m

mask = penumbral | partial | total

t = t[mask]
penumbral = penumbral[mask]
partial = partial[mask]
total = total[mask]

print(t.utc_strftime())
print(0 + penumbral + partial + total)

It produces a vector of times at which lunar eclipses occurred, and then a rating of how total the eclipse is:

['2019-01-21 05:12:51 UTC', '2019-07-16 21:31:27 UTC']
[3 2]

Its eclipse times are within 3 seconds of the times given in the huge table of lunar ephemerides at NASA:

https://eclipse.gsfc.nasa.gov/5MCLE/5MKLEcatalog.txt