Source code for galpy.df.osipkovmerrittNFWdf
# Class that implements the anisotropic spherical NFW DF with radially
# varying anisotropy of the Osipkov-Merritt type
import numpy
from ..util import conversion
from ..potential import NFWPotential
from .osipkovmerrittdf import _osipkovmerrittdf
from .isotropicNFWdf import isotropicNFWdf
_COEFFS= numpy.array([-0.9958957901383353, 4.2905266124525259,
-7.6069046709185919,7.0313234865878806,
-3.6920719890718035, 0.8313023634615980,
-0.2179687331774083, -0.0408426627412238,
0.0802975743915827])
[docs]class osipkovmerrittNFWdf(_osipkovmerrittdf):
"""Class that implements the anisotropic spherical NFW DF with radially varying anisotropy of the Osipkov-Merritt type
.. math::
\\beta(r) = \\frac{1}{1+r_a^2/r^2}
with :math:`r_a` the anistropy radius.
"""
[docs] def __init__(self,pot=None,ra=1.4,rmax=1e4,ro=None,vo=None):
"""
NAME:
__init__
PURPOSE:
Initialize a NFW DF with Osipkov-Merritt anisotropy
INPUT:
pot - NFW potential which determines the DF
ra - anisotropy radius (can be a Quantity)
rmax= (1e4) maximum radius to consider (can be Quantity); set to numpy.inf to evaluate NFW w/o cut-off
ro=, vo= galpy unit parameters
OUTPUT:
None
HISTORY:
2020-11-12 - Written - Bovy (UofT)
"""
assert isinstance(pot,NFWPotential), \
'pot= must be potential.NFWPotential'
_osipkovmerrittdf.__init__(self,pot=pot,ra=ra,rmax=rmax,ro=ro,vo=vo)
self._Qtildemax= pot._amp/pot.a
self._Qtildemin= -pot(self._rmax,0,use_physical=False)/self._Qtildemax
self._a2overra2= self._pot.a**2./self._ra2
self._fQnorm= self._a2overra2/(4.*numpy.pi)/pot.a**1.5/pot._amp**0.5
# Initialize isotropic version to use as part of fQ
self._idf= isotropicNFWdf(pot=pot,rmax=rmax,ro=ro,vo=vo)
def fQ(self,Q):
"""
NAME:
fQ
PURPOSE
Calculate the f(Q) portion of an Osipkov-Merritt NFW distribution function
INPUT:
Q - The Osipkov-Merritt 'energy' E-L^2/[2ra^2] (can be Quantity)
OUTPUT:
fQ - The value of the f(Q) portion of the DF
HISTORY:
2021-02-09 - Written - Bovy (UofT)
"""
Qtilde= conversion.parse_energy(Q,vo=self._vo)/self._Qtildemax
out= numpy.zeros_like(Qtilde)
indx= (Qtilde > self._Qtildemin)*(Qtilde <= 1.)
# The 'ergodic' part
out[indx]= self._idf.fE(-Q[indx])
# The other part
out[indx]+= self._fQnorm*numpy.polyval(_COEFFS,Qtilde[indx])\
/(Qtilde[indx]**(2./3.)\
*(numpy.log(Qtilde[indx])/(1-Qtilde[indx]))**2.)
return out