背景をあまり読みたくない場合は、以下の更新 2 にスキップしてください。
単純な軌道シミュレーション (2 体) のモデルを実装しようとしています。
しかし、私が書いたコードを使用しようとすると、結果から生成されたプロットはかなり奇妙に見えます。
プログラムは、初期状態ベクトル (位置と速度) を使用してケプラー軌道要素を計算します。これは次の位置を計算するために使用され、次の 2 つの状態ベクトルとして返されます。
これは正常に機能しているようで、プロットを軌道面に保持している限り、それ自体で正しくプロットされます。しかし、軌道がどのように見えるか（obvs）のクールな3Dビューを見ることができるように、プロットを参照フレーム（親体）に回転させたいと思います。
現時点では、バグは軌道面の 2 つの状態ベクトルから参照フレームへの回転に変換する方法にあると思われます。このドキュメントのステップ 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)、より妥当なプロットが得られます (ただし、明らかに 1 つの回転が失われ、中心が少しずれています)。
そして、それをさらに に減らすと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 のシステムからのデータを使用して、ISS からの特定の状態ベクトル (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 軌道がかなり円形であることはわかっていますが)。そのため、問題の原因が何であるかについて少し困惑しています。