Source code for galpy.potential.DoubleExponentialDiskPotential

###############################################################################
#   DoubleExponentialDiskPotential.py: class that implements the double
#                                      exponential disk potential
#
#                                      rho(R,z) = rho_0 e^-R/h_R e^-|z|/h_z
###############################################################################
import numpy
from scipy import special
from ..util import conversion
from .Potential import Potential, check_potential_inputs_not_arrays

def _de_psi(t):
    return t*numpy.tanh(numpy.pi/2.*numpy.sinh(t))
def _de_psiprime(t):
    return (numpy.sinh(numpy.pi*numpy.sinh(t))
            +numpy.pi*t*numpy.cosh(t))/(numpy.cosh(numpy.pi*numpy.sinh(t))+1)

[docs]class DoubleExponentialDiskPotential(Potential): """Class that implements the double exponential disk potential .. math:: \\rho(R,z) = \\mathrm{amp}\\,\\exp\\left(-R/h_R-|z|/h_z\\right) """
[docs] def __init__(self,amp=1.,hr=1./3.,hz=1./16.,normalize=False, ro=None,vo=None, de_h=1e-3,de_n=10000): """ NAME: __init__ PURPOSE: initialize a double-exponential disk potential INPUT: amp - amplitude to be applied to the potential (default: 1); can be a Quantity with units of mass density or Gxmass density hr - disk scale-length (can be Quantity) hz - scale-height (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. de_h= (1e-3) step used in numerical integration (use 1000 for a lower accuracy version that is typically still high accuracy enough, but faster) de_b= (10000) number of points used in numerical integration ro=, vo= distance and velocity scales for translation into internal units (default from configuration file) OUTPUT: DoubleExponentialDiskPotential object HISTORY: 2010-04-16 - Written - Bovy (NYU) 2013-01-01 - Re-implemented using faster integration techniques - Bovy (IAS) 2020-12-24 - Re-implemented again using more accurate integration techniques for Bessel integrals - Bovy (UofT) """ Potential.__init__(self,amp=amp,ro=ro,vo=vo,amp_units='density') hr= conversion.parse_length(hr,ro=self._ro) hz= conversion.parse_length(hz,ro=self._ro) self.hasC= True self.hasC_dens= True self._hr= hr self._scale= self._hr self._hz= hz self._alpha= 1./self._hr self._beta= 1./self._hz self._zforceNotSetUp= True #We have not calculated a typical Kz yet # For double-exponential formula self._de_h= de_h self._de_n= de_n self._de_j0zeros= special.jn_zeros(0,self._de_n)/numpy.pi self._de_j1zeros= special.jn_zeros(1,self._de_n)/numpy.pi self._de_j0_xs= numpy.pi/self._de_h\ *_de_psi(self._de_h*self._de_j0zeros) self._de_j0_weights= 2./(numpy.pi*self._de_j0zeros\ *special.j1(numpy.pi*self._de_j0zeros)**2.)\ *special.j0(self._de_j0_xs)\ *_de_psiprime(self._de_h*self._de_j0zeros) self._de_j1_xs= numpy.pi/self._de_h\ *_de_psi(self._de_h*self._de_j1zeros) self._de_j1_weights= 2./(numpy.pi*self._de_j1zeros\ *special.jv(2,numpy.pi*self._de_j1zeros)**2.)\ *special.j1(self._de_j1_xs)\ *_de_psiprime(self._de_h*self._de_j1zeros) # Potential at zero in case we want that _gamma= self._beta/self._alpha _gamma2= _gamma**2. self._pot_zero= (2.*(_gamma-1.)*numpy.sqrt(1.+_gamma2) +2.*numpy.arctanh(1./numpy.sqrt(1.+_gamma2)) -numpy.log(1.-_gamma/numpy.sqrt(1.+_gamma2)) +numpy.log(1.+_gamma/numpy.sqrt(1.+_gamma2)))\ /(2.*(1.+_gamma2)**1.5) self._pot_zero*= -4.*numpy.pi/self._alpha**2. # Normalize? if normalize or \ (isinstance(normalize,(int,float)) \ and not isinstance(normalize,bool)): #pragma: no cover self.normalize(normalize)
def _evaluate(self,R,z,phi=0.,t=0.,dR=0,dphi=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: 2010-04-16 - Written - Bovy (NYU) 2012-12-26 - New method using Gaussian quadrature between zeros - Bovy (IAS) 2020-12-24 - New method using Ogata's Bessel integral formula - Bovy (UofT) """ if isinstance(R,(float,int)): floatIn= True R= numpy.array([R]) z= numpy.array([z]) else: if isinstance(z,float): z= z*numpy.ones_like(R) floatIn= False outShape= R.shape # this code can't do arbitrary shapes R= R.flatten() z= z.flatten() fun= lambda x: (self._alpha**2.+(x/R[:,numpy.newaxis])**2.)**-1.5\ *(self._beta*numpy.exp(-x/R[:,numpy.newaxis]*numpy.fabs(z[:,numpy.newaxis])) -x/R[:,numpy.newaxis]*numpy.exp(-self._beta*numpy.fabs(z[:,numpy.newaxis])))\ /(self._beta**2.-(x/R[:,numpy.newaxis])**2.) out= -4.*numpy.pi*self._alpha/R*\ numpy.nansum(fun(self._de_j0_xs)*self._de_j0_weights, axis=1) out[(R == 0)*(z == 0)]= self._pot_zero if floatIn: return out[0] else: return numpy.reshape(out,outShape) @check_potential_inputs_not_arrays def _Rforce(self,R,z,phi=0.,t=0.): """ NAME: Rforce PURPOSE: evaluate radial force K_R (R,z) INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: K_R (R,z) HISTORY: 2010-04-16 - Written - Bovy (NYU) 2012-12-26 - New method using Gaussian quadrature between zeros - Bovy (IAS) 2020-12-24 - New method using Ogata's Bessel integral formula - Bovy (UofT) """ fun= lambda x: x*(self._alpha**2.+(x/R)**2.)**-1.5\ *(self._beta*numpy.exp(-x/R*numpy.fabs(z)) -x/R*numpy.exp(-self._beta*numpy.fabs(z)))\ /(self._beta**2.-(x/R)**2.) return -4.*numpy.pi*self._alpha/R**2.\ *numpy.nansum(fun(self._de_j1_xs)*self._de_j1_weights) @check_potential_inputs_not_arrays def _zforce(self,R,z,phi=0.,t=0.): """ NAME: zforce PURPOSE: evaluate vertical force K_z (R,z) INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: K_z (R,z) HISTORY: 2010-04-16 - Written - Bovy (NYU) 2012-12-26 - New method using Gaussian quadrature between zeros - Bovy (IAS) 2020-12-24 - New method using Ogata's Bessel integral formula - Bovy (UofT) """ fun= lambda x: (self._alpha**2.+(x/R)**2.)**-1.5*x/R\ *(numpy.exp(-x/R*numpy.fabs(z)) -numpy.exp(-self._beta*numpy.fabs(z)))\ /(self._beta**2.-(x/R)**2.) out= -4.*numpy.pi*self._alpha*self._beta/R*\ numpy.nansum(fun(self._de_j0_xs)*self._de_j0_weights) if z > 0.: return out else: return -out @check_potential_inputs_not_arrays def _R2deriv(self,R,z,phi=0.,t=0.): """ NAME: R2deriv PURPOSE: evaluate R2 derivative INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: -d K_R (R,z) d R HISTORY: 2012-12-27 - Written - Bovy (IAS) 2020-12-24 - New method using Ogata's Bessel integral formula - Bovy (UofT) """ fun= lambda x: x**2*(self._alpha**2.+(x/R)**2.)**-1.5\ *(self._beta*numpy.exp(-x/R*numpy.fabs(z)) -x/R*numpy.exp(-self._beta*numpy.fabs(z)))\ /(self._beta**2.-(x/R)**2.) return 4.*numpy.pi*self._alpha/R**3.\ *numpy.nansum(fun(self._de_j0_xs)*self._de_j0_weights -fun(self._de_j1_xs)/self._de_j1_xs\ *self._de_j1_weights) @check_potential_inputs_not_arrays def _z2deriv(self,R,z,phi=0.,t=0.): """ NAME: z2deriv PURPOSE: evaluate z2 derivative INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: -d K_Z (R,z) d Z HISTORY: 2012-12-26 - Written - Bovy (IAS) 2020-12-24 - New method using Ogata's Bessel integral formula - Bovy (UofT) """ fun= lambda x: (self._alpha**2.+(x/R)**2.)**-1.5*x/R\ *(x/R*numpy.exp(-x/R*numpy.fabs(z)) -self._beta*numpy.exp(-self._beta*numpy.fabs(z)))\ /(self._beta**2.-(x/R)**2.) return -4.*numpy.pi*self._alpha*self._beta/R*\ numpy.nansum(fun(self._de_j0_xs)*self._de_j0_weights) @check_potential_inputs_not_arrays def _Rzderiv(self,R,z,phi=0.,t=0.): """ NAME: Rzderiv PURPOSE: evaluate the mixed R,z derivative INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: d2phi/dR/dz HISTORY: 2013-08-28 - Written - Bovy (IAS) 2020-12-24 - New method using Ogata's Bessel integral formula - Bovy (UofT) """ fun= lambda x: (self._alpha**2.+(x/R)**2.)**-1.5*(x/R)**2.\ *(numpy.exp(-x/R*numpy.fabs(z)) -numpy.exp(-self._beta*numpy.fabs(z)))\ /(self._beta**2.-(x/R)**2.) out= -4.*numpy.pi*self._alpha*self._beta/R*\ numpy.nansum(fun(self._de_j1_xs)*self._de_j1_weights) if z > 0.: return out else: return -out def _dens(self,R,z,phi=0.,t=0.): """ NAME: _dens PURPOSE: evaluate the density INPUT: R - Cylindrical Galactocentric radius z - vertical height phi - azimuth t - time OUTPUT: rho (R,z) HISTORY: 2010-08-08 - Written - Bovy (NYU) """ return numpy.exp(-self._alpha*R-self._beta*numpy.fabs(z)) 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: 2018-08-19 - Written - Bovy (UofT) """ return 2.*numpy.exp(-self._alpha*R)/self._beta\ *(1.-numpy.exp(-self._beta*numpy.fabs(z)))