fix error in NFW profile sampling

pyHalo/Halos/concentration.py

@@ -24,8 +24,6 @@ def nfw_concentration(self, m, z):
         Evaluates the concentration of a halo of mass 'm' at redshift z
         :param M: halo mass [M_sun]
         :param z: halo redshift
-        :param scatter: bool; add log-normal scatter to concentration
-        :param scatter_amplitude: the amount of scatter in dex, assumes log-normal distribution
         :return:
         """
         c = self._evaluate_concentration(m, z)

pyHalo/Halos/tidal_truncation.py

@@ -199,3 +199,44 @@ def _make_params_in_bounds_tau_evaluate(self, point):
         log10mass_loss_fraction = max(-1.5, log10mass_loss_fraction)
         log10mass_loss_fraction = min(-0.01, log10mass_loss_fraction)
         return (log10c, log10mass_loss_fraction)
+
+class TruncateMeanDensity(object):
+
+    def __init__(self, lens_cosmo, median_rt_over_rs=1.0, c_power=3.0):
+        """
+
+        :param lens_cosmo:
+        :param median_rt_over_rs:
+        :param c_power:
+        """
+        self._norm = median_rt_over_rs
+        self._cpower = c_power
+        self.lens_cosmo = lens_cosmo
+        self._concentration_cdm = ConcentrationDiemerJoyce(lens_cosmo.cosmo,
+                                                           scatter=False)
+
+    def truncation_radius_halo(self, halo):
+
+        """
+        Thiis method computess the truncation radius using the class attributes of an instance of Halo
+        :param halo: an instance of halo
+        :return: the truncation radius in physical kpc
+        """
+        c_median = self._concentration_cdm.nfw_concentration(halo.mass, halo.z_eval)
+        c_actual = halo.c
+        halo_rpericenter = halo.rperi_units_r200
+        return self.truncation_radius(halo.mass, halo.z, c_median, c_actual, halo_rpericenter)
+
+    def truncation_radius(self, halo_mass, halo_redshift, c_median, c_actual, r_peri):
+        """
+
+        :param halo_mass:
+        :param halo_redshift:
+        :param c_median:
+        :param c_actual:
+        :param r_peri:
+        :return:
+        """
+        rt_over_rs = self._norm * (c_actual / c_median) ** self._cpower * (r_peri / 0.5)
+        _, rs, _ = self.lens_cosmo.NFW_params_physical(halo_mass, c_actual, halo_redshift)
+        return rs * rt_over_rs

pyHalo/Rendering/SpatialDistributions/nfw.py

@@ -1,13 +1,11 @@
 import numpy as np
 from lenstronomy.LensModel.Profiles.cnfw import CNFW
 from pyHalo.Rendering.SpatialDistributions.base import SpatialDistributionBase
-import inspect
 from pyHalo.Halos.concentration import ConcentrationDiemerJoyce
 from pyHalo.utilities import inverse_transform_sampling
+from scipy.interpolate import interp1d
 
 
-local_path = inspect.getfile(inspect.currentframe())[0:-11] + 'nfw_tables/'
-#
 class ProjectedNFW(SpatialDistributionBase):
 
     """
@@ -26,6 +24,7 @@ def __init__(self, rmax2d_arcsec, rs_arcsec, r_core_arcsec, r_200_arcsec, arcsec
         :param rmax3d: the virial radius of the host dark matter halo [kpc]
         :param r_core_parent: the core radius of the host dark matter halo [kpc]
         """
+
         self._rs_arcsec = rs_arcsec
         self._rmax2d = rmax2d_arcsec
         self._rmax3d = r_200_arcsec
@@ -60,18 +59,44 @@ def draw(self, N, z_plane=None):
         if N == 0:
             return [], [], []
 
-        x2d = np.linspace(1e-3 * self._rmax2d, self._rmax2d, 50000)
-        function_2d = self._cnfw_profile.density_2d
+        x2d = np.linspace(1e-3 * self._rmax2d, self._rmax2d, 20000)
         args = (0.0, self._rs_arcsec, 1.0, self._rcore_arcsec)
-        r2d_arcsec = inverse_transform_sampling(x2d, function_2d, args, N)
+        rho2d = self._cnfw_profile.density_2d(x2d, *args)
+        rho2d_integral = [rho2d[i] * (x2d[i+1]**2 - x2d[i]**2) for i in range(0, len(x2d)-1)]
+        function_2d = interp1d(x2d[0:-1], rho2d_integral)
+        r2d_arcsec = inverse_transform_sampling(x2d[0:-1], function_2d, (), N)
         r2d_kpc = r2d_arcsec * self._arcsec_to_kpc
         theta = np.random.uniform(0, 2*np.pi, N)
         x_kpc, y_kpc = r2d_kpc * np.cos(theta), r2d_kpc * np.sin(theta)
 
-        x3d = np.linspace(1e-3 * self._rmax2d, self._rmax3d, 50000)
-        function_3d = self._cnfw_profile.density
+        x3d = np.linspace(1e-3 * self._rmax2d, self._rmax3d, 20000)
         args = (self._rs_arcsec, 1.0, self._rcore_arcsec)
-        r3d_arcsec = inverse_transform_sampling(x3d, function_3d, args, N)
+        rho3d = self._cnfw_profile.density(x3d, *args)
+        rho3d_integral = [rho3d[i] * (x3d[i+1]**3 - x3d[i]**3) for i in range(0, len(x3d)-1)]
+        function_3d = interp1d(x3d[0:-1], rho3d_integral)
+        r3d_arcsec = inverse_transform_sampling(x3d[0:-1], function_3d, (), N)
         r3d_kpc = r3d_arcsec * self._arcsec_to_kpc
 
         return x_kpc, y_kpc, r3d_kpc
+
+    # def draw(self, N, z_plane=None):
+    #
+    #     if N == 0:
+    #         return [], [], []
+    #
+    #     x2d = np.linspace(1e-3 * self._rmax2d, self._rmax2d, 50000)
+    #     function_2d = self._cnfw_profile.density_2d
+    #     args = (0.0, self._rs_arcsec, 1.0, self._rcore_arcsec)
+    #     r2d_arcsec = inverse_transform_sampling(x2d, function_2d, args, N)
+    #     r2d_kpc = r2d_arcsec * self._arcsec_to_kpc
+    #     theta = np.random.uniform(0, 2*np.pi, N)
+    #     x_kpc, y_kpc = r2d_kpc * np.cos(theta), r2d_kpc * np.sin(theta)
+    #
+    #     x3d = np.linspace(1e-3 * self._rmax2d, self._rmax3d, 50000)
+    #     function_3d = self._cnfw_profile.density
+    #     args = (self._rs_arcsec, 1.0, self._rcore_arcsec)
+    #     r3d_arcsec = inverse_transform_sampling(x3d, function_3d, args, N)
+    #     r3d_kpc = r3d_arcsec * self._arcsec_to_kpc
+    #
+    #     return x_kpc, y_kpc, r3d_kpc
+