Source code for galpy.df.constantbetaHernquistdf
# Class that implements the anisotropic spherical Hernquist DF with constant
# beta parameter
import numpy
import scipy.special
import scipy.integrate
from ..util import conversion
from ..potential import evaluatePotentials,HernquistPotential
from .constantbetadf import _constantbetadf
[docs]class constantbetaHernquistdf(_constantbetadf):
"""Class that implements the anisotropic spherical Hernquist DF with constant beta parameter"""
[docs] def __init__(self,pot=None,beta=0,ro=None,vo=None):
"""
NAME:
__init__
PURPOSE:
Initialize a Hernquist DF with constant anisotropy
INPUT:
pot - Hernquist potential which determines the DF
beta - anisotropy parameter
OUTPUT:
None
HISTORY:
2020-07-22 - Written - Lane (UofT)
"""
assert isinstance(pot,HernquistPotential),'pot= must be potential.HernquistPotential'
_constantbetadf.__init__(self,pot=pot,beta=beta,ro=ro,vo=vo)
self._psi0= -evaluatePotentials(self._pot,0,0,use_physical=False)
self._potInf= 0.
self._GMa = self._psi0*self._pot.a**2.
# Final factor is mass to make the DF that of the mass density
self._fEnorm= (2.**self._beta/(2.*numpy.pi)**2.5)\
*scipy.special.gamma(5.-2.*self._beta)\
/scipy.special.gamma(1.-self._beta)\
/scipy.special.gamma(3.5-self._beta)\
/self._GMa**(1.5-self._beta)\
*self._psi0*self._pot.a
def fE(self,E):
"""
NAME:
fE
PURPOSE
Calculate the energy portion of a Hernquist distribution function
INPUT:
E - The energy (can be Quantity)
OUTPUT:
fE - The value of the energy portion of the DF
HISTORY:
2020-07-22 - Written
"""
Etilde= -conversion.parse_energy(E,vo=self._vo)/self._psi0
# Handle potential E outside of bounds
Etilde_out = numpy.where(numpy.logical_or(Etilde<0,Etilde>1))[0]
if len(Etilde_out)>0:
# Dummy variable now and 0 later, prevents numerical issues?
Etilde[Etilde_out]=0.5
# First check algebraic solutions, all adjusted such that DF = mass den
if self._beta == 0.: # isotropic case
sqrtEtilde= numpy.sqrt(Etilde)
fE= self._psi0*self._pot.a\
/numpy.sqrt(2.)/(2*numpy.pi)**3/self._GMa**1.5\
*sqrtEtilde/(1-Etilde)**2.\
*((1.-2.*Etilde)*(8.*Etilde**2.-8.*Etilde-3.)\
+((3.*numpy.arcsin(sqrtEtilde))\
/numpy.sqrt(Etilde*(1.-Etilde))))
elif self._beta == 0.5:
fE= (3.*Etilde**2.)/(4.*numpy.pi**3.*self._pot.a)
elif self._beta == -0.5:
fE= ((20.*Etilde**3.-20.*Etilde**4.+6.*Etilde**5.)\
/(1.-Etilde)**4)/(4.*numpy.pi**3.*self._GMa*self._pot.a)
else:
fE= self._fEnorm*numpy.power(Etilde,2.5-self._beta)*\
scipy.special.hyp2f1(5.-2.*self._beta,1.-2.*self._beta,
3.5-self._beta,Etilde)
if len(Etilde_out) > 0:
fE[Etilde_out]= 0.
return fE
def _icmf(self,ms):
'''Analytic expression for the normalized inverse cumulative mass
function. The argument ms is normalized mass fraction [0,1]'''
return self._pot.a*numpy.sqrt(ms)/(1-numpy.sqrt(ms))