update sidm preset model

pyHalo/PresetModels/mbh.py

@@ -22,6 +22,8 @@ def CDM_plus_BH(z_lens,
                  kwargs_globular_clusters=None):
     """
     Add a population of black holes to a CDM realization
+    some literature:
+    - https://arxiv.org/pdf/0709.0529 (evolution of Mbh seeds)
     :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

pyHalo/PresetModels/sidm.py

@@ -152,3 +152,34 @@ class for details
     extension = RealizationExtensions(cdm)
     sidm = extension.toSIDM(log10_mass_ranges, collapse_timescales, subhalo_time_scaling)
     return sidm
+
+def SIDM_parametric_fixedbins(z_lens, z_source, t0_mlow_75, t0_75_85, t0_85_10, subhalo_time_scaling=1.0,
+        sigma_sub=0.025, log10_sigma_sub=None, log_mlow=6., log_mhigh=10.,
+        concentration_model_subhalos='DIEMERJOYCE19', kwargs_concentration_model_subhalos={},
+        concentration_model_fieldhalos='DIEMERJOYCE19', kwargs_concentration_model_fieldhalos={},
+        truncation_model_subhalos='TRUNCATION_GALACTICUS', kwargs_truncation_model_subhalos={},
+        infall_redshift_model='HYBRID_INFALL', kwargs_infall_model={},
+        truncation_model_fieldhalos='TRUNCATION_RN', kwargs_truncation_model_fieldhalos={},
+        shmf_log_slope=-1.9, cone_opening_angle_arcsec=6., log_m_host=13.3,  r_tidal=0.25,
+        LOS_normalization=1.0, geometry_type='DOUBLE_CONE', kwargs_cosmo=None):
+    """
+    This function uses the SIDM_parametric model with a preset list of logarithmic mass bins 10^log_mlow - 10^7.5,
+    10^7.5-10^8.5, 10^8.5-10^10. The three parameters t0_mlow_75, t0_75_85, t0_85_10 set the collapse timescale in each
+    bin. For detailed documentation for the rest of the arguments, see the documentation in the SIDM_parametric function
+    :param t0_mlow_75: core collapse timescale log10(Gyr) in the mass range 10^log_mlow - 10^7.5
+    :param t0_75_85: core collapse timescale log10(Gyr) in the mass range 10^7.5 - 10^8.5
+    :param t0_85_10: core collapse timescale log10(Gyr) in the mass range 10^8.5 - 10^10
+    """
+    if log_mlow > 7.5:
+        raise Exception('to use this function log_mlow must be < 7.5')
+    log10_mass_ranges = [[log_mlow, 7.5], [7.5, 8.5], [8.5, 10.0]]
+    collapse_timescales = [t0_mlow_75, t0_75_85, t0_85_10]
+    return SIDM_parametric(z_lens, z_source, log10_mass_ranges, collapse_timescales, subhalo_time_scaling,
+        sigma_sub, log10_sigma_sub, log_mlow, log_mhigh,
+        concentration_model_subhalos, kwargs_concentration_model_subhalos,
+        concentration_model_fieldhalos, kwargs_concentration_model_fieldhalos,
+        truncation_model_subhalos, kwargs_truncation_model_subhalos,
+        infall_redshift_model, kwargs_infall_model,
+        truncation_model_fieldhalos, kwargs_truncation_model_fieldhalos,
+        shmf_log_slope, cone_opening_angle_arcsec, log_m_host,  r_tidal,
+        LOS_normalization, geometry_type, kwargs_cosmo)

pyHalo/preset_models.py

@@ -25,6 +25,9 @@ def preset_model_from_name(name):
     elif name == 'SIDM_parametric':
         from pyHalo.PresetModels.sidm import SIDM_parametric
         return SIDM_parametric
+    elif name == 'SIDM_parametric_fixedbins':
+        from pyHalo.PresetModels.sidm import SIDM_parametric_fixedbins
+        return SIDM_parametric_fixedbins
     elif name == 'ULDM':
         from pyHalo.PresetModels.uldm import ULDM
         return ULDM

tests/test_preset_models.py

@@ -96,6 +96,10 @@ def test_SIDM(self):
         _ = preset_model_from_name('SIDM_parametric')
         self._test_default_infall_model(sidm, 'hybrid')
 
+        model = preset_model_from_name('SIDM_parametric_fixedbins')
+        realization = model(0.5, 1.5, 10.0, 5.0, 1.0)
+        _ = realization.lensing_quantities()
+
     def test_ULDM(self):
 
         flucs_shape = 'ring'
@@ -123,7 +127,7 @@ def test_WDM_mixed(self):
         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,