Go to Qiita Advent Calendar 2024 Top
LoginSignup
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@yamadasuzaku(Shinya Yamada)in立教大学

Pythonで簡単に実装するブラックホール連星Cygnus X-1の軌道位相計算

Posted at

はじめに

連星系の軌道位相の計算は重要です。この記事では、ブラックホール連星 はくちょう座 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 で結果を可視化するサンプルコードです。

コードの説明

  1. 擬似的な FITS ファイルの生成:

    • create_mock_fits 関数で、観測時間のデータを仮定し、それを FITS ファイルに保存します。MJDREF は 55000.0 として仮定しています。
    • 時間データは秒単位で仮定し、TIME カラムに保存します。

  2. FITS ファイルの読み込みと軌道位相の計算:

    • calculate_orbital_phase 関数で、先ほど作成した FITS ファイルを読み込み、観測時間を MJD に変換して、与えられた軌道周期 (To) と基準時刻 (Do) を使用して軌道位相を計算します。

  3. 結果の可視化:

    • 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 ファイルの作成から、位相計算、およびその可視化までの一連の流れを確認できます。

phase_demo_cyg.png

まとめ

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

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

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

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

1
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?