【发布时间】：2016-03-29 20:02:25
【问题描述】：
如果您不想阅读太多背景信息，请跳至下面的更新 2。
我正在尝试为简单的轨道模拟（两个物体）实现一个模型。
但是，当我尝试使用我编写的代码时，结果生成的图看起来很奇怪。
程序使用初始状态向量（位置和速度）计算开普勒轨道元素，然后使用这些元素计算下一个位置，并作为接下来的两个状态向量返回。
这似乎工作正常，只要我将绘图保持在轨道平面上，它本身就可以正确绘图。但我想将绘图旋转到参考框架（父体），以便我可以看到轨道外观的酷 3D 视图（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 Horizon 系统的数据，我使用来自国际空间站的特定状态向量（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 的变化确实比我预期的要大。 （我目前在 64 位系统上使用 float128 数字运行我所有的 numpy 计算）。
此外，生成的轨道仍然看起来像上面的（非常）偏心的第一个图（尽管我知道这个 LEO ISS 轨道是非常圆形的）。所以我有点不知道问题的根源是什么。
【问题讨论】：
标签： python numpy matrix vector orbital-mechanics