# Orbital elements of solar system using Pyephem

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 ?

• You only have a position vector at the moment. To obtain the full set of Keplerian elements you also need a velocity vector. Commented Nov 3, 2016 at 10:30
• 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 Commented Nov 3, 2016 at 10:37
• 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.
– user21
Commented Nov 3, 2016 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):
'''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)

• 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
Commented Nov 11, 2016 at 0:55