【问题标题】:Elliptical orbit in vpythonvpython中的椭圆轨道
【发布时间】:2019-04-05 14:16:15
【问题描述】:

我有以下代码。此代码是模拟围绕其他物体的轨道物体,例如太阳系。当你运行它时,物体以圆形轨迹运行。

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

我想知道如何将轨道路径更改为椭圆? 我尝试了几种方法,但我无法找到任何解决方案。

谢谢。 谢谢

【问题讨论】:

  • 它们已经是椭圆的了,只是偏心率太小了,你把它看成一个圆。您需要在不更改任何其他参数的情况下(显着)修改动能。尝试将 vE 硬编码为 100km/s 这样荒谬的东西,你会看到它发生了变化
  • @vencaslac 顺便说一句,我改变了一切,但它仍然是循环的。 :(
  • 既然你在画圆圈,你就会得到圆圈。看看这个:tinyurl.com/vporbit,然后点击“View this program”查看VPython代码。
  • @user1114907 在那个页面我看到ReferenceError: glowscript_compile is not defined
  • 这很奇怪,因为您的 VPython 程序使用 VPython 7 运行。您的计算机环境是什么?你使用的是什么浏览器?你可以试试不同的电脑。或者,您可以通过从github.com/BruceSherwood/vpython-jupyter 下载这些演示程序的 VPython 7 版本,在文件夹 Demos (Jupyter) 或 Demos_no_notebook 中。

标签: python python-3.x vpython


【解决方案1】:

这似乎更像是一个物理问题,而不是编程问题。问题是您在计算速度和线性积分位置时假设每个轨道都是圆形的(例如 v * dt)。这不是你计算轨道物体轨迹的方式。

为了简单起见,我们假设所有质量都是点质量,因此没有任何奇怪的重力梯度或姿态动力学需要解释。

从那里,您可以参考这个 MIT 页面。 (http://web.mit.edu/12.004/TheLastHandout/PastHandouts/Chap03.Orbital.Dynamics.pdf) 关于轨道动力学。在第 7 页,有一个方程将您的中心体的径向位置与多个轨道参数的函数联系起来。除了轨道的偏心率之外,您似乎拥有所有参数。如果您有详细的临时数据或远点/近点信息,您可以在线查找或计算它。

从该等式中,您将在分母中看到一个 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

【讨论】:

    猜你喜欢
    • 2012-09-04
    • 1970-01-01
    • 1970-01-01
    • 2021-07-21
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多