Source code for galpy.actionAngle.actionAngleIsochroneInverse

###############################################################################
#   actionAngle: a Python module to calculate  actions, angles, and frequencies
#
#      class: actionAngleIsochroneInverse
#
#             Calculate (x,v) coordinates for the Isochrone potential from 
#             given actions-angle coordinates
#
###############################################################################
import numpy
from scipy import optimize
from ..util import conversion
from ..potential import IsochronePotential
from .actionAngleInverse import actionAngleInverse
[docs]class actionAngleIsochroneInverse(actionAngleInverse): """Inverse action-angle formalism for the isochrone potential, on the Jphi, Jtheta system of Binney & Tremaine (2008); following McGill & Binney (1990) for transformations"""
[docs] def __init__(self,*args,**kwargs): """ NAME: __init__ PURPOSE: initialize an actionAngleIsochroneInverse 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: 2017-11-14 - Started - Bovy (UofT) """ actionAngleInverse.__init__(self,*args,**kwargs) 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= conversion.parse_length(kwargs['b'],ro=self._ro) rb= numpy.sqrt(self.b**2.+1.) self.amp= (self.b+rb)**2.*rb # In case we ever decide to implement this in C... 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,jr,jphi,jz,angler,anglephi,anglez,**kwargs): """ NAME: __call__ PURPOSE: evaluate the phase-space coordinates (x,v) for a number of angles on a single torus INPUT: jr - radial action (scalar) jphi - azimuthal action (scalar) jz - vertical action (scalar) angler - radial angle (array [N]) anglephi - azimuthal angle (array [N]) anglez - vertical angle (array [N]) tol= (object-wide value) goal for |dJ|/|J| along the torus OUTPUT: [R,vR,vT,z,vz,phi] HISTORY: 2017-11-14 - Written - Bovy (UofT) """ return self._xvFreqs(jr,jphi,jz,angler,anglephi,anglez,**kwargs)[:6] def _xvFreqs(self,jr,jphi,jz,angler,anglephi,anglez,**kwargs): """ NAME: xvFreqs PURPOSE: evaluate the phase-space coordinates (x,v) for a number of angles on a single torus as well as the frequencies INPUT: jr - radial action (scalar) jphi - azimuthal action (scalar) jz - vertical action (scalar) angler - radial angle (array [N]) anglephi - azimuthal angle (array [N]) anglez - vertical angle (array [N]) OUTPUT: ([R,vR,vT,z,vz,phi],OmegaR,Omegaphi,Omegaz) HISTORY: 2017-11-15 - Written - Bovy (UofT) """ L= jz+numpy.fabs(jphi) # total angular momentum L2= L**2. sqrtfourbkL2= numpy.sqrt(L2+4.*self.b*self.amp) H= -2.*self.amp**2./(2.*jr+L+sqrtfourbkL2)**2. # Calculate the frequencies omegar= (-2.*H)**1.5/self.amp omegaz= (1.+L/sqrtfourbkL2)/2.*omegar # Start on getting the coordinates a= -self.amp/2./H-self.b ab= a+self.b e= numpy.sqrt(1.+L2/(2.*H*a**2.)) # Solve Kepler's-ish equation; ar must be between 0 and 2pi angler= (numpy.atleast_1d(angler) % (-2.*numpy.pi)) % (2.*numpy.pi) anglephi= numpy.atleast_1d(anglephi) anglez= numpy.atleast_1d(anglez) eta= numpy.empty(len(angler)) for ii,ar in enumerate(angler): try: eta[ii]= optimize.newton(lambda x: x-a*e/ab*numpy.sin(x)-ar, 0., lambda x: 1-a*e/ab*numpy.cos(x)) except RuntimeError: # Newton-Raphson did not converge, this has to work, # bc 0 <= ra < 2pi the following start x have different signs eta[ii]= optimize.brentq(lambda x: x-a*e/ab*numpy.sin(x)-ar, 0.,2.*numpy.pi) coseta= numpy.cos(eta) r= a*numpy.sqrt((1.-e*coseta)*(1.-e*coseta+2.*self.b/a)) vr= numpy.sqrt(self.amp/ab)*a*e*numpy.sin(eta)/r taneta2= numpy.tan(eta/2.) tan11= numpy.arctan(numpy.sqrt((1.+e)/(1.-e))*taneta2) tan12= numpy.arctan(\ numpy.sqrt((a*(1.+e)+2.*self.b)/(a*(1.-e)+2.*self.b))*taneta2) tan11[tan11 < 0.]+= numpy.pi tan12[tan12 < 0.]+= numpy.pi Lambdaeta= tan11+L/sqrtfourbkL2*tan12 psi= anglez-omegaz/omegar*angler+Lambdaeta lowerl= numpy.sqrt(1.-jphi**2./L2) sintheta= numpy.sin(psi)*lowerl costheta= numpy.sqrt(1.-sintheta**2.) vtheta= L*lowerl*numpy.cos(psi)/costheta/r R= r*costheta z= r*sintheta vR= vr*costheta-vtheta*sintheta vz= vr*sintheta+vtheta*costheta sinu= sintheta/costheta*jphi/L/lowerl u= numpy.arcsin(sinu) u[vtheta < 0.]= numpy.pi-u[vtheta < 0.] phi= anglephi-numpy.sign(jphi)*anglez+u # For non-inclined orbits, phi == psi phi[True^numpy.isfinite(phi)]= psi[True^numpy.isfinite(phi)] phi= phi % (2.*numpy.pi) phi[phi < 0.]+= 2.*numpy.pi return (R,vR,jphi/R,z,vz,phi, omegar,numpy.sign(jphi)*omegaz,omegaz) def _Freqs(self,jr,jphi,jz,**kwargs): """ NAME: Freqs PURPOSE: return the frequencies corresponding to a torus INPUT: jr - radial action (scalar) jphi - azimuthal action (scalar) jz - vertical action (scalar) OUTPUT: (OmegaR,Omegaphi,Omegaz) HISTORY: 2017-11-15 - Written - Bovy (UofT) """ L= jz+numpy.fabs(jphi) # total angular momentum sqrtfourbkL2= numpy.sqrt(L**2.+4.*self.b*self.amp) H= -2.*self.amp**2./(2.*jr+L+sqrtfourbkL2)**2. # Calculate the frequencies omegar= (-2.*H)**1.5/self.amp omegaz= (1.+L/sqrtfourbkL2)/2.*omegar return (omegar,numpy.sign(jphi)*omegaz,omegaz)