python - 二体轨道建模问题

Question

如果您不想阅读太多背景信息，请跳至下面的更新 2。

我正在尝试实现一个用于简单轨道模拟(两个物体)的模型。

但是，当我尝试使用我编写的代码时，从结果生成的图看起来很奇怪。

程序使用初始状态向量(位置和速度)来计算开普勒轨道元素，然后使用这些元素计算下一个位置，并作为接下来的两个状态向量返回。

这似乎工作正常，并且只要我将绘图保持在轨道平面上，它本身就可以正确绘制。但我想将绘图旋转到引用系(父体)，以便我可以看到轨道外观的酷炫 3D View (obvs)。

现在，我怀疑问题在于如何将轨道平面中的两个状态向量转换为将它们旋转到引用系。我使用 this document 第 6 步中的方程创建以下代码(但应用单独的旋转矩阵 [ copied from here ]):

from numpy import sin, cos, matrix, newaxis, asarray, squeeze, dot

def Rx(theta):
    """
    Return a rotation matrix for the X axis and angle *theta*
    """
    return matrix([
        [1, 0,           0           ],
        [0, cos(theta), -sin(theta)  ],
        [0, sin(theta), cos(theta)   ],
    ], dtype="float64")

def Rz(theta):
    """
    Return a rotation matrix for the Z axis and angle *theta*
    """
    return matrix([
        [cos(theta), -sin(theta),   0],
        [sin(theta), cos(theta),    0],
        [0,          0,             1],
    ], dtype="float64")

def rotate1(vector, O, i, w):
    # The starting value of *vector* is just a 1-dimensional numpy
    # array.
    # Transform into a column vector.
    vector = vector[:, newaxis]
    # Perform the rotation
    R = Rz(-O) * Rx(-i) * Rz(-w)
    res2 = dot(R, vector)
    # Transform back into a row vector (because that's what
    # the rest of the program uses)
    return squeeze(asarray(res2))

(作为上下文，我使用 this is the full class 作为轨道模型。)

当我根据结果绘制 X 和 Y 坐标时，我得到:

但是当我将旋转矩阵更改为R = Rz(-O) * Rx(-i)时，我得到了这个更合理的图(尽管明显缺少一次旋转，并且稍微偏离中心) ):

当我进一步将其减少到 R = Rx(-i) 时，正如人们所期望的那样，我得到了这个:

正如我所说，我相当确定不是轨道计算代码表现得很奇怪，而是旋转代码中出现了一些错误。但我不确定在哪里缩小范围，因为我对 numpy 和矩阵数学总体来说都很陌生。

<小时/> 更新:基于 stochastic's answer，我转置了矩阵( R = Rz(-O).T * Rx(-i).T * Rz(-w).T) ，但后来得到了这个情节:

这让我想知道我对屏幕坐标的转换是否在某种程度上是错误的 - 但它对我来说看起来是正确的(并且与旋转较少的更正确的绘图相同的代码)即:

def recenter(v_position, viewport_width, viewport_height):
    x, y, z = v_position
    # the size of the viewport in meters
    bounds = 20000000
    # viewport_width is the screen pixels (800)
    scale = viewport_width/bounds
    # Perform the scaling operation
    x *= scale
    y *= scale
    # recenter to screen X and Y measured from the top-left corner
    # of the viewport
    x += viewport_width/2
    y = viewport_height/2 - y
    # Cast to int, because we don't care about pixel fractions
    return int(x), int(y)

<小时/> 更新2

虽然我在随机的帮助下对方程的实现以及旋转进行了三次检查，但我仍然无法得到正确的轨道。它们看起来仍然与上面的图中基本相同。

使用来自 NASA Horizo​​n 系统的数据，我设置了一条具有来自 ISS 的特定状态向量的轨道 (2457380.183935185 = A.D. 2015-Dec-23 16:24:52.0000 (TDB))，并根据开普勒轨道元素对其进行了检查在同一时刻，产生以下结果:

inclination :
  0.900246137041
  0.900246137041
true_anomaly :
  0.11497063007
  0.0982485984565
long_of_asc_node :
  3.80727461492
  3.80727461492
eccentricity :
  0.000429082122137
  0.000501850615905
semi_major_axis :
  6778560.7037
  6779057.01374
mean_anomaly :
  0.114872215066
  0.0981501816537
argument_of_periapsis :
  0.843226618347
  0.85994864996

顶部的值是我的(计算的)值，底部的值是 NASA 的值。显然，一些浮点精度误差是可以预料的，但 mean_anomaly 和 true_anomaly 的变化确实让我觉得比我预期的要大。 (我目前正在 64 位系统上使用 float128 数字运行所有 numpy 计算)。

此外，生成的轨道仍然看起来像上面的(相当)偏心的第一个图(尽管我知道这个 LEO ISS 轨道是相当圆形的)。所以我有点困惑问题的根源可能是什么。

Answer 1

我相信您至少有两个问题。

在更仔细地研究了您正在做的轨道模拟之后(请参阅评论中的 this additional document )，我认为主要问题是最初非常合理但不真实的假设，即最终情节应该看起来像一个椭圆。一般来说不会，因为绕轨道运行的物体不一定会停留在一个平面上。

我认为另一个问题是，根据您描述的文档(见下文)，您的旋转矩阵是它们应有的转置。

转置旋转矩阵

您引用的文档没有直接指定 R_x 和 R_z 是否应该是轴的右旋或它们将相乘的矢量的右旋，尽管您可以从方程 9(或 10)中计算出来。事实证明，它们应该是轴的右旋，而不是矢量。这意味着它们应该这样定义:

    return matrix([
        [1, 0,           0           ],
        [0, cos(theta), sin(theta)  ],
        [0,-sin(theta), cos(theta)   ],
    ], dtype="float64")

而不是像这样:

    return matrix([
        [1, 0,           0           ],
        [0, cos(theta),-sin(theta)  ],
        [0, sin(theta), cos(theta)   ],
    ], dtype="float64")

我通过在纸上手工重现方程 9 发现了这一点。

• 在该方程中，查看向量r(t) 的第一个分量。
• 有两个术语:一个包含 o_x，另一个包含 o_y。
• 看看那个倍增的东西。它是:-(sin(omega)*cos(Omega)+cos(omega)*cos(i)*sin(Omega))。
• 前面的减号是关键。它来自 Rz 矩阵第一行的减号。
• 由于等式 9 中的 Omega、i 和 omega 均被否定，这意味着减号需要位于 R_z 的第二行，这意味着 R_z 表示轴的右旋，而不是向量。
• 类似地，我们可以查看最后一项的 o_y 分量，发现减号需要位于 R_x 的第二行，这意味着(谢天谢地) R_z 和 R_x 都是右手旋转轴。

您的 Rx 和 Rz 函数当前正在定义向量的右手旋转，而不是轴。

您可以通过以下任一方式解决此问题(所有三个都是等效的):

• 删除欧拉角上的减号:Rz(O) * Rx(i) * Rz(w)

• 转置旋转矩阵:Rz(-O).T * Rx(-i).T * Rz(-w).T

• 将 Rx 和 Rz 定义中的 - 符号移至第二行正弦项，如上所示