2015-12-23 25 views
1

如果您不想讀取太多背景,請跳到下面的更新2。兩體軌道建模問題

我正在嘗試實現一個簡單的軌道模擬模型(兩個主體)。

但是,當我嘗試使用我寫的代碼時,從結果生成的圖很奇怪。

該程序使用初始狀態向量(位置和速度)來計算開普勒軌道元素,然後計算下一個位置,並返回爲下兩個狀態向量。

這似乎工作正常,並且本身,只要我保持在軌道平面上的陰謀,繪圖正確。但是我想旋轉繪圖到參考框架(父體),以便我可以看到軌道看起來像什麼樣的酷3D視圖(obvs)。

現在,我懷疑這個錯誤在於我如何從軌道平面中的兩個狀態向量轉換爲旋轉到參考系。我使用的方程從this document步驟6來創建從下面的代碼(但施加個別輪替] matricies [copied from here]):(對於上下文,this is the full class我使用了軌道模型)

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從結果座標,我得到這個:

enter image description here

但是,當我旋轉矩陣更改爲R = Rz(-O) * Rx(-i) ,我得到這個更合理的情節(雖然明顯缺失的一轉,和稍微偏離中心):

enter image description here

當我進一步降低它R = Rx(-i),正如人們所期望的,我得到這樣的:

enter image description here

所以,正如我所說,我相當確信,這不是軌道計算代碼行爲怪異,而是旋轉代碼中的一些錯誤。但我不確定在哪裏縮小這個範圍,因爲我對一般的數學和矩陣數學都很陌生。


更新:基於 stochastic's answer我調換了matricies( R = Rz(-O).T * Rx(-i).T * Rz(-w).T),但後來得到這個情節:

enter image description here

這讓我想知道如果我轉換到屏幕座標,在某種程度上錯了 - - 但它對我來說看起來是正確的(並且與具有較少旋轉的更正確的圖相同的代碼),即:

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

雖然我有三重檢查我的執行方程,以及與隨機的幫助下旋轉,我仍然不能得到軌道出來的權利。他們仍然看起來基本上與上面的情節相同。我使用來自美國國家航空航天局地平線系統的數據,從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_anomalytrue_anomaly的變化的確讓我覺得比我預期的要大。 (我目前正在使用64位系統上的float128號碼運行我所有的numpy計算)。

此外,所產生的軌道仍然看起來像(相當)偏心的第1圖,上面(儘管我知道這LEO軌道ISS是相當圓形)。所以我有點難以確定問題的根源。

回答

1

我相信你至少有兩個問題。

在仔細觀察你正在進行的軌道模擬之後(見註釋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。
  • 看看multliplying o_y的東西。它是:-(sin(omega)*cos(Omega)+cos(omega)*cos(i)*sin(Omega))
  • 領先的負號是關鍵。它來自Rz矩陣第一行的減號。
  • 由於歐米茄,i和在式(9)的ω都否定的,這意味着減號需要是對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定義爲第二行正弦項,如上所示

+0

謝謝,這確實有點澄清 - 但不幸的是我仍然得到意想不到的結果。我在上面添加了更多可能會澄清事情的信息。 –

+0

另外,抱歉厚厚,但你說:「事實證明,他們是軸的旋轉,而不是矢量」,我不明確的含義。這是否意味着有一箇中間步驟,而不是直接將旋轉矩陣應用於矢量? –

+1

沒問題。每個旋轉矩陣可以被解釋爲兩種方式:它可以被認爲是在一個方向上旋轉矢量或者在相反方向上旋轉軸。這裏唯一的問題是你的旋轉矩陣定義中的正弦項在它前面應該有一個負號。我會盡力澄清一些更多的答案 – stochastic

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屬性調用之前更新的旋轉速度矢量發生;人們在閱讀代碼時可能會注意到,我巧妙地決定使裝飾器中包含的所有軌道元素值方法更清晰。

當然不過,這也意味着該方法被稱爲 - 和值重新計算 - 每一個屬性被訪問的時間。因此,第二次調用self.rotate()時,第一次調用的軌道元素的值略有不同,更重要的是,值與「當前」位置和速度狀態向量無法正確匹配100%!

所以經過幾天我的頭撞到這個bug後,我從一些犛牛剃鬚中找出了我正在做的一個重構的形式,現在它完美地工作。