Remove old Keplerian library

With thanks to its author, @hoylemd, for this early prototype!
This commit is contained in:
Brandon Rhodes 2020-07-09 19:31:46 -04:00
parent 4fb5e7a9f4
commit 43bb53036d
2 changed files with 0 additions and 251 deletions

View File

@ -1,146 +0,0 @@
"""A work in progress, that does not yet agree with JPL Horizons.
See https://github.com/brandon-rhodes/python-skyfield/issues/11
"""
from math import sin, cos
import math
from skyfield import constants
# class to represent a point in the IC reference frame
class ICRCoordinates:
def __init__(self, x, y, z):
self.x = x
self.y = y
self.z = z
def equalTo(self, other):
# TODO: override ==, and add epsilons here
return (self.x == other.x) and (self.y == other.y) and (self.z == other.z)
def __repr__(self):
return '(%s, %s, %s)' % (self.x, self.y, self.z)
def semimajorAxisToOrbitalPeriod(axis):
return (axis ** 3) ** 0.5
def orbitalPeriodToSemimajorAxis(period):
return (period ** 2) ** (1.0 / 3.0)
def convergeEccentricAnomaly(mean_anomaly, eccentricity, precision):
# calculate the delta
delta = 10 ** -precision
# normalize the mean anomaly
m = mean_anomaly % constants.tau
# set up the first guess
eccentric_anomaly = constants.tau
if eccentricity < 0.8:
eccentric_anomaly = m
# do the initial test
test = eccentric_anomaly - eccentricity * sin(m) - m
# while we're not happy with the result, and we haven't been dawdling too long
max_iterations = 30
count = 0
while ((math.fabs(test) > delta) and (count < max_iterations)):
# calculate the next guess for an eccentric anomaly
eccentric_anomaly = (
eccentric_anomaly - test /
(1.0 - eccentricity * cos(eccentric_anomaly))
)
# try it
test = eccentric_anomaly - eccentricity * sin(eccentric_anomaly) - m
# count the runs, so we don't go forever
count += 1
# convert to degrees
return eccentric_anomaly
def calculateMeanAnomaly(L, wb):
return L - wb
class KeplerianOrbit:
def __init__(
self,
semimajor_axis,
eccentricity,
inclination,
longitude_ascending,
argument_perihelion,
mean_anomaly,
epoch
):
self.semimajor_axis = semimajor_axis
self.eccentricity = eccentricity
self.inclination = inclination
self.longitude_ascending = longitude_ascending
self.argument_perihelion = argument_perihelion
self.mean_anomaly = mean_anomaly
self.epoch = epoch
def getECLCoordinatesAtTime(self, date):
# localize the orbital parameters
a = self.semimajor_axis
e = self.eccentricity
I = self.inclination
Om = self.longitude_ascending
#n = 0.230605479
n = 0.230652907
w = self.argument_perihelion
M = self.mean_anomaly
d = date.tdb - self.epoch.tdb
Om = Om / 360.0 * constants.tau
w = w / 360.0 * constants.tau
I = I / 360.0 * constants.tau
M = M / 360.0 * constants.tau
n = n / 360.0 * constants.tau
M += d * n
# calculate the mean anomaly in rads
E = convergeEccentricAnomaly(M, e, 30)
# calculate the initial primes
x_prime = a * (cos(E) - e)
y_prime = a * (1 - e ** 2.0) ** (0.5) * sin(E)
"""
http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf
x_ecl = cos(w)cos(Om)-sin(w)sin(Om)cos(I) * x_prime +
(-sin(w)cos(Om) - cos(w)sin(Om)cos(I)) * y_prime
y_ecl = cos(w)sin(Om)-sin(w)cos(Om)cos(I) * x_prime +
(-sin(w)cos(Om) - cos(w)sin(Om)cos(I)) * y_prime
z_ecl = (sin(w)sin(I)) * x_prime +
(cos(w)sin(I)) * y_prime
"""
# calculate the ecliptic coordinates
x_ecl = ((cos(w) * cos(Om) - sin(w) * sin(Om) * cos(I)) * x_prime +
(-1 * sin(w) * cos(Om) - cos(w) * sin(Om) * cos(I)) * y_prime)
y_ecl = ((cos(w) * sin(Om) + sin(w) * cos(Om) * cos(I)) * x_prime +
(-1 * sin(w) * sin(Om) + cos(w) * cos(Om) * cos(I)) * y_prime)
z_ecl = ((sin(w) * sin(I)) * x_prime + (cos(w) * sin(I)) * y_prime)
return ICRCoordinates(x_ecl, y_ecl, z_ecl)
def getICRSCoordinatesAtTime(self, date):
# J2000 obliquity
e = 23.43928 * math.pi / 180.0
# get the ecliptic coords
ecliptic = self.getECLCoordinatesonTime(date);
# calculate the equatorial (ICRF) coordinates
x_eq = ecliptic.x;
y_eq = cos(e) * ecliptic.y - sin(e) * ecliptic.z
z_eq = sin(e) * ecliptic.y + cos(e) * ecliptic.z
x_eq, y_eq, z_eq

View File

@ -1,105 +0,0 @@
"""Compare the output of Skyfield with the routines from NOVAS for keplerian orbiting bodies"""
import keplerianlib
from skyfield import api
from keplerianlib import KeplerianOrbit, ICRCoordinates
DISTANCE_EPSILON = 0.026
def test_semimajorAxisToOrbitalPeriod():
assert keplerianlib.semimajorAxisToOrbitalPeriod(1) == 1
assert keplerianlib.semimajorAxisToOrbitalPeriod(1.523679) == 1.8807896358663763
assert keplerianlib.semimajorAxisToOrbitalPeriod(4.27371348392) == 8.835031547398543
def test_orbitalPeriodToSemimajorAxis():
assert keplerianlib.orbitalPeriodToSemimajorAxis(1) == 1
assert keplerianlib.orbitalPeriodToSemimajorAxis(1.8807896358663763) == 1.523679
assert keplerianlib.orbitalPeriodToSemimajorAxis(8.835031547398543) == 4.27371348392
def test_convergeEccentricAnomaly():
test = keplerianlib.convergeEccentricAnomaly(
hoyle_8077['mean_anomaly'],
hoyle_8077['eccentricity'],
15
)
assert test == hoyle_8077['eccentric_anomaly']
def test_instantiate_8077_hoyle():
ts = api.load.timescale()
hoyle = KeplerianOrbit( hoyle_8077['semimajor_axis'],
hoyle_8077['eccentricity'],
hoyle_8077['inclination'],
hoyle_8077['longitude_ascending'],
hoyle_8077['argument_perihelion'],
hoyle_8077['mean_anomaly'],
ts.tt(jd=hoyle_8077['epoch_tt']))
assert hoyle is not None
def test_instantiate_coordinates():
coords = ICRCoordinates(x=500.25, y=10.76, z=0.1125)
assert coords is not None
def test_coordinatesEquivalence():
coords_the_first = ICRCoordinates(x=500.25, y=10.76, z=0.1125)
coords_the_second = ICRCoordinates(x=500.25, y=10.76, z=0.1125)
assert coords_the_first.equalTo(coords_the_second)
"""
data gotten from Horizons
date: 2456517.500000000 = A.D. 2013-Aug-13 00:00:00.0000 (CT)
expected coords (AU) 2.421251132790093E+00 -1.918893156489506E+00 -9.813409585464707E-02
Horizon Params
Ephemeris Type [change] : VECTORS
Target Body [change] : Asteroid 8077 Hoyle (1986 AW2)
Coordinate Origin [change] : Solar System Barycenter (SSB) [500@0]
Time Span [change] : Start=2013-08-13, Stop=2013-09-12, Step=1 d
Table Settings [change] : defaults
Display/Output [change] : default (formatted HTML)
EPOCH= 2453995.5 ! 2006-Sep-17.00 (CT) Residual RMS= .43359
EC= .2110946491840378 QR= 2.077692130214496 TP= 2454360.1855338747
OM= 135.855972529608 W= 34.4378477722205 IN= 17.25814783060462
A= 2.633639292806857 MA= 275.9015153135912 ADIST= 3.189586455399217
PER= 4.27408 N= .230605479 ANGMOM= .027287332
DAN= 2.14316 DDN= 3.04671 L= 169.07331
B= 9.658454799999999 MOID= 1.10581994 TP= 2007-Sep-16.6855338747
"""
hoyle_8077 = {
'semimajor_axis' : 2.633278254269645,
'eccentricity' : .2109947010748546,
'inclination' : 17.25945395594321,
'longitude_ascending' : 135.8512354853258,
'argument_perihelion' : 34.46503170092878,
'mean_anomaly' : 330.9918926661418,
'eccentric_anomaly' : 4.0942988262501965,
'epoch_tt' : (2007, 5, 14),
}
def test_get_8077_hoyle_ecliptic_on_dev_sprint_day_2():
ts = api.load.timescale()
hoyle = KeplerianOrbit( hoyle_8077['semimajor_axis'],
hoyle_8077['eccentricity'],
hoyle_8077['inclination'],
hoyle_8077['longitude_ascending'],
hoyle_8077['argument_perihelion'],
hoyle_8077['mean_anomaly'],
ts.tt(*hoyle_8077['epoch_tt']))
date = ts.tt(2013, 8, 13)
# print date.tt
test = hoyle.getECLCoordinatesAtTime(date)
#print test
epsilon = 2e-2
assert abs(test.x - 2.421251271197979) < epsilon
assert abs(test.y - -1.918893007049262) < epsilon
assert abs(test.z - -0.09813403009731327) < epsilon