from __future__ import print_function import numpy as np import h5py def read_cite_data_pbh_main(): # Puchwein+2019 UVB model (No need to change this) uvb_model = "P19" # zmean at which tau data needs to be read zmean = 2.4 # Location of directory with files inp_path = "./" # Name of Temperature evolution data file inp_file_T0_evol = "CITE_Temperature_Evolution_Data.hdf5" # Name of tau files for a give zmean inp_file_tau = generate_tau_file_name(zmean) # T0_type Could be "T0_best_fit" or "T0_lb" or "T0_ub" # T0_best_fit -->> P19 photo-heating rates scales to go through BEST-FIT T0 measurements # T0_lb -->> P19 photo-heating rates scales to go through LOWER 1-sigma T0 measurements # T0_lb -->> P19 photo-heating rates scales to go through UPPER 1-sigma T0 measurements T0_type = "T0_best_fit" # Could be "T0_best_fit" or "T0_lb" or "T0_ub" # pbh_model could be any model from the following list # The list of PBH models for which CITE needs to be run. Note that the # model not listed in this list, will not be generated # pbh_model_list = ["m115","m215","m315","m415","m515","m615","m715",\ # "m815","m915","m116","m216","m316","m416","m516",\ # "m616","m716","m816","m916","m117"] pbh_model = "m415" # Fraction of Primordial Black-Hole Heating could be anything from following list # log_f_pbh_list = [-10.0,-9.0,-8.0,-7.0,-6.0,-5.0,-4.0,-3.0,-2.0,-1.0,0.0] log_f_pbh = -7.0 # Generate all the hdf5 blocks (No need to change this function) block_dict = generate_hdf5_block_dictionary(pbh_model,log_f_pbh,T0_type,uvb_model) # --------------------------------------------------------------- # Code to Read Optical Depth data starts # --------------------------------------------------------------- hdf_filename_tau = inp_path + inp_file_tau z_block = block_dict["z_block"] tau_HI_block = block_dict["tau_HI_block"] z_arr = hdf_read_data(hdf_filename_tau,z_block) tau_HI_arr = hdf_read_data(hdf_filename_tau,tau_HI_block) print() print("*************************************************************************************") print("Optical Depth data reading done ...") print("Redshift array size : ",z_arr.shape) print("HI Optical depth data size : ",tau_HI_arr.shape) print("*************************************************************************************") print() # --------------------------------------------------------------- # --------------------------------------------------------------- # Code to Read T0 Evolution data starts # --------------------------------------------------------------- # Structure of T0_evol_data_arr hdf_filename_T0_evol = inp_path + inp_file_T0_evol T0_evol_block = block_dict["T0_evol_block"] T0_evol_data_arr = hdf_read_data(hdf_filename_T0_evol,T0_evol_block) print() print("*************************************************************************************") print("CITE Temperature Evolution data reading done Array Size : ",T0_evol_data_arr.shape) print() print("Redshifts can be accessed by using T0_evol_data_arr[:,0] (see comments in code)") print("T0 Evolution can be accessed by using T0_evol_data_arr[:,7] (see comments in code)") print() print("To plot T0 evolution use following code ax.plot(T0_evol_data_arr[:,0],T0_evol_data_arr[:,7])") print("*************************************************************************************") print() # Structure of T0_evol_data_arr # T0_evol_data_arr[:,0] -->> zmean # T0_evol_data_arr[:,1] -->> Gamma_HI # T0_evol_data_arr[:,2] -->> Gamma_HeI # T0_evol_data_arr[:,3] -->> Gamma_HeII # T0_evol_data_arr[:,4] -->> E_ph_HI # T0_evol_data_arr[:,5] -->> E_ph_HeI # T0_evol_data_arr[:,6] -->> E_ph_HeII # T0_evol_data_arr[:,7] -->> T0 # T0_evol_data_arr[:,8] -->> x_HI_mean # T0_evol_data_arr[:,9] -->> x_HeII_mean # --------------------------------------------------------------- def generate_tau_file_name(zmean): inp_file_tau = "cite2048_n5000_z" + str("%.3f"%round(zmean,3)) + ".hdf5" return inp_file_tau def generate_hdf5_block_dictionary(pbh_model,log_f_pbh,T0_type="T0_best_fit",uvb_model="P19"): """ No need to change this function. This function generates the hdf5 blocks names where different data is stored """ # Scaling of UVB models used to match best fit / Lower and Upper T0 measurements # No need to change the following dictionary T0_dict = {\ "T0_best_fit":{"factor_a_HI":0.8,"factor_a_HeI":0.8,"factor_a_HeII":0.90},\ "T0_lb" :{"factor_a_HI":0.6,"factor_a_HeI":0.6,"factor_a_HeII":0.75},\ "T0_ub" :{"factor_a_HI":1.2,"factor_a_HeI":1.2,"factor_a_HeII":1.10},\ } # No need to change this parameter factor_b = 0.0 pbh_model_block = "PBH_Model-" + pbh_model log_f_pbh_block = "log_f_pbh-" + str("%.3f"%log_f_pbh) factor_a_HI = T0_dict[T0_type]["factor_a_HI"] factor_a_HeI = T0_dict[T0_type]["factor_a_HeI"] factor_a_HeII = T0_dict[T0_type]["factor_a_HeII"] factor_b_HI = factor_b factor_b_HeI = factor_b factor_b_HeII = factor_b factor_a_HI_block = "Factor_a_HI-" + str("%.3f"%factor_a_HI) factor_a_HeI_block = "Factor_a_HeI-" + str("%.3f"%factor_a_HeI) factor_a_HeII_block = "Factor_a_HeII-" + str("%.3f"%factor_a_HeII) factor_a_block = factor_a_HI_block + "_" + factor_a_HeI_block + \ "_" + factor_a_HeII_block factor_b_HI_block = "Factor_b_HI-" + str("%.3f"%factor_b_HI) factor_b_HeI_block = "Factor_b_HeI-" + str("%.3f"%factor_b_HeI) factor_b_HeII_block = "Factor_b_HeII-" + str("%.3f"%factor_b_HeII) factor_b_block = factor_b_HI_block + "_" + factor_b_HeI_block + \ "_" + factor_b_HeII_block common_block = uvb_model + "/" + factor_a_block + "/" + \ factor_b_block + "/" + pbh_model_block + "/" + \ log_f_pbh_block + "/" T0_evol_block = common_block + "UVB_Temperature_Evolution" tau_HI_block = common_block + "/Optical_Depth_HI" z_block = "Redshift" block_dict = vars() return block_dict def hdf_read_data(hdf_filename,blockname,flag_exit="y"): hdf_read = h5py.File(hdf_filename,"r") if hdf_read.__contains__(blockname): # Check if blockname already exists data = np.array(hdf_read[blockname]) hdf_read.close() else: print("Following Block name does not exists ...") print("File :",hdf_filename) print("Block :",blockname) hdf_read.close() if flag_exit == "y": exit() else: data = None return data read_cite_data_pbh_main()