2

我在这里使用火星的方程式和数据 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
4

1 回答 1

0

尽管我不同意地球近日点的论点，但这可能会有所帮助。近日点经度很好。倾角是如此之小，以至于它并不像其他行星那样真正适用于地球。为 Omega 寻找价值具有挑战性。近日点是不断变化的。仅供参考，公元 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 表示球坐标。

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

于 2015-11-08T02:41:50.830 回答