Skip to content
Snippets Groups Projects
Commit 9414bba1 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Applied python code formatting script

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