diff --git a/example_notebooks/data/TNFW_example.hdf5 b/example_notebooks/data/TNFW_example.hdf5 new file mode 100644 index 00000000..1398c459 Binary files /dev/null and b/example_notebooks/data/TNFW_example.hdf5 differ diff --git a/pyHalo/Halos/HaloModels/TNFWFromParams.py b/pyHalo/Halos/HaloModels/TNFWFromParams.py new file mode 100644 index 00000000..c7bab198 --- /dev/null +++ b/pyHalo/Halos/HaloModels/TNFWFromParams.py @@ -0,0 +1,97 @@ +from pyHalo.Halos.halo_base import Halo +from lenstronomy.LensModel.Profiles.tnfw import TNFW as TNFWLenstronomy +import numpy as np +from pyHalo.Halos.tnfw_halo_util import tnfw_mass_fraction +from pyHalo.Halos.HaloModels.TNFW import TNFWSubhalo +from lenstronomy.LensModel.Profiles.tnfw import TNFW + +class TNFWFromParams(TNFWSubhalo): + """ + Creates a TNFW halo based on physical params. + """ + + KEY_RT = "r_trunc_kpc" + KEY_RS = "rs" + KEY_RHO_S = "rho_s" + KEY_RV = "rv" + + def __init__(self, mass, x_kpc, y_kpc, r3d, z, + sub_flag, lens_cosmo_instance, args, unique_tag=None): + """ + Defines a TNFW subhalo with physical params r_trunc_kpc, rs, rhos passed in the args argument + """ + + self._lens_cosmo = lens_cosmo_instance + + self._kpc_per_arcsec_at_z = self._lens_cosmo.cosmo.kpc_proper_per_asec(z) + + x = x_kpc / self._kpc_per_arcsec_at_z + + y = y_kpc / self._kpc_per_arcsec_at_z + + keys_physical = (self.KEY_RV,self.KEY_RS,self.KEY_RHO_S,self.KEY_RV,self.KEY_RT) + self._params_physical = {key:args[key] for key in keys_physical} + + self._c = self._params_physical[self.KEY_RV] / self._params_physical[self.KEY_RS] + + super(TNFWFromParams, self).__init__(mass,x,y,r3d,z,sub_flag,lens_cosmo_instance,args,None,None,unique_tag) + + def density_profile_3d(self, r, params_physical=None): + """ + Computes the 3-D density profile of the halo + :param r: distance from center of halo [kpc] + :return: the density profile in units M_sun / kpc^3 + """ + + _params = self._params_physical if params_physical is None else params_physical + + r_t = _params[self.KEY_RT] + r_s = _params[self.KEY_RS] + rho_s = _params[self.KEY_RHO_S] + + n = 1 + + x = r / r_s + tau = r_t / r_s + + density_nfw = (rho_s / ((x)*(1+x)**2)) + + return density_nfw * (tau**2 / (tau**2 + x**2))**n + + + @property + def profile_args(self): + """ + See documentation in base class (Halos/halo_base.py) + """ + if not hasattr(self, '_profile_args'): + truncation_radius_kpc = self.params_physical[self.KEY_RT] + self._profile_args = (self.c, truncation_radius_kpc) + return self._profile_args + + + @property + def lenstronomy_params(self): + """ + See documentation in base class (Halos/halo_base.py) + """ + if not hasattr(self, '_kwargs_lenstronomy'): + + r_t = self.params_physical[self.KEY_RT] + r_s = self.params_physical[self.KEY_RS] + rho_s = self.params_physical[self.KEY_RHO_S] + + Rs_angle, theta_Rs = self.nfw_physical2angle_fromNFWparams(rho_s,r_s,self.z) + + x, y = np.round(self.x, 4), np.round(self.y, 4) + + Rs_angle = np.round(Rs_angle, 10) + theta_Rs = np.round(theta_Rs, 10) + r_trunc_arcsec = r_t / self._lens_cosmo.cosmo.kpc_proper_per_asec(self.z) + + kwargs = [{'alpha_Rs': self._rescale_norm * theta_Rs, 'Rs': Rs_angle, + 'center_x': x, 'center_y': y, 'r_trunc': r_trunc_arcsec}] + + self._kwargs_lenstronomy = kwargs + + return self._kwargs_lenstronomy, None \ No newline at end of file diff --git a/pyHalo/Halos/galacticus_util/__init__.py b/pyHalo/Halos/galacticus_util/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/pyHalo/Halos/galacticus_util/galacticus_filter.py b/pyHalo/Halos/galacticus_util/galacticus_filter.py new file mode 100644 index 00000000..ff345ee0 --- /dev/null +++ b/pyHalo/Halos/galacticus_util/galacticus_filter.py @@ -0,0 +1,97 @@ +from pyHalo.Halos.galacticus_util.galacticus_util import GalacticusUtil +import numpy as np + +def nodedata_apply_filter(nodedata,filter,galacticus_util = None): + """Takes a dictionary with numpy arrays as values and apply as boolean filter to all. + Returns a dictionary with same keys, but a filter applied to all arrays.""" + + return {key:val[filter] for key,val in nodedata.items()} + +def nodedata_filter_tree(nodedata, treenum,galacticus_util = None): + """Returns a filter that excludes all nodes but nodes in the specified tree""" + gutil = GalacticusUtil() if galacticus_util is None else galacticus_util + return nodedata[gutil.PARAM_TREE_INDEX] == treenum + +def nodedata_filter_subhalos(nodedata,galacticus_util= None): + """Returns a filter that excludes all but subhalos (excludes host halos)""" + gutil = GalacticusUtil() if galacticus_util is None else galacticus_util + return nodedata[gutil.PARAM_ISOLATED] == 0 + +def nodedata_filter_halos(nodedata,galacticus_util = None): + """Returns a filter that excludes all but halo (excludes sub-halos)""" + return np.logical_not(nodedata_filter_subhalos(nodedata)) + +def nodedata_filter_range(nodedata,range,key,galacticus_util=None): + """Returns a filter that excludes nodes not within the given range for a specified parameter""" + return (nodedata[key] > range[0]) & (nodedata[key] < range[1]) + +def nodedata_filter_virialized(nodedata,galacticus_util= None): + """ + Returns a filter that excludes everything outside of the host halos virial radius + WARNING: Current implementation only works if there is only one host halo per tree, + IE we are looking at the last output from galacticus. + """ + gutil = GalacticusUtil() if galacticus_util is None else galacticus_util + + #Get radial position of halos + rvec = np.asarray((nodedata[gutil.PARAM_X],nodedata[gutil.PARAM_Y],nodedata[gutil.PARAM_Z])) + r = np.linalg.norm(rvec,axis=0) + + #Filter halos and get there virial radii + filtered_halos = nodedata_filter_halos(nodedata) + rv_halos = nodedata[gutil.PARAM_RADIUS_VIRIAL][filtered_halos] + halo_output_n = nodedata[gutil.PARAM_TREE_ORDER] + + filter_virialized = r < rv_halos[halo_output_n] + + return filter_virialized + +def nodedata_filter_r2d(nodedata,r2d_max,plane_normal, + galacticus_util= None): + + """ + Filters based on projected radii. + """ + gutil = GalacticusUtil() if galacticus_util is None else galacticus_util + + r = np.asarray((nodedata[gutil.PARAM_X],nodedata[gutil.PARAM_Y],nodedata[gutil.PARAM_Z])) + + r2d = project_r2d(*r,plane_normal) + + return r2d < r2d_max + +def project_r2d(x,y,z,plane_normal): + """ + Takes in arrays of coordinates, calculates the projected radius on the plane. + + + :param x: An array of x coordinates + :param y: An array of y coordinates + :param z: An array of z coordinates + :plane_normal: Normal vector of the plane to project onto + """ + coords = np.asarray((x,y,z)) + + #Reshape so if scalers are passed for x,y,z we do not encounter issue + if coords.ndim == 1: + coords.reshape((3,1)) + + #convert to unit normal + un = plane_normal / np.linalg.norm(plane_normal) + + #Project distance to plane + #Projected distance to plane is sqrt(r.r - (r.un)^2) + rdotr = np.linalg.norm(coords,axis=0)**2 + rdotun = np.dot(un,coords) + return np.sqrt(rdotr - rdotun**2) + + + + + + + + + + + diff --git a/pyHalo/Halos/galacticus_util/galacticus_util.py b/pyHalo/Halos/galacticus_util/galacticus_util.py new file mode 100644 index 00000000..d2bcd6a8 --- /dev/null +++ b/pyHalo/Halos/galacticus_util/galacticus_util.py @@ -0,0 +1,237 @@ +import h5py +import numpy as np +from typing import Iterable + +class GalacticusUtil(): + """ + This class provides utility for working with galacticus files. + """ + + #Names of galacticus parameters + #Position of the subhalo relative to the top level host halo + PARAM_X = "positionOrbitalX" + """The GALACTICUS output parameter for the X coordinate of the subhalo relative to the main halo""" + PARAM_Y = "positionOrbitalY" + """The GALACTICUS output parameter for the Y coordinate of the subhalo relative to the main halo""" + PARAM_Z = "positionOrbitalZ" + """The GALACTICUS output parameter for the Z coordinate of the subhalo relative to the main halo""" + + #Position of the subhalo relative to it's host halo or subhalo (not top level) + PARAM_X_REL = "satellitePositionX" + + PARAM_Y_REL = "satellitePositionY" + + PARAM_Z_REL = "satellitePositionZ" + + + PARAM_MASS_BOUND = "satelliteBoundMass" + + PARAM_MASS_INFALL = PARAM_MASS_BASIC = "basicMass" + """The infall mass of the subhalo""" + + PARAM_ISOLATED = "nodeIsIsolated" + + PARAM_HIERARCHY = "nodeHierarchyLevel" + + PARAM_RADIUS_VIRIAL = 'darkMatterOnlyRadiusVirial' + + PARAM_SPHERE_RADIUS = "spheroidRadius" + + PARAM_SPHERE_ANGULAR_MOMENTUM = "spheroidAngularMomentum" + + PARAM_SPHERE_MASS_STELLAR = "spheroidMassStellar" + + PARAM_SPHERE_MASS_GAS = "spheroidMassGas" + + PARAM_RADIUS_SCALE = "darkMatterProfileScaleRadius" + + PARAM_DENSITY_PROFILE_RADIUS = "densityProfileRadius" + + PARAM_DENSITY_PROFILE = "densityProfile" + + PARAM_Z_LAST_ISOLATED = "redshiftLastIsolated" + + PARAM_TNFW_RADIUS_TRUNCATION = "radiusTidalTruncationNFW" + PARAM_TNFW_RHO_S = "densityNormalizationTidalTruncationNFW" + + PARAM_TREE_ORDER = "custom_treeOutputOrder" + PARAM_TREE_INDEX = "custom_treeIndex" + + HDF5_GROUP_OUTPUT_PRIMARY = "Outputs" + """ + The name of the primary output group in a GALACTICUS output file. + This is the primary output group that contains nodedata at different output times + """ + + HDF5_GROUP_OUTPUT_N_PREFIX = "Output" + """The name of the prefix of the "OutputN" groups in GALACTICUS output file""" + + HDF5_GROUP_NODEDATA = "nodeData" + """The name of the hdf5 group containing nodedata.""" + + HDF5_DSET_TREECOUNT = "mergerTreeCount" + + HDF5_DSET_TREEINDEX = "mergerTreeIndex" + + HDF5_DSET_TREESTART = "mergerTreeStartIndex" + + def hdf5_read_output_indices(self,f): + """ + Returns the numbers assigned to the various galacticus outputs. + For if the Galacticus file has the following groups: Outputs/Output1, Outputs/Output10 + this function will return [1,10] + + :param f: h5py.File read from Galacticus output + """ + group_outputs = f[self.HDF5_GROUP_OUTPUT_PRIMARY] + + outputs = [] + + trim = len(self.HDF5_GROUP_OUTPUT_N_PREFIX) + for key in group_outputs.keys(): + try: + if(key[:trim]) != self.HDF5_GROUP_OUTPUT_N_PREFIX: + continue + outputs.append(int(key[trim:])) + except(ValueError): + pass + + return np.array(outputs) + + + def hdf5_access_output_n_nodedata(self,f,output_n): + """ + Returns the Outputs/OutputN/nodedata groups + + :param f: h5py.File read from Galacticus output + :param output_n: The number corresponding to the output to be read + """ + return f[self.HDF5_GROUP_OUTPUT_PRIMARY][f"{self.HDF5_GROUP_OUTPUT_N_PREFIX}{output_n}"][self.HDF5_GROUP_NODEDATA] + + + def hdf5_read_dset(self,dset): + """ + Reads a hdf5 dataset into a numpy array + """ + arr = np.zeros(dset.shape,dtype=dset.dtype) + dset.read_direct(arr) + return arr + + def hdf5_read_custom_nodedata(self,f,output_index): + """ + To make analysis more convenient, it is useful to add custom datasets to those present in the Outputs/OutputN/nodedata. + This function is used when read_nodedata_galacticus is called to create custom nodedata datasets for convenience. + + :param f: h5py.File read from Galacticus output + :param output_index: The output to read indices for. + """ + group_outputn = f[self.HDF5_GROUP_OUTPUT_PRIMARY][f"{self.HDF5_GROUP_OUTPUT_N_PREFIX}{output_index}"] + + tree_index = self.hdf5_read_dset(group_outputn[self.HDF5_DSET_TREEINDEX]) + tree_start = self.hdf5_read_dset(group_outputn[self.HDF5_DSET_TREESTART]) + + total_count = self.hdf5_read_nodecount_total(f,output_index) + + node_index = np.zeros(total_count,dtype=int) + node_order = np.zeros(total_count,dtype=int) + + for n in range(1,len(tree_start)): + start,stop = tree_start[n-1],tree_start[n] + node_order[start:stop] = n - 1 + node_index[start:stop] = tree_index[n - 1] + + node_order[stop:] = n + node_index[stop:] = tree_index[n] + + return {GalacticusUtil.PARAM_TREE_INDEX:node_index, + GalacticusUtil.PARAM_TREE_ORDER:node_order} + + def hdf5_read_nodecount_total(self,f,output_index): + """ + Returns the total number of nodes at a given output. + + :param f: h5py.File read from Galacticus output + :param output_index: The output index to read from. + """ + + group_outputn = f[self.HDF5_GROUP_OUTPUT_PRIMARY][f"{self.HDF5_GROUP_OUTPUT_N_PREFIX}{output_index}"] + + tree_count = self.hdf5_read_dset(group_outputn[self.HDF5_DSET_TREECOUNT]) + + return np.sum(tree_count) + + def read_nodedata_galacticus(self,path, output_index = None, + params_to_read = None, + nodes_only=True): + """ + Reads a galacticus output file at a given path, returns nodedata for a specified tree. + Galacticus: https://github.com/galacticusorg/galacticus + + Note: if your build of galacticus has a different output formats you can inherit the GalacticusUtil class and modify + parameter names / methods used when reading galacticus output. + + :param path: Path to nodedata + :param output_index: The index of the tree to read. If None defaults to the final output. + :param params_to_read: A list of parameters to read, if None reads all parameters. + :nodes_only: If true only reads nodedata entries that have entry for each node, if false reads all datasets in the nodedata group. + """ + with h5py.File(path,"r") as f: + return self.hdf5_read_galacticus_nodedata(f,output_index,params_to_read,nodes_only) + + + def hdf5_read_galacticus_nodedata(self,f, output_index = None, + params_to_read = None, + nodes_only=True): + """ + Reads a galacticus output given a h5py.File read from galacticus output + Galacticus: https://github.com/galacticusorg/galacticus + + Note: if your build of galacticus has a different output formats you can inherit the GalacticusUtil class and modify + parameter names / methods used when reading galacticus output. + + :param f: A h5py.File read from galacticus output. + :param output_index: The index of the tree to read. If None defaults to the final output. + :param params_to_read: A list of parameters to read, if None reads all parameters. + :nodes_only: If true only reads nodedata entries that have entry for each node, if false reads all datasets in the nodedata group. + """ + + #If no output is specified read the final output + if output_index is None: + + outputs_ns = self.hdf5_read_output_indices(f) + + #If no suitable outputs are found, return None + + if len(outputs_ns) == 0: + + return None + + output_index = np.max(outputs_ns) + + group_nodedata = self.hdf5_access_output_n_nodedata(f,output_index) + + nodedata = {} + + #Allow user to pass a single parameter to read as just a string + + if isinstance(params_to_read,str): + + params_to_read = [params_to_read] + + nodecount = self.hdf5_read_nodecount_total(f,output_index) + + #Loop through datasets in nodedata + + for key in group_nodedata.keys(): + + dset = group_nodedata[key] + + #Read dataset if it is a dataset and it is contained in the params_to_read variable or that variable is None + if isinstance(dset,h5py.Dataset) and (params_to_read is None or key in params_to_read): + if not nodes_only or (nodes_only and dset.shape[0] == nodecount): + nodedata[key] = self.hdf5_read_dset(dset) + + + custom = self.hdf5_read_custom_nodedata(f,output_index) + return nodedata | custom + diff --git a/pyHalo/PresetModels/external.py b/pyHalo/PresetModels/external.py index 0ea4837e..f4f9db72 100644 --- a/pyHalo/PresetModels/external.py +++ b/pyHalo/PresetModels/external.py @@ -1,19 +1,164 @@ +import numpy as np +from scipy.spatial.transform import Rotation from pyHalo.Halos.HaloModels.TNFWemulator import TNFWSubhaloEmulator +from pyHalo.Halos.HaloModels.TNFWFromParams import TNFWFromParams from pyHalo.PresetModels.cdm import CDM from pyHalo.single_realization import Realization +from pyHalo.Halos.galacticus_util.galacticus_util import GalacticusUtil +from pyHalo.Halos.galacticus_util.galacticus_filter import nodedata_filter_subhalos,nodedata_filter_tree,nodedata_filter_virialized,nodedata_filter_range,nodedata_apply_filter +from lenstronomy.LensModel.Profiles.tnfw import TNFW +import h5py -def galacticus_subhalos(zlens, zsource, galacticus_hdf5, mdef='TNFW', rmax_arcsec=3.0): - +def DMFromGalacticus(z_lens,z_source,galacticus_hdf5,tree_index, kwargs_cdm,mass_range,mass_range_is_bound = True, + proj_plane_normal = None,include_field_halos=True,nodedata_filter = None, + galacticus_utilities = None,galacticus_params_additional = None, proj_rotation_angles = None): """ - This module loads a realization of halos from galacticus and transforms it into a pyHalo-style Realization - :param zlens: main deflector redshift - :param zsource: source redshift - :param galacticus_hdf5: an HDF5 file the contains galacticus output - :param rmax_arcsec: the maximum radius out to which one wants to render subhalos [arcsec] - :return: + This generates a realization of halos using subhalo parameters provided from a specified tree in the galacticus file. + See https://github.com/galacticusorg/galacticus/ for information on the galacticus galaxy formation model. + + :param z_lens: main deflector redshift + :param z_source: source redshift + :param galacticus_hdf5: Galacticus output hdf5 file. If str loads file from specified path + :param tree_index: The number of the tree to create a realization from. A single galacticus file contains multiple realizations (trees). + NOTE: Trees output from galacticus are assigned a number starting at 1. + :param mass_range: Specifies mass range to include subhalos within. + :param mass_range_is_bound: If true subhalos are filtered bound mass, if false subhalos are filtered by infall mass. + :param projection_normal: Projects the coordinates of subhalos from parameters onto a plane defined with the given (3D) normal vector. + Use this to generate multiple realizations from a single galacticus tree. If is None, coordinates are projected on the x,y plane. + :param include_field_halos: If true includes field halos, if false no field halos are included. + :param nodedata_filter: Expects a callable function that has input and output: (dict[str,np.ndarray], GalacticusUtil) -> np.ndarray[bool] + ,subhalos are filtered based on the output np array. Defaults to None + :param galacticus_params: Extra parameters to read when loading in galacticus hdf5 file. + :param proj_rotation_angle: Alternative to providing proj_plane_normal argument. Expects a length 2 array: (theta, phi) + with theta and phi being the angles in spherical coordinates of the normal vector of defining the plane to project coordinates onto. + + :return: A realization from Galacticus halos """ - pass + gutil = GalacticusUtil() if galacticus_utilities is None else galacticus_utilities + + MPC_TO_KPC = 1E3 + + #Only read needed parameters to save memory. + PARAMS_TO_READ_DEF = [gutil.PARAM_X,gutil.PARAM_Y,gutil.PARAM_Z,gutil.PARAM_TNFW_RHO_S, + gutil.PARAM_TNFW_RADIUS_TRUNCATION,gutil.PARAM_RADIUS_VIRIAL, + gutil.PARAM_RADIUS_SCALE,gutil.PARAM_MASS_BOUND, + gutil.PARAM_MASS_BASIC,gutil.PARAM_ISOLATED] + + if galacticus_params_additional is None: + galacticus_params_additional = [] + else: + galacticus_params_additional = list(galacticus_params_additional) + + params_to_read = galacticus_params_additional + PARAMS_TO_READ_DEF + + # we create a realization of only line-of-sight halos by setting sigma_sub = 0.0 + # only include these halos if requested + kwargs_cdm['sigma_sub'] = 0.0 + + los_norm = 1 if include_field_halos else 0 + los_norm = kwargs_cdm.get("LOS_normalization") if not kwargs_cdm.get("LOS_normalization") is None else los_norm + + + cdm_halos_LOS = CDM(z_lens, z_source, **kwargs_cdm,LOS_normalization=los_norm) + + # get lens_cosmo class from class containing LOS objects; note that this will work even if there are no LOS halos + lens_cosmo = cdm_halos_LOS.lens_cosmo + + if isinstance(galacticus_hdf5,str): + nodedata = gutil.read_nodedata_galacticus(galacticus_hdf5,params_to_read=params_to_read) + else: + nodedata = gutil.hdf5_read_galacticus_nodedata(galacticus_hdf5,params_to_read=params_to_read) + + + #Set up for rotation of coordinates + #Specify the normal vector for the plane we are projecting onto, if user specified ensure the vector is normalized + nh = np.asarray((0,0,1)) if proj_plane_normal is None else proj_plane_normal / np.linalg.norm(proj_plane_normal) + nh_x,nh_y,nh_z = nh + + theta = np.arccos(nh_z) + phi = np.sign(nh_y) * np.arccos(nh_x/np.sqrt(nh_x**2 + nh_y**2)) if nh_x != 0 or nh_y != 0 else 0 + + if not proj_rotation_angles is None: + theta,phi = proj_rotation_angles + + + #This rotation rotation maps the coordinates such that in the new coordinates zh = nh and the x,y coordinates after rotation + #are the x-y coordinates in the plane + rotation = Rotation.from_euler("zyz",(0,theta,phi)) + + coords = np.asarray((nodedata[gutil.PARAM_X],nodedata[gutil.PARAM_Y],nodedata[gutil.PARAM_Z])) * MPC_TO_KPC + + #We would like define the x and y unit vectors, so we can project our coordinates + xh_r = rotation.apply(np.array((1,0,0))) + yh_r = rotation.apply(np.array((0,1,0))) + + kpc_per_arcsec_at_z = lens_cosmo.cosmo.kpc_proper_per_asec(z_lens) + + #Get the maximum r2d for the subhalo to be within the rendering volume + r2dmax_kpc = (kwargs_cdm["cone_opening_angle_arcsec"] / 2) * kpc_per_arcsec_at_z + + coords_2d = np.asarray((np.dot(xh_r,coords),np.dot(yh_r,coords))) + r2d_mag = np.linalg.norm(coords_2d,axis=0) + + filter_r2d = r2d_mag < r2dmax_kpc + + #Choose to filter by bound / infall mass + mass_key = gutil.PARAM_MASS_BOUND if mass_range_is_bound else gutil.PARAM_MASS_BASIC + + # Filter subhalos + # We should exclude nodes that are not valid subhalos, such as host halo nodes and nodes that are outside the virial radius. + # (Galacticus is not calibrated for nodes outside of virial radius) + + if nodedata_filter is None: + filter_subhalos = nodedata_filter_subhalos(nodedata,gutil) + filter_virialized = nodedata_filter_virialized(nodedata,gutil) + filter_mass = nodedata_filter_range(nodedata,mass_range,mass_key,gutil) + filter_tree = nodedata_filter_tree(nodedata,tree_index,gutil) + + filter_combined = filter_subhalos & filter_virialized & filter_mass & filter_tree & filter_r2d + + else: + filter_combined = nodedata_filter(nodedata,gutil) + + #Apply filter to nodedata and rvec + nodedata = nodedata_apply_filter(nodedata,filter_combined) + coords_2d = coords_2d[:,filter_combined] + r2d_mag = r2d_mag[filter_combined] + coords = coords[:,filter_combined] + r3d_mag = np.linalg.norm(coords,axis=0) + + # Get rhos_s factor of 4 comes from the this galacticus output is + # The density normalization of the underlying NFW halo at r = rs + # Multiply by 4 to get the normalization for the halo profile + rho_s = 4 * nodedata[gutil.PARAM_TNFW_RHO_S] / (MPC_TO_KPC)**3 + + + rs = nodedata[gutil.PARAM_RADIUS_SCALE] * MPC_TO_KPC + rt = nodedata[gutil.PARAM_TNFW_RADIUS_TRUNCATION] * MPC_TO_KPC + rv = nodedata[gutil.PARAM_RADIUS_VIRIAL] * MPC_TO_KPC + + halo_list = [] + #Loop thought properties of each subhalos + for n,m_infall in enumerate(nodedata[gutil.PARAM_MASS_BASIC]): + x,y = coords_2d[0][n], coords_2d[1][n] + + tnfw_args = { + TNFWFromParams.KEY_RT:rt[n], + TNFWFromParams.KEY_RS:rs[n], + TNFWFromParams.KEY_RHO_S:rho_s[n], + TNFWFromParams.KEY_RV:rv[n] + } + + + + halo_list.append(TNFWFromParams(m_infall,x,y,r3d_mag[n],z_lens,True,lens_cosmo,tnfw_args)) + + subhalos_from_params = Realization.from_halos(halo_list,lens_cosmo,kwargs_halo_model={}, + msheet_correction=False, rendering_classes=None) + + return cdm_halos_LOS.join(subhalos_from_params) + diff --git a/pyHalo/preset_models.py b/pyHalo/preset_models.py index c3263f79..6b97bd33 100644 --- a/pyHalo/preset_models.py +++ b/pyHalo/preset_models.py @@ -5,11 +5,12 @@ presented here show what each keyword argument accepted by pyHalo does. """ from pyHalo.PresetModels.cdm import CDM -from pyHalo.PresetModels.external import CDMFromEmulator +from pyHalo.PresetModels.external import CDMFromEmulator, DMFromGalacticus from pyHalo.PresetModels.sidm import SIDM_core_collapse from pyHalo.PresetModels.uldm import ULDM from pyHalo.PresetModels.wdm import WDM, WDM_mixed, WDMGeneral + __all__ = ['preset_model_from_name'] def preset_model_from_name(name): @@ -32,7 +33,26 @@ def preset_model_from_name(name): return WDM_mixed elif name == 'WDMGeneral': return WDMGeneral + elif name == "DMGalacticus": + return DMFromGalacticus else: raise Exception('preset model '+ str(name)+' not recognized!') + + + + + + + + + + + + + + + + + diff --git a/requirements.txt b/requirements.txt index 56d5bae1..22ffddaf 100644 --- a/requirements.txt +++ b/requirements.txt @@ -11,6 +11,7 @@ numpy>=1.13 astropy>=2.0 scipy>=0.19.1 colossus>=1.3.1 +h5py>=3.9.0 mcfit>=0.0.21 lenstronomy @ git+https://github.com/lenstronomy/lenstronomy mpmath