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

Question

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)
    return Answer

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


## 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(φ)))
    r = radius

    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

#################################################################################
#Only works with Radians  #######################################################
#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.

Answer 1

The z-scale of your graph is much smaller than the x- and y- scales. Matplotlib requires that you set them manually if you want them to be equal, Your graph is probably correct.

That said, you should probably sample Venus' position more often than once a month if you don't want your graph to do that spirograph thing.

Answer 2

In def Keplar2CartesianSecond(...):, these equations are questionable.

W = LongofAscendingNode

p = LongofAscendingNode + ArgPeriapsis

w = p - W

The reason is that when we substitute 'p' into next calculation, we get:

w = LongofAscendingNode + ArgPeriapsis - LongofAscendingNode

which reduces down so simply

w = ArgPeriapsis

Note, I have not gone back to check accuracy of your translation of the equations versus their source.

---

6/22/24

I have gone back to several references (listed below). The cartesian calculation is correctly done. Just the translation from whatever the source(s) were was misleading.

Argument of Periapsis is [(lower case) omega] in scientific equations. Here, w is used because it looks almost like [l/c omega].

Since ArgPeriapsis is a formal function parameter, there is no need to do the calculation of p (aka LongofPeriapsis).

REF:

space.stackexchange.com/questions/8911

ssd.jpl.nasa.gov/planets/approx_pos.html