确定天空中的月食答案

【问题标题】：Determining lunar eclipse in skyfield确定天空中的月食
【发布时间】：2021-02-15 20:42:38
【问题描述】：

我得到了一个 UTC 日期列表，所有时间都转换为 00:00。

我想确定某一天（即过去 24 小时）是否发生（月食）日食

考虑到python sn-p

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 小时内是否发生了日食

检查第一个角度是否足够接近 180（简单）
检查二度是否足够接近 0（不是那么简单？）

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

因此，我的问题是

有人可以提供一个函数来确定在给定的日期 t 是否发生月食吗？

编辑。编辑此问题以反映下面 cmets 中留下的 Brandon Rhodes 的反馈。

【问题讨论】：

如果你询问月球在“t 减一天”的位置，第一个问题——它需要你知道月球在天空中行进的速度，它会变化——会消失吗？检查该位置是否在 180° 的另一侧？
@BrandonRhodes 谢谢，效果很好。类似的逻辑也可以应用于日食吗？两个平面的角度是否单调？如果不是，那么确定日食的明智方法是什么？
——我从来没有写过任何日食发现例程，唉（这就是为什么我只添加了一个评论而没有尝试回答你的问题）。我知道它涉及地球阴影在太空中的 3D 形状，但是由于从月球的角度来看，阴影可以近似为一对圆锥体（本影和半影），所以图表似乎总是将阴影近似为一对月球穿过的圆圈——穿过这两个圆锥的平面截面。所以我猜他们会测量月球相对于远处那些阴影圈的位置？
@BrandonRhodes 感谢您的总体想法。这告诉我我没有足够的经验来自己编码。因此，我只是要重构这个问题并增加一点赏金。
由于日子在倒计时，我会在评论中询问我的答案是否接近满足您的赏金条款，或者是否还有任何关键问题悬而未决。谢谢！

标签： python astronomy skyfield

【解决方案1】：

我刚刚浏览了天文年历补充说明的第 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, 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]

它的日食时间与美国宇航局巨大的月球星历表中给出的时间相差不到 3 秒：

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

【讨论】：