在Skyfield中确定月食

Question

我被给定了一个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之内确定日食是否发生：

1. 检查第一个角度是否足够接近180（简单）
2. 检查第二个角度是否足够接近0（不太容易？）

根据评论中的答案，解决第二个问题并不是一个简单的测试角度是否接近于0。

因此，我的问题是：

有人能提供一个函数来确定给定日期t是否发生了月食吗？

编辑。 这个问题已经被编辑以反映以下Brandon Rhodes在下面留下的反馈。

- Jernej

---

如果你询问“t减一天”时的月球位置，并检查该位置是否在180度的另一侧，那么第一个问题——需要你知道月球在天空中移动的速度（这是变化的）——是否会消失？ - Brandon Rhodes

@BrandonRhodes 谢谢，那很完美。类似的逻辑是否也适用于日食呢？这两个平面的角度是否是单调的？如果不是，有什么合理的方法来确定日食吗？ - Jernej

很遗憾，我从未编写过任何查找日食的程序（这就是为什么我只添加了一条评论而没有尝试回答你的问题）。我知道它涉及到地球阴影在空间中的三维形状，但由于从月球的角度来看，阴影可以近似为一对圆锥体（本影和半影），因此图表似乎总是将阴影近似为一对圆形，月球穿过这两个圆形的平面截面。所以我猜他们测量月球相对于其距离处的那些阴影圆的位置？ - Brandon Rhodes

@BrandonRhodes 感谢您提供的一般想法。这告诉我我没有足够的经验来自己编写代码。因此，我只是要重构这个问题并增加一点赏金。 - Jernej

由于时间正在倒计时，我在此评论中询问我的答案是否接近满足您的悬赏条件，或者是否仍然存在任何重要问题。谢谢！ - Brandon Rhodes

Answer 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]

它的月食时间与NASA的大型月球历表中给出的时间相差不超过3秒：

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

- Brandon Rhodes