Source code for galpy.potential.BurkertPotential
###############################################################################
# BurkertPotential.py: Potential with a Burkert density
###############################################################################
import numpy
from scipy import special
from .SphericalPotential import SphericalPotential
from ..util import conversion
[docs]class BurkertPotential(SphericalPotential):
"""BurkertPotential.py: Potential with a Burkert density
.. math::
\\rho(r) = \\frac{\\mathrm{amp}}{(1+r/a)\\,(1+[r/a]^2)}
"""
[docs] def __init__(self,amp=1.,a=2.,normalize=False,
ro=None,vo=None):
"""
NAME:
__init__
PURPOSE:
initialize a Burkert-density potential
INPUT:
amp - amplitude to be applied to the potential (default: 1); can be a Quantity with units of mass density or Gxmass density
a = scale radius (can be Quantity)
normalize - if True, normalize such that vc(1.,0.)=1., or, if
given as a number, such that the force is this fraction
of the force necessary to make vc(1.,0.)=1.
ro=, vo= distance and velocity scales for translation into internal units (default from configuration file)
OUTPUT:
(none)
HISTORY:
2013-04-10 - Written - Bovy (IAS)
2020-03-30 - Re-implemented using SphericalPotential - Bovy (UofT)
"""
SphericalPotential.__init__(self,amp=amp,ro=ro,vo=vo,
amp_units='density')
a= conversion.parse_length(a,ro=self._ro,vo=self._vo)
self.a=a
self._scale= self.a
if normalize or \
(isinstance(normalize,(int,float)) \
and not isinstance(normalize,bool)): #pragma: no cover
self.normalize(normalize)
self.hasC= True
self.hasC_dxdv= True
self.hasC_dens= True
return None
def _revaluate(self,r,t=0.):
"""Potential as a function of r and time"""
x= r/self.a
return -self.a**2.*numpy.pi*(-numpy.pi/x+2.*(1./x+1)*numpy.arctan(1/x)
+(1./x+1)*numpy.log((1.+1./x)**2./(1.+1/x**2.))
+special.xlogy(2./x,1.+x**2.))
#Previous way, not stable as r -> infty
#return -self.a**2.*numpy.pi/x*(-numpy.pi+2.*(1.+x)*numpy.arctan(1/x)
# +2.*(1.+x)*numpy.log(1.+x)
# +(1.-x)*numpy.log(1.+x**2.))
def _rforce(self,r,t=0.):
x= r/self.a
return self.a*numpy.pi/x**2.*(numpy.pi-2.*numpy.arctan(1./x)
-2.*numpy.log(1.+x)-numpy.log(1.+x**2.))
def _r2deriv(self,r,t=0.):
x= r/self.a
return 4.*numpy.pi/(1.+x**2.)/(1.+x)+2.*self._rforce(r)/x/self.a
def _rdens(self,r,t=0.):
x= r/self.a
return 1./(1.+x)/(1.+x**2.)
def _surfdens(self,R,z,phi=0.,t=0.):
"""
NAME:
_surfdens
PURPOSE:
evaluate the surface density for this potential
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
OUTPUT:
the surface density
HISTORY:
2018-08-19 - Written - Bovy (UofT)
"""
r= numpy.sqrt(R**2.+z**2.)
x= r/self.a
Rpa= numpy.sqrt(R**2.+self.a**2.)
Rma= numpy.sqrt(R**2.-self.a**2.+0j)
if Rma == 0:
za= z/self.a
return self.a**2./2.*((2.-2.*numpy.sqrt(za**2.+1)
+numpy.sqrt(2.)*za\
*numpy.arctan(za/numpy.sqrt(2.)))/z
+numpy.sqrt(2*za**2.+2.)\
*numpy.arctanh(za/numpy.sqrt(2.*(za**2.+1)))
/numpy.sqrt(self.a**2.+z**2.))
else:
return self.a**2.*(numpy.arctan(z/x/Rma)/Rma
+numpy.arctanh(z/x/Rpa)/Rpa
-numpy.arctan(z/Rma)/Rma
+numpy.arctan(z/Rpa)/Rpa).real