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

# はじめに

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

https://colab.research.google.com/drive/1bfl14byMAhNK_VTPRew5XR9Jk4OiqODx?usp=sharing

には、[Brocksopp et al. (1999b)](https://ui.adsabs.harvard.edu/abs/1999A%26A...343..861B/abstract) の MJD, HJD, 
[Lasala et al. 1998](https://academic.oup.com/mnras/article/301/1/285/1058242?login=false) の HJD の計算例を載せておきました。この記事では、MJDで計算する最も簡易的な方法を載せています。Cyg X-1 h発見から40年以上経ちますが、現在のところ軌道周期が時間変化したという観測証拠はありません。

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

## Cyg X-1 の軌道情報

参考にする文献は [Brocksopp et al. (1999b)](https://ui.adsabs.harvard.edu/abs/1999A%26A...343..861B/abstract) です。この文献では以下の値が使用されています：

- **軌道周期 (Porb)**: 5.599829 ± 0.000016 日
- **基準時刻 (T0)**: 41874.207 ± 0.009 MJD（O型星が inferior conjunction のとき、つまり、観測者-O星-BHの順番に並ぶ時が phase = 0 に相当する）

https://ui.adsabs.harvard.edu/abs/1999A%26A...343..861B/abstract

:::note info

## JD（Julian Date）をMJDや日付に変換する方法の確認

[Brocksopp et al. (1999b)](https://ui.adsabs.harvard.edu/abs/1999A%26A...343..861B/abstract)

には、

- **基準時刻 (T0)**: 2441874.707 ± 0.009 JD 

として、JD で記載されています。

MJDとJDの変換式は、**MJD=JD−2400000.5** になります。

## JD 2441874.707 を MJD 変換してみる

```python
from astropy.time import Time
# JDの値を指定
jd_value = 2441874.707
t = Time(jd_value, format='jd')
```

## JD → MJD に変換

```python
print("MJD =", t.mjd)
```

出力：
```
MJD = 41874.207
```

> JD 2441874.707 は MJD 41874.207 に対応！

---

## JD → UTC日付（datetime形式）に変換

```python
print("UTC datetime =", t.to_datetime())
```

出力：
```
UTC datetime = 1979-11-05 04:58:04.447999
```

> JD 2441874.707 は **1979年11月5日 04:58:04.45 (UTC)**！

:::


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

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

```python: 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")
```

### 実行例

```bash
%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` ライブラリを使用します。

```python
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 までの範囲に収められ、周期的な変動が確認できるはずです。

``` python:
#!/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, 1000000, 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.mjd + 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 ファイルの作成から、位相計算、およびその可視化までの一連の流れを確認できます。

![スクリーンショット 2025-06-10 14.16.54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/b2868e57-e915-4266-8153-6233a43d4881.png)


### XRISMの観測時期のフェーズの計算

https://github.com/yamadasuzaku/rksysoft/blob/main/astro/cygx1/orbital_phase/astro_cygx1_calc_orbphase.py

を実行すると、

``` bash:
$ ./astro_cygx1_calc_orbphase.py 
Compare orbital phase computed in 3 slightly different methods:
----------------------------------------------------------------------
Brocksopp MJD    | 2024-04-07 17:55:48.99 → Phase = 0.66179
Brocksopp MJD    | 2024-04-10 11:48:20.93 → Phase = 0.15195
Brocksopp HJD    | 2024-04-07 17:55:48.99 → Phase = 0.67867
Brocksopp HJD    | 2024-04-10 11:48:20.93 → Phase = 0.16886
LaSala HJD       | 2024-04-07 17:55:48.99 → Phase = 0.66410
LaSala HJD       | 2024-04-10 11:48:20.93 → Phase = 0.15429
----------------------------------------------------------------------
The orbital phases calculated from these different methods differ only slightly,
typically less than ~0.01 in phase, meaning they are consistent for most purposes.
```

のように３通りの計算結果が一発で出ます

# まとめ

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

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

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

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