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.

Acquire

Would the first problem — which would require you to know the speed at which the Moon travels across the sky, which varies — go away if you asked for the Moon’s position at “t minus one day” and checked whether that position was on the other side of 180°? – Spree

@BrandonRhodes Thanks, that works perfectly. Can a similar logic be applied for eclipses as well? Is the angle of the two planes behaving monotonically? If not, what's a sensible way to determine eclipses? – Acquire

— I have never written any eclipse-finding routines, alas (which is why I added merely a comment and didn’t attempt an answer to your question). I know that it involves the 3D shape of the Earth’s shadow in space, but that since the shadow can be approximated as a pair of cones (umbral and penumbral) from the Moon’s point of view, diagrams always seem to approximate the shadow as a pair of circles that the Moon passes through — a plane section through those two cones. So I guess they measure the Moon’s position relative to those circles of shadow at its distance? – Spree

@BrandonRhodes Thanks for the general idea. This tells me I am not experienced enough to code it myself. Therefore, I'm simply going to refactor this question and bump up a little bounty. – Acquire

Since the days are counting down, I’ll ask here in a comment whether my answer gets close to meeting the terms of your bounty, or if it leaves any crucial questions still open. Thanks! – Spree

@BrandonRhodes Yes, so far it looks like exactly what I've wanted. I'll accept the answer before bounty's expiry as I'd like to make one extensive test before that. – Acquire

I will be interested to hear how the test comes out! – Spree

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

Spree