Source code for galpy.potential.PowerSphericalPotentialwCutoff

###############################################################################
#   PowerSphericalPotentialwCutoff.py: spherical power-law potential w/ cutoff
#
#                                     amp
#                          rho(r)= ---------   e^{-(r/rc)^2}
#                                   r^\alpha
###############################################################################
import numpy
from scipy import special
from ..util import conversion
from .Potential import Potential, kms_to_kpcGyrDecorator
[docs]class PowerSphericalPotentialwCutoff(Potential): """Class that implements spherical potentials that are derived from power-law density models .. math:: \\rho(r) = \\mathrm{amp}\,\\left(\\frac{r_1}{r}\\right)^\\alpha\\,\\exp\\left(-(r/rc)^2\\right) """
[docs] def __init__(self,amp=1.,alpha=1.,rc=1.,normalize=False,r1=1., ro=None,vo=None): """ NAME: __init__ PURPOSE: initialize a power-law-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 alpha= inner power rc= cut-off radius (can be Quantity) r1= (1.) reference radius for amplitude (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-06-28 - Written - Bovy (IAS) """ Potential.__init__(self,amp=amp,ro=ro,vo=vo,amp_units='density') r1= conversion.parse_length(r1,ro=self._ro) rc= conversion.parse_length(rc,ro=self._ro) self.alpha= alpha # Back to old definition self._amp*= r1**self.alpha self.rc= rc self._scale= self.rc 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 self._nemo_accname= 'PowSphwCut'
def _evaluate(self,R,z,phi=0.,t=0.): """ NAME: _evaluate PURPOSE: evaluate the potential at R,z INPUT: R - Galactocentric cylindrical radius z - vertical height phi - azimuth t - time OUTPUT: Phi(R,z) HISTORY: 2013-06-28 - Started - Bovy (IAS) """ r= numpy.sqrt(R**2.+z**2.) out= 2.*numpy.pi*self.rc**(3.-self.alpha)*(1/self.rc*special.gamma(1.-self.alpha/2.)*special.gammainc(1.-self.alpha/2.,(r/self.rc)**2.)-special.gamma(1.5-self.alpha/2.)*special.gammainc(1.5-self.alpha/2.,(r/self.rc)**2.)/r) if isinstance(r,(float,int)): if r == 0: return 0. else: return out else: out[r==0]= 0. return out def _Rforce(self,R,z,phi=0.,t=0.): """ NAME: _Rforce PURPOSE: evaluate the radial force for this potential INPUT: R - Galactocentric cylindrical radius z - vertical height phi - azimuth t - time OUTPUT: the radial force HISTORY: 2013-06-26 - Written - Bovy (IAS) """ r= numpy.sqrt(R*R+z*z) return -self._mass(r)*R/r**3. def _zforce(self,R,z,phi=0.,t=0.): """ NAME: _zforce PURPOSE: evaluate the vertical force for this potential INPUT: R - Galactocentric cylindrical radius z - vertical height phi - azimuth t - time OUTPUT: the vertical force HISTORY: 2013-06-26 - Written - Bovy (IAS) """ r= numpy.sqrt(R*R+z*z) return -self._mass(r)*z/r**3. def _R2deriv(self,R,z,phi=0.,t=0.): """ NAME: _Rderiv PURPOSE: evaluate the second radial derivative for this potential INPUT: R - Galactocentric cylindrical radius z - vertical height phi - azimuth t - time OUTPUT: the second radial derivative HISTORY: 2013-06-28 - Written - Bovy (IAS) """ r= numpy.sqrt(R*R+z*z) return 4.*numpy.pi*r**(-2.-self.alpha)*numpy.exp(-(r/self.rc)**2.)*R**2.\ +self._mass(r)/r**5.*(z**2.-2.*R**2.) def _z2deriv(self,R,z,phi=0.,t=0.): """ NAME: _z2deriv PURPOSE: evaluate the second vertical derivative for this potential INPUT: R - Galactocentric cylindrical radius z - vertical height phi - azimuth t- time OUTPUT: the second vertical derivative HISTORY: 2013-06-28 - Written - Bovy (IAS) """ r= numpy.sqrt(R*R+z*z) return 4.*numpy.pi*r**(-2.-self.alpha)*numpy.exp(-(r/self.rc)**2.)*z**2.\ +self._mass(r)/r**5.*(R**2.-2.*z**2.) def _Rzderiv(self,R,z,phi=0.,t=0.): """ NAME: _Rzderiv PURPOSE: evaluate the mixed R,z derivative for this potential INPUT: R - Galactocentric cylindrical radius z - vertical height phi - azimuth t- time OUTPUT: d2phi/dR/dz HISTORY: 2013-08-28 - Written - Bovy (IAS) """ r= numpy.sqrt(R*R+z*z) return R*z*(4.*numpy.pi*r**(-2.-self.alpha)*numpy.exp(-(r/self.rc)**2.) -3.*self._mass(r)/r**5.) def _dens(self,R,z,phi=0.,t=0.): """ NAME: _dens PURPOSE: evaluate the density for this potential INPUT: R - Galactocentric cylindrical radius z - vertical height phi - azimuth t - time OUTPUT: the density HISTORY: 2013-06-28 - Written - Bovy (IAS) """ r= numpy.sqrt(R**2.+z**2.) return 1./r**self.alpha*numpy.exp(-(r/self.rc)**2.) def _mass(self,R,z=None,t=0.): """ NAME: _mass PURPOSE: evaluate the mass within R for this potential INPUT: R - Galactocentric cylindrical radius z - vertical height t - time OUTPUT: the mass enclosed HISTORY: 2013-XX-XX - Written - Bovy (IAS) 2021-04-07 - Switched to hypergeometric function equivalent to incomplete gamma function for better behavior at alpha < -3 - Bovy (UofT) """ if z is not None: raise AttributeError # use general implementation return 2.*numpy.pi*R**(3.-self.alpha)/(1.5-self.alpha/2.)\ *special.hyp1f1(1.5-self.alpha/2.,2.5-self.alpha/2., -(R/self.rc)**2.) @kms_to_kpcGyrDecorator def _nemo_accpars(self,vo,ro): """ NAME: _nemo_accpars PURPOSE: return the accpars potential parameters for use of this potential with NEMO INPUT: vo - velocity unit in km/s ro - length unit in kpc OUTPUT: accpars string HISTORY: 2014-12-18 - Written - Bovy (IAS) """ ampl= self._amp*vo**2.*ro**(self.alpha-2.) return "0,%s,%s,%s" % (ampl,self.alpha,self.rc*ro)