# coding=utf8
# Copyright (c) MPE/IR-Submm Group. See LICENSE.rst for license information.
#
# Halo mass models for DysmalPy
from __future__ import (absolute_import, division, print_function,
unicode_literals)
# Standard library
import abc
import logging
import copy
# Third party imports
import numpy as np
import scipy.special as scp_spec
import scipy.optimize as scp_opt
import astropy.constants as apy_con
import astropy.units as u
import astropy.cosmology as apy_cosmo
# Local imports
from .model_set import ModelSet
from .base import MassModel
from dysmalpy.parameters import DysmalParameter
__all__ = ['NFW', 'TwoPowerHalo', 'Burkert', 'Einasto', 'DekelZhao', 'LinearNFW']
# DEFAULT COSMOLOGY
_default_cosmo = apy_cosmo.FlatLambdaCDM(H0=70., Om0=0.3)
# CONSTANTS
G = apy_con.G
Msun = apy_con.M_sun
pc = apy_con.pc
g_pc_per_Msun_kmssq = G.to(u.pc / u.Msun * (u.km / u.s) ** 2).value
# # +++++++++++++++++++++++++++++
# # TEMP:
# G = 6.67e-11 * u.m**3 / u.kg / (u.s**2) #(unit='m3 / (kg s2)')
# Msun = 2e30 * u.kg
# pc = 3e16 * u.m
# # +++++++++++++++++++++++++++++
# LOGGER SETTINGS
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger('DysmalPy')
logger.setLevel(logging.INFO)
import warnings
warnings.filterwarnings("ignore")
_dict_lmvir_fac_test_z = {'zarr': np.arange(0.5, 2.75, 0.25),
'fdmarr': np.array([1.e-3, 1.e-2, 0.1, 0.25, 0.5, 0.75, 0.9, 1-1.e-2, 1-1.e-3]),
'facarr': np.array([[ 1.87027726e-01, 6.83525784e-01, 1.70740388e+00,
2.43045482e+00, 3.27806760e+00, 4.19108979e+00,
5.13078742e+00, 7.20681515e+00, 9.21406239e+00],
[ 1.42061858e-01, 6.00007080e-01, 1.53906815e+00,
2.21927776e+00, 3.03914007e+00, 3.93936616e+00,
4.87420697e+00, 6.94791177e+00, 8.95494591e+00],
[ 1.00551748e-01, 5.24281952e-01, 1.38498363e+00,
2.02073465e+00, 2.80851657e+00, 3.69274241e+00,
4.62126978e+00, 6.69190395e+00, 8.69865476e+00],
[ 6.31519829e-02, 4.57268141e-01, 1.24826995e+00,
1.84043178e+00, 2.59324412e+00, 3.45844497e+00,
4.37910433e+00, 6.44582962e+00, 8.45221694e+00],
[ 2.97794509e-02, 3.98491637e-01, 1.12869340e+00,
1.67974956e+00, 2.39616708e+00, 3.23967562e+00,
4.15086381e+00, 6.21275987e+00, 8.21869377e+00],
[ 5.36103016e-05, 3.46978470e-01, 1.02457472e+00,
1.53784436e+00, 2.21772071e+00, 3.03735532e+00,
3.93747250e+00, 5.99354714e+00, 7.99892789e+00],
[-2.64841174e-02, 3.01676546e-01, 9.33823097e-01,
1.41290431e+00, 2.05708386e+00, 2.85121748e+00,
3.73871828e+00, 5.78791915e+00, 7.79263749e+00],
[-5.02759808e-02, 2.61622377e-01, 8.54407920e-01,
1.30283654e+00, 1.91285041e+00, 2.68042250e+00,
3.55386264e+00, 5.59509028e+00, 7.59902746e+00],
[-7.17141754e-02, 2.25990008e-01, 7.84535250e-01,
1.20560418e+00, 1.78340279e+00, 2.52388554e+00,
3.38196112e+00, 5.41408290e+00, 7.41711095e+00]])}
class DarkMatterHalo(MassModel):
r"""
Base model for dark matter halos
Parameters
----------
mvirial : float
Virial mass
conc : float
Concentration parameter
fdm : float
Dark matter fraction
"""
# Standard parameters for a dark matter halo profile
mvirial = DysmalParameter(default=1.0, bounds=(5, 20))
fdm = DysmalParameter(default=-99.9, fixed=True, bounds=(0,1))
_subtype = 'dark_matter'
def __init__(self, z=0, cosmo=_default_cosmo, **kwargs):
self._z = z
self._cosmo = cosmo
self._set_hz()
super(DarkMatterHalo, self).__init__(**kwargs)
def __setstate__(self, state):
super(DarkMatterHalo, self).__setstate__(state)
# quick test if necessary to migrate:
if '_z' in state.keys():
pass
else:
migrate_keys = ['z', 'cosmo']
for mkey in migrate_keys:
if (mkey in state.keys()) and ('_{}'.format(mkey) not in state.keys()):
self.__dict__['_{}'.format(mkey)] = state[mkey]
del self.__dict__[mkey]
if '_hz' not in state.keys():
self._set_hz()
@property
def z(self):
return self._z
@z.setter
def z(self, value):
if value < 0:
raise ValueError("Redshift can't be negative!")
self._z = value
# Reset hz:
self._set_hz()
@property
def cosmo(self):
return self._cosmo
@cosmo.setter
def cosmo(self, new_cosmo):
if not isinstance(new_cosmo, apy_cosmo.FLRW):
raise TypeError("Cosmology must be an astropy.cosmology.FLRW "
"instance.")
if new_cosmo is None:
self._cosmo = _default_cosmo
self._cosmo = new_cosmo
# Reset hz:
self._set_hz()
def _set_hz(self):
self._hz = self.cosmo.H(self.z).value
def calc_rvir(self):
r"""
Calculate the virial radius based on virial mass and redshift
Returns
-------
rvir : float
Virial radius
Notes
-----
Formula:
.. math::
M_{\rm vir} = 100 \frac{H(z)^2 R_{\rm vir}^3}{G}
This is based on Mo, Mao, & White (1998) [1]_ which defines the virial
radius as the radius where the mean mass density is :math:`200\rho_{\rm crit}`.
:math:`\rho_{\rm crit}` is the critical density for closure at redshift, :math:`z`.
"""
# hz = self.cosmo.H(self.z).value
# rvir = ((10 ** self.mvirial * (g_pc_per_Msun_kmssq * 1e-3) /
# (10 * hz * 1e-3) ** 2) ** (1. / 3.))
rvir = ((10 ** self.mvirial * (g_pc_per_Msun_kmssq * 1e-3) /
(10 * self._hz * 1e-3) ** 2) ** (1. / 3.))
return rvir
@abc.abstractmethod
def calc_rho0(self, *args, **kwargs):
"""
Method to calculate the scale density
"""
def velocity_profile(self, r, model):
"""
Calculate velocity profile, including any adiabatic contraction
"""
if model.kinematic_options.adiabatic_contract:
raise NotImplementedError("Adiabatic contraction not currently supported!")
else:
return self.circular_velocity(r)
def calc_mvirial_from_fdm(self, baryons, r_fdm, adiabatic_contract=False):
"""
Calculate virial mass given dark matter fraction and baryonic distribution
Parameters
----------
baryons : `~dysmalpy.models.MassModel` or dictionary
Model component representing the baryons (assumed to be light emitting),
or dictionary containing a list of the baryon components (baryons['components'])
and a list of whether the baryon components are light emitting or not (baryons['light'])
r_fdm : float
Radius at which the dark matter fraction is determined
Returns
-------
mvirial : float
Virial mass in logarithmic solar units
Notes
-----
This uses the current value of `fdm` together with
the input baryon distribution to calculate the inferred `mvirial`.
"""
if (self.fdm.value > self.bounds['fdm'][1]) | \
((self.fdm.value < self.bounds['fdm'][0])):
mvirial = np.NaN
elif (self.fdm.value == 1.):
mvirial = np.inf
elif (self.fdm.value == 0.):
mvirial = -np.inf #-5. # as a small but finite value
elif (self.fdm.value < 1.e-10):
mvirial = -np.inf
elif (r_fdm < 0.):
mvirial = np.NaN
else:
if isinstance(baryons, dict):
vsqr_bar_re = 0
bar_mtot_linear = 0
for bcmp in baryons['components']:
vsqr_bar_re += bcmp.vcirc_sq(r_fdm)
bar_mtot_linear += 10**bcmp.total_mass.value
bar_mtot = np.log10(bar_mtot_linear)
else:
vsqr_bar_re = baryons.vcirc_sq(r_fdm)
bar_mtot = baryons.total_mass.value
vsqr_dm_re_target = vsqr_bar_re / (1./self.fdm.value - 1)
if not np.isfinite(vsqr_dm_re_target):
mvirial = np.NaN
else:
#mtest = np.arange(-5, 50, 1.0)
short_mtest = False
try:
if ((bar_mtot >= 8.) & (bar_mtot <=13.)):
whminz = np.argmin(np.abs(_dict_lmvir_fac_test_z['zarr']-self.z))
whminfdm = np.argmin(np.abs(_dict_lmvir_fac_test_z['fdmarr']-self.fdm.value))
fac_lmvir = _dict_lmvir_fac_test_z['facarr'][whminz,whminfdm]
rough_mvir = fac_lmvir - np.log10(1./self.fdm.value-1)+np.log10(0.5)+bar_mtot
mtest = np.arange(rough_mvir-1., rough_mvir+1.5, 0.5)
mtest = np.append(-5., mtest)
mtest = np.append(mtest, 50.)
short_mtest = True
else:
mtest = np.arange(-5, 50, 1.0)
except:
mtest = np.arange(-5, 50, 1.0)
vtest = np.array([self._minfunc_vdm_mvir_from_fdm(m, vsqr_dm_re_target,
r_fdm, baryons, adiabatic_contract) for m in mtest])
try:
a = mtest[vtest < 0][-1]
b = mtest[vtest > 0][0]
# # TEST
# if adiabatic_contract:
# a_noAC = mtest[vtest_noAC < 0][-1]
# b_noAC = mtest[vtest_noAC > 0][0]
except:
print("adiabatic_contract={}".format(adiabatic_contract))
print("fdm={}".format(self.fdm.value))
print("r_fdm={}".format(r_fdm))
print(mtest, vtest)
raise ValueError
# ------------------------------------------------------------------
# Do a quick check to make sure it was inside the inner range:
# If not, redo calc with better range:
if short_mtest & ((a == mtest[0]) | (b == mtest[-1])):
mtest_orig = mtest[:]
if (a == mtest[0]):
mtest = np.arange(np.ceil(mtest_orig[1])-2., np.ceil(mtest_orig[1])+1., 1.0)
mtest = np.append(np.arange(np.ceil(mtest_orig[1])-22., np.ceil(mtest_orig[1])-2., 2.0), mtest)
mtest = np.append(-5., mtest)
elif (b == mtest[-1]):
mtest = np.arange(np.floor(mtest_orig[-2]), np.floor(mtest_orig[-2])+3., 1.0)
mtest = np.append(mtest, np.arange(np.floor(mtest_orig[-2])+3., np.floor(mtest_orig[-2])+23., 2.0))
mtest = np.append(mtest, 50.)
vtest = np.array([self._minfunc_vdm_mvir_from_fdm(m, vsqr_dm_re_target, r_fdm,
baryons, adiabatic_contract) for m in mtest])
try:
a = mtest[vtest < 0][-1]
b = mtest[vtest > 0][0]
except:
print("adiabatic_contract={}".format(adiabatic_contract))
print("fdm={}".format(self.fdm.value))
print("r_fdm={}".format(r_fdm))
print(mtest, vtest)
raise ValueError
# ------------------------------------------------------------------
# Run optimizer:
mvirial = scp_opt.brentq(self._minfunc_vdm_mvir_from_fdm, a, b, args=(vsqr_dm_re_target, r_fdm,
baryons, adiabatic_contract))
return mvirial
def _minfunc_vdm_mvir_from_fdm(self, mvirial, vsqtarget, r_fdm, bary, adiabatic_contract):
halotmp = self.copy()
halotmp.__setattr__('mvirial', mvirial)
if adiabatic_contract:
modtmp = ModelSet()
if isinstance(bary, dict):
for bcmp,b_light in zip(bary['components'], bary['light']):
modtmp.add_component(bcmp, light=b_light)
else:
modtmp.add_component(bary, light=True)
modtmp.add_component(halotmp)
modtmp.kinematic_options.adiabatic_contract = True
modtmp.kinematic_options.adiabatic_contract_modify_small_values = True
modtmp._update_tied_parameters()
vc_sq, vc_sq_dm = modtmp.vcirc_sq(r_fdm, compute_dm=True)
return vc_sq_dm - vsqtarget
else:
return halotmp.vcirc_sq(r_fdm)- vsqtarget
[docs]
class NFW(DarkMatterHalo):
r"""
Dark matter halo following an NFW profile
Parameters
----------
mvirial : float
Virial mass in logarithmic solar units
conc : float
Concentration parameter
fdm : float
Dark matter fraction
z : float
Redshift
cosmo : `~astropy.cosmology` object
The cosmology to use for modelling.
If this model component will be attached to a `~dysmalpy.galaxy.Galaxy` make sure
the respective cosmologies are the same. Default is
`~astropy.cosmology.FlatLambdaCDM` with H0=70., and Om0=0.3.
Notes
-----
Model formula:
The mass density follows Navarro, Frenk, & White (1995) [1]_:
.. math::
\rho = \frac{\rho_0}{(r/r_s)(1 + r/r_s)^2}
:math:`r_s` is the scale radius defined as :math:`r_{\rm vir}/c`.
:math:`\rho_0` is the normalization parameter.
References
----------
.. [1] https://ui.adsabs.harvard.edu/abs/1995MNRAS.275..720N/abstract
"""
mvirial = DysmalParameter(default=1.0, bounds=(5, 20))
conc = DysmalParameter(default=5.0, bounds=(2, 20))
fdm = DysmalParameter(default=-99.9, fixed=True, bounds=(0,1))
def __init__(self, z=0, cosmo=_default_cosmo, **kwargs):
super(NFW, self).__init__(z=z, cosmo=cosmo, **kwargs)
[docs]
def evaluate(self, r, mvirial, conc, fdm):
"""Mass density as a function of radius"""
rvirial = self.calc_rvir()
rho0 = self.calc_rho0(rvirial=rvirial)
rs = rvirial / self.conc
return rho0 / (r / rs * (1 + r / rs) ** 2)
[docs]
def enclosed_mass(self, r):
"""
Enclosed mass as a function of radius
Parameters
----------
r : float or array
Radius or radii in kpc
Returns
-------
menc : float or array
Enclosed mass in solar units
"""
rvirial = self.calc_rvir()
rho0 = self.calc_rho0(rvirial=rvirial)
rs = rvirial/self.conc
aa = 4.*np.pi*rho0*rvirial**3/self.conc**3
# For very small r, bb can be negative.
bb = np.abs(np.log((rs + r)/rs) - r/(rs + r))
return aa*bb
[docs]
def calc_rho0(self, rvirial=None):
r"""
Normalization of the density distribution
Returns
-------
rho0 : float
Mass density normalization in :math:`M_{\odot}/\rm{kpc}^3`
"""
if rvirial is None:
rvirial = self.calc_rvir()
aa = 10**self.mvirial/(4.*np.pi*rvirial**3)*self.conc**3
bb = 1./(np.log(1.+self.conc) - (self.conc/(1.+self.conc)))
return aa * bb
[docs]
class TwoPowerHalo(DarkMatterHalo):
r"""
Two power law density model for a dark matter halo
Parameters
----------
mvirial : float
Virial mass in logarithmic solar units
conc : float
Concentration parameter
alpha : float
Power law index at small radii
beta : float
Power law index at large radii
fdm : float
Dark matter fraction
z : float
Redshift
cosmo : `~astropy.cosmology` object
The cosmology to use for modelling.
If this model component will be attached to a `~dysmalpy.galaxy.Galaxy` make sure
the respective cosmologies are the same. Default is
`~astropy.cosmology.FlatLambdaCDM` with H0=70., and Om0=0.3.
Notes
-----
Model formula:
The mass density follows Equation 2.64 of Binney & Tremaine (2008) [1]_:
.. math::
\rho=\frac{\rho_0}{(r/r_s)^\alpha(1 + r/r_s)^{\beta - \alpha}}
:math:`r_s` is the scale radius and defined as :math:`r_\mathrm{vir}/c` where
:math:`r_\mathrm{vir}` is the virial radius and :math:`c` is the concentration
parameter. :math:`\rho_0` is the normalization parameter.
References
----------
.. [1] https://ui.adsabs.harvard.edu/abs/2008gady.book.....B/abstract
"""
# Powerlaw slopes for the density model
mvirial = DysmalParameter(default=1.0, bounds=(5, 20))
conc = DysmalParameter(default=5.0, bounds=(2, 20))
alpha = DysmalParameter(default=1.0)
beta = DysmalParameter(default=3.0)
fdm = DysmalParameter(default=-99.9, fixed=True, bounds=(0,1))
_subtype = 'dark_matter'
def __init__(self, z=0, cosmo=_default_cosmo, **kwargs):
super(TwoPowerHalo, self).__init__(z=z, cosmo=cosmo, **kwargs)
[docs]
def evaluate(self, r, mvirial, conc, alpha, beta, fdm):
""" Mass density for the TwoPowerHalo"""
rvirial = self.calc_rvir()
rho0 = self.calc_rho0()
rs = rvirial / self.conc
return rho0 / ((r/rs)**alpha * (1 + r/rs)**(beta - alpha))
[docs]
def enclosed_mass(self, r):
"""
Enclosed mass as a function of radius
Parameters
----------
r : float or array
Radius or radii in kpc
Returns
-------
menc : float or array
Enclosed mass in solar units
"""
rvirial = self.calc_rvir()
rs = rvirial/self.conc
aa = 10**self.mvirial*(r/rvirial)**(3 - self.alpha)
bb = (scp_spec.hyp2f1(3-self.alpha, self.beta-self.alpha, 4-self.alpha, -r/rs) /
scp_spec.hyp2f1(3 - self.alpha, self.beta - self.alpha, 4 - self.alpha, -self.conc))
return aa*bb
[docs]
def calc_rho0(self, rvirial=None):
r"""
Normalization of the density distribution
Returns
-------
rho0 : float
Mass density normalization in :math:`M_{\odot}/\rm{kpc}^3`
"""
if rvirial is None:
rvirial = self.calc_rvir()
rs = rvirial/self.conc
aa = -10**self.mvirial/(4*np.pi*self.conc**(3-self.alpha)*rs**3)
bb = (self.alpha - 3) / scp_spec.hyp2f1(3-self.alpha, self.beta-self.alpha, 4-self.alpha, -self.conc)
return aa*bb
[docs]
def calc_alpha_from_fdm(self, baryons, r_fdm):
"""
Calculate alpha given dark matter fraction and baryonic distribution
Parameters
----------
baryons : `~dysmalpy.models.MassModel` or dictionary
Model component representing the baryons (assumed to be light emitting),
or dictionary containing a list of the baryon components (baryons['components'])
and a list of whether the baryon components are light emitting or not (baryons['light'])
r_fdm : float
Radius at which the dark matter fraction is determined
Returns
-------
alpha : float
alpha value
Notes
-----
This uses the current values of `fdm`, `mvirial`, and `beta` together with
the input baryon distribution to calculate the necessary value of `alpha`.
"""
if (self.fdm.value > self.bounds['fdm'][1]) | \
((self.fdm.value < self.bounds['fdm'][0])):
alpha = np.NaN
else:
if isinstance(baryons, dict):
vsqr_bar_re = 0
for bcmp in baryons['components']:
vsqr_bar_re += bcmp.vcirc_sq(r_fdm)
else:
vsqr_bar_re = baryons.vcirc_sq(r_fdm)
vsqr_dm_re_target = vsqr_bar_re / (1./self.fdm - 1)
# alphtest = np.arange(-50, 50, 1.)
alphtest = np.arange(0., 5., 1.0)
alphtest = np.append(-5., alphtest)
alphtest = np.append(alphtest, 50.)
vtest = np.array([self._minfunc_vdm_alpha_from_fdm(alph, vsqr_dm_re_target, self.mvirial, self.conc,
self.beta, self.z, r_fdm) for alph in alphtest])
try:
a = alphtest[vtest < 0][-1]
try:
b = alphtest[vtest > 0][0]
except:
a = alphtest[-2] # Even if not perfect, force in case of no convergence...
b = alphtest[-1]
except:
a = alphtest[0] # Even if not perfect, force in case of no convergence...
b = alphtest[1]
alpha = scp_opt.brentq(self._minfunc_vdm_alpha_from_fdm, a, b, args=(vsqr_dm_re_target, self.mvirial, self.conc,
self.beta, self.z, r_fdm))
return alpha
def _minfunc_vdm_alpha_from_fdm(self, alpha, vsqtarget, mass, conc, beta, z, r_fdm):
halo = TwoPowerHalo(mvirial=mass, conc=conc, alpha=alpha, beta=beta, z=z)
return halo.vcirc_sq(r_fdm)- vsqtarget
[docs]
class Burkert(DarkMatterHalo):
r"""
Dark matter halo following a Burkert profile
Parameters
----------
mvirial : float
Virial mass in logarithmic solar units
rB : float
Size of the dark matter core in kpc
fdm : float
Dark matter fraction
z : float
Redshift
cosmo : `~astropy.cosmology` object
The cosmology to use for modelling.
If this model component will be attached to a `~dysmalpy.galaxy.Galaxy` make sure
the respective cosmologies are the same. Default is
`~astropy.cosmology.FlatLambdaCDM` with H0=70., and Om0=0.3.
Notes
-----
Model formula:
The mass density follows Burkert (1995) [1]_:
.. math::
\rho=\frac{\rho_0}{(1 + r/r_B)(1 + (r/r_B)^2)}
References
----------
.. [1] https://ui.adsabs.harvard.edu/abs/1995ApJ...447L..25B/abstract
"""
mvirial = DysmalParameter(default=1.0, bounds=(5, 20))
rB = DysmalParameter(default=1.0)
fdm = DysmalParameter(default=-99.9, fixed=True, bounds=(0,1))
_subtype = 'dark_matter'
def __init__(self, z=0, cosmo=_default_cosmo, **kwargs):
super(Burkert, self).__init__(z=z, cosmo=cosmo, **kwargs)
[docs]
def evaluate(self, r, mvirial, rB, fdm):
"""Mass density as a function of radius"""
rho0 = self.calc_rho0()
return rho0 / ((1 + r/rB) * (1 + (r/rB)**2))
[docs]
def I(self, r):
Ival = 0.25 * (np.log(r**2 + self.rB**2) + 2.*np.log(r + self.rB)
- 2.*np.arctan(r/self.rB) - 4.*np.log(self.rB))
return Ival
[docs]
def enclosed_mass(self, r):
"""
Enclosed mass as a function of radius
Parameters
----------
r : float or array
Radius or radii in kpc
Returns
-------
menc : float or array
Enclosed mass in solar units
"""
rvir = self.calc_rvir()
Irvir = self.I(rvir)
aa = 10**self.mvirial / Irvir
bb = self.I(r)
return aa*bb
[docs]
def calc_rho0(self):
r"""
Normalization of the density distribution
Returns
-------
rho0 : float
Mass density normalization in :math:`M_{\odot}/\rm{kpc}^3`
"""
rvir = self.calc_rvir()
Irvir = self.I(rvir)
aa = 10**self.mvirial / (4*np.pi* self.rB**3)
bb = 1./Irvir
return aa*bb
[docs]
def calc_conc(self):
"""
Calculate the concentration parameter
Returns
-------
conc : float
Concentration based on the core radius, `rB`.
"""
rvir = self.calc_rvir()
conc = rvir/self.rB
self.conc = conc
return conc
[docs]
def calc_rB_from_fdm(self, baryons, r_fdm):
"""
Calculate core radius given dark matter fraction and baryonic distribution
Parameters
----------
baryons : `~dysmalpy.models.MassModel` or dictionary
Model component representing the baryons (assumed to be light emitting),
or dictionary containing a list of the baryon components (baryons['components'])
and a list of whether the baryon components are light emitting or not (baryons['light'])
r_fdm : float
Radius at which the dark matter fraction is determined
Returns
-------
rB : float
Core radius in kpc
Notes
-----
This uses the current values of `fdm`, and `mvirial` together with
the input baryon distribution to calculate the necessary value of `rB`.
"""
if (self.fdm.value > self.bounds['fdm'][1]) | \
((self.fdm.value < self.bounds['fdm'][0])):
rB = np.NaN
else:
if isinstance(baryons, dict):
vsqr_bar_re = 0
for bcmp in baryons['components']:
vsqr_bar_re += bcmp.vcirc_sq(r_fdm)
else:
vsqr_bar_re = baryons.vcirc_sq(r_fdm)
vsqr_dm_re_target = vsqr_bar_re / (1./self.fdm - 1)
rBtest = np.arange(0., 250., 5.0)
vtest = np.array([self._minfunc_vdm_rB_from_fDM(rBt, vsqr_dm_re_target, self.mvirial, self.z, r_fdm) for rBt in rBtest])
try:
a = rBtest[vtest < 0][-1]
try:
b = rBtest[vtest > 0][0]
except:
a = rBtest[0] # Even if not perfect, force in case of no convergence...
b = rBtest[1]
except:
a = rBtest[-2] # Even if not perfect, force in case of no convergence...
b = rBtest[-1]
try:
rB = scp_opt.brentq(self._minfunc_vdm_rB_from_fDM, a, b, args=(vsqr_dm_re_target, self.mvirial, self.z, r_fdm))
except:
# SOMETHING, if it's failing...
rB = np.average([a,b])
return rB
def _minfunc_vdm_rB_from_fDM(self, rB, vsqtarget, mass, z, r_fdm):
halo = Burkert(mvirial=mass, rB=rB, z=z)
return halo.vcirc_sq(r_fdm) - vsqtarget
[docs]
class Einasto(DarkMatterHalo):
r"""
Dark matter halo following an Einasto profile
Parameters
----------
mvirial : float
Virial mass in logarithmic solar units
conc : float
Concentration parameter
nEinasto : float
Inverse of the power law logarithmic slope
alphaEinasto : float
Power law logarithmic slope
fdm : float
Dark matter fraction
z : float
Redshift
cosmo : `~astropy.cosmology` object
The cosmology to use for modelling.
If this model component will be attached to a `~dysmalpy.galaxy.Galaxy` make sure
the respective cosmologies are the same. Default is
`~astropy.cosmology.FlatLambdaCDM` with H0=70., and Om0=0.3.
Einasto_param : {'None', 'nEinasto', 'alphaEinasto'}
Which parameter to leave as the free parameter. If 'None', the model
determines which parameter to use based on if `nEinasto` or `alphaEinasto`
is None. Default is 'None'
Notes
-----
Model formula following Retana-Montenegro et al (2012) [1]_:
.. math::
\rho = \rho_0 \exp\left\{-\left(\frac{r}{h}\right)^{1/n}\right\}
where :math:`h=r_s/(2n)^n` is the scale length and
:math:`r_s` is the scale radius defined as :math:`r_\mathrm{vir}/c`.
In this model only `nEinasto` or `alphaEinasto` can be free since :math:`n=1/\alpha`.
References
----------
.. [1] https://ui.adsabs.harvard.edu/abs/2012A%26A...540A..70R/abstract
"""
mvirial = DysmalParameter(default=1.0, bounds=(5, 20))
conc = DysmalParameter(default=5.0, bounds=(2, 20))
nEinasto = DysmalParameter(default=1.0)
alphaEinasto = DysmalParameter(default=-99., fixed=True)
fdm = DysmalParameter(default=-99.9, fixed=True, bounds=(0,1))
_subtype = 'dark_matter'
def __init__(self, z=0, cosmo=_default_cosmo,
Einasto_param='None', alphaEinasto=None, nEinasto=None, **kwargs):
# Check whether at least *one* of alphaEinasto and nEinasto is set:
if (alphaEinasto is None) & (nEinasto is None):
raise ValueError("Must set at least one of alphaEinasto and nEinasto!")
if (alphaEinasto is not None) & (nEinasto is not None) & (Einasto_param == 'None'):
raise ValueError("If both 'alphaEinasto' and 'nEinasto' are set, must specify which is the fit variable with 'Einasto_param'")
super(Einasto, self).__init__(z=z, cosmo=cosmo, **kwargs)
# Setup the "alternating" of whether to use nEinasto or alphaEinasto:
if (Einasto_param.lower() == 'neinasto') | (alphaEinasto is None):
self.Einasto_param = 'nEinasto'
self.alphaEinasto.fixed = False
self.alphaEinasto.tied = self.tie_alphaEinasto
elif (Einasto_param.lower() == 'alphaeinasto') | (nEinasto is None):
self.Einasto_param = 'alphaEinasto'
self.nEinasto.fixed = False
self.nEinasto.tied = self.tie_nEinasto
else:
raise ValueError("Einasto_param = {} not recognized! [options: 'nEinasto', 'alphaEinasto']".format(Einasto_param))
[docs]
def evaluate(self, r, mvirial, conc, alphaEinasto, nEinasto, fdm):
"""Mass density as a function of radius"""
if self.Einasto_param.lower() == 'alphaeinasto':
nEinasto = 1./alphaEinasto
rvirial = self.calc_rvir()
rho0 = self.calc_rho0(rvirial=rvirial)
rs = rvirial / conc
h = rs / np.power(2.*nEinasto, nEinasto)
# Return the density at a given radius:
return rho0 * np.exp(- np.power(r/h, 1./nEinasto))
# Equivalent to:
# rho0 * np.exp( - 2 * nEinasto * ( np.power(r/rs, 1./nEinasto) -1.) )
# or
# rho0 * np.exp( - 2 / alphaEinasto * ( np.power(r/rs, alphaEinasto) -1.) )
[docs]
def enclosed_mass(self, r):
"""
Enclosed mass as a function of radius
Parameters
----------
r : float or array
Radius or radii in kpc
Returns
-------
menc : float or array
Enclosed mass in solar units
"""
rvirial = self.calc_rvir()
rs = rvirial/self.conc
h = rs / np.power(2.*self.nEinasto, self.nEinasto)
rho0 = self.calc_rho0(rvirial=rvirial)
# Explicitly substituted for s = r/h before doing s^(1/nEinasto)
incomp_gam = scp_spec.gammainc(3*self.nEinasto, 2.*self.nEinasto * np.power(r/rs, 1./self.nEinasto) ) \
* scp_spec.gamma(3*self.nEinasto)
Menc = 4.*np.pi * rho0 * np.power(h, 3.) * self.nEinasto * incomp_gam
return Menc
[docs]
def calc_rho0(self, rvirial=None):
r"""
Density at the scale length
Returns
-------
rho0 : float
Mass density at the scale radius in :math:`M_{\odot}/\rm{kpc}^3`
"""
if rvirial is None:
rvirial = self.calc_rvir()
rs = rvirial/self.conc
h = rs / np.power(2.*self.nEinasto, self.nEinasto)
incomp_gam = scp_spec.gammainc(3*self.nEinasto, (2.*self.nEinasto) * \
np.power(self.conc, 1./self.nEinasto) ) \
* scp_spec.gamma(3*self.nEinasto)
rho0 = 10**self.mvirial / (4.*np.pi*self.nEinasto * np.power(h, 3.) * incomp_gam)
return rho0
[docs]
def calc_alphaEinasto_from_fdm(self, baryons, r_fdm):
"""
Calculate alpha given dark matter fraction and baryonic distribution
Parameters
----------
baryons : `~dysmalpy.models.MassModel` or dictionary
Model component representing the baryons (assumed to be light emitting),
or dictionary containing a list of the baryon components (baryons['components'])
and a list of whether the baryon components are light emitting or not (baryons['light'])
r_fdm : float
Radius at which the dark matter fraction is determined
Returns
-------
alphaEinasto : float
Power law logarithmic slope
Notes
-----
This uses the current values of `fdm`, and `mvirial` together with
the input baryon distribution to calculate the necessary value of `alphaEinasto`.
"""
nEinasto = self.calc_nEinasto_from_fdm(baryons, r_fdm)
if np.isfinite(nEinasto):
return 1./nEinasto
else:
return np.NaN
[docs]
def calc_nEinasto_from_fdm(self, baryons, r_fdm):
"""
Calculate n given the dark matter fraction and baryonic distribution
Parameters
----------
baryons : `~dysmalpy.models.MassModel` or dictionary
Model component representing the baryons (assumed to be light emitting),
or dictionary containing a list of the baryon components (baryons['components'])
and a list of whether the baryon components are light emitting or not (baryons['light'])
r_fdm : float
Radius at which the dark matter fraction is determined
Returns
-------
alphaEinasto : float
Power law logarithmic slope
Notes
-----
This uses the current values of `fdm`, and `mvirial` together with
the input baryon distribution to calculate the necessary value of `nEinasto`.
"""
if (self.fdm.value > self.bounds['fdm'][1]) | \
((self.fdm.value < self.bounds['fdm'][0])):
nEinasto = np.NaN
else:
# NOTE: have not tested this yet
if isinstance(baryons, dict):
vsqr_bar_re = 0
for bcmp in baryons['components']:
vsqr_bar_re += bcmp.vcirc_sq(r_fdm)
else:
vsqr_bar_re = baryons.vcirc_sq(r_fdm)
vsqr_dm_re_target = vsqr_bar_re / (1./self.fdm - 1)
nEinastotest = np.arange(-50, 50, 1.)
vtest = np.array([self._minfunc_vdm_nEin_from_fdm(nEinast, vsqr_dm_re_target, self.mvirial, self.conc,
self.alphaEinasto, self.z, r_fdm) for nEinast in nEinastotest])
try:
a = nEinastotest[vtest < 0][-1]
try:
b = nEinastotest[vtest > 0][0]
except:
a = nEinastotest[-2] # Even if not perfect, force in case of no convergence...
b = nEinastotest[-1]
except:
a = nEinastotest[0] # Even if not perfect, force in case of no convergence...
b = nEinastotest[1]
alpha = scp_opt.brentq(self._minfunc_vdm_nEin_from_fdm, a, b, args=(vsqr_dm_re_target,
self.mvirial, self.conc,
self.alphaEinasto, self.z, r_fdm))
return nEinasto
def _minfunc_vdm_nEin_from_fdm(self, nEinasto, vsqtarget, mass, conc, alphaEinasto, z, r_fdm):
halo = Einasto(mvirial=mass, conc=conc, nEinasto=nEinasto, alphaEinasto=alphaEinasto, z=z)
return halo.vcirc_sq(r_fdm) - vsqtarget
[docs]
def tie_nEinasto(self, model_set):
"""
Function to tie n to alpha
Parameters
----------
model_set : `ModelSet`
`ModelSet` the component is a part of and will be used in the fitting
Returns
-------
nEinasto : float
`nEinastro` given the current value of `alphaEinasto`
"""
if model_set.components['halo'].alphaEinasto.value != self.alphaEinasto:
raise ValueError
return 1./self.alphaEinasto
[docs]
def tie_alphaEinasto(self, model_set):
"""
Function to tie alpha to n
Parameters
----------
model_set : `ModelSet`
`ModelSet` the component is a part of and will be used in the fitting
Returns
-------
alphaEinasto : float
`alphaEinastro` given the current value of `nEinasto`
"""
return 1./self.nEinasto
[docs]
class DekelZhao(DarkMatterHalo):
r"""
Dekel-Zhao halo profile (Dekel et al. 2017, Freundlich et al. 2020)
Parameters
----------
mvirial : float
Virial mass in logarithmic solar units
s1 : float
Inner logarithmic slope (at resolution r1=0.01*Rvir)
c2 : float
Concentration parameter (defined relative to c, a)
fdm : float
Dark matter fraction
z : float
Redshift
cosmo : `~astropy.cosmology` object
The cosmology to use for modelling.
If this model component will be attached to a `~dysmalpy.galaxy.Galaxy` make sure
the respective cosmologies are the same. Default is
`~astropy.cosmology.FlatLambdaCDM` with H0=70., and Om0=0.3.
Notes
-----
The formula for this implementation are given in Freundlich et al. 2020. [1]_
Specifically, we use the forms where b=2, gbar=3 (see Eqns 4, 5, 14, 15)
References
----------
.. [1] https://ui.adsabs.harvard.edu/abs/2020MNRAS.499.2912F/abstract
"""
# Powerlaw slopes for the density model
mvirial = DysmalParameter(default=1.0, bounds=(5, 20))
s1 = DysmalParameter(default=1.5, bounds=(0.0, 2.0))
c2 = DysmalParameter(default=25., bounds=(0.0, 40.0))
fdm = DysmalParameter(default=-99.9, fixed=True, bounds=(0,1))
_subtype = 'dark_matter'
def __init__(self, z=0, cosmo=_default_cosmo, **kwargs):
super(DekelZhao, self).__init__(z=z, cosmo=cosmo, **kwargs)
[docs]
def evaluate(self, r, mvirial, s1, c2, fdm):
""" Mass density for the DekelZhao halo profile"""
rvirial = self.calc_rvir()
a, c = self.calc_a_c()
rhoc = self.calc_rho0(rvirial=rvirial, a=a, c=c)
rc = rvirial / c
x = r / rc
return rhoc / (np.power(x, a) * np.power((1.+np.sqrt(x)), 2.*(3.5-a)))
[docs]
def enclosed_mass(self, r):
"""
Enclosed mass as a function of radius
Parameters
----------
r : float or array
Radius or radii in kpc
Returns
-------
menc : float or array
Enclosed mass in solar units
"""
mvir = 10**self.mvirial
rvirial = self.calc_rvir()
a, c = self.calc_a_c()
rc = rvirial / c
x = r / rc
mu = self.calc_mu(a=a, c=c)
return mu * mvir / (np.power(x, a-3.)*np.power((1.+np.sqrt(x)), 2.*(3.-a)))
[docs]
def calc_a_c(self):
r"""
Calculate a, c from s1, c2 for the Dekel-Zhao halo.
Returns
-------
a, c: inner asymptotic slope, concentration parameter for DZ halo
"""
#rvirial = self.calc_rvir()
#r12 = np.sqrt(0.01*rvirial/rvirial)
r12 = np.sqrt(0.01)
c12 = np.sqrt(self.c2)
a = (1.5*self.s1 - 2.*(3.5-self.s1)*r12*c12)/(1.5 - (3.5-self.s1)*r12*c12)
c = ((self.s1-2.)/((3.5-self.s1)*r12 - 1.5/c12))**2
return a, c
[docs]
def calc_rho0(self, rvirial=None, a=None, c=None):
r"""
Normalization of the density distribution, rho_c
Returns
-------
rhoc : float
Mass density normalization in :math:`M_{\odot}/\rm{kpc}^3`
"""
if (a is None) or (c is None):
a, c = self.calc_a_c()
mu = self.calc_mu(a=a, c=c)
rhovirbar = self.calc_rhovirbar(rvirial=rvirial)
rhocbar = c**3 * mu * rhovirbar
return (1.-a/3.)*rhocbar
[docs]
def calc_rhovirbar(self, rvirial=None):
r"""
Average density in the virial radius, in :math:`M_{\odot}/\rm{kpc}^3`
"""
mvir = 10**self.mvirial
if rvirial is None:
rvirial = self.calc_rvir()
rhovirbar = (3.*mvir)/(4.*np.pi*(rvirial**3))
return rhovirbar
[docs]
def calc_mu(self, a=None, c=None):
"""
Subfunction for describing DZ profile
"""
if (a is None) or (c is None):
a, c = self.calc_a_c()
mu = np.power(c, a-3.) * np.power((1.+np.sqrt(c)), 2.*(3.-a))
return mu
def _minfunc_vdm_mvir_from_fdm(self, mvirial, vsqtarget, r_fdm, bary, adiabatic_contract):
halotmp = self.copy()
halotmp.__setattr__('mvirial', mvirial)
modtmp = ModelSet()
if isinstance(bary, dict):
for bcmp,b_light in zip(bary['components'], bary['light']):
modtmp.add_component(bcmp, light=b_light)
else:
modtmp.add_component(bary, light=True)
modtmp.add_component(halotmp)
modtmp.kinematic_options.adiabatic_contract = adiabatic_contract
if adiabatic_contract:
modtmp.kinematic_options.adiabatic_contract_modify_small_values = True
modtmp._update_tied_parameters()
if adiabatic_contract:
vc_sq, vc_sq_dm = modtmp.vcirc_sq(r_fdm, compute_dm=True)
return vc_sq_dm - vsqtarget
else:
return modtmp.components['halo'].vcirc_sq(r_fdm) - vsqtarget
[docs]
class LinearNFW(DarkMatterHalo):
r"""
Same as `NFW` except with the virial mass in linear units
Parameters
----------
mvirial : float
Virial mass in solar units
conc : float
Concentration parameter
fdm : float
Dark matter fraction
z : float
Redshift
cosmo : `~astropy.cosmology` object
The cosmology to use for modelling.
If this model component will be attached to a `~dysmalpy.galaxy.Galaxy` make sure
the respective cosmologies are the same. Default is
`~astropy.cosmology.FlatLambdaCDM` with H0=70., and Om0=0.3.
Notes
-----
Model formula:
The mass density follows Navarro, Frenk, & White (1995) [1]_:
.. math::
\rho = \frac{\rho_0}{(r/r_s)(1 + r/r_s)^2}
:math:`r_s` is the scale radius defined as :math:`r_{\rm vir}/c`.
:math:`\rho_0` is the normalization parameter.
References
----------
.. [1] https://ui.adsabs.harvard.edu/abs/1995MNRAS.275..720N/abstract
"""
mvirial = DysmalParameter(default=1.e1, bounds=(1.e5, 1.e20))
conc = DysmalParameter(default=5.0, bounds=(2, 20))
fdm = DysmalParameter(default=-99.9, fixed=True, bounds=(0,1))
def __init__(self, z=0, cosmo=_default_cosmo, **kwargs):
super(LinearNFW, self).__init__(z=z, cosmo=cosmo, **kwargs)
[docs]
def evaluate(self, r, mvirial, conc, fdm):
"""Mass density as a function of radius"""
rvirial = self.calc_rvir()
rho0 = self.calc_rho0(rvirial=rvirial)
rs = rvirial / self.conc
return rho0 / (r / rs * (1 + r / rs) ** 2)
[docs]
def enclosed_mass(self, r):
"""
Enclosed mass as a function of radius
Parameters
----------
r : float or array
Radius or radii in kpc
Returns
-------
menc : float or array
Enclosed mass in solar units
"""
rvirial = self.calc_rvir()
rho0 = self.calc_rho0(rvirial=rvirial)
rs = rvirial/self.conc
aa = 4.*np.pi*rho0*rvirial**3/self.conc**3
# For very small r, bb can be negative.
bb = np.abs(np.log((rs + r)/rs) - r/(rs + r))
return aa*bb
[docs]
def calc_rho0(self, rvirial=None):
r"""
Normalization of the density distribution
Returns
-------
rho0 : float
Mass density normalization in :math:`M_{\odot}/\rm{kpc}^3`
"""
if rvirial is None:
rvirial = self.calc_rvir()
aa = self.mvirial/(4.*np.pi*rvirial**3)*self.conc**3
bb = 1./(np.log(1.+self.conc) - (self.conc/(1.+self.conc)))
return aa * bb
[docs]
def calc_rvir(self):
r"""
Calculate the virial radius based on virial mass and redshift
Returns
-------
rvir : float
Virial radius
Notes
-----
Formula:
.. math::
M_{\rm vir} = 100 \frac{H(z)^2 R_{\rm vir}^3}{G}
This is based on Mo, Mao, & White (1998) [1]_ which defines the virial
radius as the radius where the mean mass density is :math:`200\rho_{\rm crit}`.
:math:`\rho_{\rm crit}` is the critical density for closure at redshift, :math:`z`.
"""
# hz = self.cosmo.H(self.z).value
# rvir = ((self.mvirial * (g_pc_per_Msun_kmssq * 1e-3) /
# (10 * hz * 1e-3) ** 2) ** (1. / 3.))
rvir = ((self.mvirial * (g_pc_per_Msun_kmssq * 1e-3) /
(10 * self._hz * 1e-3) ** 2) ** (1. / 3.))
return rvir