2016年11月25日金曜日

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

さて、前回に引き続き関数電卓で行う衛星位置の計算です。

前回、「次回はECIまで計算するよ」と書きましたが、ECIの前の段階までに変更となりました。ということでECI(からBLHまで)は次回となります。今回は軌道面上での2次元位置を計算するところで終了となります。


まず、その前に、いくつかの前提を説明しておきます。

1. 使用する電卓
使用する電卓はもちろん関数電卓ですが、今回はCASIO fx-JP900というモデルを使用します。変数が最大10個ほど(A-Fとx,yにMとAns, preAnsが)使用できます。
他の電卓には無い機能をつかったりしているかもしれません。その点はご了承下さい。
以降、このシリーズではキーの操作はすべてこの電卓のものです。

2. 精度
今回、衛星の位置は直交座標時の距離で40km、緯度経度時の角度で5度程度の誤差が発生するはずです。天体望遠鏡で撮影したい場合やレーザー通信を行いたい場合は問題が有る精度ですが、そもそもそんなことには使わないと思うので、問題ないこととします。


ということで、実際の計算に入ろうと思います。必用な値は前回書いておきましたので、別ウインドウでその数字をいつでも見れるようにしておくと便利だと思います。




0) 下準備

以降、電卓のメモリをガンガン上書きしていきますから、必用な数字が有るならそのあたりは各自うまくどうにかして下さい。
また紛らわしくなるのですべての変数を0クリアしておくと良いでしょう。

それから、メモリEにBSTARを記録しておきます。BSTARは結構あちこちで使うので一々入力するのは大変なためです。

[0.53814][SHIFT][log][-4][STO][cos]
入力式に"→E"が増え、計算結果に"5.3814x10-5"のような表示になれば正常です。
[SHIFT][log]は10を基数とした指数を入力するコマンド、[STO][cos]は変数Eにストアする、というコマンドです。本来であれば[SHIFT][10■]や[STO][E]のようになりますが、この電卓の別機能表示は見づらいので、標準機能の文字で書いていきます。

1) 現在の時間を調べる

前回、観測日時は2016年11月22日 21時53分15秒(JST)だよ、と書きました。が、EpochTimeに月時分秒が無いように、観測日時も年と日に変換したいと思います。

まず月日時分秒を1月1日からの経過日数に変換します。
1月:31日まで、2月:28日まで、3月:31日まで、4月:30日までというように、前月までの日数をすべて加算していきます。
1月から10月までをすべて加算すると304日ですが、今年はうるう年ですから、2月が29日までなのでそれを含めて305日です。
さらに前日までの日数22日間を加算し、327日となりました。

さらに、これに時, 分, 秒を加算します。時/24+分/1440+秒/86400を加え、327.9119792日となりました。ただこの計算の元になった時間はJST(日本標準時)ですから、UT(世界標準時)に変換する必要があります。ということでJSTの時差、9/24を引いて327.5369792日となります。

ただし、衛星の計算に必用な日時は、Epochからの経過日数ですから、ここからEpochの日数を引きます。最終的には0.9835008567日となりました。
この数字は頻繁に使うので、とりあえずメモリyに記録しておきます。
またこの0.98-日のことはΔtと呼ばれます。

2) 軌道長半径αを計算する

軌道長半径というのは軌道の楕円に外接する円の半径です。ISSはほぼ円軌道ですから、ほぼ軌道半径に相当します。

軌道長半径を計算するには平均運動Mが必用ですが、平均運動は絶えず変化するため、それを補正する必要があります。そしてその後で実際の軌道長半径を計算します。


まずMmをメモリDに登録し、その後αをメモリDに登録します。Ansメモリを使っても良いのですが、αの計算で入力ミスをするとまたMmの計算からやり直しになってしまうため、Ans, PreAnsはあまり使わないことをおすすめします。

M2とΔtはメモリにありますから、それを利用してMmを計算しましょう。
[15.53693589][+][ALPHA][cos][ALPHA][S←→D][STO][sin]
[(][2.975537][x10x][15][÷][(][4][SHIFT][x10x][x2][ALPHA][sin][x2][)][)][x■][1][÷][3][STO][sin]

Mm = 15.53698882, α = 6784.07495となります。

ちなみにαから地球半径を引くと406kmとなり、もっともらしいISSの高度となります。が、前述したようにαは一方の焦点を地球中心とする楕円に外接する円の半径ですから、原点は地球の中心ではありません。充分に離心率の小さい軌道であれば高度として使えますが、基本的にα-r=高度という計算は行うべきではありません。

3) 平均近点角Mを計算する




[278.5578][÷][360][+][15.53693589][ALPHA][S←→D][+][(][ALPHA][cos][÷][2][)][ALPHA][S←→D][x2][STO][°'"]
[360][(][ALPHA][°'"][-][ALPHA][+][ALPHA][°'"][)][)][STO][°'"]

M = 16.05438745 (rev)
M = 19.57948232 (deg)

Mの結果はメモリBに登録しておきます。ただ、Mはrevとして結果が出ていますから、一旦小数部を取り出して360をかけることにより度単位に変換します。本来関数電卓には余りを計算する機能がありますが、少なくともfx-JP900ではまともに使えないので、このようにint(小数部を切り捨てて整数化する)関数を使用しています。

4) 離心近点角Eを計算する

離心近点角Eを計算するにはニュートン・ラフソン法という方法を使うのが一般的です。関数電卓にはそのためのソルブ機能というものもあります。ただ、ソルブ機能を使ったことがない人に文字だけで使い方を説明するのは非常に難解なため、今回は人間が1回1回計算することにします。念のため、ソルブ機能を使える人のために書いておくと、Eは以下の式で求めることができます。

さて、手動で計算する方法です。
まずメモリBをメモリCにコピーします。
[ALPHA][°'"][STO][x-1]
次に以下の式を入力します
[ALPHA][°'"][+][0.0006168][sin][ALPHA][x-1][STO][x-1]
その後何回か[=]を押して計算を繰り返します。5回程度繰り返せば十分です。

E = 19.57968902

手動での繰り返し計算とソルブでは答えに若干の違いがあります。0.1未満程度のはずなので、気にしないで次に進みましょう。

5) 2次元位置を計算する

いよいよ位置情報の計算に進みます。といってもまずは2次元位置、つまり軌道面上での直交座標です。


[ALPHA][sin][cos][ALPHA][x-1][)][-][0.0006168][ALPHA][sin][STO][(-)]
[ALPHA][sin][√■][1][-][0.0006168][x2](→)[sin][ALPHA][x-1][STO][°'"]

U = 6387.610279
V = 2273.462452


前述したとおり、今回は軌道面上での2次元座標を計算したところで終了です。現在メモリA, Bに入っている値が軌道面上でのX, Y座標となります。ということで、今回はここまで。

さて、次回予告です。
次はECIへの変換、ECEFへの変換、そしてBLHへの変換を行います。



0 件のコメント:

コメントを投稿