问题描述
我获得了以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小时之内是否发生了日食
- 检查第一个角度是否足够接近180度(容易)
- 检查第二学位是否足够接近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]
日食时间在美国宇航局巨大的星历表表中给出的时间的三秒之内: