離心率eは以下の式で表されます。基本的にはここに書かれているとおりにやります。以下は現在時刻の衛星位置を求める手順です。図1の平面上での地球を原点とした衛星のXY座標を求めるためには、図のVk(真近点角)を求めればよいのですが、そのためにはEk(離心近点角)を求めることが必要です。衛星の楕円軌道上の地球を原点とした2次元座標は............(3)となります。Rkは地球と衛星との間の長さで、以下で計算できます。............(4)ここでVkは離心近点角Ekを用いると以下のように書けます。............(5)したがって、............(6)となります。Ekはどうやって求めるかというとここでケプラーの方程式が登場します。............(7)Mkは平均近点角といってアルマナックに含まれる情報から計算できます。Anomary)、Tkは現在時刻、N0は平均運動といって以下の式で求まります。3.986005E+14です。Applicabilityを使って計算します。Applicabilityは現在のGPS週の始まり(日曜日の午前0時)からの経過秒なので、以下の式で求めることができます。現在時刻TkをGPSエポック(1980年1月6日0時0分0秒)からの経過秒で表した値です。Applicability)です。GPS週番号に1024足しているのは、1999年8月にGPS週が1023を越えて一度リセットされているためです。以上から(8)の式を解いてEkを求めます。(8)はニュートン法を使って解きます。Mkとして上式を繰り返し計算しある程度Eの値が収束したところで、その値をEkの値とします。(4)からVkを求めて(3)から座標を求めれば楕円軌道平面上の2次元座標系における座標値X’,Y’を求めることができます。楕円軌道座標をECEF座標に変換する次に求めた楕円軌道面上の2次元座標を地球中心の3次元座標(ECEF:Earth-Centered,Earth-Fixed)に変換します。ECEF座標系は地球の中心を原点とし、北極方向をZ軸、グリニッジ子午線方向をX軸とする3次元の右手座標系です。衛星の楕円軌道は天球に対して固定であるのに対して、ECEF座標系は地球の自転といっしょに回転しています。楕円軌道面とECEFの赤道面(XY平面)との関係から座標系の回転を行ってECEF座標を求めます。ここのところがわかりにくいので、別記事を書きました。Perigee)ωの値を加えた角度Pkを求めます。Vkの代わりにPkを使ったX’,Y’を求めます。(Inclination)だけX軸まわりに回転させ、さらにΩkだけZ軸まわりに回転させればOKです。アルマナックに含まれていますが、Ωkは以下の式で計算します。............(15)Ωeは地球の自転速度で7.2921151467E-05です。Ascention)です。Applicability)です。以上から、ECEF座標系のX,Y,Zは以下の回転行列で求めることができます。............(16)これを展開すると以下のようになります。............(17)ECEF座標をENU座標に変換する最後にENU座標に変換します。Up)座標は地表面のある点を原点として天頂方向をZ軸、東方向をX軸、北方向をY軸とした座標系です。ECEFからENUへの変換は座標系の回転と原点の移動でできます。・Z軸中心にλ度回転する(Z軸の正方向に向かって右ねじの方向に回す。・Z軸中心に90度回転する・原点を移動させる式で書くと以下のようになります。............(18)ここまでできればあとはAzimuth(方位角:北から時計回りの角度)とElevation(仰角)は簡単な計算で求めることができます。自分で書いたプログラムで以上の計算をしてみて、QZレーダーという衛星位置表示アプリケーションと同じ数値を得ることができましたので、間違いはないと思います。(上記の説明に間違いはあるかもしれませんが)http://qz-vision.jaxa.jp/USE/ja/qz_radarアルマナックから衛星の位置を求める方法を調べるにあたっては、以下のサイトを何度も何度も参照しました。