Merge pull request #67 from dangilman/sidm_evolution

fix gcs

pyHalo/Halos/HaloModels/PTMass.py → pyHalo/Halos/HaloModels/blackhole.py

@@ -1,7 +1,7 @@
 import numpy as np
 from pyHalo.Halos.halo_base import Halo
 
-class PTMass(Halo):
+class BlackHole(Halo):
     """
     Class that defines a point mass object in the lens model
     """
@@ -14,8 +14,9 @@ def __init__(self, mass, x, y, r3d, z,
         self._truncation_class = truncation_class
         self._concentration_clss = concentration_class
         mdef = 'PT_MASS'
-        super(PTMass, self).__init__(mass, x, y, r3d, mdef, z, sub_flag,
-                                            lens_cosmo_instance, args, unique_tag)
+        super(BlackHole, self).__init__(mass, x, y, r3d, mdef, z, sub_flag,
+                                        lens_cosmo_instance, args, unique_tag,
+                                        fixed_position=True)
 
     @property
     def lenstronomy_ID(self):

pyHalo/Halos/HaloModels/powerlaw.py

@@ -8,7 +8,8 @@ class PowerLawSubhalo(Halo):
     The base class for a halo modeled as a power law profile
     """
     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):
         """
         See documentation in base class (Halos/halo_base.py)
         """
@@ -139,20 +140,28 @@ def __init__(self, mass, x, y, z, lens_cosmo_instance, args, unique_tag):
         self._prof = SPLCORE()
         self._lens_cosmo = lens_cosmo_instance
         mdef = 'SPL_CORE'
-        super(GlobularCluster, self).__init__(mass, x, y, None, mdef, z, True,
-                                              lens_cosmo_instance, args, unique_tag)
+        is_subhalo = False
+        super(GlobularCluster, self).__init__(mass, x, y, None, mdef, z, is_subhalo,
+                                              lens_cosmo_instance, args, unique_tag,
+                                              fixed_position=True)
 
     @property
     def lenstronomy_params(self):
         """
         See documentation in base class (Halos/halo_base.py)
         """
         if not hasattr(self, '_lenstronomy_args'):
-
-            profile_args = self.profile_args
-            rho0 = profile_args['rho0']
-            r_core_arcsec = profile_args['r_core_arcsec']
-            gamma = profile_args['gamma']
+            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']
+            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,
                                        'r_core': r_core_arcsec}]
@@ -164,17 +173,16 @@ def profile_args(self):
         See documentation in base class (Halos/halo_base.py)
         """
         if not hasattr(self, '_profile_args'):
-            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']
-            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 = gc_size_lightyear_0 * (self.mass / 10 ** 5) ** (1 / 3) * kpc_per_lightyear
-            r_core_arcsec = r_core_fraction * gc_size
-            rho0 = self.mass / self._prof.mass_3d(gc_size, sigma_crit_arcsec, r_core_arcsec, gamma)
-            self._profile_args = {'rho0': rho0, 'gc_size': gc_size, 'gamma': gamma, 'r_core_arcsec': r_core_arcsec}
+            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
+            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_arcsec': r_core_kpc}
         return self._profile_args
 
     @property

pyHalo/realization_extensions.py

@@ -64,7 +64,9 @@ def add_globular_clusters(self, log10_mgc_mean, log10_mgc_sigma, rendering_radiu
             m = mfunc.draw()
             z = self._realization.lens_cosmo.z_lens
             uniform_spatial_distribution = Uniform(rendering_radius_arcsec, self._realization.geometry)
-            x, y = uniform_spatial_distribution.draw(len(m), z, x_center, y_center)
+            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):

pyHalo/single_realization.py

@@ -6,7 +6,7 @@
 from pyHalo.Halos.HaloModels.generalized_nfw import GeneralNFWFieldHalo, GeneralNFWSubhalo
 from pyHalo.Halos.HaloModels.powerlaw import PowerLawFieldHalo, PowerLawSubhalo
 from pyHalo.Halos.HaloModels.PsuedoJaffe import PJaffeSubhalo
-from pyHalo.Halos.HaloModels.PTMass import PTMass
+from pyHalo.Halos.HaloModels.blackhole import BlackHole
 from pyHalo.Halos.HaloModels.ULDM import ULDMFieldHalo, ULDMSubhalo
 from pyHalo.Halos.HaloModels.NFW_core_trunc import TNFWCFieldHaloSIDM, TNFWCSubhaloSIDM
 from pyHalo.Halos.HaloModels.gaussianhalo import GaussianHalo
@@ -667,7 +667,7 @@ def _load_halo_model(mass, x, y, r3d, mdef, z, is_subhalo,
             else:
                 model = TNFWCFieldHaloSIDM
         elif mdef == 'PT_MASS':
-            model = PTMass
+            model = BlackHole
         elif mdef == 'PJAFFE':
             model = PJaffeSubhalo
         elif mdef == 'ULDM':

tests/test_halos/test_globular_clusters.py

@@ -36,17 +36,34 @@ def test_lenstronomy_args(self):
         profile = GlobularCluster(mass, 0.0, 0.0, self.zhalo, self.lens_cosmo,
                                   args, 1)
         lenstronomy_args, _ = profile.lenstronomy_params
-        profile_lenstronomy = SPLCORE()
+
+    def test_mass(self):
+
+        logM = 5.0
+        mass = 10 ** logM
+        args = {'gamma': 2.5,
+                'r_core_fraction': 0.05,
+                'gc_size_lightyear': 100}
+        profile = GlobularCluster(mass, 0.0, 0.0, self.zhalo, self.lens_cosmo,
+                                  args, 1)
         profile_args = profile.profile_args
-        gc_size = profile_args['gc_size']
         rho0 = profile_args['rho0']
         r_core = profile_args['r_core_arcsec']
+        r_max = profile_args['gc_size']
         gamma = profile_args['gamma']
-        mass = profile_lenstronomy.mass_3d(gc_size, rho0, r_core, gamma)
-        sigma_crit_mpc = self.lens_cosmo.get_sigma_crit_lensing(profile.z, self.lens_cosmo.z_source)
-        kpc_per_arcsec = self.lens_cosmo.cosmo.kpc_proper_per_asec(profile.z)
+        total_mass = profile._prof.mass_3d(r_max, rho0, r_core, gamma)
+        npt.assert_almost_equal(total_mass, mass)
+
+        kpc_per_arcsec = profile.lens_cosmo.cosmo.kpc_proper_per_asec(profile.z)
+        lenstronomy_args = profile.lenstronomy_params[0]
+        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
+        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
-        npt.assert_almost_equal(logM, np.log10(mass * sigma_crit_arcsec))
+        total_mass = profile._prof.mass_3d(r_max, sigma0/rcore, rcore, gamma) * sigma_crit_arcsec
+        npt.assert_almost_equal(total_mass, mass)
 
 if __name__ == '__main__':
     pytest.main()

tests/test_halos/test_point_mass.py

@@ -1,6 +1,6 @@
 import numpy.testing as npt
 import numpy as np
-from pyHalo.Halos.HaloModels.PTMass import PTMass
+from pyHalo.Halos.HaloModels.blackhole import BlackHole
 from astropy.cosmology import FlatLambdaCDM
 from pyHalo.Halos.lens_cosmo import LensCosmo
 from pyHalo.Cosmology.cosmology import Cosmology
@@ -17,7 +17,7 @@ def setup_method(self):
         self.truncation_class = None
         self.concentration_class = None
         kwargs_profile = {}
-        self.ptmass = PTMass(10**9, 1.0, 1.0, None, self.zhalo, False, self.lens_cosmo, kwargs_profile, None, None, 1.0)
+        self.ptmass = BlackHole(10 ** 9, 1.0, 1.0, None, self.zhalo, False, self.lens_cosmo, kwargs_profile, None, None, 1.0)
 
     def test_lenstronomy_ID(self):
 
@@ -38,8 +38,8 @@ def test_thetaE(self):
         is_subhalo = False
         kwargs_profile = {}
         unique_tag = 0.1
-        ptmass = PTMass(m, x, y, r3d, self.zhalo, is_subhalo, self.lens_cosmo, kwargs_profile,
-                                      self.truncation_class, self.concentration_class, unique_tag)
+        ptmass = BlackHole(m, x, y, r3d, self.zhalo, is_subhalo, self.lens_cosmo, kwargs_profile,
+                           self.truncation_class, self.concentration_class, unique_tag)
         kwargs, _ = ptmass.lenstronomy_params
         theta_E = kwargs[0]['theta_E']
         theta_E_theory = self.lens_cosmo.point_mass_factor_z(ptmass.z) * np.sqrt(10**8)