いよいよ大ずめです。前回「アルマナックデータによるGPS衛星の位置計算其の2」で離心近点角Eをニュートン法により算出しました。今回は、クライマックスです。

さてさて、ここまでくれば、後は気が楽です。

まず、前に「航法メッセージによる GPS衛星の位置の求め方 其の3」に記載した、平面軌道上のGPS衛星位置(x',y')を求める公式を利用して、2次元面での、GPS衛星位置を求めます。

これはもう簡単ですよね。

x' = a*cos(E) - a*e
y' = a*sqrt(1-e²)*sine(E)

これで算出できます。ここで、お分かりだとは思いますが念のため、aは長半径、eは離心率、Eは離心近点角になります。
あとは、アルマナックデータおよび、前回求めた離心近点角Eを当てはめて解くだけです。解くと以下のようになります。

x' = -11120449.437278792
y' = -24439569.638433773


さてさて、最後に、上記2次元座標を3次元座標に変換してやりますがこれも、「航法メッセージによる GPS衛星の位置の求め方 其の3」で記述した、回転行列を利用します。

まずは、行列を解きます。すると、

X     =  x' * cos(Ωk) * cos(ω) - x' * sin(Ωk) * sin(ω) * cos(i)
           - y' * cos(Ωk) * sin(ω) - y' * sin(Ωk) * cos(ω) * cos(i)

Y     =  x' * sin(Ωk) * cos(ω) + x' * cos(Ωk) * sin(ω) * cos(i)
      　　- y' * sin(Ωk) * sin(ω) + y' * cos(Ωk) * cos(ω) * cos(i)

Z     =  x' * sin(Ωk) * sin(i)　+ y' * cos(Ωk) * sin(i)


となります。少々見にくいですが、こんな感じになります。ここでΩkは前回求めた、現在時刻による昇交点経度、i、ωはどちらもアルマナックデータより取得できそれぞれ、軌道傾斜角、近地点引数になります。考え方のおさらいをしますと。

「3次元の軌道面は、2次元の軌道面をz 軸回りにω 回転させ，続いてx 軸まわりにi，z 軸回りにΩ 回転させたものである」

でしたよね。

さて、実際にときますと、以下のようになります。

X = -26312458.28576335
Y = -3442369.7618006254
Z = 4094091.9850355694

どうでしたか？とうとう3次元空間における、GPS衛星の位置情報が求まりました。これは紛れも無い解です。しかしながら、GPS受信機という意味では不十分です。といいますのは、航法メッセージの中にはアルマナックデータだけでなく、エフェメリスデータというものが入っており、こちらのほうが正確な情報を持っているのですが、残念ながら、このやり方では、エフェメリスデータに含まれる補正データがうまく扱えません。よって、新たな角度、真近点角を利用した方法に置き換える必要があります。え！じゃー最初からそれをおしえろよ！と言われるかもしれませんが、概念的には今までの内容を理解していれば、問題ありませんし、入門としては、今回の方法から入ったほうが分かりやすいと思います。

機会があれば、エフェメリスデータによるGPS衛星の3次元空間上での位置の求め方を公開したいと思います。

[0回]