NASAの公開データに触れて宇宙を描こう。〜SPICEの解説〜

これはなに?

この記事はeeic Advent Calendar2017の23日目として書かれたものです。
qiita.com

ソースコードやデータを公開する動きは世界的に強くなっています。 実は航空宇宙分野にも公開されているデータは数多く存在します。 今回はその中の一つであるSPICEについて解説を行い、手軽に軌道を描いてみようと思います。

僕が軌道をいい感じに描いて楽しいって言っているだけの記事になってしまうかもしれません。すみません。

扱うミッションは

です。 Cassini Grand Finaleから読んでも大丈夫なようにしているので、 その前は読み飛ばしていただいても構わないです。

描画したものをみてるだけでも楽しいと思うので忙しい方は画像だけでも!

はじめに

下の写真をみてください。 f:id:aerith7:20171221143139j:plain

Hayabusaの影がItokawaに写っているのがわかるかと思います。

ここから小惑星の大きさを推測する事を考えてみましょう。これだけではHayabusaの影とItokawaの大きさの比率しかわかりません。まずはHayabusaの大きさを知る必要がありそうです。次にどれくらいイトカワから離れているかも知る必要がありそうです。太陽との距離・位置もなければ計算できなさそうです。

他には

  • 何処で写したのか
  • 何時に写したのか
  • どの方向に写したのか
  • どのようなカメラでとったのか
  • そもそも太陽は何処にいるのか
  • カメラは機体のどの位置についていたのか

ざっとこのような情報が必要です。

特に地球では、重力方向にz軸をとって座標にするのが高校物理などでも一般的ですが、 宇宙では何も目印がありません。全てが回転し、全てが動いています(太陽ですらも!)。 さらに機体の内部クロックと実際の時間は違います。世界標準時日本標準時も考慮する必要があります。

このように、観測した実際のデータから何か解析を行うときには、その補助となるデータが必要不可欠となります。 そのようなデータはAncillary Dataと呼ばれています。

Ancillary Dataは生データであり、多岐に渡ります。 このままにしておくと、フォーマット・時間軸・単位・エンコードが異なる大量のデータを管理する羽目になってしまいます。 データが管理しにくいと、ミスの温床になる・労力がかかる・データを利用してもらいにくくなるなどの問題が生じます。

そこで、これらを全てまとめる仕組みとしてSPICEが開発されました。

SPICEとは

このSPICE、NASA内部のNavigation Ancillary Information Facility(NAIF)によって開発されました。

まずSPICEは共有データ・個別データの二種類に大別されます。 共有データはミッションの種類によらないデータ(天体の情報・うるう秒情報・座標系変換など)が格納されています。 また、個別データには各ミッションで必要なデータ(探査機の軌道情報・設置機器の情報)が格納されています。

SPICEの名前は扱うデータの頭文字から取られており、以下の5種類に分類されます。

  • Spacecraft : 人工衛星の軌道情報(NAIFでSPKと呼ばれるファイルはこれ)
  • Planet : 惑星・小惑星・彗星の軌道情報(NAIFでPckと呼ばれるファイルはこれ)
  • Instrument : 機器の視野などの情報(NAIFでIKと呼ばれるファイルはこれ)
  • C-matrix : 科学機器の方向情報(NAIFでCKと呼ばれるファイルはこれ)
  • Events : これは滅多に使用されません(NAIFでEKと呼ばれるファイルはこれ)

そしてこれらの保存されているファイル形式のことをSPICEkernelと呼びます。

このSPICEから具体的な軌道情報を計算するためには航空宇宙に関する豊富な知識が必要とされます。 これを解消するために、NAIFからツールキットが配布されています。

ツールキット

SPICEデータを扱うためのツールキットにはベクトルや行列計算などの関数から、光行差補正に関する関数まで含まれています。

当初ANSI FORTRAN 77のみで実装されていましたが、C、IDL、MATLABJava Native Interfaceでも現在は使用可能となっています。 今回はその中のMATLAB版(MICEと呼ばれます)に絞って解説を進めていきたいと思います。

インストール

インストール手順は以下の通りです。 すでにコンパイルされた状態ですので、ダウンロードしてきて解凍したものにMATLAB側からPATHを通すだけですぐに使えるようになります。 手順は以下の通りです。

  1. サイトからOSに対応するリンクを選択した後、mice.tar.gzをダウンロード、解凍します
  2. 解凍するとmiceフォルダが出てくるので好きな位置におきます(今回はMATLABディレクトリ直下におきます)
  3. MATLABセッションで以下の呪文を唱えます
addpath('/path_to_mice/mice/src/mice/')  
addpath('/path_to_mice/mice/lib/')

4.最後に以下のコマンドを打ってみて動くようなら完了です

cspice_tkvrsn('toolkit')

ツールキットはここに置かれています。 直下のディレクトリに、インストール方法からチュートリアルまで一通り揃っているのでぜひ参照してみてください。

解析の手順

SPICEからの解析は大まかに以下の手順で進みます。

  1. 何を解析するか決める
  2. 対応するkernelファイルをダウンロードしてくる
  3. cspice_furnsh関数を使用してMATLABセッションに呼び出す
  4. 実際に計算する

問題になるのは2のプロセスです。

なんと

どのkernelファイルを落としてくるかは自分で決める

必要があります。関数を叩いて自動で必要なkernelが落とされてくるようにはなっていないのです。 なので自在に使えるようになるためにはある程度の経験が必要になります。 今回の例では何を落とすかは事前に探しておきましたので、それを使いましょう。

Kernelファイルのダウンロード

NAIF所蔵のkernelファイルは全てここに集約されています。 惑星情報・うるう秒情報などのミッションに左右されない共有データはgeneric_kernelsディレクトリを参照しましょう。

ダウンロードは何を使用しても問題ないですが、とりあえずwgetを使います。

wget -r -np --cut-dirs=1 https://naif.jpl.nasa.gov/pub/naif/APOLLO/kernels/fk/apollo_naif_ids.tf

のようにしてディレクトリ構造を維持したままダウンロードしてくる方が後々楽になると思います。 kernelファイルが増えすぎてどれがどれだかわからなくなってしまうからです。

ここからは具体的に軌道を描いてみます。

Cassini Grand Finale

f:id:aerith7:20171222003438j:plain

Cassiniは1977年にNASAESAによって開発された土星探査機で、以下のような実績を残した偉大な探査機です。

  • 飛行距離 : 79億km
  • 発見した土星の衛星 : 6つ
  • 撮影した写真 : 453000枚
  • データ量 : 635GB
  • 論文数 : 3948本

そのCassiniは2017年9月17日に木星の大気圏に突入してその生涯を終えることとなりました。 ここではその最後の足跡を追ってみたいと思います。 まずは必要なkernelファイルを列挙します。

共有データ

個別データ

以上のkernelファイルを元に、まずは

  • 2017年07月27日00:00:00から2017年09月15日00:00:00まで
  • J2000座標系で
    • SPICEでは様々な座標系を扱いますが、全ての座標の基準としてJ2000という座標系を使用しています(解説はこちら)。
    • 特に必要性がなければとりあえずJ2000を選択しましょう。
  • 土星中心からみた
  • Cassiniの軌道

を描画してみようと思います。

% 形状・自転状況
cspice_furnsh('/naif/generic_kernels/pck/pck00010.tpc');
% 惑星の軌道データ
cspice_furnsh('/naif/generic_kernels/spk/planets/de432s.bsp');
% 土星の衛星(Titanなどenceladus)
cspice_furnsh('/naif/generic_kernels/spk/satellites/sat393.bsp');
% 時間
cspice_furnsh('/naif/generic_kernels/lsk/naif0011.tls');
%Cassiniの2017年の203-258日の軌道データ
cspice_furnsh('/naif/CASSINI/kernels/spk/170807BP_SK_17203_17258.bsp');
%Cassiniの2017年の224-258日の軌道データ
cspice_furnsh('/naif/CASSINI/kernels/spk/171003BP_SK_17224_17258.bsp');

このように一通りのkernelファイルを呼び出します。 cspice_furish関数の引数は各kernelファイルへのPATHです。

次に以下のように計算を行います。 軌道を計算するだけならば、以下のコードの変数を変えるだけで全て事足りますので参考にしていただければ嬉しいです。

% 開始時間、終了時間を設定します
et    = cspice_str2et( {'2017-07-27T00:00:00', '2017-09-15T10:30:00'} );
STEP = 100000;

% 100000等分します
times = (0:STEP-1) * ( et(2) - et(1) )/STEP + et(1);

%ここで座標を計算します。
ptarg = mice_spkpos( 'Cassini', times, 'J2000', 'NONE', 'SATURN BARYCENTER' );
pos   = [ptarg.pos];
pos2   = [ptarg2.pos];
pos3   = [ptarg3.pos];
x = pos(1,:);
y = pos(2,:);
z = pos(3,:);

whitebg([25/255 25/255 25/255])
plot3(x,y,z,'Color',[67/255 135/255 233/255],'LineWidth',0.5)
set(gcf, 'InvertHardCopy', 'off')

mice_spkpos( 'Cassini', times, 'J2000', 'NONE', 'SATURN BARYCENTER' )が肝です。

  • times : ある時間内で
  • J2000 : J2000座標系で
  • NONE : 光学補正を行うことなく
  • SATURN BARYCENTER : 土星中心からの
  • Cassini : カッシーニの位置

を計算しています。
計算した結果は以下の図の通りです。

f:id:aerith7:20171222000517j:plain

うーん。。。あってるかどうかよくわかりませんね。。。

そこで以下のオブジェクトを追加していきます。

  • 土星(クリーム色)
    • 半径58,232kmの球を描画しました。
  • 土星の環(クリーム色)
    • cspice_pxform( 'IAU_SATURN','J2000', et ); とすることで、ある時刻における2つの座標系の変換行列を得ることができます。
    • 変換行列から赤道面ベクトルを算出し、そこに最小半径から最大半径に至るまで円を描くことで近似しました。
  • 最後の軌道に移る3ヶ月前の軌道(白)
    • 2016年12月31日00:00:00から2017年04月25日23:00:00までの軌道を同じようにして描画しました。
  • 土星の衛星Titan(オレンジ色)
    • mice_spkposの第一引数を'Cassini'から'Titan'にし、同じように作業することで描画しました。

追加した結果が以下の通りです。

f:id:aerith7:20171221235345j:plain

おぉ。。。それっぽくなってきました。

ところでこの青い軌道、土星ギリギリじゃないですか? この距離からだとくっついてしまっているのかどうかほとんど判別がつきませんし、まるで土星と環の隙間に飛び込んでいるかのようです。 拡大してみましょう。

f:id:aerith7:20171221235304j:plain

おぉ。。。すごい。。確かに環の内側を通っていますね。。 同時に白い軌道も土星の環の外側にかなり接近した軌道であることがわかると思います。

真横からみてみましょう。 f:id:aerith7:20171222003209j:plain この軌道がTitanの軌道のちょうど内側に入るように設計されているのがわかりますでしょうか。 軌道設計の細かさがよくわかると思います。本当にすごいです。。

視点を変えてみます。 f:id:aerith7:20171221235143j:plain Cassiniは22回に渡ってこの軌道で周回しましたが、22回の周回中にTitanの重力を利用しながら少しずつ移動して最後の周回軌道に調整していきました。 青い軌道がずれていることからそれがよくわかるかと思います。

楽しくなってきたので土星の衛星Enceladusも回してみて、
土星の回転座標系も追加してみました。 f:id:aerith7:20171222003228j:plain

Hayabusa

HayabusaのSPICEデータもNAIFサイト内で公開されています。 他の日本の衛星に関してはすでにかぐやが公開し、あかつきも順次公開する予定だそうです。 ここではItokawaからみたHayabusaの着陸軌道(1回目のリハーサル時)の様子をみてみましょう。

使用したkernelファイルは以下の通りです。 共有データ

固有データ

  • hayabusa.tsc
  • HAYABUSA_HP.TF
  • itokawa_1989-2010.bsp
  • hayabusa_itokawarendezvous_v01.bsp

時間は2005年11月03日10:00:00から2005年11月04日10:48:00まで(1回目のリハーサルが行われた時刻です)、

mice_spkpos関数は

mice_spkpos( 'HAYABUSA', times, 'J2000' , 'NONE', 'ITOKAWA' );

のように変更して描画を行いました。 f:id:aerith7:20171222142936j:plain Itokawaの直上700mまで接近したのち、上昇していることがよくわかります。 この時高度700mで自律航法の誤差が規定値を上回ったため、12時30分に地上の指令によってリハーサルを中止したのでした。

追記 : ところでこのItokawa、少しこだわっています。 ここで公開されている多面体モデルをお借りしました。 かなり綺麗に描画できたので(色々こだわってたのに大きさの都合上あまり意味なくて悲しかったので)こちらの方も載せておきます。 f:id:aerith7:20171221235006j:plain

Apollo 15

やっぱりApolloが好きなので挑戦してみます。 aerith7.hatenablog.com

Apollo15号は7月31日に着陸船が分離し、その間司令船は周回し続けています。
そのときに取られた軌道データがSPICEとして管理されています。

APOLLOディレクトリ内部を見ると、 以下の画像のようにApollo15号とApollo16号のデータのみあるのがわかります。

f:id:aerith7:20171222212409p:plain

この中のapollo15-1.cmtを見てみましょう。ここにbspファイルに記録されているデータがいつからいつのものなのかが書かれています。

f:id:aerith7:20171222212945p:plain

1971年7月31日20時から30分間データが途切れていることがわかります。このようなデータもあるのでしっかりコメントは読みましょう(僕はこれを見落としていて、いくらロードしてもうまくいかないので本当に悩んでいました。。)。

使用したkernelファイルは以下の通りです。

共通データ

  • pck00010.tpc
  • de432s.bsp
  • naif0011.tls

固有データ

時間は1971年07月30日01:01:00から1971年07月30日20:00:00、
mices_pkpos関数は

mice_spkpos( 'APOLLO15', times, 'J2000', 'NONE', 'MOON' );

のようにして描画を行いました。結果は以下の通りです。

f:id:aerith7:20171222220216j:plain

周回軌道しっかり回っていますね! せっかくなので月も浮かべてみました(半径1737kmの球に近似しています)。
月の大きさに対する周回軌道の高度の比がよくわかると思います。

その他有用なもの

他国もこのフォーマットに習い、SPICEデータを公開しています。

  • ESA : Rosettaなどなど
  • Russia : Phobos Sample Return Missionを中心に扱っています
  • JAXA : セレーネやHAYABUSAのデータを公開しています。

サードパーティ製ですが、Pythonによるツールキットも公開されています。 github.com

国・年代に跨って幅広いデータを扱っているので命名規則やフォーマットが微妙に異なっている部分があります。必ずreadmeを読み込みましょう。

最後に

  • PC上で土星木星などの惑星や探査機の軌跡まで描けるの楽しい!
  • カッシーニってこんなギリギリの軌道を飛んでたんですね・・・
  • kernelファイルさえ落としてくれば数行で軌道が描けます。
  • 適当なkernelファイルを探してくるのが非常に大変でした。
  • もっと突っ込んだ話が聞きたかった方には少し内容が薄かったかもしれません。すみません。それはまたの機会にということでご勘弁を。。

「宇宙に出た後何をするか」という分野に関しては電気電子の知識が必須になってきますし、将来的にそのような人材が求められているように思います。eeicにおいて宇宙電気電子の分野に進む人はごく少数ですが、 少しでも宇宙に興味を持ってくれる人がいたらそれはとっても嬉しいなって 。

追記 : 最後に1950年から2050年まで100年分の太陽系の軌道を描いてみました!スケール感。。 f:id:aerith7:20171223180634j:plain

参考文献