tgas-to-galactocentric.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from astropy.io import fits\n",
    "from astropy import coordinates as coords\n",
    "from astropy import units as u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Converting to Galactocentric coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "- Bovy 2011 http://search.proquest.com/docview/883388353\n",
    "- Poleski 2013 https://arxiv.org/abs/1306.2945\n",
    "- astropy galactocentric coordinates http://docs.astropy.org/en/stable/coordinates/galactocentric.html#coordinates-galactocentric\n",
    "\n",
    "Uses [astropy](astropy.org) coordinates module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tgas = fits.getdata('/Users/semyeong/data/gaia/tgas_source/stacked_tgas.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_distance(parallax, parallax_error):\n",
    "    \"\"\"\n",
    "    Return the distance [kpc] point estimate with the Lutz-Kelker correction\n",
    "    \n",
    "    parallax : float, in mas\n",
    "    parallax_error : float, in mas\n",
    "    \"\"\"\n",
    "    snr = parallax / parallax_error\n",
    "    pnew = parallax * (0.5 + 0.5*np.sqrt(1 - 16./snr**2))\n",
    "    return 1./pnew"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ra dec parallax pmra pmdec = 153.808788479 -65.9406859114 2.03657951417 -22.6175107522 7.4909732812\n",
      "parallax S/N = 7.9900621273\n",
      "distance = 0.526374345342\n"
     ]
    }
   ],
   "source": [
    "source_id = 5245499633699505280\n",
    "\n",
    "i = where(tgas['source_id'] == source_id)[0][0]\n",
    "ra = tgas['ra'][i]\n",
    "dec = tgas['dec'][i]\n",
    "parallax = tgas['parallax'][i]\n",
    "parallax_error = tgas['parallax_error'][i]\n",
    "distance = get_distance(parallax, parallax_error)\n",
    "pmra = tgas['pmra'][i]\n",
    "pmdec = tgas['pmdec'][i]\n",
    "print('ra dec parallax pmra pmdec =', ra, dec, parallax, pmra, pmdec)\n",
    "print('parallax S/N =', parallax/parallax_error)\n",
    "print('distance =', distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/semyeong/anaconda2/envs/py35/lib/python3.5/site-packages/astropy/coordinates/builtin_frames/galactocentric.py:176: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "  xyz = R.dot(xyz.reshape(xyz.shape[0], np.prod(xyz.shape[1:]))).reshape(orig_shape)\n",
      "/Users/semyeong/anaconda2/envs/py35/lib/python3.5/site-packages/astropy/coordinates/builtin_frames/galactocentric.py:184: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "  xyz = R.dot(xyz.reshape(xyz.shape[0], np.prod(xyz.shape[1:]))).reshape(orig_shape)\n"
     ]
    }
   ],
   "source": [
    "c = coords.SkyCoord(ra*u.deg, dec*u.deg, distance=distance*u.kpc, frame='icrs')\n",
    "# Galactic (l,b)\n",
    "cg = c.transform_to(coords.Galactic)\n",
    "# Galactocentric (x,y,z)\n",
    "cgc = c.transform_to(coords.Galactocentric(galcen_distance=8*u.kpc, z_sun=0.*u.kpc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "287d54m33.245s -7d44m07.3956s\n",
      "-7.839607924897606 kpc -0.49631106663705576 kpc -0.07084955053980296 kpc\n"
     ]
    }
   ],
   "source": [
    "print(cg.l, cg.b)\n",
    "print(cgc.x, cgc.y, cgc.z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-22.8871032655 -6.62170502536\n"
     ]
    }
   ],
   "source": [
    "# find ra, dec of NGP\n",
    "ngp = coords.SkyCoord(0*u.deg, 90*u.deg, frame='galactic')\n",
    "ngp_icrs = ngp.transform_to(coords.ICRS)\n",
    "ra_ngp, dec_ngp = ngp_icrs.ra.to(u.rad).value, ngp_icrs.dec.to(u.rad).value\n",
    "\n",
    "c1 = sin(dec_ngp)*cos(deg2rad(dec)) - cos(dec_ngp)*sin(deg2rad(dec))*cos(deg2rad(ra)-ra_ngp)\n",
    "c2 = sin(deg2rad(ra)-ra_ngp)*cos(dec_ngp)\n",
    "norm = sqrt(c1**2 + c2**2)\n",
    "c1 = c1 / norm\n",
    "c2 = c2 / norm\n",
    "A = array([[c1, c2], [-c2, c1]])\n",
    "\n",
    "pml, pmb = einsum('ij,j->i', A, [pmra,pmdec])\n",
    "print(pml, pmb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "R1 = array([[cos(cg.l), -sin(cg.l), 0],\n",
    "            [sin(cg.l),  cos(cg.l), 0],\n",
    "            [0, 0, 1]])\n",
    "R2 = array([[cos(cg.b), 0, -sin(cg.b)],\n",
    "            [0, 1, 0],\n",
    "            [sin(cg.b), 0, cos(cg.b)]])\n",
    "R = R1.dot(R2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 1 mas/yr (as/yr) at 1 kpc (pc) is 4.74 km/s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vr = 0.\n",
    "vl = distance*pml*4.74\n",
    "vb = distance*pmb*4.74"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-55.02050572, -15.4439638 , -16.37091024])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "einsum('ij,j->i', R, [vr, vl, vb])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using galpy\n",
    "If you have [galpy](galpy.readthedocs.io) installed, it has the following routines already implemented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from galpy.util import bovy_coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-22.887103265465381, -6.6217050253646068)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# epoch=None for ICRS\n",
    "bovy_coords.pmrapmdec_to_pmllpmbb(pmra, pmdec, ra, dec, degree=True, epoch=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-55.025961339867202, -15.445495166367138, -16.372533516785733)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bovy_coords.vrpmllpmbb_to_vxvyvz(0, pml, pmb, cg.l.to(u.rad).value, cg.b.to(u.rad).value, distance,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}

You could add a block showing how to do it with Gala :)