python - 我与 skyfield 的木星/土星合相计算与维基百科不匹配

Question

使用 Python-Skyfield 计算即将到来的木星和土星合相。

Wikipedia Great conjunction times (1800 to 2100)

使用赤经:

• 2020 年 12 月 21 日 13:22:00 UTC - 维基百科。
• 2020 年 12 月 21 日 13:34:33 UTC - 我的计算。

使用黄经:

• 2020 年 12 月 21 日 18:37:31 UTC - 维基百科
• 2020 年 12 月 21 日 18:20:40 UTC - 我的计算。

from skyfield.api import load, tau, pi
from skyfield.almanac import find_discrete

planets = load('de421.bsp')
sun = planets['sun']
earth = planets['earth']
jupiter = planets['jupiter barycenter']
saturn = planets['saturn barycenter']

ts = load.timescale(builtin=True)

def longitude_difference(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()
    _, lon1, _ = s.ecliptic_latlon()
    _, lon2, _ = j.ecliptic_latlon()
    return (lon1.degrees - lon2.degrees) > 0

def longitude_difference1(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()

    jRa, _, _ = j.radec()
    sRa, _, _ = s.radec()
    return (sRa._degrees - jRa._degrees) > 0


longitude_difference.rough_period = 300.0
longitude_difference1.rough_period = 300.0

print()
print("Great conjunction in ecliptic longitude:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), longitude_difference)
for ti in t:
    print(t.utc_jpl())

print()
print("Great conjunction in right ascension:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), longitude_difference1)
for ti in t:
    print(t.utc_jpl())

我是 Skyfield 的新手，非常感谢任何帮助。

Answer 1

你的猜想代码看起来很好，我怀疑它的结果比维基百科的要好得多——查看版本历史，甚至不清楚他们的数字从何而来，而且无法轻松进行未归因的天文计算在不知道他们从哪些星历表和软件中得出的情况下进行双重检查。

我在下面附上了你的求解器的一个稍微改进的版本。以下是我推荐的调整:

• 在计算坐标时传递 epoch='date'，因为合相和冲日传统上是根据它们发生的年份的天球定义的，而不是标准的 J2000 天球。
• 在对日期进行任何计算之前，强制将它们置于从零到整圈(360 度或 24 小时)的范围内。否则，每当赤经或经度之一恰好越过 0°/0h 并改变符号时，您将看到减法结果跳跃 ±360°/24h。这会给你“幻影合相”，其中行星并没有真正交换位置，而只是改变了它们返回的角度的符号。 (例如:从 -69° 跳到 291° 真的是天上没有任何运动。)
• 请注意，我们的两个脚本还应该找到木星和土星何时彼此穿过天空，因为它们的差异符号也应该在该点翻转。
• 如果您追踪到的任何来源都引用了最接近时刻或行星之间的角度，我已将其添加进去。

这是输出，与您的非常接近，但再次与那些无法解释的未记入维基百科的旧数字不太一致:

Great conjunction in ecliptic longitude:
['A.D. 2020-Dec-21 18:20:37.5144 UT']

Great conjunction in right ascension:
['A.D. 2020-Dec-21 13:32:04.9486 UT']

Great conjunction at closest approach:
A.D. 2020-Dec-21 18:21:00.2161 UT - 0.1018 degrees

和脚本:

from skyfield.api import load
from skyfield.searchlib import find_discrete, find_minima

planets = load('de421.bsp')
sun = planets['sun']
earth = planets['earth']
jupiter = planets['jupiter barycenter']
saturn = planets['saturn barycenter']

ts = load.timescale(builtin=True)

def longitude_difference(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()
    _, lon1, _ = s.ecliptic_latlon(epoch='date')
    _, lon2, _ = j.ecliptic_latlon(epoch='date')
    return (lon1.degrees - lon2.degrees) % 360.0 > 180.0

def ra_difference(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()

    jRa, _, _ = j.radec(epoch='date')
    sRa, _, _ = s.radec(epoch='date')
    return (sRa.hours - jRa.hours) % 24.0 > 12.0

def separation(t):
    e = earth.at(t)
    j = e.observe(jupiter).apparent()
    s = e.observe(saturn).apparent()

    return j.separation_from(s).degrees

longitude_difference.step_days = 30.0
ra_difference.step_days = 30.0
separation.step_days = 30.0

print()
print("Great conjunction in ecliptic longitude:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), longitude_difference)
for ti in t:
    print(t.utc_jpl())

print()
print("Great conjunction in right ascension:")
t, b = find_discrete(ts.utc(2020), ts.utc(2021), ra_difference)
for ti in t:
    print(t.utc_jpl())

print()
print("Great conjunction at closest approach:")
t, b = find_minima(ts.utc(2020, 6), ts.utc(2021), separation)
for ti, bi in zip(t, b):
    print('{} - {:.4f} degrees'.format(ti.utc_jpl(), bi))