我有以下代码。此代码是围绕其他物体（例如，太阳系。在运行时，对象沿圆形轨迹运动。
import math
from vpython import *
lamp = local_light(pos=vector(0,0,0), color=color.yellow)
# Data in units according to the International System of Units
G = 6.67 * math.pow(10,-11)
# Mass of the Earth
ME = 5.973 * math.pow(10,24)
# Mass of the Moon
MM = 7.347 * math.pow(10,22)
# Mass of the Mars
MMa = 6.39 * math.pow(10,23)
# Mass of the Sun
MS = 1.989 * math.pow(10,30)
# Radius Earth-Moon
REM = 384400000
# Radius Sun-Earth
RSE = 149600000000
RMS = 227900000000
# Force Earth-Moon
FEM = G*(ME*MM)/math.pow(REM,2)
# Force Earth-Sun
FES = G*(MS*ME)/math.pow(RSE,2)
# Force Mars-Sun
FEMa = G*(MMa*MS)/math.pow(RMS,2)
# Angular velocity of the Moon with respect to the Earth (rad/s)
wM = math.sqrt(FEM/(MM * REM))
# Velocity v of the Moon (m/s)
vM = wM * REM
print("Angular velocity of the Moon with respect to the Earth: ",wM," rad/s")
print("Velocity v of the Moon: ",vM/1000," km/s")
# Angular velocity of the Earth with respect to the Sun(rad/s)
wE = math.sqrt(FES/(ME * RSE))
# Angular velocity of the Mars with respect to the Sun(rad/s)
wMa = math.sqrt(FEMa/(MMa * RMS))
# Velocity v of the Earth (m/s)
vE = wE * RSE
# Velocity v of the Earth (m/s)
vMa = wMa * RMS
print("Angular velocity of the Earth with respect to the Sun: ",wE," rad/s")
print("Velocity v of the Earth: ",vE/1000," km/s")
# Initial angular position
theta0 = 0
# Position at each time
def positionMoon(t):
theta = theta0 + wM * t
return theta
def positionMars(t):
theta = theta0 + wMa * t
return theta
def positionEarth(t):
theta = theta0 + wE * t
return theta
def fromDaysToS(d):
s = d*24*60*60
return s
def fromStoDays(s):
d = s/60/60/24
return d
def fromDaysToh(d):
h = d * 24
return h
# Graphical parameters
print("\nSimulation Earth-Moon-Sun motion\n")
days = 365
seconds = fromDaysToS(days)
print("Days: ",days)
print("Seconds: ",seconds)
v = vector(384,0,0)
E = sphere(pos = vector(1500,0,0), color = color.blue, radius = 60, make_trail=True)
Ma = sphere(pos = vector(2300,0,0), color = color.orange, radius = 30, make_trail=True)
M = sphere(pos = E.pos + v, color = color.white,radius = 10, make_trail=True)
S = sphere(pos = vector(0,0,0), color = color.yellow, radius=700)
t = 0
thetaTerra1 = 0
dt = 5000
dthetaE = positionEarth(t+dt)- positionEarth(t)
dthetaM = positionMoon(t+dt) - positionMoon(t)
dthetaMa = positionMars(t+dt) - positionMars(t)
print("delta t:",dt,"seconds. Days:",fromStoDays(dt),"hours:",fromDaysToh(fromStoDays(dt)),sep=" ")
print("Variation angular position of the Earth:",dthetaE,"rad/s that's to say",degrees(dthetaE),"degrees",sep=" ")
print("Variation angular position of the Moon:",dthetaM,"rad/s that's to say",degrees(dthetaM),"degrees",sep=" ")
while t < seconds:
rate(500)
thetaEarth = positionEarth(t+dt)- positionEarth(t)
thetaMoon = positionMoon(t+dt) - positionMoon(t)
thetaMars = positionMars(t+dt) - positionMars(t)
# Rotation only around z axis (0,0,1)
E.pos = rotate(E.pos,angle=thetaEarth,axis=vector(0,1,0))
Ma.pos = rotate(Ma.pos,angle=thetaMars,axis=vector(0,1,0))
v = rotate(v,angle=thetaMoon,axis=vector(0,1,0))
M.pos = E.pos + v
t += dt
我想知道如何将轨道转换为椭圆形吗？
我尝试了几种方法，但是找不到任何解决方案。
谢谢。
谢谢
最佳答案
与编程问题相比，这似乎更像是物理问题。问题是，您在计算速度和线性积分位置（例如v * dt）时假设每个轨道都是圆形的。这不是计算轨道物体的轨迹的方式。
为了简单起见，我们将假定所有质量都是点质量，因此没有任何奇怪的重力梯度或姿态动力学可以解释。
从那里，您可以参考此MIT页面。 （http://web.mit.edu/12.004/TheLastHandout/PastHandouts/Chap03.Orbital.Dynamics.pdf）在轨道动力学上。在第7页上，存在一个方程，该方程将您的中心体的径向位置与多个轨道参数的函数相关联。似乎您拥有除轨道偏心率之外的所有参数。如果您有详细的短暂数据或Apapsis / periapsis信息，则可以在线查找或进行计算。
从该方程式中，您将看到分母中的phi-phi_0项。通俗地讲，这就是卫星的真正异常。您可以从0到360迭代这个真实的异常参数来找到您的径向距离，而无需花费时间，而从真实的异常，倾斜度，到上升节点的直角以及橄榄石的参数，您可以在以下位置找到3D笛卡尔坐标一个特定的真实异常。
从真正的异常情况来看是不那么琐碎的。您需要先找到偏心异常，然后再找到每个偏心异常步骤的均值异常。现在，您的平均数异常随时间而变。您可以在用v * dt计算位置的“节点”之间线性插值。您可以使用vis-viva方程计算速度，而dt将是计算出的时间步长之间的差。
在每个时间步上，您都可以在python程序中更新卫星的位置，它会正确绘制您的轨迹。
有关真实异常的更多信息，维基百科对此有一个很好的描述：https://en.wikipedia.org/wiki/True_anomaly
有关轨道元素（将其从径向位置转换为笛卡尔坐标所必需的）的更多信息：https://en.wikipedia.org/wiki/Orbital_elements
关于python - VPython中的椭圆轨道，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/55537402/