はじめに

連星系の軌道位相の計算は重要です。この記事では、ブラックホール連星はくちょう座 X-1 (Cygnus X-1)の計算方法を例にとって、簡単に実行できるPythonのコードも紹介しておきます。

軌道位相を計算する際、時間の単位を揃えることが非常に重要です。例えば、MJD（修正ユリウス日）を使う場合は、すべての計算をMJDで統一しましょう。

Cyg X-1 の軌道情報

参考にする文献は Brocksopp et al. (1999b) です。この文献では以下の値が使用されています：

軌道周期 (Porb): 5.599829 ± 0.000016 日
基準時刻 (T0): 41874.207 ± 0.009 MJD（O型星が inferior conjunction のとき）

Python を使った軌道位相の計算

ここでは、任意の日付からCyg X-1 の軌道位相を計算する方法を紹介します。コードは以下の通りです。

def obs_phase()

from astropy.time import Time

def obs_phase():
    # 軌道周期と基準時刻
    To = 5.599829
    Do = 41874.207

    while True:
        try:
            # ユーザーから観測日時を入力してもらう
            ut = str(input("OBS(UT) (例: 2023-11-04 14:10:21): ")).strip()  # 入力前後の空白を取り除く
            t = Time(ut)  # Time オブジェクトに変換

            # 正しいフォーマットであれば位相を計算して表示
            phase = ((t.mjd - Do) % To) / To
            print(f"軌道位相: {phase}")
            break

        except ValueError:
            # フォーマットが間違っていた場合、再入力を求める
            print("日付フォーマットが正しくありません。以下の形式で再入力してください。")
            print("例: 2023-11-04 14:10:21")

実行例

%python
>>> 上の def obs_phase() を定義

>>> obs_phase()
OBS(UT) (例: 2023-11-04 14:10:21): 2023-11-03 21:45:39.00
軌道位相: 0.8323026272559432
>>> obs_phase()
OBS(UT) (例: 2023-11-04 14:10:21): 2023-11-03 23:49:56.06
軌道位相: 0.8477153416113522

>>> obs_phase()
OBS(UT) (例: 2023-11-04 14:10:21): 2023-11-04 14:10:21.02
軌道位相: 0.9544170125214977
>>> obs_phase()
OBS(UT) (例: 2023-11-04 14:10:21): 2023-11-04 21:21:19.99
軌道位相: 0.007863936394253326

>>> obs_phase()
OBS(UT) (例: 2023-11-04 14:10:21): 2024-04-07 17:55:48.99
軌道位相: 0.6617940997243086
>>> obs_phase()
OBS(UT) (例: 2023-11-04 14:10:21): 2024-04-10 11:48:20.93    
軌道位相: 0.15195444575832176

FITSファイルから軌道位相を計算する

次に、FITSファイルの観測データから直接軌道位相を計算する方法です。astropy ライブラリを使用します。

from astropy.io import fits
from astropy.time import Time

# FITSファイルを開く
hdu = fits.open("fitsファイル名")

# 観測時間列を取得 (MJDREF からの秒数)
time = hdu[1].data["TIME"]

# MJDREF (基準時刻) を取得
mjdref = hdu[0].header["MJDREF"]

# 基準時刻を astropy の Time オブジェクトに変換
t_mjdref = Time(mjdref, format='mjd')

# 観測時間を MJD に変換
t_mjd = Time(t_mjdref.jd + time / 86400., format='mjd')

# 軌道位相の計算
To = 5.599829  # 軌道周期
Do = 41874.207 # 基準時刻 (MJD)

phase = []
for onetime in t_mjd.mjd:
    # 位相を計算してリストに追加
    onephase = (onetime - Do) % To / To # オフセットを引いて、軌道周期で余りを求めて、軌道周期で割るので、0-1 になる。

    phase.append(onephase)

位相の図の描き方

天文学では、位相の図を描く際に、横軸を2周期分表示することがよくあります。これは、一周分を繰り返してプロットするためです。図を描く際には、この点を意識すると見やすくなります。

擬似的な FITS ファイルの生成と Phase 計算の例

以下は、擬似的な FITS ファイルを astropy を使って生成し、その後に与えられたコードを使って軌道位相を計算し、matplotlib で結果を可視化するサンプルコードです。

コードの説明

擬似的な FITS ファイルの生成:
- create_mock_fits 関数で、観測時間のデータを仮定し、それを FITS ファイルに保存します。MJDREF は 55000.0 として仮定しています。
- 時間データは秒単位で仮定し、TIME カラムに保存します。
FITS ファイルの読み込みと軌道位相の計算:
- calculate_orbital_phase 関数で、先ほど作成した FITS ファイルを読み込み、観測時間を MJD に変換して、与えられた軌道周期 (To) と基準時刻 (Do) を使用して軌道位相を計算します。
結果の可視化:
- matplotlib を使って、観測時刻（MJD）に対する軌道位相をプロットします。位相は 0 から 1 までの範囲に収められ、周期的な変動が確認できるはずです。

#!/usr/bin/env python

import numpy as np
from astropy.io import fits
from astropy.time import Time
import matplotlib.pyplot as plt

# 擬似的な FITS ファイルを生成する関数
def create_mock_fits(filename="mock_fits.fits"):
    # 観測時間データ（MJDREFからの秒数で仮定）
    nrows = 100  # 行数（観測点の数）
    time_data = np.linspace(0, 500000, nrows)  # 仮の時間データ（秒）

    # MJDREF (基準時刻) を仮定
    mjdref = 55000.0  # 仮の基準MJD

    # FITSヘッダーとデータを作成
    col1 = fits.Column(name='TIME', format='E', array=time_data)
    cols = fits.ColDefs([col1])
    hdu1 = fits.BinTableHDU.from_columns(cols)

    # プライマリヘッダーにMJDREFを追加
    hdu0 = fits.PrimaryHDU()
    hdu0.header['MJDREF'] = mjdref

    # FITSファイルを作成
    hdul = fits.HDUList([hdu0, hdu1])
    hdul.writeto(filename, overwrite=True)

# FITSファイルを生成
create_mock_fits()

# FITSファイルを読み込んで位相を計算する関数
def calculate_orbital_phase(fits_filename):
    # FITSファイルを開く
    hdu = fits.open(fits_filename)

    # 観測時間列を取得 (MJDREF からの秒数)
    time = hdu[1].data["TIME"]

    # MJDREF (基準時刻) を取得
    mjdref = hdu[0].header["MJDREF"]

    # 基準時刻を astropy の Time オブジェクトに変換
    t_mjdref = Time(mjdref, format='mjd')

    # 観測時間を MJD に変換
    t_mjd = Time(t_mjdref.jd + time / 86400., format='mjd')

    # 軌道位相の計算
    To = 5.599829  # 軌道周期 (日)
    Do = 41874.207 # 基準時刻 (MJD)

    phase = []
    for onetime in t_mjd.mjd:
        # 位相を計算してリストに追加
        onephase = (onetime - Do) % To / To  # 0-1 の範囲に収める
        phase.append(onephase)

    # 計算結果をプロット
    plt.plot(t_mjd.mjd, phase, 'o-', label="orbital phase")
    plt.xlabel("MJD")
    plt.ylabel("Phase (0-1)")
    plt.title("Orbital Phase of Cyg X-1")
    plt.legend()
    plt.savefig("phase_demo_cyg.png")    
    plt.show()

# 軌道位相の計算とプロット
calculate_orbital_phase("mock_fits.fits")

これにより、擬似データを使って FITS ファイルの作成から、位相計算、およびその可視化までの一連の流れを確認できます。

まとめ

この記事では、連星系の軌道位相の計算方法について、ブラックホール連星 Cygnus X-1 を例に、詳細なPythonコードを使用して解説しました。具体的には、以下のポイントに注目しました。

軌道位相の基本的な計算方法: astropyライブラリを使用して、任意の日付からCygnus X-1の軌道位相を計算する方法を紹介しました。
FITSファイルからの計算: FITSファイルに保存された観測データから軌道位相を計算する方法も解説し、天文学的データ解析の実践的な例を示しました。
擬似的なFITSファイルの生成と可視化: 自分でFITSファイルを生成し、そのデータをもとに位相を計算し、matplotlibを使ってグラフ化する方法も提供しました。

これらのコードを使えば、Cygnus X-1 の軌道位相の計算方法を理解し、実際のデータ解析に応用できるでしょう。

中性子星の自転位相や、軌道位相が変化する連星系の場合などは、もう少し複雑になりますし、バリセントリック補正(e.g., barycorr )など、時間の扱いに注意を払う必要が出てきます。