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.