Skyfield 测试 2 星修复数学时出现巨大错误

Question

我正在尝试使用 Skyfield python lib 和 Celestial 编程网站 https://www.celestialprogramming.com/snippets/twoStarFix.html 进行 2 星修复数学的简单测试 因此，我正在使用 Skyfield 在时间 t 获取对来自依巴谷的 2 颗恒星的想象观测数据。然后我使用天体编程网站的数学计算回观察者的地理位置。我遇到了一个大错误，例如：

数据：49N，25E 结果：48.966578719047284、25.25305554808433

问题出在哪里？代码中是否有错误？感谢您的帮助！

代码本身：

# Greg Miller ([email protected]) 2023 http://www.celestialprogramming.com/

from skyfield.api import Topos, N, E, W, wgs84, load, Star
from skyfield.data import hipparcos
from skyfield.positionlib import position_of_radec
from math import pi, sin, cos, acos, atan2, sqrt

rad = pi/180.0

# Angles returned in radians
def rectToPolar(x, y, z):
    r = sqrt(x * x + y * y + z * z);
    lon = atan2(y, x);
    lat = acos(z / r);
    lat = pi/2 - lat # Make sure lat is in range +/-90deg
    return [lat, lon, r]

# star[0]=GP Longitude, star[1]=Declination, star[2]=Zenith Distance
def sightReduce(star1, star2):
    # From "An Analytical Solution of the Two Star Sight Problem of Celestial Navigation." 
    # James A. Van Allen

    a1 = cos(star1[0])*cos(star1[1])
    a2 = cos(star2[0])*cos(star2[1])

    b1=sin(star1[0])*cos(star1[1])
    b2=sin(star2[0])*cos(star2[1])

    c1=sin(star1[1])
    c2=sin(star2[1])

    p1=cos(star1[2])
    p2=cos(star2[2])

    # Eq 3
    A=a1*b2 - a2*b1
    B=b2*c1 - b1*c2
    C=b2*p1 - b1*p2
    D=a1*c2 - a2*c1
    E=b1*c2 - b2*c1
    F=c2*p1 - c1*p2

    # Eq 4 substituded in to eq 5
    # x^2(1+A^2/B^2 + D^2/E^2) + x(-2AC/B^2 - 2FD/E^2) + C^2/B^2 + F^2/E^2 - 1

    # Quadratic solve
    a=(1+(A*A)/(B*B)+(D*D)/(E*E))
    b=((-2*A*C)/(B*B) - (2*F*D)/(E*E))
    c=(C*C)/(B*B) + (F*F)/(E*E) - 1

    x1=(-b+sqrt(b*b - 4*a*c))/(2*a)
    x2=(-b-sqrt(b*b - 4*a*c))/(2*a)

    # eq  4
    y1=(F-D*x1)/E
    y2=(F-D*x2)/E

    z1=(C-A*x1)/B
    z2=(C-A*x2)/B

    s1=rectToPolar(x1,y1,z1)
    s2=rectToPolar(x2,y2,z2)

    return [s1[0]/rad, s1[1]/rad, s2[0]/rad, s2[1]/rad]

with load.open(hipparcos.URL) as f:
    df = hipparcos.load_dataframe(f)
planets = load('de421.bsp') # Load the JPL ephemeris DE421 (covers 1900-2050).
earth = planets['earth']

origLat, origLng = 49, 25

(txtLat, txtLatS) = (-origLat, "S") if (origLat < 0) else (origLat, "N")
(txtLng, txtLngS) = (-origLng, "W") if (origLng < 0) else (origLng, "E")
te = Topos(str(txtLat)+" "+txtLatS, str(txtLng)+" "+txtLngS)# elevation_m=300

tx = 1710621420+27
t = load.timescale().utc(1970, 1, 1, 0, 0, tx)
print(t.utc_strftime())
body = [4427, 67301]

def getHip(body):
    star = Star.from_dataframe(df.loc[body])
    pps = (earth + te).at(t).observe(star).apparent()
    alt0, az0, distance = pps.altaz()
    pps = (earth).at(t).observe(star)
    ra0, dec0, distance = pps.radec('date')
    xppq = position_of_radec(ra0.radians*12/pi, dec0.radians*180/pi, t=t, center=399)
    xq2 = wgs84.geographic_position_of(xppq)
    return [xq2.longitude.radians, dec0.radians, pi/2-alt0.radians] # [0]=GP Longitude, [1]=Declination, [2]=Zenith Distance

print(origLat, origLng)
print(sightReduce(getHip(body[0]), getHip(body[1])))

已经尝试过这个Python代码示例，但出现了一个大错误。

Answer 1

听起来您在 Skyfield 中测试两星修复的数学时遇到了重大错误。在使用天文导航或类似应用时，准确性至关重要。以下是解决此问题的方法：

1. 检查输入数据：确保您提供给 Skyfield 的输入数据准确无误。这包括恒星或天体的坐标、观测时间以及任何其他相关参数。

2. 检查代码逻辑： 仔细检查用于计算 Skyfield 中两星定位的代码逻辑。查找数学计算或算法中任何潜在的错误或不准确之处。

3. 比较结果：将从 Skyfield 获得的结果与已知值或其他来源的结果进行比较。这可以帮助识别差异并查明错误根源。

4. 调试：使用调试技术逐步执行代码并跟踪计算每一步的变量值。这可以帮助您确定错误发生的位置以及原因。

5. 查阅文档： 请参阅 Skyfield 提供的文档和资源，了解如何正确实施两星修正计算。请注意提到的任何具体要求或注意事项。

6. 寻求社区支持：如果您仍然无法解决错误，请考虑联系 Skyfield 社区或论坛寻求帮助。其他用户或开发人员可能也遇到过类似的问题，并且可以提供有价值的见解或解决方案。

7. 报告错误：如果您认为自己已发现 Skyfield 中的错误，请考虑向开发人员报告。提供有关您遇到的错误的详细信息，包括任何相关的代码片段、输入数据以及预期结果与实际结果。

通过遵循这些步骤并仔细检查您的代码和计算，您应该能够解决 Skyfield 中的两星修复数学中的错误并获得准确的结果。