# Keplerian Orbital Elements to Cartesian Coordinates issue with inclination, and possibly other stuff

this is the first time I am asking for help on this site, so if I format badly, or need to add/change something, let me know.

I am currently trying to make a game using python (python powered board game is the easiest way to explain it), which is about colonizing space from 1950-2100. In order to properly calculate things like travel time, I need to calculate distance, and so I need xyz position, and so I need to calculate the position of the object, which requires the Keplerian orbital elements. I have gone through this site, and although it seems others have been satisfied with their answers, my implementation of those answers has resulted in some clearly wrong orbits. The best way to explain it to to attach an image: Red = Mars, Blue = Earth, Green = Venus

I have used the numbers from Wikipedia, made them all into years, AU, and Radians, but for some reason only Earth's orbit is fine. I have used two different translation formulas from this website, and they both give almost exactly the same thing, so I don't feel like its that, but just in case I will be attaching my full code.

full_kep2cart.py:

import math

## EDIT - Orbital Period is from here: https://phas.ubc.ca/~newhouse/p210/orbits/cometreport.pdf (page 7, Values into P) ; I can't currently find the Mean Anomaly equation ##
def f_Mean_Anomaly(TimeSincePerihelion, Semi_Major_Axis):
OrbitalPeriod = math.sqrt(Semi_Major_Axis**3)
MeanAnomaly = (2 * math.pi * TimeSincePerihelion) / OrbitalPeriod
return MeanAnomaly

def fApoapsis(Semi_Major_Axis, Eccentricity):
Answer = Semi_Major_Axis * (1 + Eccentricity)

def fPeriapsis (Semi_Major_Axis, Eccentricity):
Answer = Semi_Major_Axis * (1 - Eccentricity)

## EDIT - From here: https://space.stackexchange.com/a/59600/48051 ##
def TrueAnomaly(MeanAnomaly, eccentricity) :
# Place M in [0,2*pi)
MeanAnomaly -= math.floor(MeanAnomaly/(2*math.pi))*(2*math.pi)

# Solve Kepler's equation for eccentric anomaly
if eccentricity > 0.5 :
Eccentric_Anomaly = math.pi
else :
Eccentric_Anomaly = MeanAnomaly
while True :
delta = (MeanAnomaly - (Eccentric_Anomaly - eccentricity*math.sin(Eccentric_Anomaly))) / (1.0 - eccentricity*math.cos(Eccentric_Anomaly))
Eccentric_Anomaly += delta
if abs(delta) <= 1e-15*Eccentric_Anomaly : break

# Solve for true anomaly given the eccentric anomaly
True_Anomaly = 2.0*math.atan(math.sqrt((1.0-eccentricity)/(1.0+eccentricity))*math.tan(Eccentric_Anomaly/2))
if True_Anomaly < 0.0 : True_Anomaly += (2*math.pi)

return True_Anomaly

## EDIT - Also from here: https://space.stackexchange.com/a/59600/48051 ##
def EccentricAnomaly(MeanAnomaly, eccentricity) :
# Place M in [0,2*pi)
MeanAnomaly -= math.floor(MeanAnomaly/(2*math.pi))*(2*math.pi)

# Solve Kepler's equation for eccentric anomaly
if eccentricity > 0.5 :
Eccentric_Anomaly = math.pi
else :
Eccentric_Anomaly = MeanAnomaly
while True :
delta = (MeanAnomaly - (Eccentric_Anomaly - eccentricity*math.sin(Eccentric_Anomaly))) / (1.0 - eccentricity*math.cos(Eccentric_Anomaly))
Eccentric_Anomaly += delta
if abs(delta) <= 1e-15*Eccentric_Anomaly : break

# Solve for true anomaly given the eccentric anomaly
nu = 2.0*math.atan(math.sqrt((1.0+eccentricity)/(1.0-eccentricity))*math.tan(0.5*Eccentric_Anomaly))
if nu < 0.0 : nu += (2*math.pi)

return Eccentric_Anomaly

## EDIT - I switched this and the next one around, so that the second implementation was after the first. This first one was taken from this google doc: https://web.archive.org/web/20170810015111/http://ccar.colorado.edu/asen5070/handouts/kep2cart_2002.doc ; Radius equation is from https://phas.ubc.ca/~newhouse/p210/orbits/cometreport.pdf page 10 Solving for r ##
def Keplar2Cartesian(True_Anomaly, Semi_Major_Axis, Eccentricity, Longitude_Of_TheAscending_Node, Argument_Periapsis, Inclination):
φ = True_Anomaly
a = Semi_Major_Axis
e = Eccentricity
Ω = Longitude_Of_TheAscending_Node
ω = Argument_Periapsis
i = Inclination

radius = a * ((1 - e**2)/(1 + e * math.cos(φ)))

X_Value = r * (math.cos(Ω) * math.cos(ω+φ) - math.sin(Ω) * math.sin(ω+φ) * math.cos(i))
Y_Value = r * (math.sin(Ω) * math.cos(ω+φ) + math.cos(Ω) * math.sin(ω+φ) * math.cos(i))
Z_Value = r * (math.sin(ω+φ) * math.sin(i))
return X_Value, Y_Value, Z_Value

## EDIT - I should first note that this and the equation above give very similar results. I basically just did a translation from whatever language its in to python: https://space.stackexchange.com/a/8915/48051
def Keplar2CartesianSecond(ArgPeriapsis, LongofAscendingNode, EccentricAnomaly, SemiMajorAxis, Eccentricity, Inclination):
W = LongofAscendingNode
E = EccentricAnomaly
a = SemiMajorAxis
e = Eccentricity
i = Inclination

p = LongofAscendingNode + ArgPeriapsis

w = p - W

P = a * (math.cos(E) - e)
Q = a * math.sin(E) * math.sqrt(1 - e **2)

# rotate by argument of periapsis
x = math.cos(w) * P - math.sin(w) * Q
y = math.sin(w) * P + math.cos(w) * Q
# rotate by inclination
z = math.sin(i) * y
y = math.cos(i) * y
# rotate by longitude of ascending node
xtemp = x
x = math.cos(W) * xtemp - math.sin(W) * y
y = math.sin(W) * xtemp + math.cos(W) * y

return x, y, z



Kep2Cart_SeriesOfOutputs.py:

import full_kep2cart as k2c
import matplotlib.pyplot as plt

#################################################################################
#Input:                   #######################################################
#################################################################################
Semi_Major_Axis = 0.723332
Eccentricity_In = 0.006772
LongitudeOfTheAscendingNode = 1.3383184704292521
ArgumentPeriapsis = 0.9579065066645678
Inclination_In = 0.05924659772234911

Time_Change_Factor = 12

#how many pixels per input unit
Scale_Factor = 50

TimesToRepeat = 2000

X_Values = []
Y_Values = []
Z_Values = []

for x in range (1, (TimesToRepeat + 1)):

Time_Since_Perihelion = x * (1/Time_Change_Factor)            #use same unit as Orbital Period

Mean_Anomaly = k2c.f_Mean_Anomaly(Time_Since_Perihelion, Semi_Major_Axis)
TrueAnomaly_In = k2c.TrueAnomaly(Mean_Anomaly, Eccentricity_In)
EccentricAnomaly_In = k2c.EccentricAnomaly(Mean_Anomaly, Eccentricity_In)

X, Y, Z = k2c.Keplar2Cartesian(TrueAnomaly_In, Semi_Major_Axis, Eccentricity_In, LongitudeOfTheAscendingNode, ArgumentPeriapsis, Inclination_In)

X_Values.append(X)
Y_Values.append(Y)
Z_Values.append(Z)

############################################

fig = plt.figure()

ax = plt.axes(projection="3d")

ax.plot3D(X_Values, Y_Values, Z_Values, color = 'green')

plt.show()


I am not very good with this stuff, so I'm not sure what to do next. If someone can figure out whats wrong, or even offer up an alternative method for estimating travel time, that would be great.

• Please clearly state what your question is. Is it "please debug my code?" Aug 10, 2022 at 11:24
• Also consider adding links to the pages you got the equations from. Aug 10, 2022 at 12:42
• @OrganicMarble I'm not entirely sure what I need answered. I guess its A: is there an easier way to do this (since it doesn't need to be %100 correct), or B: does anyone know if my code is bad, or if something else is the issue? Aug 10, 2022 at 20:33
• @GregMiller I've tried with at least these: space.stackexchange.com/a/19335/48051 , phas.ubc.ca/~newhouse/p210/orbits/cometreport.pdf , space.stackexchange.com/a/8915/48051 , space.stackexchange.com/a/59600/48051 . I'll go through and edit in some of these into the appropriate places for better readability, however overall its hard to say what comes from where since I've looked at a lot of different stuff and I probably have small bits and pieces of stuff all over the place. Aug 10, 2022 at 21:03
• The first broken thing is "OrbitalPeriod = math.sqrt(Semi_Major_Axis**3)". Orbital period has units of time, and semimajor axis has units of length, so you need a conversion factor with the proper units. Some research will tell you that the factor in question depends on the masses of the two objects (though the smaller can safely be ignored for artificial satellites, it can be too big to ignore for natural planets and moons). Aug 10, 2022 at 21:22