setting the possibility to vary the energy injection as we want

.gitignore

@@ -9,4 +9,6 @@ build
 __pycache__
 pyhyrec.egg-info
 
-dist
+dist
+
+example_pmf.ipynb

README.md

@@ -32,7 +32,7 @@ In this fork we have implement a lightweight python wrapper for HYREC-2 using cy
   pip install hyrec
   ```
 
-## How to run HYREC-2
+### How to run HYREC-2
 
 An example of notebook using the python wrapper is given [here](https://github.com/gaetanfacchinetti/HYREC-2/blob/master/example.ipynb).
 

example.py

@@ -2,7 +2,55 @@
 import numpy as np
 import pyhyrec as pyhy
 
-inj_params = pyhy.init_INPUT_INJ_PARAMS(0., 0., 0., 0., 0., 0., 0., 0, 0., 1., 0.)
-cosmo_params = pyhy.init_INPUT_COSMOPARAMS(6.735837e-01, 2.7255, 0.0494142797907188, 0.31242079216478097, 0., -1, 0, 3.046, 1.0, 0.06, 0., 0., 0.245, 1., 1.)
+cosmo = pyhy.HyRecCosmoParams()
+inj  = pyhy.HyRecInjectionParams()
 
-z, xe, Tm = pyhy.call_run_hyrec(cosmo_params, inj_params)
+z, xe, Tm = pyhy.call_run_hyrec(cosmo(), inj())
+
+# load data for comparisons to the precomputed default case
+data = np.loadtxt("./tests/output_xe.dat")
+
+z_c = np.flip(data[:, 0])
+xe_c = np.flip(data[:, 1])
+Tm_c = np.flip(data[:, 2])
+
+
+# create a cosmo params and injection params object with Primordial Black Holes
+inj_pbh  = pyhy.HyRecInjectionParams({'fpbh' : 1.0, 'Mpbh': 1e+3})
+z_pbh, xe_pbh, Tm_pbh = pyhy.call_run_hyrec(cosmo(), inj_pbh())
+
+
+# plotting xe
+fig = plt.figure()
+ax = fig.gca()
+
+ax.plot(z, xe, '-b', label='without pbh')
+ax.plot(z_pbh, xe_pbh, '-r', label='with pbh')
+ax.plot(z_c, xe_c, 'c-.',  label='without pbh (default)')
+
+ax.set_xscale("log")
+ax.set_yscale("log")
+ax.set_ylabel('xe')
+ax.set_xlabel('z')
+ax.legend()
+
+fig.savefig("../test_xe.pdf")
+
+
+# plotting Tm
+fig = plt.figure()
+ax = fig.gca()
+
+ax.plot(z, Tm, '-b', label='without pbh')
+ax.plot(z_pbh, Tm_pbh, '-r', label='with pbh')
+ax.plot(z_c, Tm_c, 'c-.',  label='without pbh (default)')
+ax.plot(z, cosmo.T0 * (1+z), color='k', linestyle=':', label='CMB temp.')
+
+ax.set_xscale("log")
+ax.set_yscale("log")
+ax.set_ylabel('Tm [K]')
+ax.set_xlabel('z')
+ax.legend()
+
+
+fig.savefig("../test_TK.pdf")

pyproject.toml

@@ -4,6 +4,7 @@
 requires = ["setuptools>=61.0.0", "wheel", "cython"]
 build-backend = "setuptools.build_meta"
 
+
 [project]
 name = "pyhyrec"
 version = "1.0.0"
@@ -34,4 +35,5 @@ PythonVersion = "https://github.com/gaetanfacchinetti/HYREC-2"
 
 [tool.pytest.ini_options]
 minversion = "6.0"
-testpaths = ["tests"]
+testpaths = ["tests"]
+

setup.py

@@ -1,28 +1,39 @@
 import subprocess
 import os
 
-# go to where the project needs to be compiled
-os.chdir("./src/pyhyrec/src/")
+from setuptools import setup, Extension
+from Cython.Build import cythonize
 
-try:
-    subprocess.run(['make'], check=True)
-    print("make executed successfully.")
+def main():
 
-    # Run the 'make' command with the desired target
-    subprocess.run(['make', 'clean'], check=True)
-    print("make clean executed successfully.")
+    # go to where the project needs to be compiled
+
+    os.chdir("./src/pyhyrec/src/")
 
-except subprocess.CalledProcessError as e:
-    print("Error executing Makefile:", e)
+    try:
+
+        if os.path.exists("libhyrec.so"):
+            os.remove("libhyrec.so")
 
-# go back to the main project folder
-os.chdir("../../../")
+        subprocess.run(['make'], check=True)
+        print("make executed successfully.")
 
-from setuptools import setup, Extension
-from Cython.Build import cythonize
+        # Run the 'make' command with the desired target
+        subprocess.run(['make', 'clean'], check=True)
+        print("make clean executed successfully.")
+
+    except subprocess.CalledProcessError as e:
+        print("Error executing Makefile:", e)
 
-cythonize("./src/pyhyrec/wrapperhyrec.pyx", force = True)
+    # go back to the main project folder
+    os.chdir("../../../")
+
+    #extensions = [Extension("pyhyrec.src.hyrec", sources = ["./src/pyhyrec/src/hyrec.c"])]
+    extensions = [Extension("pyhyrec.wrapperhyrec", sources = ["./src/pyhyrec/wrapperhyrec.c"], libraries = ["hyrec"], library_dirs=["./src/pyhyrec/src/"], build_temp="src/pyhyrec/")]
+    cythonize("./src/pyhyrec/wrapperhyrec.pyx", force = True)
+
+    setup(ext_modules=extensions, package_dir={'' : 'src'},  packages=['pyhyrec'], package_data= {'pyhyrec' :['data/*']}, include_package_data=True)
 
-extensions = [Extension("pyhyrec.wrapperhyrec", sources = ["./src/pyhyrec/wrapperhyrec.c"], libraries = ["hyrec"], library_dirs=["./src/pyhyrec/src/"], build_temp="src/pyhyrec/")]
 
-setup(ext_modules=extensions, package_dir={'': 'src'},  packages=['pyhyrec'], package_data= {'pyhyrec' :['data/*']}, include_package_data=True)
+if __name__ == "__main__":
+    main()

src/pyhyrec.egg-info/PKG-INFO

@@ -728,7 +728,7 @@ In this fork we have implement a lightweight python wrapper for HYREC-2 using cy
   pip install hyrec
   ```
 
-## How to run HYREC-2
+### How to run HYREC-2
 
 An example of notebook using the python wrapper is given [here](https://github.com/gaetanfacchinetti/HYREC-2/blob/master/example.ipynb).
 

src/pyhyrec.egg-info/SOURCES.txt

@@ -4,11 +4,13 @@ README.md
 pyproject.toml
 setup.py
 ./src/pyhyrec/wrapperhyrec.c
+./src/pyhyrec/src/hyrec.c
 src/pyhyrec/.DS_Store
 src/pyhyrec/__init__.py
 src/pyhyrec/__main__.py
 src/pyhyrec/params.py
 src/pyhyrec/wrapperhyrec.c
+src/pyhyrec/wrapperhyrec.cpython-310-darwin.so
 src/pyhyrec/wrapperhyrec.pyx
 src/pyhyrec.egg-info/PKG-INFO
 src/pyhyrec.egg-info/SOURCES.txt
@@ -21,19 +23,4 @@ src/pyhyrec/data/__init__.py
 src/pyhyrec/data/fit_swift.dat
 src/pyhyrec/data/two_photon_tables.dat
 src/pyhyrec/data/two_photon_tables_hires.dat
-src/pyhyrec/src/Makefile
-src/pyhyrec/src/energy_injection.c
-src/pyhyrec/src/energy_injection.h
-src/pyhyrec/src/helium.c
-src/pyhyrec/src/helium.h
-src/pyhyrec/src/history.c
-src/pyhyrec/src/history.h
-src/pyhyrec/src/hydrogen.c
-src/pyhyrec/src/hydrogen.h
-src/pyhyrec/src/hyrec
-src/pyhyrec/src/hyrec.c
-src/pyhyrec/src/hyrectools.c
-src/pyhyrec/src/hyrectools.h
-src/pyhyrec/src/input.dat
-src/pyhyrec/src/libhyrec.so
 tests/test.py

src/pyhyrec/__init__.py

@@ -1,3 +1,8 @@
-from .wrapperhyrec import (init_INPUT_COSMOPARAMS, init_INPUT_INJ_PARAMS, call_run_hyrec)
+from .wrapperhyrec import (
+    init_INPUT_COSMOPARAMS, 
+    init_INPUT_INJ_PARAMS, 
+    call_run_hyrec,
+    call_decay_rate_pmf_turbulences,
+    )
 
 from .params import HyRecCosmoParams, HyRecInjectionParams

src/pyhyrec/src/energy_injection.c

@@ -66,24 +66,25 @@ double dEdtdV_DM_ann(double z, INJ_PARAMS *params){
 
 /***************************************************************************************
 Effect of primordial magnetic fields
-According to Chluba et al. 2015
+According to Kunze and Komatzu 2014
 ***************************************************************************************/
 
-// Energy injection rate due to turbulences
-// H must be in 1/s for a result in eV / cm^3 / s 
-double gamma_turbulences(double z, double H, double B0, double nB)
+/* Mode for a magnetic field with variance sigB0 (nG) on scale lambda (Mpc), normalisation sig, and power index nB */
+double pmf_mode(double sigB0, double nB, double sig, double lambda)
+{
+  return 2 * M_PI * pow(pow(sigB0 / sig, 2) * pow(lambda / 2.0 / M_PI, 3.0 + nB), -1.0/(5.0+nB));
+}
+
+// Energy injection rate due to turbulences in units of Hubble time
+double decay_rate_pmf_turbulences(double z, double tdti, double nB)
 {
   double zi = 1088;
 
   if (z > zi)
     return 0;
 
   double m = 2.0*(nB+3.0)/(nB + 5.0);
-  double kd = 286.91 * B0;
-  double titd = 14.8 / B0 / kd;
-  double rhoB = 9.5e-8 * square(B0) * 0.26 * pow(1+z, 4);
-
-  return 3.0*m/2.0 * pow(log(1+titd), m)/pow(log(1+titd) + 1.5 * log((1+zi)/(1+z)), m+1) * H * rhoB;
+  return 3.0*m/2.0 * pow(log(1+tdti), m)/pow(log(1+tdti) + 1.5 * log((1+zi)/(1+z)), m+1);
 }
 
 
@@ -233,49 +234,4 @@ double dEdtdV_inj(double z, double xe, double Tgas, INJ_PARAMS *params){
 }
 
 
-/**********************************************************************************
-Energy *deposition* rate per unit volume
-Essentially use a very simple ODE solution
-**********************************************************************************/
-
-void update_dEdtdV_dep(double z_out, double dlna, double xe, double Tgas,
-		       double nH, double H, INJ_PARAMS *params, double *dEdtdV_dep) {
-
-  double inj  = dEdtdV_inj(z_out, xe, Tgas, params);
-
-  if (params->on_the_spot == 1) *dEdtdV_dep = inj;
-
-  // Else put in your favorite recipe to translate from injection to deposition
-  // Here I assume injected photon spectrum Compton cools at rate dE/dt = - 0.1 n_h c sigma_T E
-  // This is valid for E_photon ~ MeV or so.
-
-  else { // c sigma_T = 2e-14 (cgs)
-    *dEdtdV_dep = (*dEdtdV_dep *exp(-7.*dlna) + 2e-15* dlna*nH/H *inj)
-                 /(1.+ 2e-15 *dlna*nH/H);                              
-  }       
-}
-
-/*******************************************************************************
-Fraction of energy deposited in the form of heat, ionization and excitations
-*******************************************************************************/
-
-double chi_heat(double xe) { 
-  return (1.+2.*xe)/3.; // Approximate value of Chen & Kamionkowski 2004 
-
-  // fit by Vivian Poulin of columns 1 and 2 in Table V of Galli et al. 2013
-  // overall coefficient slightly changed by YAH so it reaches exactly 1 at xe = 1.  
-  // return (xe < 1.? 1.-pow(1.-pow(xe,0.300134),1.51035) : 1.);
-}
-
-double chi_ion(double xe) {
-  return (1.-xe)/3.; // Approximate value of Chen & Kamionkowski 2004 
-
-  // fit by Vivian Poulin of columns 1 and 4 in Table V of Galli et al. 2013
-  // return 0.369202*pow(1.-pow(xe,0.463929),1.70237); 
-}
-
-double chi_exc(double xe) {
-  return 1. - chi_ion(xe) - chi_heat(xe);
-}
-
 

src/pyhyrec/src/energy_injection.h

@@ -22,11 +22,12 @@ typedef struct {
   int on_the_spot;            /* if set to 1 assume energy deposition rate = injection rate */
                               /* Otherwise solves for deposition given injection with simple recipe */
 
-  double ion, exclya;
+  double ion, exclya, dEdtdV_heat;   /* Adding the possibility to have a heating decorrelated from ion and exclya */
 
 } INJ_PARAMS;
 
-void update_dEdtdV_dep(double z_out, double dlna, double xe, double Tgas,
-                       double nH, double H, INJ_PARAMS *params, double *dEdtdV_dep);
+
+double dEdtdV_inj(double z, double xe, double Tgas, INJ_PARAMS *params);
+double decay_rate_pmf_turbulences(double z, double tdti, double nB); /* decay rate of PMF turbulences*/
 
 #endif