Source code for galpy.potential.EllipticalDiskPotential

###############################################################################
#   EllipticalDiskPotential: Kuijken & Tremaine (1994)'s elliptical disk 
#   potential
###############################################################################
import numpy
from ..util import conversion
from .planarPotential import planarPotential
_degtorad= numpy.pi/180.
[docs]class EllipticalDiskPotential(planarPotential): """Class that implements the Elliptical disk potential of Kuijken & Tremaine (1994) .. math:: \\Phi(R,\\phi) = \\mathrm{amp}\\,\\phi_0\\,\\left(\\frac{R}{R_1}\\right)^p\\,\\cos\\left(2\\,(\\phi-\\phi_b)\\right) This potential can be grown between :math:`t_{\mathrm{form}}` and :math:`t_{\mathrm{form}}+T_{\mathrm{steady}}` in a similar way as DehnenBarPotential, but times are given directly in galpy time units """
[docs] def __init__(self,amp=1.,phib=25.*_degtorad, p=1.,twophio=0.01,r1=1., tform=None,tsteady=None, cp=None,sp=None,ro=None,vo=None): """ NAME: __init__ PURPOSE: initialize an Elliptical disk potential phi(R,phi) = phio (R/Ro)^p cos[2(phi-phib)] INPUT: amp= amplitude to be applied to the potential (default: 1.), see twophio below tform= start of growth (to smoothly grow this potential (can be Quantity) tsteady= time delay at which the perturbation is fully grown (default: 2.; can be Quantity) p= power-law index of the phi(R) = (R/Ro)^p part r1= (1.) normalization radius for the amplitude (can be Quantity) Either: a) phib= angle (in rad; default=25 degree; or can be Quantity) twophio= potential perturbation (in terms of 2phio/vo^2 if vo=1 at Ro=1; can be Quantity with units of velocity-squared) b) cp, sp= twophio * cos(2phib), twophio * sin(2phib) (can be Quantity with units of velocity-squared) OUTPUT: (none) HISTORY: 2011-10-19 - Started - Bovy (IAS) """ planarPotential.__init__(self,amp=amp,ro=ro,vo=vo) phib= conversion.parse_angle(phib) r1= conversion.parse_length(r1,ro=self._ro) tform= conversion.parse_time(tform,ro=self._ro,vo=self._vo) tsteady= conversion.parse_time(tsteady,ro=self._ro,vo=self._vo) twophio= conversion.parse_energy(twophio,vo=self._vo) cp= conversion.parse_energy(cp,vo=self._vo) sp= conversion.parse_energy(sp,vo=self._vo) # Back to old definition self._amp/= r1**p self.hasC= True self.hasC_dxdv= True if cp is None or sp is None: self._phib= phib self._twophio= twophio else: self._twophio= numpy.sqrt(cp*cp+sp*sp) self._phib= numpy.arctan2(sp,cp)/2. self._p= p if not tform is None: self._tform= tform else: self._tform= None if not tsteady is None: self._tsteady= self._tform+tsteady else: if self._tform is None: self._tsteady= None else: self._tsteady= self._tform+2.
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) """ #Calculate relevant time if not self._tform is None: 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: #fully on smooth= 1. else: smooth= 1. return smooth*self._twophio/2.*R**self._p\ *numpy.cos(2.*(phi-self._phib)) 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) """ #Calculate relevant time if not self._tform is None: 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: #fully on smooth= 1. else: smooth= 1. return -smooth*self._p*self._twophio/2.*R**(self._p-1.)\ *numpy.cos(2.*(phi-self._phib)) 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) """ #Calculate relevant time if not self._tform is None: 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: #fully on smooth= 1. else: smooth= 1. return smooth*self._twophio*R**self._p*numpy.sin(2.*(phi-self._phib)) def _R2deriv(self,R,phi=0.,t=0.): #Calculate relevant time if not self._tform is None: 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: #fully on smooth= 1. else: smooth= 1. return smooth*self._p*(self._p-1.)/2.*self._twophio*R**(self._p-2.)\ *numpy.cos(2.*(phi-self._phib)) def _phi2deriv(self,R,phi=0.,t=0.): #Calculate relevant time if not self._tform is None: 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: #perturbation is fully on smooth= 1. else: smooth= 1. return -2.*smooth*self._twophio*R**self._p*numpy.cos(2.*(phi-self._phib)) def _Rphideriv(self,R,phi=0.,t=0.): #Calculate relevant time if not self._tform is None: 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: #perturbation is fully on smooth= 1. else: smooth= 1. return -smooth*self._p*self._twophio*R**(self._p-1.)*numpy.sin(2.*(phi-self._phib)) def tform(self): #pragma: no cover """ NAME: tform PURPOSE: return formation time of the perturbation INPUT: (none) OUTPUT: tform in normalized units HISTORY: 2011-10-19 - Written - Bovy (IAS) """ return self._tform