确定天空中的月食

如何解决确定天空中的月食

我获得了以UTC列出的日期列表，所有时间都强制转换为00:00。

我想确定在特定的一天（即过去的24小时）内是否发生了（月食）

考虑python代码段

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

我假设人们能够确定在 t 时间之内24小时之内是否发生了日食

1. 检查第一个角度是否足够接近180度（容易）
2. 检查第二学位是否足够接近0（不是很麻烦吗？）

现在，据评论中的答案表明，仅通过测试角度是否接近0来解决第二个问题并不是那么简单。

因此，我的问题是

有人可以提供确定在给定的t天是否发生月食的功能吗？

编辑。此问题经过编辑，以反映以下评论中剩余的Brandon Rhodes的反馈。

解决方法

我刚刚阅读了《天文历书解释性补充说明》第11.2.3节，并尝试将其转换为Skyfield Python代码。这是我想出的：

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)

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)

它产生月食发生的时间的向量，然后对月食的总量进行评级：

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

日食时间在美国宇航局巨大的星历表表中给出的时间的三秒之内：

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