from __future__ import division, print_function
import os
import copy
import pickle
import numpy
from scipy import integrate
from ..util import bovy_plot as plot
from ..util import config
from ..util.bovy_conversion import physical_conversion,\
potential_physical_input, freq_in_Gyr, physical_compatible
from .Potential import Potential, PotentialError, lindbladR, flatten
from .DissipativeForce import _isDissipative
from .plotRotcurve import plotRotcurve
from .plotEscapecurve import _INF, plotEscapecurve
_APY_LOADED= True
try:
from astropy import units
except ImportError:
_APY_LOADED= False
class planarPotential(object):
"""Class representing 2D (R,\phi) potentials"""
def __init__(self,amp=1.,ro=None,vo=None):
self._amp= amp
self.dim= 2
self.isNonAxi= True #Gets reset by planarAxiPotential
self.isRZ= False
self.hasC= False
self.hasC_dxdv= False
self.hasC_dens= False
# Parse ro and vo
if ro is None:
self._ro= config.__config__.getfloat('normalization','ro')
self._roSet= False
else:
if _APY_LOADED and isinstance(ro,units.Quantity):
ro= ro.to(units.kpc).value
self._ro= ro
self._roSet= True
if vo is None:
self._vo= config.__config__.getfloat('normalization','vo')
self._voSet= False
else:
if _APY_LOADED and isinstance(vo,units.Quantity):
vo= vo.to(units.km/units.s).value
self._vo= vo
self._voSet= True
return None
[docs] def __mul__(self,b):
"""
NAME:
__mul__
PURPOSE:
Multiply a planarPotential's amplitude by a number
INPUT:
b - number
OUTPUT:
New instance with amplitude = (old amplitude) x b
HISTORY:
2019-01-27 - Written - Bovy (UofT)
"""
if not isinstance(b,(int,float)):
raise TypeError("Can only multiply a planarPotential instance with a number")
out= copy.deepcopy(self)
out._amp*= b
return out
# Similar functions
__rmul__= __mul__
def __div__(self,b): return self.__mul__(1./b)
__truediv__= __div__
[docs] def __add__(self,b):
"""
NAME:
__add__
PURPOSE:
Add planarPotential instances together to create a multi-component potential (e.g., pot= pot1+pot2+pot3)
INPUT:
b - planarPotential instance or a list thereof
OUTPUT:
List of planarPotential instances that represents the combined potential
HISTORY:
2019-01-27 - Written - Bovy (UofT)
"""
from ..potential import flatten as flatten_pot
if not isinstance(flatten_pot([b])[0],(Potential,planarPotential)):
raise TypeError("""Can only combine galpy Potential"""
"""/planarPotential objects with """
"""other such objects or lists thereof""")
assert physical_compatible(self,b), \
"""Physical unit conversion parameters (ro,vo) are not """\
"""compatible between potentials to be combined"""
if isinstance(b,list):
return [self]+b
else:
return [self,b]
# Define separately to keep order
def __radd__(self,b):
from ..potential import flatten as flatten_pot
if not isinstance(flatten_pot([b])[0],(Potential,planarPotential)):
raise TypeError("""Can only combine galpy Force objects with """
"""other Force objects or lists thereof""")
assert physical_compatible(self,b), \
"""Physical unit conversion parameters (ro,vo) are not """\
"""compatible between potentials to be combined"""
# If we get here, b has to be a list
return b+[self]
[docs] def turn_physical_off(self):
"""
NAME:
turn_physical_off
PURPOSE:
turn off automatic returning of outputs in physical units
INPUT:
(none)
OUTPUT:
(none)
HISTORY:
2016-01-30 - Written - Bovy (UofT)
"""
self._roSet= False
self._voSet= False
return None
[docs] def turn_physical_on(self,ro=None,vo=None):
"""
NAME:
turn_physical_on
PURPOSE:
turn on automatic returning of outputs in physical units
INPUT:
ro= reference distance (kpc; can be Quantity)
vo= reference velocity (km/s; can be Quantity)
OUTPUT:
(none)
HISTORY:
2016-01-30 - Written - Bovy (UofT)
2020-04-22 - Don't turn on a parameter when it is False - Bovy (UofT)
"""
if not ro is False: self._roSet= True
if not vo is False: self._voSet= True
if not ro is None and ro:
if _APY_LOADED and isinstance(ro,units.Quantity):
ro= ro.to(units.kpc).value
self._ro= ro
if not vo is None and vo:
if _APY_LOADED and isinstance(vo,units.Quantity):
vo= vo.to(units.km/units.s).value
self._vo= vo
return None
[docs] @potential_physical_input
@physical_conversion('energy',pop=True)
def __call__(self,R,phi=0.,t=0.,dR=0,dphi=0):
"""
NAME:
__call__
PURPOSE:
evaluate the potential
INPUT:
R - Cylindrica radius (can be Quantity)
phi= azimuth (optional; can be Quantity)
t= time (optional; can be Quantity)
OUTPUT:
Phi(R(,phi,t)))
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return self._call_nodecorator(R,phi=phi,t=t,dR=dR,dphi=dphi)
def _call_nodecorator(self,R,phi=0.,t=0.,dR=0,dphi=0):
# Separate, so it can be used during orbit integration
if dR == 0 and dphi == 0:
try:
return self._amp*self._evaluate(R,phi=phi,t=t)
except AttributeError: #pragma: no cover
raise PotentialError("'_evaluate' function not implemented for this potential")
elif dR == 1 and dphi == 0:
return -self.Rforce(R,phi=phi,t=t,use_physical=False)
elif dR == 0 and dphi == 1:
return -self.phiforce(R,phi=phi,t=t,use_physical=False)
elif dR == 2 and dphi == 0:
return self.R2deriv(R,phi=phi,t=t,use_physical=False)
elif dR == 0 and dphi == 2:
return self.phi2deriv(R,phi=phi,t=t,use_physical=False)
elif dR == 1 and dphi == 1:
return self.Rphideriv(R,phi=phi,t=t,use_physical=False)
elif dR != 0 or dphi != 0:
raise NotImplementedError('Higher-order derivatives not implemented for this potential')
[docs] @potential_physical_input
@physical_conversion('force',pop=True)
def Rforce(self,R,phi=0.,t=0.):
"""
NAME:
Rforce
PURPOSE:
evaluate the radial force
INPUT:
R - Cylindrical radius (can be Quantity)
phi= azimuth (optional; can be Quantity)
t= time (optional; can be Quantity)
OUTPUT:
F_R(R,(\phi,t)))
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return self._Rforce_nodecorator(R,phi=phi,t=t)
def _Rforce_nodecorator(self,R,phi=0.,t=0.):
# Separate, so it can be used during orbit integration
try:
return self._amp*self._Rforce(R,phi=phi,t=t)
except AttributeError: #pragma: no cover
raise PotentialError("'_Rforce' function not implemented for this potential")
[docs] @potential_physical_input
@physical_conversion('force',pop=True)
def phiforce(self,R,phi=0.,t=0.):
"""
NAME:
phiforce
PURPOSE:
evaluate the phi force
INPUT:
R - Cylindrical radius (can be Quantity)
phi= azimuth (optional; can be Quantity)
t= time (optional; can be Quantity)
OUTPUT:
F_phi(R,(phi,t)))
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return self._phiforce_nodecorator(R,phi=phi,t=t)
def _phiforce_nodecorator(self,R,phi=0.,t=0.):
# Separate, so it can be used during orbit integration
try:
return self._amp*self._phiforce(R,phi=phi,t=t)
except AttributeError: #pragma: no cover
raise PotentialError("'_phiforce' function not implemented for this potential")
@potential_physical_input
@physical_conversion('forcederivative',pop=True)
def R2deriv(self,R,phi=0.,t=0.):
"""
NAME:
R2deriv
PURPOSE:
evaluate the second radial derivative
INPUT:
R - Cylindrical radius (can be Quantity)
phi= azimuth (optional; can be Quantity)
t= time (optional; can be Quantity)
OUTPUT:
d2phi/dR2
HISTORY:
2011-10-09 - Written - Bovy (IAS)
"""
try:
return self._amp*self._R2deriv(R,phi=phi,t=t)
except AttributeError: #pragma: no cover
raise PotentialError("'_R2deriv' function not implemented for this potential")
@potential_physical_input
@physical_conversion('forcederivative',pop=True)
def phi2deriv(self,R,phi=0.,t=0.):
"""
NAME:
phi2deriv
PURPOSE:
evaluate the second azimuthal derivative
INPUT:
R - Cylindrical radius (can be Quantity)
phi= azimuth (optional; can be Quantity)
t= time (optional; can be Quantity)
OUTPUT:
d2phi/daz2
HISTORY:
2014-04-06 - Written - Bovy (IAS)
"""
try:
return self._amp*self._phi2deriv(R,phi=phi,t=t)
except AttributeError: #pragma: no cover
raise PotentialError("'_phi2deriv' function not implemented for this potential")
@potential_physical_input
@physical_conversion('forcederivative',pop=True)
def Rphideriv(self,R,phi=0.,t=0.):
"""
NAME:
Rphideriv
PURPOSE:
evaluate the mixed radial and azimuthal derivative
INPUT:
R - Cylindrical radius (can be Quantity)
phi= azimuth (optional can be Quantity)
t= time (optional; can be Quantity)
OUTPUT:
d2phi/dR d az
HISTORY:
2014-05-21 - Written - Bovy (IAS)
"""
try:
return self._amp*self._Rphideriv(R,phi=phi,t=t)
except AttributeError: #pragma: no cover
raise PotentialError("'_Rphideriv' function not implemented for this potential")
def plot(self,*args,**kwargs):
"""
NAME:
plot
PURPOSE:
plot the potential
INPUT:
Rrange - range (can be Quantity)
grid - number of points to plot
savefilename - save to or restore from this savefile (pickle)
+bovy_plot(*args,**kwargs)
OUTPUT:
plot to output device
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return plotplanarPotentials(self,*args,**kwargs)
class planarAxiPotential(planarPotential):
"""Class representing axisymmetric planar potentials"""
def __init__(self,amp=1.,ro=None,vo=None):
planarPotential.__init__(self,amp=amp,ro=ro,vo=vo)
self.isNonAxi= False
return None
def _phiforce(self,R,phi=0.,t=0.):
return 0.
def _phi2deriv(self,R,phi=0.,t=0.): #pragma: no cover
"""
NAME:
_phi2deriv
PURPOSE:
evaluate the second azimuthal derivative for this potential
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
OUTPUT:
the second azimuthal derivative
HISTORY:
2011-10-17 - Written - Bovy (IAS)
"""
return 0.
def _Rphideriv(self,R,phi=0.,t=0.): #pragma: no cover
"""
NAME:
_Rphideriv
PURPOSE:
evaluate the radial+azimuthal derivative for this potential
INPUT:
R - Galactocentric cylindrical radius
z - vertical height
phi - azimuth
t - time
OUTPUT:
the radial+azimuthal derivative
HISTORY:
2011-10-17 - Written - Bovy (IAS)
"""
return 0.
@potential_physical_input
@physical_conversion('velocity',pop=True)
def vcirc(self,R,phi=None,t=0.):
"""
NAME:
vcirc
PURPOSE:
calculate the circular velocity at R in potential Pot
INPUT:
Pot - Potential instance or list of such instances
R - Galactocentric radius (can be Quantity)
phi= (None) azimuth to use for non-axisymmetric potentials
t - time (optional; can be Quantity)
OUTPUT:
circular rotation velocity
HISTORY:
2011-10-09 - Written - Bovy (IAS)
2016-06-15 - Added phi= keyword for non-axisymmetric potential - Bovy (UofT)
"""
return numpy.sqrt(R*-self.Rforce(R,phi=phi,t=t,use_physical=False))
@potential_physical_input
@physical_conversion('frequency',pop=True)
def omegac(self,R,t=0.):
"""
NAME:
omegac
PURPOSE:
calculate the circular angular speed at R in potential Pot
INPUT:
Pot - Potential instance or list of such instances
R - Galactocentric radius (can be Quantity)
t - time (optional; can be Quantity)
OUTPUT:
circular angular speed
HISTORY:
2011-10-09 - Written - Bovy (IAS)
"""
return numpy.sqrt(-self.Rforce(R,t=t,use_physical=False)/R)
@potential_physical_input
@physical_conversion('frequency',pop=True)
def epifreq(self,R,t=0.):
"""
NAME:
epifreq
PURPOSE:
calculate the epicycle frequency at R in this potential
INPUT:
R - Galactocentric radius (can be Quantity)
t - time (optional; can be Quantity)
OUTPUT:
epicycle frequency
HISTORY:
2011-10-09 - Written - Bovy (IAS)
"""
return numpy.sqrt(self.R2deriv(R,t=t,use_physical=False)
-3./R*self.Rforce(R,t=t,use_physical=False))
@physical_conversion('position',pop=True)
def lindbladR(self,OmegaP,m=2,t=0.,**kwargs):
"""
NAME:
lindbladR
PURPOSE:
calculate the radius of a Lindblad resonance
INPUT:
OmegaP - pattern speed (can be Quantity)
m= order of the resonance (as in m(O-Op)=kappa (negative m for outer)
use m='corotation' for corotation
+scipy.optimize.brentq xtol,rtol,maxiter kwargs
t - time (optional; can be Quantity)
OUTPUT:
radius of Linblad resonance, None if there is no resonance
HISTORY:
2011-10-09 - Written - Bovy (IAS)
"""
if _APY_LOADED and isinstance(OmegaP,units.Quantity):
OmegaP= OmegaP.to(1/units.Gyr).value/freq_in_Gyr(self._vo,self._ro)
return lindbladR(self,OmegaP,m=m,t=t,use_physical=False,**kwargs)
@potential_physical_input
@physical_conversion('velocity',pop=True)
def vesc(self,R,t=0.):
"""
NAME:
vesc
PURPOSE:
calculate the escape velocity at R for potential Pot
INPUT:
Pot - Potential instances or list thereof
R - Galactocentric radius (can be Quantity)
t - time (optional; can be Quantity)
OUTPUT:
escape velocity
HISTORY:
2011-10-09 - Written - Bovy (IAS)
"""
return numpy.sqrt(2.*(self(_INF,t=t,use_physical=False)
-self(R,t=t,use_physical=False)))
[docs] def plotRotcurve(self,*args,**kwargs):
"""
NAME:
plotRotcurve
PURPOSE:
plot the rotation curve for this potential
INPUT:
Rrange - range (can be Quantity)
grid - number of points to plot
savefilename - save to or restore from this savefile (pickle)
+bovy_plot(*args,**kwargs)
OUTPUT:
plot to output device
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return plotRotcurve(self,*args,**kwargs)
[docs] def plotEscapecurve(self,*args,**kwargs):
"""
NAME:
plotEscapecurve
PURPOSE:
plot the escape velocity curve for this potential
INPUT:
Rrange - range (can be Quantity)
grid - number of points to plot
savefilename - save to or restore from this savefile (pickle)
+bovy_plot(*args,**kwargs)
OUTPUT:
plot to output device
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return plotEscapecurve(self,*args,**kwargs)
class planarPotentialFromRZPotential(planarAxiPotential):
"""Class that represents an axisymmetic planar potential derived from a
RZPotential"""
def __init__(self,RZPot):
"""
NAME:
__init__
PURPOSE:
Initialize
INPUT:
RZPot - RZPotential instance
OUTPUT:
planarAxiPotential instance
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
planarAxiPotential.__init__(self,amp=1.,ro=RZPot._ro,vo=RZPot._vo)
# Also transfer roSet and voSet
self._roSet= RZPot._roSet
self._voSet= RZPot._voSet
self._Pot= RZPot
self.hasC= RZPot.hasC
self.hasC_dxdv= RZPot.hasC_dxdv
self.hasC_dens= RZPot.hasC_dens
return None
def _evaluate(self,R,phi=0.,t=0.):
"""
NAME:
_evaluate
PURPOSE:
evaluate the potential
INPUT:
R
phi
t
OUTPUT:
Pot(R(,\phi,t))
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return self._Pot(R,0.,t=t,use_physical=False)
def _Rforce(self,R,phi=0.,t=0.):
"""
NAME:
_Rforce
PURPOSE:
evaluate the radial force
INPUT:
R
phi
t
OUTPUT:
F_R(R(,\phi,t))
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return self._Pot.Rforce(R,0.,t=t,use_physical=False)
def _R2deriv(self,R,phi=0.,t=0.):
"""
NAME:
_R2deriv
PURPOSE:
evaluate the second radial derivative
INPUT:
R
phi
t
OUTPUT:
d2phi/dR2
HISTORY:
2011-10-09 - Written - Bovy (IAS)
"""
return self._Pot.R2deriv(R,0.,t=t,use_physical=False)
[docs]def RZToplanarPotential(RZPot):
"""
NAME:
RZToplanarPotential
PURPOSE:
convert an RZPotential to a planarPotential in the mid-plane (z=0)
INPUT:
RZPot - RZPotential instance or list of such instances (existing planarPotential instances are just copied to the output)
OUTPUT:
planarPotential instance(s)
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
RZPot= flatten(RZPot)
if _isDissipative(RZPot):
raise NotImplementedError("Converting dissipative forces to 2D potentials is currently not supported")
if isinstance(RZPot,list):
out= []
for pot in RZPot:
if isinstance(pot,planarPotential) and not pot.isNonAxi:
out.append(pot)
elif isinstance(pot,Potential) and not pot.isNonAxi:
out.append(planarPotentialFromRZPotential(pot))
else:
raise PotentialError("Input to 'RZToplanarPotential' is neither an RZPotential-instance or a list of such instances")
return out
elif isinstance(RZPot,Potential) and not RZPot.isNonAxi:
return planarPotentialFromRZPotential(RZPot)
elif isinstance(RZPot,planarPotential) and not RZPot.isNonAxi:
return RZPot
else:
raise PotentialError("Input to 'RZToplanarPotential' is neither an RZPotential-instance or a list of such instances")
class planarPotentialFromFullPotential(planarPotential):
"""Class that represents a planar potential derived from a non-axisymmetric
3D potential"""
def __init__(self,Pot):
"""
NAME:
__init__
PURPOSE:
Initialize
INPUT:
Pot - Potential instance
OUTPUT:
planarPotential instance
HISTORY:
2016-06-02 - Written - Bovy (UofT)
"""
planarPotential.__init__(self,amp=1.,ro=Pot._ro,vo=Pot._vo)
# Also transfer roSet and voSet
self._roSet= Pot._roSet
self._voSet= Pot._voSet
self._Pot= Pot
self.hasC= Pot.hasC
self.hasC_dxdv= Pot.hasC_dxdv
self.hasC_dens= Pot.hasC_dens
return None
def _evaluate(self,R,phi=0.,t=0.):
"""
NAME:
_evaluate
PURPOSE:
evaluate the potential
INPUT:
R
phi
t
OUTPUT:
Pot(R(,\phi,t))
HISTORY:
2016-06-02 - Written - Bovy (UofT)
"""
return self._Pot(R,0.,phi=phi,t=t,use_physical=False)
def _Rforce(self,R,phi=0.,t=0.):
"""
NAME:
_Rforce
PURPOSE:
evaluate the radial force
INPUT:
R
phi
t
OUTPUT:
F_R(R(,\phi,t))
HISTORY:
2016-06-02 - Written - Bovy (UofT)
"""
return self._Pot.Rforce(R,0.,phi=phi,t=t,use_physical=False)
def _phiforce(self,R,phi=0.,t=0.):
"""
NAME:
_phiforce
PURPOSE:
evaluate the azimuthal force
INPUT:
R
phi
t
OUTPUT:
F_phi(R(,\phi,t))
HISTORY:
2016-06-02 - Written - Bovy (UofT)
"""
return self._Pot.phiforce(R,0.,phi=phi,t=t,use_physical=False)
def _R2deriv(self,R,phi=0.,t=0.):
"""
NAME:
_R2deriv
PURPOSE:
evaluate the second radial derivative
INPUT:
R
phi
t
OUTPUT:
d2phi/dR2
HISTORY:
2016-06-02 - Written - Bovy (UofT)
"""
return self._Pot.R2deriv(R,0.,phi=phi,t=t,use_physical=False)
def _phi2deriv(self,R,phi=0.,t=0.):
"""
NAME:
_phi2deriv
PURPOSE:
evaluate the second azimuthal derivative
INPUT:
R
phi
t
OUTPUT:
d2phi/dphi2
HISTORY:
2016-06-02 - Written - Bovy (UofT)
"""
return self._Pot.phi2deriv(R,0.,phi=phi,t=t,use_physical=False)
def _Rphideriv(self,R,phi=0.,t=0.):
"""
NAME:
_Rphideriv
PURPOSE:
evaluate the mixed radial-azimuthal derivative
INPUT:
R
phi
t
OUTPUT:
d2phi/dRdphi
HISTORY:
2016-06-02 - Written - Bovy (UofT)
"""
return self._Pot.Rphideriv(R,0.,phi=phi,t=t,use_physical=False)
def OmegaP(self):
"""
NAME:
OmegaP
PURPOSE:
return the pattern speed
INPUT:
(none)
OUTPUT:
pattern speed
HISTORY:
2016-05-31 - Written - Bovy (UofT)
"""
return self._Pot.OmegaP()
[docs]def toPlanarPotential(Pot):
"""
NAME:
toPlanarPotential
PURPOSE:
convert an Potential to a planarPotential in the mid-plane (z=0)
INPUT:
Pot - Potential instance or list of such instances (existing planarPotential instances are just copied to the output)
OUTPUT:
planarPotential instance(s)
HISTORY:
2016-06-11 - Written - Bovy (UofT)
"""
Pot= flatten(Pot)
if _isDissipative(Pot):
raise NotImplementedError("Converting dissipative forces to 2D potentials is currently not supported")
elif isinstance(Pot,list):
out= []
for pot in Pot:
if isinstance(pot,planarPotential):
out.append(pot)
elif isinstance(pot,Potential) and pot.isNonAxi:
out.append(planarPotentialFromFullPotential(pot))
elif isinstance(pot,Potential):
out.append(planarPotentialFromRZPotential(pot))
else:
raise PotentialError("Input to 'toPlanarPotential' is neither an Potential-instance or a list of such instances")
return out
elif isinstance(Pot,Potential) and Pot.isNonAxi:
return planarPotentialFromFullPotential(Pot)
elif isinstance(Pot,Potential):
return planarPotentialFromRZPotential(Pot)
elif isinstance(Pot,planarPotential):
return Pot
else:
raise PotentialError("Input to 'toPlanarPotential' is neither an Potential-instance or a list of such instances")
[docs]@potential_physical_input
@physical_conversion('energy',pop=True)
def evaluateplanarPotentials(Pot,R,phi=None,t=0.,dR=0,dphi=0):
"""
NAME:
evaluateplanarPotentials
PURPOSE:
evaluate a (list of) planarPotential instance(s)
INPUT:
Pot - (list of) planarPotential instance(s)
R - Cylindrical radius (can be Quantity)
phi= azimuth (optional; can be Quantity)
t= time (optional; can be Quantity)
dR=, dphi= if set to non-zero integers, return the dR,dphi't derivative instead
OUTPUT:
Phi(R(,phi,t))
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return _evaluateplanarPotentials(Pot,R,phi=phi,t=t,dR=dR,dphi=dphi)
def _evaluateplanarPotentials(Pot,R,phi=None,t=0.,dR=0,dphi=0):
from .Potential import _isNonAxi
isList= isinstance(Pot,list)
nonAxi= _isNonAxi(Pot)
if nonAxi and phi is None:
raise PotentialError("The (list of) planarPotential instances is non-axisymmetric, but you did not provide phi")
if isList and numpy.all([isinstance(p,planarPotential) for p in Pot]):
sum= 0.
for pot in Pot:
if nonAxi:
sum+= pot._call_nodecorator(R,phi=phi,t=t,dR=dR,dphi=dphi)
else:
sum+= pot._call_nodecorator(R,t=t,dR=dR,dphi=dphi)
return sum
elif isinstance(Pot,planarPotential):
if nonAxi:
return Pot._call_nodecorator(R,phi=phi,t=t,dR=dR,dphi=dphi)
else:
return Pot._call_nodecorator(R,t=t,dR=dR,dphi=dphi)
else: #pragma: no cover
raise PotentialError("Input to 'evaluatePotentials' is neither a Potential-instance or a list of such instances")
[docs]@potential_physical_input
@physical_conversion('force',pop=True)
def evaluateplanarRforces(Pot,R,phi=None,t=0.):
"""
NAME:
evaluateplanarRforces
PURPOSE:
evaluate the Rforce of a (list of) planarPotential instance(s)
INPUT:
Pot - (list of) planarPotential instance(s)
R - Cylindrical radius (can be Quantity)
phi= azimuth (optional can be Quantity)
t= time (optional; can be Quantity)
OUTPUT:
F_R(R(,phi,t))
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return _evaluateplanarRforces(Pot,R,phi=phi,t=t)
def _evaluateplanarRforces(Pot,R,phi=None,t=0.):
"""Raw, undecorated function for internal use"""
from .Potential import _isNonAxi
isList= isinstance(Pot,list)
nonAxi= _isNonAxi(Pot)
if nonAxi and phi is None:
raise PotentialError("The (list of) planarPotential instances is non-axisymmetric, but you did not provide phi")
if isinstance(Pot,list) \
and numpy.all([isinstance(p,planarPotential) for p in Pot]):
sum= 0.
for pot in Pot:
if nonAxi:
sum+= pot._Rforce_nodecorator(R,phi=phi,t=t)
else:
sum+= pot._Rforce_nodecorator(R,t=t)
return sum
elif isinstance(Pot,planarPotential):
if nonAxi:
return Pot._Rforce_nodecorator(R,phi=phi,t=t)
else:
return Pot._Rforce_nodecorator(R,t=t)
else: #pragma: no cover
raise PotentialError("Input to 'evaluatePotentials' is neither a Potential-instance or a list of such instances")
[docs]@potential_physical_input
@physical_conversion('force',pop=True)
def evaluateplanarphiforces(Pot,R,phi=None,t=0.):
"""
NAME:
evaluateplanarphiforces
PURPOSE:
evaluate the phiforce of a (list of) planarPotential instance(s)
INPUT:
Pot - (list of) planarPotential instance(s)
R - Cylindrical radius (can be Quantity)
phi= azimuth (optional; can be Quantity)
t= time (optional; can be Quantity)
OUTPUT:
F_phi(R(,phi,t))
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
return _evaluateplanarphiforces(Pot,R,phi=phi,t=t)
def _evaluateplanarphiforces(Pot,R,phi=None,t=0.):
from .Potential import _isNonAxi
isList= isinstance(Pot,list)
nonAxi= _isNonAxi(Pot)
if nonAxi and phi is None:
raise PotentialError("The (list of) planarPotential instances is non-axisymmetric, but you did not provide phi")
if isinstance(Pot,list) \
and numpy.all([isinstance(p,planarPotential) for p in Pot]):
sum= 0.
for pot in Pot:
if nonAxi:
sum+= pot._phiforce_nodecorator(R,phi=phi,t=t)
else:
sum+= pot._phiforce_nodecorator(R,t=t)
return sum
elif isinstance(Pot,planarPotential):
if nonAxi:
return Pot._phiforce_nodecorator(R,phi=phi,t=t)
else:
return Pot._phiforce_nodecorator(R,t=t)
else: #pragma: no cover
raise PotentialError("Input to 'evaluatePotentials' is neither a Potential-instance or a list of such instances")
[docs]@potential_physical_input
@physical_conversion('forcederivative',pop=True)
def evaluateplanarR2derivs(Pot,R,phi=None,t=0.):
"""
NAME:
evaluateplanarR2derivs
PURPOSE:
evaluate the second radial derivative of a (list of) planarPotential instance(s)
INPUT:
Pot - (list of) planarPotential instance(s)
R - Cylindrical radius (can be Quantity)
phi= azimuth (optional; can be Quantity)
t= time (optional; can be Quantity)
OUTPUT:
F_R(R(,phi,t))
HISTORY:
2010-10-09 - Written - Bovy (IAS)
"""
from .Potential import _isNonAxi
isList= isinstance(Pot,list)
nonAxi= _isNonAxi(Pot)
if nonAxi and phi is None:
raise PotentialError("The (list of) planarPotential instances is non-axisymmetric, but you did not provide phi")
if isinstance(Pot,list) \
and numpy.all([isinstance(p,planarPotential) for p in Pot]):
sum= 0.
for pot in Pot:
if nonAxi:
sum+= pot.R2deriv(R,phi=phi,t=t,use_physical=False)
else:
sum+= pot.R2deriv(R,t=t,use_physical=False)
return sum
elif isinstance(Pot,planarPotential):
if nonAxi:
return Pot.R2deriv(R,phi=phi,t=t,use_physical=False)
else:
return Pot.R2deriv(R,t=t,use_physical=False)
else: #pragma: no cover
raise PotentialError("Input to 'evaluatePotentials' is neither a Potential-instance or a list of such instances")
[docs]def LinShuReductionFactor(axiPot,R,sigmar,nonaxiPot=None,
k=None,m=None,OmegaP=None):
"""
NAME:
LinShuReductionFactor
PURPOSE:
Calculate the Lin & Shu (1966) reduction factor: the reduced linear response of a kinematically-warm stellar disk to a perturbation
INPUT:
axiPot - The background, axisymmetric potential
R - Cylindrical radius (can be Quantity)
sigmar - radial velocity dispersion of the population (can be Quantity)
Then either provide:
1) m= m in the perturbation's m x phi (number of arms for a spiral)
k= wavenumber (see Binney & Tremaine 2008)
OmegaP= pattern speed (can be Quantity)
2) nonaxiPot= a non-axisymmetric Potential instance (such as SteadyLogSpiralPotential) that has functions that return OmegaP, m, and wavenumber
OUTPUT:
reduction factor
HISTORY:
2014-08-23 - Written - Bovy (IAS)
"""
axiPot= flatten(axiPot)
from ..potential import omegac, epifreq
if nonaxiPot is None and (OmegaP is None or k is None or m is None):
raise IOError("Need to specify either nonaxiPot= or m=, k=, OmegaP= for LinShuReductionFactor")
elif not nonaxiPot is None:
OmegaP= nonaxiPot.OmegaP()
k= nonaxiPot.wavenumber(R)
m= nonaxiPot.m()
tepif= epifreq(axiPot,R)
s= m*(OmegaP-omegac(axiPot,R))/tepif
chi= sigmar**2.*k**2./tepif**2.
return (1.-s**2.)/numpy.sin(numpy.pi*s)\
*integrate.quad(lambda t: numpy.exp(-chi*(1.+numpy.cos(t)))\
*numpy.sin(s*t)*numpy.sin(t),
0.,numpy.pi)[0]
[docs]def plotplanarPotentials(Pot,*args,**kwargs):
"""
NAME:
plotplanarPotentials
PURPOSE:
plot a planar potential
INPUT:
Rrange - range (can be Quantity)
xrange, yrange - if relevant (can be Quantity)
grid, gridx, gridy - number of points to plot
savefilename - save to or restore from this savefile (pickle)
ncontours - number of contours to plot (if applicable)
+bovy_plot(*args,**kwargs) or bovy_dens2d(**kwargs)
OUTPUT:
plot to output device
HISTORY:
2010-07-13 - Written - Bovy (NYU)
"""
Pot= flatten(Pot)
Rrange= kwargs.pop('Rrange',[0.01,5.])
xrange= kwargs.pop('xrange',[-5.,5.])
yrange= kwargs.pop('yrange',[-5.,5.])
if _APY_LOADED:
if hasattr(Pot,'_ro'):
tro= Pot._ro
else:
tro= Pot[0]._ro
if isinstance(Rrange[0],units.Quantity):
Rrange[0]= Rrange[0].to(units.kpc).value/tro
if isinstance(Rrange[1],units.Quantity):
Rrange[1]= Rrange[1].to(units.kpc).value/tro
if isinstance(xrange[0],units.Quantity):
xrange[0]= xrange[0].to(units.kpc).value/tro
if isinstance(xrange[1],units.Quantity):
xrange[1]= xrange[1].to(units.kpc).value/tro
if isinstance(yrange[0],units.Quantity):
yrange[0]= yrange[0].to(units.kpc).value/tro
if isinstance(yrange[1],units.Quantity):
yrange[1]= yrange[1].to(units.kpc).value/tro
grid= kwargs.pop('grid',100)
gridx= kwargs.pop('gridx',100)
gridy= kwargs.pop('gridy',gridx)
savefilename= kwargs.pop('savefilename',None)
isList= isinstance(Pot,list)
nonAxi= ((isList and Pot[0].isNonAxi) or (not isList and Pot.isNonAxi))
if not savefilename is None and os.path.exists(savefilename):
print("Restoring savefile "+savefilename+" ...")
savefile= open(savefilename,'rb')
potR= pickle.load(savefile)
if nonAxi:
xs= pickle.load(savefile)
ys= pickle.load(savefile)
else:
Rs= pickle.load(savefile)
savefile.close()
else:
if nonAxi:
xs= numpy.linspace(xrange[0],xrange[1],gridx)
ys= numpy.linspace(yrange[0],yrange[1],gridy)
potR= numpy.zeros((gridx,gridy))
for ii in range(gridx):
for jj in range(gridy):
thisR= numpy.sqrt(xs[ii]**2.+ys[jj]**2.)
if xs[ii] >= 0.:
thisphi= numpy.arcsin(ys[jj]/thisR)
else:
thisphi= -numpy.arcsin(ys[jj]/thisR)+numpy.pi
potR[ii,jj]= evaluateplanarPotentials(Pot,thisR,
phi=thisphi,
use_physical=False)
else:
Rs= numpy.linspace(Rrange[0],Rrange[1],grid)
potR= numpy.zeros(grid)
for ii in range(grid):
potR[ii]= evaluateplanarPotentials(Pot,Rs[ii],
use_physical=False)
if not savefilename is None:
print("Writing planar savefile "+savefilename+" ...")
savefile= open(savefilename,'wb')
pickle.dump(potR,savefile)
if nonAxi:
pickle.dump(xs,savefile)
pickle.dump(ys,savefile)
else:
pickle.dump(Rs,savefile)
savefile.close()
if nonAxi:
if not 'orogin' in kwargs:
kwargs['origin']= 'lower'
if not 'cmap' in kwargs:
kwargs['cmap']= 'gist_yarg'
if not 'contours' in kwargs:
kwargs['contours']= True
if not 'xlabel' in kwargs:
kwargs['xlabel']= r"$x / R_0$"
if not 'ylabel' in kwargs:
kwargs['ylabel']= "$y / R_0$"
if not 'aspect' in kwargs:
kwargs['aspect']= 1.
if not 'cntrls' in kwargs:
kwargs['cntrls']= '-'
ncontours= kwargs.pop('ncontours',10)
if not 'levels' in kwargs:
kwargs['levels']= numpy.linspace(numpy.nanmin(potR),numpy.nanmax(potR),ncontours)
return plot.bovy_dens2d(potR.T,
xrange=xrange,
yrange=yrange,**kwargs)
else:
kwargs['xlabel']=r"$R/R_0$"
kwargs['ylabel']=r"$\Phi(R)$"
kwargs['xrange']=Rrange
return plot.bovy_plot(Rs,potR,*args,**kwargs)