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

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

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

ここでは，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]:
# デフォルト設定
config = {
    'axes.labelsize': 'x-large',
    'mathtext.fontset': 'cm'
}
plt.rcParams.update(config)

$\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, 100)
In [6]:
tl = []
for Omega in Oml:
    tl.append(t(Omega))
    
Tl = []
for Omega in Oml:
    Tl.append(T(Omega))

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

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

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

plt.plot(Oml, Tl, label = r"$\Omega_{\Lambda} = 1 - \Omega_{\rm m}$", 
         color = "red", linewidth = 2)
plt.plot(Oml, tl, label = r"$\Omega_{\Lambda} = 0$", 
         color = "black", linewidth = 2)

# 大目盛に格子線
plt.grid(which="major", linestyle = 'dotted')

# 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(r"$\Omega_{\rm m}$")
plt.ylabel(r"$H_0\ t_0$")

# 凡例の表示
plt.legend(prop={"size":"x-large"});