光度とエディントン光度比を計算し可視化する方法(備忘録)

はじめに

天文学において、天体の光度とエディントン光度比は、天体の物理的性質を理解するために重要な指標です。本記事では、これらの量を計算し、結果を可視化するPythonコードを紹介します。このコードでは、与えられたフラックス、質量、距離に基づいて光度とエディントン光度を計算し、その比をグラフで表示します。

使い方

• フラックス(erg/cm^2/s)
• mass(太陽質量単位)
• 距離(kpc)

をプログラムに与えると、計算結果の表示と簡単な図が生成されます。

• python astro_util_calc_luminosity.py フラックス(erg/cm^2/s) mass(太陽質量単位) 距離(kpc) --outfile 出力ファイル名

astro_util_calc_luminosity.py は、

に置いています。

Cyg X-1 の例

$python astro_util_calc_luminosity.py 4.8538e-8 21.2 2.2 --outfile newcyg.png
Calculated luminosity: 2.811e+37 erg/s
Luminosity to Eddington luminosity ratio: 1.054e-02
Eddington luminosity for input mass: 2.668e+39 erg/s

赤い点が入力したフラックス、質量、距離に相当する光度とエディントン比になります。

1. 光度とエディントン光度とは？

光度 (Luminosity)

光度は天体が単位時間あたりに放射するエネルギーの総量です。光度は天文学では cgs単位系で、erg/s（エルグ毎秒）で表され、天体の明るさやエネルギー放出の強さを示します。

光度は次の式で求められます：

L = 4 \pi d^2 F

ここで：

• $L$ は光度（erg/s）、
• $d$ は天体からの距離（cm）、
• $F$ は天体から放射されるフラックス（erg/cm²/s）です。

エディントン光度 (Eddington Luminosity)

エディントン光度は、天体がその質量に対して放射できる最大の光度です。天体の質量が大きいほど、エディントン光度も大きくなります。エディントン光度は次の式で計算されます：

L_{\text{Edd}} = \frac{4 \pi G M m_p c}{\sigma_T}

ここで：

• $L_{\text{Edd}}$ はエディントン光度（erg/s）、
• $G$ は万有引力定数、
• $M$ は天体の質量（g）、
• $m_p$ は陽子質量（g）、
• $c$ は光速（cm/s）、
• $\sigma_T$ はトムソン散乱断面積（cm²）です。

光度とエディントン光度比

光度とエディントン光度の比は、天体の放射がエディントン制限をどれだけ超えているかを示します。この比が1より大きい場合、天体はエディントン光度を超えていることになります。

\frac{L}{L_{\text{Edd}}}

2. 必要なライブラリのインポート

以下のPythonコードでは、必要なライブラリをインポートし、ユーザーから入力を受け取って計算を行います。

import numpy as np
import matplotlib.pyplot as plt
import math
import argparse

• numpy は数値計算用のライブラリで、特に対数スケールやデータ処理に使います。
• matplotlib.pyplot はデータを可視化するためのライブラリです。
• math は数学的な定数や関数（例えば、円周率π）を使用するためにインポートします。
• argparse はコマンドライン引数を処理するために使用します。

3. 光度とエディントン光度比を計算する関数

次に、光度とエディントン光度を計算する関数を定義します。フラックス、質量、距離を入力として、光度とエディントン光度比を計算します。

def calculate_luminosity_ratio(flux, distance_kpc, mass_solar):
    # 距離の単位変換
    distance_cm = distance_kpc * 3.086e21  # kpc → cm

    # 光度の計算
    luminosity = 4 * math.pi * distance_cm**2 * flux

    # エディントン光度の計算
    mass = mass_solar * 1.989e33  # 太陽質量をグラム単位に変換
    L_edd = (4 * math.pi * 6.67430e-8 * mass * 1.6726e-24 * 3e10) / 6.652e-25  # エディントン光度

    # 光度とエディントン光度の比を計算
    luminosity_ratio = luminosity / L_edd
    return luminosity, luminosity_ratio, L_edd

説明：

• distance_kpc は距離（kpc単位）で、これをcmに変換します。
• 光度 $ L $ はフラックス $ F $ と距離 $ d $ を使って計算されます。
• エディントン光度 $ L_{\text{Edd}} $ は質量と定数を使って計算されます。

4. メイン関数：引数の受け取りと計算、プロット

ユーザーからフラックス、質量、距離を受け取り、光度とエディントン光度比を計算して可視化します。

コード全文

astro_util_calc_luminosity.py

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'serif'
import math
import argparse

# 定数の設定
G = 6.67430e-8  # 万有引力定数 (cm^3/g/s^2)
m_p = 1.6726e-24  # 陽子質量 (g)
c = 3e10  # 光速 (cm/s)
sigma_T = 6.652e-25  # トムソン散乱断面積 (cm^2)
solar_mass = 1.989e33  # 太陽質量 (g)

# 関数：光度とエディントン光度の比を計算
def calculate_luminosity_ratio(flux, distance_kpc, mass_solar):
    # 距離の単位変換
    distance_cm = distance_kpc * 3.086e21  # kpc → cm

    # 光度の計算
    luminosity = 4 * math.pi * distance_cm**2 * flux

    # エディントン光度の計算
    mass = mass_solar * solar_mass  # 質量 (g)
    L_edd = (4 * math.pi * G * mass * m_p * c) / sigma_T  # エディントン光度

    # 光度とエディントン光度の比を計算
    luminosity_ratio = luminosity / L_edd
    return luminosity, luminosity_ratio, L_edd

# argparse を使用して引数を処理
def main():
    parser = argparse.ArgumentParser(description='Calculate luminosity and Eddington luminosity ratio.')
    parser.add_argument('flux', type=float, help='Flux in erg/cm^2/s')
    parser.add_argument('mass', type=float, help='Mass in solar masses')
    parser.add_argument('distance', type=float, help='Distance in kpc')
    parser.add_argument('--flux_min', type=float, default=-10, help='Minimum flux range (default: -10)')
    parser.add_argument('--flux_max', type=float, default=-6, help='Maximum flux range (default: -6)')
    parser.add_argument('--outfile', type=str, default="mass_lumi_limiedd.png", help='output file name')

    
    args = parser.parse_args()

    # 引数から値を取得
    flux = args.flux
    mass_solar = args.mass
    distance_kpc = args.distance
    flux_min = args.flux_min
    flux_max = args.flux_max
    outfile = args.outfile

    # 光度とエディントン光度の比を計算
    user_luminosity, user_luminosity_ratio, user_L_edd = calculate_luminosity_ratio(flux, distance_kpc, mass_solar)

    # 結果を表示
    print(f"Calculated luminosity: {user_luminosity:.3e} erg/s")
    print(f"Luminosity to Eddington luminosity ratio: {user_luminosity_ratio:.3e}")
    print(f"Eddington luminosity for input mass: {user_L_edd:.3e} erg/s")

    # 可視化
    flux_values = np.logspace(flux_min, flux_max, 100)  # フラックスの範囲を変更可能に
    luminosities = []
    luminosity_ratios = []

    for flux_val in flux_values:
        luminosity, luminosity_ratio, _ = calculate_luminosity_ratio(flux_val, distance_kpc, mass_solar)
        luminosities.append(luminosity)
        luminosity_ratios.append(luminosity_ratio)

    # 可視化
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))

    # 光度のプロット
    ax1.plot(flux_values, luminosities, label='Luminosity (erg/s)', color='blue')
    ax1.plot(flux, user_luminosity, "o", label='(USER) Luminosity (erg/s)', color='red')
    ax1.set_xscale('log')
    ax1.set_yscale('log')
    ax1.set_xlabel(r'Flux (erg/cm$^2$/s)')
    ax1.set_ylabel('Luminosity (erg/s)')
    ax1.set_title('Luminosity vs Flux')
    ax1.legend()

    # エディントン光度比のプロット
    ax2.plot(flux_values, luminosity_ratios, label=r'$L / L_{\rm{Edd}}$', color='blue')
    ax2.plot(flux, user_luminosity_ratio, "o", label=r'(USER) $L / L_{\rm{Edd}}$', color='red')

    ax2.set_xscale('log')
    ax2.set_yscale('log')
    ax2.set_xlabel(r'Flux (erg/cm$^2$/s)')
    ax2.set_ylabel(r'$L / L_{\rm{Ledd}}$')
    ax2.set_title(r'Luminosity to Eddington Luminosity Ratio')
    ax2.legend()

    # 上部余白を調整し、テキストを配置
    fig.subplots_adjust(top=0.8)  # 上部の余白を大きく設定

    # 入力パラメータのテキスト表示
    input_text_str = (f"Input Parameters:\n"
                      f"Flux: {flux} erg/cm²/s\n"
                      f"Mass: {mass_solar}" r"$M_{\odot}$\n"
                      f"Distance: {distance_kpc} kpc")

    # 計算結果のテキスト表示
    output_text_str = (f"Calculated Results:\n"
                       f"Luminosity: {user_luminosity:.3e} erg/s\n"
                       f"Eddington Luminosity: {user_L_edd:.3e} erg/s\n"
                       f"Luminosity / Eddington Luminosity Ratio: {user_luminosity_ratio:.3e}")

    # テキストを図の上部に表示
    plt.figtext(0.05, 0.87, input_text_str, transform=fig.transFigure, fontsize=8,
                verticalalignment='bottom', horizontalalignment='left', color='black',
                bbox=dict(facecolor='white', alpha=0.8))

    plt.figtext(0.55, 0.87, output_text_str, transform=fig.transFigure, fontsize=8,
                verticalalignment='bottom', horizontalalignment='left', color='black',
                bbox=dict(facecolor='white', alpha=0.8))

#    plt.tight_layout()
    plt.savefig(outfile)
    plt.show()

if __name__ == '__main__':
    main()

まとめ

このコードは、与えられたフラックス、質量、距離に基づいて光度とエディントン光度を計算し、その結果を対数スケールで可視化します。また、テキストを図の上部に表示して、計算結果を簡単に確認できるようにしています。