Chandra軌道データを3Dアニメーションで可視化してみた

はじめに

Chandra X線天文衛星の太陽系重心時刻（Barycentric Julian Date, BJD）などを計算する際に欠かせないのが軌道情報です。Chandra衛星の軌道情報はeph1ファイルに軌道データが格納されており、その中身には時刻ごとの衛星の位置と速度が記録されています。

この記事では、その軌道情報がどのようなものかを直感的に理解するためにアニメーションで可視化してみました。Pythonとmatplotlib、astropyを使って、Chandraの一周分の軌道を3DアニメーションGIFとして可視化出力しています。

Chandraの軌道とは

Chandra X線天文衛星は、NASAが1999年に打ち上げた高空間分解能（0.5秒角）の撮像能力を持つ、現在も運用中の人工衛星です。

特徴的なのはその「高楕円軌道（highly elliptical orbit）」です。近地点が約1万6千km、遠地点が最大で約13万km以上にも達し、軌道周期は約64時間。この軌道により、地球の影響を最小限に抑えつつ、長時間にわたる観測が可能になります。

今回のアニメーションでも、速度が近地点で速く、遠地点でゆっくりになる様子が確認できます。

Chandraの周回軌道実装コード

カラム情報の確認

まずは衛星ごとに必ずカラムを確認するようにします。これは単位が衛星ごとで違う可能性があるからです。

カラムの確認

from astropy.io import fits

filename = "orbitf077976300N001_eph1.fits"
with fits.open(filename) as hdul:
    print(repr(hdul[1].columns))

出力結果

ColDefs(
    name = 'Time'; format = 'D'; unit = 's'
    name = 'X'; format = 'D'; unit = 'm'
    name = 'Y'; format = 'D'; unit = 'm'
    name = 'Z'; format = 'D'; unit = 'm'
    name = 'Vx'; format = 'D'; unit = 'm/s'
    name = 'Vy'; format = 'D'; unit = 'm/s'
    name = 'Vz'; format = 'D'; unit = 'm/s'
)

• Time: 時刻情報（MET (Mission Elapsed Time)で、Chandraの場合は、1998-01-01 00:00:00 (TT)）
• x, y, z: 衛星の位置
• Vx, Vy, Vz: 衛星の速度

このように時刻と3次元の位置と速度が決まれば運動を記述することができます。これらの情報を使って実際にアニメーションを作成してそのスケールを体感します。

3Dアニメーション表示

Chandraの周回軌道

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
from astropy.io import fits

# FITSファイル読み込み
# 例: orbitf077976300N001_eph1.fits
filename = "orbitf077976300N001_eph1.fits"
with fits.open(filename) as hdul:
    data = hdul[1].data

x = data['X']
y = data['Y']
z = data['Z']
vx = data['VX']
vy = data['VY']
vz = data['VZ']
t = data['TIME'] - data['TIME'][0]  # Elapsed Time（秒）

# 全データの1/7（約1周）だけを使用
portion = 1 / 7
max_index = int(len(t) * portion)

# そこから5点に1点だけを抽出（軽量化のため）
step = 5
x = x[:max_index:step]
y = y[:max_index:step]
z = z[:max_index:step]
vx = vx[:max_index:step]
vy = vy[:max_index:step]
vz = vz[:max_index:step]
t = t[:max_index:step]

# 速度の大きさを計算
speed = np.sqrt(vx**2 + vy**2 + vz**2)
n_points = len(x)

# 描画範囲の自動設定
max_range = max(np.ptp(x), np.ptp(y), np.ptp(z)) / 2
x_mid = (np.max(x) + np.min(x)) / 2
y_mid = (np.max(y) + np.min(y)) / 2
z_mid = (np.max(z) + np.min(z)) / 2

# プロット準備
fig = plt.figure(figsize=(5, 5), dpi=72)
ax = fig.add_subplot(111, projection='3d')
fig.subplots_adjust(left=0.08, right=0.88, bottom=0.08, top=0.92)

# 地球の球体を描画
earth_radius = 6371*1e3
u = np.linspace(0, 2 * np.pi, 50)
v = np.linspace(0, np.pi, 50)
xe = earth_radius * np.outer(np.cos(u), np.sin(v))
ye = earth_radius * np.outer(np.sin(u), np.sin(v))
ze = earth_radius * np.outer(np.ones(np.size(u)), np.cos(v))
earth_point = ax.plot_surface(xe, ye, ze, color='green', alpha=0.5)


orbit_line, = ax.plot([], [], [], lw=2, color='blue')
satellite_point, = ax.plot([], [], [], 'ro')

ax.set_xlim([x_mid - max_range, x_mid + max_range])
ax.set_ylim([y_mid - max_range, y_mid + max_range])
ax.set_zlim([z_mid - max_range, z_mid + max_range])
ax.set_box_aspect([1, 1, 1])
ax.set_xlabel('X [m]')
ax.set_ylabel('Y [m]')
ax.set_zlabel('Z [m]')
ax.set_title('Chandra Orbit Animation')

ax.view_init(elev=30, azim=60)

# 時刻＋速度のテキスト表示
info_text = ax.text2D(0.05, 0.95, "", transform=ax.transAxes)

# 初期化関数
def init():
    orbit_line.set_data([], [])
    orbit_line.set_3d_properties([])
    satellite_point.set_data([], [])
    satellite_point.set_3d_properties([])
    info_text.set_text("")
    return orbit_line, satellite_point, earth_point, info_text

# フレーム更新関数
def update(frame):
    orbit_line.set_data(x[:frame], y[:frame])
    orbit_line.set_3d_properties(z[:frame])
    satellite_point.set_data([x[frame]], [y[frame]])
    satellite_point.set_3d_properties([z[frame]])
    info_text.set_text(
        f"Elapsed Time: {t[frame]*1e-3:.0f} ks\nVelocity: {speed[frame]:.2f} m/s"
    )
    return orbit_line, satellite_point, earth_point, info_text

# アニメーション作成・保存
ani = FuncAnimation(fig, update, frames=n_points, init_func=init,
                    interval=30, blit=False)

ani.save("chandra_orbit_with_earth.gif", writer="pillow", fps=30)
print("GIF saved as chandra_orbit_with_earth.gif")

plt.show()

出力結果

赤がChandara衛星、青がその軌道、緑が地球です。衛星以外は実寸スケールで描画しています。

説明

この例では、元の軌道データには約7周期分の情報が含まれていました。そこで、その中から約1周期分だけを抽出し、さらに5点ごとに間引くことで、軽量かつスムーズなアニメーションを実現しました。

補足

衛星が軌道上を動きながら観測するため、観測された光子が衛星に届くタイミングにわずかなずれが生じます。このようなずれは観測結果に影響を与えるため、データ解析では注意が必要です。

これを補正するために用いられるのが、太陽系の重心を基準としたBJD時刻です。Chandra衛星のデータにおけるBJDへの時刻補正は、公式解析ツールであるCIAOで行うことができますので、詳しくはCIAOのドキュメントをご参照ください。

• CIAO axbary

• Chandraの軌道情報一覧（現在利用可能なすべてのeph1ファイル）

まとめ

Chandraの軌道情報は、BJD時刻の正確な補正に欠かせないものですが、数値だけを見ていてもなかなかイメージができないため今回可視化してみました。この記事がBJD時刻を求める方のモチベーションになれば幸いです。