###############################################################################
# DehnenBarPotential: Dehnen (2000)'s bar potential
###############################################################################
import numpy
from ..util import conversion
from .Potential import Potential
_degtorad= numpy.pi/180.
[docs]class DehnenBarPotential(Potential):
"""Class that implements the Dehnen bar potential (`Dehnen 2000 <http://adsabs.harvard.edu/abs/2000AJ....119..800D>`__), generalized to 3D following `Monari et al. (2016) <http://adsabs.harvard.edu/abs/2016MNRAS.461.3835M>`__
.. math::
\\Phi(R,z,\\phi) = A_b(t)\\,\\cos\\left(2\\,(\\phi-\\Omega_b\\,t)\\right))\\,\\left(\\frac{R}{r}\\right)^2\\,\\times \\begin{cases}
-(R_b/r)^3\\,, & \\text{for}\\ r \\geq R_b\\\\
(r/R_b)^3-2\\,, & \\text{for}\\ r\\leq R_b.
\\end{cases}
where :math:`r^2 = R^2+z^2` is the spherical radius and
.. math::
A_b(t) = A_f\\,\\left(\\frac{3}{16}\\xi^5-\\frac{5}{8}\\xi^3+\\frac{15}{16}\\xi+\\frac{1}{2}\\right)\\,, \\xi = 2\\frac{t/T_b-t_\\mathrm{form}}{T_\mathrm{steady}}-1\\,,\ \mathrm{if}\ t_\\mathrm{form} \\leq \\frac{t}{T_b} \\leq t_\\mathrm{form}+T_\\mathrm{steady}
and
.. math::
A_b(t) = \\begin{cases}
0\\,, & \\frac{t}{T_b} < t_\mathrm{form}\\\\
A_f\\,, & \\frac{t}{T_b} > t_\mathrm{form}+T_\mathrm{steady}
\\end{cases}
where
.. math::
T_b = \\frac{2\pi}{\\Omega_b}
is the bar period and the strength can also be specified using :math:`\\alpha`
.. math::
\\alpha = 3\\,\\frac{A_f}{v_0^2}\\,\\left(\\frac{R_b}{r_0}\\right)^3\,.
"""
normalize= property() # turn off normalize
[docs] def __init__(self,amp=1.,omegab=None,rb=None,chi=0.8,
rolr=0.9,barphi=25.*_degtorad,
tform=-4.,tsteady=None,beta=0.,
alpha=0.01,Af=None,ro=None,vo=None):
"""
NAME:
__init__
PURPOSE:
initialize a Dehnen bar potential
INPUT:
amp - amplitude to be applied to the potential (default:
1., see alpha or Ab below)
barphi - angle between sun-GC line and the bar's major axis
(in rad; default=25 degree; or can be Quantity))
tform - start of bar growth / bar period (default: -4)
tsteady - time from tform at which the bar is fully grown / bar period (default: -tform/2, st the perturbation is fully grown at tform/2)
Either provide:
a) rolr - radius of the Outer Lindblad Resonance for a
circular orbit (can be Quantity)
chi - fraction R_bar / R_CR (corotation radius of bar)
alpha - relative bar strength (default: 0.01)
beta - power law index of rotation curve (to
calculate OLR, etc.)
b) omegab - rotation speed of the bar (can be Quantity)
rb - bar radius (can be Quantity)
Af - bar strength (can be Quantity)
OUTPUT:
(none)
HISTORY:
2010-11-24 - Started - Bovy (NYU)
2017-06-23 - Converted to 3D following Monari et al. (2016) - Bovy (UofT/CCA)
"""
Potential.__init__(self,amp=amp,ro=ro,vo=vo)
barphi= conversion.parse_angle(barphi)
rolr= conversion.parse_length(rolr,ro=self._ro)
rb= conversion.parse_length(rb,ro=self._ro)
omegab= conversion.parse_frequency(omegab,ro=self._ro,vo=self._vo)
Af= conversion.parse_energy(Af,vo=self._vo)
self.hasC= True
self.hasC_dxdv= True
self.isNonAxi= True
self._barphi= barphi
if omegab is None:
self._rolr= rolr
self._chi= chi
self._beta= beta
#Calculate omegab and rb
self._omegab= 1./((self._rolr**(1.-self._beta))/(1.+numpy.sqrt((1.+self._beta)/2.)))
self._rb= self._chi*self._omegab**(1./(self._beta-1.))
self._alpha= alpha
self._af= self._alpha/3./self._rb**3.
else:
self._omegab= omegab
self._rb= rb
self._af= Af
self._tb= 2.*numpy.pi/self._omegab
self._tform= tform*self._tb
if tsteady is None:
self._tsteady= self._tform/2.
else:
self._tsteady= self._tform+tsteady*self._tb
def _smooth(self,t):
if isinstance(t,numpy.ndarray):
smooth=numpy.ones(len(t))
indx=(t < self._tform)
smooth[indx]=0.
indx=(t < self._tsteady) * (t >= self._tform)
deltat=t[indx]-self._tform
xi= 2.*deltat/(self._tsteady-self._tform)-1.
smooth[indx]= (3./16.*xi**5.-5./8*xi**3.+15./16.*xi+.5)
else:
if t < self._tform:
smooth= 0.
elif t < self._tsteady:
deltat= t-self._tform
xi= 2.*deltat/(self._tsteady-self._tform)-1.
smooth= (3./16.*xi**5.-5./8*xi**3.+15./16.*xi+.5)
else: #bar is fully on
smooth= 1.
return smooth
def _evaluate(self,R,z,phi=0.,t=0.):
"""
NAME:
_evaluate
PURPOSE:
evaluate the potential at R,phi,t
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
OUTPUT:
Phi(R,z,phi,t)
HISTORY:
2010-11-24 - Started - Bovy (NYU)
"""
#Calculate relevant time
smooth=self._smooth(t)
r2= R**2.+z**2.
r= numpy.sqrt(r2)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= ((r[indx]/self._rb)**3.-2.)\
*numpy.divide(R[indx]**2.,r2[indx],numpy.ones_like(R[indx]),
where=R[indx]!=0)
indx=numpy.invert(indx)
out[indx]= -(self._rb/r[indx])**3.*1./(1.+z[indx]**2./R[indx]**2.)
out*=self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))
return out
else:
if r == 0:
return -2.*self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))
elif r <= self._rb:
return self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))\
*((r/self._rb)**3.-2.)*R**2./r2
else:
return -self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3.\
*1./(1.+z**2./R**2.)
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:
2010-11-24 - Written - Bovy (NYU)
"""
#Calculate relevant time
smooth=self._smooth(t)
r= numpy.sqrt(R**2.+z**2.)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= -((r[indx]/self._rb)**3.*R[indx]*(3.*R[indx]**2.+2.*z[indx]**2.)-4.*R[indx]*z[indx]**2.)/r[indx]**4.
indx= numpy.invert(indx)
out[indx]= -(self._rb/r[indx])**3.*R[indx]/r[indx]**4.*(3.*R[indx]**2.-2.*z[indx]**2.)
out*=self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))
return out
else:
if r <= self._rb:
return -self._af*smooth*numpy.cos(2.*(phi-self._omegab*t
-self._barphi))\
*((r/self._rb)**3.*R*(3.*R**2.+2.*z**2.)-4.*R*z**2.)/r**4.
else:
return -self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3.*R/r**4.*(3.*R**2.-2.*z**2.)
def _phiforce(self,R,z,phi=0.,t=0.):
"""
NAME:
_phiforce
PURPOSE:
evaluate the azimuthal force for this potential
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
OUTPUT:
the azimuthal force
HISTORY:
2010-11-24 - Written - Bovy (NYU)
"""
#Calculate relevant time
smooth=self._smooth(t)
r2= R**2.+z**2.
r= numpy.sqrt(r2)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= ((r[indx]/self._rb)**3.-2.)*R[indx]**2./r2[indx]
indx=numpy.invert(indx)
out[indx]= -(self._rb/r[indx])**3.*R[indx]**2./r2[indx]
out*=2.*self._af*smooth*numpy.sin(2.*(phi-self._omegab*t-self._barphi))
return out
else:
if r <= self._rb:
return 2.*self._af*smooth*numpy.sin(2.*(phi-self._omegab*t-
self._barphi))\
*((r/self._rb)**3.-2.)*R**2./r2
else:
return -2.*self._af*smooth*numpy.sin(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3.*R**2./r2
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:
2017-06-23 - Written - Bovy (NYU)
"""
#Calculate relevant time
smooth=self._smooth(t)
r= numpy.sqrt(R**2.+z**2.)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= -((r[indx]/self._rb)**3.+4.)*R[indx]**2.*z[indx]/r[indx]**4.
indx=numpy.invert(indx)
out[indx]= -5.*(self._rb/r[indx])**3.*R[indx]**2.*z[indx]/r[indx]**4.
out*=self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))
return out
else:
if r <= self._rb:
return -self._af*smooth*numpy.cos(2.*(phi-self._omegab*t
-self._barphi))\
*((r/self._rb)**3.+4.)*R**2.*z/r**4.
else:
return -5.*self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3.*R**2.*z/r**4.
def _R2deriv(self,R,z,phi=0.,t=0.):
#Calculate relevant time
smooth=self._smooth(t)
r= numpy.sqrt(R**2.+z**2.)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= ((r[indx]/self._rb)**3.*((9.*R[indx]**2.+2.*z[indx]**2.)/r[indx]**4.
-R[indx]**2./r[indx]**6.*(3.*R[indx]**2.+2.*z[indx]**2.))\
+4.*z[indx]**2./r[indx]**6.*(4.*R[indx]**2.-r[indx]**2.))
indx=numpy.invert(indx)
out[indx]= (self._rb/r[indx])**3./r[indx]**6.*((r[indx]**2.-7.*R[indx]**2.)*(3.*R[indx]**2.-2.*z[indx]**2.)\
+6.*R[indx]**2.*r[indx]**2.)
out*=self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))
return out
else:
if r <= self._rb:
return self._af*smooth*numpy.cos(2.*(phi-self._omegab*t
-self._barphi))\
*((r/self._rb)**3.*((9.*R**2.+2.*z**2.)/r**4.
-R**2./r**6.*(3.*R**2.+2.*z**2.))\
+4.*z**2./r**6.*(4.*R**2.-r**2.))
else:
return self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3./r**6.*((r**2.-7.*R**2.)*(3.*R**2.-2.*z**2.)\
+6.*R**2.*r**2.)
def _phi2deriv(self,R,z,phi=0.,t=0.):
#Calculate relevant time
smooth=self._smooth(t)
r= numpy.sqrt(R**2.+z**2.)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= -((r[indx]/self._rb)**3.-2.)*R[indx]**2./r[indx]**2.
indx=numpy.invert(indx)
out[indx]= (self._rb/r[indx])**3.*R[indx]**2./r[indx]**2.
out*=4.*self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))
return out
else:
if r <= self._rb:
return -4.*self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-
self._barphi))\
*((r/self._rb)**3.-2.)*R**2./r**2.
else:
return 4.*self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3.*R**2./r**2.
def _Rphideriv(self,R,z,phi=0.,t=0.):
#Calculate relevant time
smooth=self._smooth(t)
r= numpy.sqrt(R**2.+z**2.)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= ((r[indx]/self._rb)**3.*R[indx]*(3.*R[indx]**2.+2.*z[indx]**2.)-4.*R[indx]*z[indx]**2.)/r[indx]**4.
indx=numpy.invert(indx)
out[indx]= (self._rb/r[indx])**3.*R[indx]/r[indx]**4.*(3.*R[indx]**2.-2.*z[indx]**2.)
out*=-2.*self._af*smooth*numpy.sin(2.*(phi-self._omegab*t-self._barphi))
return out
else:
if r <= self._rb:
return -2.*self._af*smooth*numpy.sin(2.*(phi-self._omegab*t
-self._barphi))\
*((r/self._rb)**3.*R*(3.*R**2.+2.*z**2.)-4.*R*z**2.)/r**4.
else:
return -2.*self._af*smooth*numpy.sin(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3.*R/r**4.*(3.*R**2.-2.*z**2.)
def _z2deriv(self,R,z,phi=0.,t=0.):
#Calculate relevant time
smooth=self._smooth(t)
r= numpy.sqrt(R**2.+z**2.)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= R[indx]**2./r[indx]**6.*((r[indx]/self._rb)**3.*(r[indx]**2.-z[indx]**2.)
+4.*(r[indx]**2.-4.*z[indx]**2.))
indx=numpy.invert(indx)
out[indx]=5.*(self._rb/r[indx])**3.*R[indx]**2./r[indx]**6.*(r[indx]**2.-7.*z[indx]**2.)
out*=self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))
return out
else:
if r <= self._rb:
return self._af*smooth*numpy.cos(2.*(phi-self._omegab*t
-self._barphi))\
*R**2./r**6.*((r/self._rb)**3.*(r**2.-z**2.)
+4.*(r**2.-4.*z**2.))
else:
return 5.*self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3.*R**2./r**6.*(r**2.-7.*z**2.)
def _Rzderiv(self,R,z,phi=0.,t=0.):
#Calculate relevant time
smooth=self._smooth(t)
r= numpy.sqrt(R**2.+z**2.)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= R[indx]*z[indx]/r[indx]**6.*((r[indx]/self._rb)**3.*(2.*r[indx]**2.-R[indx]**2.)
+8.*(r[indx]**2.-2.*R[indx]**2.))
indx=numpy.invert(indx)
out[indx]= 5.*(self._rb/r[indx])**3.*R[indx]*z[indx]/r[indx]**6.*(2.*r[indx]**2.-7.*R[indx]**2.)
out*=self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-self._barphi))
return out
else:
if r <= self._rb:
return self._af*smooth*numpy.cos(2.*(phi-self._omegab*t
-self._barphi))\
*R*z/r**6.*((r/self._rb)**3.*(2.*r**2.-R**2.)
+8.*(r**2.-2.*R**2.))
else:
return 5.*self._af*smooth*numpy.cos(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3.*R*z/r**6.*(2.*r**2.-7.*R**2.)
def _phizderiv(self,R,z,phi=0.,t=0.):
"""
NAME:
_phizderiv
PURPOSE:
evaluate the mixed azimuthal, vertical derivative for this potential
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
OUTPUT:
the mixed azimuthal, vertical derivative
HISTORY:
2021-04-30 - Written - Bovy (UofT)
"""
#Calculate relevant time
smooth=self._smooth(t)
r= numpy.sqrt(R**2.+z**2.)
if isinstance(r,numpy.ndarray):
if not isinstance(R,numpy.ndarray):
R=numpy.repeat(R,len(r))
if not isinstance(z,numpy.ndarray):
z=numpy.repeat(z,len(r))
out=numpy.empty(len(r))
indx= r <= self._rb
out[indx]= -((r[indx]/self._rb)**3.+4.)*R[indx]**2.*z[indx]/r[indx]**4.
indx=numpy.invert(indx)
out[indx]= -5.*(self._rb/r[indx])**3.*R[indx]**2.*z[indx]/r[indx]**4.
out*=self._af*smooth*numpy.sin(2.*(phi-self._omegab*t-self._barphi))
return 2.*out
else:
if r <= self._rb:
return -2*self._af*smooth*numpy.sin(2.*(phi-self._omegab*t
-self._barphi))\
*((r/self._rb)**3.+4.)*R**2.*z/r**4.
else:
return -10.*self._af*smooth*numpy.sin(2.*(phi-self._omegab*t-
self._barphi))\
*(self._rb/r)**3.*R**2.*z/r**4.
def tform(self): #pragma: no cover
"""
NAME:
tform
PURPOSE:
return formation time of the bar
INPUT:
(none)
OUTPUT:
tform in normalized units
HISTORY:
2011-03-08 - Written - Bovy (NYU)
"""
return self._tform
def OmegaP(self):
"""
NAME:
OmegaP
PURPOSE:
return the pattern speed
INPUT:
(none)
OUTPUT:
pattern speed
HISTORY:
2011-10-10 - Written - Bovy (IAS)
"""
return self._omegab