Two-body 轨道建模问题
Two-body orbit modelling problems
如果您不想阅读太多背景知识，请跳至下面的更新 2。
我正在尝试实现一个简单轨道模拟模型（双体）。
但是，当我尝试使用我编写的代码时，结果生成的图表看起来很奇怪。
程序使用初始状态向量（位置和速度）来计算开普勒轨道元素，这些元素用于计算下一个位置，并作为下两个状态向量返回。
这似乎工作正常，并且只要我将绘图保持在轨道平面上，它本身就能正确绘图。但我想将绘图旋转到参考系（母体），这样我就可以看到一个很酷的轨道外观的 3D 视图（obvs）。
现在，我怀疑错误在于我如何从轨道平面中的两个状态向量转换为将它们旋转到参考系。我正在使用 this document to create the following code from (but applying individual roation matricies [copied from here]):
步骤 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))
（对于上下文，this is the full class 我正在使用轨道模型。）
当我根据结果绘制 X 和 Y 坐标时，我得到：
但是当我将旋转矩阵更改为R = Rz(-O) * Rx(-i)时，我得到了这个更合理的图（虽然显然缺少一个旋转，并且稍微off-center）：
当我将它进一步减少到 R = Rx(-i) 时，正如人们所期望的那样，我得到了这个：
所以正如我所说，我相当确定不是轨道计算代码表现异常，而是旋转代码中的一些错误。但我不确定在哪里缩小范围，因为我对一般的 numpy 和矩阵数学都很陌生。
更新：基于我转置了矩阵（R = Rz(-O).T * Rx(-i).T * Rz(-w).T），但后来得到了这个图：
这让我想知道我转换到屏幕坐标是否有问题——但它对我来说看起来是正确的（并且与 more-correct 图的代码相同，但旋转较少）即：
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
尽管我已经 triple-checked 实现了方程式，并且在随机变量的帮助下进行了旋转，但我仍然无法得到正确的轨道。它们仍然与上图中的基本相同。
使用来自 NASA Horizon 系统的数据，我建立了一个轨道，其中包含来自国际空间站的特定状态向量 (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 的变化确实让我感到震惊，因为它比我预期的要大。 （我目前 运行 我所有的 numpy 计算都是在 64 位系统上使用 float128 数字）。
此外，由此产生的轨道看起来仍然像上面的（相当）偏心的第一个图（即使我知道这个 LEO ISS 轨道是很圆的）。所以我对问题的根源可能是什么感到有点困惑。
我相信你至少有两个问题。
在更仔细地研究了您正在进行的轨道模拟之后（参见评论中的 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。
- 看看乘积的东西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定义中的-符号移至第二行正弦项，如上图
我要把 stochastic 的回答标记为正确，因为 a) 他的帮助值得加分，b) 他的建议基本上是正确的。
然而，奇怪情节的来源实际上最终是链接中的这些行 Orbit class:
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() 时轨道元素的值与第一次调用略有不同，更重要的是，这些值与 "current" 位置和速度状态不 100% 正确匹配矢量！
因此，经过几天的努力解决这个错误后，我以重构的形式进行了一些 yak 剃须，现在一切正常。
如果您不想阅读太多背景知识，请跳至下面的更新 2。
我正在尝试实现一个简单轨道模拟模型（双体）。
但是，当我尝试使用我编写的代码时，结果生成的图表看起来很奇怪。
程序使用初始状态向量（位置和速度）来计算开普勒轨道元素，这些元素用于计算下一个位置，并作为下两个状态向量返回。
这似乎工作正常，并且只要我将绘图保持在轨道平面上，它本身就能正确绘图。但我想将绘图旋转到参考系（母体），这样我就可以看到一个很酷的轨道外观的 3D 视图（obvs）。
现在，我怀疑错误在于我如何从轨道平面中的两个状态向量转换为将它们旋转到参考系。我正在使用 this document to create the following code from (but applying individual roation matricies [copied from here]):
步骤 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))
（对于上下文，this is the full class 我正在使用轨道模型。）
当我根据结果绘制 X 和 Y 坐标时，我得到：
但是当我将旋转矩阵更改为R = Rz(-O) * Rx(-i)时，我得到了这个更合理的图（虽然显然缺少一个旋转，并且稍微off-center）：
当我将它进一步减少到 R = Rx(-i) 时，正如人们所期望的那样，我得到了这个：
所以正如我所说，我相当确定不是轨道计算代码表现异常，而是旋转代码中的一些错误。但我不确定在哪里缩小范围，因为我对一般的 numpy 和矩阵数学都很陌生。
更新：基于
R = Rz(-O).T * Rx(-i).T * Rz(-w).T），但后来得到了这个图：
这让我想知道我转换到屏幕坐标是否有问题——但它对我来说看起来是正确的（并且与 more-correct 图的代码相同，但旋转较少）即：
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
尽管我已经 triple-checked 实现了方程式，并且在随机变量的帮助下进行了旋转，但我仍然无法得到正确的轨道。它们仍然与上图中的基本相同。
使用来自 NASA Horizon 系统的数据，我建立了一个轨道，其中包含来自国际空间站的特定状态向量 (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 的变化确实让我感到震惊，因为它比我预期的要大。 （我目前 运行 我所有的 numpy 计算都是在 64 位系统上使用 float128 数字）。
此外，由此产生的轨道看起来仍然像上面的（相当）偏心的第一个图（即使我知道这个 LEO ISS 轨道是很圆的）。所以我对问题的根源可能是什么感到有点困惑。
我相信你至少有两个问题。
在更仔细地研究了您正在进行的轨道模拟之后（参见评论中的 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。
- 看看乘积的东西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定义中的
-符号移至第二行正弦项，如上图
我要把 stochastic 的回答标记为正确，因为 a) 他的帮助值得加分，b) 他的建议基本上是正确的。
然而，奇怪情节的来源实际上最终是链接中的这些行 Orbit class:
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() 时轨道元素的值与第一次调用略有不同，更重要的是，这些值与 "current" 位置和速度状态不 100% 正确匹配矢量！
因此，经过几天的努力解决这个错误后，我以重构的形式进行了一些 yak 剃须，现在一切正常。