Source code for galpy.actionAngle.actionAngleIsochrone

###############################################################################
#   actionAngle: a Python module to calculate  actions, angles, and frequencies
#
#      class: actionAngleIsochrone
#
#             Calculate actions-angle coordinates for the Isochrone potential
#
#      methods:
#             __call__: returns (jr,lz,jz)
#             actionsFreqs: returns (jr,lz,jz,Or,Op,Oz)
#             actionsFreqsAngles: returns (jr,lz,jz,Or,Op,Oz,ar,ap,az)
#
###############################################################################
import copy
import warnings
import numpy
from .actionAngle import actionAngle
from ..potential import IsochronePotential
from ..util import galpyWarning
_APY_LOADED= True
try:
    from astropy import units
except ImportError:
    _APY_LOADED= False
[docs]class actionAngleIsochrone(actionAngle): """Action-angle formalism for the isochrone potential, on the Jphi, Jtheta system of Binney & Tremaine (2008)"""
[docs] def __init__(self,*args,**kwargs): """ NAME: __init__ PURPOSE: initialize an actionAngleIsochrone object INPUT: Either: b= scale parameter of the isochrone parameter (can be Quantity) ip= instance of a IsochronePotential ro= distance from vantage point to GC (kpc; can be Quantity) vo= circular velocity at ro (km/s; can be Quantity) OUTPUT: instance HISTORY: 2013-09-08 - Written - Bovy (IAS) """ actionAngle.__init__(self, ro=kwargs.get('ro',None),vo=kwargs.get('vo',None)) if not 'b' in kwargs and not 'ip' in kwargs: #pragma: no cover raise IOError("Must specify b= for actionAngleIsochrone") if 'ip' in kwargs: ip= kwargs['ip'] if not isinstance(ip,IsochronePotential): #pragma: no cover raise IOError("'Provided ip= does not appear to be an instance of an IsochronePotential") # Check the units self._pot= ip self._check_consistent_units() self.b= ip.b self.amp= ip._amp else: self.b= kwargs['b'] if _APY_LOADED and isinstance(self.b,units.Quantity): self.b= self.b.to(units.kpc).value/self._ro rb= numpy.sqrt(self.b**2.+1.) self.amp= (self.b+rb)**2.*rb self._c= False ext_loaded= False if ext_loaded and (('c' in kwargs and kwargs['c']) or not 'c' in kwargs): #pragma: no cover self._c= True else: self._c= False if not self._c: self._ip= IsochronePotential(amp=self.amp,b=self.b) #Define _pot, because some functions that use actionAngle instances need this self._pot= IsochronePotential(amp=self.amp,b=self.b) # Check the units self._check_consistent_units() return None
def _evaluate(self,*args,**kwargs): """ NAME: __call__ (_evaluate) PURPOSE: evaluate the actions (jr,lz,jz) INPUT: Either: a) R,vR,vT,z,vz[,phi]: 1) floats: phase-space value for single object (phi is optional) (each can be a Quantity) 2) numpy.ndarray: [N] phase-space values for N objects (each can be a Quantity) b) Orbit instance: initial condition used if that's it, orbit(t) if there is a time given as well as the second argument OUTPUT: (jr,lz,jz) HISTORY: 2013-09-08 - Written - Bovy (IAS) """ if len(args) == 5: #R,vR.vT, z, vz R,vR,vT, z, vz= args elif len(args) == 6: #R,vR.vT, z, vz, phi R,vR,vT, z, vz, phi= args else: self._parse_eval_args(*args) R= self._eval_R vR= self._eval_vR vT= self._eval_vT z= self._eval_z vz= self._eval_vz if isinstance(R,float): R= numpy.array([R]) vR= numpy.array([vR]) vT= numpy.array([vT]) z= numpy.array([z]) vz= numpy.array([vz]) if self._c: #pragma: no cover pass else: Lz= R*vT Lx= -z*vT Ly= z*vR-R*vz L2= Lx*Lx+Ly*Ly+Lz*Lz E= self._ip(R,z)+vR**2./2.+vT**2./2.+vz**2./2. L= numpy.sqrt(L2) #Actions Jphi= Lz Jz= L-numpy.fabs(Lz) Jr= self.amp/numpy.sqrt(-2.*E)\ -0.5*(L+numpy.sqrt((L2+4.*self.amp*self.b))) return (Jr,Jphi,Jz) def _actionsFreqs(self,*args,**kwargs): """ NAME: actionsFreqs (_actionsFreqs) PURPOSE: evaluate the actions and frequencies (jr,lz,jz,Omegar,Omegaphi,Omegaz) INPUT: Either: a) R,vR,vT,z,vz[,phi]: 1) floats: phase-space value for single object (phi is optional) (each can be a Quantity) 2) numpy.ndarray: [N] phase-space values for N objects (each can be a Quantity) b) Orbit instance: initial condition used if that's it, orbit(t) if there is a time given as well as the second argument OUTPUT: (jr,lz,jz,Omegar,Omegaphi,Omegaz) HISTORY: 2013-09-08 - Written - Bovy (IAS) """ if len(args) == 5: #R,vR.vT, z, vz R,vR,vT, z, vz= args elif len(args) == 6: #R,vR.vT, z, vz, phi R,vR,vT, z, vz, phi= args else: self._parse_eval_args(*args) R= self._eval_R vR= self._eval_vR vT= self._eval_vT z= self._eval_z vz= self._eval_vz if isinstance(R,float): R= numpy.array([R]) vR= numpy.array([vR]) vT= numpy.array([vT]) z= numpy.array([z]) vz= numpy.array([vz]) if self._c: #pragma: no cover pass else: Lz= R*vT Lx= -z*vT Ly= z*vR-R*vz L2= Lx*Lx+Ly*Ly+Lz*Lz E= self._ip(R,z)+vR**2./2.+vT**2./2.+vz**2./2. L= numpy.sqrt(L2) #Actions Jphi= Lz Jz= L-numpy.fabs(Lz) Jr= self.amp/numpy.sqrt(-2.*E)\ -0.5*(L+numpy.sqrt((L2+4.*self.amp*self.b))) #Frequencies Omegar= (-2.*E)**1.5/self.amp Omegaz= 0.5*(1.+L/numpy.sqrt(L2+4.*self.amp*self.b))*Omegar Omegaphi= copy.copy(Omegaz) indx= Lz < 0. Omegaphi[indx]*= -1. return (Jr,Jphi,Jz,Omegar,Omegaphi,Omegaz) def _actionsFreqsAngles(self,*args,**kwargs): """ NAME: actionsFreqsAngles (_actionsFreqsAngles) PURPOSE: evaluate the actions, frequencies, and angles (jr,lz,jz,Omegar,Omegaphi,Omegaz,angler,anglephi,anglez) INPUT: Either: a) R,vR,vT,z,vz[,phi]: 1) floats: phase-space value for single object (phi is optional) (each can be a Quantity) 2) numpy.ndarray: [N] phase-space values for N objects (each can be a Quantity) b) Orbit instance: initial condition used if that's it, orbit(t) if there is a time given as well as the second argument OUTPUT: (jr,lz,jz,Omegar,Omegaphi,Omegaz,angler,anglephi,anglez) HISTORY: 2013-09-08 - Written - Bovy (IAS) """ if len(args) == 5: #R,vR.vT, z, vz pragma: no cover raise IOError("You need to provide phi when calculating angles") elif len(args) == 6: #R,vR.vT, z, vz, phi R,vR,vT, z, vz, phi= args else: self._parse_eval_args(*args) R= self._eval_R vR= self._eval_vR vT= self._eval_vT z= self._eval_z vz= self._eval_vz phi= self._eval_phi if isinstance(R,float): R= numpy.array([R]) vR= numpy.array([vR]) vT= numpy.array([vT]) z= numpy.array([z]) vz= numpy.array([vz]) phi= numpy.array([phi]) if self._c: #pragma: no cover pass else: Lz= R*vT Lx= -z*vT Ly= z*vR-R*vz L2= Lx*Lx+Ly*Ly+Lz*Lz E= self._ip(R,z)+vR**2./2.+vT**2./2.+vz**2./2. L= numpy.sqrt(L2) #Actions Jphi= Lz Jz= L-numpy.fabs(Lz) Jr= self.amp/numpy.sqrt(-2.*E)\ -0.5*(L+numpy.sqrt((L2+4.*self.amp*self.b))) #Frequencies Omegar= (-2.*E)**1.5/self.amp Omegaz= 0.5*(1.+L/numpy.sqrt(L2+4.*self.amp*self.b))*Omegar Omegaphi= copy.copy(Omegaz) indx= Lz < 0. Omegaphi[indx]*= -1. #Angles c= -self.amp/2./E-self.b e2= 1.-L2/self.amp/c*(1.+self.b/c) e= numpy.sqrt(e2) if self.b == 0.: coseta= 1/e*(1.-numpy.sqrt(R**2.+z**2.)/c) else: s= 1.+numpy.sqrt(1.+(R**2.+z**2.)/self.b**2.) coseta= 1/e*(1.-self.b/c*(s-2.)) pindx= (coseta > 1.)*(coseta < (1.+10.**-7.)) coseta[pindx]= 1. pindx= (coseta < -1.)*(coseta > (-1.-10.**-7.)) coseta[pindx]= -1. eta= numpy.arccos(coseta) costheta= z/numpy.sqrt(R**2.+z**2.) sintheta= R/numpy.sqrt(R**2.+z**2.) vrindx= (vR*sintheta+vz*costheta) < 0. eta[vrindx]= 2.*numpy.pi-eta[vrindx] angler= eta-e*c/(c+self.b)*numpy.sin(eta) tan11= numpy.arctan(numpy.sqrt((1.+e)/(1.-e))*numpy.tan(0.5*eta)) tan12= numpy.arctan(numpy.sqrt((1.+e+2.*self.b/c)/(1.-e+2.*self.b/c))*numpy.tan(0.5*eta)) vzindx= (-vz*sintheta+vR*costheta) > 0. tan11[tan11 < 0.]+= numpy.pi tan12[tan12 < 0.]+= numpy.pi pindx= (Lz/L > 1.)*(Lz/L < (1.+10.**-7.)) Lz[pindx]= L[pindx] pindx= (Lz/L < -1.)*(Lz/L > (-1.-10.**-7.)) Lz[pindx]= -L[pindx] i= numpy.arccos(Lz/L) sinpsi= costheta/numpy.sin(i) pindx= (sinpsi > 1.)*(sinpsi < (1.+10.**-7.)) sinpsi[pindx]= 1. pindx= (sinpsi < -1.)*(sinpsi > (-1.-10.**-7.)) sinpsi[pindx]= -1. psi= numpy.arcsin(sinpsi) psi[vzindx]= numpy.pi-psi[vzindx] psi= psi % (2.*numpy.pi) anglez= psi+Omegaz/Omegar*angler\ -tan11-1./numpy.sqrt(1.+4*self.amp*self.b/L2)*tan12 sinu= z/R/numpy.tan(i) pindx= (sinu > 1.)*(sinu < (1.+10.**-7.)) sinu[pindx]= 1. pindx= (sinu < -1.)*(sinu > (-1.-10.**-7.)) sinu[pindx]= -1. u= numpy.arcsin(sinu) u[vzindx]= numpy.pi-u[vzindx] Omega= phi-u anglephi= Omega anglephi[indx]-= anglez[indx] anglephi[True^indx]+= anglez[True^indx] angler= angler % (2.*numpy.pi) anglephi= anglephi % (2.*numpy.pi) anglez= anglez % (2.*numpy.pi) return (Jr,Jphi,Jz,Omegar,Omegaphi,Omegaz,angler,anglephi,anglez) def _EccZmaxRperiRap(self,*args,**kwargs): """ NAME: _EccZmaxRperiRap PURPOSE: evaluate the eccentricity, maximum height above the plane, peri- and apocenter for an isochrone potential INPUT: Either: a) R,vR,vT,z,vz[,phi]: 1) floats: phase-space value for single object (phi is optional) (each can be a Quantity) 2) numpy.ndarray: [N] phase-space values for N objects (each can be a Quantity) b) Orbit instance: initial condition used if that's it, orbit(t) if there is a time given as well as the second argument OUTPUT: (e,zmax,rperi,rap) HISTORY: 2017-12-22 - Written - Bovy (UofT) """ if len(args) == 5: #R,vR.vT, z, vz pragma: no cover R,vR,vT, z, vz= args elif len(args) == 6: #R,vR.vT, z, vz, phi R,vR,vT, z, vz, phi= args else: self._parse_eval_args(*args) R= self._eval_R vR= self._eval_vR vT= self._eval_vT z= self._eval_z vz= self._eval_vz if isinstance(R,float): R= numpy.array([R]) vR= numpy.array([vR]) vT= numpy.array([vT]) z= numpy.array([z]) vz= numpy.array([vz]) if self._c: #pragma: no cover pass else: Lz= R*vT Lx= -z*vT Ly= z*vR-R*vz L2= Lx*Lx+Ly*Ly+Lz*Lz E= self._ip(R,z)+vR**2./2.+vT**2./2.+vz**2./2. if self.b == 0: warnings.warn("zmax for point-mass (b=0) isochrone potential is only approximate, because it assumes that zmax is attained at rap, which is not necessarily the case",galpyWarning) a= -self.amp/2./E me2= L2/self.amp/a e= numpy.sqrt(1.-me2) rperi= a*(1.-e) rap= a*(1.+e) else: smin= 0.5*((2.*E-self.amp/self.b)\ +numpy.sqrt((2.*E-self.amp/self.b)**2. +2.*E*(4.*self.amp/self.b+L2/self.b**2.)))/E smax= 2.-self.amp/E/self.b-smin rperi= smin*numpy.sqrt(1.-2./smin)*self.b rap= smax*numpy.sqrt(1.-2./smax)*self.b return ((rap-rperi)/(rap+rperi),rap*numpy.sqrt(1.-Lz**2./L2), rperi,rap)