Source code for galpy.potential.HomogeneousSpherePotential

###############################################################################
#   HomogeneousSpherePotential.py: The potential of a homogeneous sphere
###############################################################################
import numpy
from ..util import conversion
from .Potential import Potential
[docs]class HomogeneousSpherePotential(Potential): """Class that implements the homogeneous sphere potential for :math:`\\rho(r) = \\rho_0 = \\mathrm{constant}` for all :math:`r < R` and zero otherwise. The potential is given by .. math:: \\Phi(r) = \\mathrm{amp}\\times\\left\\{\\begin{array}{lr} (r^2-3R^2), & \\text{for } r < R\\\\ -\\frac{2R^3}{r} & \\text{for } r \\geq R \\end{array}\\right. We have that :math:`\\rho_0 = 3\\,\\mathrm{amp}/[2\\pi G]`. """
[docs] def __init__(self,amp=1.,R=1.1,normalize=False, ro=None,vo=None): """ NAME: __init__ PURPOSE: initialize a homogeneous sphere potential INPUT: amp= amplitude to be applied to the potential (default: 1); can be a Quantity with units of mass density or Gxmass density R= size of the sphere (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: 2019-12-20 - Written - Bovy (UofT) """ Potential.__init__(self,amp=amp,ro=ro,vo=vo,amp_units='density') R= conversion.parse_length(R,ro=self._ro) self.R= R self._R2= self.R**2. self._R3= self.R**3. 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
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: 2019-12-20 - Written - Bovy (UofT) """ r2= R**2.+z**2. if r2 < self._R2: return r2-3.*self._R2 else: return -2.*self._R3/numpy.sqrt(r2) 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: 2019-12-20 - Written - Bovy (UofT) """ r2= R**2.+z**2. if r2 < self._R2: return -2.*R else: return -2.*self._R3*R/r2**1.5 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: 2019-12-20 - Written - Bovy (UofT) """ r2= R**2.+z**2. if r2 < self._R2: return -2.*z else: return -2.*self._R3*z/r2**1.5 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: 2019-12-20 - Written - Bovy (UofT) """ r2= R**2.+z**2. if r2 < self._R2: return 2. else: return 2.*self._R3/r2**1.5-6.*self._R3*R**2./r2**2.5 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: 2019-12-20 - Written - Bovy (UofT) """ r2= R**2.+z**2. if r2 < self._R2: return 2. else: return 2.*self._R3/r2**1.5-6.*self._R3*z**2./r2**2.5 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: 2019-12-20 - Written - Bovy (UofT) """ r2= R**2.+z**2. if r2 < self._R2: return 0. else: return -6.*self._R3*R*z/r2**2.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: 2019-12-20 - Written - Bovy (UofT) """ r2= R**2.+z**2. if r2 < self._R2: return 1.5/numpy.pi else: return 0.