2016年12月5日月曜日

Three.jsの軌道表示

とりあえず暫定的にThree.jsで軌道を表示するやつをJS Orbitの下に置いときました。
http://www13.plala.or.jp/rian/jsOrbit/three/?autosize

軸とかラベルとか前部表示だとかなり見づらいです。とは言え不要な情報というのはあまりないので、せいぜいECEFのラベルを消すくらいでしょうか。

フォントの関係でラベルはすべて英語です。間違ってても文句言わないで下さい。そもそも軌道の表示自体が間違ってるかもしれませんが(そのためのページタイトル「ThreeJSの練習」です?)。

黒いところでマウスを動かすと視点を動かすことができるので、見やすい視点を探してみて下さい。画面をクリックすればその位置で視点を固定し、もう一度クリックすれば固定を解除します。

視線固定を解除していて、カーソルがブラウザの外にあるときは自動的に視線を回転させます。フルスクリーンで勝手に回しておけばちょっとした癒やし効果がありそうです(個人の感想です)。


軽く見方を説明しておきます。

まず赤、緑、青の軸はそれぞれX, Y, Zの向きを示しています。

ECIはX+が春分点を向き、ECEIはX+が経度0度を向いています。ECI, ECEFはZ+が北極方向、Yはそれらに直行する軸です。

OrbSという軸はOrbit Surfaceの略で、軌道面になります。

中くらいの白い円・楕円は衛星軌道を示しています。

外側の大きなオレンジの円弧は軌道面に至る回転を示しています。Ω Right Ascensionは昇交点を、i Inclinationは軌道傾斜角を、ω Argument of Perigeeは近点角を示しています。



Three.jsをダウンロードしてきたのがおよそ3日前ですから、初めて触ってからたった3日くらいでこれくらいのモノが作れるようになります。なれてれば半日もかからないと思います。Webブラウザでマルチプラットフォームで動くものがこんなに簡単に作れるんですからすごい時代ですね(ウチのモバイル端末はWebGLに未対応なのでマルチプラットフォームの動作確認はしてないですけども)。

Three.jsはObject3Dで入れ子にできますが、JavaScriptでこれを書くのはかなり面倒です。僕はあんまりGUIツールは好きじゃないのですが、このあたりは木構造で簡単にThree.jsのコードを生成できるようなソフトウェアがあればすごい便利だと思います。そういうのがあれば今回のようなテクスチャのないモノなら10分や20分で作れるだろうなぁ。


とりあえず、後は緯度経度の矢印を作ったりとか、そのあたりでしょうか。
本当は昇交点だったり傾斜角だったりを少しずつ増やしていくアニメーションとかも作ってみたいんですが、リアルタイムに書き換えるのをどうすれば良いのかわからないのでもうちょっとThree.jsを調べてからになると思います。メモリリーク怖い。

2016年12月4日日曜日

軌道面



Three.jsで作ってる衛星軌道の概要を表示するやつ、とりあえずECI, ECEF, Orbit Surfaceの情報を表示できるようになってきました。ということで経過報告。といってもECIとECEFはジャマなので消してますが。

青の矢印が地球の中心(ECI, ECEF, 軌道面すべての原点)です。その周りの円と、長半径が同じ楕円は地球の大きさを示しています。
その外側に有る円と楕円は軌道の大きさです。楕円が実際の衛星軌道、外側の円は衛星軌道に外接する真円です。

赤が衛星軌道に外接する真円の原点で、緑の矢印は衛星の位置によって移動しますが、緑の矢印の部分にあるオレンジと赤の線は常に直角です。
水色の矢印で示した、オレンジの線と楕円の交点に衛星がいます(ラベルの位置が不適切なので真円のところに衛星がいるように見えてしまいますが)。


他の言語で直接3Dを作るよりはThree.jsを使ったほうが簡単なんでしょうけども、それでも結構面倒です。特にpositionやquaternionが思い通りに設定できないのが面倒。いちいちObject3Dで入れ子にしたりとかしてて、なんかもうちょっとうまくできないものかなぁと。

メモ:離心率から短半径と超半径の比率を計算

楕円のパラメータには超半径a, 短半径b, 焦点c, 離心率e, 扁平率f等がある(らしい)。abからcやeやfを計算する式はちょっと探せばいくらでも出て来るが、eから計算する方法はあんまり見つからなかった。そういう計算を必要とする人がいないからどこにも書いてないのかもしれないけど、せっかく苦労して調べたのでメモしておく。


たったこれだけの計算式を調べるのにだいぶ時間を使ってしまった。こういうところで色々不足してるのが思い知らされる。

2行軌道要素だと離心率が与えられるし、長半径を計算するのも簡単だけど、短半径がないと楕円を作れない、というところでこういう計算が必要になる。

/*
ぜんぜん関係ないけど、いまディスカバリーの怪しい伝説でテレビゲームを扱ってる。フルーツニンジャとかDOOMとか。このDOOMみたいなのやったら楽しいだろうなぁ。夏休みとか人集まるような時期にお化け屋敷みたいなふうにやってみたい。「超人的な狙撃手」でやってた、紙で作った部屋でやるサバゲも楽しそうだよね。ちゃんとした壁は邪魔者でしかないから、どっかのだだっ広いところに垂木とかでフレームを作って、模造紙で壁を作る、そんな感じでやるならエアガンを撃てる部屋を探す必要もないし、広いことについては定評のある北海道だから、壁の有る部屋を作るのは簡単だろう。ただ天井もないから、雨とか降ったら最悪なことになる。そうでなくても前日に雨が降るだけでも問題。ほんとは単管とかで2階フロアとか作って、ラペリングとかもやってみたいけど、安全確保が大変だからお化け屋敷みたいな人を集めるイベントでは無理だろうなぁ。まぁこのあたりは準備中にスタッフだけで遊ぶとかもできるだろうし。やるにしたって7月か8月になるだろうから、もうちょっとじっくり考えてみよう。
*/

2016年12月3日土曜日

衛星の軌道面とか軸の話

文字だけで軌道の座標とか説明するのはかなりつらいので、試しにCADで書いてみた。


まだ無理があるかな、という感じ。そもそも3次元を2次元画像に押し込むのが無理がある。
マウスで視点をグリグリ動かすとなんとなく3次元形状を把握できるようになるので、試しにThree.jsで表示してみた。


「5分でできるThree.js入門」とか言っておきながらThree.js本体をダウンロードするのに30分かかるじゃねぇか!!とか毒づきながらも、丸1日でここまでできるんだから楽なもんだなぁとも思ったり。


今回は離心率が約0.49、平均運動が7.4くらい、軌道傾斜角が63くらいの衛星を表示してみた。
勝手な思い込みでQZSは結構離心率が大きい印象だったけど、実際には0.075くらいしかない。これくらい小さい離心率でもあんな動きになるんだねぇ。


最終的にはTLEを渡したら表示できるような感じになると面白いかなーとか考えながらも。もうちょっと色々表示を追加したい。いくらライブラリにまとまって簡単とはいえ、やはりそのままだと文字の表示とかが大変なのでこのあたりはうまくラッパー関数を作ってやる必要があるのかな。

2016年11月30日水曜日

第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等を使えば簡単に知ることができます。
しかし、「宇宙?人工衛星?なんか難しそう。分かんないや」と思っていても、実は電卓で計算できるんだよ、というところから宇宙を身近に感じる人がいるかな、と思ってこのシリーズを書いてみました。

2016年11月27日日曜日

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

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

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


2016年11月25日金曜日

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

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

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


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

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

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


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