{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Transforming positions and velocities to and from a Galactocentric frame\n\n\nThis document shows a few examples of how to use and customize the\n`~astropy.coordinates.Galactocentric` frame to transform Heliocentric sky\npositions, distance, proper motions, and radial velocities to a Galactocentric,\nCartesian frame, and the same in reverse.\n\nThe main configurable parameters of the `~astropy.coordinates.Galactocentric`\nframe control the position and velocity of the solar system barycenter within\nthe Galaxy. These are specified by setting the ICRS coordinates of the\nGalactic center, the distance to the Galactic center (the sun-galactic center\nline is always assumed to be the x-axis of the Galactocentric frame), and the\nCartesian 3-velocity of the sun in the Galactocentric frame. We'll first\ndemonstrate how to customize these values, then show how to set the solar motion\ninstead by inputting the proper motion of Sgr A*.\n\nNote that, for brevity, we may refer to the solar system barycenter as just \"the\nsun\" in the examples below.\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, set up numpy,\nmatplotlib, and use a nicer set of plot parameters:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom astropy.visualization import astropy_mpl_style\nplt.style.use(astropy_mpl_style)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Import the necessary astropy subpackages\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import astropy.coordinates as coord\nimport astropy.units as u"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's first define a barycentric coordinate and velocity in the ICRS frame.\nWe'll use the data for the star HD 39881 from the `Simbad\n<simbad.harvard.edu/simbad/>`_ database:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "c1 = coord.ICRS(ra=89.014303*u.degree, dec=13.924912*u.degree,\n                distance=(37.59*u.mas).to(u.pc, u.parallax()),\n                pm_ra_cosdec=372.72*u.mas/u.yr,\n                pm_dec=-483.69*u.mas/u.yr,\n                radial_velocity=0.37*u.km/u.s)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is a high proper-motion star; suppose we'd like to transform its position\nand velocity to a Galactocentric frame to see if it has a large 3D velocity\nas well. To use the Astropy default solar position and motion parameters, we\ncan simply do:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gc1 = c1.transform_to(coord.Galactocentric)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "From here, we can access the components of the resulting\n`~astropy.coordinates.Galactocentric` instance to see the 3D Cartesian\nvelocity components:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(gc1.v_x, gc1.v_y, gc1.v_z)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The default parameters for the `~astropy.coordinates.Galactocentric` frame\nare detailed in the linked documentation, but we can modify the most commonly\nchanges values using the keywords ``galcen_distance``, ``galcen_v_sun``, and\n``z_sun`` which set the sun-Galactic center distance, the 3D velocity vector\nof the sun, and the height of the sun above the Galactic midplane,\nrespectively. The velocity of the sun must be specified as a\n`~astropy.coordinates.CartesianDifferential` instance, as in the example\nbelow. Note that, as with the positions, the Galactocentric frame is a\nright-handed system - the x-axis is positive towards the Galactic center, so\n``v_x`` is opposite of the Galactocentric radial velocity:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "v_sun = coord.CartesianDifferential([11.1, 244, 7.25]*u.km/u.s)\ngc_frame = coord.Galactocentric(galcen_distance=8*u.kpc,\n                                galcen_v_sun=v_sun,\n                                z_sun=0*u.pc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can then transform to this frame instead, with our custom parameters:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gc2 = c1.transform_to(gc_frame)\nprint(gc2.v_x, gc2.v_y, gc2.v_z)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It's sometimes useful to specify the solar motion using the `proper motion\nof Sgr A* <https://arxiv.org/abs/astro-ph/0408107>`_ instead of Cartesian\nvelocity components. With an assumed distance, we can convert proper motion\ncomponents to Cartesian velocity components using `astropy.units`:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "galcen_distance = 8*u.kpc\npm_gal_sgrA = [-6.379, -0.202] * u.mas/u.yr # from Reid & Brunthaler 2004\nvy, vz = -(galcen_distance * pm_gal_sgrA).to(u.km/u.s, u.dimensionless_angles())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We still have to assume a line-of-sight velocity for the Galactic center,\nwhich we will again take to be 11 km/s:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "vx = 11.1 * u.km/u.s\n\ngc_frame2 = coord.Galactocentric(galcen_distance=galcen_distance,\n                                 galcen_v_sun=coord.CartesianDifferential(vx, vy, vz),\n                                 z_sun=0*u.pc)\ngc3 = c1.transform_to(gc_frame2)\nprint(gc3.v_x, gc3.v_y, gc3.v_z)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The transformations also work in the opposite direction. This can be useful\nfor transforming simulated or theoretical data to observable quantities. As\nan example, we'll generate 4 theoretical circular orbits at different\nGalactocentric radii with the same circular velocity, and transform them to\nHeliocentric coordinates:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ring_distances = np.arange(10, 25+1, 5) * u.kpc\ncirc_velocity = 220 * u.km/u.s\n\nphi_grid = np.linspace(90, 270, 512) * u.degree # grid of azimuths\nring_rep = coord.CylindricalRepresentation(\n    rho=ring_distances[:,np.newaxis],\n    phi=phi_grid[np.newaxis],\n    z=np.zeros_like(ring_distances)[:,np.newaxis])\n\nangular_velocity = (-circ_velocity / ring_distances).to(u.mas/u.yr,\n                                                        u.dimensionless_angles())\nring_dif = coord.CylindricalDifferential(\n    d_rho=np.zeros(phi_grid.shape)[np.newaxis]*u.km/u.s,\n    d_phi=angular_velocity[:,np.newaxis],\n    d_z=np.zeros(phi_grid.shape)[np.newaxis]*u.km/u.s\n)\n\nring_rep = ring_rep.with_differentials(ring_dif)\ngc_rings = coord.Galactocentric(ring_rep)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, let's visualize the geometry in Galactocentric coordinates. Here are\nthe positions and velocities of the rings; note that in the velocity plot,\nthe velocities of the 4 rings are identical and thus overlaid under the same\ncurve:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig,axes = plt.subplots(1, 2, figsize=(12,6))\n\n# Positions\naxes[0].plot(gc_rings.x.T, gc_rings.y.T, marker='None', linewidth=3)\naxes[0].text(-8., 0, r'$\\odot$', fontsize=20)\n\naxes[0].set_xlim(-30, 30)\naxes[0].set_ylim(-30, 30)\n\naxes[0].set_xlabel('$x$ [kpc]')\naxes[0].set_ylabel('$y$ [kpc]')\n\n# Velocities\naxes[1].plot(gc_rings.v_x.T, gc_rings.v_y.T, marker='None', linewidth=3)\n\naxes[1].set_xlim(-250, 250)\naxes[1].set_ylim(-250, 250)\n\naxes[1].set_xlabel('$v_x$ [{0}]'.format((u.km/u.s).to_string(\"latex_inline\")))\naxes[1].set_ylabel('$v_y$ [{0}]'.format((u.km/u.s).to_string(\"latex_inline\")))\n\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we can transform to Galactic coordinates and visualize the rings in\nobservable coordinates:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "gal_rings = gc_rings.transform_to(coord.Galactic)\n\nfig,ax = plt.subplots(1, 1, figsize=(8,6))\nfor i in range(len(ring_distances)):\n    ax.plot(gal_rings[i].l.degree, gal_rings[i].pm_l_cosb.value,\n            label=str(ring_distances[i]), marker='None', linewidth=3)\n\nax.set_xlim(360, 0)\n\nax.set_xlabel('$l$ [deg]')\nax.set_ylabel(r'$\\mu_l \\, \\cos b$ [{0}]'.format((u.mas/u.yr).to_string('latex_inline')))\n\nax.legend()"
      ]
    }
  ],
  "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
}