Merge 220ffb9 into 2c65d9e

.travis.yml

@@ -23,9 +23,11 @@ env:
         # overidden underneath. They are defined here in order to save having
         # to repeat them for all configurations.
         - NUMPY_VERSION=stable
-        - ASTROPY_VERSION=stable
+        # - ASTROPY_VERSION=stable # HACK: to use mhvk's branch
         - SETUP_CMD='test'
         - CONDA_DEPENDENCIES='cython jinja2 scipy matplotlib pyyaml h5py sympy qt'
+        # HACK: to use mhvk's branch
+        - PIP_DEPENDENCIES='git+https://github.com/mhvk/astropy.git@representation-differentials'
 
     matrix:
         # Make sure that egg_info works without dependencies

docs/coordinates/index.rst

@@ -29,7 +29,8 @@ The core functions in this subpackage provide support to:
 - 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.
+- 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
@@ -39,38 +40,45 @@ 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::
+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)
+    >>> 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`` ::
+    >>> vhel = coord.SphericalDifferential(d_lon=1.5 * u.mas/u.yr,
+    ...                                    d_lat=-1.7 * u.mas/u.yr,
+    ...                                    d_distance=151.1 * u.km/u.s)
+    >>> gc.vhel_to_gal(c.icrs, vhel)
+    <CartesianDifferential (d_x, d_y, d_z) in km / s
+        (-134.86592639,  229.212658,  51.24398172)>
+
+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, \mu_\delta)`. If we instead passed in a
+coordinate in the Galactic frame, the components are assumed to be
+:math:`(\mu_l, \mu_b)` ::
+
+    >>> gc.vhel_to_gal(c.galactic, vhel)
+    <CartesianDifferential (d_x, d_y, d_z) in km / s
+        (-137.60840819,  232.4106025,  40.73809195)>
+
+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,
+    >>> gc.vhel_to_gal(c.icrs, vhel, 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>
+    <CartesianDifferential (d_x, d_y, d_z) in km / s
+        (-144.95589471,  221.962658,  43.77717535)>
 
 The inverse transformations are also available, with the function
 :func:`~gala.coordinates.vgal_to_hel`. Here, because the input coordinate is passed
@@ -79,27 +87,21 @@ 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>
+    >>> c = xyz.transform_to(coord.ICRS)
+    >>> gc.vgal_to_hel(c, vxyz) # doctest: +FLOAT_CMP
+    <SphericalDifferential (d_lon, d_lat, d_distance) in (mas / yr, mas / yr, km / s)
+        ( 0.30661983, -0.88817441, -29.71790625)>
 
 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>
+instead be :math:`(\mu_l, \mu_b)` ::
 
-All of these functions also work on arrays of coordinates and velocities, e.g.::
+    >>> c = xyz.transform_to(coord.Galactic)
+    >>> gc.vgal_to_hel(c, vxyz) # doctest: +FLOAT_CMP
+    <SphericalDifferential (d_lon, d_lat, d_distance) in (mas / yr, mas / yr, km / s)
+        (-2.58898658,  0.07942092, -137.0530746)>
 
-    >>> 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
-    ...
+All of these functions also work on arrays of many coordinates and velocities as
+well.
 
 Using gala.coordinates
 ======================

docs/dynamics/actionangle.rst

@@ -37,12 +37,15 @@ been executed::
 
     >>> import astropy.coordinates as coord
     >>> import astropy.units as u
-    >>> import matplotlib.pyplot as pl
+    >>> import matplotlib.pyplot as plt
     >>> import numpy as np
     >>> import gala.potential as gp
     >>> import gala.dynamics as gd
     >>> from gala.units import galactic
 
+For many more options for action calculation, see
+`tact <https://github.com/jls713/tact>`_.
+
 .. _tube-axisymmetric:
 
 A tube orbit in an axisymmetric potential
@@ -74,42 +77,34 @@ We first create a potential and set up our initial conditions::
 
     >>> pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, q1=1., q2=1., q3=0.9, r_h=0,
     ...                               units=galactic)
-    >>> w0 = gd.CartesianPhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
-    ...                                     vel=[75, 150, 50.]*u.km/u.s)
+    >>> w0 = gd.PhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
+    ...                            vel=[75, 150, 50.]*u.km/u.s)
 
 We will now integrate the orbit and plot it in the meridional plane::
 
     >>> w = pot.integrate_orbit(w0, dt=0.5, n_steps=10000)
-    >>> cyl_pos, cyl_vel = w.represent_as(coord.CylindricalRepresentation)
-    >>> fig,ax = pl.subplots(1,1,figsize=(6,6))
-    >>> ax.plot(cyl_pos.rho.to(u.kpc), cyl_pos.z.to(u.kpc),
-    ...         marker=None, linestyle='-') # doctest: +SKIP
-    >>> ax.set_xlabel("R [kpc]") # doctest: +SKIP
-    >>> ax.set_ylabel("z [kpc]") # doctest: +SKIP
+    >>> cyl = w.represent_as('cylindrical')
+    >>> fig = cyl.plot(['rho', 'z'], linestyle='-')
 
 .. plot::
     :align: center
 
     import astropy.coordinates as coord
     import astropy.units as u
-    import matplotlib.pyplot as pl
+    import matplotlib.pyplot as plt
     import numpy as np
     import gala.potential as gp
     import gala.dynamics as gd
     from gala.units import galactic
 
     pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, q1=1., q2=1., q3=0.9, r_h=0,
                                   units=galactic)
-    w0 = gd.CartesianPhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
-                                        vel=[75, 150, 50.]*u.km/u.s)
+    w0 = gd.PhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
+                               vel=[75, 150, 50.]*u.km/u.s)
 
     w = pot.integrate_orbit(w0, dt=0.5, n_steps=10000)
-    cyl_pos, cyl_vel = w.represent_as(coord.CylindricalRepresentation)
-    fig,ax = pl.subplots(1,1,figsize=(6,6))
-    ax.plot(cyl_pos.rho.to(u.kpc).value, cyl_pos.z.to(u.kpc).value,
-            marker=None, linestyle='-')
-    ax.set_xlabel("R [kpc]")
-    ax.set_ylabel("z [kpc]")
+    cyl = w.represent_as('cylindrical')
+    cyl.plot(['rho', 'z'], linestyle='-')
 
 To solve for the actions in the true potential, we first compute the actions in
 a "toy" potential -- a potential in which we can compute the actions and angles
@@ -132,7 +127,7 @@ angles would be perfectly straight lines with slope equal to the frequencies.
 Instead, the orbit is wobbly in the toy potential angles::
 
     >>> toy_actions,toy_angles,toy_freqs = toy_potential.action_angle(w)
-    >>> fig,ax = pl.subplots(1,1,figsize=(5,5))
+    >>> fig,ax = plt.subplots(1,1,figsize=(5,5))
     >>> ax.plot(toy_angles[0], toy_angles[2], linestyle='none', marker=',') # doctest: +SKIP
     >>> ax.set_xlim(0,2*np.pi) # doctest: +SKIP
     >>> ax.set_ylim(0,2*np.pi) # doctest: +SKIP
@@ -144,21 +139,21 @@ Instead, the orbit is wobbly in the toy potential angles::
 
     import astropy.coordinates as coord
     import astropy.units as u
-    import matplotlib.pyplot as pl
+    import matplotlib.pyplot as plt
     import numpy as np
     import gala.potential as gp
     import gala.dynamics as gd
     from gala.units import galactic
 
     pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, q1=1., q2=1., q3=0.9, r_h=0,
                                   units=galactic)
-    w0 = gd.CartesianPhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
-                                        vel=[75, 150, 50.]*u.km/u.s)
+    w0 = gd.PhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
+                               vel=[75, 150, 50.]*u.km/u.s)
 
     w = pot.integrate_orbit(w0, dt=0.5, n_steps=10000)
     toy_potential = gd.fit_isochrone(w)
     actions,angles,freqs = toy_potential.action_angle(w)
-    fig,ax = pl.subplots(1,1,figsize=(5,5))
+    fig,ax = plt.subplots(1,1,figsize=(5,5))
     ax.plot(angles[0], angles[2], linestyle='none', marker=',')
     ax.set_xlim(0,2*np.pi)
     ax.set_ylim(0,2*np.pi)
@@ -169,7 +164,7 @@ Instead, the orbit is wobbly in the toy potential angles::
 This can also be seen in the value of the action variables, which are not
 time-independent in the toy potential::
 
-    >>> fig,ax = pl.subplots(1,1)
+    >>> fig,ax = plt.subplots(1,1)
     >>> ax.plot(w.t, toy_actions[0], marker=None) # doctest: +SKIP
     >>> ax.set_xlabel(r"$t$ [Myr]") # doctest: +SKIP
     >>> ax.set_ylabel(r"$J_1$ [rad]") # doctest: +SKIP
@@ -179,21 +174,21 @@ time-independent in the toy potential::
 
     import astropy.coordinates as coord
     import astropy.units as u
-    import matplotlib.pyplot as pl
+    import matplotlib.pyplot as plt
     import numpy as np
     import gala.potential as gp
     import gala.dynamics as gd
     from gala.units import galactic
 
     pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, q1=1., q2=1., q3=0.9, r_h=0,
                                   units=galactic)
-    w0 = gd.CartesianPhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
-                                        vel=[75, 150, 50.]*u.km/u.s)
+    w0 = gd.PhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
+                               vel=[75, 150, 50.]*u.km/u.s)
 
     w = pot.integrate_orbit(w0, dt=0.5, n_steps=10000)
     toy_potential = gd.fit_isochrone(w)
     actions,angles,freqs = toy_potential.action_angle(w)
-    fig,ax = pl.subplots(1,1)
+    fig,ax = plt.subplots(1,1)
     ax.plot(w.t, actions[0].to(u.km/u.s*u.kpc*u.Msun), marker=None)
     ax.set_xlabel(r"$t$ [Myr]")
     ax.set_ylabel(r"$J_1$ [kpc ${\rm M}_\odot$ km/s]")
@@ -222,7 +217,7 @@ actions computed using this machinery::
     >>> act_correction = nvecs.T[...,None] * result['Sn'][None,:,None] * np.cos(nvecs.dot(toy_angles))[None]
     >>> action_approx = toy_actions - 2*np.sum(act_correction, axis=1)*u.kpc**2/u.Myr*u.Msun
     >>>
-    >>> fig,ax = pl.subplots(1,1)
+    >>> fig,ax = plt.subplots(1,1)
     >>> ax.plot(w.t, toy_actions[0].to(u.km/u.s*u.kpc*u.Msun), marker=None, label='$J_1$') # doctest: +SKIP
     >>> ax.plot(w.t, action_approx[0].to(u.km/u.s*u.kpc*u.Msun), marker=None, label="$J_1'$") # doctest: +SKIP
     >>> ax.set_xlabel(r"$t$ [Myr]") # doctest: +SKIP
@@ -234,16 +229,16 @@ actions computed using this machinery::
 
     import astropy.coordinates as coord
     import astropy.units as u
-    import matplotlib.pyplot as pl
+    import matplotlib.pyplot as plt
     import numpy as np
     import gala.potential as gp
     import gala.dynamics as gd
     from gala.units import galactic
 
     pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, q1=1., q2=1., q3=0.9, r_h=0,
                                   units=galactic)
-    w0 = gd.CartesianPhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
-                                        vel=[75, 150, 50.]*u.km/u.s)
+    w0 = gd.PhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
+                               vel=[75, 150, 50.]*u.km/u.s)
 
     w = pot.integrate_orbit(w0, dt=0.5, n_steps=10000)
     toy_potential = gd.fit_isochrone(w)
@@ -252,7 +247,7 @@ actions computed using this machinery::
     nvecs = gd.generate_n_vectors(8, dx=1, dy=2, dz=2)
     act_correction = nvecs.T[...,None] * result['Sn'][None,:,None] * np.cos(nvecs.dot(toy_angles))[None]
     action_approx = toy_actions - 2*np.sum(act_correction, axis=1)*u.kpc**2/u.Myr*u.Msun
-    fig,ax = pl.subplots(1,1)
+    fig,ax = plt.subplots(1,1)
     ax.plot(w.t, toy_actions[0].to(u.km/u.s*u.kpc*u.Msun), marker=None, label='$J_1$')
     ax.plot(w.t, action_approx[0].to(u.km/u.s*u.kpc*u.Msun), marker=None, label="$J_1'$")
     ax.set_xlabel(r"$t$ [Myr]")
@@ -298,7 +293,7 @@ and the same initial conditions as above:
 
     import astropy.coordinates as coord
     import astropy.units as u
-    import matplotlib.pyplot as pl
+    import matplotlib.pyplot as plt
     import numpy as np
     import gala.potential as gp
     import gala.dynamics as gd
@@ -309,8 +304,8 @@ and the same initial conditions as above:
                                   units=galactic)
 
     # define initial conditions
-    w0 = gd.CartesianPhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
-                                        vel=[75, 150, 50.]*u.km/u.s)
+    w0 = gd.PhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
+                               vel=[75, 150, 50.]*u.km/u.s)
 
     # integrate orbit
     w = pot.integrate_orbit(w0, dt=0.5, n_steps=10000)
@@ -330,7 +325,7 @@ and the same initial conditions as above:
     act_correction = nvecs.T[...,None] * result['Sn'][None,:,None] * np.cos(nvecs.dot(toy_angles))[None]
     action_approx = toy_actions - 2*np.sum(act_correction, axis=1)*u.kpc**2/u.Myr*u.Msun
 
-    fig,axes = pl.subplots(3,1,figsize=(6,14))
+    fig,axes = plt.subplots(3,1,figsize=(6,14))
 
     for i,ax in enumerate(axes):
         ax.plot(w.t, toy_actions[i].to(u.km/u.s*u.kpc*u.Msun), marker=None, label='$J_{}$'.format(i+1))