using astro cosmo and constants

models/galaxy_model/mergerrate.py

@@ -3,13 +3,15 @@ import numpy.random as nr
 import scipy.optimize as so
 import scipy.integrate as si
 import matplotlib.pyplot as plt
+from astropy import constants as const
+import astropy.units as u
 
 #cosmological functions begin
 
 def invE(z):
-	OmegaM = 0.3
-	Omegak = 0.
-	OmegaLambda = 0.7
+	OmegaM = cosmo.Om0 #0.3
+	Omegak = cosmo.Ok0 #0.
+	OmegaLambda = cosmo.Ode0 #0.7
 	return 1./np.sqrt(OmegaM*(1.+z)**3.+Omegak*(1.+z)**2.+OmegaLambda)
 
 def dtdz(z):
@@ -19,13 +21,13 @@ def dtdz(z):
 	return t0/(1.+z)*invE(z)
 
 def dVcdz(z):
-	c = 3.e8
-	H0 = 70.e3
+	c = const.c.value  # 3.e8
+	H0 = cosmo.H(0).value*1.e3 # 70.e3 -> check with Siyuan
 	return 4.*np.pi*c/H0*invE(z)*DM(z)*DM(z)
 
 def DM(z):
-	c = 3.e8
-	H0 = 70.e3
+	c = const.c.value   # 3.e8
+	H0 = cosmo.H(0).value*1.e3 # 70.e3 -> check with Siyuan
 	return c/H0*si.quad(invE, 0., z)[0]
 
 def mchirpq(q):
@@ -49,11 +51,11 @@ def mfractionp(m):
 #cosmological functions end
 #spectrum computation functions begin
 
-Gu = 6.673848e-11
-cu = 2.997925e8
-Msun = 1.98855e30
-pc = 3.0856776e16
-Gyr = 3.15576e16
+Gu = const.G.value                # 6.673848e-11
+cu = const.G.value                # 2.997925e8
+Msun = const.M_sun.value          # 1.98855e30
+pc = 1.*u.parsec.to(u.m)          # 3.0856776e16
+Gyr = 1.e9 * 1.*u.year.to(u.s)    # 3.15576e16
 G = Gu*Msun/pc**3.0
 c = cu/pc
 
@@ -401,8 +403,8 @@ class mergerrate(object):
 	def hdrop(self,n0,f,fbin):
 		if fbin is None:
 			return n0
-		c = 3.e8
-		G = 1.33e20
+		c = const.c.value   # 3.e8
+		G = 1.33e20 # should we use const.G.value here too (not sure what the unit is here)
 		fhigh = f + fbin
 		flow = f - fbin
 		Mcbh = np.linspace(5,11,30)