# Source code for galpy.potential.EllipticalDiskPotential

###############################################################################
#   EllipticalDiskPotential: Kuijken & Tremaine (1994)'s elliptical disk
#   potential
###############################################################################
import numpy
from ..util import conversion
from .planarPotential import planarPotential
[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

"""
p=1.,twophio=0.01,r1=1.,
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)
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
else:
if self._tform is None: self._tsteady= None

def _evaluate(self,R,phi=0.,t=0.):
"""
NAME:
_evaluate
PURPOSE:
evaluate the potential at R,phi,t
INPUT:
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.
deltat= t-self._tform
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:
phi - azimuth
t - time
OUTPUT:
HISTORY:
2011-10-19 - Written - Bovy (IAS)
"""
#Calculate relevant time
if not self._tform is None:
if t < self._tform:
smooth= 0.
deltat= t-self._tform
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:
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.
deltat= t-self._tform
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.
deltat= t-self._tform
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.
deltat= t-self._tform
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.
deltat= t-self._tform
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