Orbital elements of solar system using Pyephem

Question

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 ?

Answer 1

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):
        '''view-source:http://orbitsimulator.com/formulas/OrbitalElements.html'''
        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
            else:
                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

            return(a,E,omega,Omega,i,M)
        except ZeroDivisionError:
            #meaningless, but avoids crash
            return(0,0,0,0,0,0)