.. include:: references.txt .. _orbits-in-detail: ***************************************************** Orbit and phase-space position objects in more detail ***************************************************** For the examples below the following imports have already been executed:: >>> import astropy.units as u >>> import numpy as np >>> import gala.potential as gp >>> import gala.dynamics as gd >>> from astropy.coordinates import (CylindricalRepresentation, ... CylindricalDifferential) >>> from gala.units import galactic >>> np.random.seed(42) We will also set the default Astropy Galactocentric frame parameters to the values adopted in Astropy v4.0: >>> import astropy.coordinates as coord >>> _ = coord.galactocentric_frame_defaults.set('v4.0') Introduction ============ The `astropy.units` subpackage is excellent for working with numbers and associated units, but dynamical quantities often contain many quantities with mixed units. An example is a position in phase-space, which may contain some quantities with length units and some quantities with velocity or momentum units. The |psp| and |orb| classes are designed to work with these data structures and provide a consistent API for visualizing and computing further dynamical quantities. Click these shortcuts to jump to a section below, or start reading below: * :ref:`phase-space-position` * :ref:`orbit` .. _phase-space-position: Phase-space Positions ===================== The |psp| class provides an interface for representing full phase-space positions--coordinate positions and momenta (velocities). This class is useful as a container for initial conditions and for transforming phase-space positions to new coordinate representations or reference frames. The easiest way to create a |psp| object is to pass in a pair of `~astropy.units.Quantity` objects that represent the Cartesian position and velocity vectors:: >>> gd.PhaseSpacePosition(pos=[4., 8., 15.] * u.kpc, ... vel=[-150., 50., 15.] * u.km/u.s) By default, passing in `~astropy.units.Quantity`'s are interpreted as Cartesian coordinates and velocities. This works with arrays of positions and velocities as well:: >>> x = np.arange(24).reshape(3, 8) >>> v = np.arange(24).reshape(3, 8) >>> w = gd.PhaseSpacePosition(pos=x * u.kpc, ... vel=v * u.km/u.s) >>> w This is interpreted as 8, 6-dimensional phase-space positions. The class internally stores the positions and velocities as `~astropy.coordinates.BaseRepresentation` and `~astropy.coordinates.BaseDifferential` subclasses; in this case, `~astropy.coordinates.CartesianRepresentation` and `~astropy.coordinates.CartesianDifferential`:: >>> w.pos >>> w.vel All of the components of these classes are mapped to attributes of the phase-space position class for convenience, but with more user-friendly names. These mappings are defined in the class definition of `~gala.dynamics.PhaseSpacePosition`. For example, to access the ``x`` component of the position and the ``v_x`` component of the velocity:: >>> w.x # doctest: +FLOAT_CMP >>> w.v_x # doctest: +FLOAT_CMP The default representation is Cartesian, but the class can also be instantiated with representation objects instead of `~astropy.units.Quantity`'s -- this is useful for creating |psp| or |orb| instances from non-Cartesian representations of the position and velocity:: >>> pos = CylindricalRepresentation(rho=np.linspace(1., 4, 4) * u.kpc, ... phi=np.linspace(0, np.pi, 4) * u.rad, ... z=np.linspace(-1, 1., 4) * u.kpc) >>> vel = CylindricalDifferential(d_rho=np.linspace(100, 150, 4) * u.km/u.s, ... d_phi=np.linspace(-1, 1, 4) * u.rad/u.Myr, ... d_z=np.linspace(-15, 15., 4) * u.km/u.s) >>> w = gd.PhaseSpacePosition(pos=pos, vel=vel) >>> w >>> w.rho We can easily transform the full phase-space vector to new representations or coordinate frames. These transformations use the :mod:`astropy.coordinates` |astropyrep|_:: >>> cart = w.represent_as('cartesian') >>> cart.x >>> sph = w.represent_as('spherical') >>> sph.distance There is also support for transforming the positions and velocities (assumed to be in a `~astropy.coordinates.Galactocentric` frame) to any of the other coordinate frames. For example, to transform to :class:`~astropy.coordinates.Galactic` coordinates:: >>> from astropy.coordinates import Galactic >>> gal_c = w.to_coord_frame(Galactic) >>> gal_c # doctest: +FLOAT_CMP We can easily plot projections of the phase-space positions using the `~gala.dynamics.PhaseSpacePosition.plot` method:: >>> np.random.seed(42) >>> x = np.random.uniform(-10, 10, size=(3,128)) >>> v = np.random.uniform(-200, 200, size=(3,128)) >>> w = gd.PhaseSpacePosition(pos=x * u.kpc, ... vel=v * u.km/u.s) >>> fig = w.plot() # doctest: +SKIP .. plot:: :align: center :context: close-figs import astropy.units as u import numpy as np import gala.dynamics as gd np.random.seed(42) x = np.random.uniform(-10,10,size=(3,128)) v = np.random.uniform(-200,200,size=(3,128)) w = gd.PhaseSpacePosition(pos=x*u.kpc, vel=v*u.km/u.s) fig = w.plot() This is a thin wrapper around the `~gala.dynamics.plot_projections` function and any keyword arguments are passed through to that function:: >>> fig = w.plot(components=['x', 'v_z'], color='r', ... facecolor='none', marker='o', s=20, alpha=0.5) # doctest: +SKIP .. plot:: :align: center :context: close-figs :width: 60% fig = w.plot(components=['x', 'v_z'], color='r', facecolor='none', marker='o', s=20, alpha=0.5) .. _orbit: Orbits ====== The |orb| class inherits much of the functionality from |psp| (described above) and adds some additional features that are useful for time-series orbits. An |orb| instance is initialized like the |psp|--with arrays of positions and velocities-- but usually also requires specifying a time array as well. Also, the extra axes in these arrays hold special meaning for the |orb| class. The position and velocity arrays passed to |psp| can have arbitrary numbers of dimensions as long as the 0th axis specifies the dimensionality. For the |orb| class, the 0th axis remains the axis of dimensionality, but the 1st axis now is always assumed to be the time axis. For example, an input position with shape ``(2,128)`` to a |psp| represents 128 independent 2D positions, but to a |orb| it represents a single orbit's positions at 128 times:: >>> t = np.linspace(0, 100, 128) * u.Myr >>> Om = 1E-1 * u.rad / u.Myr >>> pos = np.vstack((5*np.cos(Om*t), np.sin(Om*t))).value * u.kpc >>> vel = np.vstack((-5*np.sin(Om*t), np.cos(Om*t))).value * u.kpc/u.Myr >>> orbit = gd.Orbit(pos=pos, vel=vel) >>> orbit To create a single object that contains multiple orbits, the input position object should have 3 axes. The last axis (``axis=2``) specifies the number of orbits. So, an input position with shape ``(2,128,16)`` would represent 16, 2D orbits, each with the same 128 times:: >>> t = np.linspace(0, 100, 128) * u.Myr >>> Om = np.random.uniform(size=16) * u.rad / u.Myr >>> angle = Om[None] * t[:, None] >>> pos = np.stack((5*np.cos(angle), np.sin(angle))).value * u.kpc >>> vel = np.stack((-5*np.sin(angle), np.cos(angle))).value * u.kpc/u.Myr >>> orbit = gd.Orbit(pos=pos, vel=vel) >>> orbit To make full use of the orbit functionality, you must also pass in an array with the time values and an instance of a `~gala.potential.potential.PotentialBase` subclass that represents the potential that the orbit was integrated in:: >>> pot = gp.PlummerPotential(m=1E10, b=1., units=galactic) >>> orbit = gd.Orbit(pos=pos*u.kpc, vel=vel*u.km/u.s, ... t=t*u.Myr, potential=pot) (note, in this case ``pos`` and ``vel`` were not generated from integrating an orbit in the potential ``pot``!). However, most of the time you won't need to create |orb| objects from scratch! They are returned from any of the numerical integration routines provided in `gala`. For example, they are returned by the `~gala.potential.potential.PotentialBase.integrate_orbit` method of potential objects and will automatically contain the ``time`` array and ``potential`` object. For example:: >>> pot = gp.PlummerPotential(m=1E10 * u.Msun, b=1. * u.kpc, units=galactic) >>> w0 = gd.PhaseSpacePosition(pos=[10.,0,0] * u.kpc, ... vel=[0.,75,0] * u.km/u.s) >>> orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=1., n_steps=5000) >>> orbit >>> orbit.t >>> orbit.potential Just like for |psp|, we can quickly visualize an orbit using the `~gala.dynamics.Orbit.plot` method:: >>> fig = orbit.plot() # doctest: +SKIP .. plot:: :align: center :context: close-figs import astropy.units as u import gala.dynamics as gd import gala.potential as gp from gala.units import galactic pot = gp.PlummerPotential(m=1E10 * u.Msun, b=1. * u.kpc, units=galactic) w0 = gd.PhaseSpacePosition(pos=[2.,0,0] * u.kpc, vel=[0.,75,15] * u.km/u.s) orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=1., n_steps=5000) fig = orbit.plot() Again, this is a thin wrapper around the `~gala.dynamics.plot_projections` function and any keyword arguments are passed through to that function:: >>> fig = orbit.plot(linewidth=4., alpha=0.5, color='r') # doctest: +SKIP .. plot:: :align: center :context: close-figs fig = orbit.plot(linewidth=4., alpha=0.5, color='r') Alternatively, for three-dimensional orbits, we can visualize the orbit using the 3D projection capabilities in `matplotlib`:: >>> fig = orbit.plot_3d(alpha=0.5, color='k') # doctest: +SKIP .. plot:: :align: center :context: close-figs :width: 60% fig = orbit.plot_3d(alpha=0.5, color='k') We can also quickly create an animation of the progression of an orbit using the `~gala.dynamics.Orbit.animate` method, which animated projections of the orbit:: >>> fig, anim = orbit[:1000].animate(stride=10) # doctest: +SKIP .. raw:: html The animate method acts like `~gala.dynamics.Orbit.plot`, in that it works for any coordinate representation (Cartesian, cylindrical, etc.) and supports only animating subsets of the phase-space components. For example, to make an animation of an orbit in cylindrical coordinates, showing the orbit proress in the R,z meridional plane:: >>> fig, anim = orbit[:1000].cylindrical.animate(components=['rho', 'z'], # doctest: +SKIP ... stride=10) .. raw:: html We can also quickly compute quantities like the angular momentum, and estimates for the pericenter, apocenter, eccentricity of the orbit. Estimates for the latter few get better with smaller timesteps:: >>> orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.1, n_steps=100000) >>> np.mean(orbit.angular_momentum(), axis=1) # doctest: +FLOAT_CMP >>> orbit.eccentricity() # doctest: +FLOAT_CMP >>> orbit.pericenter() # doctest: +FLOAT_CMP >>> orbit.apocenter() # doctest: +FLOAT_CMP More information ================ Internally, both of the above classes rely on the Astropy representation transformation framework (i.e. the subclasses of `~astropy.coordinates.BaseRepresentation` and `~astropy.coordinates.BaseDifferential`). However, at present these classes only support 3D positions and differentials (velocities). The |psp| and |orb| classes both support arbitrary numbers of dimensions and, when relevant, rely on custom subclasses of the representation classes to handle such cases. See the :ref:`nd-representations` page for more information about these classes.