2016年11月27日日曜日

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

前回に続き衛星の位置計算です。

前回は軌道面上の2次元位置を計算しました。今回はBLH(緯度経度高度)の計算までを行います。




1) 近地点引数ωと昇交点赤経Ωを計算する

まずは軌道面上の2次元直交座標をECI(地球中心慣性座標系)に変換しますが、それには近地点引数ωと昇交点赤経Ωが必要になります。これらはTLEで与えられていますが、値は変化するため、現在の値を計算する必要があります。


[238.1568][+][(][(][31.32][(][2][-][2.5][sin][51.6454][)][x2][)][)][÷][(][SHIFT][x10x][(][ALPHA][sin][÷][6371][)][x■][3.5](→)[)][)][ALPHA][S←→D][STO][cos]

[353.2810][-][(][(][31.32][cos][51.6454][)][)][÷][(][SHIFT][x10x][(][ALPHA][sin][÷][6371][)][x■][3.5](→)[)][)][ALPHA][S←→D][STO][sin]

E = 241.797566
D = 348.3976662

2) 3次元位置を計算する(ECI)

いよいよ一番面倒なところに入ってきます。ベクトルの回転です。まずECIで3回、ECEFで1回、計4回の回転を行います。

とりあえず計算式です。


一応X, Y, Zの3種類ですが、今回の計算で使用するのはXとZの2つのみです。

最初にZ軸に近地点引数ωで回転させます。次にX軸に軌道傾斜角iで回転させます。最後に再びZ軸で昇交点赤経Ωで回転させます。

前回、メモリAにX軸の、メモリBにY軸の結果が入っています。
またメモリCはZ軸を割り当てる予定ですが、現在メモリCには離心近点角Eが入っています。ということでまずはこれを0クリアします。
[0][STO][x-1]


最初に近地点引数ωでZ周りの回転を行います。
[ALPHA][(-)][cos][ALPHA][cos][)][-][ALPHA][°'"][sin][ALPHA][cos][STO][)]

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

[ALPHA][)][STO][(-)]
まずA(X)を計算しますが、計算前のA(X)は次の計算で必用ですから、一旦別のメモリにストアします。今回はxを使用しました。次にメモリB(Y)を計算します。これは直接メモリB(Y)にストアします。最後にxに保存したX軸をメモリAに戻します。現在の値はA = -1015.1445, B = -6703.7062, C = 0 となっています。


次に軌道傾斜角iでX周りの回転を行います。ただ4回も軌道傾斜角を入力するのは面倒なので、メモリEに軌道傾斜角を設定して、その後で回転を行います。
[51.6454][STO][cos]

[ALPHA][°'"][cos][ALPHA][cos][)][-][ALPHA][x-1][sin][ALPHA][cos][STO][)]

[ALPHA][°'"][sin][ALPHA][cos][)][+][ALPHA][x-1][cos][ALPHA][cos][STO][x-1]

[ALPHA][)][STO][°'"]
現在の値はA = -1015.1445, B = -4159.828, C = -5256.9485です。


最後に昇交点赤経ΩでZ周りの回転を行います。
[ALPHA][(-)][cos][ALPHA][sin][)][-][ALPHA][°'"][sin][ALPHA][sin][STO][)]

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

[ALPHA][)][STO][(-)]
現在の値はA = -1831.0177, B = -3870.6668, C = -5256.9485です。


お疲れ様でした。やっとメモリにECI形式の3次元座標が入りました。

さて、ちょっと休憩です。とりあえず現在のところ精度がどれくらいなのか確認してみましょう。
前回のTLEをPreviSatに読み込ませてECIの座標を確認してみます。位置はX = -1829.828, Y = -3875.740, Z = -5259.548 となりました。直交座標ではわかりづらいので、距離を計算してみました。結果、PreviSatの位置と関数電卓で計算した位置のズレは5.82km程でした。ECIの原点から衛星の位置までは6785kmですから、ズレは0.1%にも満たない程度です。意外と精度よく計算できるでしょ?

3) グリニッジ恒星時を計算する

ECIからECEFへ変換するにはグリニッジ恒星時の計算が必要になります。


年Y, 月M, 日D, 時h, 分m, 秒sですが、1月と2月は年-1, 月+12を使用します。また日本標準時で計算している場合はTJDから9/24を引きます。

[ALPHA][+][365.25][×][2016][)][+][ALPHA][+][2016][÷][400][)][-][ALPHA][+][2016][÷][100][)][+][ALPHA][+][30.59][(][11][-][2][)][)][+][22][+][1721088.5][+][21][÷][24][+][53][÷][1440][+][15][÷][86400][-][2440000.5][-][9][÷][24][STO][sin]

D = 17714.53698


[1.0027379094][ALPHA][sin][+][0.671262][STO][sin]

D = 17763.70904

得られたGはrev単位なので、小数部に-360をかけて度単位にします。

[-360][(][ALPHA][sin][-][ALPHA][+][ALPHA][sin][)][)][STO][sin]

D = -255.2538498

4) 3次元位置を計算する(ECEF)

次にZ周りの回転を行います。

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

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

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

A = 4209.240625, B = -785.481799, C = -5256.9485

これで終了です。いままでに比べれば遥かに簡単ですね。
ということでこれでECEF(地球中心地球固定座標系)になりました。「地球固定」というくらいですから、地球に同期した座標となります。ですがほとんどの人が使い慣れているのはBLH, 緯度経度高度でしょう。ということで次はBLHに変換します。

5) 3次元位置を計算する(BLH)

ECEFから大雑把な緯度経度を計算するのは比較的簡単です(正確な緯度経度高度を計算するのは難解です)。
一般的に大雑把な緯度経度を計算するにはatan(tan-1)関数を使用します。しかしこれは象限を考えたりと面倒なので、関数電卓のpol関数を使用します。三角関数で言うとAtan2となりますが、polとAtan2ではxとyを逆に与える場合があるので注意してください。

[SHIFT][+][ALPHA][(-)][SHIFT][)][ALPHA][°'"][=]

x = 4281.902416
y = -10.57032552

これによりメモリxにXYの長さが、メモリyに経度が入ります。
次の計算でyが上書きされてしまいますから、とりあえずメモリFに退避させておきましょう。

[ALPHA][S←→D][STO][tan]

次に緯度を計算します。

[SHIFT][+][ALPHA][)][SHIFT][)][ALPHA][x-1][=]

x = 6780.132491
y = -50.83641222

これによりxにXYZの長さが、メモリyに緯度が入ります。
xに入っている値はECEFのsqrt(X^2+Y^2+Z^2)と同じ値です。つまりECEF原点からの距離です。ということで地表からの距離を計算しておきましょう。

[ALPHA][)][-][6371][STO][)]

x = 409.132491


これで大雑把な緯度, 経度, 高度の値が出ました。現在、衛星は南緯50.836412度, 西経10.570325度, 高度409.132491kmにいるようです。
先ほどと同じようにPreviSatで計算した緯度経度高度と比較してみましょう。
PreviSatで得た緯度経度高度は次のとおりです
 緯度:南緯51.00度 経度:西経10.53度
 高度:419.5km

差は緯度で0.16度、経度で0.04度、高度で10.37kmとなりました。緯度経度はかなり正確に出ているとわかります。高度は誤差がすこし大きいですが、これは地球の断面が楕円に対し、今回の計算では真円と仮定して計算していることが原因として考えられます。WGS84から得た緯度51度付近での地球の直径はおよそ6361.5kmです。一方今回地球の直径として使用した値は6371kmです。その差は9.5kmであり、誤差の10.37kmとかなり近い値です。この分を差し引くと、誤差は0.9kmほどになります。
どちらにしろ、関数電卓で人間が計算した割には良い精度だと思いませんか?


予定では第3回の今回で終了でしたが、もう1回続けようと思っています。次回は2点間の位置から方位や仰角を計算してみようと思っています。せっかく衛星の位置を計算したんですから、どの方向に見えるかも計算したいですよね。

0 件のコメント:

コメントを投稿