Calculating the x,y,z position of solar bodies using ruby

Question

I am using the equations and data for Mars here http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf

and the solution for the eccentric anomaly kepler equation given here at the top of page four http://murison.alpheratz.net/dynamics/twobody/KeplerIterations_summary.pdf

And checking the output by modifying the date in the get_centuries_past to the following dates and looking at page E-7 for the actual x,y,z coordinates of Mars (sample data below, but link for the curious: http://books.google.com/books/about/Astronomical_Almanac_for_the_Year_2013_a.html?id=7fl_-DLwJ8YC)

date 2456320.5 is 2013, 1, 28 and should output

x = 1.283762
y = -0.450111
z = -0.241123

date 2456357.5 is 2013, 3, 6 and should output

x = 1.300366
y = 0.533593
z = 0.209626

date 2456539.500000 is 2013, 9, 4 and should output

x = - 0.325604
y = 1.418110
z = 0.659236

I tested mean anomaly equation and it was fine. However, I cannot get a good set of x,y,z coordinates. I have been tweaking my kepler and coordinate function but cannot get them to match the tables in the astronomical almanac.

Any suggestions or advice on solving the positions of the stars is greatly appreciated. The code below can be put in a .rb file and running it on the command line will output the x,y,z values.

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

Douglas G. Allen · Answer 1 · 2015-11-20T20:50:52.677

This may be helpful although I don't agree with arguments of perihelion for Earth. Longitude of perihelion is fine. The inclination is so small that it doesn't really apply for Earth as it does for other planets. And finding values for Omega is challenging. The perihelion is ever changing. FYI in 1248 AD it coincided with the winter solstice.

Firstly IAU has free SOFA C and FORTRAN lib's with standardized astronomical functions. Matric tables are included in certain routines so you don't have to go look them up.

But if you are so inclined to use old school methods then this site has what you need http://www.stjarnhimlen.se/comp/tutorial.html

NOVA C and JAVA, MICA, JPL catalogs, Jean Meeus book, AA USNO and a host of others beside wikipedia have tons of information. It looks like you want rectangular values so I think Paul Schlyter can help you out.

SOFA has these as well but documentation on how to use them is not going to teach the techniques. It will take a lot of research to understand them.

It looks like you are using Ruby and there is a wrapper gem for the SOFA lib named Celes. Just gem install celes and you'll have it.

Try looking at the fundamental arguments which all begin with fa:

** iauFal03 mean anomaly of the Moon ** iauFaf03 mean argument of the latitude of the Moon ** iauFaom03 mean longitude of the Moon's ascending node ** iauFame03 mean longitude of Mercury ** iauFave03 mean longitude of Venus ** iauFae03 mean longitude of Earth ** iauFama03 mean longitude of Mars ** iauFaju03 mean longitude of Jupiter ** iauFasa03 mean longitude of Saturn ** iauFaur03 mean longitude of Uranus ** iauFapa03 general accumulated precession in longitude

Have fun!

Edit Update:

Two functions in this gem that will give you heliocentric and barycentric x,y,z for Earth.

p is position and v is velocity. h is heliocentric and b is barycentric.

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

sc = Celes.pv2s(pvh)

sc means spherical coordinates.

As you can see, all you need to provide is a JD time value. So lots of good stuff in that gem and the SOFA C code. I have yet to learn how to use them all.

Calculating the x,y,z position of solar bodies using ruby

1 Answers1