Estoy usando el módulo de pitón astrométrico Pyephem, y me gustaría obtener los elementos orbitales (keplerianos) para los planetas del sistema solar.
Los únicos valores que encontré son la latitud heliocéntrica, la longitud y la distancia al sol. ¿Hay alguna forma de calcular los parámetros orbitales en función de esos valores? ¿Me perdí una función en Pyephem?
Hay un sitio web: http://orbitsimulator.com/formulas/OrbitalElements.html que tiene un programa javascript para convertir vectores de estado en elementos orbitales y viceversa.
La fuente de ese sitio web se puede convertir a python: Aquí hay una Body
clase que tiene un método para calcular elementos orbitales basado en ese sitio web. El método toma un argumento, un cuerpo que es el Principal (es decir, el Sol)
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)
SE - deja de despedir a los buenos
james k
usuario21