Return to 宇宙論パラメータと宇宙年齢

補足：Python で宇宙年齢のグラフを描く

宇宙年齢の表式の導出については，以下のページを参照。

ここでは，Python の matplotlib.pyplot.plot() を使って宇宙年齢のグラフを描いてみる。

In [1]:
# NumPy も使います
import numpy as np
# Matplotlib でグラフを描きます
import matplotlib.pyplot as plt

# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
In [2]:
# savefig した svg が 600 x 450 になるように 72 で割った値にしてみる
plt.rcParams["figure.figsize"] = (8.3334, 6.25)
# font size の設定
plt.rcParams["font.size"] = 12

$\Omega_{\Lambda} = 0$ の場合

$$H_0 t_0 = -\frac{1}{\Omega_{\rm m} -1}+\frac{\Omega_{\rm m}}{(\Omega_{\rm m}-1)^{\frac{3}{2}} }
\tan^{-1}\sqrt{\Omega_{\rm m}-1} \quad \mbox{for}\ \ \Omega_{\rm m} > 1$$$$H_0 t_0 = \frac{1}{1-\Omega_{\rm m}}-\frac{\Omega_{\rm m}}{(1-\Omega_{\rm m})^{\frac{3}{2}} }
\tanh^{-1}\sqrt{1-\Omega_{\rm m}} \quad \mbox{for}\ \ \Omega_{\rm m} < 1$$

$\Omega_{\rm m} \rightarrow \Omega$ として…

独り言：あぁそう。NumPy は $\tanh^{-1} x$ を np.arctanh(x) としているんだぁ。逆双曲線関数は arc と書くべきではないんだよなぁ。

In [3]:
# Omega > 1
def t1(Omega):
    return -1/(Omega-1) + Omega/((Omega-1)*np.sqrt(Omega-1)) * np.arctan(np.sqrt(Omega-1))
# Omega < 1
def t2(Omega):
    return 1/(1-Omega) - Omega/((1-Omega)*np.sqrt(1-Omega)) * np.arctanh(np.sqrt(1-Omega))
# Omega = 1
def t0(Omega):
    return 2/3

def t(Omega):
    if Omega > 1:
        return t1(Omega)
    elif Omega == 1:
        return t0(Omega)
    else:
        return t2(Omega)

$k = 0$ の場合

$$H_0 t_0 = \frac{2}{3\sqrt{\Omega_{\rm m} -1}}\tan^{-1} \sqrt{\Omega_{\rm m} -1} \quad \mbox{for}\ \ \Omega_{\rm m} > 1$$$$H_0 t_0 = \frac{2}{3\sqrt{1-\Omega_{\rm m} }}\tanh^{-1} \sqrt{1-\Omega_{\rm m} } \quad \mbox{for}\ \ \Omega_{\rm m} < 1$$

$\Omega_{\rm m} \rightarrow \Omega$ として…

In [4]:
# Omega > 1
def T1(Omega):
    return 2/(3*np.sqrt(Omega-1)) * np.arctan(np.sqrt(Omega-1))
# Omega < 1
def T2(Omega):
    return 2/(3*np.sqrt(1-Omega)) * np.arctanh(np.sqrt(1-Omega))

def T(Omega):
    if Omega > 1:
        return T1(Omega)
    elif Omega == 1:
        return t0(Omega)
    else:
        return T2(Omega)

横軸・縦軸の数値リストを用意する

In [5]:
Oml = np.linspace(0.001, 2, 200)
In [6]:
tl = []
for Omega in Oml:
    tmp = t(Omega)
    tl.append(tmp)
    
Tl = []
for Omega in Oml:
    TMP = T(Omega)
    Tl.append(TMP)

オプションを設定して plot する

In [7]:
# tick の向きを内側に
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'

#グラフの上下左右に目盛線を付ける。
f = plt.figure()
ax = f.add_subplot(111)
ax.yaxis.set_ticks_position('both')
ax.xaxis.set_ticks_position('both')

# 大目盛に格子線
plt.grid(which="major", color="lightgray", ls="--", linewidth=0.5)

# x の大目盛を 0.2 刻みに
plt.xticks(np.arange(0.2, 2.1, 0.2))
# y の大目盛を 0.1 刻みに
plt.yticks(np.arange(0.5, 1.6, 0.1))
# 小目盛も表示
plt.minorticks_on()

# 横軸の表示範囲
plt.xlim([0, 2])
plt.ylim([0.5, 1.5])
# タイトル
plt.title("宇宙年齢の密度パラメータ依存性")
# 軸のラベル
plt.xlabel("$Ω_m$")
plt.ylabel("$H_0\ t_0$")

key = "$Ω_Λ = 1 - Ω_m$"
plt.plot(Oml, Tl, label = key, color = "red", linewidth = 2)

key = "$Ω_Λ = 0$"
plt.plot(Oml, tl, label = key, color = "black", linewidth = 2)
# 凡例の表示
plt.legend();