Spectralcoord_target_projection_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy import units as u\n",
    "from astropy import time\n",
    "from astropy.coordinates import SkyCoord, EarthLocation, ICRS, GCRS, Galactic, CartesianDifferential, SphericalRepresentation, RadialDifferential, get_body_barycentric_posvel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from astropy.visualization import quantity_support\n",
    "quantity_support()\n",
    "\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "keck = EarthLocation.of_site('keck')\n",
    "obstime = time.Time('2018-12-13 9:00')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<GCRS Coordinate (obstime=2018-12-13 09:00:00.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, m)\n",
       "    (61.24879773, 19.65742965, 6379854.71495458)\n",
       " (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)\n",
       "    (4.46861165e+11, 7.59822165e+08, 2.14743051e-08)>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "observer_gcrs = keck.get_gcrs(obstime)\n",
    "observer_gcrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<CartesianDifferential (d_x, d_y, d_z) in km / s\n",
       "    (-0.38410488, 0.21045191, 0.00070133)>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "observer_gcrs.velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Jupiter (in the solar system):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos, vel = get_body_barycentric_posvel('jupiter', obstime)\n",
    "jupiter_icrs = ICRS(pos.with_differentials(CartesianDifferential(vel.xyz)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<GCRS Coordinate (obstime=2018-12-13 09:00:00.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, AU)\n",
       "    (245.80154929, -20.92556215, 6.30709173)\n",
       " (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)\n",
       "    (2.86074921e+08, -48407081.05369898, -7.22056426)>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jupiter_gcrs = jupiter_icrs.transform_to(observer_gcrs)\n",
    "jupiter_gcrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<CartesianDifferential (d_x, d_y, d_z) in km / s\n",
       "    (41.61529311, -8.55978957, -3.9750781)>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jupiter_gcrs.velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the LOS velocity from observer to jupiter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<CartesianRepresentation (x, y, z) [dimensionless]\n",
       "    (-0.3828641, -0.85197163, -0.35715462)>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_pos = jupiter_gcrs.data.without_differentials() - observer_gcrs.data.without_differentials()\n",
    "pos_hat = d_pos/d_pos.norm()\n",
    "pos_hat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-7.1880966 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity -7.18809664 km / s>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_vel = jupiter_gcrs.velocity - observer_gcrs.velocity\n",
    "\n",
    "rv = np.sum(d_vel.d_xyz* pos_hat.xyz, axis=0)\n",
    "rv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-7.2205643 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity -7.22056426 km / s>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jupiter_gcrs.radial_velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not the same!  Which makes sense on reflection - there's a difference between the radial velocity from the origin of the *frame* and the radial velocity from the *observer*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alpha cen "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<SkyCoord (ICRS): (ra, dec) in deg\n",
       "    (219.90085, -60.83561944)>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SkyCoord.from_name('alpha cen')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rough estimates\n",
    "αcen = ICRS(ra=219.90085*u.deg, dec=-60.83562*u.deg, \n",
    "              distance=4.37*u.lightyear, radial_velocity=-18.*u.km/u.s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<GCRS Coordinate (obstime=2018-12-13 09:00:00.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, lyr)\n",
       "    (219.89232033, -60.83197963, 4.37001073)\n",
       " (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)\n",
       "    (83197.10620566, 52735.6421058, -26.32813814)>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "αcen_gcrs = αcen.transform_to(observer_gcrs)\n",
    "αcen_gcrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<CartesianDifferential (d_x, d_y, d_z) in km / s\n",
       "    (124.34698135, -584.78772197, 186.23582661)>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "αcen_gcrs.velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$384.13206 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity 384.1320584 km / s>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_vel = αcen_gcrs.velocity - observer_gcrs.velocity\n",
    "rv = np.sum(d_vel.d_xyz* pos_hat.xyz, axis=0)\n",
    "rv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Those are very wrong.  Probably due to the transformation and its numerical precision issues"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Where does it go wrong?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$276363.51 \\; \\mathrm{AU}$"
      ],
      "text/plain": [
       "<Distance 276363.50685824 AU>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "αcen.distance.to(u.au)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$[10,~11.23324,~12.618569,~14.174742,~15.922828,~17.886495,~20.09233,~22.570197,~25.353645,~28.480359,~31.992671,~35.938137,~40.370173,~45.348785,~50.94138,~57.223677,~64.280731,~72.20809,~81.113083,~91.116276,~102.3531,~114.9757,~129.15497,~145.08288,~162.97508,~183.07383,~205.65123,~231.01297,~259.50242,~291.50531,~327.45492,~367.83798,~413.20124,~464.15888,~521.40083,~585.70208,~657.93322,~739.0722,~830.21757,~932.60335,~1047.6158,~1176.812,~1321.9411,~1484.9683,~1668.1005,~1873.8174,~2104.9041,~2364.4894,~2656.0878,~2983.6472,~3351.6027,~3764.9358,~4229.2429,~4750.8102,~5336.6992,~5994.8425,~6734.1507,~7564.6333,~8497.5344,~9545.4846,~10722.672,~12045.035,~13530.478,~15199.111,~17073.526,~19179.103,~21544.347,~24201.283,~27185.882,~30538.555,~34304.693,~38535.286,~43287.613,~48626.016,~54622.772,~61359.073,~68926.121,~77426.368,~86974.9,~97700.996,~109749.88,~123284.67,~138488.64,~155567.61,~174752.84,~196304.07,~220513.07,~247707.64,~278255.94,~312571.58,~351119.17,~394420.61,~443062.15,~497702.36,~559081.02,~628029.14,~705480.23,~792482.9,~890215.09,~1000000] \\; \\mathrm{AU}$"
      ],
      "text/plain": [
       "<Quantity [1.00000000e+01, 1.12332403e+01, 1.26185688e+01, 1.41747416e+01,\n",
       "           1.59228279e+01, 1.78864953e+01, 2.00923300e+01, 2.25701972e+01,\n",
       "           2.53536449e+01, 2.84803587e+01, 3.19926714e+01, 3.59381366e+01,\n",
       "           4.03701726e+01, 4.53487851e+01, 5.09413801e+01, 5.72236766e+01,\n",
       "           6.42807312e+01, 7.22080902e+01, 8.11130831e+01, 9.11162756e+01,\n",
       "           1.02353102e+02, 1.14975700e+02, 1.29154967e+02, 1.45082878e+02,\n",
       "           1.62975083e+02, 1.83073828e+02, 2.05651231e+02, 2.31012970e+02,\n",
       "           2.59502421e+02, 2.91505306e+02, 3.27454916e+02, 3.67837977e+02,\n",
       "           4.13201240e+02, 4.64158883e+02, 5.21400829e+02, 5.85702082e+02,\n",
       "           6.57933225e+02, 7.39072203e+02, 8.30217568e+02, 9.32603347e+02,\n",
       "           1.04761575e+03, 1.17681195e+03, 1.32194115e+03, 1.48496826e+03,\n",
       "           1.66810054e+03, 1.87381742e+03, 2.10490414e+03, 2.36448941e+03,\n",
       "           2.65608778e+03, 2.98364724e+03, 3.35160265e+03, 3.76493581e+03,\n",
       "           4.22924287e+03, 4.75081016e+03, 5.33669923e+03, 5.99484250e+03,\n",
       "           6.73415066e+03, 7.56463328e+03, 8.49753436e+03, 9.54548457e+03,\n",
       "           1.07226722e+04, 1.20450354e+04, 1.35304777e+04, 1.51991108e+04,\n",
       "           1.70735265e+04, 1.91791026e+04, 2.15443469e+04, 2.42012826e+04,\n",
       "           2.71858824e+04, 3.05385551e+04, 3.43046929e+04, 3.85352859e+04,\n",
       "           4.32876128e+04, 4.86260158e+04, 5.46227722e+04, 6.13590727e+04,\n",
       "           6.89261210e+04, 7.74263683e+04, 8.69749003e+04, 9.77009957e+04,\n",
       "           1.09749877e+05, 1.23284674e+05, 1.38488637e+05, 1.55567614e+05,\n",
       "           1.74752840e+05, 1.96304065e+05, 2.20513074e+05, 2.47707636e+05,\n",
       "           2.78255940e+05, 3.12571585e+05, 3.51119173e+05, 3.94420606e+05,\n",
       "           4.43062146e+05, 4.97702356e+05, 5.59081018e+05, 6.28029144e+05,\n",
       "           7.05480231e+05, 7.92482898e+05, 8.90215085e+05, 1.00000000e+06] AU>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances_to_try = np.logspace(1, 6, 100)*u.au\n",
    "distances_to_try"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x113301710>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEQCAYAAACa+vIpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl81dWd//HXJxv7HjYhiAji1go8QFymwlRU1CIdtXWpLdQF22Jbfr9qi47LaO2UPtT5qeMygkWqWKGDtMAUscAUbIeqLGUwKCgEAgECWTCBJGS75/fHvcEQsn2Te7/fe5P38/HII/eee+45n0MgH8453/s95pxDRESkpZKCDkBERBKbEomIiLSKEomIiLSKEomIiLSKEomIiLSKEomIiLSKEomIiLSKEomIiLSKEomIiLSKEomIiLRKStAB+CE9Pd0NHTo06DBEpI2oqKgAIC0tLeBIYmvz5s35zrm+TdVrF4lk6NChbNq0KegwRKSNmDhxIgDr1q0LNI5YM7Ps5tRrF4lERCSaHn744aBDiCtKJCIiHk2aNCnoEOKKNttFRDzKysoiKysr6DDihmYkIiIe3XnnnUDb3yNpLiUSERGPHn/88aBDiCtKJCIiHk2YMCHoEOKK9khERDzauXMnO3fuDDqMJv33jsN8dvhYzPtRIhER8ejee+/l3nvvDTqMRjnn+PGirbz+t2Z9FKRVtLQlIuLRv/7rvwYdQpOOllZy7EQVZ/bpHPO+lEhERDy67LLLgg6hSdkFJQAM7dMl5n1paUtExKPMzEwyMzODDqNR2QWlAAxN14xERCTu3HfffUB8f45kb0EJZjC4lxKJiEjceeqpp4IOoUnZBaUM7N6RjqnJMe9LiURExKNx48YFHUKTsgtKONOH/RHQHomIiGdbt25l69atQYfRqOyCUl/2R0AzEhERz2bNmgXE7x5J8YlKCkoqGNLbnxmJEomIiEfPPvts0CE0al/NFVs+fIYElEhERDwbNWpU0CE0qubSX+2RiIjEqY0bN7Jx48agw2jQ3siHEf34VDtoRiIi4tkDDzwAxO8eSXZBCeldO9Clgz+/4pVIREQ8euGFF4IOoVHZBaW+7Y+AEomIiGcXXnhh0CE0KruglMuHp/vWn/ZIREQ82rBhAxs2bAg6jHqVVVSTW3xCMxIRkXj20EMPAfG5R7KvMHzF1hAlEhGR+PXKK68EHUKD/Lx9fA3flrbMLMPM/mxmn5jZdjP7caS8t5mtNrPPIt97RcrNzJ43s11mts3MxtRqa1qk/mdmNs2vMYiIAIwcOZKRI0cGHUa9Tt4+vi0mEqAK+Ilz7jzgEmCmmZ0PzAbWOudGAGsjzwGuBUZEvmYAL0M48QCPAeOBi4HHapKPiIgf1q9fz/r164MOo157C0ro2TmVHp1TfevTt6Ut59wh4FDk8TEz+wQYBEwFJkaq/QZYB/wsUv66c84B75tZTzMbGKm72jlXCGBmq4HJwFt+jUVE2rfHHnsMiM89kuyCUs7s7d/+CAS0R2JmQ4HRwAdA/0iSwTl3yMz6RaoNAvbXeltOpKyhchERX8yfPz/oEBqUXVjC6Ax/F2l8v/zXzLoCbwOznHPFjVWtp8w1Ul63nxlmtsnMNuXl5bUsWBGRegwbNoxhw4YFHcZpyquqOXC0zNdLf8HnRGJmqYSTyJvOuaWR4sORJSsi349EynOAjFpvHwwcbKT8FM65uc65sc65sX379o3uQESkXVuzZg1r1qwJOozTZBeUEnJwdr+uvvbr51VbBvwa+MQ592+1XloO1Fx5NQ1YVqv8O5Grty4BiiJLYO8CV5tZr8gm+9WRMhERXzz55JM8+eSTQYdxmt1HjgNwdl9/E4mfeySXA98GPjKzmqPFHgLmAL8zs7uAfcA3Iq+tBK4DdgGlwHcBnHOFZvZzoObWm0/UbLyLiPjhjTfeCDqEeu3OCyeSYX39u/QX/L1q66/Uv78BcGU99R0ws4G25gPxu9slIm1aRkZG05UCsOvIcQb17ETnNH+vo9K9tkREPFq1ahWrVq0KOozT7M4r8X02ArpFioiIZ3PmzAFg8uTJAUfyBeccu/OO882x/s+WlEhERDxatGhR0CGcJrf4BKUV1b5fsQVKJCIing0YMCDoEE6z6+QVW/4vbWmPRETEoxUrVrBixYqgwzhFzaW/wzUjERGJf8888wwAU6ZMCTiSL+zOK6FbxxT6du3ge99KJCIiHi1ZsiToEE6zO+84Z/ftSviz3/5SIhER8Sg93b/z0Jtrd95x/mF4MLeD0h6JiIhHS5cuZenSpU1X9EnxiUoOF5cHsj8CmpGIiHj2/PPPA3DjjTcGHElYVl74eN0grtgCJRIREc+WLVvWdCUfnbxZo2YkIiKJoUePHkGHcIrdecdJSTKG+HwyYg3tkYiIeLR48WIWL14cdBgn7TpynDP7dCY1OZhf6ZqRiIh49PLLLwNwyy23BBxJ2O6844FttIMSiYiIZytXrgw6hJMqq0NkF5RyzQXB3bZFiURExKPOnYPZi6jP3vwSqkIu0BmJ9khERDxauHAhCxcuDDoMAHbkHgNg5IBugcWgGYmIiEevvvoqAHfccUfAkcCO3GKSk0x7JCIiiWT16tVBh3DSztxjDEvvQoeU5MBi0NKWiIhHqamppKamBh0GEF7aCnJZC5RIREQ8W7BgAQsWLAg6DI6XV5FztIxzlUhERBJLvCSSnSc32rsHGof2SEREPFq3bl3QIQBfJBLNSEREpEV25hbTJS2ZQT07BRqHEomIiEfz5s1j3rx5QYfBjtxjnDOgG0lJ/p+KWJsSiYiIR/Fw00bnHDsPHwt8WQu0RyIi4tmaNWuCDoHDxeV8XlrJyP7BJxLNSEREEtCO3GIg+Cu2QIlERMSzl156iZdeeinQGOLlii1QIhER8WzFihWsWLEi0Bh25h6jf/cO9OqSFmgcoD0SERHP3nnnnaBDiNwaJfhlLdCMREQk4VRVh9iVdzwulrVAiURExLPnnnuO5557LrD+9+SXUFEViosrtkCJRETEs7Vr17J27drA+v/oQBEAFw7qEVgMtWmPRETEo+XLlwfaf+aBYjqmJnF23y6BxlFDMxIRkQSTeaCI8wd2JyU5Pn6F+xaFmc03syNmllmr7F/M7ICZbY18XVfrtQfNbJeZ7TSza2qVT46U7TKz2X7FLyJS4+mnn+bpp58OpO9QyLH9YFHcLGuBv0tbC4AXgNfrlP8/59wpPxEzOx+4FbgAOANYY2bnRF5+EbgKyAE2mtly59zHsQxcRKS2v/3tb4H1vaeghJKK6vaZSJxz75nZ0GZWnwoscs6VA3vMbBdwceS1Xc65LAAzWxSpq0QiIr55++23A+s7M7LR/qU4SiTxsMB2n5ltiyx99YqUDQL216qTEylrqFxEpF34KKeItJQkhvfrGnQoJwWdSF4GzgZGAYeAZyLl9d1c3zVSfhozm2Fmm8xsU15eXjRiFREBYM6cOcyZMyeQvjMPFnHewO6kxslGOwR8+a9z7nDNYzObB/xX5GkOkFGr6mDgYORxQ+V1254LzAUYO3ZsvclGRKQltm7dGki/oZBj+4Fipo4+I5D+GxJoIjGzgc65Q5Gn/wTUXNG1HPitmf0b4c32EcCHhGckI8zsLOAA4Q352/2NWkTau0WLFgXSb3ZhKcfKq7jwjPjZHwEfE4mZvQVMBNLNLAd4DJhoZqMIL0/tBe4FcM5tN7PfEd5ErwJmOueqI+3cB7wLJAPznXPb/RqDiEiQMuPsE+01/Lxq67Z6in/dSP1fAL+op3wlsDKKoYmIePLzn/8cgEceecTXfjMPFJGWnMQ5cXKPrRq6RYqIiEc7d+4MpN+PDhQxckA30lLiZ6MdlEhERDxbuHCh730658g8UMT1X46vjXYI/vJfERFphv2FZRSfqIqrDyLWUCIREfHo0Ucf5dFHH/W1z7/vPwrAlwfHXyLR0paIiEf79+9vulKU/X3f53ROS46bUxFrUyIREfHotdde873PLfuO8uXBPeLm1vG1xV9EIiJyihOV1Xx8sJgxQ3o1XTkASiQiIh49+OCDPPjgg771ty2niKqQY3ScJhItbYmIeFRQUOBrf1v2hTfaRw/p6Wu/zaVEIiLi0dy5c33tb0v2Uc7s05n0rh187be5tLQlIhLHnHNs2fd53O6PgBKJiIhn999/P/fff78vfeUcLSP/eDlj4nRZC1qQSMzsZ7EIREQkUZSVlVFWVuZLX1/sj8TvjKTJPZLI7dxPPiV8muGvYhaRiEice/HFF33ra0v20bj9IGKN5my2Fzvn7q55YmYvxzAeERGpZcu+z+P2g4g1mhNZ3TNB/jkWgYiIJIpZs2Yxa9asmPdTVlHNJ4fi94OINZpMJM65PQBmlh55XhjroEREBLblfE5VyMV9IvHyOZL5wA2xCkREJFE8++yzvvTz4Z5CzGDs0PhOJF4W3SxmUYiIyGne31PAuQO607NzWtChNMpLInExi0JEJIHMnDmTmTNnxrSP8qpqNmcf5ZJhvWPaTzR4WdrSjEREBOjUqVPM+9iWU8SJyhCXDOsT875ay0si8e9WlyIicezpp5+OeR/v7y7ADMafFf8zkmYvbTnnMs3sG2bWDcDMHjazpWY2JnbhiYi0T4myPwLeb5HyiHPumJn9A3AN8BtAH1AUkXZlxowZzJgxI2btJ9L+CHhPJNWR79cDLzvnlgHxny5FRKKoT58+9OkTu72Lmv2RSxNgfwS8n0dywMxeASYBvzKzDugOwiLSzvzyl7+Mafs1+yMXJ8D+CHhPAt8E3gUmO+c+B3oDD0Q9KhGRduz9PQWclyD7I+BxRuKcKwWW1np+CDgU7aBEROLZd7/7XQBee+21qLddsz9y+8VnRr3tWNFRuyIiHmVkZMSs7f/dX/P5kcRY1gIlEhERz5544omYtf3XXfkkGYw/KzE22sFjIjGzsYRvI39m5L0GOOfcl2MQm4hIu7P+0zxGZfSkR+fUoENpNq8zkjcJb65/BISiH46ISPy74447AFi4cGFU2y0sqWBbzufMuvKcqLYba14TSZ5zbnlMIhERSRAjR46MSbt/+SwP52DCyL4xaT9WvCaSx8zsVWAtUF5T6Jxb2vBbRETalkceeSQm7a7/NI9enVP50qAeMWk/Vrwmku8C5wKpfLG05ah1SbCIiHgXCjne+zSfr4zoS3JSYt1s3Wsiucg596WYRCIikiBuvfVWABYtWhS1Nj/JLSb/eDlXnJNYy1rg/ZPt75vZ+S3pyMzmm9kRM8usVdbbzFab2WeR770i5WZmz5vZLjPbVvsOw2Y2LVL/MzOb1pJYRERaY9SoUYwaNSqqba7/NA+AK0akR7VdP3idkfwDMM3M9hDeI/Fy+e8C4AXg9Vpls4G1zrk5ZjY78vxnwLXAiMjXeMJ3GB5vZr2Bx4CxhJfUNpvZcufcUY/jEBFpsdmzZ0e9zfU78zh/YHf6de8Y9bZjzWsi+TaQV6fshua80Tn3npkNrVM8FZgYefwbYB3hRDIVeN055wjPgnqa2cBI3dXOuUIAM1sNTAbe8jgOEZG4cexEJZuzj3LPFcOCDqVFvC5tzQW6O+eynXPZwKXAHa3ov3/kfl019+3qFykfBOyvVS8nUtZQuYiIb2666SZuuummqLW3YXcBVSHHhATcHwHvM5KbgSVm9i3Cy1zfAa6OelT1nw/vGik/vQGzGcAMgCFDhkQvMhFp9y699NKotrfm48N065jCmCG9otquX7ze/TfLzG4F/kB4ZnC1c66sFf0fNrOBzrlDkaWrI5HyHKD2XdEGAwcj5RPrlK9rINa5hGdQjB07tt5kIyLSEvfff3/U2qqqDrHmk8N89dx+pKUk5vFOzYrazD6KXD21DVhC+BySocAHkbKWWg7UXHk1DVhWq/w7kau3LgGKIktf7wJXm1mvyBVeV0fKREQS0qbsoxwtreSaCwYEHUqLNXdG8rXWdmRmbxGeTaSbWQ7hq6/mAL8zs7uAfcA3ItVXAtcBu4BSwh+ExDlXaGY/BzZG6j1Rs/EuIuKXG24IX2O0fHnr7xj17vZc0lKSEnZ/BJqZSCIb663inLutgZeurKeuA2Y20M58YH5r4xERaakrrzzt11aLOOf40/bDfGV4Ol06JO6pHokbuYhIQH784x9HpZ3tB4s58HkZP75yRFTaC0pi7uyIiLQBf9qeS5LBlef1a7pyHFMiERHx6Nprr+Xaa69tdTvvbj/MuKG96dO1QxSiCo6WtkREPJoyZUqr29ibX8LOw8d49Gstun1hXFEiERHx6Ac/+EGr23h3ey4AV53fv9VtBU1LWyIiAVi29SAXZfQko3fnoENpNSUSERGPJk2axKRJk1r8/s8OH+PjQ8V8fdQZUYwqOFraEhHx6JZbbmnV+/+w9QBJBtd/eWCUIgqWEomIiEf33HNPi9/rnGPZ1oNcPjydft0S7+yR+mhpS0TER1v2HSXnaBlfH9V2TsBQIhER8WjixIlMnDixRe/9w98P0jE1iWsuTNybNNalpS0REY+mT5/eovdVVof440eHmHRef7om8L216mo7IxER8UlLE8lfPsujsKSiTS1rgZa2REQ8q6yspLKy0vP73t5ygF6dU7kigW8ZXx8lEhERj6666iquuuoqT+8pOF7On7bncuOYwQl7EmJDtLQlIuLR3Xff7fk9b2/JobLacdvFGU1XTjBKJCIiHt1xxx2e6jvneOvD/Ywb2ovh/brFKKrgtK35lYiID0pLSyktLW12/fezCtmTX8JtFw+JYVTB0YxERMSj6667DoB169Y1q/5bH+6je8cUrvtS27glSl1KJCIiHn3/+99vdt3CkgpWZeZy+/ghdExNjmFUwVEiERHxyMtNG5duyaGiOtRml7VAeyQiIp4VFRVRVFTUZL3qkOP1v2Uz9sxejBzQ9jbZa2hGIiLi0dSpU4Gm90hWf5zLvsJSHrz2XB+iCo4SiYiIRz/60Y+aVW/eX/aQ0bsTV1/Qdm7QWB8lEhERj2688cYm62zZd5TN2Ud5bMr5JCeZD1EFR3skIiIe5efnk5+f32idX/9lD906pvDNsW3vk+x1aUYiIuLRzTffDDS8R7K/sJR3Mg9xzxXD6NKGbhffkLY/QhGRKPvJT37S6Ovz/2cPSWZMv2yoPwEFTIlERMSjKVOmNPjakWMneOvDfdww6gwG9ujkY1TB0R6JiIhHubm55Obm1vvaK+uzqKx2/PCrI3yOKjiakYiIeHTrrbcCp++RHDl2goXvZ/P1UYM4K71LAJEFQ4lERMSj2bNn11v+H+uyqAo5fvjV4T5HFCwlEhERjyZPnnxa2ZHiE7z5QTb/NHoQQ9vRbAS0RyIi4tn+/fvZv3//KWUvrdvdLmcjoBmJiIhn3/72t4Ev9kiyC0p484Nsbh4zmDP7tK/ZCCiRiIh49vDDD5/y/Jcrd5CanMRPrj4noIiCFReJxMz2AseAaqDKOTfWzHoDi4GhwF7gm865o2ZmwHPAdUApMN05tyWIuEWkfZo0adLJx+9nFbBqey4/ueoc+nXvGGBUwYmnPZJ/dM6Ncs6NjTyfDax1zo0A1kaeA1wLjIh8zQBe9j1SEWnXsrKyyMrKIhRyPPnHjzmjR0fuuWJY0GEFJi5mJA2YCkyMPP4NsA74WaT8deecA943s55mNtA5dyiQKEWk3bnzzjsBuO+ZhWQeKObZW0a12WN0myNeEokD/mRmDnjFOTcX6F+THJxzh8ysX6TuIKD25RI5kTIlEhHxxeOPP87x8ioeX7WDizJ6csNFZwQdUqDiZWnrcufcGMLLVjPN7IpG6tZ3Y393WiWzGWa2ycw25eXlRStOEREmTJjA/5T0I/94OU9OvZCkNn7eSFPiIpE45w5Gvh8Bfg9cDBw2s4EAke9HItVzgNo3+B8MHKynzbnOubHOubF9+/aNZfgi0s7859oP+M3KDdx5+Vl8aXCPoMMJXOCJxMy6mFm3msfA1UAmsByYFqk2DVgWebwc+I6FXQIUaX9ERPxSXlXNjBn3cnzty/yfq9rn5b51xcMeSX/g9+GrekkBfuucW2VmG4HfmdldwD7gG5H6Kwlf+ruL8OW/3/U/ZBFpr1768246XHoHD19/Xrs4tKo5Av9TcM5lARfVU14AXFlPuQNm+hCaiMgpMg8U8dK6XXzz+iu577bRQYcTNwJf2hIRSQSlFVX86K2/06dLB24dDpmZmUGHFDcCn5GIiCSCn//Xx+wpKOHNu8fz0F2Nn9ne3iiRiIg0YVVmLm99uJ/vTTiby85O56mnngo6pLiiRCIi0oj9haXMXrqNLw3qwf+NXKU1bty4gKOKL9ojERFpQGlFFTPe2Ewo5Pj320aTlhL+lbl161a2bt0acHTxQzMSEZF6OOf46ZJt7Mgt5rXp40459XDWrFmA9khqKJGIiNRj7ntZ/Ne2Q/xs8rlMHNnvlNeeffbZgKKKT0okIiJ1rMrM5VerdnD9lwfyvQmn3x5+1KhRAUQVv7RHIiJSywdZBfxo0d+5KKMnT998EZG7bpxi48aNbNy4MYDo4pNmJCIiETtyi7n79U1k9OrE/Gnj6JRW/xkjDzzwAKA9khpKJCIiwN78EqbN/5AuaSm8ftd4enVJa7DuCy+84GNk8U+JRETavay849w+7wMqqkIsmnEpg3p2arT+hRde6FNkiUGJRETatV1HjnP7vPepDjnemnEJIwd0a/I9GzZsAOCyyy6LdXgJQYlERNqtjw8W8535HwKwaMYljOjfdBIBeOihhwDtkdRQIhGRdmn9p3n8YOFmundK5Y27xjO8X9dmv/eVV16JYWSJR4lERNqdxRv38dDvMzmnfzdemz6OAT06enr/yJEjYxRZYlIiEZF2o7I6xJx3dvDrv+7hinP68uLto+nWMdVzO+vXrwdgwoQJ0Q4xISmRiEi7cKT4BDN/u4WNe48y/bKh/PP155Ga3LLPZD/22GOA9khqKJGISJv318/ymbV4KyXlVTx36yimjhrUqvbmz58fpcjaBiUSEWmzyiqq+dWqHSzYsJfh/bry23vGc04zr8xqzLBhp99/qz1TIhGRNmnT3kJ+umQbWfkl3Hn5Wfx08kg6ptZ/yxOv1qxZA8CkSZOi0l6iUyIRkTal4Hg5c97ZwX9uzmFQz0789u7xXDY8Pap9PPnkk4ASSQ0lEhFpEyqqQrz14T7+bfWnlJRX8b0JZ/OjK4fTOS36v+beeOONqLeZyJRIRCShOef440eHeOrdnWQXlHLpsD48MfWCZn9KvSUyMjJi1nYiUiIRkYRUHXKsyszlxT/v4uNDxYyMfLhw4si+9Z4hEk2rVq0CYPLkyTHtJ1EokYhIQimrqOYPWw8w770ssvJLGNa3C8984yK+PnoQyUmxTSA15syZAyiR1FAiEZGEkF1Qwpsf7GPxxv0UlVVy/sDuvHj7GCZfOMC3BFJj0aJFvvYX75RIRCRuHTtRyTsf5bJkcw4f7i0kOcmYfMEApl02lHFDe8V8CashAwYMCKTfeKVEIiJxpfhEJWs/OczKj3JZ/2keFVUhhvXtwgPXjOSmMYM932AxFlasWAHAlClTAo4kPiiRiEigQiHHp0eOsX5nHn/eeYRNe49SFXIM6N6Rb40fwpSLzmB0Rs/AZh/1eeaZZwAlkhpKJCLiq6rqEDtyj7Fl31Hezyrg/axCCksqADh3QDfuuWIYk87rz+iMniT5vPfRXEuWLAk6hLiiRCIiMVNZHWJPfgnbDxaReaCYzANFbMspoqyyGoCBPToy8Zy+XHJ2H74yIp2BPRo/Kz1epKdH95PyiU6JRERaxTlH3vFy9heWsTe/hD2Rr08PH2NPfglVIQdAh5Qkzh3YnVvGZTB6SE/GDOnF4F6d4mrJqrmWLl0KwI033hhwJPFBiURE6hUKOYrKKikoqaCwpIK8Y+XkHy8n71g5h4pOkFtcxqHPT3Dg8zLKq0In35ecZGT06sTwft246vz+jOjflQvO6MGw9C6ktPD8j3jz/PPPA0okNZRIRBJUdchRWR2isjpEVXX4cUV1iIqqEOVV4e8nKqspj3w/URWirKKKsopqSiqqKa2ooqS8mpLyKo6dqOJ4eRXFJyopKgt/FZdVEplMnCI5yejXrQMDenTk3IHdmHR+fwb36sTgXp0Y2qcLGb07t/jAqESxbNmyoEOIK0okjfjkUDHz3ssKOgwA6vn3HAjnvEXSnNoNNdnQexuLwTXwxNV64twXfZ5WflpMrlbdcN+uVt2aWMLP3cm2az8OnXxP+HvIhR9Xh1ydx+Hn1aEvntd8rwo5qqvD36tCIapCrsE/t+ZKMuiSlkLXjil065hC1w4p9O6SxlnpXejeMZWenVPp3SWN3l3S6NU5jb7dOpDetQO9u6T5/gHAeNOjR4+gQ4grCZtIzGwy8ByQDLzqnJsT7T6KyirZmF0Y7WZbzIiPf7xel7SbU72hdfIG39tIo7Vfqt3uqeVf/HnW7brmPVa7rtW0YZH3hl+wk22F32e12jaDpCRIsaRTy8xIqvmeVPexkWyQlGQkW/h5SrKRkhR+PTU5ieSk8POUpCRSU4zUpCRSk43UlCRSk5LokJpEWnISaSlJdExNpkNKEh1SkumUlkSntBQ6pSbTOS1cnoj7E/Fg8eLFANxyyy0BRxIfzOv/MOOBmSUDnwJXATnARuA259zH9dUfO3as27Rpk48RikhbNnHiRKDtn9luZpudc2ObqpeoM5KLgV3OuSwAM1sETAXqTSQiItG0cuXKoEOIK4maSAYB+2s9zwHGBxSLiLQznTt3DjqEuJKol1bUt7B7yhqdmc0ws01mtikvL8+nsESkPVi4cCELFy4MOoy4kaiJJAeofUTZYOBg7QrOubnOubHOubF9+/b1NTgRadteffVVXn311aDDiBuJurS1ERhhZmcBB4BbgduDDUlE2ovVq1cHHUJcSchE4pyrMrP7gHcJX/473zm3PeCwRKSdSE1NDTqEuJKQiQTAObcS0KUTIuK7BQsWADB9+vRA44gXibpHIiISmAULFpxMJpKgH0j0yszygGygB1BU66XGntc8TgfyoxRK3f5aWreh1+orb84Y674WxJibqheLMdd+rDG3XLT+Xjf2elP/duuWtYfFo/8pAAAEZUlEQVQxN/dn3poxn+mca/pqJedcu/kC5jb3ec1jYFOs+m9p3YZeq6+8OWOMhzE3VS8WY67zWGOO8XijOeam/gzaw5g9/Myj9ne7oa/2trS1wsPzuq/Fov+W1m3otfrKvYwxyDE3VS8WY47FeL2021bGHK2/14293py/q+1tzH7//mpQu1jaag0z2+Saca+ZtkRjbh805vbBjzG3txlJS8wNOoAAaMztg8bcPsR8zJqRiIhIq2hGIiIiraJEIiIiraJEIiIiraJE4pGZDTOzX5vZkqBj8YuZfd3M5pnZMjO7Ouh4/GBm55nZf5jZEjP7ftDx+MXMupjZZjP7WtCx+MHMJprZXyI/64lBxxNrZpZkZr8ws383s2nRaleJBDCz+WZ2xMwy65RPNrOdZrbLzGYDOOeynHN3BRNp9Hgc8x+cc/cA04GEPaTa45g/cc59D/gmkLCXi3oZc8TPgN/5G2V0eRyzA44DHQkfT5FwPI53KuGDASuJ5nhj/YnHRPgCrgDGAJm1ypKB3cAwIA34X+D8Wq8vCTruAMb8DDAm6Nj9GjNwA7ABuD3o2P0YMzCJ8JEM04GvBR27T2NOirzeH3gz6Nh9GO9s4N5Inaj9DtOMBHDOvQcU1ik+eS68c64CqDkXvk3wMmYL+xXwjnNui9+xRovXn7Nzbrlz7jLgW/5GGj0ex/yPwCWEz/a5x8wS8veDlzE750KR148CHXwMM2o8/oxzCI8VoDpaMSTsbeR9UO+58GbWB/gFMNrMHnTO/TKQ6GKj3jEDPyT8v9UeZjbcOfcfQQQXIw39nCcCNxL+5dLWjiuod8zOufsAzGw6kF/rl2xb0NDP+UbgGqAn8EIQgcVIQ/+WnwP+3cy+ArwXrc6USBpW77nwzrkC4Ht+B+OThsb8PPC838H4pKExrwPW+RuKb+od88kHzi3wLxTfNPRzXgos9TsYHzQ03lIg6nu8CTl19UmT58K3QRqzxtxWtbcx+zpeJZKGnTwX3szSCG9CLg84pljTmDXmtqq9jdnX8SqRAGb2FvA3YKSZ5ZjZXc65KqDmXPhPgN+5NnQuvMasMaMxt4kxx8N4ddNGERFpFc1IRESkVZRIRESkVZRIRESkVZRIRESkVZRIRESkVZRIRESkVZRIRESkVZRIRESkVZRIRHxkZv9kZs7Mzo08H1rPgUT/Ymb3BxOhiHdKJCL+ug34K+F7H4m0CUokIj4xs67A5YRv461EIm2GEomIf74OrHLOfQoUmtmYoAMSiQYlEhH/3Eb4yFMi32+j1oFSdehuqpIwdEKiiA8iRzR/FbjQzByQTDhZPAH0qlO9N7DH3whFWk4zEhF/3Ay87pw70zk31DmXQThZjAIOmdmVAGbWG5hMeENeJCEokYj44zbg93XK3gZuB74DPGxmW4H/Bh53zu32OT6RFtPBViIi0iqakYiISKsokYiISKsokYiISKsokYiISKsokYiISKsokYiISKsokYiISKsokYiISKv8f9uD0sbbO7QuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fake_αcens = ICRS(ra=αcen.ra, dec=αcen.dec, \n",
    "                  radial_velocity=[αcen.radial_velocity]*len(distances_to_try), \n",
    "                  distance=distances_to_try)\n",
    "fake_αcens_gcrs = fake_αcens.transform_to(observer_gcrs)\n",
    "plt.semilogx(distances_to_try, np.sqrt(np.sum(fake_αcens_gcrs.velocity.d_xyz**2, axis=0)))\n",
    "plt.axvline(1*u.pc, c='k', ls=':')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmm. That seems pretty close for the onset of the effect?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fix by converting observer to ICRS? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<CartesianDifferential (d_x, d_y, d_z) in km / s\n",
       "     (-0.38410488, 0.21045191, 0.00070133)>,\n",
       " <CartesianDifferential (d_x, d_y, d_z) in km / s\n",
       "     (-30.31439059, 4.35992313, 1.80095324)>)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "observer_icrs = observer_gcrs.transform_to(ICRS)\n",
    "observer_gcrs.velocity, observer_icrs.velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "~30 km/s is about right, as that's roughly the orbital speed of the earth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$[0.0011716038,~0.0076441094,~0.00077090226] \\; \\mathrm{\\frac{cm}{s}}$"
      ],
      "text/plain": [
       "<Quantity [0.0011716 , 0.00764411, 0.0007709 ] cm / s>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "roundtrip_delta_v = observer_icrs.transform_to(observer_gcrs).velocity - observer_gcrs.velocity\n",
    "roundtrip_delta_v.d_xyz.to('cm/s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The velocity round-trips at much better than the cm/s level, which is roughly the limit of the best RV observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ICRS Coordinate: (ra, dec, distance) in (deg, deg, lyr)\n",
       "    (219.90085, -60.83562, 4.37)\n",
       " (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)\n",
       "    (0., 0., -18.)>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "αcen_icrs = αcen.transform_to(observer_icrs)\n",
    "αcen_icrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-26.397423 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity -26.39742346 km / s>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_pos = (αcen_icrs.data.without_differentials() - observer_icrs.data.without_differentials()).to_cartesian()\n",
    "pos_hat = d_pos/d_pos.norm()\n",
    "\n",
    "d_vel = αcen_icrs.velocity - observer_icrs.velocity\n",
    "\n",
    "rv = np.sum(d_vel.d_xyz* pos_hat.xyz, axis=0)\n",
    "rv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aha, that seems much more plausible.  Lets just double-check the answer depends on time (as it should because the earth's velocity vector moves over time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-26.397423 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity -26.39742346 km / s>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "observer_gcrs2 = keck.get_gcrs(obstime)\n",
    "observer_icrs2 = observer_gcrs2.transform_to(ICRS)\n",
    "\n",
    "d_pos = (αcen_icrs.data.without_differentials() - observer_icrs2.data.without_differentials()).to_cartesian()\n",
    "pos_hat = d_pos/d_pos.norm()\n",
    "\n",
    "d_vel = αcen_icrs.velocity - observer_icrs2.velocity\n",
    "\n",
    "rv = np.sum(d_vel.d_xyz* pos_hat.xyz, axis=0)\n",
    "rv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1248cc160>]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEMCAYAAADTfFGvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd4lFXC/vHvSSehhJBQQhISSkKvkY4VBQsitrV3XVfdYtdl1V1dddXfru+uq7i6umv3VRFFFLEiICK9EyCEQOgJLY308/sjw27kDWVS5swk9+e65iLzzEzmzpOEO087x1hrEREROVFBrgOIiEhgUXGIiIhXVBwiIuIVFYeIiHhFxSEiIl5RcYiIiFdUHCIi4hUVh4iIeEXFISIiXglxHaAxxMbG2uTkZNcxREQCypIlS/KstXHHe16TLI7k5GQWL17sOoaISEAxxmw5kedpV5WIiHhFxSEiIl5RcYiIiFdUHCIi4hUVh4iIeEXFISIiXlFxiIiIV1QcIlIvq7Yd5PUfstm4u4CjTUV9qKySldsOUFFZ5dtw0iia5AWAItL4SsorefbLDbw8N4sqT190jm7BmB6x9OzYitQOrQgNCeLDpduZsWIHBaUVxLeJ4LKhSVx2UiLtW0e4/QKkzszR/kIIZOnp6VZXjos0Dmstczfm8fvpa8jKK+LyoUncODqZRdn7+TZjDwuz93GguPw/z28RGszZ/ToyLCWGGSt3MndjHqHBhnvHpXHzmK4YYxx+NVKTMWaJtTb9uM9TcYjI0VhrKa+0BAcZggx8u34Pf/s6k+U5B0ho24KnLurPqO6x/+c1uYWlbNxdyIHick5OjaVVROh/Hs/OK+JPMzP4fM0uxvZqz/+7ZADRkWG+/tKkFioOFYdIveTsK+YXby1h9fb8nyzvHN2CX5zajUvSEwgPCa7T57bW8tr8bB7/bB2xLcNJT46hVUQI0S1CObd/J/rEt2mIL0G8pOJQcYjU2fxNedz+1lIqqyzXjUoh2BgqrSW5XSQTBsQTGtww59WsyDnAE5+tY3d+CQUlFRw8VE5FlWVcnw78+oxUese3bpD3kROj4lBxiNTJWz9u4eGP15ASG8XL16STEhvls/c+eKicf32/mVfmbaagpILfnduLm8Z09dn7N3cnWhw6HVdE/uO9RTlMnraaU1LjmHbbSJ+WBkCbFqH8Zmwq8+4/nXF9OvD4Z+uYvX6PTzPI8ak4RASAWWt28cCHKzk5NY4XrxrykwPavtamRSjP/mwgaR1a8at3lpGdV+Qsi/xfKg4R4fvMPH75zjIGJEbz4lWDCQtx/19DZFgIL1+TTlCQ4ZY3FlNYWuE6kni4/+kQESfyS8p5d+FWLnlxPlf+80dS2kXxr+tOIjLMf64LToyJ5O+XDyZzTyEXT5nP8pwDriMJKg6RZml+Zh5jnvqWBz5cxd6iMu4dl8a7twz3y+spRveI5ZVrT+JAcTkXvvA9j36yliJtfTils6pEmpl3F27ldx+tJiU2iqcu7s+gxOiAuHo7v6Scpz/P4M0FW4lvE8HDE3ozrk/HgMgeKHRWlYj8hLWWJz9bxwMfrmJU91im3jaSwUltA+Y/3tYRofzxgn58cOsIWrcI5dY3l3LdvxaRs6/YdbRmR8Uh0kxMW7adf8zJ4qrhSbxybTqtHZ41VR/pyTHM+OVoHj6vN0u27OfqV37kUFml61jNSkAUhzFmvDFmvTEm0xjzgOs8IoFmT0EJf/hkLeld2vLo+X0JaaArv10JCQ7ihtEpvHTNELL3FvOXL9e7jtSs+P1PjzEmGHgeOBvoDVxujOntNpVI4LDW8tBHqzlUXslTF/cnKCgwdk2diJHdYrlyWBKvzNvM0q37XcdpNvy+OIChQKa1NstaWwa8C0x0nEkkYHy6aiez1uzmrjNT6RbX0nWcBvfA2T3p2DqC+z5YSWmFdln5QiAUR2cgp8b9bZ5lP2GMucUYs9gYszg3N9dn4UT8WVZuIY98vIYBCW24aXSK6ziNolVEKE9e1J/MPYX8z1cbXcdpFgKhOGrbrv4/5xBba1+y1qZba9Pj4uJ8EEvEvy3Zso+LpswH4M+XDgj44xrHckpqHJcPTWTK7E18sWaX6zhNXiD8JG0DEmvcTwB2OMoiEhA+X72TK17+kejIMD68bSTd27dyHanRPTKhDwMS2nDn/y5n4+4C13GatEAojkVAD2NMijEmDLgMmO44k4jfem1+Nr94ayl94lsz9Rcj6dLOtyPcuhIRGsyLVw+hRVgIN7++mIM1pq+VhuX3xWGtrQDuAGYB64D3rLVr3KYS8T/WWp6ZlcEj09cwtlcH3rppODFR/jeESGPq1KYFL141mO0HDnHLG4s5UFzmOlKT5PfFAWCt/cxam2qt7Watfdx1HhF/U1FZxX0frOT5bzdx+dAkplw5mBZhdZvWNdClJ8fwzMUDWLb1AOc9N4/V2w+6jtTkBERxiMix/Xt+Nu8v2cavz+jBE5MC/wK/+rpgUGfeu3UElVWWi6bM58Ol21xHalKa90+XSBNwoLiM577J5OTUOO48MzVgxp5qbAMTo5nxy9EMSorm3g9WkrEr33WkJkPFIRLgnvsmk4KSciaf08t1FL/TrmU4U64cQuuIEH43bTVVVU1vNHAXVBwiASw7r4jXf8jm0vRE0jo2/VNu66JtVBgPnt2LxVv2M1W7rBqEikMkgD09K4PQ4CDuOjPVdRS/dvGQBIZ0acuTMzN0plUDUHGIBKgfs/by2apd/PzkbrRvHeE6jl8LCjL88YK+HDxUztOzNJJufak4RAJQYWkF936wkoS2Lbj55KY5BlVD69WpNdeOSOadhVtZu0MHyutDxSESgP44Yy05+4t59mcDiQwLcR0nYPz6jB60jgjlyZnrXEcJaCoOkQDz5drdvLsoh1tP6cZJyTGu4wSUNpGh/PL07szdmMd3GzSKdl2pOEQCSG5BKQ9MXUmvTq25c6wOiNfF1SO6kBjTgic/W0elTs+tExWHSIA4UFzGdf9aSEFpBf/zs4GEhejXty7CQ4K5f3xPMnYV6PTcOtJPnkgAOFBcxpX//JGNewp56eohumajns7t14mBidH8+Yv17C0sdR0n4Kg4RPzckaVxalp715ECnjGGRyf24UBxOTf8exFFpRWuIwUUFYeIn3v44zVs3K3SaGj9E6L5+xWDWbX9IL94ayllFVWuIwUMFYeIH1u17SDTV+zg5pNTVBqN4MzeHXjywn7M2ZDL/VNXYq0Olp8IFYeIH3vq8wzaRoby81O6uY7SZP3spCTuPjOVacu2M23ZdtdxAoKKQ8RPzdmQy7zMPO44vfqiNWk8t5/WnUFJ0Tz+6TpNOXsCVBwifqiqyvKnmRkktG3BVcOTXMdp8g6PZbW/uIxnvshwHcfvqThE/ND0FTtYuzOfe8elER7SPKeA9bU+8W24dmQyb/24leU5B1zH8WsqDhE/U1ZRxV++3EDvTq2Z0D/edZxm5a4zU4lrGc7vPlqlq8qPQcUh4mfeW5zD1n3F3DsujaAgTQPrS60iQpl8bi9Wb89nxsodruP4LRWHiB8pKa/kuW82kt6lLaemxbmO0yxN6B9P19goXp6bpdNzj0LFIeJH3vhhC7vzS7lnXBrGaGvDhaAgw01jurJ6ez4Lsva5juOXVBwifqKgpJwXZmcypkcsw7u2cx2nWbtwcGfaRYXx8tws11H8kopDxE/847ss9heXc89Zaa6jNHsRocFcMyKZbzL2sHF3ges4fkfFIeIHvlq7m+dnZ3LhoM4MSIx2HUeonrcjPCSIf87d7DqK31FxiDi2flcBv353GX3j2/D4pH6u44hHTFQYl6QnMG3Zdvbkl7iO41dUHCIO7S0s5cbXFhEVHsLL16TTIkwX+/mTm8d0BeCR6Wt0hlUNKg4RR6y13P3+CvYUlPLSNel0bBPhOpIcoUu7KH5zZg9mrt7FjJU7XcfxGyoOEUdmb8hl9vpc7huXxkAd1/Bbt4zpyoDEaB7+eDW5BZotEFQcIk5UVFbx+KfrSG4XyTUjkl3HkWMICQ7iz5f0p6iskt99tEq7rFBxiDjxzqIcMvcU8uA5vQgL0a+hv+vevhV3n5nKrDW7+US7rFQcIr6WX1LOs19uYFhKDGf17uA6jpygmzy7rP4wfQ37i8pcx3FKxSHiY89/m8n+4jIeOq+3hhUJIMFBhj9d2I+Dh8r546frXMdxSsUh4kMHi8t5ff4WJg6Ip2/nNq7jiJd6dWrNrad0Y+rSbczdmOs6jjN+URzGmEuMMWuMMVXGmPQjHnvQGJNpjFlvjBnnKqNIQ3h30VYOlVdy88ldXUeROrrj9O50jY3it9NWUVxW4TqOE35RHMBq4EJgTs2FxpjewGVAH2A88IIxRldISUCqqKzitfnZDO8aQ594bW0EqojQYJ68sB85+w7x6rzmORyJXxSHtXadtXZ9LQ9NBN611pZaazcDmcBQ36YTaRiz1uxmx8ESbhiV4jqK1NOwru0Y3T2Wt37cSkVlles4PucXxXEMnYGcGve3eZaJBJxXv99MUkwkZ/TSmVRNwdUjurDzYAlfrdvjOorP+aw4jDFfGWNW13KbeKyX1bKs1qtvjDG3GGMWG2MW5+Y234NW4p+W5xxgyZb9XDcymWBNB9sknNGzPfFtInhjQbbrKD4X4qs3staOrcPLtgGJNe4nALVOBGytfQl4CSA9PV2XdopfeXXeZlqGh3BJeoLrKNJAQoKDuHJ4F56ZtZ7MPYV0b9/SdSSf8fddVdOBy4wx4caYFKAHsNBxJhGvLM7exycrd3DFsCRaRYS6jiMN6NL0REKDDW8u2OI6ik/5RXEYYyYZY7YBI4BPjTGzAKy1a4D3gLXA58Dt1tpKd0lFvFNSXsl9U1cS36YFvzqjh+s40sDiWoVzTr9OTF2yjaLS5nNqrl8Uh7V2mrU2wVobbq3tYK0dV+Oxx6213ay1adbamS5zinjrr19vJCu3iCcv7EfLcJ/tGRYfumZEFwpKK/hw2XbXUXzGL4pDpClavf0gL83J4pIhCZycGuc6jjSSwUltGZgYzZRvMymtaB47RFQcIo2gorKKez9YSUxUGL87t7frONKIjDHcc1YaOw6W8O7CnOO/oAlQcYg0gn/Pz2bdznwePb8PbSJ1QLypG9W9HcNSYvj7t5kcKmv6Wx0qDpEGtvPgIZ79cgOnpcUxvm9H13HEB4wx3H1WGrkFpc3iug4Vh0gD++OMdVRUWf5wfl8Nm96MDE2J4eTUOKbM3kRhEz/DSsUh0oC+25DLp6t2csdp3UlqF+k6jvjY3Wemsr+4nJfmZLmO0qhUHCINpLyyikc+Xk3X2ChuOUXDpjdHAxKjOa9/J16cvYnMPQWu4zQaFYdIA/lq7W6y9xZz/9k9CQ/R6P/N1SMT+hAZHsx9H6yksqppjn6k4hBpIG8s2ELn6BaM1ei3zVpcq3AePq83S7ce4I0fsl3HaRQqDpEGkLmngPmb9nLFsCSNfitMGtSZk1PjeHrWerbtL3Ydp8GpOEQawJsLthIWHMTPTko8/pOlyTPG8MSkvgA8+VmG4zQNT8UhUk/FZRVMXbKNc/p1JLZluOs44icS2kZy9fAufL5mFzsPHnIdp0GpOETq6ePlOygoreDqEV1cRxE/c9XwLlRZy9s/bnUdpUGpOETqwVrLGz9soVen1gxOaus6jviZxJhIzujZnncWbm1SAyCqOETqYc7GPNbuzOeaEV10lbjU6uoRyeQVljFz1S7XURqMikOkjqqqLE/NzCAxpgUXDdaUsFK7Md1jSYmN4vUfsl1HaTAqDpE6+mTlDtbuzOfuM9MIC9GvktQuKMhw9fAuLN16gFXbDrqO0yC8/mk3xtzfGEFEAklZRRV//mIDvTq15vwB8a7jiJ+7aEgCkWHBTWar47jFYYx5r8btfeAmH+QS8WvvLtrK1n3F3Dc+jSBd8CfH0aZFKJMGdebjFTvYW1jqOk69ncgWR7619lLP7RLgq8YOJeLPCksr+NvXGxmWEsOpmhJWTtD1o5Ipq6jinYWBf2ruiRTH40fcn9wYQUQCxZ9mrmNvURkPntNLZ1LJCevevhVjesTy+g9bKKuoch2nXo5bHNbazQDGmFjP/X2NHUrEX83PzOPNBVu5cVQKAxOjXceRAHPDqBT2FJQyc/VO11HqxZuD4682WgqRAFBUWsF9U1eSEhvF3WeluY4jAeiU1Di6xkbx6vfZrqPUizfFoW1yadae+jyD7QcO8czF/WkRpvk2xHtBQYZrRyazIucAS7fudx2nzrwpjqY5I4nICViRc4DXf9jCDaNSSE+OcR1HAthFQxJoFR7CK3M3u45SZ9riEDkBr8zbTKvwEO48M9V1FAlwLcNDuG5UMp+u2sn8zDzXcerEm+J4sNFSiPix3fklfLZqJ5eelEjL8BDXcaQJuP207iS3i+TBaasoKQ+8wQ9PuDistauNMZcYY1oBGGN+Z4z50BgzuPHiibj31o9bqbSWazRsujSQiNBgnpjUjy17i/nr1xtdx/Gat0OOPGStLTDGjAbGAa8BUxo+loh/KK2o5O0ft3B6Wnu6tItyHUeakJHdY7lkSAIvzclizY7AGsPK2+I4vE11LjDFWvsxENawkUT8x2erdpJXWMa1I5NdR5EmaPK5vWgbGcrkaauxNnDOP/K2OLYbY/4BXAp8ZowJr8PnEAkY//4+m25xUYzpEes6ijRB0ZFh3H1WGstzDjAvgA6Ue/uf/qXALGC8tfYAEAPc2+CpRPzA0q37WbHtINeOTNbQItJoLhzcmfatwnnxu02uo5wwr4rDWltsrf3QWrvRc3+ntfaLxokm4tY/52bROiKECzVJkzSi8JBgbhydwveZe1m57YDrOCdEu5lEarFlbxGfr97FlcO76BRcaXRXDEuiVURIwGx1qDhEavHKvM0EBxmu10Fx8YFWEaFcPbwLM1fvIiu30HWc4/KqOIwx6caYacaYpcaYlcaYVcaYlY0VTsSF/UVlvLc4hwsGdqZ96wjXcaSZuH5UCqHBQbw8N8t1lOPydovjLeBfwEXABOA8z7/1Yox5xhiT4SmjacaY6BqPPWiMyTTGrDfGjKvve4kczxsLtlBSXsXNJ3d1HUWakbhW4Vw8JIGpS7azr6jMdZxj8rY4cq210621m621Ww7fGiDHl0Bfa21/YAOe4U2MMb2By4A+wHjgBWOMhiWVRlNSXslr87M5LS2O1A6tXMeRZuaaEV0oq6ziw6XbXEc5Jm+L4xFjzD+NMZcbYy48fKtvCGvtF9baCs/dBcDh01gmAu9aa0s9E0plAkPr+34iR/PGD1vYW1TGLSd3cx1FmqGeHVszKCmadxZu9esLAr0tjuuBgVT/9T+B/+6uakg3ADM9H3cGcmo8ts2zTKTBbcot5P99sZ7Te7ZneFcNnS5uXD40iU25RSzK9t/5Orw9z3CAtbZfXd7IGPMV0LGWhyZ7hi7BGDMZqKD6WArUPpR7rTVsjLkFuAUgKSmpLhGlGausstzz/goiQoP504X9dMGfOHNe/0489sla3lm4laEp/vkHjLdbHAs8xx28Zq0da63tW8vtcGlcS/XWy5X2v9to24DEGp8mAdhxlM//krU23VqbHhcXV5eI0oy9PDeLZVsP8OjEPjqTSpyKDAvhgkGd+XTVTg4U++dBcm+LYzSw3HOGU4OdjmuMGQ/cD5xvrS2u8dB04DJjTLgxJgXoASys7/uJ1LRxdwF/+WID4/p04PwB8a7jiHD50CTKKqr4cOl211Fq5e2uqquB3COWnd8AOf4OhANfenYRLLDW3mqtXWOMeQ9YS/UurNuttYE364n4tac+X09EaBCPT9IuKvEPveNbMyAxmrcXbuX6Uf43Vpq3WxwvAa1rnIY7AriqviGstd2ttYnW2oGe2601HnvcWtvNWptmrZ15rM8j4q2MXfl8tW43149KIbZluOs4Iv9x7YguZO4pZNoy/9vq8LY4LgZeM8b0MsbcDNwOnNXwsUR84/lvNxEVFsz1o5JdRxH5iQsGdmZAQhuenJlBQUm56zg/4e3ouFlUX5A3leoSOctaG1hTV4l4bM4r4tOVO7hqeBeiIzUfmfiXoCDDHyb2JbeglOe+yXQd5ydO6BiHMWYVPz0NNgYIBn40xuC54lskoEyZnUlIcBA3jklxHUWkVgMTo7k0PYFX523m0vREurdv6ToScOIHxxv6Ij8Rp7YfOMSHS7dzxbAk2rfS6bfiv+4b35OZq3fxh0/W8PoNQ/3iQPkJ7aqqOS5VbbfGDinS0F6ZuxmAWzSQofi52Jbh3Dk2lbkb8/h2/R7XcQDNxyHNUFFpBe8vzuGcfp1IaBvpOo7IcV09ogvJ7SJ5auZ6Kqvcj2Gl4pBmZ9qy7RSUVnDtyC6uo4ickNDgIO4d15P1uwuY6gcj56o4pFmx1vL6D9n07dyawUltXccROWHn9OvIgMRonv1yAyXlbq+DVnFIs7Igax8bdhdyzQj/uxpX5FiMMTx4dk92HizhX99nO82i4pBm5fUfsomODNWYVBKQhndtx+k92/PC7EwOFru7KFDFIc3GjgOH+GLtbn52UiIRoZpIUgLTXWemUlBSwYfL3B3rUHFIs/Ha/GystVw1TAfFJXD17dyGAYluZwlUcUizsHZHPq/M28wFAzuTGKNTcCWwXX5SIht2F7J06wEn76/ikCavvLKKe95fQXRkGA+dV6d5yET8yoQB8USFBfPOwq1O3l/FIU3elNmbWLsznz9e0Je2URrMUAJfVHgI5w/szIyVOzh4yPcHyVUc0qRl7MrnuW82MmFAPOP71jblvUhgumJoEiXlVUxf7vv5OlQc0mRVVVkemLqK1hGh/OH8Pq7jiDSofglt6BPfmrd+9P1BchWHNFkfLtvO8pwD/PacXsRoF5U0QZcPTSJjVwGLt+z36fuqOKRJyi8p508zMxiUFM2kQZ1dxxFpFBcO7kxcq3Ce/jzDp1sdKg5pkp77eiN7i0r5w/l9CArS0CLSNEWGhfCbsT1YlL2fr9f5bsh1FYc0OZl7CvnX99n8LD2R/gnRruOINKpL0xPpGhvFU59nUFFZ5ZP3VHFIk/PkZ+toERbMPePSXEcRaXShwUHcNz6NjXsKfTbkuorjCPuKylxHkHrYlFvI1xl7uGl0V2JbhruOI+IT4/p0ZFBSNM9+uZFDZY0/5LqKo4ZnZmVw3t/mkl/ibtRJqZ83fthCWHAQVwxLch1FxGeqh1zvxa78EqavaPzrOlQcNYzt1YFd+SU8+dk611GkDgpKynl/cQ7n9u9EXCttbUjzMjQlhvd+PoJL0xMb/b1UHDUMSmrLzSd35Z2FOczZkOs6jnhp6pJtFJVVct3IZNdRRJwYmhLjkwnKVBxHuHNsKt3bt+SBqSu1yyqAVFVZXvthCwMToxmQqDOpRBqTiuMIEaHBPHNxf3bll/DEp9plFSjmbMxlc14R149Kdh1FpMlTcdRiUFJbbhrTlXcX5bB+V4HrOHICXv0+m7hW4Zzdt5PrKCJNnorjKG49pRsRoUG8Om+z6yhyHN9k7GbOhlxuHJ1CWIh+pEUam37LjiImKowLBycwbfl28gpLXceRozhUVsnDH6+he/uW3DAqxXUckWZBxXEMN4xKoayiircWuJllS47v799uZNv+Qzw2sa+2NkR8RL9px9C9fUtOTYvjjQVbKK1o/KsxxTuZewp4aU4WFw7qzIhu7VzHEWk2VBzHcePoFPIKS5m+fIfrKFJDVZXldx+tpkVoML89t5frOCLNiorjOEZ3jyWtQytembfZ57NsydFN+W4TC7L2MfncXhqTSsTHVBzHYYzhxjEpZOwqYPZ6XU3uDxZl7+MvX25gwoB4nwyvICI/5RfFYYx5zBiz0hiz3BjzhTEm3rPcGGP+ZozJ9Dw+2EW+SYM6kxQTyZ+/XK+tDsf2FZXxy7eXkdC2BU9M6uuT4RVE5Kf8ojiAZ6y1/a21A4EZwMOe5WcDPTy3W4ApLsKFBgfxqzN6sHp7PrPW7HYRQQBrLfd9sIJ9RWU8f8VgWkWEuo4k0iz5RXFYa/Nr3I0CDv9ZPxF43VZbAEQbY5xcGnzBwHi6xkbx7JcbqKrSVocLC7L28dW6Pdx9Vip9O7dxHUek2fKL4gAwxjxujMkBruS/WxydgZwaT9vmWeZzIcFB/HpsD9bvLuDTVTtdRGj2XpidSWzLcK7V6LciTvmsOIwxXxljVtdymwhgrZ1srU0E3gLuOPyyWj5VrX/uG2NuMcYsNsYszs1tnIPYE/rHk9qhJf/z1Qafze0r1VbkHGDuxjxuGpNCRGiw6zgizZrPisNaO9Za27eW28dHPPVt4CLPx9uAmqfNJAC1XlBhrX3JWpturU2Pi4tr+C8ACAoy3HVmGptyi3h61vpGeQ+p3QuzM2kdEcKVmtlPxDm/2FVljOlR4+75QIbn4+nANZ6zq4YDB621TvcTje/bkWtGdOGlOVm8vzjn+C+Qetu4u4BZa3Zz3chkHRAX8QMhrgN4/MkYkwZUAVuAWz3LPwPOATKBYuB6N/F+6qHzerMpt5DJ01aTEhtFenKM60hN2pTZm2gRGsx1GsRQxC/4xRaHtfYiz26r/tbaCdba7Z7l1lp7u7W2m7W2n7V2seusUH167gtXDKFz2xb8/I0l7M4vcR2pycrcU8jHK3Zw+dAkYqLCXMcREfykOAJRm8hQXr4mnYKSCp7+XMc7GoO1lsdmrCUyLJjbTuvmOo6IeKg46qF7+5ZcPyqZD5dtY/X2g67jNDnfZOzhuw25/PqMHhqPSsSPqDjq6bbTuhPdIpTHP12n4UgaUGlFJY/NWEu3uChdtyHiZ1Qc9dSmRSi/GZvKD1l7+XrdHtdxmox/f59N9t5iHjqvN6HB+jEV8Sf6jWwAVwxLomtcFE/MXEe5Lgyst8w9hTz3TSZn9GzPqWntXccRkSOoOBpAaHAQD57di6zcIt7TtR31sim3kMtfXkBEaBCPTOjjOo6I1ELF0UDG9mrPkC5tee7rTErKNc1sXWzOK+KKlxdQVWV55+bhJLWLdB1JRGqh4mggxhjuOSuNXfklvPXjVtdxAs6e/BIuf2kB5ZWWt28eTo8OrVxHEpGjUHE0oBHd2jGqezsnjTBAAAANKklEQVSmzM6kqLTCdZyA8uiMtewrLuPNG4eR1lGlIeLPVBwN7O6z0sgrLOPf87NdRwkYczfmMmPlTm47tRu941u7jiMix6HiaGCDk9pyRs/2/OO7TRw8VO46jt8rKa/k4Y/XkNwukltP0dXhIoFAxdEI7jorlYLSCu5+bzllFTo991j+8V0Wm/OKeHRiX82zIRIgVByNoE98Gx6d2Jev1u3hV+8s07UdR7FlbxHPz87k3P6dODm1ceZQEZGGp+JoJFcP78LD5/Xm8zW7uOu9FZoxsBaPzVhHSJDhoXN7u44iIl7wl/k4mqQbRqdQXlnFkzMzaNMihD9e0M91JL8xe/0evlq3m/vH96RjmwjXcUTECyqORvbzU7qxt6iMl+Zk0Se+DZcP1dSnZRVVPDpjLcntIrlhdLLrOCLiJe2q8oH7x/dkTI9YHv54NUu27HMdx7nX5meTlVvEwxN6Ex6iA+IigUbF4QPBQYa/Xz6Y+OgW3PrmUnYdbL4zBu46WMJfv97IaWlxnN6zg+s4IlIHKg4fOTxjYFFpBb+dtsp1HCdWbTvIpBe+p7LK8tB5OiAuEqhUHD6U2qEVvz6jx39mtmtOPly6jYtfnE+QMbx/6wi6xrV0HUlE6kjF4WPXjUqmS7tI/jhjbbM5Rff5bzO5670VDE5qyye/HE3fzm1cRxKRelBx+Fh4SDC/PacXG/cU8vbCpj+K7psLtvDMrPVcMDCe128cSkxUmOtIIlJPKg4HzurdgZHd2vGXLzdwsLjpjmf1yYodPPTxasb2as8zlwzQFLAiTYR+kx0wxvDQeb3JP1TO5I9WNcldVvMz87jrveWc1CWGv18xWKUh0oTot9mRXp1ac8+4NGas3Mmtby5pUrMGFpZWcM/7K+jSLop/XpeuwQtFmhgVh0O3ndqdxyb24euMPVz9yo9NZhj2pz/PYGd+CU9f3J/WEaGu44hIA1NxOHb1iGT+dtkgluccYMJz85i3Mc91pHpZlL2PNxZs4fqRKQxOaus6jog0AhWHH5gwIJ63bx5OcJDhqld+5J73V3CguMx1LK+VlFdy/9SVdI5uwT3jUl3HEZFGouLwEyclxzDz12O47dRuTFu2nYumzKcwgOYtLygp594PVpKVW8QTk/oRGabxM0WaKhWHH4kIDea+8T157fqhbM4rYvK0VVhrXcc6rvmZeYz/n7l8unIH945L06RMIk2c/iz0Q6N7xHLn2FT+/OUGhqbEcOWwLq4j1Sort5Apszfx/pJtdI2NYuovRjJIxzVEmjwVh5+6/bTuLNqynz98spYBCdF+NUxH5p4Cnv1qI5+t2klYcBA3j0nhrjPTaBGm025FmgPtqvJTQUGGZy8dQExkGDe9tpjV2w+6jgTAptxCLn7xB+asz+UXp3Rj3v2nM/nc3ioNkWZExeHH2rUM59XrTsIYuPjF+XyyYofTPHmFpVz/r0UEG8OnvxrDfeN7Etcq3GkmEfE97aryc73jWzP9jtH84s0l/PKdZczbmMfApGgS20aS2rEl7Vv5Zr7ukvJKbnptMXsKSnjn5uEktYv0yfuKiP9RcQSAuFbhvH3zcB6dsYZ3F+bwv4tzAAgLDuLhCb25clgSxphGe/+cfcU8+OEqVmw7wJQrh+gAuEgzZ/zpdE9jzD3AM0CctTbPVP9v+FfgHKAYuM5au/R4nyc9Pd0uXry4ccM6UlFZxa78EnL2HeIfczYxe30u5w+I54kL+9EyvP5/B+QVllJYUkHLiBBCggyvzNvMP+ZkEWwMvz+/Nz87KakBvgoR8UfGmCXW2vTjPc9vtjiMMYnAmUDNSSrOBnp4bsOAKZ5/m62Q4CAS2kaS0DaSYSkxTPluE3/+Yj2rtx/kr5cNol9C3c6+Ki6r4LlvMvnn3CzKK3/6x8TEgfE8cHZPOrVp0RBfgogEOL8pDuBZ4D7g4xrLJgKv2+rNogXGmGhjTCdr7U4nCf1MUJDh9tO6MzipLXf+73ImvfA9d56Zyq2ndANg3c58MnYVEBkWTHRkKLEtw+nRvuVPdmtZa/l89S4em7GWHQdLuGhwAiO7taOorIKi0kqGpsQwpIt2TYnIf/lFcRhjzge2W2tXHLGvvjOQU+P+Ns+y/1McxphbgFsAkpKa1+6UEd3a8flvxjD5o9U8M2s97y/OIa+wrNYhS1I7tOS6kSlMHBjPNxl7eP7bTDJ2FdCzYyv+evkgTkqOcfAViEgg8dkxDmPMV0DHWh6aDPwWOMtae9AYkw2ke45xfAo8aa2d5/kcXwP3WWuXHOu9mvIxjmOx1vLR8u28szCH1A4tOSk5hr6d21BWUcWB4nI25xXx5oItrN2ZT3CQobLK0i0uittO7c7EgfGEaLIlkWbN745xWGvH1rbcGNMPSAEOb20kAEuNMUOp3sJIrPH0BMDtxQx+zBjDpEEJTBqUUOvjI7q14/KhiSzcvI9PV+1kWEo7xvftSHBQ452RJSJNj/NdVdbaVUD7w/eP2OKYDtxhjHmX6oPiB3V8o36MMQzr2o5hXdu5jiIiAcp5cRzHZ1SfiptJ9em417uNIyIiflcc1trkGh9b4HZ3aURE5Eg6GioiIl5RcYiIiFdUHCIi4hUVh4iIeEXFISIiXlFxiIiIV/xqWPWGYozJBbb48C1jgTwfvl9dBEJGUM6GFgg5AyEjNI+cXay1ccd7UpMsDl8zxiw+kfFdXAqEjKCcDS0QcgZCRlDOmrSrSkREvKLiEBERr6g4GsZLrgOcgEDICMrZ0AIhZyBkBOX8Dx3jEBERr2iLQ0REvKLi8IIx5jFjzEpjzHJjzBfGmHjP8lONMQc9y5cbYx6u8Zrxxpj1xphMY8wDjnMaY8zfPFlWGmMG13jNtcaYjZ7btT7K+YwxJsOTZZoxJtqzPNkYc6jG+nyxxmuGGGNWeb6Gv5kj5hr2VUbPYw96cqw3xoyrsdzF9/wSY8waY0yVMSa9xnK/WZfHyul5zG/WZ433/r0xZnuN9XfO8fK64tP1ZK3V7QRvQOsaH/8KeNHz8anAjFqeHwxsAroCYcAKoLfDnOcAMwEDDAd+9CyPAbI8/7b1fNzWBznPAkI8Hz8FPOX5OBlYfZTXLARGeL6GmcDZjjL29nw/w6mewXKT5/vt6nveC0gDZlM9Edrh5X6zLo+T06/WZ41cvwfuqWV5rXl9lauWPD5dT9ri8IK1Nr/G3SjgeAeIhgKZ1tosa20Z8C4wsbHyHXaMnBOB1221BUC0MaYTMA740lq7z1q7H/gSGO+DnF9Yays8dxdQPTXwUXmytrbW/mCrf1teBy5wlHEi8K61ttRau5nqycaG4u57vs5au/5En+9iXcIxc/rV+jwBR8vrik/Xk4rDS8aYx40xOcCVwMM1HhphjFlhjJlpjOnjWdYZyKnxnG2eZa5yHi2Ps5w13ED1X72HpRhjlhljvjPGjPEs6+zJdpivc9bM6M/r8kj+uC6P5M/r8w7PrspXjTFtPcv8IVdNPs3jdzMAumaM+QroWMtDk621H1trJwOTjTEPAncAjwBLqb5Uv9CzD/QjoAfVuwCO1CCnsdUx59HyOMvpec5koAJ4y/PYTiDJWrvXGDME+MhTxo2Ss44Zj5altj/GfLYua+HTdVmPnD5fn/9542PkBaYAj3ne8zHgz1T/AdFo66+OfJpHxXEEa+3YE3zq28CnwCM1dw1Zaz8zxrxgjImluvUTa7wmAdjhKucx8myj+jhNzeWz6x2S4+f0HIg/DzjDs8sEa20pUOr5eIkxZhOQ6slZc3dWg6zPumTk2N9b19/zmq/x6bqsa04crM/DTjSvMeZlYIbnbqP9bteRb/O4OpgTiDegR42Pfwl84Pm4I/+9JmYosJXqvwBCqD7QnMJ/D1j1cZjzXH56cHyhZ3kMsJnqA+NtPR/H+CDneGAtEHfE8jg8BxqpPti3/XAeYJEn++EDuuc4ytiHnx4czaL6AKWT73mNXLP56UFnv1mXx8npr+uzU42P76T6uMZR8/oqVy05fbqenHyRgXoDpgKrgZXAJ0Bnz/I7gDWeb9YCYGSN15wDbKD6jIfJjnMa4HlPllVH/OLeQPUBvkzgeh/lzKR6v+xyz+3w2V8X1VifS4EJNV6T7vnaNgF/x1PYvs7oeWyyJ8d6apyR5Oh7PonqvzpLgd3ALH9bl8fK6W/rs8Z7v+H5XVkJTD+iSGrN6+rmy/WkK8dFRMQrOqtKRES8ouIQERGvqDhERMQrKg4REfGKikNERLyi4hAREa+oOERExCsqDhEHPPM83OM6h0hdqDhERMQrKg4RHzHGTPbM0PYV1ZMZiQQkjY4r4gOeIc0vAwZR/Xu3FFjiNJRIHak4RHxjDDDNWlsMYIyZ7jiPSJ1pV5WI72hEUWkSVBwivjEHmGSMaWGMaQVMcB1IpK60q0rEB6y1S40x/0v1nB5bgLmOI4nUmebjEBERr2hXlYiIeEXFISIiXlFxiIiIV1QcIiLiFRWHiIh4RcUhIiJeUXGIiIhXVBwiIuKV/w/WGEKNS7pNRgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "observer_gcrs2 = keck.get_gcrs(obstime - np.linspace(0, 365, 100)*u.day)\n",
    "observer_icrs2 = observer_gcrs2.transform_to(ICRS)\n",
    "\n",
    "d_pos = (αcen_icrs.data.without_differentials() - observer_icrs2.data.without_differentials()).to_cartesian()\n",
    "pos_hat = d_pos/d_pos.norm()\n",
    "\n",
    "d_vel = αcen_icrs.velocity - observer_icrs2.velocity\n",
    "\n",
    "rv = np.sum(d_vel.d_xyz * pos_hat.xyz, axis=0)\n",
    "plt.plot((observer_gcrs2.obstime-observer_gcrs2.obstime[0]).to(u.day), rv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Consider M31 (a very distant object):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)\n",
       "    (10.6847, 41.269, 710.)\n",
       " (radial_velocity) in km / s\n",
       "    (-300.,)>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# these are all just rough estimates, but right orders of magnitude for M31\n",
    "m31 = ICRS(ra=10.6847*u.deg, dec=41.269*u.deg, \n",
    "              distance=710*u.kpc, radial_velocity=-300*u.km/u.s)\n",
    "m31"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more complete testing, lets also try representing in *Galactic* coordinates, which are a rotation of ICRS but should have the same RV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Galactic Coordinate: (l, b, distance) in (deg, deg, kpc)\n",
       "    (121.1743321, -21.57305874, 710.)\n",
       " (radial_velocity) in km / s\n",
       "    (-300.,)>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m31 = m31.transform_to(Galactic)\n",
    "m31"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<GCRS Coordinate (obstime=2018-12-13 09:00:00.000, obsgeoloc=(0., 0., 0.) m, obsgeovel=(0., 0., 0.) m / s): (ra, dec, distance) in (deg, deg, kpc)\n",
       "    (10.68714893, 41.27286941, 710.)\n",
       " (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)\n",
       "    (-114874.28172291, -11517.00802253, 1878.85362296)>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m31_in_gcrs = m31.transform_to(observer_gcrs)\n",
    "m31_in_gcrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<CartesianDifferential (d_x, d_y, d_z) in km / s\n",
       "    (96828071.72300434, -3.7518773e+08, -29132262.93832655)>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m31_in_gcrs.velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again this is horribly wrong.  But does the observer-icrs transform fix it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$-279.40557 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity -279.40557273 km / s>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m31_in_obs = m31.transform_to(observer_icrs)\n",
    "\n",
    "d_pos = (m31_in_obs.data.without_differentials() - observer_icrs.data.without_differentials()).to_cartesian()\n",
    "pos_hat = d_pos/d_pos.norm()\n",
    "\n",
    "d_vel = m31_in_obs.velocity - observer_icrs.velocity\n",
    "\n",
    "rv = np.sum(d_vel.d_xyz* pos_hat.xyz, axis=0)\n",
    "rv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aha!  That now looks right - it's M31 with an added component from the Earth's rotation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we transform the observer to ICRS (or one of the systems that aren’t FiniteDifference connected), and then transform the target to that, all is well.  \n",
    "\n",
    "The long-term problem is that this is not strictly right in a GR sense.  Probably it’s mostly-ignorable for practical use cases, since the earth and sun are effectively in the same place for any cosmological considerations.  But for careful work especially inside the Solar System that might not be sufficient."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (astro36)",
   "language": "python",
   "name": "astro36"
  },
  "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": 2
}