# Written by David Sims in Python 3.4.3, released to the public domain.
# Download Python 3.4.3 from
https://www.python.org/downloads/release/python-343/
# This program is a modification of ephem, intended to track the Juno spacecraft in orbit around Jupiter.
import math
AU = 1.495978707e11
pi = 3.1415926535897932384626433832795
dr = pi/180.0
GM = 1.26686534e17
toi = 2457627.9966
# Update Juno's orbital elements relative to Jupiter Body Center (500@599)
# from JPL Horizons (
http://ssd.jpl.nasa.gov/horizons.cgi)
# sma = semimajor axis, AU
sma = 0.02740009798134841
# ecc = eccentricity
ecc = 0.9816687162279770
# tpp = time of perjove passage, Julian date
tpp = 2457628.036427238025
sma = sma*AU
# The results will be vectors in the orbit's canonical coordinate system.
# There will be no conversion to ecliptic coordinates.
P = (pi/43200)*math.sqrt(sma**3/GM)
m0 = (toi-tpp)/P
m = 2*pi*(m0-int(m0))
u = m + (ecc-ecc**3/8+ecc**5/192)*math.sin(m)
u = u + (ecc*ecc/2-ecc**4/6)*math.sin(2*m)
u = u + (3*ecc**3/8-27*ecc**5/128)*math.sin(3*m)
u = u + (ecc**4/3)*math.sin(4*m)
U = 999.9
# Replace the four underscores with four spaces, where necessary.
while abs(u-U)>1.0e-14:
____U = u
____F0 = U-ecc*math.sin(U)-m
____F1 = 1-ecc*math.cos(U)
____F2 = ecc*math.sin(U)
____F3 = ecc*math.cos(U)
____D1 = -F0/F1
____D2 = -F0/(F1+D1*F2/2)
____D3 = -F0/(F1+D1*F2/2+D2*D2*F3/6)
____u = U+D3
if u<0:
____u = u+2*pi
x = sma*(math.cos(u)-ecc)
y = sma*math.sin(u)*math.sqrt(1-ecc*ecc)
r = math.sqrt(x*x+y*y)
q = math.atan(y/x)
if x<0:
____q=q+pi
if x>0 and y<0:
____q=q+2.0*pi
k = math.sqrt(GM/(sma*(1-ecc*ecc)))
Vx = -k*math.sin(q)
Vy = k*(ecc+math.cos(q))
V = math.sqrt(Vx*Vx+Vy*Vy)
print('toi {:15.7f}'.format(toi),'JD')
print('x {:15.3f}'.format(x),'meters')
print('y {:15.3f}'.format(y),'meters')
print('r {:15.3f}'.format(r),'meters')
print('Vx {:15.9f}'.format(Vx),'m/s')
print('Vy {:15.9f}'.format(Vy),'m/s')
print('V {:15.9f}'.format(V),'m/s')
print('T.A. {:15.11f}'.format(q/dr),'degrees')
print('E.A. {:15.11f}'.format(u/dr),'degrees')
print('P {:15.11f}'.format(P),'days')
keypress = input('Press return to exit program.')
