cosmology module

src/pyhyrec.egg-info/SOURCES.txt

@@ -8,6 +8,7 @@ setup.py
 src/pyhyrec/.DS_Store
 src/pyhyrec/__init__.py
 src/pyhyrec/__main__.py
+src/pyhyrec/cosmology.py
 src/pyhyrec/params.py
 src/pyhyrec/wrapperhyrec.c
 src/pyhyrec/wrapperhyrec.cpython-310-darwin.so

src/pyhyrec/__init__.py

@@ -7,4 +7,19 @@
     call_dEdtdV_heat_ambipolar_pmf,
     )
 
-from .params import HyRecCosmoParams, HyRecInjectionParams
+from .params import HyRecCosmoParams, HyRecInjectionParams
+
+from .cosmology import (
+    hubble_rate, 
+    hubble_factor, 
+    z_rec, 
+    visibility_function, 
+    optical_depth, 
+    compute_z_rec, 
+    acoustic_damping_scale, 
+    n_baryons, 
+    rho_baryons, 
+    rho_radiation, 
+    rho_gamma,
+    compute_acoustic_damping_scale,
+    )

src/pyhyrec/cosmology.py

@@ -0,0 +1,168 @@
+########################
+# Some simple functions implemented to compute simple cosmological quantities
+# from the result of HYREC runs - without invoking involved codes like CLASS
+########################
+
+import numpy as np
+from scipy import integrate
+
+from .wrapperhyrec import compute_hubble_rate, call_run_hyrec
+from .params import HyRecInjectionParams, HyRecCosmoParams
+
+
+# define usefull units and conversion rates
+_MPC_TO_M_      = 3.08567758128e+22
+_MSUN_TO_KG_    = 2.0e+30
+_KG_TO_EV_      = 5.60958860380445e+35
+_KM_TO_MPC_     = 1.0/_MPC_TO_M_ * 1e+3
+_M_TO_MPC_     =  1.0/_MPC_TO_M_
+_MASS_HYDROGEN_ = 0.93878299831e+9 # in eV
+_MASS_PROTON_   = 0.938272e+9 # in eV
+_MASS_HELIUM_   = 3.728400e+9 # in eV
+_SIGMA_THOMSON_ = 6.6524616e-29 # in m^2
+_C_LIGHT_       = 299792458 # in m / s
+
+
+def hubble_rate(z, cosmo = HyRecCosmoParams()):
+    """
+        hubble rate in s^{-1} directly computed from the C-code
+
+    Parameters:
+    -----------
+    - z: np.ndarray, float
+        redshift
+    - cosmo: HyRecCosmoParams, optional
+        cosmology
+
+    Returns:
+    --------
+    hubble rate in s^{-1}
+    """
+
+    _IS_A_LIST = True
+    if not isinstance(z, (list, np.ndarray)):
+        _IS_A_LIST = False
+        z = np.array([z])
+
+    res = np.zeros(len(z))
+
+    for iz, _z in enumerate(z):
+        res[iz] = compute_hubble_rate(_z, cosmo(), HyRecInjectionParams()())
+
+    if _IS_A_LIST == True:
+        return res
+
+    return res[0]
+
+
+def hubble_factor(z, cosmo=HyRecCosmoParams()):
+    """
+        hubble factor H(z)/H_0 (dimensionless)
+
+    Parameters:
+    -----------
+    - z: np.ndarray, float
+        redshift
+    - cosmo: HyRecCosmoParams, optional
+        cosmology
+
+    Returns:
+    --------
+    hubble factor H(z)/H0
+    """
+    return hubble_rate(z, cosmo) / 3.2407792896393e-18 / cosmo.h
+
+
+def rho_baryons(cosmo = HyRecCosmoParams()):
+    return 2.7754e+11 * cosmo.Omega_b * cosmo.h**2 * _MSUN_TO_KG_ * _KG_TO_EV_ / _MPC_TO_M_**3 # in eV / m^3
+
+def n_baryons(cosmo = HyRecCosmoParams()):
+    return rho_baryons(cosmo) / _MASS_PROTON_ / (1 + cosmo.YHe / 4 * (_MASS_HELIUM_/_MASS_HYDROGEN_ -1)) # in 1/m^3
+
+def rho_gamma(cosmo = HyRecCosmoParams()):
+    omega_gamma = 4.48162687719e-7 * cosmo.T0**4 / cosmo.h**2
+    return 2.7754e+11 * omega_gamma * cosmo.h**2 * _MSUN_TO_KG_ * _KG_TO_EV_ / _MPC_TO_M_**3 # in eV / m^3
+
+def rho_radiation(cosmo):
+    neff_array  = [3.046, 2.0328, 1.0196, 0.00641]
+    m_neutrinos = [cosmo.mnu1, cosmo.mnu2, cosmo.mnu3]
+    n_ur = neff_array[np.count_nonzero(m_neutrinos)]
+    omega_r = 4.48162687719e-7 * cosmo.T0**4 * (1. + 0.227107317660239 * n_ur) / cosmo.h**2
+
+    return 2.7754e+11 * omega_r * cosmo.h**2 * _MSUN_TO_KG_ * _KG_TO_EV_ / _MPC_TO_M_**3 # in eV / m^3
+
+def optical_depth(z, xe, cosmo = HyRecCosmoParams()):
+    """
+        optical depth assuming no reionization (dimensionless)
+
+    Parameters:
+    -----------
+    - z: np.ndarray
+        redshift
+    - xe: np.ndarray
+        free electron fraction
+    - cosmo: HyRecCosmoParams, optional
+        cosmology
+
+    Returns:
+    --------
+    optical depth for each value of z
+    """
+
+    e_z = hubble_factor(z)
+
+    # fast trapezoid integration scheme
+    integrand = xe * (1+z)**2 / e_z
+    trapz = (integrand[1:] + integrand[:-1])/2.0
+    dz = np.diff(z)
+
+    res = np.zeros(len(z))
+    for i in range(len(z)-1):
+        res[i+1] = res[i] + trapz[i] * dz[i]
+
+    pref = _C_LIGHT_ * _SIGMA_THOMSON_ * n_baryons(cosmo) / (cosmo.h * 3.2407792896393e-18)
+    return pref * res
+
+
+def visibility_function(z, xe, cosmo):
+    pref = _C_LIGHT_ * _SIGMA_THOMSON_ * n_baryons(cosmo) / (cosmo.h * 3.2407792896393e-18)
+    return pref * (1 + z)**2 / hubble_factor(z) * xe * np.exp(-optical_depth(z, xe, cosmo))
+
+def z_rec(z, xe, cosmo):
+    vis = visibility_function(z, xe, cosmo)
+    return z[np.argmax(vis)]
+
+def compute_z_rec(cosmo = HyRecCosmoParams()):
+    z, xe, _ = call_run_hyrec(cosmo(),  HyRecInjectionParams()(), zmax = 8000, zmin = 500, nz = 40000)
+    return z_rec(z, xe, cosmo)
+
+def acoustic_damping_scale(z, xe, cosmo):
+
+    # get the value of the recombination redshift
+    z_reco = z_rec(z, xe, cosmo)
+    iz_min = np.argmin(np.abs(z - z_reco)) 
+
+    # restrict the redshift and free electron fraction 
+    # to the interesting interval
+    z_int  = z[iz_min:]
+    xe_int = xe[iz_min:]
+
+    # get the baryon and photon densities
+    rho_b  = rho_baryons(cosmo) * (1+z_int)**3
+    rho_g  = rho_gamma(cosmo) * (1+z_int)**4
+    r_arr = 3./4. * rho_b / rho_g
+
+    # evaluate the Thomson scattering length
+    l_scatt = 1.0/(_SIGMA_THOMSON_ * xe_int * n_baryons(cosmo) * (1+z_int)**3 ) * _M_TO_MPC_ # in Mpc
+
+    # hubble rate in 1/Mpc
+    h_z = hubble_rate(z_int) /_M_TO_MPC_ / _C_LIGHT_ # in 1/Mpc
+
+    # damping length square
+    ld2 = integrate.trapezoid((1+z_int)/(1+ r_arr)/h_z * l_scatt * (16.0/15.0 + r_arr**2/(1+r_arr)), z_int)/6.0
+
+    return np.sqrt(1/ld2)
+
+def compute_acoustic_damping_scale(cosmo = HyRecCosmoParams()):
+    z, xe, _ = call_run_hyrec(cosmo(),  HyRecInjectionParams()(), zmax = 10000, zmin = 500, nz = 40000)
+    return acoustic_damping_scale(z, xe, cosmo)

src/pyhyrec/src/energy_injection.c

@@ -122,15 +122,24 @@ double fit_Lorentz_force_average(double x)
 }
 
 
-/* Energy injection rate due to ambipolar diffusion (result is in eV / cm^3 / s) */
+/* 
+  Energy injection rate due to ambipolar diffusion (result is in eV / cm^3 / s) 
+  - z       : redshift
+  - xe      : free electron fraction
+  - Tgas    : Temperature of the gas (in K)
+  - obh2    : \Omega_{\rm b} h^2, baryon abundance times hubble parameter squared
+  - sigma_A : alfven magnetic scale (in nG, ~30 nG fixed by the value of the Aflven mode)
+  - sigma_B : \sigma_{B, 0} (in nG),  standard deviation of the PMF spectrum on 1 Mpc scale today
+  - nB      : power index of the PMF spectrum
+*/
 double dEdtdV_heat_ambipolar_pmf(double z, double xe, double Tgas, double obh2, double sigmaA, double sigmaB, double nB, double smooth_z)
 {
 
   double zi = 1088;
 
   double gamma_AD = 6.49e-10 * pow(Tgas, 0.375) / (2.0 * mH); // in cm^3 * clight^2 / s / eV
-  double rho_b    = obh2 * 10539.850865068418; // in eV / clight^2 / cm^3
-  double eta_AD = (1.0-xe)/xe / (4*M_PI) / rho_b / rho_b / gamma_AD; // in s * clight^2 / eV * cm^3  
+  double rho_b    = obh2 * 10539.850865068418 * pow(1+z, 3); // in eV / clight^2 / cm^3
+  double eta_AD = (1.0-xe)/xe / rho_b / rho_b / gamma_AD; // in s * clight^2 / eV * cm^3  
 
   /* 
     Note that, assuming Helium ionization history similar to Hydrogen ionization history,
@@ -140,13 +149,16 @@ double dEdtdV_heat_ambipolar_pmf(double z, double xe, double Tgas, double obh2,
     rho_n = rho_HI + rho_HeI and rho_+ = rho_HII + rho_HeII (negleting second ionization)
    */
 
-  // prefactor sigma(k_A)^4 kA^2
-  double pref = 1.0502650402891526e-49 * pow(sigmaB / sigmaA, 2.0/(5.0+nB)) * pow(sigmaA, 4) * 4 * M_PI * M_PI; // in nG^4 / cm^2
+  // prefactor sigma(k_A)^4 kA^2 in convinient units
+  double sA4kA2 = 1.0502650402891526e-49 * pow(sigmaB / sigmaA, 2.0/(5.0+nB)) * pow(sigmaA, 4) * 4 * M_PI * M_PI; // in nG^4 / cm^2
 
   // mu_0 in convinient units
   double mu_0 = 4 * M_PI * 1e+19 * 1.7826619216278999e-34; // in nG^2 * cm * clight^2 * s^2 / eV 
 
-  double res =  pref * pow(1+z, 10) * eta_AD / mu_0 / mu_0 * fit_Lorentz_force_average(nB + 3.0); // in eV / s^3 / cm / clight^2 
+  // result
+  double res =  sA4kA2 * pow(1+z, 10) * eta_AD / mu_0 / mu_0 * fit_Lorentz_force_average(nB + 3.0); // in eV / s^3 / cm / clight^2 
+
+  // result in the correct output units (devide by speed of light factors)
   res = res / pow(299792458e+3, 2); // in eV / cm^3 / s 
 
   if (smooth_z > 0)

src/pyhyrec/src/history.c

@@ -133,13 +133,6 @@ void init_hyrec(REC_COSMOPARAMS * param, INPUT_COSMOPARAMS cosmo_params, INPUT_I
 HYREC_DATA * run_hyrec(INPUT_COSMOPARAMS cosmo_params, INPUT_INJ_PARAMS inj_params, double zmax, double zmin, char * table_location)
 {
 
-  //memset(&sa, 0, sizeof(struct sigaction));
-  //sigemptyset(&sa.sa_mask);
-  //sa.sa_sigaction = segfault_sigaction;
-  //sa.sa_flags   = SA_SIGINFO;
-
-  //sigaction(SIGSEGV, &sa, NULL);
-  //printf("We are here, %e \n", inj_params.fpbh);
 
   HYREC_DATA * data = malloc(sizeof(*data));
 
@@ -151,6 +144,16 @@ HYREC_DATA * run_hyrec(INPUT_COSMOPARAMS cosmo_params, INPUT_INJ_PARAMS inj_para
   return data;
 }
 
+double compute_Hubble_rate(double z, INPUT_COSMOPARAMS cosmo_params, INPUT_INJ_PARAMS inj_params)
+{
+
+  REC_COSMOPARAMS * cosmo  = (REC_COSMOPARAMS *) malloc(sizeof(REC_COSMOPARAMS));
+  cosmo->inj_params = (INJ_PARAMS *)  malloc(sizeof(INJ_PARAMS));
+  init_hyrec(cosmo, cosmo_params, inj_params);
+  double res = rec_HubbleRate(cosmo, z);
+  hyrec_free_cosmo(cosmo);
+  return res;
+}
 
 
 /*************************************************************************************
@@ -950,6 +953,11 @@ void hyrec_allocate(HYREC_DATA * data, double zmax, double zmin) {
 }
 
 
+void hyrec_free_cosmo(REC_COSMOPARAMS *cosmo){
+  free(cosmo->inj_params);
+  free(cosmo);
+}
+
 void hyrec_free(HYREC_DATA *data) {
   free_atomic(data->atomic);
   free(data->cosmo->inj_params);

src/pyhyrec/src/history.h

@@ -68,9 +68,11 @@ char* rec_build_history(HYREC_DATA *data, int model, double *hubble_array);
 
 void init_hyrec(REC_COSMOPARAMS * param, INPUT_COSMOPARAMS cosmo_params, INPUT_INJ_PARAMS injection_params);
 HYREC_DATA * run_hyrec(INPUT_COSMOPARAMS global_cosmo, INPUT_INJ_PARAMS inj_params, double zmax, double zmin, char * table_location);
+double compute_Hubble_rate(double z, INPUT_COSMOPARAMS cosmo_params, INPUT_INJ_PARAMS inj_params);
 //----------------------------------------------------//
 
 void hyrec_allocate(HYREC_DATA *data, double zmax, double zmin);
+void hyrec_free_cosmo(REC_COSMOPARAMS *cosmo);
 void hyrec_free(HYREC_DATA *data);
 void hyrec_compute(HYREC_DATA *data, int model);
 double hyrec_xe(double z, HYREC_DATA *data);