0
$\begingroup$

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: enter image description here 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.

Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ Please clearly state what your question is. Is it "please debug my code?" $\endgroup$
    – Organic Marble
    Aug 10, 2022 at 11:24
  • $\begingroup$ Also consider adding links to the pages you got the equations from. $\endgroup$
    – Greg Miller
    Aug 10, 2022 at 12:42
  • $\begingroup$ @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? $\endgroup$
    – finnnosam
    Aug 10, 2022 at 20:33
  • $\begingroup$ @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. $\endgroup$
    – finnnosam
    Aug 10, 2022 at 21:03
  • 1
    $\begingroup$ 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). $\endgroup$
    – Ryan C
    Aug 10, 2022 at 21:22

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ oh, you are absolutely correct, after doing what you said it looks fine. Sorry for bringing this problem here then. Thank you. $\endgroup$
    – finnnosam
    Aug 11, 2022 at 0:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.