Orbit#

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

Bases: PhaseSpacePosition

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:
posrepresentation, quantity_like, or 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.

veldifferential, quantity_like, or 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.

tarray_like, Quantity (optional)

Array of times. If a numpy array (e.g., has no units), this will be stored as a dimensionless Quantity.

hamiltonianHamiltonian (optional)

The Hamiltonian that the orbit was integrated in.

copybool, optional

If True, the input arrays are copied. If False, the input data is referenced directly (if possible). Default is True.

Represents phase-space positions, i.e. positions and conjugate momenta (velocities).

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

If passing in representation objects, the default representation is taken to be the class that is passed in.

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 axis (0) has special meaning:

  • axis=0 is the coordinate dimension (e.g., x, y, z for Cartesian)

So if the input position array, pos, has shape pos.shape = (3, 100), this would represent 100 3D positions (pos[0] is x, pos[1] is y, etc.). The same is true for velocity.

Parameters:
posrepresentation, quantity_like, or 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.

veldifferential, quantity_like, or 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.

frameFrameBase (optional)

The reference frame of the input phase-space positions.

copybool (optional)

If True, the input position and velocity data is copied. If False, the input data is referenced directly (if possible). Default is True.

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])

Align the circulation axis with the z-axis for tube orbits.

angular_momentum()

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

animate([components, units, stride, ...])

Animate an orbit or collection of orbits.

apocenter([return_times, func, approximate])

Estimate the apocenter(s) of the orbit by identifying local maxima in the spherical radius, fitting a parabola around these local maxima and then solving this parabola to find the apocenter(s).

circulation()

Determine which axes the orbit circulates around.

eccentricity(**kw)

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

energy([hamiltonian])

The total energy per unit mass:

estimate_period()

Estimate the period of the orbit in each dimension.

from_galpy_orbit(galpy_orbit)

Create a Gala PhaseSpacePosition or Orbit instance from a galpy.Orbit instance.

from_hdf5(f)

Load an object from an HDF5 file.

from_w(w[, units, copy])

Create a PhaseSpacePosition from a single array of positions and velocities.

get_components(which)

Get the component name dictionary for the desired object.

guiding_radius([potential, t])

Compute the guiding-center radius

kinetic_energy()

The kinetic energy per unit mass:

orbit_gen()

Generator for iterating over each orbit.

pericenter([return_times, func, approximate])

Estimate the pericenter(s) of the orbit by identifying local minima in the spherical radius, fitting a parabola around these local minima and then solving this parabola to find the pericenter(s).

plot([components, units, auto_aspect])

Plot the positions in all projections.

plot_3d([components, units, auto_aspect, ...])

Plot the specified 3D components.

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.

reshape(new_shape)

Reshape the underlying position and velocity arrays.

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])

Transform to a different reference frame.

to_galpy_orbit([ro, vo])

Convert this object to a galpy.Orbit instance.

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, approximate])

Estimate the maximum z height of the orbit by identifying local maxima in the absolute value of the z position, fitting a parabola around these local maxima and then solving this parabola to find the maximum z height.

Attributes Documentation

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

Align the circulation axis with the z-axis for tube orbits.

If the input orbit is a tube orbit (circulates around one axis), this method rotates the coordinate system so that circulation occurs around the z-axis. This is useful for standardizing orbit orientations when computing actions or comparing orbits.

Parameters:
circulationndarray, optional

The circulation pattern of the orbit. If not specified, this is computed using the circulation() method.

Returns:
orbOrbit

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}}

See Array shapes for more information about the shapes of input and output objects.

Returns:
LQuantity

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()
<Quantity [0.        ,0.        ,6.28318531] AU2 / yr>
animate(components=None, units=None, stride=1, segment_nsteps=10, underplot_full_orbit=True, show_time=True, marker_style=None, segment_style=None, FuncAnimation_kwargs=None, orbit_plot_kwargs=None, axes=None)[source]#

Animate an orbit or collection of orbits.

Parameters:
componentsiterable (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'].

unitsUnitBase, iterable, gala.units.UnitSystem (optional)

A single unit or list of units to display the components in.

strideint (optional)

How often to draw a new frame, in terms of orbit timesteps.

segment_nstepsint (optional)

How many timesteps to draw in an orbit segment trailing the timestep marker. Set this to 0 or None to disable.

underplot_full_orbitbool (optional)

Controls whether to under-plot the full orbit as a thin line.

show_timebool (optional)

Controls whether to show a label of the current timestep

marker_styledict or list of dict (optional)

Matplotlib style arguments passed to matplotlib.pyplot.plot that control the plot style of the timestep marker. If a single dict is passed then the marker_style is applied to all orbits. If a list of dicts is passed then each dict will be applied to each orbit.

segment_styledict (optional)

Matplotlib style arguments passed to matplotlib.pyplot.plot that control the plot style of the orbit segment. If a single dict is passed then the segment_style is applied to all orbits. If a list of dicts is passed then each dict will be applied to each orbit.

FuncAnimation_kwargsdict (optional)

Keyword arguments passed through to matplotlib.animation.FuncAnimation.

orbit_plot_kwargsdict (optional)

Keyword arguments passed through to gala.dynamics.Orbit.plot.

axesmatplotlib.axes.Axes (optional)

Where to draw the orbit.

Returns:
figmatplotlib.figure.Figure
animmatplotlib.animation.FuncAnimation
apocenter(return_times=False, func=<function mean>, approximate=False)[source]#

Estimate the apocenter(s) of the orbit by identifying local maxima in the spherical radius, fitting a parabola around these local maxima and then solving this parabola to find the apocenter(s).

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:
funcfunc (optional)

A function to evaluate on all of the identified apocenter times.

return_timesbool (optional)

Also return the apocenter times.

approximatebool (optional)

Compute an approximate apocenter by skipping interpolation.

Returns:
apofloat, ndarray

Either a single number or an array of apocenters.

timesndarray (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.

This method checks whether there is a change of sign of the angular momentum about each axis to determine circulation patterns. For example, a tube orbit will circulate around one axis, while a box orbit will not circulate around any axis.

Returns:
circulationndarray

An integer array indicating circulation about each axis. For each axis: 1 indicates circulation, 0 indicates no circulation. For a single orbit, returns a 1D array of shape (ndim,). For multiple orbits, returns a 2D array of shape (ndim, norbits).

Notes

This method works by checking whether the angular momentum about each axis changes sign or becomes very small during the orbit integration.

Examples

For a single 3D orbit:

  • Box and boxlet orbits: [0, 0, 0]

  • z-axis (short-axis) tube orbit: [0, 0, 1]

  • x-axis (long-axis) tube orbit: [1, 0, 0]

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.

Returns:
eccfloat

The orbital eccentricity.

energy(hamiltonian=None)[source]#

The total energy per unit mass:

Parameters:
hamiltoniangala.potential.Hamiltonian, gala.potential.PotentialBase instance

The Hamiltonian object to evaluate the energy. If a potential is passed in, this assumes a static reference frame.

Returns:
EQuantity

The total energy.

estimate_period()[source]#

Estimate the period of the orbit in each dimension.

Returns:
periodsQTable

The estimated orbital periods for each phase-space component, for each orbit.

Examples

>>> from gala.potential import MilkyWayPotential2022
>>> pot = MilkyWayPotential2022()

Compute an orbit and estimate the orbital period in each Cartesian component:

>>> orbit = pot.integrate_orbit([8., 0, 0, 0, 0.18, 0], dt=1., n_steps=4000)
>>> P_xyz = orbit.estimate_period()
>>> P_xyz
<QTable length=1>
        x                  y                  z
       Myr                Myr                Myr
     float64            float64            float64
------------------ ------------------ -----------------
176.02380952380952 176.07034632034632 56.43902691511387

Or, to estimate the period in cylindrical radius:

>>> orbit.cylindrical.estimate_period()["rho"]
<Quantity [121.09375] Myr>
classmethod from_galpy_orbit(galpy_orbit)[source]#

Create a Gala PhaseSpacePosition or Orbit instance from a galpy.Orbit instance.

Parameters:
galpy_orbitgalpy.orbit.Orbit
Returns:
orbitOrbit
classmethod from_hdf5(f)[source]#

Load an object from an HDF5 file.

Requires h5py.

Parameters:
fstr, h5py.File

Either the filename or an open HDF5 file.

classmethod from_w(w, units=None, copy=True, **kwargs)#

Create a PhaseSpacePosition from a single array of positions and velocities.

Parameters:
warray_like

The array of phase-space positions. Should have shape (2*ndim, ...) where the first ndim rows contain positions and the last ndim rows contain velocities.

unitsUnitSystem, optional

The unit system that the input position+velocity array, w, is represented in. If not provided, the array is assumed to be dimensionless.

copybool, optional

If True, the input array is copied. If False, the input data is referenced directly (if possible). Default is True.

**kwargs

Additional keyword arguments passed to the class initializer.

Returns:
objPhaseSpacePosition

A new PhaseSpacePosition instance created from the input array.

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

Can either be 'pos' or 'vel' to get the components for the position or velocity object.

guiding_radius(potential=None, t=0.0, **root_kwargs)[source]#

Compute the guiding-center radius

Parameters:
potentialgala.potential.PotentialBase subclass instance (optional)

The potential to compute the guiding radius in.

tquantity-like (optional)

Time.

**root_kwargs

Any additional keyword arguments are passed to root.

Returns:
RgQuantity

Guiding-center radius.

kinetic_energy()#

The kinetic energy per unit mass:

\[E_K = \frac{1}{2} \, |\boldsymbol{v}|^2\]
Returns:
EQuantity

The kinetic energy.

orbit_gen()[source]#

Generator for iterating over each orbit.

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

Estimate the pericenter(s) of the orbit by identifying local minima in the spherical radius, fitting a parabola around these local minima and then solving this parabola to find the pericenter(s).

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:
funcfunc (optional)

A function to evaluate on all of the identified pericenter times.

return_timesbool (optional)

Also return the pericenter times.

approximatebool (optional)

Compute an approximate pericenter by skipping interpolation.

Returns:
perifloat, ndarray

Either a single number or an array of pericenters.

timesndarray (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:
componentsiterable (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'].

unitsUnitBase, iterable, gala.units.UnitSystem (optional)

A single unit or list of units to display the components in.

auto_aspectbool (optional)

Automatically enforce an equal aspect ratio.

relative_tobool (optional)

Plot the values relative to this value or values.

autolimbool (optional)

Automatically set the plot limits to be something sensible.

axesarray_like (optional)

Array of matplotlib Axes objects.

subplots_kwargsdict (optional)

Dictionary of kwargs passed to subplots().

labelsiterable (optional)

List or iterable of axis labels as strings. They should correspond to the dimensions of the input orbit.

plot_functioncallable() (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.

Returns:
figFigure
plot_3d(components=None, units=None, auto_aspect=True, subplots_kwargs=None, **kwargs)[source]#

Plot the specified 3D components.

Parameters:
componentsiterable (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'].

unitsUnitBase, iterable, gala.units.UnitSystem (optional)

A single unit or list of units to display the components in.

auto_aspectbool (optional)

Automatically enforce an equal aspect ratio.

axmatplotlib.axes.Axes

The matplotlib Axes object to draw on.

subplots_kwargsdict (optional)

Dictionary of kwargs passed to subplots().

labelsiterable (optional)

List or iterable of axis labels as strings. They should correspond to the dimensions of the input orbit.

plot_functionstr (optional)

The matplotlib plot function to use. By default, this is ‘plot’ but can also be, e.g., ‘scatter’.

**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.

Returns:
figFigure
potential_energy(potential=None)[source]#

The potential energy per unit mass:

\[E_\Phi = \Phi(\boldsymbol{q})\]
Returns:
EQuantity

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_posBaseRepresentation

The type of representation to generate. Must be a class (not an instance), or the string name of the representation class.

new_velBaseDifferential (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.

Returns:
new_orbitgala.dynamics.Orbit
reshape(new_shape)[source]#

Reshape the underlying position and velocity arrays.

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

The frame instance specifying the desired output frame. For example, ICRS.

galactocentric_frameGalactocentric

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.

Returns:
cBaseCoordinateFrame

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

Transform to a different reference frame.

Parameters:
framegala.potential.CFrameBase

The frame to transform to.

current_framegala.potential.CFrameBase (optional)

If the Orbit has no associated Hamiltonian, this specifies the current frame of the orbit.

Returns:
orbitgala.dynamics.Orbit

The orbit in the new reference frame.

to_galpy_orbit(ro=None, vo=None)[source]#

Convert this object to a galpy.Orbit instance.

Parameters:
roastropy.units.Quantity or astropy.units.UnitBase

“Natural” length unit.

voastropy.units.Quantity or astropy.units.UnitBase

“Natural” velocity unit.

Returns:
galpy_orbitgalpy.orbit.Orbit
to_hdf5(f)[source]#

Serialize this object to an HDF5 file.

Requires h5py.

Parameters:
fstr, 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:
unitsUnitSystem (optional)

The unit system to represent the position and velocity in before combining into the full array.

Returns:
wndarray

A numpy array of all positions and velocities, without units. Will have shape (2*ndim, ...).

zmax(return_times=False, func=<function mean>, approximate=False)[source]#

Estimate the maximum z height of the orbit by identifying local maxima in the absolute value of the z position, fitting a parabola around these local maxima and then solving this parabola to find the maximum z height.

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:
funcfunc (optional)

A function to evaluate on all of the identified z maximum times.

return_timesbool (optional)

Also return the times of maximum.

approximatebool (optional)

Compute approximate values by skipping interpolation.

Returns:
zsfloat, ndarray

Either a single number or an array of maximum z heights.

timesndarray (optional, see return_times)

If return_times=True, also returns an array of the apocenter times.