使用Vpython模拟日地月轨道答案

【问题标题】：Using Vpython to simulate sun-earth-moon orbit使用Vpython模拟日地月轨道
【发布时间】：2015-03-12 10:41:00
【问题描述】：

我正在尝试使用牛顿第二定律来模拟日-地-月系统。我可以让月球和地球绕太阳平行旋转。但是月球不绕地球旋转。

我知道月亮和太阳之间应该有引力，月球和地球之间也应该有引力。所以我把加速度加在一起。但在 3D 可视化中，似乎没有地球对月球的影响。

谁能检查我在运动部分的代码：

def movePlanets(self):
    rate(200)

    G = (6.673e-11) 
    dt =  12*3600 # half day

    for p in self.planets:

        p.moveTo((p.getXPos() + dt * p.getXVel()),
                 (p.getYPos() + dt * p.getYVel()),
                 (p.getZPos() + dt * p.getZVel()))


        rx = self.thesun.getXPos() - p.getXPos()

        ry = self.thesun.getYPos() - p.getYPos()
        rz = self.thesun.getZPos() - p.getZPos()

        r = math.sqrt(rx**2 + ry**2 + rz**2)

        accx = G * self.thesun.getMass()*rx/r**3
        accy = G * self.thesun.getMass()*ry/r**3
        accz = G * self.thesun.getMass()*rz/r**3


        for pTwo in self.planets:
            if(pTwo != p):
                rx = pTwo.getXPos() - p.getXPos()
                ry = pTwo.getYPos() - p.getYPos()
                rz = pTwo.getZPos() - p.getZPos()
                r = math.sqrt(rx**2 + ry**2 + rz**2)

                accx = accx + (G * pTwo.getMass()*rx/r**3)
                accy = accy + (G * pTwo.getMass()*ry/r**3)
                accz = accz + (G * pTwo.getMass()*rz/r**3)

        p.setXVel(p.getXVel() + dt * accx)
        p.setYVel(p.getYVel() + dt * accy)
        p.setZVel(p.getZVel() + dt * accz)

非常感谢！

【问题讨论】：

标签： simulation vpython

【解决方案1】：

我已经完成了这个模拟，我认为你没有使用“足够”的物理，这使得代码更容易理解。老实说，我无法理解您的代码。这是方法（正如我的物理学教授教给我的）：

创建 3 个球体对象，分别为地球、太阳和月亮
每个对象都有以下属性：质量、位置、速度和合力

在主循环中：

while True:
    for each "planet":
        calculate net force (use the gravity formula)
        update momentum: P = P + net_force*dt
        update velocity: v = P/mass
        update position: pos = pos + v*dt

注意：以下变量（属性）应为矢量类型，以便使用矢量运算符（加法、减法和标量乘法）进行计算：位置、速度和 net_force

【讨论】：

感谢您的回答！我发现了我的代码的错误。这是因为我在月球上使用的比例因子。

【解决方案2】：

我认为可以轻松添加第三个“星”：

from visual import *

G = 6.7e-11

giant = sphere(pos=(-1e11,0,0), radius=2e10, color=color.red,
               make_trail=True, interval=10)
giant.mass = 2e30
giant.p = vector(0, 0, -1e4) * giant.mass

dwarf = sphere(pos=(1.5e11,0,0), radius=1e10, color=color.yellow,
               make_trail=True, interval=10)
dwarf.mass = 1e30
dwarf.p = -giant.p

dt = 1e5

while True:
  rate(200)

  dist = dwarf.pos - giant.pos
  force = G * giant.mass * dwarf.mass * dist / mag(dist)**3
  giant.p = giant.p + force*dt
  dwarf.p = dwarf.p - force*dt

  for star in [giant, dwarf]:
    star.pos = star.pos + star.p/star.mass * dt

希望这篇文章对你有帮助！

【讨论】：