tidy

pyhmcode/constants.py

@@ -8,7 +8,7 @@
 Mpc = 3.0857e16*1e6  # Mpc [m]
 
 # Cosmology
-H0 = 100. # Hubble parameter today in h [km/s/Mpc]
-G_cosmological = G*Sun_mass/(Mpc*1e3**2) # Gravitational constant [(Msun/h)^-1 (km/s)^2 (Mpc/h)] ~4.301e-9 (1e3**2 m -> km)
+H0 = 100.                                      # Hubble parameter today in h [km/s/Mpc]
+G_cosmological = G*Sun_mass/(Mpc*1e3**2)       # Gravitational constant [(Msun/h)^-1 (km/s)^2 (Mpc/h)] ~4.301e-9 (1e3**2 m -> km)
 rho_critical = 3.*H0**2/(8.*pi*G_cosmological) # Critical density [(Msun/h) (Mpc/h)^-3] ~2.775e11 (Msun/h)/(Mpc/h)^3
-neutrino_constant = 93.1 # Neutrino mass required to close universe [eV]
+neutrino_constant = 93.1                       # Neutrino mass required to close universe [eV]

pyhmcode/cosmology.py

@@ -12,7 +12,7 @@
 # Parameters
 xmin_Tk = 1e-5    # Scale at which to switch to Taylor expansion approximation in tophat Fourier functions
 eps_sigmaR = 1e-4 # Accuracy of the sigmaR integration
-eps_sigmaV = 1e-4 # Accuracy of the sigmaV integration NOTE: Seems to fail with higher accuracy (1e-4; when zs=[0] only?!?)
+eps_sigmaV = 1e-4 # Accuracy of the sigmaV integration
 eps_neff = 1e-4   # Accuracy of the neff integration (d ln sigma^2/d ln k)
 
 ### Backgroud ###
@@ -147,9 +147,9 @@ def sigmaV(R:float, Pk:callable, kmin=0., kmax=np.inf, eps=eps_sigmaV, transform
     '''
     Integration to get the 1D RMS in the linear displacement field
     TODO: This generates a warning sometimes, and is slow, there must be a cleverer way to integrate here.
-    Unless eps_sigmaV > 1e-3 the integration fails for z=0 sometimes, but not after being called for
-    z > 0. I really don't understand this, but it's annoying and should be fixed (kmin /= 0. helps).
-    I should look at how CAMB deals with these type of integrals (e.g., sigmaR).
+    TODO: Unless eps_sigmaV > 1e-3 the integration fails for z=0 sometimes, but not after being called for z > 0.
+    TODO: I really don't understand this, but it's annoying and should be fixed (kmin /= 0. helps).
+    TODO: I should look at how CAMB deals with these type of integrals (e.g., sigmaR).
     args:
         Pk: Function of k to evaluate the linear power spectrum
         eps: Integration accuracy

pyhmcode/linear_growth.py

@@ -68,13 +68,14 @@ def get_growth_interpolator(CAMB_results:camb.CAMBdata, LCDM=False) -> callable:
     Solve the linear growth ODE and returns an interpolating function for the solution
     LCDM = True forces w = -1 and imposes flatness by modifying the dark-energy density
     TODO: w dependence for initial conditions; f here is correct for w=0 only
+    TODO: Could use d_init = a(1+(w-1)/(w(6w-5))*(Om_w/Om_m)*a**-3w) at early times with w = w(a<<1)
     '''
     from scipy.integrate import solve_ivp
     from scipy.interpolate import interp1d as interp
     na = 129 # Number of scale factors used to construct interpolator
     a = np.linspace(a_init, 1., na)
     f = 1.-_Omega_m(a_init, CAMB_results, LCDM=LCDM) # Early mass density
-    d_init = a_init**(1.-3.*f/5.) # Initial condition (~ a_init; but f factor accounts for EDE-ish)
+    d_init = a_init**(1.-3.*f/5.)            # Initial condition (~ a_init; but f factor accounts for EDE-ish)
     v_init = (1.-3.*f/5.)*a_init**(-3.*f/5.) # Initial condition (~ 1; but f factor accounts for EDE-ish)
     y0 = (d_init, v_init)
     def fun(a, y):