ruby - 使用 ruby​​ 计算太阳体的 x、y、z 位置

Question

我在这里使用火星的方程式和数据 http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf

以及第四页顶部给出的偏心反常开普勒方程的解 http://murison.alpheratz.net/dynamics/twobody/KeplerIterations_summary.pdf

并通过将 get_centuries_past 中的日期修改为以下日期并查看第 E-7 页以获取火星的实际 x、y、z 坐标来检查输出(下面的示例数据，但好奇的链接: http://books.google.com/books/about/Astronomical_Almanac_for_the_Year_2013_a.html?id=7fl_-DLwJ8YC)

日期 2456320.5 是 2013, 1, 28 应该输出

x = 1.283762
y = -0.450111
z = -0.241123

日期 2456357.5 是 2013, 3, 6 应该输出

x = 1.300366
y = 0.533593
z = 0.209626

日期 2456539.500000 是 2013, 9, 4 应该输出

x = - 0.325604
y = 1.418110
z = 0.659236

我测试了均值异常方程，结果很好。但是，我无法获得一组好的 x、y、z 坐标。我一直在调整我的开普勒和坐标函数，但无法使它们与天文年历中的表格相匹配。

非常感谢任何关于解决星星位置的建议或意见。下面的代码可以放在 .rb 文件中，在命令行上运行它会输出 x、y、z 值。

def get_centuries_past_j2000()
        #second number is from DateTime.new(2000,1,1,12).amjd.to_f - 1 the modified julian date for the J2000 Epoch
        #Date.today.jd.to_f - 51544.5
        (DateTime.new(2013,1,28).amjd.to_f - 51544.5)/36525
end

class Planet
    attr_accessor :semi_major_axis, :semi_major_axis_delta, :eccentricity, :eccentricity_delta,
    :inclination, :inclination_delta, :mean_longitude, :mean_longitude_delta, :longitude_of_perihelion, 
    :longitude_of_perihelion_delta, :longitude_of_ascending_node, :longitude_of_ascending_node_delta, :time_delta
def initialize(semi_major_axis, semi_major_axis_delta, eccentricity, eccentricity_delta,
inclination, inclination_delta, mean_longitude, mean_longitude_delta, longitude_of_perihelion, 
longitude_of_perihelion_delta, longitude_of_ascending_node, longitude_of_ascending_node_delta, time_delta)
    @semi_major_axis = semi_major_axis + (semi_major_axis_delta * time_delta)
    @eccentricity = eccentricity + (eccentricity_delta * time_delta)
    @inclination = inclination + (inclination_delta * time_delta)
    @mean_longitude = mean_longitude + (mean_longitude_delta * time_delta)
    @longitude_of_perihelion = longitude_of_perihelion + (longitude_of_perihelion_delta * time_delta)
    @longitude_of_ascending_node = longitude_of_ascending_node + (longitude_of_ascending_node_delta * time_delta)
    @argument_of_perhelion = @longitude_of_perihelion - @longitude_of_ascending_node
end

def mean_anomaly
    ((@mean_longitude - @longitude_of_perihelion)%360).round(8)
end

def eccentric_anomaly
    mod_mean_anomaly = mean_anomaly
    if mod_mean_anomaly > 180
        mod_mean_anomaly = mod_mean_anomaly - 360
    elsif mod_mean_anomaly < -180
        mod_mean_anomaly = mod_mean_anomaly + 360
    end
    e34 = @eccentricity**2
    e35 = @eccentricity*e34
    e33 = Math.cos(mod_mean_anomaly*Math::PI/180)
    mod_mean_anomaly + (-0.5 * e35 + @eccentricity + (e34 + 1.5 * e33 * e35) * e33) * Math.sin(mod_mean_anomaly*Math::PI/180)
end

def J2000_ecliptic_plane
    x_prime = @semi_major_axis * (Math.cos(eccentric_anomaly*Math::PI/180) - @eccentricity)
    y_prime = @semi_major_axis * Math.sqrt(1-@eccentricity**2) * Math.sin(eccentric_anomaly*Math::PI/180)
    z_prime = 0
    x = x_prime * (Math.cos(@argument_of_perhelion*Math::PI/180) * Math.cos(@longitude_of_ascending_node*Math::PI/180) - Math.sin(@argument_of_perhelion * Math::PI/180) * Math.sin(@longitude_of_ascending_node * Math::PI/180) * Math.cos(@inclination * Math::PI/180)) + y_prime * (-Math.sin(@argument_of_perhelion* Math::PI/180) * Math.cos(@longitude_of_ascending_node * Math::PI/180) - Math.cos(@argument_of_perhelion * Math::PI/180) * Math.sin(@longitude_of_ascending_node * Math::PI/180) * Math.cos(@inclination * Math::PI/180))
    y = x_prime * (Math.cos(@argument_of_perhelion*Math::PI/180) * Math.sin(@longitude_of_ascending_node*Math::PI/180) + Math.sin(@argument_of_perhelion * Math::PI/180) * Math.cos(@longitude_of_ascending_node * Math::PI/180) * Math.cos(@inclination * Math::PI/180)) + y_prime * (-Math.sin(@argument_of_perhelion* Math::PI/180) * Math.sin(@longitude_of_ascending_node * Math::PI/180) + Math.cos(@argument_of_perhelion * Math::PI/180) * Math.cos(@longitude_of_ascending_node * Math::PI/180) * Math.cos(@inclination * Math::PI/180))
    z = x_prime * Math.sin(@argument_of_perhelion*Math::PI/180) * Math.sin(@inclination*Math::PI/180) + y_prime * Math.cos(@argument_of_perhelion*Math::PI/180) * Math.sin(@inclination*Math::PI/180)
    return x, y, z
end
end

time = get_centuries_past_j2000
mars = Planet.new(1.52371034, 0.00001847, 0.09339410, 0.00007882, 1.84969142, -0.00813131, -4.553443205, 19140.30268499, -23.94362959, 0.44441088, 49.55952891, -0.29257343, time)
puts time
puts mars.mean_anomaly
puts mars.eccentric_anomaly
puts mars.J2000_ecliptic_plane

Answer 1

虽然我不同意地球近日点的论点，但这可能会有所帮助。近日点的经度很好。倾角是如此之小，以至于它并不真正适用于地球，就像它适用于其他行星一样。为 Omega 寻找值(value)具有挑战性。近日点在不断变化。仅供引用，公元 1248 年恰逢冬至。

首先，IAU 有免费的 SOFA C 和 FORTRAN 库，带有标准化的天文函数。矩阵表包含在某些例程中，因此您不必去查找它们。

但是，如果您非常倾向于使用老派方法，那么这个网站可以满足您的需求 http://www.stjarnhimlen.se/comp/tutorial.html

NOVA C 和 JAVA、MICA、JPL 目录、Jean Meeus 的书、AA USNO 和除维基百科之外的许多其他书籍都有大量信息。看起来你想要矩形值，所以我认为 Paul Schlyter 可以帮助你。

SOFA 也有这些，但是关于如何使用它们的文档不会教授这些技术。需要大量研究才能理解它们。

看起来您正在使用 Ruby，并且有一个名为 Celes 的 SOFA 库的包装器 gem。只需 gem install celes 即可。

尝试查看所有以 fa 开头的基本参数:

** iauFal03 表示月球异常** iauFaf03 表示月球纬度的参数** iauFaom03 表示月球升交点的经度** iauFame03 表示水星经度** iauFave03 是金星的经度** iauFae03 表示地球经度** iauFama03 表示火星的经度** iauFaju03 表示木星的经度** iauFasa03 土星的平均经度** iauFaur03 天王星的平均经度** iauFapa03 一般累积经度进动

玩得开心!

编辑更新:

此 gem 中的两个函数将为您提供地球的日心和重心 x、y、z。

p 是位置，v 是速度。h 是日心，b 是质心。

pvh = Celes.epv00(jd_now, 0.0)[0]
pvb = Celes.epv00(jd_now, 0.0)[1]

sc = Celes.pv2s(pvh)

sc 表示球坐标。

如您所见，您只需提供一个京东时间值即可。那个 gem 和 SOFA C 代码中有很多好东西。我还没有学会如何使用它们。