3. 実践―座標変換

ひたすら行列計算

前章で得られたUVW座標のECEF座標への変換は結局

となるのでした。θGは、観測時刻における、春分点からみたグリニッジ子午線までの時差（グリニッジ恒星時という）を調べる必要があります。理科年表などで基準となる過去のグリニッジ恒星時を調べ、そこからの地球の自転回数を計算し…とやれば求まるのですが、面倒なので今回は直接計算してくれるサイトを使ってみました。

それによると時差は12:31:10.504、角度にして187.79377 [deg] になります。

vect_sat_3dが衛星のECEF座標で、

(-1315.5624, -5457.8393, 4313.6127)

となりました。さてECEF座標と緯度経度との関係は

でした。これを逆にとくと、

となります。Φ'が地心緯度、Φが地理緯度です。南半球（z'<0）の場合には地心緯度・地理緯度ともに負符号にする必要がありますので、そこだけ注意です。今回は北半球なので大丈夫。

結果は

sat_lon:[-103.55211004]

sat_lat:[37.72282564]

sat_alt:[711.57534636]

と出ました。というわけでTerra衛星は37.72283N, 103.55211Wの位置におり、楕円体高は711 kmのようです。周回衛星の高度として妥当な感じがします。

観測者の気持ちになる

観測者の緯度経度（38.50486N, 115.69041W）および楕円体高を与えてECEF座標に変換します。楕円体高のかわりに観測値の標高（1435 m）を与えてもほとんど変わらないとは思いますが、正しくは観測値のジオイド高を足す必要があります。ジオイド高はどこかからモデルを拾ってきて調べてもいいのですが、例によって適当なサイトがあったので調べてもらうと、-23 mくらいのようでした。つまりHs = 1435-23=1412 [m] です。では。

今回はECEFへの座標変換までをクラスに突っ込んでみました。

なぜ衛星の座標変換のときもそうしなかったのかについてですが、気分の問題です。「本筋と関係ない瑣末な処理」と感じたときクラスに突っ込んだり関数化したりしてmain部分から不可視化すると読みやすいような気がする、という程度です。ソフト開発をしているわけでもないですし、自分の好きに書けばよいと思います。

最後に、

によってENU座標に変換すれば完了です。

結果は

sat_azimuth:[91.0089892]

sat_zenith:[62.67530722]

方位角91°、天頂角63°になりました。東の空に衛星が飛んでいたのですね。

全部の処理をまとめると以下のような感じです：

おまけ

ところでこれは内緒なんですが、これまでの日々のことは一旦忘れて

などとしますとよくわからない数値が出てきます。無事に車輪が再発明できたみたいで何よりです。