new SimpleMC

py/BaseCosmology.py → Cosmo/BaseCosmology.py

@@ -2,53 +2,56 @@
 # parameterization of the equation of state or densities or anything.
 # However, it does know about Hubble's constant at z=0 OR the prefactor
 # c/(H0*rd) which should be fit for in the case of "rd agnostic" fits.
-# That is why you should let it declare those parameterd based on its settings
+# That is why you should let it declare those parameters based on its settings
 ##
 # However, to get the angular diameter distance you need to pass it
-# its Curvature parameter (Omega_k basically), so you need to update
-# it
+# its Curvature parameter (Omega_k basically), so you need to update it
+#
 
 
-from scipy import *
-from RadiationAndNeutrinos import *
-from Parameter import *
-from ParamDefs import *
+
+import scipy as sp
+from scipy import constants
 import scipy.integrate as intg
-from scipy.interpolate import griddata
-from scipy.integrate import odeint
-import CosmoApprox as CA
+
+from ParamDefs import h_par, Pr_par
 
 
 class BaseCosmology:
-    # speed of light
-    c_ = 299792.45
+    # speed of light in km s^-1
+    c_ = constants.c/1000.
 
     def __init__(self, h=h_par.value):
-        self.Curv = 0
-        self.rd = 149.50
-        self.h = h
+        self.Curv    = 0
+        self.rd      = 149.50
+        self.h       = h
         self.prefact = Pr_par.value
         self.varyPrefactor = False
         BaseCosmology.updateParams(self, [])
 
     def setCurvature(self, R):
         self.Curv = R
 
+
     def setrd(self, rd):
         self.rd = rd
 
+
     def setVaryPrefactor(self, T=True):
         self.varyPrefactor = T
 
+
     def setPrefactor(self, p):
         self.prefact = p
 
+
     def prefactor(self):
         if self.varyPrefactor:
             return self.prefact
         else:
             return self.c_/(self.rd*self.h*100)
 
+
     def freeParameters(self):
         if (self.varyPrefactor):
             Pr_par.setValue(self.prefact)
@@ -58,14 +61,17 @@ def freeParameters(self):
             l = [h_par]
         return l
 
+
     def printFreeParameters(self):
         print("Free parameters:")
         self.printParameters(self.freeParameters())
 
+
     def printParameters(self, params):
         for p in params:
             print(p.name, '=', p.value, '+/-', p.error)
 
+
     def updateParams(self, pars):
         for p in pars:
             if p.name == "h":
@@ -74,26 +80,32 @@ def updateParams(self, pars):
                 self.setPrefactor(p.value)
                 # h shouldn't matter here
                 # we do not want it to enter secondarily through
-                # say neutrinos, so let's keep it
-                # sane
+                # say neutrinos, so let's keep it sane
+                #
                 # self.h=p.value*self.rd*100/self.c_
         return True
 
+
     def prior_loglike(self):
         return 0
 
+
     # this is relative hsquared as a function of a
     ## i.e. H(z)^2/H(z=0)^2
     def RHSquared_a(self, a):
         print("You should not instatiate BaseCosmology")
-        error("BAD")
+        print("BAD")
+        return 0
+
 
     def Hinv_z(self, z):
-        return 1/sqrt(self.RHSquared_a(1.0/(1+z)))
+        return 1./sp.sqrt(self.RHSquared_a(1.0/(1+z)))
+
 
     # @autojit
     def DistIntegrand_a(self, a):
-        return 1/sqrt(self.RHSquared_a(a))/a**2
+        return 1./sp.sqrt(self.RHSquared_a(a))/a**2
+
 
     # @autojit
     def Da_z(self, z):
@@ -105,47 +117,51 @@ def Da_z(self, z):
         if self.Curv == 0:
             return r
         elif (self.Curv > 0):
-            q = sqrt(self.Curv)
+            q = sp.sqrt(self.Curv)
             # someone check this eq
             # Pure ADD has a 1+z fact, but have
             # comoving one
-            return sinh(r*q)/(q)
+            return sp.sinh(r*q)/(q)
         else:
-            q = sqrt(-self.Curv)
-            return sin(r*q)/(q)
+            q = sp.sqrt(-self.Curv)
+            return sp.sin(r*q)/(q)
+
 
     # D_a / rd
     def DaOverrd(self, z):
         return self.prefactor()*self.Da_z(z)
-    # H^{-1} / rd
 
+
+    # H^{-1} / rd
     def HIOverrd(self, z):
         return self.prefactor()*self.Hinv_z(z)
-    # Dv / rd
 
+
+    # Dv / rd
     def DVOverrd(self, z):
         return self.prefactor()*(self.Da_z(z)**(2./3.)*(z*self.Hinv_z(z))**(1./3.))
 
+
     # distance modulus
     def distance_modulus(self, z):
-
         # I think this should also work with varyPrefactor as long as BAO is there too
         # assert(not self.varyPrefactor)
 
         # note that our Da_z is comoving, so we're only
         # multilpyting with a single (1+z) factor
-        return 5*log10(self.Da_z(z)*(1+z))
+        return 5*sp.log10(self.Da_z(z)*(1+z))
 
-    # returns the growth factor as a function of redshift
 
+    # returns the growth factor as a function of redshift
     def GrowthIntegrand_a(self, a):
-        return 1/(self.RHSquared_a(a)*a*a)**(1.5)
+        return 1./(self.RHSquared_a(a)*a*a)**(1.5)
+
 
     def growth(self, z):
         # Equation 7.80 from Doddie
         af = 1/(1.+z)
         r = intg.quad(self.GrowthIntegrand_a, 1e-7, af)
-        gr = sqrt(self.RHSquared_a(af))*r[0]  # assume precision is ok
+        gr = sp.sqrt(self.RHSquared_a(af))*r[0]  # assume precision is ok
         # If we have Omega_m, let's normalize that way
         if hasattr(self, "Om"):
             gr *= 5/2.*self.Om

py/CosmoApprox.py → Cosmo/CosmoApprox.py

@@ -2,6 +2,7 @@
 # and rd calculators in the second section
 ##
 
+import sys
 import math as N
 Tcmb = 2.7255
 
@@ -72,14 +73,15 @@ def rd_anderson_approx(obh2, ocbh2, onuh2, Nnu):
     if (abs(Nnu-3) > 0.1):
         print("ERROR, cannot use anderson approx with Nnu")
         print("Nnu=", Nnu)
+        sys.exit(1)
     return 55.234 / (ocbh2**0.2538 * obh2**0.1278 * (1+onuh2)**0.3794)
 
 
 def rd_cuesta_approx(obh2, ocbh2, onuh2, Nnu):
     if (abs(Nnu-3) > 0.1):
         print("ERROR, Tony Cuesta says: 'not in this ceral box.'")
         print("Nnu=", Nnu)
-        stop()
+        sys.exit(1)
     return 55.154 / (ocbh2**0.25351 * (obh2)**0.12807 * N.exp(
         (onuh2+0.0006)**2.0/0.1176**2))
 
@@ -92,5 +94,5 @@ def rd_EH_approx(obh2, ocbh2, onuh2, Nnu):
     if (abs(Nnu-3) > 0.1):
         print("ERROR, cannot use EH approx with Nnu.")
         print("Nnu=", Nnu)
-        stop()
+        sys.exit(1)
     return self.CA.soundhorizon_eh(ocbh2, obh2, Nnu)

py/NuDensity.py → Cosmo/NuDensity.py

@@ -3,40 +3,44 @@
 # of neutrino energy densities.
 #
 
-from scipy import *
+
+import scipy as sp
+from scipy import constants as ct
+from scipy.special import zeta
 from scipy.integrate import quad
 from scipy.interpolate import interp1d
-from CosmoApprox import *
+import CosmoApprox as CA
 #from numba import autojit
 
 
 class NuIntegral:
     def __init__(self):
         "this integrates the stupid integral"
         print("Initalizing nu density look up table...", end=' ')
-        rat = 10**(arange(-4, 5, 0.1))
+        rat = 10**(sp.arange(-4, 5, 0.1))
         intg = []
         for r in rat:
             # min below is to supress the stupid overlow warning
-            res = quad(lambda x: sqrt(x*x+r**2) /
-                       (exp(min(x, 400))+1.0)*x**2, 0, 1000)
+            res = quad(lambda x: sp.sqrt(x*x+r**2) /
+                       (sp.exp(min(x, 400))+1.0)*x**2, 0, 1000)
             intg.append(res[0]/(1+r))
-        intg = array(intg)
+        intg = sp.array(intg)
         intg *= 7/8./intg[0]  # the right normalization
 
-        self.interpolator = interp1d(log(rat), intg)
+        self.interpolator = interp1d(sp.log(rat), intg)
         # type this into maple
-        # evalf(45*Zeta(3)/(2*Pi^4));
-        self.int_infty = 0.2776566337
+        # evalf(45*Zeta(3)/(2*Pi^4));  0.2776566337
+        self.int_infty = 45*zeta(3)/(2*ct.pi**4)
         print("Done")  # self.int_infty,intg[-1]*(1+r)/r
 
+
     def SevenEights(self, mnuOT):
         if (mnuOT < 1e-4):
             return 7/8.
         elif (mnuOT > 1e4):
-            return self.int_infty*mnuOT  # I don't think this every matters
+            return self.int_infty*mnuOT  # I don't think this ever matters
         else:
-            return self.interpolator(log(mnuOT))*(1+mnuOT)
+            return self.interpolator(sp.log(mnuOT))*(1+mnuOT)
 
 
 class ZeroNuDensity:
@@ -63,45 +67,49 @@ def __init__(self, TCMB, Nnu=3.046, mnu=0.06, degenerate=False):
         # one neutrino worth of radiation at the nominal temperature, or heated on?
         # See CAMB notes Eq. 4-7.
 
-        self.gfact = (3.046/3.0)
+        self.gfact    = (3.046/3.0)
         self.gfact_o4 = self.gfact**(0.25)
         # ideal neutrino temp
-        self.Tnu0 = (4./11.)**(1./3.)*TCMB
+        self.Tnu0     = (4./11.)**(1./3.)*TCMB
         # actual neutrino temp
-        self.Tnu = self.Tnu0*self.gfact_o4
+        self.Tnu      = self.Tnu0*self.gfact_o4
         # same for prefactors
-        self.prefix0 = 4.48130979e-7 * TCMB**4 * ((4./11.)**(4./3.))
-        self.prefix = self.prefix0*self.gfact
+        self.prefix0  = 4.48130979e-7 * TCMB**4 * ((4./11.)**(4./3.))
+        self.prefix   = self.prefix0*self.gfact
 
         self.degenerate = degenerate
         self.set_mnuone_()
 
+
     def set_mnuone_(self):
         if self.degenerate:
-            self.mnuone = self.mnu_/self.Nnu_
+            self.mnuone  = self.mnu_/self.Nnu_
         else:
-            self.mnuone = self.mnu_*1.0
+            self.mnuone  = self.mnu_*1.0
         self.omnuh2today = self.rho(1)
 
+
     def setMnu(self, mnu):
         self.mnu_ = mnu
         self.set_mnuone_()
 
+
     def setNnu(self, Nnu):
         self.Nnu_ = Nnu
         self.set_mnuone_()
 
-    # @autojit
 
+    # @autojit
     def rho(self, a):
         # This returns the density at a normalized so that
         # we get nuh2 at a=0
-        # (1 eV) / (Boltzmann constant * 1 kelvin) =
-        # 11 604.5193
+        # (1 eV) / (Boltzmann constant * 1 kelvin) = 11 604.5193
+
         if (self.mnuone == 0):
             return self.Nnu_*7/8.*self.prefix0/a**4
 
-        mnuOT = self.mnuone/(self.Tnu/a)*11604.5193
+        mnuOT = self.mnuone/(self.Tnu/a)*(1./ct.value(u'Boltzmann constant in eV/K')) #11604.5193
+
         # Here for massive we use 1*prefix (accounting for 1.015 in Tnu)
         # For massles we use Neff*prefix0 (so we account
         if self.degenerate:
@@ -110,16 +118,17 @@ def rho(self, a):
             return ((self.I.SevenEights(mnuOT)*self.prefix+(self.Nnu_-1.015)*7/8.*self.prefix0))/a**4
 
 
+
 if __name__ == "__main__":
-    A = NuDensity(Tcmb, 3.046, 0.06)
+    A = NuDensity(CA.Tcmb, 3.046, 0.06)
     print(A.omnuh2today, '=including massless neutrinos')
-    B = NuDensity(Tcmb, 2.030, 0.00)
+    B = NuDensity(CA.Tcmb, 2.030, 0.00)
     print(B.omnuh2today, end=' ')
     print(A.omnuh2today-B.omnuh2today, '=excluding massless neutrinos')
 
-    A = NuDensity(Tcmb, 1.015, 60.)
+    A = NuDensity(CA.Tcmb, 1.015, 60.)
     print(A.omnuh2today/1000., '=assuming very cold')
-    A = NuDensity(Tcmb, 1.015, 0.06)
+    A = NuDensity(CA.Tcmb, 1.015, 0.06)
 
     print(A.omnuh2today, '=assuming real temperature')