<p>There is a website: <a href="http://orbitsimulator.com/formulas/OrbitalElements.html" rel="nofollow noreferrer">http://orbitsimulator.com/formulas/OrbitalElements.html</a> that has a javascript program for converting state vectors to orbital elements and back.</p>

<p>The source of that website is possible to convert to python: Here is a <code>Body</code> 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)</p>

<pre><code>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&lt;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)
</code></pre>

CC BY-SA 3.0
2016-11-11T00:55:00.747Z
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?
uhoh
CC BY-SA 3.0
2016-11-04T15:24:45.600Z
2016-11-04T15:31:44.713Z
2016-11-04T15:31:44.713Z
James K
<p>I'm using the astrometric python module Pyephem, and I would like to get the orbital (keplerian) elements for the  solar system planets. </p>

<p>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 ? </p>

CC BY-SA 3.0
2016-11-03T10:30:48.913Z
You only have a position vector at the moment. To obtain the full set of Keplerian elements you also need a velocity vector.
SE - stop firing the good guys
CC BY-SA 3.0
2016-11-03T10:37:40.900Z
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
James K
<planet><solar-system><orbital-elements><python>
Orbital elements of solar system using Pyephem
Covich
# Orbital elements of solar system using Pyephem

(Asked by: username_0 on 2016-11-03)

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 ?

**Tags:** planet, solar-system, orbital-elements, python

> **Comments:**
> **username_1**: You only have a position vector at the moment. To obtain the full set of Keplerian elements you also need a velocity vector.  
> **username_2**: 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  

---

## 1 answers

### Answer by username_2 on 2016-11-04. Score: 3

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)

> **Comments:**
> **username_3**: 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?  

---

