[WIP] Simplify PhaseSpacePosition and Orbit classes

.rtd-environment.yml

@@ -4,7 +4,7 @@ channels:
   - astropy
 
 dependencies:
-  - astropy
+  - astropy>=2
   - astropy-helpers
   - numpy
   - cython

CHANGES.md

@@ -1,27 +1,41 @@
 0.2 (unreleased)
-----------------
+================
 
+New Features
+------------
 - Added a new potential class for the Satoh density (Satoh 1980).
 - Added support for Leapfrog integration when generating mock stellar streams.
 - Added new colormaps and defaults for the matplotlib style.
 - Added support for non-inertial reference frames and implemented a constant
   rotating reference frame.
-- Added a new class - `Hamiltonian` - for storing potentials with reference
-  frames and should be used for orbit integration.
-- Added a new mock stream argument to output orbits of all of the mock stream
-  star particles to an HDF5 file.
+- Added a new class - ``Hamiltonian`` - for storing potentials with reference
+  frames. This should be used for easy orbit integration instead of the
+  potential classes.
+- Added a new argument to the mock stream generation functions t output orbits
+  of all of the mock stream star particles to an HDF5 file.
 - Cleaned up and simplified the process of subclassing a C-implemented
   gravitational potential.
 - Gravitational potential class instances can now be composed by just adding the
   instances.
 
+API-breaking changes
+--------------------
+- ``CartesianPhaseSpacePosition`` and ``CartesianOrbit`` are deprecated. Use
+  ``PhaseSpacePosition`` and ``Orbit`` with a Cartesian representation instead.
+- Overhauled the storage of position and velocity information on
+  ``PhaseSpacePosition`` and ``Orbit`` classes. This uses new features in
+  Astropy 2.0 that allow attaching "differential" classes to representation
+  classes for storing velocity information. ``.pos`` and ``.vel`` no longer
+  point to arrays of Cartesian coordinates, but now instead point to
+  astropy.coordinates representation and differential objects, respectively.
+
 0.1.1 (2016-05-20)
-------------------
+==================
 
 - Removed debug statement.
 - Added 'Why' page to documentation.
 
 0.1.0 (2016-05-19)
-------------------
+==================
 
 - Initial release.

docs/coordinates/index.rst

@@ -7,10 +7,13 @@ Coordinate Systems (`gala.coordinates`)
 Introduction
 ============
 
-The `~gala.coordinates` subpackage provides functions for converting
-velocities between various astronomical coordinate frames and systems.
-This subpackage also provides :mod:`astropy.coordinates` frame classes
-for coordinate sytems defined by the Sagittarius and Orphan streams.
+The `~gala.coordinates` subpackage primarily provides specialty
+:mod:`astropy.coordinates` frame classes for coordinate systems defined by the
+Sagittarius, Orphan, and GD1 streams. It also contains functions for converting
+velocities between various astronomical coordinate frames and systems, but these
+are now deprecated with v0.2 because Astropy v2.0 `now supports velocity
+transformations
+<http://docs.astropy.org/en/stable/coordinates/velocities.html>`_.
 
 For the examples below the following imports have already been executed::
 
@@ -19,97 +22,184 @@ For the examples below the following imports have already been executed::
     >>> import astropy.units as u
     >>> import gala.coordinates as gc
 
-Getting Started
-===============
-
-The core functions in this subpackage provide support to:
+Stellar stream coordinate frames
+================================
+
+`Gala` provides Astropy coordinate frame classes for
+transforming to Sagittarius, Orphan, and GD1 stream coordinates (as defined in the
+references below). These classes behave like the built-in astropy coordinates
+frames (e.g., :class:`~astropy.coordinates.ICRS` or
+:class:`~astropy.coordinates.Galactic`) and can be transformed to and from
+other astropy coordinate frames. For example, to convert a set of
+`~astropy.coordinates.ICRS` (RA, Dec) coordinates to the
+`~gala.coordinates.Sagittarius` frame::
+
+    >>> c = coord.ICRS(ra=100.68458*u.degree, dec=41.26917*u.degree)
+    >>> sgr = c.transform_to(gc.Sagittarius)
+    >>> (sgr.Lambda, sgr.Beta) # doctest: +FLOAT_CMP
+    (<Longitude 179.58511053544734 deg>, <Latitude -12.558450192162654 deg>)
+
+Or, to transform from `~gala.coordinates.Sagittarius` coordinates to the
+`~astropy.coordinates.Galactic` frame::
+
+    >>> sgr = gc.Sagittarius(Lambda=156.342*u.degree, Beta=1.1*u.degree)
+    >>> c = sgr.transform_to(coord.Galactic)
+    >>> (c.l, c.b) # doctest: +FLOAT_CMP
+    (<Longitude 182.5922090437946 deg>, <Latitude -9.539692094685893 deg>)
+
+These transformations also handle velocities so that proper motion components
+can be transformed between the systems. For example, to transform from
+`~gala.coordinates.GD1` proper motions to `~astropy.coordinates.Galactic` proper
+motions::
+
+    >>> gd1 = gc.GD1(phi1=-35.00*u.degree, phi2=0*u.degree,
+    ...              pm_phi1_cosphi2=-12.20*u.mas/u.yr,
+    ...              pm_phi2=-3.10*u.mas/u.yr)
+    >>> gd1.transform_to(coord.Galactic) # doctest: +FLOAT_CMP
+    <Galactic Coordinate: (l, b) in deg
+        ( 181.28968151,  54.84972806)
+     (pm_l_cosb, pm_b) in mas / yr
+        ( 12.03209393, -3.69847479)>
+
+As with the other Astropy coordinate frames, with a full specification of the 3D
+position and velocity, we can transform to a
+`~astropy.coordinates.Galactocentric` frame::
+
+    >>> gd1 = gc.GD1(phi1=-35.00*u.degree, phi2=0.04*u.degree,
+    ...              distance=7.83*u.kpc,
+    ...              pm_phi1_cosphi2=-12.20*u.mas/u.yr,
+    ...              pm_phi2=-3.10*u.mas/u.yr,
+    ...              radial_velocity=-32*u.km/u.s)
+    >>> gd1.transform_to(coord.Galactocentric) # doctest: +FLOAT_CMP
+    <Galactocentric Coordinate (galcen_coord=<ICRS Coordinate: (ra, dec) in deg
+        ( 266.4051, -28.936175)>, galcen_distance=8.3 kpc, galcen_v_sun=( 11.1,  232.24,  7.25) km / s, z_sun=27.0 pc, roll=0.0 deg): (x, y, z) in kpc
+        (-12.78977138, -0.09870921,  6.44110283)
+     (v_x, v_y, v_z) in km / s
+        (-73.01933674, -216.37648654, -97.60065189)>
+
+References
+----------
+
+* `A 2MASS All-Sky View of the Sagittarius Dwarf Galaxy: I. Morphology of the
+  Sagittarius Core and Tidal Arms <http://arxiv.org/abs/astro-ph/0304198>`_
+* `The Orbit of the Orphan Stream <http://arxiv.org/abs/1001.0576>`_
+* `Constraining the Milky Way potential with a 6-D phase-space map of the GD-1
+  stellar stream <https://arxiv.org/abs/0907.1085>`_
+
+Moving from ``gala.coordinates`` to ``astropy.coordinates``
+===========================================================
+
+``Gala`` previously supported the following transformations:
 
 - Convert Galactocentric, cartesian velocities to heliocentric proper motion
   and radial velocities.
 - Convert proper motions and radial velocities to Galactocentric, cartesian
   velocities.
 - Convert proper motions from/to ICRS to/from Galactic.
-- Convert radial velocities from/to the Galactic Standard of Rest (GSR) to/from a barycentric frame.
-
-These functions work naturally with the :mod:`astropy.units` and
-:mod:`astropy.coordinates` subpackages. Handling positional transformations
-is already supported by :mod:`astropy.coordinates` and new to Astropy v1.0 is a
-:class:`~astropy.coordinates.Galactocentric` reference frame. However, there is
-currently no support for transforming velocities in Astropy. The functions below
-attempt to bridge that gap as a temporary solution until support is added
-(planned for v1.2).
-
-For example, to convert a spherical, heliocentric velocity (proper motion and radial
-velocity) in an ICRS frame to a Galactocentric, cartesian velocity, we first have
-to define an Astropy coordinate to specify the position of the object::
-
-    >>> c = coord.SkyCoord(ra=100.68458*u.deg, dec=41.26917*u.deg, distance=1.1*u.kpc)
-
-Then pass this object in to the heliocentric to galactocentric conversion
-function, :func:`~gala.coordinates.vhel_to_gal`::
-
-    >>> pm = [1.5, -1.7] * u.mas/u.yr
-    >>> rv = 151.1 * u.km/u.s
-    >>> gc.vhel_to_gal(c.icrs, pm=pm, rv=rv)
-    <Quantity [-134.44094022, 228.42957796,  52.97041271] km / s>
-
-Because the input coordinate is given in the ICRS frame, the function assumes that
-the proper motion is also in this frame, e.g., that the proper motion components are
-:math:`(\mu_\alpha\cos\delta, \mu_\delta)`. If we instead passed in a coordinate in
-the Galactic frame, the components are assumed to be :math:`(\mu_l\cos b, \mu_b)` ::
-
-    >>> gc.vhel_to_gal(c.galactic, pm=pm, rv=rv)
-    <Quantity [-137.63839061, 232.10966761,  40.73818895] km / s>
-
-The velocity transformation functions allow specifying the circular velocity at the Sun
-(``vcirc``) and a 3-vector specifying the Sun's velocity with respect to the local
-standard of rest (``vlsr``). Further customization of the Sun's location can be made via
-the :class:`~astropy.coordinates.Galactocentric` frame attributes and passed in with the
-keyword argument ``galactocentric_frame`` ::
-
-    >>> frame = coord.Galactocentric(z_sun=10.*u.pc, galcen_distance=8.3*u.kpc)
-    >>> gc.vhel_to_gal(c.icrs, pm=pm, rv=rv, galactocentric_frame=frame,
-    ...                vcirc=218*u.km/u.s, vlsr=[0.,0.,0.]*u.km/u.s)
-    <Quantity [-144.5344455 , 221.17957796,  45.50447318] km / s>
-
-The inverse transformations are also available, with the function
-:func:`~gala.coordinates.vgal_to_hel`. Here, because the input coordinate is passed
-in after being transformed to the ICRS frame, the output proper motions will also be
-given in the ICRS frame :math:`(\mu_\alpha\cos\delta, \mu_\delta)`::
-
-    >>> xyz = coord.Galactocentric([11., 15, 25] * u.kpc)
-    >>> vxyz = [121., 150., -75.] * u.km/u.s
-    >>> pm_ra,pm_dec,vr = gc.vgal_to_hel(xyz.transform_to(coord.ICRS), vxyz)
-    >>> pm_ra # doctest: +FLOAT_CMP
-    <Quantity 0.2834641666390529 mas / yr>
-    >>> pm_dec # doctest: +FLOAT_CMP
-    <Quantity -0.888174413651107 mas / yr>
-    >>> vr # doctest: +FLOAT_CMP
-    <Quantity -29.71790624810498 km / s>
-
-Passing in coordinates in the Galactic frame means that the output proper motions will
-instead be :math:`(\mu_l\cos b, \mu_b)` ::
-
-    >>> pm_l,pm_b,vr = gc.vgal_to_hel(xyz.transform_to(coord.Galactic), vxyz)
-    >>> pm_l # doctest: +FLOAT_CMP
-    <Quantity -0.7713637315333076 mas / yr>
-
-All of these functions also work on arrays of coordinates and velocities, e.g.::
-
-    >>> xyz = coord.Galactocentric(np.random.uniform(-20,20,size=(3,10)) * u.kpc)
-    >>> vxyz = np.random.uniform(-150,150,size=(3,10)) * u.km/u.s
-    >>> gc.vgal_to_hel(xyz.transform_to(coord.ICRS), vxyz) # doctest: +SKIP
-    ...
-
-Using gala.coordinates
-======================
-More details are provided in the linked pages below:
-
-.. toctree::
-   :maxdepth: 1
-
-   propermotion
-   streamframes
+- Convert radial velocities from/to the Galactic Standard of Rest (GSR) to/from
+  a barycentric frame.
+
+Below, we'll show examples of how to do each of these transformations using
+``astropy.coordinates`` instead.
+
+Convert Galactocentric (or simulated) to Heliocentric coordinates and velocities
+--------------------------------------------------------------------------------
+
+Let's assume we have a set of Cartesian positions and velocities that we
+assume are in a Galactocentric (i.e. Milky Way-centric) frame::
+
+    >>> xyz = [[  6.47945349, -34.4772621 ],
+    ...        [-17.77019357,  31.11681441],
+    ...        [-25.32101661,   3.54273331]] * u.kpc
+    >>> vxyz = [[-184.32952533,   76.94467146],
+    ...         [ -19.19375771,  -98.15603007],
+    ...         [  22.49696323,   26.35110153]] * u.km/u.s
+
+These could be, for example, the output from a simulation or sampling from a
+distribution function. We want to transform these coordinates to a Helio- or
+Barycentric coordinate frame to compute observable quantities like proper
+motions and radial velocities. This is now supported in Astropy v2.0. To start,
+we have to define the Galactocentric frame. That is, we have to define the sun's
+position and velocity within the assumed Galactocentric frame. We do this by
+setting frame attributes of the `astropy.coordinates.Galactocentric` frame. This
+is a right-handed coordinate system with defaults for the position of the
+Galactic center in ICRS coordinates, the sun-galactic-center distance, the
+height of the sun above the Galactic midplane, and the solar velocity vector. We
+can modify all of these parameters, but for the sake of example we'll just
+change the distance and solar velocity (full velocity of the sun, including the
+circular velocity and motion with respect to the local standard of rest)::
+
+    >>> v_sun = coord.CartesianDifferential([10, 250., 7] * u.km/u.s)
+    >>> gc = coord.Galactocentric(x=xyz[0], y=xyz[1], z=xyz[2],
+    ...                           v_x=vxyz[0], v_y=vxyz[1], v_z=vxyz[2],
+    ...                           galcen_distance=8*u.kpc,
+    ...                           galcen_v_sun=v_sun)
+
+To transform to any other `astropy.coordinates` frame, we use the
+``transform_to()`` method::
+
+    >>> icrs = gc.transform_to(coord.ICRS)
+
+From the new frame, we can access the sky positions, Barycentric distances,
+proper motion components, and Barycentric radial velocity. For example::
+
+    >>> icrs.distance # doctest: +FLOAT_CMP
+    <Distance [ 34.17535632, 41.00810092] kpc>
+    >>> icrs.pm_ra_cosdec # doctest: +FLOAT_CMP
+    <Quantity [ 1.78185606, 0.92782761] mas / yr>
+    >>> icrs.radial_velocity # doctest: +FLOAT_CMP
+    <Quantity [  46.14528132,-305.74417355] km / s>
+
+Convert Heliocentric (observed) to Galactocentric coordinates and velocities
+----------------------------------------------------------------------------
+
+To transform from Heliocentric coordinates to Galactocentric coordinates, we
+also have to define the sun's position and velocity within the assumed
+Galactocentric frame. For the example below, we'll use the same frame parameter
+values as above, but note that here we don't pass data in to the frame class::
+
+    >>> v_sun = coord.CartesianDifferential([10, 250., 7] * u.km/u.s)
+    >>> gc = coord.Galactocentric(galcen_distance=8*u.kpc,
+    ...                           galcen_v_sun=v_sun)
+
+Now we need some Heliocentric coordinates to transform. The Galactocentric
+transformation requires full 3D position and velocity information, so we'll
+have to specify a sky position, distance, proper motion components, and a
+radial velocity. Let's start with coordinates in the ICRS frame::
+
+    >>> icrs = coord.ICRS(ra=11.23*u.degree, dec=58.13*u.degree,
+    ...                   distance=21.34*u.pc,
+    ...                   pm_ra_cosdec=-55.89*u.mas/u.yr, pm_dec=71*u.mas/u.yr,
+    ...                   radial_velocity=210.43*u.km/u.s)
+
+We again use the ``transform_to()`` method to do the transformation::
+
+    >>> icrs.transform_to(gc) # doctest: +FLOAT_CMP
+    <Galactocentric Coordinate (galcen_coord=<ICRS Coordinate: (ra, dec) in deg
+        ( 266.4051, -28.936175)>, galcen_distance=8.0 kpc, galcen_v_sun=( 10.,  250.,  7.) km / s, z_sun=27.0 pc, roll=0.0 deg): (x, y, z) in pc
+        (-8011.25186648,  18.02217595,  25.27812036)
+     (v_x, v_y, v_z) in km / s
+        (-97.06296832,  431.12942421, -2.69495881)>
+
+Convert proper motions between the ``ICRS`` and ``Galactic`` frames
+-------------------------------------------------------------------
+
+The above Galactocentric coordinate transformations require full 3D position and
+velocity information. However, transforming proper motion components between
+different Barycentric coordinate frames is just a rotation, and can therefore be
+done with just sky position and proper motions. For example, to convert from
+ICRS proper motions to Galactic proper motions::
+
+    >>> icrs = coord.ICRS(ra=11.23*u.degree, dec=58.13*u.degree,
+    ...                   pm_ra_cosdec=-55.89*u.mas/u.yr, pm_dec=71*u.mas/u.yr)
+    >>> gal = icrs.transform_to(coord.Galactic)
+    >>> gal # doctest: +FLOAT_CMP
+    <Galactic Coordinate: (l, b) in deg
+        ( 122.06871373, -4.73082278)
+     (pm_l_cosb, pm_b) in mas / yr
+        (-54.0689922,  72.39638239)>
+    >>> gal.pm_l_cosb # doctest: +FLOAT_CMP
+    <Quantity -54.06899219513397 mas / yr>
 
 .. _gala-coordinates-api:
 
@@ -118,3 +208,11 @@ API
 
 .. automodapi:: gala.coordinates
     :no-inheritance-diagram:
+    :skip: pm_gal_to_icrs
+    :skip: pm_icrs_to_gal
+    :skip: transform_proper_motion
+    :skip: vgal_to_hel
+    :skip: vgsr_to_vhel
+    :skip: vhel_to_gal
+    :skip: vhel_to_vgsr
+    :skip: Quaternion