growth function updates

simonsobs · Aug 12, 2024 · 75c92c3 · 75c92c3
1 parent 630c877
commit 75c92c3
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 10 deletions.
diff --git a/examples/growth.py b/examples/growth.py
@@ -26,20 +26,22 @@
 avals = 1./(1+zs)
 Dmd = avals
 
+pl = io.Plotter('rCL',xlabel='$z$',ylabel=r'$\Delta D(z) / D(z)$',xyscale='loglin')
+
 for engine in ['camb','class']:
-    pl = io.Plotter('rCL',xlabel='$z$',ylabel=r'$\Delta D(z) / D(z)$',xyscale='loglin')
 
-    for model in models:
+    for i,model in enumerate(models):
         params = model[1]
         h = hcosmo.Cosmology(params,accuracy='low',engine=engine)
         label = model[0]
         dapprox = h.D_growth(avals,exact=False)
         dexact = h.D_growth(avals,exact=True)
         ddiff = (dapprox-dexact)/dexact
-        pl.add(zs,ddiff,label=label)
+        pl.add(zs,ddiff,label=label if engine=='camb' else None,ls={'camb':'-','class':'--'}[engine],color=[f'C{x}' for x in range(len(models))][i])
         print(label)
-    pl.hline(y=0)
-    pl._ax.set_ylim(-0.2,0.2)
-    pl.done(f'growth_{engine}.png')
+pl.hline(y=0)
+pl.legend('outside')
+pl._ax.set_ylim(-0.25,0.6)
+pl.done(f'growth.png')
 
 
diff --git a/hmvec/cosmology.py b/hmvec/cosmology.py
@@ -227,9 +227,9 @@ def get_sigma8(self,zs,exact=False):
                 return np.vectorize(lambda x : self._class_results.sigma(8./self.h,x))(zs)
         else: return np.sqrt(self.get_sigma2_R(8./self.params['H0']*100.,zs))
 
-    def D_growth_exact_arbitrary_norm(self,a):
+    def D_growth_exact_arbitrary_norm(self,a,k_camb=1e-5):
         if self.engine=='camb':
-            deltakz = self._camb_results.get_redshift_evolution(1e-5, a2z(a), ['delta_cdm']) #index: z,0
+            deltakz = self._camb_results.get_redshift_evolution(k_camb, a2z(a), ['delta_cdm']) #index: z,0
             D = deltakz[:,0]
         elif self.engine=='class':
             D = np.vectorize(self._class_results.scale_independent_growth_factor)(a2z(a))
@@ -247,6 +247,7 @@ def D_growth_approx(self,a):
         omm0 = self.omm0
         oml0 = self.oml0
         x = (oml0/omm0)**(1./3.) * a
+        # Exact analytic integral
         Dovera = np.sqrt(1.+x**3.)*(hyp2f1(5/6.,3/2.,11/6.,-x**3.))
         # An approximation to this integral
         # oma = 1./(1+x**3)
@@ -255,9 +256,9 @@ def D_growth_approx(self,a):
         return Dovera*a
 
 
-    def D_growth(self, a,type="anorm",exact=False):
+    def D_growth(self, a,type="anorm",exact=False,k_camb=1e-5):
         if exact:
-            Dfunc = self.D_growth_exact_arbitrary_norm
+            Dfunc = lambda a: self.D_growth_exact_arbitrary_norm(a,k_camb=k_camb)
         else:
             Dfunc = self.D_growth_approx
         Dtoday = Dfunc(1)