我的 Jupiter/Saturn 与 skyfield 的结合计算与维基百科不匹配

My Jupiter/Saturn conjunction calculation with skyfield doesn't match wikipedia

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

Wikipedia Great conjunction times (1800 to 2100)

使用赤经:

使用黄经:

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 的新手,非常感谢您的帮助。

我已尽力避免对这个问题有基于意见的答案的可能性,并在互联网上进行了查找。发现很难找到任何我可以信任的相关信息,所以我列举了这些帖子(维基百科除外):

timeanddate 表示确切时间是 12 月 21 日 18:20 UTC,这与您计算的一样

winstars 表示行星最接近角度的时间为 18:25 UTC 并且他们提到合相将发生在 13:30 UTC,不知道是不是第一次

不确定 this 的相关性如何,但这里的合相据说在格林威治标准时间 17:32 时为 6.2 度,因此 18: 32 协调世界时

我能找到的最相关的来源是 in the sky, where the time was estimated exactly to 13:24 UTC., based on calculations on data by Jet Propulsion Laboratory - source code can be checked here (c)

您可以看到,大多数情况下并没有同时使用这两种计算方式,而且时间各不相同。原因是在此类计算中,您需要非常长的浮点数才能获得最佳精度。由于您受到所用机器的限制,精度并不完美。如@bad_coder所建议,你可能在Astronomy stack exchange.

中得到更好的答案

你的猜想代码看起来很好,我怀疑它的结果比维基百科的结果好得多——看版本历史,甚至不清楚他们的数字从何而来,未归属的天文计算也不能'在不知道他们从哪个星历表和软件中导出它们的情况下,很容易被双重检查。

我在下方附上了您的求解器的略微改进版本。以下是我推荐的调整:

  • 在计算坐标时通过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))