"""
This module contains the `NFWProfile` class,
which is used to model the spatial distribution of mass and/or galaxies
inside dark matter halos according to the fitting function introduced in
Navarry, Frenk and White (1995), `arXiv:9508025 <http://arxiv.org/abs/astro-ph/9508025/>`_.
a sub-class of `~halotools.empirical_models.AnalyticDensityProf`.
"""
from __future__ import division, print_function, absolute_import, unicode_literals
import numpy as np
from .conc_mass import direct_from_halo_catalog, dutton_maccio14
from .kernels import nfw_dimensionless_mass_density, nfw_cumulative_mass_PDF
from .kernels import standalone_mc_generate_nfw_radial_positions
from ...profile_model_template import AnalyticDensityProf
from ..... import model_defaults
from ......sim_manager import sim_defaults
__all__ = ("NFWProfile",)
__author__ = ("Andrew Hearin", "Benedikt Diemer")
[docs]
class NFWProfile(AnalyticDensityProf):
r"""Model for the spatial distribution of mass
and/or galaxies residing in an NFW halo profile,
based on Navarro, Frenk and White (1995),
`arXiv:9508025 <http://arxiv.org/abs/astro-ph/9508025/>`_.
For a review of the mathematics underlying the NFW profile,
including descriptions of how the relevant equations are
implemented in the Halotools code base, see :ref:`nfw_profile_tutorial`.
"""
def __init__(
self,
cosmology=sim_defaults.default_cosmology,
redshift=sim_defaults.default_redshift,
mdef=model_defaults.halo_mass_definition,
conc_mass_model=model_defaults.conc_mass_model,
concentration_key=model_defaults.concentration_key,
halo_boundary_key=None,
**kwargs
):
r"""
Parameters
----------
cosmology : object, optional
Instance of an astropy `~astropy.cosmology`.
Default cosmology is set in
`~halotools.sim_manager.sim_defaults`.
redshift : float, optional
Default is set in `~halotools.sim_manager.sim_defaults`.
mdef: str, optional
String specifying the halo mass definition, e.g., 'vir' or '200m'.
Default is set in `~halotools.empirical_models.model_defaults`.
halo_boundary_key : str, optional
Default behavior is to use the column associated with the input mdef.
conc_mass_model : string or callable, optional
Specifies the function used to model the relation between
NFW concentration and halo mass.
Can either be a custom-built callable function,
or one of the following strings:
``dutton_maccio14``, ``direct_from_halo_catalog``.
concentration_key : string, optional
Column name of the halo catalog storing NFW concentration.
This argument is only relevant when ``conc_mass_model``
is set to ``direct_from_halo_catalog``. In such a case,
the default value is ``halo_nfw_conc``,
which is consistent with all halo catalogs provided by Halotools
but may differ from the convention adopted in custom halo catalogs.
Examples
--------
>>> nfw = NFWProfile()
"""
AnalyticDensityProf.__init__(
self, cosmology, redshift, mdef, halo_boundary_key=halo_boundary_key
)
self.gal_prof_param_keys = ["conc_NFWmodel"]
self.halo_prof_param_keys = ["conc_NFWmodel"]
self.publications = ["arXiv:9611107", "arXiv:0002395"]
self._initialize_conc_mass_behavior(
conc_mass_model, concentration_key=concentration_key
)
def _initialize_conc_mass_behavior(self, conc_mass_model, **kwargs):
if conc_mass_model == "direct_from_halo_catalog":
self.concentration_key = kwargs.get(
"concentration_key", model_defaults.concentration_key
)
self.conc_mass_model = conc_mass_model
[docs]
def conc_NFWmodel(self, *args, **kwargs):
r"""NFW concentration as a function of halo mass.
Parameters
----------
prim_haloprop : array, optional
Array storing the mass-like variable, e.g., ``halo_mvir``.
If ``prim_haloprop`` is not passed,
then ``table`` keyword argument must be passed.
table : object, optional
`~astropy.table.Table` storing the halo catalog.
If your NFW model is based on the virial definition,
then ``halo_mvir`` must appear in the input table,
and likewise for other halo mass definitions.
If ``table`` is not passed,
then ``prim_haloprop`` keyword argument must be passed.
Returns
-------
conc : array_like
Concentrations of the input halos.
Note that concentrations will be clipped to their min/max permitted
values set in the `~halotools.empirical_models.model_defaults` module.
The purpose of this clipping is to ensure stable results during
mock galaxy population. Due to this clipping,
the behavior of the `conc_NFWmodel` function
is different from the concentration-mass relation that underlies it.
Examples
---------
In the examples below, we'll demonstrate the various ways to use the
`~halotools.empirical_models.NFWProfile.conc_NFWmodel` function, depending
on the initial choice for the ``conc_mass_model``.
>>> fake_masses = np.logspace(12, 15, 10)
If you use the ``direct_from_halo_catalog`` option, you must pass a
``table`` argument storing a `~astropy.table.Table` with a column name
for the halo mass that is consistent with your chosen halo mass definition:
>>> from astropy.table import Table
>>> nfw = NFWProfile(conc_mass_model='direct_from_halo_catalog', mdef='vir')
>>> fake_conc = np.zeros_like(fake_masses) + 5.
>>> fake_halo_table = Table({'halo_mvir': fake_masses, 'halo_nfw_conc': fake_conc})
>>> model_conc = nfw.conc_NFWmodel(table=fake_halo_table)
In case your halo catalog uses a different keyname from the Halotools
default ``halo_nfw_conc``:
>>> nfw = NFWProfile(conc_mass_model='direct_from_halo_catalog', mdef='vir', concentration_key='my_conc_keyname')
>>> fake_halo_table = Table({'halo_mvir': fake_masses, 'my_conc_keyname': fake_conc})
>>> model_conc = nfw.conc_NFWmodel(table=fake_halo_table)
One of the available options provided by Halotools is ``dutton_maccio14``.
With this option, you can either pass in a ``table`` argument, or alternatively
an array of masses via the ``prim_haloprop`` argument:
>>> nfw = NFWProfile(conc_mass_model='dutton_maccio14')
>>> fake_halo_table = Table({'halo_mvir': fake_masses, 'halo_nfw_conc': fake_conc})
>>> model_conc = nfw.conc_NFWmodel(table=fake_halo_table)
>>> model_conc = nfw.conc_NFWmodel(prim_haloprop=fake_masses)
Finally, you may also have chosen to define your own concentration-mass relation.
If so, your function must at a minimum accept a ``table`` keyword argument.
Below we give a trivial example of using the identity function:
>>> def identity_func(*args, **kwargs): return kwargs['table']['halo_mvir']
>>> nfw = NFWProfile(conc_mass_model=identity_func, mdef='vir')
>>> fake_halo_table = Table({'halo_mvir': fake_masses})
>>> model_conc = nfw.conc_NFWmodel(table=fake_halo_table)
"""
if self.conc_mass_model == "direct_from_halo_catalog":
try:
table = kwargs["table"]
except KeyError:
msg = (
"Must pass ``table`` argument to the ``conc_NFWmodel`` function\n"
"when ``conc_mass_model`` is set to ``direct_from_halo_catalog``\n"
)
raise KeyError(msg)
result = direct_from_halo_catalog(
table=table, concentration_key=self.concentration_key
)
elif self.conc_mass_model == "dutton_maccio14":
msg = (
"Must either pass a ``prim_haloprop`` argument, \n"
"or a ``table`` argument with an astropy Table that has the ``{0}`` key"
)
try:
mass = kwargs["table"][self.prim_haloprop_key]
except:
try:
mass = kwargs["prim_haloprop"]
except:
raise KeyError(msg.format(self.prim_haloprop_key))
result = dutton_maccio14(mass, self.redshift)
else:
result = self.conc_mass_model(*args, **kwargs)
cmin = model_defaults.min_permitted_conc
cmax = model_defaults.max_permitted_conc
result = np.where(result < cmin, cmin, result)
result = np.where(result > cmax, cmax, result)
return result
[docs]
def dimensionless_mass_density(self, scaled_radius, conc):
r"""
Physical density of the NFW halo scaled by the density threshold of the mass definition.
The `dimensionless_mass_density` is defined as
:math:`\tilde{\rho}_{\rm prof}(\tilde{r}) \equiv \rho_{\rm prof}(\tilde{r}) / \rho_{\rm thresh}`,
where :math:`\tilde{r}\equiv r/R_{\Delta}`.
For an NFW halo,
:math:`\tilde{\rho}_{\rm NFW}(\tilde{r}, c) = \frac{c^{3}/3g(c)}{c\tilde{r}(1 + c\tilde{r})^{2}},`
where :math:`g(x) \equiv \log(1+x) - x / (1+x)` is computed using the `g` function.
The quantity :math:`\rho_{\rm thresh}` is a function of
the halo mass definition, cosmology and redshift,
and is computed via the
`~halotools.empirical_models.profile_helpers.density_threshold` function.
The quantity :math:`\rho_{\rm prof}` is the physical mass density of the
halo profile and is computed via the `mass_density` function.
See :ref:`nfw_spatial_profile_derivations` for a derivation of this expression.
Parameters
-----------
scaled_radius : array_like
Halo-centric distance *r* scaled by the halo boundary :math:`R_{\Delta}`, so that
:math:`0 <= \tilde{r} \equiv r/R_{\Delta} <= 1`. Can be a scalar or numpy array.
conc : array_like
Value of the halo concentration. Can either be a scalar, or a numpy array
of the same dimension as the input ``scaled_radius``.
Returns
-------
dimensionless_density: array_like
Dimensionless density of a dark matter halo
at the input ``scaled_radius``, normalized by the
`~halotools.empirical_models.profile_helpers.density_threshold`
:math:`\rho_{\rm thresh}` for the
halo mass definition, cosmology, and redshift.
Result is an array of the dimension as the input ``scaled_radius``.
"""
return nfw_dimensionless_mass_density(scaled_radius, conc)
[docs]
def mass_density(self, radius, mass, conc):
r"""
Physical density of the halo at the input radius,
given in units of :math:`h^{3}/{\rm Mpc}^{3}`.
Parameters
-----------
radius : array_like
Halo-centric distance in Mpc/h units; can be a scalar or numpy array
mass : array_like
Total mass of the halo; can be a scalar or numpy array of the same
dimension as the input ``radius``.
conc : array_like
Value of the halo concentration. Can either be a scalar, or a numpy array
of the same dimension as the input ``radius``.
Returns
-------
density: array_like
Physical density of a dark matter halo of the input ``mass``
at the input ``radius``. Result is an array of the
dimension as the input ``radius``, reported in units of :math:`h^{3}/Mpc^{3}`.
Examples
--------
>>> model = NFWProfile()
>>> Npts = 100
>>> radius = np.logspace(-2, -1, Npts)
>>> mass = np.zeros(Npts) + 1e12
>>> conc = 5
>>> result = model.mass_density(radius, mass, conc)
>>> concarr = np.linspace(1, 100, Npts)
>>> result = model.mass_density(radius, mass, concarr)
Notes
------
See :ref:`halo_profile_definitions` for derivations and implementation details.
"""
return AnalyticDensityProf.mass_density(self, radius, mass, conc)
[docs]
def cumulative_mass_PDF(self, scaled_radius, conc):
r"""
Analytical result for the fraction of the total mass
enclosed within dimensionless radius of an NFW halo,
:math:`P_{\rm NFW}(<\tilde{r}) \equiv M_{\Delta}(<\tilde{r}) / M_{\Delta} = g(c\tilde{r})/g(\tilde{r}),`
where :math:`g(x) \equiv \int_{0}^{x}dy\frac{y}{(1+y)^{2}} = \log(1+x) - x / (1+x)` is computed
using `g`, and where :math:`\tilde{r} \equiv r / R_{\Delta}`.
See :ref:`nfw_cumulative_mass_pdf_derivation` for a derivation of this expression.
Parameters
-------------
scaled_radius : array_like
Halo-centric distance *r* scaled by the halo boundary :math:`R_{\Delta}`, so that
:math:`0 <= \tilde{r} \equiv r/R_{\Delta} <= 1`. Can be a scalar or numpy array.
conc : array_like
Value of the halo concentration. Can either be a scalar, or a numpy array
of the same dimension as the input ``scaled_radius``.
Returns
-------------
p: array_like
The fraction of the total mass enclosed
within the input ``scaled_radius``, in :math:`M_{\odot}/h`;
has the same dimensions as the input ``scaled_radius``.
Examples
--------
>>> model = NFWProfile()
>>> Npts = 100
>>> scaled_radius = np.logspace(-2, 0, Npts)
>>> conc = 5
>>> result = model.cumulative_mass_PDF(scaled_radius, conc)
>>> concarr = np.linspace(1, 100, Npts)
>>> result = model.cumulative_mass_PDF(scaled_radius, concarr)
"""
scaled_radius = np.where(scaled_radius > 1, 1, scaled_radius)
scaled_radius = np.where(scaled_radius < 0, 0, scaled_radius)
return nfw_cumulative_mass_PDF(scaled_radius, conc)
[docs]
def enclosed_mass(self, radius, total_mass, conc):
r"""
The mass enclosed within the input radius,
:math:`M(<r) = 4\pi\int_{0}^{r}dr'r'^{2}\rho(r)`.
Parameters
-----------
radius : array_like
Halo-centric distance in Mpc/h units; can be a scalar or numpy array
total_mass : array_like
Total mass of the halo; can be a scalar or numpy array of the same
dimension as the input ``radius``.
conc : array_like
Value of the halo concentration. Can either be a scalar, or a numpy array
of the same dimension as the input ``radius``.
Returns
----------
enclosed_mass: array_like
The mass enclosed within radius r, in :math:`M_{\odot}/h`;
has the same dimensions as the input ``radius``.
Examples
--------
>>> model = NFWProfile()
>>> Npts = 100
>>> radius = np.logspace(-2, -1, Npts)
>>> total_mass = np.zeros(Npts) + 1e12
>>> conc = 5
>>> result = model.enclosed_mass(radius, total_mass, conc)
>>> concarr = np.linspace(1, 100, Npts)
>>> result = model.enclosed_mass(radius, total_mass, concarr)
Notes
------
See :ref:`halo_profile_definitions` for derivations and implementation details.
"""
return AnalyticDensityProf.enclosed_mass(self, radius, total_mass, conc)
[docs]
def virial_velocity(self, total_mass):
r"""The circular velocity evaluated at the halo boundary,
:math:`V_{\rm vir} \equiv \sqrt{GM_{\rm halo}/R_{\rm halo}}`.
Parameters
--------------
total_mass : array_like
Total mass of the halo; can be a scalar or numpy array.
Returns
--------
vvir : array_like
Virial velocity in km/s.
Examples
--------
>>> model = NFWProfile()
>>> Npts = 100
>>> mass_array = np.logspace(11, 15, Npts)
>>> vvir_array = model.virial_velocity(mass_array)
Notes
------
See :ref:`halo_profile_definitions` for derivations and implementation details.
"""
return AnalyticDensityProf.virial_velocity(self, total_mass)
[docs]
def circular_velocity(self, radius, total_mass, conc):
r"""
The circular velocity, :math:`V_{\rm cir} \equiv \sqrt{GM(<r)/r}`,
as a function of halo-centric distance r.
Parameters
--------------
radius : array_like
Halo-centric distance in Mpc/h units; can be a scalar or numpy array
total_mass : array_like
Total mass of the halo; can be a scalar or numpy array of the same
dimension as the input ``radius``.
conc : array_like
Value of the halo concentration. Can either be a scalar, or a numpy array
of the same dimension as the input ``radius``.
Returns
----------
vc: array_like
The circular velocity in km/s; has the same dimensions as the input ``radius``.
Examples
--------
>>> model = NFWProfile()
>>> Npts = 100
>>> radius = np.logspace(-2, -1, Npts)
>>> total_mass = np.zeros(Npts) + 1e12
>>> conc = 5
>>> result = model.circular_velocity(radius, total_mass, conc)
>>> concarr = np.linspace(1, 100, Npts)
>>> result = model.circular_velocity(radius, total_mass, concarr)
Notes
------
See :ref:`halo_profile_definitions` for derivations and implementation details.
"""
return AnalyticDensityProf.circular_velocity(self, radius, total_mass, conc)
[docs]
def rmax(self, total_mass, conc):
r"""Radius at which the halo attains its maximum circular velocity,
:math:`R_{\rm max}^{\rm NFW} = 2.16258R_{\Delta}/c`.
Parameters
----------
total_mass: array_like
Total halo mass in :math:`M_{\odot}/h`; can be a number or a numpy array.
conc : array_like
Value of the halo concentration. Can either be a scalar, or a numpy array
of the same dimension as the input ``total_mass``.
Returns
--------
rmax : array_like
:math:`R_{\rm max}` in Mpc/h.
Examples
--------
>>> model = NFWProfile()
>>> Npts = 100
>>> total_mass = np.zeros(Npts) + 1e12
>>> conc = 5
>>> result = model.rmax(total_mass, conc)
>>> concarr = np.linspace(1, 100, Npts)
>>> result = model.rmax(total_mass, concarr)
Notes
------
See :ref:`halo_profile_definitions` for derivations and implementation details.
"""
halo_radius = self.halo_mass_to_halo_radius(total_mass)
scale_radius = halo_radius / conc
return 2.16258 * scale_radius
[docs]
def vmax(self, total_mass, conc):
r"""Maximum circular velocity of the halo profile,
:math:`V_{\rm max}^{\rm NFW} = V_{\rm cir}^{\rm NFW}(r = 2.16258R_{\Delta}/c)`.
Parameters
----------
total_mass: array_like
Total halo mass in :math:`M_{\odot}/h`; can be a number or a numpy array.
conc : array_like
Value of the halo concentration. Can either be a scalar, or a numpy array
of the same dimension as the input ``total_mass``.
Returns
--------
vmax : array_like
:math:`V_{\rm max}` in km/s.
Examples
--------
>>> model = NFWProfile()
>>> Npts = 100
>>> total_mass = np.zeros(Npts) + 1e12
>>> conc = 5
>>> result = model.vmax(total_mass, conc)
>>> concarr = np.linspace(1, 100, Npts)
>>> result = model.vmax(total_mass, concarr)
Notes
------
See :ref:`halo_profile_definitions` for derivations and implementation details.
"""
Rmax = self.rmax(total_mass, conc)
vmax = self.circular_velocity(Rmax, total_mass, conc)
return vmax
[docs]
def halo_mass_to_halo_radius(self, total_mass):
r"""
Spherical overdensity radius as a function of the input mass.
Note that this function is independent of the form of the density profile.
Parameters
----------
total_mass: array_like
Total halo mass in :math:`M_{\odot}/h`; can be a number or a numpy array.
Returns
-------
radius : array_like
Radius of the halo in Mpc/h units.
Will have the same dimension as the input ``total_mass``.
Examples
--------
>>> model = NFWProfile()
>>> halo_radius = model.halo_mass_to_halo_radius(1e13)
"""
return AnalyticDensityProf.halo_mass_to_halo_radius(self, total_mass)
[docs]
def halo_radius_to_halo_mass(self, radius):
r"""
Spherical overdensity mass as a function of the input radius.
Note that this function is independent of the form of the density profile.
Parameters
------------
radius : array_like
Radius of the halo in Mpc/h units; can be a number or a numpy array.
Returns
----------
total_mass: array_like
Total halo mass in :math:`M_{\odot}/h`.
Will have the same dimension as the input ``radius``.
Examples
--------
>>> model = NFWProfile()
>>> halo_mass = model.halo_mass_to_halo_radius(500.)
"""
return AnalyticDensityProf.halo_radius_to_halo_mass(self, radius)
[docs]
def mc_generate_nfw_radial_positions(self, **kwargs):
r"""Return a Monte Carlo realization of points in an NFW profile.
See :ref:`monte_carlo_nfw_spatial_profile` for a discussion of this technique.
Parameters
-----------
num_pts : int, optional
Number of points in the Monte Carlo realization of the profile.
Default is 1e4.
conc : float, optional
Concentration of the NFW profile being realized.
Default is 5.
halo_mass : float, optional
Total mass of the halo whose profile is being realized,
used to define the halo boundary for the mass definition
bound to the NFWProfile instance as ``mdef``.
If ``halo_mass`` is unspecified,
keyword argument ``halo_radius`` must be specified.
halo_radius : float, optional
Physical boundary of the halo whose profile is being realized
in units of Mpc/h.
If ``halo_radius`` is unspecified,
keyword argument ``halo_mass`` must be specified, in which case the
outer boundary of the halo will be determined according to the mass definition
bound to the NFWProfile instance as ``mdef``.
seed : int, optional
Random number seed used in the Monte Carlo realization.
Default is None, which will produce stochastic results.
Returns
--------
radial_positions : array_like
Numpy array storing a Monte Carlo realization of the halo profile.
All values will lie strictly between 0 and the halo boundary.
Examples
---------
>>> nfw = NFWProfile()
>>> radial_positions = nfw.mc_generate_nfw_radial_positions(halo_mass = 1e12, conc = 10)
>>> radial_positions = nfw.mc_generate_nfw_radial_positions(halo_radius = 0.25)
"""
kwargs["mdef"] = self.mdef
kwargs["cosmology"] = self.cosmology
kwargs["redshift"] = self.redshift
return standalone_mc_generate_nfw_radial_positions(**kwargs)