###############################################################################
# ChandrasekharDynamicalFrictionForce: Class that implements the
# Chandrasekhar dynamical friction
###############################################################################
import copy
import hashlib
import numpy
from scipy import special, interpolate
from ..util import conversion
from .DissipativeForce import DissipativeForce
from .Potential import evaluateDensities, _check_c
from .Potential import flatten as flatten_pot
_INVSQRTTWO= 1./numpy.sqrt(2.)
_INVSQRTPI= 1./numpy.sqrt(numpy.pi)
[docs]class ChandrasekharDynamicalFrictionForce(DissipativeForce):
"""Class that implements the Chandrasekhar dynamical friction force
.. math::
\\mathbf{F}(\\mathbf{x},\\mathbf{v}) = -2\\pi\\,[G\\,M]\\,[G\\,\\rho(\\mathbf{x})]\\,\\ln[1+\\Lambda^2] \\,\\left[\\mathrm{erf}(X)-\\frac{2X}{\\sqrt{\\pi}}\\exp\\left(-X^2\\right)\\right]\\,\\frac{\\mathbf{v}}{|\\mathbf{v}|^3}\\,
on a mass (e.g., a satellite galaxy or a black hole) :math:`M` at position :math:`\\mathbf{x}` moving at velocity :math:`\\mathbf{v}` through a background density :math:`\\rho`. The quantity :math:`X` is the usual :math:`X=|\\mathbf{v}|/[\\sqrt{2}\\sigma_r(r)`. The factor :math:`\\Lambda` that goes into the Coulomb logarithm is taken to be
.. math::
\\Lambda = \\frac{r/\\gamma}{\\mathrm{max}\\left(r_{\\mathrm{hm}},GM/|\\mathbf{v}|^2\\right)}\\,,
where :math:`\\gamma` is a constant. This :math:`\\gamma` should be the absolute value of the logarithmic slope of the density :math:`\\gamma = |\\mathrm{d} \\ln \\rho / \\mathrm{d} \\ln r|`, although for :math:`\\gamma<1` it is advisable to set :math:`\\gamma=1`. Implementation here roughly follows `2016MNRAS.463..858P <http://adsabs.harvard.edu/abs/2016MNRAS.463..858P>`__ and earlier work.
"""
[docs] def __init__(self,amp=1.,GMs=.1,gamma=1.,rhm=0.,
dens=None,sigmar=None,
const_lnLambda=False,minr=0.0001,maxr=25.,nr=501,
ro=None,vo=None):
"""
NAME:
__init__
PURPOSE:
initialize a Chandrasekhar Dynamical Friction force
INPUT:
amp - amplitude to be applied to the potential (default: 1)
GMs - satellite mass; can be a Quantity with units of mass or Gxmass; can be adjusted after initialization by setting obj.GMs= where obj is your ChandrasekharDynamicalFrictionForce instance (note that the mass of the satellite can *not* be changed simply by multiplying the instance by a number, because he mass is not only used as an amplitude)
rhm - half-mass radius of the satellite (set to zero for a black hole; can be a Quantity); can be adjusted after initialization by setting obj.rhm= where obj is your ChandrasekharDynamicalFrictionForce instance
gamma - Free-parameter in :math:`\\Lambda`
dens - Potential instance or list thereof that represents the density [default: LogarithmicHaloPotential(normalize=1.,q=1.)]
sigmar= (None) function that gives the velocity dispersion as a function of r (has to be in natural units!); if None, computed from the dens potential using the spherical Jeans equation (in galpy.df.jeans) assuming zero anisotropy; if set to a lambda function, *the object cannot be pickled* (so set it to a real function)
cont_lnLambda= (False) if set to a number, use a constant ln(Lambda) instead with this value
minr= (0.0001) minimum r at which to apply dynamical friction: at r < minr, friction is set to zero (can be a Quantity)
Interpolation:
maxr= (25) maximum r for which sigmar gets interpolated; for best performance set this to the maximum r you will consider (can be a Quantity)
nr= (501) number of radii to use in the interpolation of sigmar
You can check that sigmar is interpolated correctly by comparing the methods sigmar [the interpolated version] and sigmar_orig [the original or directly computed version]
OUTPUT:
(none)
HISTORY:
2011-12-26 - Started - Bovy (NYU)
2018-03-18 - Re-started: updated to r dependent Lambda form and integrated into galpy framework - Bovy (UofT)
2018-07-23 - Calculate sigmar from the Jeans equation and interpolate it; allow GMs and rhm to be set on the fly - Bovy (UofT)
"""
DissipativeForce.__init__(self,amp=amp*GMs,ro=ro,vo=vo,
amp_units='mass')
rhm= conversion.parse_length(rhm,ro=self._ro)
minr= conversion.parse_length(minr,ro=self._ro)
maxr= conversion.parse_length(maxr,ro=self._ro)
self._gamma= gamma
self._ms= self._amp/amp # from handling in __init__ above, should be ms in galpy units
self._rhm= rhm
self._minr= minr
self._maxr= maxr
self._dens_kwarg= dens # for pickling
self._sigmar_kwarg= sigmar # for pickling
# Parse density
if dens is None:
from .LogarithmicHaloPotential import LogarithmicHaloPotential
dens= LogarithmicHaloPotential(normalize=1.,q=1.)
if sigmar is None: # we know this solution!
sigmar= lambda x: _INVSQRTTWO
dens= flatten_pot(dens)
self._dens_pot= dens
self._dens=\
lambda R,z,phi=0.,t=0.: evaluateDensities(self._dens_pot,
R,z,phi=phi,t=t,
use_physical=False)
if sigmar is None:
from ..df import jeans
sigmar= lambda x: jeans.sigmar(self._dens_pot,x,beta=0.,
use_physical=False)
self._sigmar_rs_4interp= numpy.linspace(self._minr,self._maxr,nr)
self._sigmars_4interp= numpy.array([sigmar(x)
for x in self._sigmar_rs_4interp])
if numpy.any(numpy.isnan(self._sigmars_4interp)):
# Check for case where density is zero, in that case, just
# paint in the nearest neighbor for the interpolation
# (doesn't matter in the end, because force = 0 when dens = 0)
nanrs_indx= numpy.isnan(self._sigmars_4interp)
if numpy.all(numpy.array([self._dens(r*_INVSQRTTWO,r*_INVSQRTTWO)
for r in
self._sigmar_rs_4interp[nanrs_indx]])
== 0.):
self._sigmars_4interp[nanrs_indx]= interpolate.interp1d(\
self._sigmar_rs_4interp[True^nanrs_indx],
self._sigmars_4interp[True^nanrs_indx],
kind="nearest",fill_value="extrapolate")\
(self._sigmar_rs_4interp[nanrs_indx])
self.sigmar_orig= sigmar
self.sigmar= interpolate.InterpolatedUnivariateSpline(\
self._sigmar_rs_4interp,self._sigmars_4interp,k=3)
if const_lnLambda:
self._lnLambda= const_lnLambda
else:
self._lnLambda= False
self._amp*= 4.*numpy.pi
self._force_hash= None
self.hasC= _check_c(self._dens_pot,dens=True)
return None
def GMs(self,gms):
gms= conversion.parse_mass(gms,ro=self._ro,vo=self._vo)
self._amp*= gms/self._ms
self._ms= gms
# Reset the hash
self._force_hash= None
return None
GMs= property(None,GMs)
def rhm(self,new_rhm):
self._rhm= conversion.parse_length(new_rhm,ro=self._ro)
# Reset the hash
self._force_hash= None
return None
rhm= property(None,rhm)
def lnLambda(self,r,v):
"""
NAME:
lnLambda
PURPOSE:
evaluate the Coulomb logarithm :math:`\\ln \\Lambda`
INPUT:
r - spherical radius (natural units)
v - current velocity in cylindrical coordinates (natural units)
OUTPUT:
Coulomb logarithm
HISTORY:
2018-03-18 - Started - Bovy (UofT)
"""
if self._lnLambda:
lnLambda= self._lnLambda
else:
GMvs= self._ms/v**2.
if GMvs < self._rhm:
Lambda= r/self._gamma/self._rhm
else:
Lambda= r/self._gamma/GMvs
lnLambda= 0.5*numpy.log(1.+Lambda**2.)
return lnLambda
def _calc_force(self,R,phi,z,v,t):
r= numpy.sqrt(R**2.+z**2.)
if r < self._minr:
self._cached_force= 0.
else:
vs= numpy.sqrt(v[0]**2.+v[1]**2.+v[2]**2.)
if r > self._maxr:
sr= self.sigmar_orig(r)
else:
sr= self.sigmar(r)
X= vs*_INVSQRTTWO/sr
Xfactor= special.erf(X)-2.*X*_INVSQRTPI*numpy.exp(-X**2.)
lnLambda= self.lnLambda(r,vs)
self._cached_force=\
-self._dens(R,z,phi=phi,t=t)/vs**3.\
*Xfactor*lnLambda
def _Rforce(self,R,z,phi=0.,t=0.,v=None):
"""
NAME:
_Rforce
PURPOSE:
evaluate the radial force for this potential
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
v= current velocity in cylindrical coordinates
OUTPUT:
the radial force
HISTORY:
2018-03-18 - Started - Bovy (UofT)
"""
new_hash= hashlib.md5(numpy.array([R,phi,z,v[0],v[1],v[2],t]))\
.hexdigest()
if new_hash != self._force_hash:
self._calc_force(R,phi,z,v,t)
return self._cached_force*v[0]
def _phiforce(self,R,z,phi=0.,t=0.,v=None):
"""
NAME:
_phiforce
PURPOSE:
evaluate the azimuthal force for this potential
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
v= current velocity in cylindrical coordinates
OUTPUT:
the azimuthal force
HISTORY:
2018-03-18 - Started - Bovy (UofT)
"""
new_hash= hashlib.md5(numpy.array([R,phi,z,v[0],v[1],v[2],t]))\
.hexdigest()
if new_hash != self._force_hash:
self._calc_force(R,phi,z,v,t)
return self._cached_force*v[1]*R
def _zforce(self,R,z,phi=0.,t=0.,v=None):
"""
NAME:
_zforce
PURPOSE:
evaluate the vertical force for this potential
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
v= current velocity in cylindrical coordinates
OUTPUT:
the vertical force
HISTORY:
2018-03-18 - Started - Bovy (UofT)
"""
new_hash= hashlib.md5(numpy.array([R,phi,z,v[0],v[1],v[2],t]))\
.hexdigest()
if new_hash != self._force_hash:
self._calc_force(R,phi,z,v,t)
return self._cached_force*v[2]
# Pickling functions
def __getstate__(self):
pdict= copy.copy(self.__dict__)
# rm lambda function
del pdict['_dens']
if self._sigmar_kwarg is None:
# because an object set up with sigmar = user-provided function
# cannot typically be picked, disallow this explicitly
# (so if it can, everything should be fine; if not, pickling error)
del pdict['sigmar_orig']
return pdict
def __setstate__(self,pdict):
self.__dict__= pdict
# Re-setup _dens
self._dens=\
lambda R,z,phi=0.,t=0.: evaluateDensities(self._dens_pot,
R,z,phi=phi,t=t,
use_physical=False)
# Re-setup sigmar_orig
if self._dens_kwarg is None and self._sigmar_kwarg is None:
self.sigmar_orig= lambda x: _INVSQRTTWO
else:
from ..df import jeans
self.sigmar_orig= lambda x: jeans.sigmar(self._dens_pot,x,beta=0.,
use_physical=False)
return None
