使用 skyfield 的脚本的无效结果
Invalid results with script that uses skyfield
我正在探索 Brandon Rhodes 出色的软件 Skyfield 的可能性。我制作了一个脚本来计算随机对象之间的 Right Ascension 连词。我使用以下脚本：
from skyfield import almanac
from skyfield.searchlib import find_maxima, find_minima, find_discrete
from skyfield.api import Star, load
from datetime import datetime, date,timedelta
import pytz
planets = load('de430t.bsp')
earth = planets['earth']
x = [
['Aldebaran',[4, 35, 55.2],[16, 30, 33]],
['Regulus',[10, 8, 22.3],[11, 58, 2]],
['Pollux',[7, 45, 18.9],[28, 1, 34]],
['Antares',[16, 29, 24.4],[-26, 25, 55]],
]
ts = load.timescale(builtin=True)
t = ts.now()
tzn = 'Europe/Amsterdam'
tz = pytz.timezone(tzn)
now = datetime(2020, 1, 1, 0, 0, 0)
t0 = ts.utc(tz.localize(now))
t1 = ts.utc(tz.localize(now) + timedelta(days=+365))
def difference(t):
e = earth.at(t)
ra11, dec11, distance = e.observe(object).radec()
ra12, dec12, distance2 = e.observe(planets[target1]).radec()
diff = ra11.hours - ra12.hours
return diff >= 0
difference.rough_period = 1.0
for count in range (len(x)):
object = Star(ra_hours=(x[count][1][0], x[count][1][1], x[count][1][2]),dec_degrees=(x[count][2][0], x[count][2][1], x[count][2][2]))
target1 = 'venus'
t, b = find_discrete(t0, t1, difference)
if len(t) > 0:
print (f"{x[count][0]} and {target1}")
for i, (a, b) in enumerate(zip(t, b)):
e = earth.at(a)
ra11, dec11, distance = e.observe(object).radec('date')
ra12, dec12, distance2 = e.observe(planets[target1]).radec('date')
print (f"Diff: {ra11.hours - ra12.hours}, ra_{x[count][0]}: {ra11.hours}, ra_{target1}: {ra12.hours}")
print(f"{a.utc_iso()},{dec11._degrees - dec12._degrees}")
print ("")
我相信这是在计算两个对象具有相同 RA 的时间实例。
不幸的是我得到了这些结果：
Aldebaran and venus
Diff: 4.600705365706519, ra_Aldebaran: 4.6176105612629375, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,16.962942031863825
Diff: -0.0014205156698450239, ra_Aldebaran: 4.617748605529588, ra_venus: 4.619169121199433
2020-04-17T20:25:49Z,-10.136618155596008
Diff: -0.000670093218860579, ra_Aldebaran: 4.617892691238783, ra_venus: 4.618562784457644
2020-06-08T07:56:08Z,-4.921187478324768
Diff: -0.0001286749609095139, ra_Aldebaran: 4.618000948409605, ra_venus: 4.618129623370515
2020-07-12T06:44:16Z,-0.962286810883981
Regulus and venus
Diff: 10.140247344361857, ra_Regulus: 10.157152539918275, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,12.28436748615726
Diff: 5.852858068422506e-06, ra_Regulus: 10.157702949500333, ra_venus: 10.157697096642265
2020-10-02T23:42:45Z,0.0903429562716429
Pollux and venus
Diff: 7.758742204719277, ra_Pollux: 7.7756474002756955, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,28.391501492522927
Diff: 0.001226225278400328, ra_Pollux: 7.776229220287632, ra_venus: 7.775002995009232
2020-09-01T16:39:22Z,8.682000412217121
Antares and venus
Diff: 16.493491164600684, ra_Antares: 16.510396360157102, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,-26.059118330110437
Diff: 0.000832014094154232, ra_Antares: 16.51126040187071, ra_venus: 16.510428387776557
2020-12-23T00:34:39Z,-5.652225571050259
以“Diff”开头的行是监控输出有效性的行。
Diff 代表计算出的 RA 差异。它应该接近于零。其他两个值是两个对象的赤经。它们应该非常相似。
第二行是我想要的结果，它是计算的时间和对象之间的度数距离。
正如您所看到的，出于某种原因，对于每组对象，我得到了时间实例的无效结果：2020-02-07T21:20:06Z，并且该实例的差值肯定不接近于零。如果我将对象 venus 更改为 moon 它会变得更糟，因为每秒结果都是无效的。
我根据 Skychart / Cartes du Ciel 软件和那些结帐检查了其他结果。
我不知道这里出了什么问题。有人可以帮我吗？
好问题！我应该在 https://rhodesmill.org/skyfield/searches.html 中添加一个新部分来解释在减去两个经度或赤经时看到的这种常见行为。解开这个谜团的关键是观察角度差在输出中作为幻影合相出现的某一时刻发生了什么。我附上了一个脚本，它为您在金星和毕宿五之间打印的第一个事件打印此内容：
2020-02-07 Difference: -19.33880215224481 Venus RA: 23.93746881891146
2020-02-08 Difference: 4.5908654129343995 Venus RA: 0.007801253732248779
角度差在 -19.3 和 4.6 之间跳跃，这应该会立即引起我们的怀疑，因为它们只是完全相同角度的两个不同名称！您可以通过将 24.0 添加到 -19.3 来确认这一点，您将得到一个非常接近 4.6 的角度（给出或考虑金星在一天内完成的实际运动）。
为什么结果会在两个别名之间跳跃，因为天空中的 angular 差异完全相同？
答案就在上面打印的第二个事实中：金星的 RA。不连续恰好发生在金星恰好穿过 0h 赤经的那一刻！尽管 23.93746881891146 和 0.007801253732248779 几乎是相同的角度，但它们相差 24.0，因为它们横跨天空中我们重新命名赤经的位置。
我下面的脚本也显示了一个说明情况的情节：
你可以看到，在顶部的图中，正是在金星将其 RA 重置回零的确切时刻，RA 差异从一个别名到另一个别名跳跃了 24.0 小时，用于相同的四和 - RA 半小时的差异。
解决方案？
Angular 差异需要限制在 [-12h, +12h) 之类的范围内，以强制为每个可能的 angular 值选择一个首选别名。在 Python 中，正如您在下面的脚本中所见，典型的操作是：
(ra1.hours - ra2.hours + 12.0) % 24.0 - 12.0
这在上面的第二个图中显示为“改进的差异”，它不仅正确地隐藏了 2 月 7 日的不连续性，不再让它看起来像一个事件，而且它现在正确地识别了金星和Aldebaran 在 2020 年底（在图的右边缘）出现，之前仅显示为超过 -12.0 的差异，但现在作为 angular 差异的关键时刻闪耀。
最后，此脚本会检查反对意见并将其从搜索结果中过滤掉。您还会发现对您的 Python 代码进行一些可能的调整，您可能会考虑在 Python 中继续编码。这里是：
from skyfield.searchlib import find_maxima, find_minima, find_discrete
from skyfield.api import Star, load
from datetime import datetime, date,timedelta
import pytz
planets = load('de430t.bsp')
earth = planets['earth']
stars = [
['Aldebaran', (4, 35, 55.2), (16, 30, 33)],
['Regulus', (10, 8, 22.3), (11, 58, 2)],
['Pollux', (7, 45, 18.9), (28, 1, 34)],
['Antares', (16, 29, 24.4), (-26, 25, 55)],
]
ts = load.timescale(builtin=True)
t0 = ts.utc(2020, 1, 1)
t1 = ts.utc(2021, 1, 1)
# Exploring the first bad result for Venus and Aldebaran.
star = Star(ra_hours=stars[0][1], dec_degrees=stars[0][2])
for utc in (2020, 2, 7), (2020, 2, 8):
t = ts.utc(*utc)
ra1, _, _ = planets['earth'].at(t).observe(star).radec()
ra2, _, _ = planets['earth'].at(t).observe(planets['venus']).radec()
print(t.utc_strftime('%Y-%m-%d'),
'Difference:', ra1.hours - ra2.hours,
'Venus RA:', ra2.hours)
# Plot showing how to protect an angular difference against discontinuity.
import matplotlib.pyplot as plt
t = ts.utc(2020, 1, range(365))
e = planets['earth'].at(t)
star = Star(ra_hours=stars[0][1], dec_degrees=stars[0][2])
ra1, _, _, = e.observe(star).radec()
ra2, _, _, = e.observe(planets['venus']).radec()
fig, (ax, ax2) = plt.subplots(2, 1)
ax.plot(t.J, ra2.hours, label='Venus RA')
ax.plot(t.J, ra1.hours - ra2.hours, label='RA difference')
ax.set(xlabel='Year', ylabel='Hours')
ax.grid()
ax.legend()
ax2.plot(t.J, ra2.hours, label='Venus RA')
ax2.plot(t.J, (ra1.hours - ra2.hours + 12.0) % 24.0 - 12.0,
label='Improved difference')
ax2.set(xlabel='Year', ylabel='Hours')
ax2.grid()
ax2.legend()
fig.savefig('tmp.png')
# So we need to force the difference into the range [-12 hours .. +12 hours]
def difference(t):
e = earth.at(t)
ra11, dec11, distance = e.observe(object).radec()
ra12, dec12, distance2 = e.observe(planets[target1]).radec()
diff = (ra11.hours - ra12.hours + 12.0) % 24.0 - 12.0
return diff >= 0
difference.rough_period = 1.0
for name, ra_hms, dec_dms in stars:
object = Star(ra_hours=ra_hms, dec_degrees=dec_dms)
target1 = 'venus'
t, b = find_discrete(t0, t1, difference)
if len(t) == 0:
break
print (f"{name} and {target1}")
for a, b in zip(t, b):
e = earth.at(a)
ra1, dec1, _ = e.observe(object).radec('date')
ra2, dec2, _ = e.observe(planets[target1]).radec('date')
if abs(ra1.hours - ra2.hours) > 6.0:
continue # ignore oppositions
print(f"Diff: {ra1.hours - ra2.hours:.4f}, ra_{name}: {ra1.hours}, ra_{target1}: {ra2.hours}")
print(f"{a.utc_iso()},{dec1._degrees - dec2._degrees}")
print()
现在打印的事件是：
Aldebaran and venus
Diff: -0.0014, ra_Aldebaran: 4.617748605529591, ra_venus: 4.619169121320681
2020-04-17T20:25:49Z,-10.136618155920797
Diff: -0.0007, ra_Aldebaran: 4.617892691238804, ra_venus: 4.618562784294269
2020-06-08T07:56:08Z,-4.921187477025711
Diff: -0.0001, ra_Aldebaran: 4.618000948409604, ra_venus: 4.618129623302464
2020-07-12T06:44:16Z,-0.9622868107603999
Regulus and venus
Diff: 0.0000, ra_Regulus: 10.157702949500333, ra_venus: 10.157697096607553
2020-10-02T23:42:45Z,0.09034295610918264
Pollux and venus
Diff: 0.0012, ra_Pollux: 7.776229220287636, ra_venus: 7.775002995218732
2020-09-01T16:39:22Z,8.68200041253418
Antares and venus
Diff: 0.0008, ra_Antares: 16.511260401870718, ra_venus: 16.51042838802181
2020-12-23T00:34:39Z,-5.652225570418828
我相信这可以解决并纠正您的问题！
我正在探索 Brandon Rhodes 出色的软件 Skyfield 的可能性。我制作了一个脚本来计算随机对象之间的 Right Ascension 连词。我使用以下脚本：
from skyfield import almanac
from skyfield.searchlib import find_maxima, find_minima, find_discrete
from skyfield.api import Star, load
from datetime import datetime, date,timedelta
import pytz
planets = load('de430t.bsp')
earth = planets['earth']
x = [
['Aldebaran',[4, 35, 55.2],[16, 30, 33]],
['Regulus',[10, 8, 22.3],[11, 58, 2]],
['Pollux',[7, 45, 18.9],[28, 1, 34]],
['Antares',[16, 29, 24.4],[-26, 25, 55]],
]
ts = load.timescale(builtin=True)
t = ts.now()
tzn = 'Europe/Amsterdam'
tz = pytz.timezone(tzn)
now = datetime(2020, 1, 1, 0, 0, 0)
t0 = ts.utc(tz.localize(now))
t1 = ts.utc(tz.localize(now) + timedelta(days=+365))
def difference(t):
e = earth.at(t)
ra11, dec11, distance = e.observe(object).radec()
ra12, dec12, distance2 = e.observe(planets[target1]).radec()
diff = ra11.hours - ra12.hours
return diff >= 0
difference.rough_period = 1.0
for count in range (len(x)):
object = Star(ra_hours=(x[count][1][0], x[count][1][1], x[count][1][2]),dec_degrees=(x[count][2][0], x[count][2][1], x[count][2][2]))
target1 = 'venus'
t, b = find_discrete(t0, t1, difference)
if len(t) > 0:
print (f"{x[count][0]} and {target1}")
for i, (a, b) in enumerate(zip(t, b)):
e = earth.at(a)
ra11, dec11, distance = e.observe(object).radec('date')
ra12, dec12, distance2 = e.observe(planets[target1]).radec('date')
print (f"Diff: {ra11.hours - ra12.hours}, ra_{x[count][0]}: {ra11.hours}, ra_{target1}: {ra12.hours}")
print(f"{a.utc_iso()},{dec11._degrees - dec12._degrees}")
print ("")
我相信这是在计算两个对象具有相同 RA 的时间实例。
不幸的是我得到了这些结果：
Aldebaran and venus
Diff: 4.600705365706519, ra_Aldebaran: 4.6176105612629375, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,16.962942031863825
Diff: -0.0014205156698450239, ra_Aldebaran: 4.617748605529588, ra_venus: 4.619169121199433
2020-04-17T20:25:49Z,-10.136618155596008
Diff: -0.000670093218860579, ra_Aldebaran: 4.617892691238783, ra_venus: 4.618562784457644
2020-06-08T07:56:08Z,-4.921187478324768
Diff: -0.0001286749609095139, ra_Aldebaran: 4.618000948409605, ra_venus: 4.618129623370515
2020-07-12T06:44:16Z,-0.962286810883981
Regulus and venus
Diff: 10.140247344361857, ra_Regulus: 10.157152539918275, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,12.28436748615726
Diff: 5.852858068422506e-06, ra_Regulus: 10.157702949500333, ra_venus: 10.157697096642265
2020-10-02T23:42:45Z,0.0903429562716429
Pollux and venus
Diff: 7.758742204719277, ra_Pollux: 7.7756474002756955, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,28.391501492522927
Diff: 0.001226225278400328, ra_Pollux: 7.776229220287632, ra_venus: 7.775002995009232
2020-09-01T16:39:22Z,8.682000412217121
Antares and venus
Diff: 16.493491164600684, ra_Antares: 16.510396360157102, ra_venus: 0.016905195556418524
2020-02-07T21:20:06Z,-26.059118330110437
Diff: 0.000832014094154232, ra_Antares: 16.51126040187071, ra_venus: 16.510428387776557
2020-12-23T00:34:39Z,-5.652225571050259
以“Diff”开头的行是监控输出有效性的行。 Diff 代表计算出的 RA 差异。它应该接近于零。其他两个值是两个对象的赤经。它们应该非常相似。 第二行是我想要的结果，它是计算的时间和对象之间的度数距离。 正如您所看到的，出于某种原因，对于每组对象，我得到了时间实例的无效结果：2020-02-07T21:20:06Z，并且该实例的差值肯定不接近于零。如果我将对象 venus 更改为 moon 它会变得更糟，因为每秒结果都是无效的。 我根据 Skychart / Cartes du Ciel 软件和那些结帐检查了其他结果。
我不知道这里出了什么问题。有人可以帮我吗？
好问题！我应该在 https://rhodesmill.org/skyfield/searches.html 中添加一个新部分来解释在减去两个经度或赤经时看到的这种常见行为。解开这个谜团的关键是观察角度差在输出中作为幻影合相出现的某一时刻发生了什么。我附上了一个脚本，它为您在金星和毕宿五之间打印的第一个事件打印此内容：
2020-02-07 Difference: -19.33880215224481 Venus RA: 23.93746881891146
2020-02-08 Difference: 4.5908654129343995 Venus RA: 0.007801253732248779
角度差在 -19.3 和 4.6 之间跳跃，这应该会立即引起我们的怀疑，因为它们只是完全相同角度的两个不同名称！您可以通过将 24.0 添加到 -19.3 来确认这一点，您将得到一个非常接近 4.6 的角度（给出或考虑金星在一天内完成的实际运动）。
为什么结果会在两个别名之间跳跃，因为天空中的 angular 差异完全相同？
答案就在上面打印的第二个事实中：金星的 RA。不连续恰好发生在金星恰好穿过 0h 赤经的那一刻！尽管 23.93746881891146 和 0.007801253732248779 几乎是相同的角度，但它们相差 24.0，因为它们横跨天空中我们重新命名赤经的位置。
我下面的脚本也显示了一个说明情况的情节：
你可以看到，在顶部的图中，正是在金星将其 RA 重置回零的确切时刻，RA 差异从一个别名到另一个别名跳跃了 24.0 小时，用于相同的四和 - RA 半小时的差异。
解决方案？
Angular 差异需要限制在 [-12h, +12h) 之类的范围内，以强制为每个可能的 angular 值选择一个首选别名。在 Python 中，正如您在下面的脚本中所见，典型的操作是：
(ra1.hours - ra2.hours + 12.0) % 24.0 - 12.0
这在上面的第二个图中显示为“改进的差异”，它不仅正确地隐藏了 2 月 7 日的不连续性，不再让它看起来像一个事件，而且它现在正确地识别了金星和Aldebaran 在 2020 年底（在图的右边缘）出现，之前仅显示为超过 -12.0 的差异，但现在作为 angular 差异的关键时刻闪耀。
最后，此脚本会检查反对意见并将其从搜索结果中过滤掉。您还会发现对您的 Python 代码进行一些可能的调整，您可能会考虑在 Python 中继续编码。这里是：
from skyfield.searchlib import find_maxima, find_minima, find_discrete
from skyfield.api import Star, load
from datetime import datetime, date,timedelta
import pytz
planets = load('de430t.bsp')
earth = planets['earth']
stars = [
['Aldebaran', (4, 35, 55.2), (16, 30, 33)],
['Regulus', (10, 8, 22.3), (11, 58, 2)],
['Pollux', (7, 45, 18.9), (28, 1, 34)],
['Antares', (16, 29, 24.4), (-26, 25, 55)],
]
ts = load.timescale(builtin=True)
t0 = ts.utc(2020, 1, 1)
t1 = ts.utc(2021, 1, 1)
# Exploring the first bad result for Venus and Aldebaran.
star = Star(ra_hours=stars[0][1], dec_degrees=stars[0][2])
for utc in (2020, 2, 7), (2020, 2, 8):
t = ts.utc(*utc)
ra1, _, _ = planets['earth'].at(t).observe(star).radec()
ra2, _, _ = planets['earth'].at(t).observe(planets['venus']).radec()
print(t.utc_strftime('%Y-%m-%d'),
'Difference:', ra1.hours - ra2.hours,
'Venus RA:', ra2.hours)
# Plot showing how to protect an angular difference against discontinuity.
import matplotlib.pyplot as plt
t = ts.utc(2020, 1, range(365))
e = planets['earth'].at(t)
star = Star(ra_hours=stars[0][1], dec_degrees=stars[0][2])
ra1, _, _, = e.observe(star).radec()
ra2, _, _, = e.observe(planets['venus']).radec()
fig, (ax, ax2) = plt.subplots(2, 1)
ax.plot(t.J, ra2.hours, label='Venus RA')
ax.plot(t.J, ra1.hours - ra2.hours, label='RA difference')
ax.set(xlabel='Year', ylabel='Hours')
ax.grid()
ax.legend()
ax2.plot(t.J, ra2.hours, label='Venus RA')
ax2.plot(t.J, (ra1.hours - ra2.hours + 12.0) % 24.0 - 12.0,
label='Improved difference')
ax2.set(xlabel='Year', ylabel='Hours')
ax2.grid()
ax2.legend()
fig.savefig('tmp.png')
# So we need to force the difference into the range [-12 hours .. +12 hours]
def difference(t):
e = earth.at(t)
ra11, dec11, distance = e.observe(object).radec()
ra12, dec12, distance2 = e.observe(planets[target1]).radec()
diff = (ra11.hours - ra12.hours + 12.0) % 24.0 - 12.0
return diff >= 0
difference.rough_period = 1.0
for name, ra_hms, dec_dms in stars:
object = Star(ra_hours=ra_hms, dec_degrees=dec_dms)
target1 = 'venus'
t, b = find_discrete(t0, t1, difference)
if len(t) == 0:
break
print (f"{name} and {target1}")
for a, b in zip(t, b):
e = earth.at(a)
ra1, dec1, _ = e.observe(object).radec('date')
ra2, dec2, _ = e.observe(planets[target1]).radec('date')
if abs(ra1.hours - ra2.hours) > 6.0:
continue # ignore oppositions
print(f"Diff: {ra1.hours - ra2.hours:.4f}, ra_{name}: {ra1.hours}, ra_{target1}: {ra2.hours}")
print(f"{a.utc_iso()},{dec1._degrees - dec2._degrees}")
print()
现在打印的事件是：
Aldebaran and venus
Diff: -0.0014, ra_Aldebaran: 4.617748605529591, ra_venus: 4.619169121320681
2020-04-17T20:25:49Z,-10.136618155920797
Diff: -0.0007, ra_Aldebaran: 4.617892691238804, ra_venus: 4.618562784294269
2020-06-08T07:56:08Z,-4.921187477025711
Diff: -0.0001, ra_Aldebaran: 4.618000948409604, ra_venus: 4.618129623302464
2020-07-12T06:44:16Z,-0.9622868107603999
Regulus and venus
Diff: 0.0000, ra_Regulus: 10.157702949500333, ra_venus: 10.157697096607553
2020-10-02T23:42:45Z,0.09034295610918264
Pollux and venus
Diff: 0.0012, ra_Pollux: 7.776229220287636, ra_venus: 7.775002995218732
2020-09-01T16:39:22Z,8.68200041253418
Antares and venus
Diff: 0.0008, ra_Antares: 16.511260401870718, ra_venus: 16.51042838802181
2020-12-23T00:34:39Z,-5.652225570418828
我相信这可以解决并纠正您的问题！