Skip to content

Instantly share code, notes, and snippets.

@croxis
Last active December 21, 2015 00:48
Show Gist options
  • Save croxis/6222986 to your computer and use it in GitHub Desktop.
Save croxis/6222986 to your computer and use it in GitHub Desktop.
Compute the x,y,z position of a body relative to its parent, using precalculated orbital parameters, as a function of time
from math import sin, cos, radians, degrees, sqrt, atan2
#Imported from a json/yaml/whatever config file
#This is Earth's orbit. Time is number of Julian2000 days
planet_config = {
"a": 1.00000018,
"e": eval("lambda d: 0.01673163 - 3.661e-07 * d"),
"w": eval("lambda d: 108.04266274 + 0.0031795260 * d"),
"i": eval("lambda d: -0.00054346 + -0.0001337178 * d"),
"M": eval("lambda d: -2.4631431299999917"),
"N": eval("lambda d: -5.11260389 + -0.0024123856 * d")
}
# Time is represented in days
time = 0.0
def compute_position(config, time):
'''Returns x, y, z coordinates of a body relative to it's parent.
The units returned are the units used in the orbital parameters.
If using nasa/jpl data computing a planet will return units in AU, moons in km.
Scale as needed. Plenty of optimization can be done still for production code.'''
M = radians(config['M'](time))
w = radians(config['w'](time))
i = radians(config['i'](time))
N = radians(config['N'](time))
a = config['a']
e = config['e'](time)
# Compute eccentric anomaly
E = M + e * sin(M) * (1.0 + e * cos(M))
if degrees(E) > 0.05:
E = computeE(E, M, e)
# http://stjarnhimlen.se/comp/tutorial.html
# Compute distance and true anomaly
xv = a * (cos(E) - e)
yv = a * (sqrt(1.0 - e * e) * sin(E))
v = atan2(yv, xv)
r = sqrt(xv * xv + yv * yv)
xh = r * (cos(N) * cos(v + w) - sin(N) * sin(v + w) * cos(i))
yh = r * (sin(N) * cos(v + w) + cos(N) * sin(v + w) * cos(i))
zh = r * (sin(v + w) * sin(i))
position = (xh, yh, zh)
return position
def computeE(E0, M, e):
'''Iterative function for a higher accuracy of E'''
E1 = E0 - (E0 - e * sin(E0) - M) / (1 - e * cos(E0))
if abs(abs(degrees(E1)) - abs(degrees(E0))) > 0.001:
E1 = self.computeE(E1, M, e)
return E1
if __name__ == '__main__':
'''For the sake of the demo time will be sped up'''
while True:
print('{:10.4f}'.format(time), compute_position(planet_config, time))
time += 0.01
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment