I'm using the astrometric python module Pyephem, and I would like to get the orbital (keplerian) elements for the solar system planets.

The only values I found is the heliocentric latitude, longitude, and the distance to the sun. Is there a way to compte the orbital parameters based on thoses values ? Did I missed a function in Pyephem ?

  • 1
    $\begingroup$ You only have a position vector at the moment. To obtain the full set of Keplerian elements you also need a velocity vector. $\endgroup$ – SE - stop firing the good guys Nov 3 '16 at 10:30
  • 2
    $\begingroup$ While possible, to use the position and velocity to calculate orbital elements, this is working backwards. pyephem uses the Keplerian elements to calculate the position. The orbital elements are fixed, and can be read from the Horizons site $\endgroup$ – James K Nov 3 '16 at 10:37
  • 2
    $\begingroup$ You might try looking at skyfield, but it appears even that doesn't have orbital elements directly. I know the ephemeris files don't have orbital elements either, but you might try looking at VSOP 2013. The orbital elements are fairly constant, but do change slowly due to perturbations. $\endgroup$ – user21 Nov 3 '16 at 13:24

There is a website: http://orbitsimulator.com/formulas/OrbitalElements.html that has a javascript program for converting state vectors to orbital elements and back.

The source of that website is possible to convert to python: Here is a Body class that has a method for calculating orbital elements based on that website. The method takes one argument, a body the is the Principal (ie the Sun)

G = 0.0002946        # in units of seconds, AU and solar masses.
class Body:
    """A body has attributes r and v, which are its position and 
velocity in cartesian coordinates and a mass. implied units are 
solar masses, AU and seconds.""" 
    def __init__(self,r,v,mass):
        self.r = np.array(r,dtype="float")
        self.v = np.array(v,dtype="float")
        self.mass = mass
        self.GM = self.mass*G

    def orbital_elements(self,principal):
        mu = G*(principal.mass+self.mass)
        # calculate relative position,velocity
        r = self.r - principal.r
        v = self.v - principal.v
        try: #catch division by zero
            R = np.linalg.norm(r)
            V = np.linalg.norm(v)
            a = 1/(2/R - V**2/mu)  # semi major axis

            h = np.cross(r,v)
            H = np.linalg.norm(h)

            P = H**2/mu
            Q = np.dot(r,v)

            E= np.sqrt(1-P/a)  #eccentricity

            e = [1-R/a,Q/np.sqrt(a*mu)]
            i = np.arccos(h[2]/H)
            Omega = 0
            if i!=0: 
                Omega = np.arctan2(h[0],-h[1]) #Longitude of acending node

            ta = [H**2/(R*mu) -1,H*Q/(R*mu)]
            TA = np.arctan2(ta[1],ta[0])
            Cw = (r[0]*np.cos(Omega)+r[1]*np.sin(Omega))/R

            if i==0 or i==np.pi:
                Sw = (r[1]*np.cos(Omega) - r[0]*np.sin(Omega))/R
                Sw = r[2]/(R*np.sin(i))

            omega = np.arctan2(Sw,Cw) - TA  #argument of periapsis
            if omega<0: omega += 2*np.pi

            u = np.arctan2(e[1],e[0]) 
            M = u - E*np.sin(u) # mean anomaly

        except ZeroDivisionError:
            #meaningless, but avoids crash
| improve this answer | |
  • 1
    $\begingroup$ I've been enjoying this answer - but I'm confused by the reduced mass. In this formulation, a given input state vector $(\vec{r}, \vec{v})$ gives different answers depending on both masses. With respect to what coordinate frame is the state vector assumed to be measured - the frame of the principal or that of the center of mass? To which orbit do the returned elements apply - the orbit of the Body in the frame of the principal, or the center of mass, or is it the one-body orbit in the central field? $\endgroup$ – uhoh Nov 11 '16 at 0:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.