2. 実践―楕円上での位置決定
軌道パラメータの取得
前章で軌道計算の流れを概説しましたが、具体的に各種変数を計算するための情報は、北アメリカ航空宇宙防衛司令部（NORAD）というかっこいい名前の組織のページに載っています。伝統的に衛星の軌道パラメータは2行軌道要素形式（TLE）というフォーマットでまとめられています。
ためしに、NASAが運用している地球観測衛星であるTerraのTLEを見てみたいと思います。あなたは2016/07/01のUTC 17:50:20（現地サマータイム時間10:50:20）あたりにアメリカ・ネバダ州の砂漠（38.50486N, 115.69041W, 標高1435 m）を彷徨っていたとします。TLEは日々更新されているので、space trackなどのページを使って観測時刻に近いTLEを取得しましょう：
1 25994U 99068A 16183.78487350 .00000065 00000-0 24461-4 0 9994
2 25994 98.2080 257.8785 0001614 88.6402 271.4973 14.57112151879645
スペースの個数にも意味があるので勝手に削除しないようにしてください。1行目の3項目がTLEの元期で、2016年の1月1日からかぞえて183.78487350 [day] の意味です。軌道情報はだいたい2行目にあって、軌道傾斜角iが98.2080 [deg]、昇交点赤経Ωが257.8785 [deg]、離心率εが0.0001614、近地点引数ωが88.6402 [deg]です。
楕円上での位置決定
離心率はわかったので、あと離心近点角Eと軌道の長径aがわかれば楕円上での衛星の位置がわかります。まず離心近点角E [rad]ですが、以下の関係を使って計算することができます:
M≡2πt/T [rad]は平均近点角といいます。tは近点通過後の経過時間 [day]、Tは公転周期 [day/rev]です。U軸が近点方向なので、UV平面において、ある時間経過後に「衛星が半径aの円運動をしていたとしたらその位置にいるはずの角度（位相）」とでも解釈できます。楕円軌道上の実際の位置を単に補助円上に投影した離心近点角Eとの意味の違いに注意してください。
元期における平均近点角M0およびその一階時間微分M1・二階時間微分M2はTLEで与えられておりまして、上の例ですとM0 = 271.4973 [deg]、M1 = 14.57112151 [rev/day]、M2 = 0.00000065*2 [rev/day^2]です。角度の単位がバラバラなのと、M2はTLE1行4項目のパラメータに2をかける必要があることに注意してください。というわけで元期からの経過時間をΔt [day] とすれば、
となります。ここでは単位を度（deg）に統一していますがラジアン（rad）でも回転（rev）でもいいと思います。なおEとMの関係ですが非線形のためニュートン・ラフソン法などで数値的に解きましょう。
いっぽうの長径a [km]ですが、ケプラーの第三法則によれば
です。μは地球の重力定数で2.975537×10^15 [km^3/day^2]です。公転周期TはTLEに直接は記載されていませんが、Tとは一回転にかかる日数のことですから、一日あたりの回転数（平均運動と呼ばれます）Mm [rev/day] の逆数をとれば求まります。元期での平均運動がM1ですので、観測時刻での平均運動は
となります。
コーディング例
軌道パラメータは分かっているので、あと必要なのは観測時刻のUTC [day] くらいで、2016/07/01の17:50:20をDOYで表すと183.743287037 [day] です。ではやりましょう。
軌道関係のパラメータがとても多いのでクラスを作ってまとめることにしました。
まあパラメータをまとめるだけなら辞書やリストでもいいんですが、なにかと使いまわしもききますし、一緒にこまごました処理（角度の変換など）も突っ込んでおくほうが楽なので。あとコードにおいてsat['e']とかsat[4]などと書くよりsat.eと書いたほうが眼にうるさくないというのもあります。また、地球関連のパラメータもついでにクラスを作って用意しておきました。重力定数（μ=GM）は先述のとおりですが、回転楕円体（WGS84に準拠）の長径 aE = 6378137 [m]、扁平率 fE = 1/298.257223563も与えておきました。
36-38行目でscipyモジュールを使ってEを計算しています。けっきょく軌道長径 a = 7080.6515 km, 離心近点角 E = 0.93128351 [rad] となりました。地球を原点とする衛星のUVW座標（楕円上での位置）はvect_sat_elpに格納され、
(4224.6226, 5681.4199, 0)
となりました（有効数字は適当に8桁くらいにして表示）。