Я использую астрометрический модуль Python Pyephem и хочу получить орбитальные (кеплеровские) элементы для планет Солнечной системы.
Единственные значения, которые я нашел, это гелиоцентрическая широта, долгота и расстояние до солнца. Есть ли способ вычислить параметры орбиты на основе этих значений? Я пропустил функцию в Pyephem?
Существует веб-сайт: http://orbitsimulator.com/formulas/OrbitalElements.html , на котором есть программа javascript для преобразования векторов состояния в элементы орбиты и обратно.
Источник этого веб-сайта можно преобразовать в python: вот Body
класс, в котором есть метод расчета элементов орбиты на основе этого веб-сайта. Метод принимает один аргумент, телом является Руководитель (т.е. Солнце)
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 - хватит стрелять в хороших парней
Джеймс К.
пользователь21