第4回:関数電卓で行う衛星位置の計算

さて、「関数電卓で行う衛星位置の計算」シリーズは前回で一応の完結をみました。が、前回最後に書いたとおり、せっかく衛星の位置がわかったんだからどの方向に見えるかも計算してみたいと思います。

しかし前回計算した南緯51度, 西経10度というのは日本から見ればほとんど地球の裏側であり、さすがにこれの方位や仰角を計算しても面白くないでしょう。
ということで、前回、前々回の復習も兼ねてもう一度衛星の位置を計算してみましょう。確実にISSが見える位置にいてほしいので、とりあえずJAXAのWebサイトで可視な時間を調べておきました。

軌道要素は以下を使用して下さい。

ISS (ZARYA)             
1 25544U 98067A   16331.50278528  .00003330  00000-0  58332-4 0  9995
2 25544  51.6438 328.6268 0006073 257.1648 241.1942 15.53732614 30231

観測日時：2016/12/04 17:01:30 JST

まずは第1回を参考に各要素を取り出してみて下さい。その後、第2回、第3回の手順のとおりに計算し、ECEF座標を計算して下さい。BLHまで計算する必要はありません。

まずはECIです。答えを書いておくと、僕が計算した結果はX = 4549, Y = -2653, Z = 4280といったあたりです。PreviSatの結果はX = 4543, Y = -2639, Z = 4287といったあたりです。だいたい17kmくらいの差がありますが、誤差の範囲だと考えましょう。

充分に問題のない数字であれば、ECEFまで計算を行って下さい。何百kmもズレているようであればどこかおかしいはずなので、もう一度計算してみて下さい。

参考までに、ECEFはおよそX = -3776, Y = +3671, Z = +4280あたりになります。

観測地点の決定

さて、いよいよ方位と仰角の計算に入ります。が、その前にまずは観測地点が何処かを決めておく必要があります。とりあえず北緯35.71, 東経139.81を使用して計算しようと思います(おおよそ東京スカイツリーのあたりです)。

緯度経度は今後頻繁に使うので、メモリに保存しておきます。といってもA, B, Cは衛星のECEFを保存していますし、x, yは今後使用しますから、D, E, Fのいずれかになります。とりあえずDに緯度、Eに経度を保存しておきましょう。

[35.71][STO][sin]

[139.81][STO][cos]

観測地点の位置ECEFを求め衛星との差をとる(ECEF)

まずは衛星と観測地点の差を計算します。衛星のECEFから観測地点のECEFを引いた差を求めることになります。
本来は観測地点のECEFを計算した後に差を求めるのがわかりやすいのですが、現在自由に使えるメモリが限られているので、観測地点のECEFを計算しつつ、同時に衛星との差も計算していきます。

緯度経度からECEFを求めるにはRec関数を使用します。Rec関数は半径rと角度θからx, y要素を計算するもので、sinとcosを同時に計算する事ができる関数です。sinとcosを個別に計算しても良いのですが、せっかくRec関数があるのでそれを使用します。

まずは緯度からECEFのZを求めます。
[SHIFT][-][6371][SHIFT][)][ALPHA][sin][=]
x = 5173.135203
y = 3718.643997

yがECEFのZとなるので、衛星のECEF、メモリCから引いた差をメモリCに入れておきます。
[ALPHA][x-1][-][ALPHA][S←→D][STO][x-1]
C = 561.235594

次に経度からECEFのXとYを求めます。
[SHIFT][-][ALPHA][)][SHIFT][)][ALPHA][cos][=]

x = -3951.802836
y = 3338.350217

メモリx, yはそれぞれECEFのX, Yです。Zと同じように衛星との差をそれぞれのメモリに入れておきましょう。
[ALPHA][(-)][-][ALPHA][)][STO][(-)]

[ALPHA][°'"][-][ALPHA][S←→D][STO][°'"]
A = 175.3075856
B = 332.5032403

ENUを求める

ENUは適当な日本語訳が見当たりませんでした。近い概念としてはEND, 局地水平座標系がありますが、最後の1文字が違う通り、1軸の解釈が違います。ENUはEast, North, Upの略で、東にいくら、北にいくら、上にいくら、という情報です。対してENDはEast, North, Downの略で、東にいくら、北にいくら、下にいくら、という情報です。
ベクトルの減算と回転で得られる値はENUですから、ここではENUを使います。

さて、ENUの計算ですが、これは「ベクトルの減算と回転」により計算します(他の方法もありますが、関数電卓向きではないのでここでは扱いません)。
ベクトルの減算はすでに先程済ませたので、あとは回転を行うのみです。ベクトルの回転は第3回で行ったECIやECEFの計算と同じです。ただしその時はX軸とZ軸の回転を使いましたが、今回はY軸とZ軸の回転です。

まず最初にZ軸で-Lon分を回転させ、その次にY軸でLat分を回転させます。

LonはメモリEに入っていますが、ZはLonの逆に回転させる必要がありますから、まずはLonを反転させます。
[-][ALPHA][cos][STO][tan]

[ALPHA][(-)][cos][ALPHA][tan][)][-][ALPHA][°'"][sin][ALPHA][tan][STO][)]

[ALPHA][(-)][sin][ALPHA][tan][)][+][ALPHA][°'"][cos][ALPHA][tan][STO][°'"]

[ALPHA][)][STO][(-)]

次にLatの回転を行います。
[ALPHA][x-1][cos][ALPHA][sin][)][-][ALPHA][(-)][sin][ALPHA][sin][STO][)]

[ALPHA][x-1][sin][ALPHA][sin][)][+][ALPHA][(-)][cos][ALPHA][sin][STO][(-)]

[ALPHA][)][STO][x-1]
A = 393.07277, B = -367.13236, C = 408. 636965

以上でENUの計算が終了です。ただし軸の取り方がちょっと直感とは違うと思います。A(X)がUp, B(Y)がEast, C(Z)がNorthの方向になります。

方位と仰角を求める

お疲れ様でした、あとは方位と仰角、ついでに距離の計算を残すのみです。これはECEFからBLHを計算したときとほぼおなじ計算です。

まずは方位を計算します。
[SHIFT][+][ALPHA][x-1][SHIFT][)][ALPHA][°'"][=]
x = 549.3362752
y = -41.93752609

yが方位ですが、負数になっているので360に方位を加算して最終的な方位とします。とりあえずメモリDにでも入れておきましょう。
[360][+][ALPHA][S←→D][STO][sin]
D = 318.0624739

次に仰角を計算します。
[SHIFT][+][ALPHA][)][SHIFT][)][ALPHA][(-)][=]
x = 675.4824545
y = 35.5852838

yが仰角なので、メモリEあたりに入れておきましょう。
[ALPHA][S←→D][STO][cos]

そしてxに入っている数字が彼我の距離となります。ついでなのでメモリFあたりに入れておきましょう。
[ALPHA][)][STO][tan]

ということで、メモリD, E, Fにそれぞれ方位、仰角、距離が入りました。それぞれ318°, 35°, 675kmとなりました。前例に従いPreviSatの値を例に出すと、それぞれ321°, 34°, 688km と言ったところです。方位で3°ほど、仰角で1°ほど、距離で10数kmほどズレていますが、肉眼で観察しようと言う程度であれば数度の誤差は問題ありません。距離に関しては何割もズレていなければ大して問題になりません。

最後に

以上で全4回に渡って続いた「関数電卓で行う衛星位置の計算」シリーズは終了となります。
実際のところ、既知の衛星の位置を電卓で計算してもあんまり使い所はありません。例えばISSが見れる日時を知りたければJAXAのWebサイトで探したほうが楽ですし確実です。衛星がどこにいるかを調べたい場合でも、PreviSat等を使えば簡単に知ることができます。
しかし、「宇宙？人工衛星？なんか難しそう。分かんないや」と思っていても、実は電卓で計算できるんだよ、というところから宇宙を身近に感じる人がいるかな、と思ってこのシリーズを書いてみました。