Compute the x,y,z position of a body relative to its parent, using precalculated orbital parameters, as a function of time

orbit.py

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