diff --git a/tools/create_perturb_file.py b/tools/create_perturb_file.py
index 801551912cdbb6083053312d354d20b193c1642c..52f32ca1235eff5a263593aba0109d2c85c84f2a 100755
--- a/tools/create_perturb_file.py
+++ b/tools/create_perturb_file.py
@@ -7,20 +7,20 @@ import numpy as np
import h5py
from classy import Class
-#Name of the perturbations file to be create
+# Name of the perturbations file to be create
fname = "perturbations.hdf5"
-#Open the file
+# Open the file
f = h5py.File("perturbations.hdf5", mode="w")
-#Cosmological parameters
+# Cosmological parameters
h = 0.681
Omega_b = 0.0486
Omega_cdm = 0.2560110606
A_s = 2.0993736148e-09
n_s = 0.967
-#Neutrino and radiation parameters
+# Neutrino and radiation parameters
T_cmb = 2.7255
T_ncdm = 0.71611
N_ur = 2.0308
@@ -28,43 +28,44 @@ N_ncdm = 1
deg_ncdm = [1]
m_ncdm = [0.06]
-#Maximum wavenumber and redshift
-kmax = 30.
+# Maximum wavenumber and redshift
+kmax = 30.0
zmax = 1e3
amin = 1.0 / (zmax + 1)
-#CLASS output distance unit
+# CLASS output distance unit
Mpc_cgs = 3.085677581282e24
-#CLASS parameters
+# CLASS parameters
params = {
- "h": h,
- "Omega_b": Omega_b,
- "Omega_cdm": Omega_cdm,
- "T_cmb": T_cmb,
- "N_ncdm": N_ncdm,
- "N_ur": N_ur,
- "T_ncdm": T_ncdm,
- "deg_ncdm": "".join(str(x)+"," for x in deg_ncdm)[:-1],
- "m_ncdm": "".join(str(x)+"," for x in m_ncdm)[:-1],
- "A_s": A_s,
- "n_s": n_s,
- "output": "dTk, vTk",
- "z_max_pk": zmax,
- "P_k_max_1/Mpc": kmax,
- "reio_parametrization": "reio_none",
- "YHe": "BBN",
- "k_output_values": kmax,
- "k_per_decade_for_pk": 100}
+ "h": h,
+ "Omega_b": Omega_b,
+ "Omega_cdm": Omega_cdm,
+ "T_cmb": T_cmb,
+ "N_ncdm": N_ncdm,
+ "N_ur": N_ur,
+ "T_ncdm": T_ncdm,
+ "deg_ncdm": "".join(str(x) + "," for x in deg_ncdm)[:-1],
+ "m_ncdm": "".join(str(x) + "," for x in m_ncdm)[:-1],
+ "A_s": A_s,
+ "n_s": n_s,
+ "output": "dTk, vTk",
+ "z_max_pk": zmax,
+ "P_k_max_1/Mpc": kmax,
+ "reio_parametrization": "reio_none",
+ "YHe": "BBN",
+ "k_output_values": kmax,
+ "k_per_decade_for_pk": 100,
+}
print("Running CLASS")
-#Run CLASS
+# Run CLASS
model = Class()
model.set(params)
model.compute()
-#Extract wavenumbers and prepare redshifts
+# Extract wavenumbers and prepare redshifts
k = model.get_transfer(0)["k (h/Mpc)"] * h
a = np.exp(np.arange(0, np.log(amin), -0.01))[::-1]
z = 1.0 / a - 1.0
@@ -78,29 +79,29 @@ keys = model.get_transfer(0).keys()
print("Available transfer functions:")
print(keys)
-#Prepare dictionary
+# Prepare dictionary
pt = {}
for key in keys:
- pt[key] = np.zeros((nz, nk))
+ pt[key] = np.zeros((nz, nk))
-#Extract transfer functions
+# Extract transfer functions
for i in range(nz):
- pti = model.get_transfer(z[i])
- for key in pti:
- pt[key][i,:] = pti[key]
+ pti = model.get_transfer(z[i])
+ for key in pti:
+ pt[key][i, :] = pti[key]
-#Export the perturbations file
+# Export the perturbations file
f.create_group("Functions")
f["Redshifts"] = z
f["Wavenumbers"] = k
f.create_group("Units")
f["Units"].attrs["Unit length in cgs (U_L)"] = Mpc_cgs
-#Write the perturbations
+# Write the perturbations
for key in keys:
- f["Functions/" + key.replace("/", "\\")] = pt[key]
+ f["Functions/" + key.replace("/", "\\")] = pt[key]
-#Close the file
+# Close the file
f.close()
print("Done.")
diff --git a/tools/spawn_neutrinos.py b/tools/spawn_neutrinos.py
index 77449b2d851932aa10f967495e1aa52cc07b223e..4f454a0acf8ed46d3f86a8774964f23f2bda70c3 100755
--- a/tools/spawn_neutrinos.py
+++ b/tools/spawn_neutrinos.py
@@ -12,9 +12,9 @@ import sys
# Usage: ./spawn_neutrinos.py filename
# Constants
-Mpc_cgs = 3.08567758e+24
-Default_N_nu_per100Mpc = 72 # 72^3 neutrinos for a (100 Mpc)^3 box
-Default_nr_neutrinos_per_Mpc3 = (Default_N_nu_per100Mpc / 100.)**3
+Mpc_cgs = 3.08567758e24
+Default_N_nu_per100Mpc = 72 # 72^3 neutrinos for a (100 Mpc)^3 box
+Default_nr_neutrinos_per_Mpc3 = (Default_N_nu_per100Mpc / 100.0) ** 3
# Read command line arguments
if len(sys.argv) <= 1 or sys.argv[1] == "--help" or sys.argv[1] == "-h":
@@ -29,13 +29,13 @@ print("")
# Check the unit system
if "Units" in f.keys() and "Unit length in cgs (U_L)" in f["Units"].attrs.keys():
- Length_Unit = f["Units"].attrs["Unit length in cgs (U_L)"] / Mpc_cgs # Mpc
+ Length_Unit = f["Units"].attrs["Unit length in cgs (U_L)"] / Mpc_cgs # Mpc
else:
- Length_Unit = 1.0 # Mpc
+ Length_Unit = 1.0 # Mpc
# Extract the box dimensions and volume
-L = f["Header"].attrs["BoxSize"] / Length_Unit # Mpc
-V = L**3 if np.isscalar(L) else np.product(L) # Mpc^3
+L = f["Header"].attrs["BoxSize"] / Length_Unit # Mpc
+V = L ** 3 if np.isscalar(L) else np.product(L) # Mpc^3
if not np.isscalar(L) and len(L) != 3:
raise ValueError("Box dimensions are not cubic")
@@ -47,20 +47,31 @@ if nparts[6] != 0 or "PartType6" in f.keys():
raise IOError("This file already has neutrinos.")
# Compute the default number of neutrinos (round to nearest cubic number)
-Default_N_nu = round((Default_nr_neutrinos_per_Mpc3 * V)**(1./3.))
-Default_Nr_neutrinos = int(Default_N_nu**3)
+Default_N_nu = round((Default_nr_neutrinos_per_Mpc3 * V) ** (1.0 / 3.0))
+Default_Nr_neutrinos = int(Default_N_nu ** 3)
print("The box dimensions are " + str(L) + " Mpc.")
-print("The default number of neutrinos is " +
- "%g" % Default_N_nu_per100Mpc + "^3 per (100 Mpc)^3.")
-print("The default number of neutrinos is " +
- "%g" % Default_N_nu + "^3 = " + str(Default_Nr_neutrinos) + ".")
+print(
+ "The default number of neutrinos is "
+ + "%g" % Default_N_nu_per100Mpc
+ + "^3 per (100 Mpc)^3."
+)
+print(
+ "The default number of neutrinos is "
+ + "%g" % Default_N_nu
+ + "^3 = "
+ + str(Default_Nr_neutrinos)
+ + "."
+)
print("")
-#Request the number of neutrino particles to be spawned (with default option)
-Nr_neutrinos = int(input("Enter the number of neutrinos (default " +
- "%d" % Default_Nr_neutrinos + "): ")
- or "%d" % Default_Nr_neutrinos)
+# Request the number of neutrino particles to be spawned (with default option)
+Nr_neutrinos = int(
+ input(
+ "Enter the number of neutrinos (default " + "%d" % Default_Nr_neutrinos + "): "
+ )
+ or "%d" % Default_Nr_neutrinos
+)
nparts[6] = Nr_neutrinos
@@ -80,7 +91,7 @@ print("The first particle ID of the first neutrino will be: " + str(firstID))
print("")
confirm = input("Enter y to confirm: ")
-if (confirm != "y"):
+if confirm != "y":
print("Not confirmed. Done for now.")
exit(0)