行星倾角，用python追踪旋转轴

Question

我正试着做一个简单的模拟一个被月球环绕的行星。到目前为止，我有一个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()

Answer 1

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

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

1. 两个物体的运动质量

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

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

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

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

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

根据数字的不同，有些影响会很小，可以忽略，但有些影响不会（特别是在惯性+共振的情况下）。你知道吗

旋转方程类似于你已经得到的位置/速度/加速度。它被称为Newton D'Alembert integration/physics，但是您需要实现transform matrices来实现这一点。你知道吗

参见一些相关的QAs：

• Rigid Body Physics Rotations
• realistic n-body solar system

注意提高积分精度的最后一个环节，因为现在这个环节非常糟糕，无论dt有多小，轨道都会变形，因为重力矢量会随时间而变化，但现在每个积分的影响方向都是错误的（只有在dt间隔开始时才是正确的）迭代。你知道吗

正如你所看到的，这是很多东西来处理和大多数模拟程序，我看到不这样做（包括我的）。。。他们用章动常数和进动常数来模拟。你知道吗