"""
Analytic transformations to action-angle coordinates.
"""
# Third-party
import numpy as np
from astropy.constants import G
import astropy.coordinates as coord
import astropy.units as u
# Project
from ..potential import Hamiltonian, PotentialBase, IsochronePotential
__all__ = ['isochrone_to_aa', 'harmonic_oscillator_to_aa']
[docs]def isochrone_to_aa(w, potential):
"""
Transform the input cartesian position and velocity to action-angle
coordinates in the Isochrone potential. See Section 3.5.2 in
Binney & Tremaine (2008), and be aware of the errata entry for
Eq. 3.225.
This transformation is analytic and can be used as a "toy potential"
in the Sanders & Binney (2014) formalism for computing action-angle
coordinates in any potential.
.. note::
This function is included as a method of the
:class:`~gala.potential.IsochronePotential` and it is recommended
to call :meth:`~gala.potential.IsochronePotential.phase_space()`
instead.
Parameters
----------
w : :class:`gala.dynamics.PhaseSpacePosition`, :class:`gala.dynamics.Orbit`
potential : :class:`gala.potential.IsochronePotential`, dict
An instance of the potential to use for computing the transformation
to angle-action coordinates. Or, a dictionary of parameters used to
define an :class:`gala.potential.IsochronePotential` instance.
Returns
-------
actions : :class:`numpy.ndarray`
An array of actions computed from the input positions and velocities.
angles : :class:`numpy.ndarray`
An array of angles computed from the input positions and velocities.
freqs : :class:`numpy.ndarray`
An array of frequencies computed from the input positions and velocities.
"""
if not isinstance(potential, PotentialBase):
potential = IsochronePotential(**potential)
usys = potential.units
GM = (G*potential.parameters['m']).decompose(usys).value
b = potential.parameters['b'].decompose(usys).value
E = w.energy(Hamiltonian(potential)).decompose(usys).value
E = np.squeeze(E)
if np.any(E > 0.):
raise ValueError("Unbound particle. (E = {})".format(E))
# convert position, velocity to spherical polar coordinates
w_sph = w.represent_as(coord.PhysicsSphericalRepresentation)
r, phi, theta = map(np.squeeze, [w_sph.r.decompose(usys).value,
w_sph.phi.radian,
w_sph.theta.radian])
ang_unit = u.radian/usys['time']
vr, phi_dot, theta_dot = map(np.squeeze,
[w_sph.radial_velocity.decompose(usys).value,
w_sph.pm_phi.to(ang_unit).value,
w_sph.pm_theta.to(ang_unit).value])
vtheta = r*theta_dot
# ----------------------------
# Compute the actions
# ----------------------------
L_vec = np.squeeze(w.angular_momentum().decompose(usys).value)
Lz = L_vec[2]
L = np.linalg.norm(L_vec, axis=0)
# Radial action
Jr = GM / np.sqrt(-2*E) - 0.5*(L + np.sqrt(L*L + 4*GM*b))
# compute the three action variables
actions = np.array([Jr, Lz, L - np.abs(Lz)]) # Jr, Jphi, Jtheta
# ----------------------------
# Angles
# ----------------------------
c = GM / (-2*E) - b
e = np.sqrt(1 - L*L*(1 + b/c) / GM / c)
# Compute theta_r using eta
tmp1 = r*vr / np.sqrt(-2.*E)
tmp2 = b + c - np.sqrt(b*b + r*r)
eta = np.arctan2(tmp1, tmp2)
thetar = eta - e*c*np.sin(eta) / (c + b) # same as theta3
# Compute theta_z
psi = np.arctan2(np.cos(theta), -np.sin(theta)*r*vtheta/L)
psi[np.abs(vtheta) <= 1e-10] = np.pi/2. # blows up for small vtheta
omega_th = 0.5 * (1 + L/np.sqrt(L*L + 4*GM*b))
a = np.sqrt((1+e) / (1-e))
ap = np.sqrt((1 + e + 2*b/c) / (1 - e + 2*b/c))
def F(x, y):
z = np.zeros_like(x)
ix = y > np.pi/2.
z[ix] = np.pi/2. - np.arctan(np.tan(np.pi/2.-0.5*y[ix])/x[ix])
ix = y < -np.pi/2.
z[ix] = -np.pi/2. + np.arctan(np.tan(np.pi/2.+0.5*y[ix])/x[ix])
ix = (y <= np.pi/2) & (y >= -np.pi/2)
z[ix] = np.arctan(x[ix]*np.tan(0.5*y[ix]))
return z
A = omega_th*thetar - F(a, eta) - F(ap, eta)/np.sqrt(1 + 4*GM*b/L/L)
thetaz = psi + A
LR = Lz/L
sinu = (LR/np.sqrt(1.-LR*LR)/np.tan(theta))
uu = np.arcsin(sinu)
uu[sinu > 1.] = np.pi/2.
uu[sinu < -1.] = -np.pi/2.
uu[vtheta > 0.] = np.pi - uu[vtheta > 0.]
thetap = phi - uu + np.sign(Lz)*thetaz
angles = np.array([thetar, thetap, thetaz])
angles = angles % (2*np.pi)
# ----------------------------
# Frequencies
# ----------------------------
freqs = np.zeros_like(actions)
omega_r = GM**2 / (Jr + 0.5*(L + np.sqrt(L*L + 4*GM*b)))**3
freqs[0] = omega_r
freqs[1] = np.sign(actions[1]) * omega_th * omega_r
freqs[2] = omega_th * omega_r
a_unit = (1 * usys['angular momentum'] / usys['mass']).decompose(usys).unit
f_unit = (1 * usys['angular speed']).decompose(usys).unit
return actions*a_unit, angles*u.radian, freqs*f_unit
def isochrone_to_xv(actions, angles, potential):
"""
Transform the input actions and angles to ordinary phase space (position
and velocity) in cartesian coordinates. See Section 3.5.2 in
Binney & Tremaine (2008), and be aware of the errata entry for
Eq. 3.225.
.. note::
This function is included as a method of the
:class:`~gala.potential.IsochronePotential`
and it is recommended to call
:meth:`~gala.potential.IsochronePotential.action_angle()` instead.
Parameters
----------
actions : array_like
Action variables. Must have shape ``(3, N)`` or ``(3,)``.
angles : array_like
Angle variables. Must have shape ``(3, N)`` or ``(3,)``.
Should be in radians.
potential : :class:`gala.potential.IsochronePotential`
An instance of the potential to use for computing the transformation
to angle-action coordinates.
Returns
-------
x : :class:`numpy.ndarray`
An array of cartesian positions computed from the input
angles and actions.
v : :class:`numpy.ndarray`
An array of cartesian velocities computed from the input
angles and actions.
"""
raise NotImplementedError("Implementation not supported until working with "
"angle-action variables has a better API.")
# actions = atleast_2d(actions, insert_axis=1).copy()
# angles = atleast_2d(angles, insert_axis=1).copy()
# usys = potential.units
# GM = (G*potential.parameters['m']).decompose(usys).value
# b = potential.parameters['b'].decompose(usys).value
# # actions
# Jr = actions[0]
# Lz = actions[1]
# L = actions[2] + np.abs(Lz)
# # angles
# theta_r, theta_phi, theta_theta = angles
# # get longitude of ascending node
# theta_1 = theta_phi - np.sign(Lz)*theta_theta
# Omega = theta_1
# # Ly = -np.cos(Omega) * np.sqrt(L**2 - Lz**2)
# # Lx = np.sqrt(L**2 - Ly**2 - Lz**2)
# cosi = Lz/L
# sini = np.sqrt(1 - cosi**2)
# # Hamiltonian (energy)
# H = -2. * GM**2 / (2.*Jr + L + np.sqrt(4.*b*GM + L**2))**2
# if np.any(H > 0.):
# raise ValueError("Unbound particle. (E = {})".format(H))
# # Eq. 3.240
# c = -GM / (2.*H) - b
# e = np.sqrt(1 - L*L*(1 + b/c) / GM / c)
# # solve for eta
# theta_3 = theta_r
# eta_func = lambda x: x - e*c/(b+c)*np.sin(x) - theta_3
# eta_func_prime = lambda x: 1 - e*c/(b+c)*np.cos(x)
# # use newton's method to find roots
# niter = 100
# eta = np.ones_like(theta_3)*np.pi/2.
# for i in range(niter):
# eta -= eta_func(eta)/eta_func_prime(eta)
# # TODO: when to do this???
# eta -= 2*np.pi
# r = c*np.sqrt((1-e*np.cos(eta)) * (1-e*np.cos(eta) + 2*b/c))
# vr = np.sqrt(GM/(b+c))*(c*e*np.sin(eta))/r
# theta_2 = theta_theta
# Omega_23 = 0.5*(1 + L / np.sqrt(L**2 + 4*GM*b))
# a = np.sqrt((1+e) / (1-e))
# ap = np.sqrt((1 + e + 2*b/c) / (1 - e + 2*b/c))
# def F(x, y):
# z = np.zeros_like(x)
# ix = y>np.pi/2.
# z[ix] = np.pi/2. - np.arctan(np.tan(np.pi/2.-0.5*y[ix])/x[ix])
# ix = y<-np.pi/2.
# z[ix] = -np.pi/2. + np.arctan(np.tan(np.pi/2.+0.5*y[ix])/x[ix])
# ix = (y<=np.pi/2) & (y>=-np.pi/2)
# z[ix] = np.arctan(x[ix]*np.tan(0.5*y[ix]))
# return z
# theta_2[Lz < 0] -= 2*np.pi
# theta_3 -= 2*np.pi
# A = Omega_23*theta_3 - F(a, eta) - F(ap, eta)/np.sqrt(1 + 4*GM*b/L/L)
# psi = theta_2 - A
# # theta
# theta = np.arccos(np.sin(psi)*sini)
# vtheta = L*sini*np.cos(psi)/np.cos(theta)
# vtheta = -L*sini*np.cos(psi)/np.sin(theta)/r
# vphi = Lz / (r*np.sin(theta))
# d_phi = vphi / (r*np.sin(theta))
# d_theta = vtheta / r
# # phi
# sinu = np.sin(psi)*cosi/np.sin(theta)
# uu = np.arcsin(sinu)
# uu[sinu > 1.] = np.pi/2.
# uu[sinu < -1.] = -np.pi/2.
# uu[vtheta > 0.] = np.pi - uu[vtheta > 0.]
# sinu = cosi/sini * np.cos(theta)/np.sin(theta)
# phi = (uu + Omega) % (2*np.pi)
# # We now need to convert from spherical polar coord to cart. coord.
# pos = coord.PhysicsSphericalRepresentation(r=r*u.dimensionless_unscaled,
# phi=phi*u.rad, theta=theta*u.rad)
# pos = pos.represent_as(coord.CartesianRepresentation)
# x = pos.xyz.value
# vel = coord.PhysicsSphericalDifferential(d_phi=d_phi, d_theta=d_theta, d_r=vr)
# v = vel.represent_as(coord.CartesianDifferential, base=pos).d_xyz.value
# return x, v
[docs]def harmonic_oscillator_to_aa(w, potential):
"""
Transform the input cartesian position and velocity to action-angle
coordinates for the Harmonic Oscillator potential.
This transformation is analytic and can be used as a "toy potential"
in the Sanders & Binney (2014) formalism for computing action-angle
coordinates in any potential.
.. note::
This function is included as a method of the
:class:`~gala.potential.HarmonicOscillatorPotential`
and it is recommended to call
:meth:`~gala.potential.HarmonicOscillatorPotential.action_angle()`
instead.
Parameters
----------
w : :class:`gala.dynamics.PhaseSpacePosition`, :class:`gala.dynamics.Orbit`
potential : Potential
"""
usys = potential.units
if usys is not None:
x = w.xyz.decompose(usys).value
v = w.v_xyz.decompose(usys).value
else:
x = w.xyz.value
v = w.v_xyz.value
_new_omega_shape = (3,) + tuple([1]*(len(x.shape)-1))
# compute actions -- just energy (hamiltonian) over frequency
if usys is None:
usys = []
try:
omega = potential.parameters['omega'].reshape(_new_omega_shape).decompose(usys).value
except AttributeError: # not a Quantity
omega = potential.parameters['omega'].reshape(_new_omega_shape)
action = (v**2 + (omega*x)**2)/(2.*omega)
angle = np.arctan(-v / omega / x)
angle[x == 0] = -np.sign(v[x == 0])*np.pi/2.
angle[x < 0] += np.pi
freq = potential.parameters['omega'].decompose(usys).value
if usys is not None and usys:
a_unit = (1*usys['angular momentum']/usys['mass']).decompose(usys).unit
f_unit = (1*usys['angular speed']).decompose(usys).unit
return action*a_unit, (angle % (2.*np.pi))*u.radian, freq*f_unit
else:
return action*u.one, (angle % (2.*np.pi))*u.one, freq*u.one
def harmonic_oscillator_to_xv(actions, angles, potential):
"""
Transform the input action-angle coordinates to cartesian
position and velocity for the Harmonic Oscillator potential.
.. note::
This function is included as a method of the
:class:`~gala.potential.HarmonicOscillatorPotential`
and it is recommended to call
:meth:`~gala.potential.HarmonicOscillatorPotential.phase_space()` instead.
Parameters
----------
actions : array_like
angles : array_like
potential : Potential
"""
raise NotImplementedError("Implementation not supported until working with "
"angle-action variables has a better API.")
# TODO: bug in below...
# omega = potential.parameters['omega'].decompose(potential.units).value
# x = np.sqrt(2*actions/omega[None]) * np.sin(angles)
# v = np.sqrt(2*actions*omega[None]) * np.cos(angles)
# return x, v