{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n================================================================\nConvert a radial velocity to the Galactic Standard of Rest (GSR)\n================================================================\n\nRadial or line-of-sight velocities of sources are often reported in a\nHeliocentric or Solar-system barycentric reference frame. A common\ntransformation incorporates the projection of the Sun's motion along the\nline-of-sight to the target, hence transforming it to a Galactic rest frame\ninstead (sometimes referred to as the Galactic Standard of Rest, GSR). This\ntransformation depends on the assumptions about the orientation of the Galactic\nframe relative to the bary- or Heliocentric frame. It also depends on the\nassumed solar velocity vector. Here we'll demonstrate how to perform this\ntransformation using a sky position and barycentric radial-velocity.\n\n-------------------\n\n*By: Adrian Price-Whelan*\n\n*License: BSD*\n\n-------------------\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make print work the same in all versions of Python and import the required\nAstropy packages:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import astropy.units as u\nimport astropy.coordinates as coord"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For this example, let's work with the coordinates and barycentric radial\nvelocity of the star HD 155967, as obtained from\n`Simbad <http://simbad.harvard.edu/simbad/>`_:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "icrs = coord.ICRS(ra=258.58356362*u.deg, dec=14.55255619*u.deg,\n                  radial_velocity=-16.1*u.km/u.s)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We next need to decide on the velocity of the Sun in the assumed GSR frame.\nWe'll use the same velocity vector as used in the\n`~astropy.coordinates.Galactocentric` frame, and convert it to a\n`~astropy.coordinates.CartesianRepresentation` object using the\n``.to_cartesian()`` method of the\n`~astropy.coordinates.CartesianDifferential` object ``galcen_v_sun``:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "v_sun = coord.Galactocentric.galcen_v_sun.to_cartesian()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now need to get a unit vector in the assumed Galactic frame from the sky\nposition in the ICRS frame above. We'll use this unit vector to project the\nsolar velocity onto the line-of-sight:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gal = icrs.transform_to(coord.Galactic)\ncart_data = gal.data.to_cartesian()\nunit_vector = cart_data / cart_data.norm()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we project the solar velocity using this unit vector:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "v_proj = v_sun.dot(unit_vector)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we add the projection of the solar velocity to the radial velocity\nto get a GSR radial velocity:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rv_gsr = icrs.radial_velocity + v_proj\nprint(rv_gsr)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We could wrap this in a function so we can control the solar velocity and\nre-use the above code:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def rv_to_gsr(c, v_sun=None):\n    \"\"\"Transform a barycentric radial velocity to the Galactic Standard of Rest\n    (GSR).\n\n    The input radial velocity must be passed in as a\n\n    Parameters\n    ----------\n    c : `~astropy.coordinates.BaseCoordinateFrame` subclass instance\n        The radial velocity, associated with a sky coordinates, to be\n        transformed.\n    v_sun : `~astropy.units.Quantity` (optional)\n        The 3D velocity of the solar system barycenter in the GSR frame.\n        Defaults to the same solar motion as in the\n        `~astropy.coordinates.Galactocentric` frame.\n\n    Returns\n    -------\n    v_gsr : `~astropy.units.Quantity`\n        The input radial velocity transformed to a GSR frame.\n\n    \"\"\"\n    if v_sun is None:\n        v_sun = coord.Galactocentric.galcen_v_sun.to_cartesian()\n\n    gal = c.transform_to(coord.Galactic)\n    cart_data = gal.data.to_cartesian()\n    unit_vector = cart_data / cart_data.norm()\n\n    v_proj = v_sun.dot(unit_vector)\n\n    return c.radial_velocity + v_proj\n\nrv_gsr = rv_to_gsr(icrs)\nprint(rv_gsr)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}