add mbh seeds

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):

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

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!')
 

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):
         """

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)

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)