Source code for galpy.potential.AnyAxisymmetricRazorThinDiskPotential

###############################################################################
#   AnyAxisymmetricRazorThinDiskPotential.py: class that implements the
#                                             potential of an arbitrary
#                                             axisymmetric, razor-thin disk
###############################################################################
import numpy
from scipy import integrate, special
from .Potential import Potential, check_potential_inputs_not_arrays, \
    _APY_LOADED
from ..util import conversion
if _APY_LOADED:
    from astropy import units
[docs]class AnyAxisymmetricRazorThinDiskPotential(Potential): """Class that implements the potential of an arbitrary axisymmetric, razor-thin disk with surface density :math:`\Sigma(R)`"""
[docs] def __init__(self,surfdens=lambda R: 1.5*numpy.exp(-R/0.5),amp=1., normalize=False,ro=None,vo=None): """ NAME: __init__ PURPOSE: Initialize the potential of an arbitrary axisymmetric disk INPUT: surfdens= (1.5 e^[-R/0.3]) function of a single variable that gives the surface density as a function of radius (can return a Quantity) amp= (1.) amplitude to be applied to the potential 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: AnyAxisymmetricRazorThinDiskPotential object HISTORY: 2021-01-04 - Written - Bovy (UofT) """ Potential.__init__(self,amp=amp,ro=ro,vo=vo) # Parse surface density: does it have units? does it expect them? if _APY_LOADED: _sdens_unit_input= False try: surfdens(1) except (units.UnitConversionError,units.UnitTypeError): _sdens_unit_input= True _sdens_unit_output= False if _sdens_unit_input: try: surfdens(1.*units.kpc).to(units.Msun/units.pc**2) except (AttributeError,units.UnitConversionError): pass else: _sdens_unit_output= True else: try: surfdens(1.).to(units.Msun/units.pc**2) except (AttributeError,units.UnitConversionError): pass else: _sdens_unit_output= True if _sdens_unit_input and _sdens_unit_output: self._sdens= lambda R: conversion.parse_surfdens(\ surfdens(R*self._ro*units.kpc), ro=self._ro,vo=self._vo) elif _sdens_unit_input: self._sdens= lambda R: surfdens(R*self._ro*units.kpc) elif _sdens_unit_output: self._sdens= lambda R: conversion.parse_surfdens(surfdens(R), ro=self._ro, vo=self._vo) if not hasattr(self,'_sdens'): # unitless self._sdens= surfdens # The potential at zero, in case it's asked for self._pot_zero= -2.*numpy.pi*integrate.quad(lambda a: self._sdens(a), 0,numpy.inf)[0] if normalize or \ (isinstance(normalize,(int,float)) \ and not isinstance(normalize,bool)): #pragma: no cover self.normalize(normalize)
@check_potential_inputs_not_arrays def _evaluate(self,R,z,phi=0.,t=0.): """ NAME: _evaluate PURPOSE: evaluate the potential at (R,z) INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: potential at (R,z) HISTORY: 2021-01-04 - Written - Bovy (UofT) """ if R == 0 and z == 0: return self._pot_zero elif numpy.isinf(R**2+z**2): return 0. potint= lambda a: a*self._sdens(a)\ /numpy.sqrt((R+a)**2.+z**2.)*special.ellipk(4*R*a/((R+a)**2.+z**2.)) return -4*(integrate.quad(potint,0,2*R,points=[R])[0] +integrate.quad(potint,2*R,numpy.inf)[0]) @check_potential_inputs_not_arrays def _Rforce(self,R,z,phi=0.,t=0.): """ NAME: _Rforce PURPOSE: evaluate the radial force at (R,z) INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: F_R at (R,z) HISTORY: 2021-01-04 - Written - Bovy (UofT) """ R2= R**2 z2= z**2 def rforceint(a): a2= a**2 aRz= (a+R)**2.+z2 faRoveraRz= 4*a*R/aRz return a*self._sdens(a)\ *((a2-R2+z2)*special.ellipe(faRoveraRz) -((a-R)**2+z2)*special.ellipk(faRoveraRz))\ /R/((a-R)**2+z2)/numpy.sqrt(aRz) return 2*(integrate.quad(rforceint,0,2*R,points=[R])[0] +integrate.quad(rforceint,2*R,numpy.inf)[0]) @check_potential_inputs_not_arrays def _zforce(self,R,z,phi=0.,t=0.): """ NAME: _zforce PURPOSE: evaluate the vertical force at (R,z) INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: F_z at (R,z) HISTORY: 2021-01-04 - Written - Bovy (UofT) """ if z == 0: return 0. z2= z**2 def zforceint(a): aRz= (a+R)**2.+z2 faRoveraRz= 4*a*R/aRz return a*self._sdens(a)\ *special.ellipe(faRoveraRz)/((a-R)**2+z2)/numpy.sqrt(aRz) return -4*z*(integrate.quad(zforceint,0,2*R,points=[R])[0] +integrate.quad(zforceint,2*R,numpy.inf)[0]) @check_potential_inputs_not_arrays def _R2deriv(self,R,z,phi=0.,t=0.): """ NAME: _R2deriv PURPOSE: evaluate the 2nd radial derivative at (R,z) INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: d2 Phi / dR2 at (R,z) HISTORY: 2021-01-04 - Written - Bovy (UofT) """ R2= R**2 z2= z**2 def r2derivint(a): a2= a**2 aRz= (a+R)**2.+z2 faRoveraRz= 4*a*R/aRz return a*self._sdens(a)\ *(-(((a2-3.*R2)*(a2-R2)**2+(3.*a2**2+2.*a2*R2+3.*R2**2)*z2 +(3.*a2+7.*R2)*z**4+z**6)*special.ellipe(faRoveraRz)) +((a-R)**2+z2)*((a2-R2)**2+2.*(a2+2.*R2)*z2+z**4) *special.ellipk(faRoveraRz))\ /(2.*R2*((a-R)**2+z2)**2*((a+R)**2+z2)**1.5) return -4*(integrate.quad(r2derivint,0,2*R,points=[R])[0] +integrate.quad(r2derivint,2*R,numpy.inf)[0]) @check_potential_inputs_not_arrays def _z2deriv(self,R,z,phi=0.,t=0.): """ NAME: _z2deriv PURPOSE: evaluate the 2nd vertical derivative at (R,z) INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: d2 Phi / dz2 at (R,z) HISTORY: 2021-01-04 - Written - Bovy (UofT) """ R2= R**2 z2= z**2 def z2derivint(a): a2= a**2 aRz= (a+R)**2.+z2 faRoveraRz= 4*a*R/aRz return a*self._sdens(a)\ *(-(((a2-R2)**2-2.*(a2+R2)*z2-3.*z**4)*special.ellipe(faRoveraRz)) -z2*((a-R)**2+z2)*special.ellipk(faRoveraRz))\ /(((a-R)**2+z2)**2*((a+R)**2+z2)**1.5) return -4*(integrate.quad(z2derivint,0,2*R,points=[R])[0] +integrate.quad(z2derivint,2*R,numpy.inf)[0]) @check_potential_inputs_not_arrays def _Rzderiv(self,R,z,phi=0.,t=0.): """ NAME: _Rzderiv PURPOSE: evaluate the mixed radial, vertical derivative at (R,z) INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: d2 Phi / dRdz at (R,z) HISTORY: 2021-01-04 - Written - Bovy (UofT) """ R2= R**2 z2= z**2 def rzderivint(a): a2= a**2 aRz= (a+R)**2.+z2 faRoveraRz= 4*a*R/aRz return a*self._sdens(a)\ *(-((a**4-7.*R**4-6.*R2*z2+z**4+2.*a2*(3.*R2+z2)) *special.ellipe(faRoveraRz)) +((a-R)**2+z**2)*(a2-R2+z2)*special.ellipk(faRoveraRz))\ /R/((a-R)**2+z2)**2/((a+R)**2+z2)**1.5 return -2*z*(integrate.quad(rzderivint,0,2*R,points=[R])[0] +integrate.quad(rzderivint,2*R,numpy.inf)[0]) def _surfdens(self,R,z,phi=0.,t=0.): """ NAME: _surfdens PURPOSE: evaluate the surface density INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: Sigma (R,z) HISTORY: 2021-01-04 - Written - Bovy (UofT) """ return self._sdens(R)