0

我正在尝试对一个被月球环绕的行星进行简单的模拟。到目前为止，我有一个解决行星和月球轨道的 2 体问题。现在我想给行星添加一个固定的旋转轴，看看它是如何受到月球的影响的。知道如何通过使用 python 来完成吗？

可以使用以下代码运行两个主体问题：

import pylab
import math
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# Set Constants
G = 6.67e-11
AU = 1.5e11
daysec = 24.0*60*60

Ma =5.972e24   # Planet mass in Kg
Mb = 7.348e22   # Moon mass in Kg


gravconst = G*Ma*Mb

# Set up starting conditions

# Planet
xa = 0.0
ya = 0.0
za = 0.0

xva = 0.0
yva = 0.0
zva = 0.0

# Moon
xb = 384400000
yb = 0.0
zb = 0.0

xvb = 0.0
yvb = 1000.0
zvb = 0.0

# Time steps
t = 0.0
dt = 0.01*daysec

# Coordinate lists
xalist = []
yalist = []

xblist = []
yblist = []

zalist = []
zblist = []

# Loop
while t < 100.0*daysec:
    # Compute Force
    rx = xb-xa
    ry = yb-ya
    rz = zb-za

    modr3 = (rx**2+ry**2+rz**2)**1.5

    fx = -gravconst*rx/modr3
    fy = -gravconst*ry/modr3
    fz = -gravconst*rz/modr3

    # Update quantities
    # Moon
    xvb += fx*dt/Mb
    yvb += fy*dt/Mb
    zvb += fz*dt/Mb

    xb += xvb*dt
    yb += yvb*dt
    zb += zvb*dt

    # Planet
    xva += -fx*dt/Ma
    yva += -fy*dt/Ma
    zva += -fz*dt/Ma

    xa += xva*dt
    ya += yva*dt
    za += zva*dt

    t += dt

    # Saving them in lists
    xalist.append(xa)
    yalist.append(ya)
    zalist.append(za)

    xblist.append(xb)
    yblist.append(yb)
    zblist.append(zb)


xalist[:] = [x / 1e6 for x in xalist]
yalist[:] = [x / 1e6 for x in yalist]
zalist[:] = [x / 1e6 for x in zalist]

xblist[:] = [x / 1e6 for x in xblist]
yblist[:] = [x / 1e6 for x in yblist]
zblist[:] = [x / 1e6 for x in zblist]

#Creating the point to represent the planet at the origin (not to       scale),
plt.scatter(0,0,s=200,color='blue')
plt.annotate('Planet', xy=(-45,-50))
plt.scatter(xblist[0],0,s=100,color='grey')
plt.annotate('Mond', xy=(xblist[0]-45,-50))

# Plotting
pylab.plot(xalist, yalist, "-g")
pylab.plot(xblist, yblist, "-r")
plt.axhline(0, color='black')
plt.axvline(0, color='black')
pylab.axis("equal")
pylab.xlabel("X (Mio. Meter)")
pylab.ylabel("Y (Mio. Meter)")
pylab.show()
4

1 回答 1

0

不是答案，因为我不是这方面的专家，而只是一些提示（以评论的形式无法阅读）

您要添加的内容非常复杂，因为您需要考虑：

  1. 两个物体的移动质量

    所以你需要具有任何移动质量（如海洋、岩浆、旋转核心等）的物体的“接触面高度”，这样你就可以计算每次的真实重心。此外，您还需要对移动的质量本身施加力（这不仅由重力和旋转驱动，而且还由共振和主要惯性驱动）不要忘记地球周围也有核心和岩浆，而不仅仅是海洋，所以你需要考虑至少考虑 3 个表面...

  2. 两个物体的质量分布不均匀

    因此您可以计算重心，以及相对于实际旋转轴的旋转的二次质量惯性

  3. 行星/月亮通常至少有 3 个身体问题，而不仅仅是 2 个

    因为当地的恒星也对月球轨道影响很大......

根据数量的不同，一些效果会非常小以至于可以丢弃，但有些则不能（尤其是在惯性+共振到位的情况下）。

旋转方程类似于您已经得到的位置/速度/加速度。它称为Newton D'Alembert 积分/物理，但您需要实现变换矩阵才能做到这一点。

查看一些相关的质量保证：

请注意最后一个链接以提高您的积分的准确性，因为现在它非常糟糕，并且dt由于重力矢量随时间变化但现在您以错误的方向影响它（即仅在 dt 间隔开始时是正确的）对于每个积分迭代。

正如你所看到的，有很多东西需要处理，而我看到的大多数模拟程序都没有这样做（包括我的）......他们通过章动和进动常数来伪造它。

于 2019-04-05T08:31:13.047 回答