update globular cluster model

pyHalo/Halos/HaloModels/powerlaw.py

@@ -130,11 +130,10 @@ def __init__(self, mass, x, y, z, lens_cosmo_instance, args, unique_tag):
         :param y: y coordinate [arcsec]
         :param z: redshift
         :param lens_cosmo_instance: an instance of LensCosmo class
-        :param args: keyword arguments for the profile, must contain 'gamma', 'r_core_fraction', 'gc_size_lightyear'
-        which are the logarithmic profile slope outside the core, the core size in units of the GC size, and the gc size
-        in light years. The size and mass are related by
-        gc_size = gc_size_lightyear (m / 10^5)^1/3
-        The mass is the total mass inside a sphere with radius gc_size
+        :param args: keyword arguments for the profile, must contain 'gamma', 'gc_size_lightyear', 'gc_concentration'
+        which are the logarithmic profile slope outside the core, the size of the GC in lightyears
+        (the mass is defined as the mass inside a sphere with this radius), and the 'concentration' which is the ratio
+         of the gc core size to total size
         :param unique_tag:
         """
         self._prof = SPLCORE()
@@ -151,16 +150,14 @@ def lenstronomy_params(self):
         See documentation in base class (Halos/halo_base.py)
         """
         if not hasattr(self, '_lenstronomy_args'):
+            kpc_per_lightyear = 0.31 * 1e-3
             kpc_per_arcsec = self._lens_cosmo.cosmo.kpc_proper_per_asec(self.z)
             gamma = self._args['gamma']
-            r_core_fraction = self._args['r_core_fraction']
-            gc_size_lightyear_0 = self._args['gc_size_lightyear']
+            gc_size_lightyear = self._args['gc_size_lightyear']
+            gc_size_arcsec = gc_size_lightyear * kpc_per_lightyear / kpc_per_arcsec
+            r_core_arcsec = gc_size_arcsec/self._args['gc_concentration']
             sigma_crit_mpc = self._lens_cosmo.get_sigma_crit_lensing(self.z, self._lens_cosmo.z_source)
             sigma_crit_arcsec = sigma_crit_mpc * (0.001 * kpc_per_arcsec) ** 2
-            kpc_per_lightyear = 0.3 * 1e-3
-            gc_size_kpc = gc_size_lightyear_0 * (self.mass / 10 ** 5) ** (1 / 3) * kpc_per_lightyear
-            gc_size_arcsec = gc_size_kpc / kpc_per_arcsec
-            r_core_arcsec = r_core_fraction * gc_size_arcsec
             rho0 = self.mass / self._prof.mass_3d(gc_size_arcsec, sigma_crit_arcsec, r_core_arcsec, gamma)
             sigma0 = rho0 * r_core_arcsec
             self._lenstronomy_args = [{'sigma0': sigma0, 'gamma': gamma, 'center_x': self.x, 'center_y': self.y,
@@ -174,15 +171,14 @@ def profile_args(self):
         """
         if not hasattr(self, '_profile_args'):
 
+            kpc_per_lightyear = 0.31 * 1e-3
             gamma = self._args['gamma']
-            r_core_fraction = self._args['r_core_fraction']
-            gc_size_lightyear_0 = self._args['gc_size_lightyear']
-            kpc_per_lightyear = 0.3 * 1e-3
-            gc_size_kpc = gc_size_lightyear_0 * (self.mass / 10 ** 5) ** (1 / 3) * kpc_per_lightyear
-            r_core_kpc = r_core_fraction * gc_size_kpc
+            gc_size_lightyear = self._args['gc_size_lightyear']
+            c = self._args['gc_concentration']
+            gc_size_kpc = gc_size_lightyear * kpc_per_lightyear
+            r_core_kpc = gc_size_kpc / c
             rho0 = self.mass / self._prof.mass_3d(gc_size_kpc, 1.0, r_core_kpc, gamma)
-            self._profile_args = {'rho0': rho0, 'gc_size': gc_size_kpc,
-                                  'gamma': gamma, 'r_core_kpc': r_core_kpc}
+            self._profile_args = (rho0, gc_size_kpc, gamma, r_core_kpc)
         return self._profile_args
 
     @property

pyHalo/PresetModels/cdm.py

@@ -168,7 +168,6 @@ class for details
     if add_globular_clusters:
         from pyHalo.realization_extensions import RealizationExtensions
         ext = RealizationExtensions(realization)
-        kwargs_globular_clusters['host_halo_mass'] = 10 ** log_m_host
         realization = ext.add_globular_clusters(**kwargs_globular_clusters)
     return realization
 

pyHalo/PresetModels/wdm.py

@@ -177,7 +177,6 @@ def WDM(z_lens, z_source, log_mc, sigma_sub=0.025, log_mlow=6., log_mhigh=10., l
     if add_globular_clusters:
         from pyHalo.realization_extensions import RealizationExtensions
         ext = RealizationExtensions(realization)
-        kwargs_globular_clusters['host_halo_mass'] = 10 ** log_m_host
         realization = ext.add_globular_clusters(**kwargs_globular_clusters)
     return realization
 
@@ -404,6 +403,5 @@ def WDMGeneral(z_lens, z_source, log_mc, dlogT_dlogk, sigma_sub=0.025, log_mlow=
     if add_globular_clusters:
         from pyHalo.realization_extensions import RealizationExtensions
         ext = RealizationExtensions(realization)
-        kwargs_globular_clusters['host_halo_mass'] = 10 ** log_m_host
         realization = ext.add_globular_clusters(**kwargs_globular_clusters)
     return realization

pyHalo/realization_extensions.py

@@ -134,43 +134,58 @@ def add_black_holes(self, log10_mass_ratio,
         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):
+    def add_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radius_arcsec, gamma_mean=2.2,
+                              gamma_sigma=0.2, gc_concentration_mean=80, gc_concentration_sigma=20,
+                              gc_size_mean=100, gc_size_sigma=10, gc_surface_mass_density=10 ** 5,
+                              center_x=0, center_y=0):
         """
-        Add globular clusters at z = z_lens to a realization sampling from a log-normal mass distribution
+        Add globular clusters at main deflector redshift following a log-normal mass distribution
         :param log10_mgc_mean: the median log10(mass) of the GC's
         :param log10_mgc_sigma: the standard deviation of the Gaussian mass function for log10(m)
         :param rendering_radius_arcsec [arcsec]: sets the area around (center_x, center_y) where the GC's appear; GC's are
         distributed uniformly in a circle centered at (center_x, center_y) with this radius
-        :param gc_profile_args: the keyword arguments for the GC mass profile, must include "gamma", "gc_size_lightyear",
-        "r_core_fraction", or the logarithmic power law slope, the size of the gc in light years, and the size of the core
-        relative to the size of the GC
-        :param galaxy_Re: galaxy half-light radius
-        :param host_halo_mass: total mass of the host DM halo
-        :param f: the fraction M_gc / M_host from https://arxiv.org/pdf/1602.01105
-        :param coordinate_center_x: center of rendering area in arcsec
-        :param coordinate_center_y: center of rendering area in arcsec
+        :param gamma_mean: the mean logarithmic power-law slope for the GCs
+        :param gamma_sigma: half the width of the slope distribution around gamma mean, assuming a uniform distribution
+        :param gc_concentration_mean: the ratio of the GC core size to the total size, where the total size is defined as
+        the radius enclosing the stated mass of the GC
+        :param gc_concentration_sigma: half the width of the distribution around gc_concentration_mean
+        :param gc_size_mean: the typical radial extend of the GC in light years
+        (the mass is defined as the mass inside this radius)
+        :param gc_size_sigma: half the width of the uniform distribution of gc size
+        :param gc_surface_mass_density: the surface mass density of GCs in units of M_sun / kpc^2
+        :param center_x: center of rendering area in arcsec
+        :param center_y: center of rendering area in arcsec
         :return: an instance of Realization that includes the GC's
         """
         if isinstance(center_x, int) or isinstance(center_x, float):
             center_x = [center_x]
             center_y = [center_y]
         assert len(center_x) == len(center_y)
         GC_realization = None
+        lens_cosmo = self._realization.lens_cosmo
+        z = self._realization.lens_cosmo.z_lens
+        kpc_per_arcsec = self._realization.lens_cosmo.cosmo.kpc_proper_per_asec(z)
+        # determine number of GCs
+        integral = np.exp(np.log(10 ** log10_mgc_mean) + np.log(10 ** log10_mgc_sigma) ** 2 / 2)
+        mass_in_gc = np.pi * gc_surface_mass_density * (rendering_radius_arcsec * kpc_per_arcsec) ** 2
+        n = int(mass_in_gc / integral)
+        mfunc = Gaussian(n, log10_mgc_mean, log10_mgc_sigma)
         for x_center, y_center in zip(center_x, center_y):
-            n = self.number_globular_clusters(log10_mgc_mean, log10_mgc_sigma, rendering_radius_arcsec, galaxy_Re,
-                                              host_halo_mass, f)
-            mfunc = Gaussian(n, log10_mgc_mean, log10_mgc_sigma)
             m = mfunc.draw()
-            z = self._realization.lens_cosmo.z_lens
             uniform_spatial_distribution = Uniform(rendering_radius_arcsec, self._realization.geometry)
             x_kpc, y_kpc = uniform_spatial_distribution.draw(len(m), z, 1.0, x_center, y_center)
             x = x_kpc / self._realization.lens_cosmo.cosmo.kpc_proper_per_asec(z)
             y = y_kpc / self._realization.lens_cosmo.cosmo.kpc_proper_per_asec(z)
-            lens_cosmo = self._realization.lens_cosmo
             GCS = []
             for (m_gc, x_center_gc, y_center_gc) in zip(m, x, y):
+                gamma = np.random.uniform(gamma_mean - gamma_sigma, gamma_mean + gamma_sigma)
+                gc_concentration = np.random.uniform(gc_concentration_mean - gc_concentration_sigma,
+                                                     gc_concentration_mean + gc_concentration_sigma)
+                gc_size = np.random.uniform(gc_size_mean - gc_size_sigma, gc_size_mean + gc_size_sigma)
+                gc_size_lightyear = gc_size * (m_gc / 10 ** 5) ** (1/3)
+                gc_profile_args = {'gamma': gamma,
+                                   'gc_size_lightyear': gc_size_lightyear,
+                                   'gc_concentration': gc_concentration}
                 unique_tag = np.random.rand()
                 profile = GlobularCluster(m_gc, x_center_gc, y_center_gc, lens_cosmo.z_lens, lens_cosmo,
                                           gc_profile_args, unique_tag)
@@ -189,23 +204,6 @@ def add_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radiu
         new_realization = self._realization.join(GC_realization)
         return new_realization
 
-    def number_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radius_arcsec,
-                              galaxy_Re=10.0, host_halo_mass=10**13.3, f=3.4e-5):
-        """
-        Compute the number of globular clusters from their mass function (see documentation in add_globular_clusters)
-        :return:
-        """
-        m_gc = f * host_halo_mass
-        area = np.pi * (6 * galaxy_Re) ** 2 # physical kpc^2
-        surface_mass_density = m_gc / area
-        z = self._realization.lens_cosmo.z_lens
-        kpc_per_arcsec = self._realization.lens_cosmo.cosmo.kpc_proper_per_asec(z)
-        # assuming log-normal mass function
-        integral = np.exp(np.log(10 ** log10_mgc_mean) + np.log(10 ** log10_mgc_sigma) ** 2 / 2)
-        mass_in_gc = np.pi * surface_mass_density * (rendering_radius_arcsec * kpc_per_arcsec) ** 2
-        n = int(mass_in_gc / integral)
-        return n
-
     def toSIDM_single_halo(self, halo, t_c, subhalo_evolution_scaling):
         """
         Transform a single NFW profile to an SIDM profile based on the halo age, its status as a subhalo or a field halo,

tests/test_halos/test_globular_clusters.py

@@ -19,8 +19,8 @@ def test_lenstronomy_ID(self):
 
         mass = 10 ** 5
         args = {'gamma': 2.5,
-                'r_core_fraction': 0.05,
-                'gc_size_lightyear': 100}
+                'gc_size_lightyear': 100,
+                'gc_concentration': 50}
         profile = GlobularCluster(mass, 0.0, 0.0, self.zhalo, self.lens_cosmo,
                                   args, 1)
         lenstronomy_ID = profile.lenstronomy_ID
@@ -31,8 +31,8 @@ def test_lenstronomy_args(self):
         logM = 5.0
         mass = 10 ** logM
         args = {'gamma': 2.5,
-                'r_core_fraction': 0.05,
-                'gc_size_lightyear': 100}
+                'gc_size_lightyear': 100,
+                'gc_concentration': 50}
         profile = GlobularCluster(mass, 0.0, 0.0, self.zhalo, self.lens_cosmo,
                                   args, 1)
         lenstronomy_args, _ = profile.lenstronomy_params
@@ -42,15 +42,15 @@ def test_mass(self):
         logM = 5.0
         mass = 10 ** logM
         args = {'gamma': 2.5,
-                'r_core_fraction': 0.05,
-                'gc_size_lightyear': 100}
+                'gc_size_lightyear': 100,
+                'gc_concentration': 50}
         profile = GlobularCluster(mass, 0.0, 0.0, self.zhalo, self.lens_cosmo,
                                   args, 1)
         profile_args = profile.profile_args
-        rho0 = profile_args['rho0']
-        r_core = profile_args['r_core_kpc']
-        r_max = profile_args['gc_size']
-        gamma = profile_args['gamma']
+        rho0 = profile_args[0]
+        r_core = profile_args[3]
+        r_max = profile_args[1]
+        gamma = profile_args[2]
         total_mass = profile._prof.mass_3d(r_max, rho0, r_core, gamma)
         npt.assert_almost_equal(total_mass, mass)
 
@@ -59,7 +59,7 @@ def test_mass(self):
         sigma0 = lenstronomy_args[0]['sigma0']
         rcore = lenstronomy_args[0]['r_core']
         gamma = lenstronomy_args[0]['gamma']
-        r_max = profile_args['gc_size'] / kpc_per_arcsec
+        r_max = profile_args[1] / kpc_per_arcsec
         sigma_crit_mpc = profile.lens_cosmo.get_sigma_crit_lensing(profile.z, profile.lens_cosmo.z_source)
         sigma_crit_arcsec = sigma_crit_mpc * (0.001 * kpc_per_arcsec) ** 2
         total_mass = profile._prof.mass_3d(r_max, sigma0/rcore, rcore, gamma) * sigma_crit_arcsec

tests/test_preset_models.py

@@ -41,18 +41,10 @@ def test_CDM(self):
         self._test_default_infall_model(cdm, 'hybrid')
         self._test_default_infall_model(cdm3, 'direct')
 
-        gc_profile_args = {'gamma': 2.5,
-                           'gc_size_lightyear': 250,
-                           'r_core_fraction': 0.05}
         kwargs_globular_clusters = {'log10_mgc_mean': 5.0,
                                     'log10_mgc_sigma': 0.5,
                                     'rendering_radius_arcsec': 0.1,
-                                    'gc_profile_args': gc_profile_args,
-                                    'galaxy_Re': 6.0,
-                                    'host_halo_mass': 10**13.3,
-                                    'f': 3.4e-5,
-                                    'center_x': [0.5, 1.0],
-                                    'center_y': [-1, 0]}
+                                    }
         cdm = CDM(0.5, 1.5,
                   add_globular_clusters=True,
                   kwargs_globular_clusters=kwargs_globular_clusters)
@@ -69,18 +61,10 @@ def test_WDM(self):
         _ = wdm.lensing_quantities()
         _ = preset_model_from_name('WDM')
         self._test_default_infall_model(wdm, 'hybrid')
-        gc_profile_args = {'gamma': 2.5,
-                           'gc_size_lightyear': 250,
-                           'r_core_fraction': 0.05}
         kwargs_globular_clusters = {'log10_mgc_mean': 5.0,
                                     'log10_mgc_sigma': 0.5,
                                     'rendering_radius_arcsec': 0.1,
-                                    'gc_profile_args': gc_profile_args,
-                                    'galaxy_Re': 6.0,
-                                    'host_halo_mass': 10 ** 13.3,
-                                    'f': 3.4e-5,
-                                    'center_x': [0.5, 1.0],
-                                    'center_y': [-1, 0]}
+                                    }
         wdm = WDM(0.5, 1.5,
                   7.7,
                   add_globular_clusters=True,
@@ -143,18 +127,10 @@ def test_WDM_general(self):
         _ = wdm.lensing_quantities()
         wdm = func(0.5, 1.5, 7.7, -2.0, truncation_model_subhalos='TRUNCATION_GALACTICUS')
         self._test_default_infall_model(wdm, 'hybrid')
-        gc_profile_args = {'gamma': 2.5,
-                           'gc_size_lightyear': 250,
-                           'r_core_fraction': 0.05}
         kwargs_globular_clusters = {'log10_mgc_mean': 5.0,
                                     'log10_mgc_sigma': 0.5,
                                     'rendering_radius_arcsec': 0.1,
-                                    'gc_profile_args': gc_profile_args,
-                                    'galaxy_Re': 6.0,
-                                    'host_halo_mass': 10 ** 13.3,
-                                    'f': 3.4e-5,
-                                    'center_x': [0.5, 1.0],
-                                    'center_y': [-1, 0]}
+                                    }
         wdm = func(0.5, 1.5, 7.7, -2.5,
                   add_globular_clusters=True,
                   kwargs_globular_clusters=kwargs_globular_clusters)