宇宙の数式を頑張って読み解いてみる(Interface 2023 7月号)

概要

※この記事は上記に基づきます

先日 2023 年 7 月号の Interface を買ったところ、付録の別冊プログラミング学園 Python 部に面白い記事が有った。
JAXA 橋本樹明先生の記事で題は「宇宙の数式」、NASA 等から提供されている各種観測データを用いて衛星軌道をトレースしてみるという内容。
知らないことばかりだったので頑張って読み解いてみる。

SPICE って何?

その昔は宇宙開発競争もあり、それぞれの宇宙開発プロジェクト毎にデータフォーマットが異なっていたらしい。
しかしながらオープンソースの潮流や各国の協調した開発プロジェクトが始まるにつれ共通したデータフォーマットの需要が生まれ、NASA の Navigation Ancillary Information Facility(NAIF) が開発した SPICE(Spacecraft Planet Instrument C - matrix Events) が事実上の標準フォーマットとなった。

• SPICE に格納されているデータ

データ名	内容
SPK	探査機や自然天体の軌道情報
PCK	天体の自転軸や形状の情報
IK	観測機器の視野、形状の情報
CK	探査機の姿勢、観測機器の取り付け方向
EK	イベント情報
FK	アンテナ取り付け角など基準となる情報
SCLK	探査機の時刻カウンタと世界時の変換情報
LSK	閏秒の情報
DSK	観測する天体の3次元情報

衛星軌道を追ってみる

SPICE を使用するツールキットは C, IDL, MATLAB, Java 等で提供されており、また Python でもサードパーティ製の spiceypy が提供されている。
今回は spiceypy を用いてイトカワへランデブーしたハヤブサ1の軌道を追う。

ハヤブサ1とは

はやぶさ（第20号科学衛星MUSES-C）は、2003年5月9日13時29分25秒（日本標準時、以下同様）に宇宙科学研究所（ISAS）が打ち上げた小惑星探査機で、ひてん、はるかに続くMUSESシリーズ3番目の工学実験機である。開発・製造はNEC東芝スペースシステムが担当した。

イオンエンジンの実証試験を行いながら2005年夏にアポロ群の小惑星（25143）イトカワへ到達し、その表面を詳しく観測してサンプル採集を試みた後、2010年6月13日22時51分、60億 kmの旅を終え地球に帰還し、大気圏に再突入した。地球重力圏外にある天体固体表面に着陸してのサンプルリターンに、世界で初めて成功した。
出典: wikipedia

コード解説

import numpy as np
import spiceypy as spice
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
#　軌道データ読み込み
#　ハヤブサ1
spice.furnsh('/content/hay_osbj_050911_051118_v1.bsp')
# 各惑星
spice.furnsh('/content/de403s.bsp')
# イトカワ
spice.furnsh('/content/sb_25143_140.bsp')

#うるう秒の取得
spice.furnsh('/content/naif0009.tls')

# 探査機IDの設定
scid = spice.bodn2c('HAYABUSA')
scname =spice.bodc2n(scid,256)
print(scname)

# ハヤブサ1の軌道データが格納されている期間を取得
coverage_window = spice.spkcov(spk='/content/hay_osbj_050911_051118_v1.bsp', idcode=scid)
print("格納期間:", spice.et2utc(coverage_window[0], "ISOC", 6), "to", spice.et2utc(coverage_window[1], "ISOC", 6))

et_start = coverage_window[0]
et_end = coverage_window[1]
et_step = 3600*24;
times = np.arange(et_start, et_end, et_step)

# 各軌道を取得
hayabusa_positions, _ = spice.spkpos('HAYABUSA', times, 'J2000', 'NONE', 'ITOKAWA')
earth_positions, _ = spice.spkpos('EARTH', times, 'J2000', 'NONE', 'SUN')
venus_positions, _ = spice.spkpos('VENUS', times, 'J2000', 'NONE', 'SUN')
mercury_positions, _ = spice.spkpos('MERCURY', times, 'J2000', 'NONE', 'SUN')
itokawa_positions, _ = spice.spkpos('ITOKAWA', times, 'J2000', 'NONE', 'SUN')

# ハヤブサ1の軌道はイトカワ基準の軌道しかデータが存在しなかったため、太陽基準のイトカワの軌道データと合成
hayabusa_to_sun = hayabusa_positions + itokawa_positions

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot(hayabusa_to_sun[:, 0], hayabusa_to_sun[:, 1], hayabusa_to_sun[:, 2], color='black',label='HYABUSA')
ax.plot(earth_positions[:, 0], earth_positions[:, 1], earth_positions[:, 2], color='blue',label='EARTH')
ax.plot(venus_positions[:, 0], venus_positions[:, 1], venus_positions[:, 2], color='red',label='VENUS')
ax.plot(mercury_positions[:, 0], mercury_positions[:, 1], mercury_positions[:, 2], color='cyan',label='MERCURY')
ax.plot(itokawa_positions[:, 0], itokawa_positions[:, 1], itokawa_positions[:, 2], color='green',label='ITOKAWA')


ax.set_xlabel('X [km]')
ax.set_ylabel('Y [km]')
ax.set_zlabel('Z [km]')

plt.legend()
plt.show()

• furnsh()

  • SPICE に格納された軌道データのカーネルファイルを読み込み、使用可能にする関数
    • ハヤブサ1, 各種天体, イトカワの軌道情報、閏秒等の時刻情報をロードしているカーネル

• bodn2c()

  • 探査機名をパラメータとして、探査機の ID に変換する関数

• bodc2n()

  • 探査機 ID をパラメータとして、探査機名に変換する関数

• spkcov()

  • 指定した軌道データに対して、どの期間の軌道データがカバーされているかを取得する関数
    • coverage_window には開始時間と終了時間が ET (ユリウス年)に基づいて格納される

• times

  • 起動計算がカバーする範囲を一日毎に区切ってリストを作成

• spkpos()

  • 指定したターゲットと基準点に基づいて相対的な位置ベクトルを返却する関数
    • イトカワに基づくハヤブサ
    • 太陽に基づく地球
    • 太陽に基づく金星
    • 太陽に基づく彗星
    • 太陽に基づくイトカワ

• hayabusa_to_sun

  • ハヤブサからイトカワへ、イトカワから太陽へのベクトルを合成することで、擬似的に太陽に対するハヤブサの軌道を取得している

描画結果

ハヤブサが描画されていない。

妙だなと思い、イトカワの描画を削除してみると...

イトカワと重なった軌道に！
それもその筈

2005年9月には小惑星「イトカワ」とランデブーした。約5か月の小惑星付近滞在中、カメラやレーダーなどによる科学観測を行った
出典: wikipedia

この期間はイトカワとほぼ同軌道におり、またイトカワとランデブーを果たしていたので、軌道計算上重なった位置を航行していたのである。

終わりに

KSP 等で楽々軌道計算をしていた身としては、自身の手で軌道計算が出来るのは面白い体験だった。
今なお航行するボイジャー1号が無事にヘリオポーズを超えて三体の世界が来たら面白いなあと思う。