# Orbit¶

class gala.dynamics.orbit.Orbit(pos, vel, t=None, hamiltonian=None, potential=None, frame=None)[source]

Represents an orbit: positions and velocities (conjugate momenta) as a function of time.

The class can be instantiated with Astropy representation objects (e.g., CartesianRepresentation), Astropy Quantity objects, or plain Numpy arrays.

If passing in Quantity or Numpy array instances for both position and velocity, they are assumed to be Cartesian. Array inputs are interpreted as dimensionless quantities. The input position and velocity objects can have an arbitrary number of (broadcastable) dimensions. For Quantity or array inputs, the first axes have special meaning:

- axis=0 is the coordinate dimension (e.g., x, y, z)
- axis=1 is the time dimension


So if the input position array, pos, has shape pos.shape = (3, 100), this would be a 3D orbit at 100 times (pos[0] is x, pos[1] is y, etc.). For representing multiple orbits, the position array could have 3 axes, e.g., it might have shape pos.shape = (3, 100, 8), where this is interpreted as a 3D position at 100 times for 8 different orbits. The same is true for velocity. The position and velocity arrays must have the same shape.

If a time argument is specified, the position and velocity arrays must have the same number of timesteps as the length of the time object:

len(t) == pos.shape[1]

Parameters: pos : BaseRepresentation, Quantity, array_like Positions. If a numpy array (e.g., has no units), this will be stored as a dimensionless Quantity. See the note above about the assumed meaning of the axes of this object. vel : BaseDifferential, Quantity, array_like Velocities. If a numpy array (e.g., has no units), this will be stored as a dimensionless Quantity. See the note above about the assumed meaning of the axes of this object. t : array_like, Quantity (optional) Array of times. If a numpy array (e.g., has no units), this will be stored as a dimensionless Quantity. hamiltonian : Hamiltonian (optional) The Hamiltonian that the orbit was integrated in.

Attributes Summary

 data hamiltonian norbits ntimes pos_components representation_mappings shape This is not the shape of the position or velocity arrays. vel_components

Methods Summary

 align_circulation_with_z([circulation]) If the input orbit is a tube orbit, this function aligns the circulation axis with the z axis and returns a copy. angular_momentum() Compute the angular momentum for the phase-space positions contained in this object. apocenter([return_times, func, …]) Estimate the apocenter(s) of the orbit by identifying local maxima in the spherical radius and interpolating between timesteps near the maxima. circulation() Determine which axes the Orbit circulates around by checking whether there is a change of sign of the angular momentum about an axis. eccentricity(**kw) Returns the eccentricity computed from the mean apocenter and mean pericenter. energy([hamiltonian]) The total energy per unit mass: estimate_period([radial]) Estimate the period of the orbit. from_hdf5(f) Load an object from an HDF5 file. from_w(w[, units]) get_components(which) Get the component name dictionary for the desired object. kinetic_energy() The kinetic energy per unit mass: orbit_gen() Generator for iterating over each orbit. pericenter([return_times, func, …]) Estimate the pericenter(s) of the orbit by identifying local minima in the spherical radius and interpolating between timesteps near the minima. plot([components, units, auto_aspect]) Plot the positions in all projections. potential_energy([potential]) The potential energy per unit mass: represent_as(new_pos[, new_vel]) Represent the position and velocity of the orbit in an alternate coordinate system. to_coord_frame(frame[, galactocentric_frame]) Transform the orbit from Galactocentric, cartesian coordinates to Heliocentric coordinates in the specified Astropy coordinate frame. to_frame(frame[, current_frame]) TODO: to_hdf5(f) Serialize this object to an HDF5 file. w([units]) This returns a single array containing the phase-space positions. zmax([return_times, func, interp_kwargs, …]) Estimate the maximum z height of the orbit by identifying local maxima in the absolute value of the z position and interpolating between timesteps near the maxima.

Attributes Documentation

data
hamiltonian
norbits
ntimes
pos_components
representation_mappings = {<class 'astropy.coordinates.representation.CartesianRepresentation'>: [RepresentationMapping(repr_name='xyz', new_name='xyz', default_unit='recommended')], <class 'astropy.coordinates.representation.SphericalCosLatDifferential'>: [RepresentationMapping(repr_name='d_lon_coslat', new_name='pm_lon_coslat', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_lat', new_name='pm_lat', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_distance', new_name='radial_velocity', default_unit='recommended')], <class 'astropy.coordinates.representation.SphericalDifferential'>: [RepresentationMapping(repr_name='d_lon', new_name='pm_lon', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_lat', new_name='pm_lat', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_distance', new_name='radial_velocity', default_unit='recommended')], <class 'astropy.coordinates.representation.PhysicsSphericalDifferential'>: [RepresentationMapping(repr_name='d_phi', new_name='pm_phi', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_theta', new_name='pm_theta', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_r', new_name='radial_velocity', default_unit='recommended')], <class 'astropy.coordinates.representation.CartesianDifferential'>: [RepresentationMapping(repr_name='d_x', new_name='v_x', default_unit='recommended'), RepresentationMapping(repr_name='d_y', new_name='v_y', default_unit='recommended'), RepresentationMapping(repr_name='d_z', new_name='v_z', default_unit='recommended'), RepresentationMapping(repr_name='d_xyz', new_name='v_xyz', default_unit='recommended')], <class 'astropy.coordinates.representation.CylindricalDifferential'>: [RepresentationMapping(repr_name='d_rho', new_name='v_rho', default_unit='recommended'), RepresentationMapping(repr_name='d_phi', new_name='pm_phi', default_unit='recommended'), RepresentationMapping(repr_name='d_z', new_name='v_z', default_unit='recommended')], <class 'gala.dynamics.representation_nd.NDCartesianRepresentation'>: [RepresentationMapping(repr_name='xyz', new_name='xyz', default_unit='recommended')], <class 'gala.dynamics.representation_nd.NDCartesianDifferential'>: [RepresentationMapping(repr_name='d_xyz', new_name='v_xyz', default_unit='recommended'), RegexRepresentationMapping(repr_name='d_x([0-9])', new_name='v_x{0}', default_unit='recommended')], <class 'astropy.coordinates.representation.UnitSphericalCosLatDifferential'>: [RepresentationMapping(repr_name='d_lon_coslat', new_name='pm_lon_coslat', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_lat', new_name='pm_lat', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_distance', new_name='radial_velocity', default_unit='recommended')], <class 'astropy.coordinates.representation.UnitSphericalDifferential'>: [RepresentationMapping(repr_name='d_lon', new_name='pm_lon', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_lat', new_name='pm_lat', default_unit=Unit("mas / yr")), RepresentationMapping(repr_name='d_distance', new_name='radial_velocity', default_unit='recommended')]}
shape

This is not the shape of the position or velocity arrays. That is accessed by doing, e.g., obj.x.shape.

vel_components

Methods Documentation

align_circulation_with_z(circulation=None)[source]

If the input orbit is a tube orbit, this function aligns the circulation axis with the z axis and returns a copy.

Parameters: circulation : array_like (optional) Array of bits that specify the axis about which the orbit circulates. If not provided, will compute this using circulation(). See that method for more information. orb : Orbit A copy of the original orbit object with circulation aligned with the z axis.
angular_momentum()

Compute the angular momentum for the phase-space positions contained in this object:

.. math::

boldsymbol{{L}} = boldsymbol{{q}} times boldsymbol{{p}}

Returns: L : Quantity Array of angular momentum vectors.

Examples

>>> import numpy as np
>>> import astropy.units as u
>>> pos = np.array([1., 0, 0]) * u.au
>>> vel = np.array([0, 2*np.pi, 0]) * u.au/u.yr
>>> w = PhaseSpacePosition(pos, vel)
>>> w.angular_momentum() # doctest: +FLOAT_CMP
<Quantity [0.        ,0.        ,6.28318531] AU2 / yr>

apocenter(return_times=False, func=<function mean>, interp_kwargs=None, minimize_kwargs=None, approximate=False)[source]

Estimate the apocenter(s) of the orbit by identifying local maxima in the spherical radius and interpolating between timesteps near the maxima.

By default, this returns the mean of all local maxima (apocenters). To get, e.g., the largest apocenter, pass in func=np.max. To get all apocenters, pass in func=None.

Parameters: func : func (optional) A function to evaluate on all of the identified apocenter times. return_times : bool (optional) Also return the apocenter times. interp_kwargs : dict (optional) Keyword arguments to be passed to scipy.interpolate.InterpolatedUnivariateSpline. minimize_kwargs : dict (optional) Keyword arguments to be passed to scipy.optimize.minimize. approximate : bool (optional) Compute an approximate apocenter by skipping interpolation. apo : float, ndarray Either a single number or an array of apocenters. times : ndarray (optional, see return_times) If return_times=True, also returns an array of the apocenter times.
circulation()[source]

Determine which axes the Orbit circulates around by checking whether there is a change of sign of the angular momentum about an axis. Returns a 2D array with ndim integers per orbit point. If a box orbit, all integers will be 0. A 1 indicates circulation about the corresponding axis.

TODO: clockwise / counterclockwise?

For example, for a single 3D orbit:

• Box and boxlet = [0,0,0]
• z-axis (short-axis) tube = [0,0,1]
• x-axis (long-axis) tube = [1,0,0]
Returns: circulation : numpy.ndarray An array that specifies whether there is circulation about any of the axes of the input orbit. For a single orbit, will return a 1D array, but for multiple orbits, the shape will be (3, norbits).
eccentricity(**kw)[source]

Returns the eccentricity computed from the mean apocenter and mean pericenter.

$e = \frac{r_{\rm apo} - r_{\rm per}}{r_{\rm apo} + r_{\rm per}}$
Parameters: **kw Any keyword arguments passed to apocenter() and pericenter(). For example, approximate=True. ecc : float The orbital eccentricity.
energy(hamiltonian=None)[source]

The total energy per unit mass:

Returns: E : Quantity The total energy.
estimate_period(radial=True)[source]

Estimate the period of the orbit. By default, computes the radial period. If radial==False, this returns period estimates for each dimension of the orbit.

Parameters: radial : bool (optional) What period to estimate. If True, estimates the radial period. If False, estimates period in each dimension, e.g., if the orbit is 3D, along x, y, and z. T : Quantity The period or periods.
classmethod from_hdf5(f)[source]

Load an object from an HDF5 file.

Requires h5py.

Parameters: f : str, h5py.File Either the filename or an open HDF5 file.
classmethod from_w(w, units=None, **kwargs)
get_components(which)

Get the component name dictionary for the desired object.

The returned dictionary maps component names on this class to component names on the desired object.

Parameters: which : str Can either be 'pos' or 'vel' to get the components for the position or velocity object.
kinetic_energy()

The kinetic energy per unit mass:

$E_K = \frac{1}{2} \, |\boldsymbol{v}|^2$
Returns: E : Quantity The kinetic energy.
orbit_gen()[source]

Generator for iterating over each orbit.

pericenter(return_times=False, func=<function mean>, interp_kwargs=None, minimize_kwargs=None, approximate=False)[source]

Estimate the pericenter(s) of the orbit by identifying local minima in the spherical radius and interpolating between timesteps near the minima.

By default, this returns the mean of all local minima (pericenters). To get, e.g., the minimum pericenter, pass in func=np.min. To get all pericenters, pass in func=None.

Parameters: func : func (optional) A function to evaluate on all of the identified pericenter times. return_times : bool (optional) Also return the pericenter times. interp_kwargs : dict (optional) Keyword arguments to be passed to scipy.interpolate.InterpolatedUnivariateSpline. minimize_kwargs : dict (optional) Keyword arguments to be passed to scipy.optimize.minimize. approximate : bool (optional) Compute an approximate pericenter by skipping interpolation. peri : float, ndarray Either a single number or an array of pericenters. times : ndarray (optional, see return_times) If return_times=True, also returns an array of the pericenter times.
plot(components=None, units=None, auto_aspect=True, **kwargs)[source]

Plot the positions in all projections. This is a wrapper around plot_projections for fast access and quick visualization. All extra keyword arguments are passed to that function (the docstring for this function is included here for convenience).

Parameters: components : iterable (optional) A list of component names (strings) to plot. By default, this is the Cartesian positions ['x', 'y', 'z']. To plot Cartesian velocities, pass in the velocity component names ['v_x', 'v_y', 'v_z']. If the representation is different, the component names will be different. For example, for a Cylindrical representation, the components are ['rho', 'phi', 'z'] and ['v_rho', 'pm_phi', 'v_z']. units : UnitBase, iterable (optional) A single unit or list of units to display the components in. auto_aspect : bool (optional) Automatically enforce an equal aspect ratio. relative_to : bool (optional) Plot the values relative to this value or values. autolim : bool (optional) Automatically set the plot limits to be something sensible. axes : array_like (optional) Array of matplotlib Axes objects. subplots_kwargs : dict (optional) Dictionary of kwargs passed to subplots(). labels : iterable (optional) List or iterable of axis labels as strings. They should correspond to the dimensions of the input orbit. plot_function : callable (optional) The matplotlib plot function to use. By default, this is scatter(), but can also be, e.g., plot(). **kwargs All other keyword arguments are passed to the plot_function. You can pass in any of the usual style kwargs like color=..., marker=..., etc. fig : Figure
potential_energy(potential=None)[source]

The potential energy per unit mass:

$E_\Phi = \Phi(\boldsymbol{q})$
Returns: E : Quantity The potential energy.
represent_as(new_pos, new_vel=None)[source]

Represent the position and velocity of the orbit in an alternate coordinate system. Supports any of the Astropy coordinates representation classes.

Parameters: new_pos : BaseRepresentation The type of representation to generate. Must be a class (not an instance), or the string name of the representation class. new_vel : BaseDifferential (optional) Class in which any velocities should be represented. Must be a class (not an instance), or the string name of the differential class. If None, uses the default differential for the new position class. new_orbit : gala.dynamics.Orbit
to_coord_frame(frame, galactocentric_frame=None, **kwargs)

Transform the orbit from Galactocentric, cartesian coordinates to Heliocentric coordinates in the specified Astropy coordinate frame.

Parameters: frame : BaseCoordinateFrame The class or frame instance specifying the desired output frame. For example, ICRS. galactocentric_frame : Galactocentric This is the assumed frame that the position and velocity of this object are in. The Galactocentric instand should have parameters specifying the position and motion of the sun in the Galactocentric frame, but no data. c : BaseCoordinateFrame An instantiated coordinate frame containing the positions and velocities from this object transformed to the specified coordinate frame.
to_frame(frame, current_frame=None, **kwargs)[source]

TODO:

Parameters: frame : gala.potential.CFrameBase The frame to transform to. current_frame : gala.potential.CFrameBase (optional) If the Orbit has no associated Hamiltonian, this specifies the current frame of the orbit. orbit : gala.dynamics.Orbit The orbit in the new reference frame.
to_hdf5(f)[source]

Serialize this object to an HDF5 file.

Requires h5py.

Parameters: f : str, h5py.File Either the filename or an open HDF5 file.
w(units=None)[source]

This returns a single array containing the phase-space positions.

Parameters: units : UnitSystem (optional) The unit system to represent the position and velocity in before combining into the full array. w : ndarray A numpy array of all positions and velocities, without units. Will have shape (2*ndim,...).
zmax(return_times=False, func=<function mean>, interp_kwargs=None, minimize_kwargs=None, approximate=False)[source]

Estimate the maximum z height of the orbit by identifying local maxima in the absolute value of the z position and interpolating between timesteps near the maxima.

By default, this returns the mean of all local maxima. To get, e.g., the largest z excursion, pass in func=np.max. To get all z maxima, pass in func=None.

Parameters: func : func (optional) A function to evaluate on all of the identified z maximum times. return_times : bool (optional) Also return the times of maximum. interp_kwargs : dict (optional) Keyword arguments to be passed to scipy.interpolate.InterpolatedUnivariateSpline. minimize_kwargs : dict (optional) Keyword arguments to be passed to scipy.optimize.minimize. approximate : bool (optional) Compute approximate values by skipping interpolation. zs : float, ndarray Either a single number or an array of maximum z heights. times : ndarray (optional, see return_times) If return_times=True, also returns an array of the apocenter times.