如果您不想阅读太多背景信息,请跳至下面的更新 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 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
比我预期的要大。(我目前正在float128
64 位系统上使用数字运行我的所有 numpy 计算)。
此外,生成的轨道仍然看起来像上面的(非常)偏心的第一个图(尽管我知道这个 LEO ISS 轨道是非常圆形的)。因此,对于问题的根源可能是什么,我感到有些困惑。