宇宙年齢を求める

astro

（知らない人には不親切な）理屈

宇宙論の教科書(Cosmological Physics, J. A. Peacock とか)に
H^2(a)=H_0^2\left[\Omega_{\Lambda}+\Omega_{m} a^{-3} +\Omega_{r}a^{-4}+(\Omega -1)a^{-2}\right],\quad(\Omega=\sum_i\Omega_i)
って書いてあるのでH(a)=\dot{a}/aであることを使えば，
dt=\frac{1}{H_0}\left[\Omega_{\Lambda}+\Omega_{m} a^{-3} +\Omega_{r}a^{-4}+(\Omega -1)a^{-2}\right]^{-1/2}da/a
とかって書き換えられますから．これをt:0\rightarrow t_0つまり，a:0\rightarrow 1まで積分するとでる．
hubble parameter h を使うと 1/H_0=9.78\times10^{9}/h {\rm year}
はやりの宇宙論パラメータは

item value
h 0.71
\Omega_{\Lambda} 0.73
\Omega_{m} 0.27
\Omega_{r} 0

なので宇宙の年齢は137億年と求まります．

実際に計算する

準備

まず scipy を使うので scipy をいれましょう．たぶん numpy もいるはずです．
せっかくなので scipy を使って積分してみたいと思います．まずは scipy.integrate.__doc__ を調べてみるといくつか積分法があります．どれを使うか考えるのも面倒なのでここでは quad (General purpose integration.) を使うことにします．でもやっぱり使い方が分からないので scipy.integrate.quad.__doc__ をプリントしてみます．どうやら quad(func,a,b) で
\int_a^b dx f(x)
を計算できるようです．楽すぎ．試しに

#!/usr/bin/env python
from scipy import integrate

def f(x):
    return x**2

print integrate.quad(f,0,1)

とかってやると普通に (0.33333333333333337, 3.7007434154171887e-15) って出ます．二番目は精度でしょう．

本番

というわけでやり方が分かったのでやってみます．

#!/usr/bin/env python
from scipy import integrate

class cosmology():
    def __init__(self,l=0.73,m=0.27,r=0,h=0.71):
        self.omegal=l
        self.omegam=m
        self.omegar=r
        self.h=h

    def dtda(self,a):
        return 9.78e9/self.h \
            *(self.omegal+self.omegam/a**3+self.omegar/a**4 \
            +(self.omegal+self.omegam+self.omegar-1)/a**2)**(-0.5)/a

    def age(self):
        return integrate.quad(self.dtda,0,1)

if __name__=="__main__":
    myuniverse=cosmology()
    print myuniverse.age()

あとからいろいろ付け足せるように class にしました．