###############################################################################
# CosmphiDiskPotential: cos(mphi) potential
###############################################################################
import numpy
from .planarPotential import planarPotential, _APY_LOADED
if _APY_LOADED:
from astropy import units
_degtorad= numpy.pi/180.
[docs]class CosmphiDiskPotential(planarPotential):
"""Class that implements the disk potential
.. math::
\\Phi(R,\\phi) = \\mathrm{amp}\\,\\phi_0\\,\\,\\cos\\left[m\\,(\\phi-\\phi_b)\\right]\\times \\begin{cases}
\\left(\\frac{R}{R_1}\\right)^p\\,, & \\text{for}\\ R \\geq R_b\\\\
\\left[2-\\left(\\frac{R_b}{R}\\right)^p\\right]\\times\\left(\\frac{R_b}{R_1}\\right)^p\\,, & \\text{for}\\ R\\leq R_b.
\\end{cases}
This potential can be grown between :math:`t_{\mathrm{form}}` and :math:`t_{\mathrm{form}}+T_{\mathrm{steady}}` in a similar way as DehnenBarPotential by wrapping it with a DehnenSmoothWrapperPotential
"""
[docs] def __init__(self,amp=1.,phib=25.*_degtorad,
p=1.,phio=0.01,m=4,r1=1.,rb=None,
cp=None,sp=None,ro=None,vo=None):
"""
NAME:
__init__
PURPOSE:
initialize an cosmphi disk potential
INPUT:
amp= amplitude to be applied to the potential (default:
1.), degenerate with phio below, but kept for overall
consistency with potentials
m= cos( m * (phi - phib) ), integer
p= power-law index of the phi(R) = (R/Ro)^p part
r1= (1.) normalization radius for the amplitude (can be Quantity); amp x phio is only the potential at (R,phi) = (r1,pib) when r1 > rb; otherwise more complicated
rb= (None) if set, break radius for power-law: potential R^p at R > Rb, R^-p at R < Rb, potential and force continuous at Rb
Either:
a) phib= angle (in rad; default=25 degree; or can be Quantity)
phio= potential perturbation (in terms of phio/vo^2 if vo=1 at Ro=1; or can be Quantity with units of velocity-squared)
b) cp, sp= m * phio * cos(m * phib), m * phio * sin(m * phib); can be Quantity with units of velocity-squared)
OUTPUT:
(none)
HISTORY:
2011-10-27 - Started - Bovy (IAS)
2017-09-16 - Added break radius rb - Bovy (UofT)
"""
planarPotential.__init__(self,amp=amp,ro=ro,vo=vo)
if _APY_LOADED and isinstance(phib,units.Quantity):
phib= phib.to(units.rad).value
if _APY_LOADED and isinstance(r1,units.Quantity):
r1= r1.to(units.kpc).value/self._ro
if _APY_LOADED and isinstance(rb,units.Quantity):
rb= rb.to(units.kpc).value/self._ro
if _APY_LOADED and isinstance(phio,units.Quantity):
phio= phio.to(units.km**2/units.s**2).value/self._vo**2.
if _APY_LOADED and isinstance(cp,units.Quantity):
cp= cp.to(units.km**2/units.s**2).value/self._vo**2.
if _APY_LOADED and isinstance(sp,units.Quantity):
sp= sp.to(units.km**2/units.s**2).value/self._vo**2.
# Back to old definition
self._r1p= r1**p
self._amp/= self._r1p
self.hasC= False
self._m= int(m) # make sure this is an int
if cp is None or sp is None:
self._phib= phib
self._mphio= phio*self._m
else:
self._mphio= numpy.sqrt(cp*cp+sp*sp)
self._phib= numpy.arctan(sp/cp)/self._m
if m < 2. and cp < 0.:
self._phib= numpy.pi+self._phib
self._p= p
if rb is None:
self._rb= 0.
self._rbp= 1. # never used, but for p < 0 general expr fails
self._rb2p= 1.
else:
self._rb= rb
self._rbp= self._rb**self._p
self._rb2p= self._rbp**2.
self._mphib= self._m*self._phib
self.hasC= True
self.hasC_dxdv= True
def _evaluate(self,R,phi=0.,t=0.):
"""
NAME:
_evaluate
PURPOSE:
evaluate the potential at R,phi,t
INPUT:
R - Galactocentric cylindrical radius
phi - azimuth
t - time
OUTPUT:
Phi(R,phi,t)
HISTORY:
2011-10-19 - Started - Bovy (IAS)
"""
if R < self._rb:
return self._mphio/self._m*numpy.cos(self._m*phi-self._mphib)\
*self._rbp*(2.*self._r1p-self._rbp/R**self._p)
else:
return self._mphio/self._m*R**self._p\
*numpy.cos(self._m*phi-self._mphib)
def _Rforce(self,R,phi=0.,t=0.):
"""
NAME:
_Rforce
PURPOSE:
evaluate the radial force for this potential
INPUT:
R - Galactocentric cylindrical radius
phi - azimuth
t - time
OUTPUT:
the radial force
HISTORY:
2011-10-19 - Written - Bovy (IAS)
"""
if R < self._rb:
return -self._p*self._mphio/self._m*self._rb2p/R**(self._p+1.)\
*numpy.cos(self._m*phi-self._mphib)
else:
return -self._p*self._mphio/self._m*R**(self._p-1.)\
*numpy.cos(self._m*phi-self._mphib)
def _phiforce(self,R,phi=0.,t=0.):
"""
NAME:
_phiforce
PURPOSE:
evaluate the azimuthal force for this potential
INPUT:
R - Galactocentric cylindrical radius
phi - azimuth
t - time
OUTPUT:
the azimuthal force
HISTORY:
2011-10-19 - Written - Bovy (IAS)
"""
if R < self._rb:
return self._mphio*numpy.sin(self._m*phi-self._mphib)\
*self._rbp*(2.*self._r1p-self._rbp/R**self._p)
else:
return self._mphio*R**self._p*numpy.sin(self._m*phi-self._mphib)
def _R2deriv(self,R,phi=0.,t=0.):
if R < self._rb:
return -self._p*(self._p+1.)*self._mphio/self._m\
*self._rb2p/R**(self._p+2.)*numpy.cos(self._m*phi-self._mphib)
else:
return self._p*(self._p-1.)/self._m*self._mphio*R**(self._p-2.)\
*numpy.cos(self._m*phi-self._mphib)
def _phi2deriv(self,R,phi=0.,t=0.):
if R < self._rb:
return -self._m*self._mphio*numpy.cos(self._m*phi-self._mphib)\
*self._rbp*(2.*self._r1p-self._rbp/R**self._p)
else:
return -self._m*self._mphio*R**self._p\
*numpy.cos(self._m*phi-self._mphib)
def _Rphideriv(self,R,phi=0.,t=0.):
if R < self._rb:
return -self._p*self._mphio/self._m*self._rb2p/R**(self._p+1.)\
*numpy.sin(self._m*phi-self._mphib)
else:
return -self._p*self._mphio*R**(self._p-1.)*\
numpy.sin(self._m*phi-self._mphib)
[docs]class LopsidedDiskPotential(CosmphiDiskPotential):
"""Class that implements the disk potential
.. math::
\\Phi(R,\\phi) = \\mathrm{amp}\\,\\phi_0\\,\\left(\\frac{R}{R_1}\\right)^p\\,\\cos\\left(\\phi-\\phi_b\\right)
Special case of CosmphiDiskPotential with m=1; see documentation for CosmphiDiskPotential
"""
[docs] def __init__(self,amp=1.,phib=25.*_degtorad,
p=1.,phio=0.01,r1=1.,
cp=None,sp=None,ro=None,vo=None):
CosmphiDiskPotential.__init__(self,
amp=amp,
phib=phib,
p=p,phio=phio,m=1.,
cp=cp,sp=sp,ro=ro,vo=vo)
self.hasC= True
self.hasC_dxdv= True