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
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.
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: