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)