Source code for galpy.potential.IsochronePotential

###############################################################################
#   IsochronePotential.py: The isochrone potential
#
#                                     - amp
#                          Phi(r)= ---------------------
#                                   b + sqrt{b^2+r^2}
###############################################################################
import numpy
from ..util import conversion
from .Potential import Potential
[docs]class IsochronePotential(Potential): """Class that implements the Isochrone potential .. math:: \\Phi(r) = -\\frac{\\mathrm{amp}}{b+\\sqrt{b^2+r^2}} with :math:`\\mathrm{amp} = GM` the total mass. """
[docs] def __init__(self,amp=1.,b=1.,normalize=False, ro=None,vo=None): """ NAME: __init__ PURPOSE: initialize an isochrone potential INPUT: amp= amplitude to be applied to the potential, the total mass (default: 1); can be a Quantity with units of mass or Gxmass b= scale radius of the isochrone potential (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-09-08 - Written - Bovy (IAS) """ Potential.__init__(self,amp=amp,ro=ro,vo=vo,amp_units='mass') b= conversion.parse_length(b,ro=self._ro) self.b= b self._scale= self.b self.b2= self.b**2. 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: 2013-09-08 - Written - Bovy (IAS) """ r2= R**2.+z**2. rb= numpy.sqrt(r2+self.b2) return -1./(self.b+rb) 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-09-08 - Written - Bovy (IAS) """ r2= R**2.+z**2. rb= numpy.sqrt(r2+self.b2) dPhidrr= -1./rb/(self.b+rb)**2. return dPhidrr*R 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-09-08 - Written - Bovy (IAS) """ r2= R**2.+z**2. rb= numpy.sqrt(r2+self.b2) dPhidrr= -1./rb/(self.b+rb)**2. return dPhidrr*z 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-09-08 - Written - Bovy (IAS) """ r2= R**2.+z**2. rb= numpy.sqrt(r2+self.b2) return -(-self.b**3.-self.b*z**2.+(2.*R**2.-z**2.-self.b**2.)*rb)/\ rb**3./(self.b+rb)**3. 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-09-08 - Written - Bovy (IAS) """ r2= R**2.+z**2. rb= numpy.sqrt(r2+self.b2) return -(-self.b**3.-self.b*R**2.-(R**2.-2.*z**2.+self.b**2.)*rb)/\ rb**3./(self.b+rb)**3. 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-09-08 - Written - Bovy (IAS) """ r2= R**2.+z**2. rb= numpy.sqrt(r2+self.b2) return -R*z*(self.b+3.*rb)/\ rb**3./(self.b+rb)**3. 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-09-08 - Written - Bovy (IAS) """ r2= R**2.+z**2. rb= numpy.sqrt(r2+self.b2) return (3.*(self.b+rb)*rb**2.-r2*(self.b+3.*rb))/\ rb**3./(self.b+rb)**3./4./numpy.pi 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) """ r2= R**2.+z**2. rb= numpy.sqrt(r2+self.b2) return self.b*((R*z)/r2-(self.b*R*z*(self.b**2+2.*R**2+z**2)) /((self.b**2+R**2)*r2*rb) +numpy.arctan(z/R)-numpy.arctan(self.b*z/R/rb))/R**3/2./numpy.pi