# 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
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:

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:

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:

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)