Elementos orbitales del sistema solar usando Pyephem

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?

Solo tienes un vector de posición en este momento. Para obtener el conjunto completo de elementos keplerianos, también necesita un vector de velocidad.
Si bien es posible usar la posición y la velocidad para calcular los elementos orbitales, esto funciona al revés. pyephem usa los elementos Keplerianos para calcular la posición. Los elementos orbitales son fijos y se pueden leer desde el sitio de Horizons
Puede intentar mirar Skyfield, pero parece que incluso eso no tiene elementos orbitales directamente. Sé que los archivos de efemérides tampoco tienen elementos orbitales, pero puede intentar mirar VSOP 2013. Los elementos orbitales son bastante constantes, pero cambian lentamente debido a las perturbaciones.

Respuestas (1)

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 Bodyclase 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)
He estado disfrutando de esta respuesta, pero estoy confundido por la masa reducida. En esta formulación, un vector de estado de entrada dado ( r , v ) da diferentes respuestas dependiendo de ambas masas. ¿Con respecto a qué marco de coordenadas se supone que se mide el vector de estado: el marco del principal o el del centro de masa? ¿A qué órbita se aplican los elementos devueltos: la órbita del cuerpo en el marco del principal, o el centro de masa, o es la órbita de un cuerpo en el campo central?