■人工衛星の理論軌道の計算式について

 先の人工衛星の打上げシミュレーションプログラム "iEarth Satellite"においては、オイラー法およびルンゲ・クッタ法による数値計算により衛星の軌道を計算しています。
 ここでは水平方向の発射に限定して、人工衛星の軌道の理論式を示すとともに"iEarth Satellite"の改良版を紹介します。改良版では数値計算による軌跡に加えて、理論軌道も表示されます。

Satellite2.gif
左図で、打上げ角α=0度の水平方向発射の場合を考えます。

ただし、途中までの計算式は任意の発射方向に適用できます。
人工衛星の軌道の理論式の導出
・人工衛星の軌道に関する運動方程式(2階の微分方程式):   x"= -GMx/r3   y"= -GMy/r3   ここで、    G : 万有引力定数    M : 地球の質量    r : 地球の中心と人工衛星の距離、r2 = x2 + y2  を初期条件[時刻t=0 で x=r0, y=0, 速度v=v0]で解きます。 ・直交座標系で表現されたこれらの式を極座標系(r,θ)に変換します。両者の関係   x = rcosθ   y = rsinθ  より、   x' = r'cosθ - rsinθθ' x" = r"cosθ - 2r'sinθθ' -rcosθ(θ')2 - rsinθθ"   y' = r'sinθ + rcosθθ' y" = r"sinθ + 2r'cosθθ' -rsinθ(θ')2 + rcosθθ" ・加速度の動径(r)方向成分:Ar   Ar = x"cosθ + y"sinθ     = r" - r(θ')2  加速度の方位角(θ)方向成分:Aθ Aθ = -x"sinθ + y"cosθ     = 2r'θ' + rθ" = (r2θ')'/r ・力の動径(r)方向成分:Fr = -GMm/r2  力の方位角(θ)方向成分:Fθ = 0 ・従って、極座標系での運動方程式は   r" - r(θ')2 = -GM/r2 ・・・(式1) (r2θ')'/r = 0  後者の式より   r2θ' = c1 (c1: 積分定数)   mr2θ' = mc1  左辺は角運動量(mrv)であり、これが一定であることから t=0 での初期値 r = r0, v = v0 を代入して  (これ以降は打上げ角度α=0度の水平方向発射の場合に限定)   mr2θ' = mr0v0   r2θ' = r0v0 ・・・(式2) ・(式2)のθ'を(式1)に代入すると、動径r方向の2階の微分方程式が得られます。   r" -(r0v0)2/r3 + GM/r2 = 0 ・変数の置換 u = 1/r を行ない、整理すると(途中の過程は省略)   d2u/dθ2 + u = 1/L ここで、L = (r0v0)2/(GM) ・この解は   u = 1/r = Acos(θ+β) + 1/L (A, β: 積分定数) θが0のときrが最小(u が最大)となるように角度を設定すればβ=0となるから   u = 1/r = Acos(θ) + 1/L   r = L / (1 + εcosθ)   ここで、 ε=AL :離心率 ・t=0で θ=0, r=r0 とすると   r0 = L / (1 + ε) --> ε = L/r0 - 1 となり、結局(打上げ角0度、すなわち水平方向に発射される)人工衛星の軌道方程式は   r = L / (1 + εcosθ)   ・・・(式3)    ここで、ε = L/r0 - 1        L = (r0v0)2/(GM) となります。この式は極座標での円錐曲線(二次曲線)を表しています。
次のシミュレーションプログラムで確認してください。


(注1)ε=0 すなわち L=r0 のときはr=一定となり、円軌道となります。
    このときの初期速度v0は
     (r0v0)2/(GM) = r0
        より、v0 = (GM/r0)1/2 : 第一宇宙速度です。
       また、ε<1 のときは楕円、ε=1 のときは放物線、ε>1 のときは双曲線となります。

(注2)任意の打上げ角度αに対する軌道方程式は次のとおりです。
     r = L / (1 + εcos(θ+β))
     ここで、L = (r0v0cosα)2/(GM)
         β = tan-1[L/(L-r0)・tanα]
         ε = (L-r0)/r0/cosβ
       
(注3)極座標系での運動方程式については、Wikipediaの「軌道(力学)」を参考にしました。
ホーム