Source code for galpy.potential.MN3ExponentialDiskPotential

###############################################################################
#   MN3ExponentialDiskPotential.py: class that implements the three Miyamoto-
#                                   Nagai approximation to a radially 
#                                   exponential disk potential of Smith et al.
#                                   2015
###############################################################################
import warnings
import numpy
from ..util import galpyWarning, conversion
from .Potential import Potential, kms_to_kpcGyrDecorator
from .MiyamotoNagaiPotential import MiyamotoNagaiPotential
[docs]class MN3ExponentialDiskPotential(Potential): """class that implements the three Miyamoto-Nagai approximation to a radially-exponential disk potential of `Smith et al. 2015 <http://adsabs.harvard.edu/abs/2015arXiv150200627S>`_ .. math:: \\rho(R,z) = \\mathrm{amp}\\,\\exp\\left(-R/h_R-|z|/h_z\\right) or .. math:: \\rho(R,z) = \\mathrm{amp}\\,\\exp\\left(-R/h_R\\right)\\mathrm{sech}^2\\left(-|z|/h_z\\right) depending on whether sech=True or not. This density is approximated using three Miyamoto-Nagai disks """
[docs] def __init__(self,amp=1.,hr=1./3.,hz=1./16., sech=False,posdens=False, normalize=False, ro=None,vo=None): """ NAME: __init__ PURPOSE: initialize a 3MN approximation to an 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) sech= (False) if True, hz is the scale height of a sech vertical profile (default is exponential vertical profile) posdens= (False) if True, allow only positive density solutions (Table 2 in Smith et al. rather than Table 1) 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: MN3ExponentialDiskPotential object HISTORY: 2015-02-07 - Written - Bovy (IAS) """ 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._hr= hr self._hz= hz self._scale= self._hr # Adjust amp for definition self._amp*= 4.*numpy.pi*self._hr**2.*self._hz # First determine b/rd if sech: self._brd= _b_sechhz(self._hz/self._hr) else: self._brd= _b_exphz(self._hz/self._hr) if self._brd < 0.: raise IOError("MN3ExponentialDiskPotential's b/Rd is negative for the given hz") # Check range if (not posdens and self._brd > 3.) \ or (posdens and self._brd > 1.35): warnings.warn("MN3ExponentialDiskPotential's b/Rd = %g is outside of the interpolation range of Smith et al. (2015)" % self._brd, galpyWarning) self._b= self._brd*self._hr # Now setup the various MN disks if posdens: self._mn3= [MiyamotoNagaiPotential(amp=_mass1_tab2(self._brd), a=_a1_tab2(self._brd)*self._hr, b=self._b), MiyamotoNagaiPotential(amp=_mass2_tab2(self._brd), a=_a2_tab2(self._brd)*self._hr, b=self._b), MiyamotoNagaiPotential(amp=_mass3_tab2(self._brd), a=_a3_tab2(self._brd)*self._hr, b=self._b)] else: self._mn3= [MiyamotoNagaiPotential(amp=_mass1_tab1(self._brd), a=_a1_tab1(self._brd)*self._hr, b=self._b), MiyamotoNagaiPotential(amp=_mass2_tab1(self._brd), a=_a2_tab1(self._brd)*self._hr, b=self._b), MiyamotoNagaiPotential(amp=_mass3_tab1(self._brd), a=_a3_tab1(self._brd)*self._hr, b=self._b)] if normalize or \ (isinstance(normalize,(int,float)) \ and not isinstance(normalize,bool)): self.normalize(normalize) self.hasC= True self.hasC_dxdv= True self._nemo_accname= 'MiyamotoNagai+MiyamotoNagai+MiyamotoNagai' return None
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: 2015-02-07 - Written - Bovy (IAS) """ return self._mn3[0](R,z,phi=phi,t=t)\ +self._mn3[1](R,z,phi=phi,t=t)\ +self._mn3[2](R,z,phi=phi,t=t) 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: 2015-02-07 - Written - Bovy (IAS) """ return self._mn3[0].Rforce(R,z,phi=phi,t=t)\ +self._mn3[1].Rforce(R,z,phi=phi,t=t)\ +self._mn3[2].Rforce(R,z,phi=phi,t=t) 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: 2015-02-07 - Written - Bovy (IAS) """ return self._mn3[0].zforce(R,z,phi=phi,t=t)\ +self._mn3[1].zforce(R,z,phi=phi,t=t)\ +self._mn3[2].zforce(R,z,phi=phi,t=t) 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: 2015-02-07 - Written - Bovy (IAS) """ return self._mn3[0].dens(R,z,phi=phi,t=t)\ +self._mn3[1].dens(R,z,phi=phi,t=t)\ +self._mn3[2].dens(R,z,phi=phi,t=t) def _R2deriv(self,R,z,phi=0.,t=0.): """ NAME: _R2deriv 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: 2015-02-07 - Written - Bovy (IAS) """ return self._mn3[0].R2deriv(R,z,phi=phi,t=t)\ +self._mn3[1].R2deriv(R,z,phi=phi,t=t)\ +self._mn3[2].R2deriv(R,z,phi=phi,t=t) 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: 2015-02-07 - Written - Bovy (IAS) """ return self._mn3[0].z2deriv(R,z,phi=phi,t=t)\ +self._mn3[1].z2deriv(R,z,phi=phi,t=t)\ +self._mn3[2].z2deriv(R,z,phi=phi,t=t) 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: 2015-02-07 - Written - Bovy (IAS) """ return self._mn3[0].Rzderiv(R,z,phi=phi,t=t)\ +self._mn3[1].Rzderiv(R,z,phi=phi,t=t)\ +self._mn3[2].Rzderiv(R,z,phi=phi,t=t) @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: 2015-02-09 - Written - Bovy (IAS) """ out= "" # Loop through the self._mn3 MN potentials for ii in range(3): if ii > 0: out+= '#' ampl= self._amp*self._mn3[ii]._amp*vo**2.*ro out+= "0,%s,%s,%s" % (ampl,self._mn3[ii]._a*ro,self._mn3[ii]._b*ro) return out
# Equations from Table 1 def _mass1_tab1(brd): return -0.0090*brd**4.+0.0640*brd**3.-0.1653*brd**2.+0.1164*brd+1.9487 def _mass2_tab1(brd): return 0.0173*brd**4.-0.0903*brd**3.+0.0877*brd**2.+0.2029*brd-1.3077 def _mass3_tab1(brd): return -0.0051*brd**4.+0.0287*brd**3.-0.0361*brd**2.-0.0544*brd+0.2242 def _a1_tab1(brd): return -0.0358*brd**4.+0.2610*brd**3.-0.6987*brd**2.-0.1193*brd+2.0074 def _a2_tab1(brd): return -0.0830*brd**4.+0.4992*brd**3.-0.7967*brd**2.-1.2966*brd+4.4441 def _a3_tab1(brd): return -0.0247*brd**4.+0.1718*brd**3.-0.4124*brd**2.-0.5944*brd+0.7333 # Equations from Table 2 def _mass1_tab2(brd): return 0.0036*brd**4.-0.0330*brd**3.+0.1117*brd**2.-0.1335*brd+0.1749 def _mass2_tab2(brd): return -0.0131*brd**4.+0.1090*brd**3.-0.3035*brd**2.+0.2921*brd-5.7976 def _mass3_tab2(brd): return -0.0048*brd**4.+0.0454*brd**3.-0.1425*brd**2.+0.1012*brd+6.7120 def _a1_tab2(brd): return -0.0158*brd**4.+0.0993*brd**3.-0.2070*brd**2.-0.7089*brd+0.6445 def _a2_tab2(brd): return -0.0319*brd**4.+0.1514*brd**3.-0.1279*brd**2.-0.9325*brd+2.6836 def _a3_tab2(brd): return -0.0326*brd**4.+0.1816*brd**3.-0.2943*brd**2.-0.6329*brd+2.3193 # Equations to go from hz to b def _b_exphz(hz): return -0.269*hz**3.+1.080*hz**2.+1.092*hz def _b_sechhz(hz): return -0.033*hz**3.+0.262*hz**2.+0.659*hz