diff --git a/pyHalo/Halos/HaloModels/blackhole.py b/pyHalo/Halos/HaloModels/blackhole.py index c2dd8e3..a35dc26 100644 --- a/pyHalo/Halos/HaloModels/blackhole.py +++ b/pyHalo/Halos/HaloModels/blackhole.py @@ -6,7 +6,9 @@ class BlackHole(Halo): Class that defines a point mass object in the lens model """ def __init__(self, mass, x, y, r3d, z, - sub_flag, lens_cosmo_instance, args, truncation_class, concentration_class, unique_tag): + sub_flag, lens_cosmo_instance, args, + truncation_class, concentration_class, unique_tag, + fixed_position=True): """ See documentation in base class (Halos/halo_base.py) """ @@ -16,7 +18,7 @@ def __init__(self, mass, x, y, r3d, z, mdef = 'PT_MASS' super(BlackHole, self).__init__(mass, x, y, r3d, mdef, z, sub_flag, lens_cosmo_instance, args, unique_tag, - fixed_position=True) + fixed_position=fixed_position) @property def lenstronomy_ID(self): diff --git a/pyHalo/PresetModels/mbh.py b/pyHalo/PresetModels/mbh.py new file mode 100644 index 0000000..f6d4a99 --- /dev/null +++ b/pyHalo/PresetModels/mbh.py @@ -0,0 +1,60 @@ +from pyHalo.PresetModels.cdm import CDM +from pyHalo.realization_extensions import RealizationExtensions + + +def CDM_plus_BH(z_lens, + z_source, + log10_mass_ratio, + log10_occupation_frac, + log10_mlow_halos_subres=5.0, + log10_min_mbh=4.5, + log10_mass_maximum=6.7, + sigma_sub=0.025, + log_mlow=6., + log_mhigh=10., + log10_sigma_sub=None, + shmf_log_slope=-1.9, + cone_opening_angle_arcsec=6., + log_m_host=13.3, + LOS_normalization=1.0, + geometry_type='DOUBLE_CONE', + add_globular_clusters=False, + kwargs_globular_clusters=None): + """ + Add a population of black holes to a CDM realization + :param z_lens: lens redshift + :param z_source: source redshift + :param log10_mass_ratio: the ratio of the black hole mass to the host halo mass + :param log10_occupation_frac: the fraction of halos with a bh in them + :param log10_mlow_halos_subres: the minimum halo mass that host black holes, can be smaller than the minimum + halo mass rendered in the model (set by log_mlow) + :param log10_min_mbh: the minimum black hole mass to render + :param log10_mass_maximum: the maximum black hole mass to render + :param sigma_sub: SHMF normalization + :param log_mlow: log10 minimum DM halo mass to render + :param log_mhigh: log10 maximum DM halo mass to render + :param log10_sigma_sub: SHMF normalization + :param shmf_log_slope: SHMF log-slope + :param cone_opening_angle_arcsec: the opening angle of the lensing volume + :param log_m_host: host halo mass + :param LOS_normalization: amplitude of the LOS HMF relative to Sheth Tormen + :param geometry_type: specifies the geometry of the rendering volume + :param add_globular_clusters: bool; include a population of globular clusters + :param kwargs_globular_clusters: keyword args for the GC population, see documentation in RealizationExtensions + :return: a realization with CDM halos plus black hole seeds + """ + + cdm = CDM(z_lens, z_source, sigma_sub, log_mlow, log_mhigh, + log10_sigma_sub, shmf_log_slope=shmf_log_slope, cone_opening_angle_arcsec=cone_opening_angle_arcsec, + log_m_host=log_m_host, LOS_normalization=LOS_normalization, geometry_type=geometry_type, + add_globular_clusters=add_globular_clusters, kwargs_globular_clusters=kwargs_globular_clusters) + ext = RealizationExtensions(cdm) + mbh = ext.add_black_holes(log10_mass_ratio, + 10**log10_occupation_frac, + log10_mlow_halos_subres, + log10_min_mbh, + log_mlow, + log10_mass_maximum, + LOS_normalization) + realization = cdm.join(mbh) + return realization diff --git a/pyHalo/preset_models.py b/pyHalo/preset_models.py index 6f5739d..02b60ee 100644 --- a/pyHalo/preset_models.py +++ b/pyHalo/preset_models.py @@ -40,6 +40,9 @@ def preset_model_from_name(name): elif name == "DMGalacticus": from pyHalo.PresetModels.external import DMFromGalacticus return DMFromGalacticus + elif name == 'CDM_plus_BH': + from pyHalo.PresetModels.mbh import CDM_plus_BH + return CDM_plus_BH else: raise Exception('preset model '+ str(name)+' not recognized!') diff --git a/pyHalo/realization_extensions.py b/pyHalo/realization_extensions.py index d5c19ec..821cda8 100644 --- a/pyHalo/realization_extensions.py +++ b/pyHalo/realization_extensions.py @@ -3,13 +3,14 @@ from pyHalo.Halos.HaloModels.generalized_nfw import GeneralNFWSubhalo, GeneralNFWFieldHalo from pyHalo.single_realization import Realization from pyHalo.Halos.HaloModels.gaussianhalo import GaussianHalo -from pyHalo.concentration_models import preset_concentration_models +from pyHalo.Halos.HaloModels.blackhole import BlackHole +from pyHalo.Rendering.MassFunctions.density_peaks import ShethTormen from pyHalo.Rendering.correlated_structure import CorrelatedStructure from pyHalo.Rendering.MassFunctions.delta_function import DeltaFunction from pyHalo.Rendering.MassFunctions.gaussian import Gaussian from pyHalo.Cosmology.geometry import Geometry from pyHalo.utilities import generate_lens_plane_redshifts, mask_annular -from pyHalo.Rendering.SpatialDistributions.uniform import Uniform +from pyHalo.Rendering.SpatialDistributions.uniform import Uniform, LensConeUniform from copy import deepcopy from scipy.interpolate import RectBivariateSpline import time @@ -34,6 +35,106 @@ def __init__(self, realization): self._realization = realization + def add_black_holes(self, log10_mass_ratio, + f, + log10_mlow_halos_subres=5.0, + log10_min_mbh=4.5, + log_mlow_halos=6.0, + log10_mass_maximum=6.7, + LOS_normalization=1.0): + """ + Add a population of black holes in the center of halos + :param log10_mass_ratio: the ratio of the black hole to the mass of the host halo + :param f: the fraction of halos with a bh seed + :param log10_mlow_halos_subres: the minimum halo mass in which to inject BH seeds; this should be lower than the + minimum halo mass used to create the realization (log_mlow_halos) + :param log10_min_mbh: the minimum bh mass rendered + :param log_mlow_halos: the minimum halo mass of explicitely rendered halos in the realization + :param log10_mass_maximum: the maximum mass of the BH seeds + :param LOS_normalization: the overal normalization of the LOS halo mass function relative to Sheth Tormen + :return: a population of black holes + """ + # first we inject seeds into rendered halos (down to log10_minimum_halo_mass) + black_hole_list = [] + for halo in self._realization.halos: + u = np.random.rand() + if u > f: + # no BH in this halo + continue + kpc_per_arcsec = self._realization.lens_cosmo.cosmo.kpc_proper_per_asec(halo.z) + x_center_halo, y_center_halo = halo.x, halo.y + m_bh = min(10**log10_mass_maximum, halo.mass * 10 ** log10_mass_ratio) + halo_scale_radius_arcsec = halo.nfw_params[1] / kpc_per_arcsec + theta = np.random.uniform(0, 2*np.pi) + costheta, sintheta = np.cos(theta), np.sin(theta) + R = np.sqrt(np.random.uniform(0, halo_scale_radius_arcsec ** 2)) + x_bh, y_bh = x_center_halo + R * costheta, y_center_halo + R * sintheta + if m_bh > 10**log10_min_mbh: + mbh = BlackHole(m_bh, + x_bh, + y_bh, + r3d=None, + z=halo.z, + sub_flag=False, + lens_cosmo_instance=halo.lens_cosmo, + args={}, + truncation_class=None, + concentration_class=None, + unique_tag=np.random.rand(), + fixed_position=False) + black_hole_list.append(mbh) + + plane_redshifts = self._realization.unique_redshifts + delta_z = [] + for i, zi in enumerate(plane_redshifts[0:-1]): + delta_z.append(plane_redshifts[i + 1] - plane_redshifts[i]) + delta_z.append(self._realization.lens_cosmo.z_source - plane_redshifts[-1]) + + for (zi, delta_zi) in zip(plane_redshifts, delta_z): + kwargs_model_subres = {'m_pivot': 10 ** 8, + 'log_mlow': log10_mlow_halos_subres, + 'log_mhigh': log_mlow_halos, + 'LOS_normalization': f * LOS_normalization, + 'delta_power_law_index': 0.0, + 'draw_poisson': True} + mfunc_sub_resolution = ShethTormen.from_redshift(zi, delta_zi, + self._realization.geometry, + kwargs_model_subres) + mass_sub_resolution = mfunc_sub_resolution.draw() + uniform_spatial_distribution = LensConeUniform(self._realization.geometry.cone_opening_angle, + self._realization.geometry) + x_kpc, y_kpc = uniform_spatial_distribution.draw(len(mass_sub_resolution), + zi) + kpc_per_arcsec = self._realization.lens_cosmo.cosmo.kpc_proper_per_asec(zi) + x = x_kpc / kpc_per_arcsec + y = y_kpc / kpc_per_arcsec + for m_bh, x_bh, y_bh in zip(mass_sub_resolution, x, y): + if m_bh > 10 ** log10_min_mbh: + mbh = BlackHole(m_bh * 10 ** log10_mass_ratio, + x_bh, + y_bh, + r3d=None, + z=zi, + sub_flag=False, + lens_cosmo_instance=self._realization.lens_cosmo, + args={}, + truncation_class=None, + concentration_class=None, + unique_tag=np.random.rand(), + fixed_position=False) + black_hole_list.append(mbh) + + mbh_realization = Realization.from_halos(black_hole_list, self._realization.lens_cosmo, + self._realization.kwargs_halo_model, + self._realization.apply_mass_sheet_correction, + self._realization.rendering_classes, + self._realization._rendering_center_x, + self._realization._rendering_center_y, + self._realization.geometry) + return mbh_realization + # now we add mbh seeds into the "background" population of DM halos we didn't explicitely model + + def add_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radius_arcsec, gc_profile_args, galaxy_Re=6.0, host_halo_mass=10**13.3, f=3.4e-5, center_x=0, center_y=0): """ diff --git a/tests/test_preset_models.py b/tests/test_preset_models.py index be82248..898890f 100644 --- a/tests/test_preset_models.py +++ b/tests/test_preset_models.py @@ -1,5 +1,6 @@ from pyHalo.PresetModels.cdm import CDM, CDMCorrelatedStructure from pyHalo.PresetModels.wdm import WDM, WDM_mixed +from pyHalo.PresetModels.mbh import CDM_plus_BH from pyHalo.PresetModels.sidm import SIDM_core_collapse, SIDM_parametric from pyHalo.PresetModels.uldm import ULDM from pyHalo.preset_models import preset_model_from_name @@ -121,6 +122,17 @@ def test_WDM_mixed(self): _ = preset_model_from_name('WDM_mixed') self._test_default_infall_model(wdm_mixed, 'hybrid') + def test_CDM_blackholes(self): + + model = preset_model_from_name('CDM_plus_BH') + cdm_bh = model(0.5, + 1.5, + -0.2, + -0.3, + sigma_sub=0.01) + _ = cdm_bh.lensing_quantities() + _ = preset_model_from_name('CDM_plus_BH') + def test_WDM_general(self): func = preset_model_from_name('WDMGeneral') wdm = func(0.5, 1.5, 7.7, -2.0) diff --git a/tests/test_realization_extensions.py b/tests/test_realization_extensions.py index 95177e8..6debcf4 100644 --- a/tests/test_realization_extensions.py +++ b/tests/test_realization_extensions.py @@ -13,6 +13,19 @@ class TestRealizationExtensions(object): + def test_black_holes(self): + + cdm = CDM(0.5, 2.0, sigma_sub=0.0, LOS_normalization=1.0) + ext = RealizationExtensions(cdm) + log10_m_ratio = 0.0 + log10_f = np.log10(0.5) + mbh = ext.add_black_holes(log10_m_ratio, + log10_f, + log10_mlow_halos_subres=6.0) + for bh, nfw in zip(mbh.halos, cdm.halos): + npt.assert_almost_equal(bh.mass/2, nfw.mass) + npt.assert_almost_equal(bh.redshift, nfw.redshift) + def test_globular_clusters(self): cdm = CDM(0.5, 2.0, sigma_sub=0.01, LOS_normalization=0.1)