punyverse/punyverse/orbit.py
Quantum a862302408 Use the faster Newton's method for eccentric anomaly.
Newton's method also converges quickly for highly eccentric orbits.
2018-08-23 21:49:38 -04:00

97 lines
3.1 KiB
Python

from math import sin, cos, tan, atan, sqrt, radians, degrees
class KeplerOrbit(object):
def __init__(self, sma, eccentricity, inclination=0, longitude=0, argument=0):
self.sma = sma
self.eccentricity = eccentricity
self.inclination = inclination
self.longitude = longitude
self.argument = argument
self.__true_anomaly_factor = sqrt((1 + eccentricity)/(1 - eccentricity))
self.__distance_factor = sma * (1 - eccentricity ** 2)
@property
def inclination(self):
return degrees(self._inclination)
@inclination.setter
def inclination(self, value):
self.inclination_radian = radians(value)
@property
def inclination_radian(self):
return self._inclination
@inclination_radian.setter
def inclination_radian(self, value):
self._inclination = value
self.__sin_inclination = sin(self._inclination)
self.__cos_inclination = cos(self._inclination)
@property
def longitude(self):
return degrees(self._longitude)
@longitude.setter
def longitude(self, value):
self.longitude_radian = radians(value)
@property
def longitude_radian(self):
return self._longitude
@longitude_radian.setter
def longitude_radian(self, value):
self._longitude = value
self.__sin_longitude = sin(self._longitude)
self.__cos_longitude = cos(self._longitude)
@property
def argument(self):
return degrees(self._argument)
@argument.setter
def argument(self, value):
self.argument_radian = radians(value)
@property
def argument_radian(self):
return self._argument
@argument_radian.setter
def argument_radian(self, value):
self._argument = value
self.__sin_argument = sin(self._argument)
self.__cos_argument = cos(self._argument)
def eccentric_anomaly(self, mean_anomaly):
e1 = 0
e2 = mean_anomaly
while abs(e1 - e2) > 0.000001:
e1, e2 = e2, e2 - ((e2 - mean_anomaly - self.eccentricity * sin(e2)) /
(1 - self.eccentricity * cos(e2)))
return e2
def true_anomaly(self, mean_anomaly):
eccentric_anomaly = self.eccentric_anomaly(mean_anomaly)
return 2 * atan(self.__true_anomaly_factor * tan(eccentric_anomaly / 2))
def orbit(self, mean_anomaly):
mean_anomaly = radians(mean_anomaly)
phi = self.true_anomaly(mean_anomaly)
r = self.__distance_factor / (1 + self.eccentricity * cos(phi))
x = r * cos(phi)
y = r * sin(phi)
z = 0
# phi = longitude, theta = inclination, psi = argument
x, y = (x * self.__cos_longitude + y * self.__sin_longitude,
-x * self.__sin_longitude + y * self.__cos_longitude)
x, z = (x * self.__cos_inclination + z * self.__sin_inclination,
-x * self.__sin_inclination + z * self.__cos_inclination)
x, y = (x * self.__cos_argument + y * self.__sin_argument,
-x * self.__sin_argument + y * self.__cos_argument)
return x, y, z