VLSR_GSR.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord, EarthLocation, UnitSphericalRepresentation, FK5\n",
    "from astropy.table import Table\n",
    "from astropy.io import fits\n",
    "\n",
    "%matplotlib inline\n",
    "from matplotlib import rcParams\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see NED's velocities for M31 at\n",
    "\n",
    "http://ned.ipac.caltech.edu/cgi-bin/objsearch?objname=m31&extend=no&hconst=73&omegam=0.27&omegav=0.73&corr_z=1&out_csys=Equatorial&out_equinox=J2000.0&obj_sort=RA+or+Longitude&of=pre_text&zv_breaker=30000.0&list_limit=5&img_stamp=YES#DerivedValues_0\n",
    "\n",
    "The starting point for M31 is \n",
    "\n",
    "``V (Heliocentric)           :    -300 +/-      4 km/s``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m31 = SkyCoord.from_name('M31')\n",
    "m31_vhelio = -300*u.km/u.s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@u.quantity_input(vhelio=u.km/u.s, uvwlsr=u.km/u.s)\n",
    "def vhelio_to_vlsr(coord, vhelio=0*u.km/u.s, lsr_definition='dynamical'):\n",
    "    \"\"\"\n",
    "    Convert heliocentric radial velocity to Local Standard of Rest radial velocity\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    coord : SkyCoord\n",
    "        The direction at which to compute the conversion factor.\n",
    "    vhelio : Quantity with velocity units\n",
    "        The heliocentric radial velocity.  Should be a scalar quantity or match the shape of \n",
    "        `coord`.\n",
    "    lsr_definition : str or Quantity with velocity units\n",
    "        The definition of LSR to assume.  This can be one of three options:\n",
    "        * 'kinematic' : 20 km/s towards 18h,+30 deg (at 1900)\n",
    "        * 'dynamical' : IAU definition of (9, 12, 7) km/sec in Galactic cartesian coordinates\n",
    "        * A length-3 vector with the U, V, and W components (velocity along the galactic \n",
    "          x, y, and z axes) to assume for the LSR. \n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    vlsr : Quantity\n",
    "        The velocity in the Local Standard of Rest frame.\n",
    "    \"\"\"\n",
    "    if lsr_definition == 'kinematic':\n",
    "        direction = SkyCoord('18h', '30d', frame='fk5', equinox='J1900')\n",
    "        velocity = 20*u.km/u.s\n",
    "        uvw_lsr = direction.galactic.cartesian.xyz * velocity\n",
    "    elif lsr_definition == 'dynamical':\n",
    "        uvw_lsr = (9, 12, 7)*u.km/u.s\n",
    "    else:\n",
    "        uvw_lsr = lsr_definition\n",
    "    \n",
    "    # the unitspherical conversion ensures that the resulting cartesian is a *unit* vector\n",
    "    usr = coord.galactic.represent_as(UnitSphericalRepresentation)\n",
    "    cart = usr.to_cartesian()\n",
    "    \n",
    "    vlsr_vector = (cart.xyz.T*uvw_lsr).T\n",
    "    \n",
    "    return vhelio + np.sum(vlsr_vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-297.35836 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity -297.35835548549045 km / s>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m31_vlsr = vhelio_to_vlsr(m31, m31_vhelio)\n",
    "m31_vlsr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For M31, NED sez:\n",
    "\n",
    "``V (Kinematic LSR)          :    -296 +/-      4 km/s ``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-295.60223 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity -295.6022272005147 km / s>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# this is the \"kinematic\" LSR that NED assumes\n",
    "vhelio_to_vlsr(m31, m31_vhelio, 'kinematic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@u.quantity_input(vlsr=u.km/u.s, vrot=u.km/u.s)\n",
    "def vlsr_to_vgsr(coord, vlsr=0*u.km/u.s, vrot=220*u.km/u.s):\n",
    "    \"\"\"\n",
    "    Convert Local Standard of Rest radial velocity to galactocentric velocity\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    coord : SkyCoord\n",
    "        The direction at which to compute the conversion factor.\n",
    "    vlsr : Quantity with velocity units\n",
    "        The heliocentric radial velocity.  Should be a scalar quantity or match the shape of \n",
    "        `coord`.\n",
    "    vrot : Quantity with velocity units\n",
    "        The value of the galactic rotation curve at the solar cycle. The default is the \n",
    "        IAU convention of 220 km/s.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    vlsr : Quantity\n",
    "        The velocity in the Local Standard of Rest frame.\n",
    "    \"\"\"\n",
    "    # this adds the projection along the l=90, b=0 axis.  Could do the whole vector approach as above\n",
    "    # but along the axis makes it trivial to do this simpler version\n",
    "    l = coord.galactic.l\n",
    "    b = coord.galactic.b\n",
    "    return vlsr + vrot * np.sin(l) * np.cos(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For M31, NED says:\n",
    "\n",
    "``V (Galactocentric GSR)     :    -122 +/-      8 km/s``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-122.31317 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity -122.31316955594377 km / s>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vlsr_to_vgsr(m31, m31_vlsr)"
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}