2

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

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

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

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

这似乎工作正常,只要我将绘图保持在轨道平面上,它本身就可以正确绘图。但我想将绘图旋转到参考框架(父体),以便我可以看到轨道外观的酷 3D 视图(obvs)。

现在,我怀疑这个错误在于我如何从轨道平面中的两个状态向量转换为将它们旋转到参考系。我正在使用本文档第 6 步中的方程式创建以下代码(但应用单独的旋转矩阵 [从此处复制]):

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))

(对于上下文,这是我用于轨道模型的完整类。)

当我从结果中绘制 X 和 Y 坐标时,我得到了这个:

在此处输入图像描述

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

在此处输入图像描述

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

在此处输入图像描述

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


更新:基于随机的回答,我转置了矩阵(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 系统的数据,我使用来自国际空间站 (2457380.183935185 = AD 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比我预期的要大。(我目前正在float12864 位系统上使用数字运行我的所有 numpy 计算)。

此外,生成的轨道仍然看起来像上面的(非常)偏心的第一个图(尽管我知道这个 LEO ISS 轨道是非常圆形的)。因此,对于问题的根源可能是什么,我感到有些困惑。

4

2 回答 2

2

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

在更仔细地查看了您正在执行的轨道模拟之后(请参阅评论中的此附加文档),我认为主要问题是最初非常合理但不真实的假设,即最终图应该看起来像一个椭圆。一般来说,它不会,因为轨道物体不一定会停留在一个平面上。

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

关于转置旋转矩阵

您引用的文档没有直接指定 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。
  • 看看乘以 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定义中的符号移到第二行正弦项,如上图

于 2015-12-23T12:14:34.247 回答
0

我会将随机指标的答案标记为正确,因为 a) 他的帮助值得加分,b) 他的建议基本上是正确的。

然而,奇怪情节的来源实际上最终是链接Orbit类中的这些行:

    self.v_position = self.rotate(v_position, self.long_of_asc_node, self.inclination, self.argument_of_periapsis)
    self.v_velocity = self.rotate(v_velocity, self.long_of_asc_node, self.inclination, self.argument_of_periapsis)

请注意,self.v_position在调用旋转速度矢量之前,该属性已更新;人们可能还会注意到,在阅读代码时,我聪明地决定将所有轨道元素值方法都包装在@property装饰器中,以使计算更加清晰。

但当然,这也意味着每次访问属性时都会调用方法——并重新计算值。因此,第二次调用self.rotate()的轨道元素值与第一次调用略有不同,更重要的是,这些值与“当前”位置和速度状态向量未 100% 正确匹配!

因此,在对这个 bug 进行了几天的努力后,我从我正在以重构的形式进行的一些牦牛剃须中弄清楚了这一点,现在一切正常。

于 2015-12-24T23:08:16.423 回答