GPS Test、GPS StatusなどAndroid端末のGPS信号受信状況を表示するアプリがいくつかあります。それらのアプリは現在のGPS衛星の位置を天空図に表示しています。そのためには現在位置からみた各衛星の方位角（Azimuth)と仰角（Elevation)が必要なのですが、これらのデータはアルマナックという衛星の軌道情報から計算で求めることができます。
自分でもやってみようと思って調べながらプログラムを書き、試行錯誤の結果、GPS衛星の位置を求めることができるようになったので、忘れないためのまとめとして書いておきます。
アルマナックというのはGPS衛星の軌道情報データで、衛星のおおよその位置を求めるために使います。今回はQZ-visionで配布しているものを使いました。Web APIでデータを取得することができます。
http://qz-vision.jaxa.jp/USE/ja/almanac
アルマナックの例を以下に示します。
******** Week 739 almanac for PRN-01 ******** ID: 01 Health: 003 Eccentricity: 2.453327179E-03 Time of Applicability(s): 147456.0000 Orbital Inclination(rad): 0.9602863543 Rate of Right Ascen(r/s): -8.0117622935E-09 SQRT(A) (m 1/2): 5153.626465 Right Ascen at Week(rad): 2.189035151E+00 Argument of Perigee(rad): 0.416642440 Mean Anom(rad): 2.076584921E+00 Af0(s): 8.678436279E-05 Af1(s/s): 3.637978807E-12 week: 739
おおよその手順は以下のようになります。
1．現在時刻に衛星が楕円軌道のどの位置にあるかを地球を原点とする楕円平面の2次元座標系で求める
2．もとめた座標を地球の中心を原点とする3次元の座標系（地球中心地球固定座標系：ECEF）に変換する
3．さらに方位角と仰角を求めたい地球上の場所を原点とする地平直交座標（ENU座標）に変換する
衛星軌道上のどの位置にあるかを求める
GPS衛星は楕円軌道を描いて地球の周りを回っています。地球は楕円の中心ではなく焦点の位置にあります。
図1
楕円は以下の式で表されます。
…………(1)
ここで、a=長半径、b=短半径です。アルマナックには aの平方根(Square Root of Semi-Major Axis)が含まれています。地球は楕円の中心からa * eだけ離れています。e は離心率といいます。離心率eは以下の式で表されます。
…………(2)
まずはGPS衛星の楕円軌道上の位置を2次元で求めるわけですが、これに必要なパラメータや関係式がInterface Control Document という文書のTable 20-IV. Elements of Coordinate Systems にのっています。基本的にはここに書かれているとおりにやります。以下は現在時刻の衛星位置を求める手順です。
図1の平面上での地球を原点とした衛星のXY座標を求めるためには、図のVk（真近点角）を求めればよいのですが、そのためにはEｋ（離心近点角）を求めることが必要です。つまり、衛星の楕円軌道上の地球を原点とした2次元座標は
…………(3)
となります。Rkは地球と衛星との間の長さで、以下で計算できます。
…………(4)
ここでVkは離心近点角Ekを用いると以下のように書けます。
…………(5)
したがって、
…………(6)
となります。
Ekはどうやって求めるかというとここでケプラーの方程式が登場します。
…………(7)
Mkは平均近点角といってアルマナックに含まれる情報から計算できます。
…………(8)
M0はアルマナックに含まれる基準時刻における平均近点角(Mean Anomary)、Tkは現在時刻、N0は平均運動といって以下の式で求まります。GMは定数で 3.986005E+14です。
…………(9)
現在時刻Tkは、アルマナックデータの基準時刻(Time of Applicability)からの経過秒でなければならないので、GPS週の値とTime of Applicabilityを使って計算します。Time of Applicabilityは現在のGPS週の始まり（日曜日の午前0時）からの経過秒なので、以下の式で求めることができます。
…………(10)
Tgps は現在時刻TkをGPSエポック(1980年1月6日0時0分0秒）からの経過秒で表した値です。また、Toeはアルマナックの基準時刻(Time of Applicability)です。week はGPS週番号です。GPS週番号に1024足しているのは、1999年8月にGPS週が1023を越えて一度リセットされているためです。
以上から(8)の式を解いてEkを求めます。(8)はニュートン法を使って解きます。
…………(11)
とおいて、
…………(12)
にあてはめると、
…………(13)
となるので、初期値 E0 を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座標を求めます。
ここのところがわかりにくいので、別記事を書きました。
まず、Vkにアルマナックの近地点引数(Argument Of Perigee)ωの値を加えた角度Pkを求めます。
…………(14)
(3) で Vkの代わりにPkを使ったX’,Y’を求めます。
次に、軌道傾斜角 i (Inclination)だけX軸まわりに回転させ、さらにΩkだけZ軸まわりに回転させればOKです。i はアルマナックに含まれていますが、Ωkは以下の式で計算します。
…………(15)
Ωeは地球の自転速度で7.2921151467E-05です。Ωdは基準時刻における昇交点黄経(Longitude of Ascending Node of Orbit Plane at Weekly Epoch)です。 Ω0は基準時刻における昇交点経度変化率(Rate of Right Ascention)です。Toeは基準時刻（Time of Applicability)です。
以上から、ECEF座標系のX,Y,Zは以下の回転行列で求めることができます。
…………(16)
これを展開すると以下のようになります。
…………(17)
ECEF座標をENU座標に変換する
最後にENU座標に変換します。ENU(East, North, Up)座標は地表面のある点を原点として天頂方向をZ軸、東方向をX軸、北方向をY軸とした座標系です。
ECEFからENUへの変換は座標系の回転と原点の移動でできます。ENU座標の原点にする地表面の位置の経度をλ度、緯度をφ度とすると、
・Z軸中心にλ度回転する（Z軸の正方向に向かって右ねじの方向に回す。グリニッジ子午線のところにあった基準子午線が原点の子午線のところになる）
・Y軸中心に(90度ｰφ度)回転する（Z軸が原点の天頂方向を向くようになる）
この時点でY軸が東、X軸が南を向いているので、
・Z軸中心に90度回転する
・原点を移動させる
式で書くと以下のようになります。
…………(18)
ここまでできればあとはAzimuth(方位角:北から時計回りの角度）とElevation(仰角）は簡単な計算で求めることができます。
自分で書いたプログラムで以上の計算をしてみて、QZレーダーという衛星位置表示アプリケーションと同じ数値を得ることができましたので、間違いはないと思います。（上記の説明に間違いはあるかもしれませんが）
http://qz-vision.jaxa.jp/USE/ja/qz_radar
アルマナックから衛星の位置を求める方法を調べるにあたっては、以下のサイトを何度も何度も参照しました。
