首页 > 解决方案 > 我与天域的木星/土星合相计算与维基百科不匹配

问题描述

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

维基百科大合相时代(1800 到 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 的新手,因此不胜感激。

标签: pythonastronomyskyfield

解决方案


你的猜想代码看起来非常好,我怀疑它的结果比维基百科的结果要好得多——查看版本历史,甚至不清楚它们的数字来自哪里,而且未归属的天文学计算不容易被翻倍——在不知道他们从哪个星历表和软件中获得它们的情况下进行了检查。

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

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

推荐阅读