modeling_vel_distrib.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sys\n",
    "sys.path.insert(0, '/astro/users/ctslater/.local/lib/python3.6/site-packages')\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.coordinates\n",
    "import astropy.units as u\n",
    "\n",
    "#plt.style.use((\"paper_plots.mplstyle\"))\n",
    "\n",
    "matplotlib.rcParams['figure.dpi'] = 80\n",
    "\n",
    "from astropy.coordinates import Galactocentric, Galactic, ICRS\n",
    "\n",
    "import healpy as hp\n",
    "\n",
    "cmap = matplotlib.cm.coolwarm.set_under('w')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ra, dec = hp.pix2ang(hp.order2nside(7), np.arange(hp.nside2npix(hp.order2nside(7))),\n",
    "#                                    lonlat=True, nest=True)\n",
    "ra, dec = 20, 30\n",
    "\n",
    "ncoords = 1000\n",
    "\n",
    "equatorial_coords = astropy.coordinates.SkyCoord(ra*u.deg, dec*u.deg, 10*u.kpc, frame='icrs')\n",
    "coords = equatorial_coords #equatorial_coords.galactic\n",
    "\n",
    "velocities = np.random.multivariate_normal([0, 0, 0], np.diag([80, 80, 80]), ncoords)\n",
    "\n",
    "gc_coord = Galactocentric(x=coords.galactocentric.x + np.zeros(velocities.shape[0]),\n",
    "                          y=coords.galactocentric.y + np.zeros(velocities.shape[0]),\n",
    "                          z=coords.galactocentric.z + np.zeros(velocities.shape[0]),\n",
    "                          v_x = velocities[:,0] * u.km/u.s,\n",
    "                          v_y = velocities[:,1] * u.km/u.s,\n",
    "                          v_z = velocities[:,2] * u.km/u.s,\n",
    "                          # galcen_v_sun=coord.CartesianDifferential([0,220,0]*u.km/u.s)\n",
    "                          )\n",
    "#print(gc_coord[0])\n",
    "galactic_coords_w_pm = gc_coord.transform_to(Galactic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAusAAAGbCAYAAAB9M95hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAQ6gAAEOoBgtMKmAAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XnUJFV9//H3F0ZAhIFBRRDQUfwpgyKDAsYYBeSICFFcQDRGMUeN0Z9xgdEMLqBodECDe0RFA54oIu4IIoosKrIo+GMbSASGABkWZcYBhGH7/v6oaqdoup95ll7u0/1+ndOnp6tu1XO7pm/3p27dqorMRJIkSVJ51hl2BSRJkiR1ZliXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViX1FVEHBAR34+I6yPizoi4JCLeEhHrNMocFxHZ4bF327oeFhEfi4jlEfHniDgzIp4++HclSdLsMWfYFZBUtEOA64B3AzcDewCfAZ5YT2u5BnhN27JL215/Enhdvc5lwHuAMyJih8y8qec1lyRpBERmDrsOkgoVEY/OzFvbph0NvAXYNDNXR8RxwM6Z+bQJ1rMVVeh/e2b+ez1tY+Ba4NjMXNyv9yBJ0mzmMBhJXbUH9drFwAbAZlNY1V7AusA3G+u+HTgZ2HcmdZQkaZQZ1iVN1XOB24BbGtO2jYiVEXFPRPw2Il7atswC4ObMvK1t+hXAU5pj4CVJ0hojM2Y9ItYHdgBuBe4fcnWkmVgXeDRwaWauHnZlmiJiZ+AfgA9lZqudXQxcCFwObEo1ROZ7EXFAZn67LjMPWNlhlSuAhwEbAau6/M25wNzGpPWAbajGvdvWNdsV296Hzd91jZhpt/WRCetUDfrCYVdC6qFdgN8MuxItEbEF8B3gAuDI1vTM/HRbuR8C5wJHAN9uzOp0gkxMMK/lYODwaVRZmk2Kau+F8Hddo2jKbX2UwvqtABdccAFbbrnlsOsiTdvy5cvZddddof5MlyAiNgF+DPwZeElm3tutbGY+EBHfAY6KiIdn5l1UPejzOhTfFLgXuHOCP380cGzj9VbAebZ1jYIS23tB/F3XyJhJWx+lsH4/wJZbbsnWW2897LpIvVDEYd+I2AD4IfAY4NmZ+cfJLNb2eimweURs1jZufXvgqsx8oNuKMnMVjSEyEdWqbesaMUW098L4u65RNOW27kldkrqKiDnAt4Adgb0z87pJLLMOsD9wed2rDnA68ADwyka5jYAXA6f0ut6SJI2KUepZl9R7n6cK1O8BNoyIv2rMu4JqaMtxwAnA1fXrtwA7A69oFczMGyPiGODIiLiP6prri+rZn+rze5AkadYyrEuayAvr56M6zNsDuIRqiMphVGe530N14syLMvMnbeUPBu4APgJsApwP7OndSyVJ6s6wLqmrzJw/iWL7TXJd9wCL64ckSZoEx6xLkiRJhTKsS5IkSYUyrEuSJEmFMqxLkiRJhTKsS5IkSYUyrEuSJEmFMqxLkiRJhTKsS5IkSYUyrEuSJEmF8g6mmrL5i0+Z0fLLluzbo5pImohtVdJkzPS7Avy+6Cd71iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCzRl2BSRJkjR98xefMuwqqI/sWZckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgrlpRslSR3N9HJwy5bs26OaSNL4smddkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZKkMRcRB0TE9yPi+oi4MyIuiYi3RMQ6beX2iYiLI+LuiPh9RLy1y/oWRcSyutyFEbH7QN6INIIM65Ik6RBgNfBu4G+B7wOfAY5sFYiIZwM/AC4CXgQcB3w2It7YXFFELAI+CnwO2Af4PfDjiNih7+9CGkFzhl0BSZI0dC/OzFsbr8+MiI2At0XE+zNzNXAYcFFmvqFR5nHAERHx1cx8ICLWB94PfCozPwEQEWcDlwLvA141sHckjQh71iVJGnNtQb3lYmADYLM6hD8f+GZbma8DWwI71a//GtgEOKGx7vuBE4F9IiJ6XHVp5BnWJUlSJ88FbgNuAbYF1gOWtpW5on5e0PZ8ZYdyGwNb9b6a0mhzGIwkSXqQiNgZ+AfgQ5l5f0TMq2etbCu6on7erH6eB6zOzLsmKHdDl785F5jbmLTFdOoujZop96x7xrgkSaMrIrYAvgNcQOME01p2WSzXUiYmmNdyMHB943HhWisrjYHpDIPxjHFJkkZQRGwC/Bj4M/CSzLy3ntXqGZ/Xtsi8tvkrgA0iYoO2cpu2levkaGCbxmOXqdVeGk3TGQbjGeOSJI2YOmD/EHgM8OzM/GNj9tXAPVRj0k9rTN++fl7a9ryA6gTVZrnbgRu7/f3MXAWsatRn6m9CGkFT7ln3jHFpPDjkTRofETEH+BawI7B3Zl7XnF93xP0ceGXboq8GlrMmmJ8L/Ak4sLHudevlTs3MiYbBSOqgVyeYNs8YfwprP2P8t0zujPGOJ6FIGohDgOuohrzdDOxBNeTtifW05pC3r1GNN30O1ZC3ezLz2NaKGkPe3ks1PO5NVEPeds3MSwf2jiR183ngxcB7gA0j4q8a866oe72PAM6JiC9TdcA9h6otvzkzH4Aq1EfER4CPRsStVO39jVTfGx4xl6ZhxmF9GGeM13/Xs8al/nLImzQ+Xlg/H9Vh3h7AWZn564jYj2rH+3VUv9Fvb+6Y1/6N6oTSt1MNqbkU2Mcdc2l6ZhTWh3jGOFS9eIevrY6SpmcSQ95uoxrytritzNepett2ojqK1nHIW0ScCBwSEeGhcWm4MnP+JMudCpy6ljIJfLx+SJqhad8UachnjINnjUvD4E1SJEkaoGn1rA/7jHHwrHFp0BzyJknS4E05rLedMf68TmeMR0TrjPFPNmZNdMb4xfW6PWO8z+YvPmXYVdAs5JA3SZKGYzo9654xLo2RXg55y8y7G+WmMuSteQLbFnhnQ0nSmJhOWPeMcWlMOORNkqThmnJY94xxaTw45E2SpOHr1U2RJI0eh7xJkjRkhnVJ3TjkTZKkITOsS+rIIW+SJA3ftG+KJEmSJKm/DOuSJElSoQzrkiRJUqEM65IkSVKhDOuSJElSoQzrkiRJUqEM65IkSVKhDOuSJElSoQzrkiRJUqEM65IkSVKhDOuSJElSoQzrkiRJUqEM65IkSVKhDOuSJElSoQzrkiRJUqEM65IkSVKhDOuSJElSoQzrkiRJUqEM65IkSVKh5gy7ApKkh5q/+JRhV0GSVADDugZupiFk2ZJ9e1QTSZKksjkMRpIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpSXbpQkSRqiUbivgpdl7h971iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSRpzEfGkiDgmIn4XEfdFxGUdypwVEdnhsV1buY0j4osR8ceIuCMifhgRjx/cu5FGy5xhV0CSJA3dU4F9gfOpOvK6deb9CljUNm1Z2+sTgGcAbwNWAUcAP4uIp2fmXb2qsDQuDOuSJOnkzPwBQEQcB+zcpdzKzDyv20oi4llUoX/fzDy1nnYpcDVwEHBMLystjQOHwUiSNOYy84EerWofYCXw48a6/wf4JVWIlzRF0wrrjm2TJGks7RYRd0bE3RFxdkQ8r23+AuCqzMy26VfU8yRN0XSHwTi2TZKk8XI28DXgv4HHUv2+/ywidsvMX9dl5lH1rLdbAWw20cojYi4wtzFpixnXWBoB0w3rjm2TJGmMZObhzdcR8SPgcuADVMNf/lK0w+LRZXrTwcDhaykjjZ1pDYNxbJs0HhzyJqmbzLwTOAV4ZmPyCqre9Xab1vMmcjSwTeOxSw+qKc16/b4azG4RcSewLtWQmQ9k5jmN+RONbXthn+smae0c8iZpItH2einwgoiItt/27et5XWXmKqrvhmrF0b5qaTz182owZwPvAPamGtKyIdUP87MbZWY0ti0itm49cGyb1A8nZ+Y2mbk/cNEE5VZm5nltj7tbMxtD3t6YmSdk5inAy4D5VN8PkmaZiHgEVbu+sDH5VKpe9Bc2ym0D/A1VL7ykKepbz7pj26TZr99D3iKiNeTN81OkIYqIDVnz2/x4YG5E7F+/PhvYjuro2feA66hOMD2EqqPsgNZ6MvP8iDgF+EpEHMKao2jLgOP7/06k0TOwmyJl5p11A96/MXkF8LgOxSc7tu3YxustePDevaTBccibNLttDpzUNq31eg/gBmB94GPAI4E7gXOBf8rMC9qW+zvgE8C/A+sBPwde4XA3aXoGfQdTx7ZJo8fLuUmzXGYu46G/0e32nuS6VgH/WD8kzdDA7mDq2DZpNGXm4Zn51cz8RWaeCOwO/C/VkLcHFe2w+GSHvF3feHgETZI0NqbVs+7YNkndOORNkqTeme4wGMe2SZqIQ94kSeqBaYV1x7ZJ6maCIW+HUQ15O60u1xry9s+DrqMkSbPFoE8wlTSLOORNkqThMqxLmohD3iRJGiLDuqSuHPImSdJwDezSjZIkSZKmxrAuSZIkFcqwLkmSJBXKsC5JkiQVyrAuSZIkFcqwLkmSJBXKsC5JkiQVyrAuSZIkFcqwLkmSJBXKsC5JkiQVas6wKyBJo2j+4lOGXQVJ0giwZ12SJEkqlGFdkiRJKpRhXZIkSSqUY9ZnGcfBSpIkjQ971iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQnlTJM06vbgx1LIl+/agJpIkSf1lz7okSZJUKMO6JEmSVCjDuiRJklQow7okSZJUKMO6JEmSVCjDuiRJklQow7okSZJUKMO6JEmSVCjDuiRJklQow7okSZJUKMO6JEmSVCjDuiRJklQow7okSZJUKMO6JEmSVCjDuiRJklQow7okSZJUKMO6JEmSVCjDuiRJklQow7okSZJUKMO6JEmSVCjDuiRJklQow7okSZJUKMO6JEmSVCjDuiRJklQow7okSZJUKMO6JEljLiKeFBHHRMTvIuK+iLisS7mDIuLKiLg7Ii6LiAM6lHlYRHwsIpZHxJ8j4syIeHr/34U0mgzrkiTpqcC+wO+BKzoViIj9geOA7wEvAs4AToyIvdqKfhL4v8BhwH7AfcAZEbFFX2oujbg5w66AJGk0zV98yozXsWzJvj2oiSbh5Mz8AUBEHAfs3KHMh4GTMvPQ+vWZEbEdcARwer3sVsA/AW/PzC/X084DrgXeCSzu55uQRtG0etY9XCZJ0ujIzAcmmh8RTwC2A05om/UNYNeIeFT9ei9gXeCbjXXfDpxM1XMvaYqmOwzGw2XSGHDHXFJtQf28tG36FUBQBflWuZsz87YO5Z4SEQ6/laZouo3m5MzcJjP3By7qUuYvh8sy88zMfAfwU6rDZcCDDpctzswvZ+ZPgZdTNfx3TrNuknrHHXNJAPPq55Vt01fUz5s1yrWXaZV7GLBRtz8QEXMjYuvWA/C7QWKaYd3DZdLYcMdcUlO2vY4O09vLdCvX7mDg+sbjwulUUBo1/Toc5eEyaQS4Yy6p1upBn9c2fdO2+Ss6lGmVuxe4c4K/cTSwTeOxy7RqKo2YfoVhD5dJ48Edc2k8tNr4grbp21P1ll/ZKLd5RGzWodxVE3UAZOaqzLyh9QBu6kG9pVmv3z+QHi6TRps75tIYyMxrqQL5gW2zXg1ckJl/qF+fDjwAvLJVICI2Al4MzPxantIY6td11puHy25uTO/14bJjG6+3wMAuDUu/d8wPn2a9JE1CRGwI7FO/fDwwtz55HODszLyV6uTwEyPiaqrzUvajGuK2d2s9mXljRBwDHBkR9wHXAYvq2Z/q/zuRRk+/wnrzcNmVjeldD5e1HR6f1OEyYFXrdUR0Kyqpf9wxl0bD5sBJbdNar/cAzsrMk+pQ/16qAP574MDMPL1tuYOBO4CPAJsA5wN7ZqbDWqRp6EtYz8xrI6J1uOx7jVkTHS47Bh50uKz54yypTO6YSyMgM5ex5kjXROWOB45fS5l7qO5U6t1KpR6YVlj3cJkkcMdckqR+m27PuofLpDHgjrkkScM1rbDu4TJpbLhjLknSEPXrBFNJI8Adc0mShsuwLkmSNE3zF3v5ePWXdw2UJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKNWfYFZAkSdJ4m7/4lBktv2zJvj2qSXnsWZckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIKZViXJEmSCmVYlyRJkxIRr4+I7PBY0lbuoIi4MiLujojLIuKAYdVZmu3m9GvFEfF64D86zDoyMxc3yh0EHArMB34PfCgzT+pXvSSA+YtPmdHyy5bs26OaSNKstDfwp8brG1v/iIj9geOAJcDpwEuBEyPiT5l5+iArKY2CvoX1Bhu0NMLcMZfG0m8z8w9d5n0YOCkzD61fnxkR2wFHUP3WS5qCQYR1G7Q0Htwxl8ZcRDwB2A54b9usbwD/ERGPmiATSOpgEGG9Ixu0NHLcMZfGx+UR8SjgOuDLwFGZeT+woJ6/tK38FUBQ/e7/stMKI2IuMLcxaYue1liapQZxgunlEXF/RFwTEYdGxLr19Mk0aEmzXGPH/IS2Wd8Adq1/8CXNDsuBw4HXAS8CTgU+Any6nj+vfl7ZttyK+nmzCdZ9MHB943FhD+orzXr97FlvNejzgQReQtWgtwLexswatHvgUnl63tMmqSyZ+RPgJ41Jp0fEXcC7IuJfm0XbFo0u05uOBo5tvN4CA7vUv7De5wYN1R744TOrpaQecMdcGm/fAhYBC1nTrucBNzfKbFo/r6CLzFwFrGq9johuRaWxMugx6z1p0LVZuQc+00sGSqVxx1wae81U3TqCtgC4sjF9e6q23pwmaRIGfVOkbg26aVINOjNXZeYNrQdwU++qKWmGvgWsy0N3zJumsmO+TeOxS4/qKKk3DgTuBy7OzGupfr8PbCvzauACLxwhTd2ge9abDfqmiGg16O81ytigpdmvZz1tHhqXyhERPwHOAC6rJ70E+Efg05nZ6jQ7jOrSrFcDPwX2A/aiuryrpCnq5x1MbdDS+HLHXBpNVwJvBLamOjr/X8A7gc+2CmTmSRGxIdWlmRdR3QTtQO+pIE1PP3vWbdDSGHDHXP000/N8li3Zt0c1EUBmvgN4xyTKHQ8c3/8aSaOvn1eDsUFL48Edc0mS+mRodzCVNBrcMZckqX8GfTUYSZIkSZNkWJckSZIKZViXJEmSCmVYlyRJkgplWJckSZIK5dVgJKmDmV7fW5KkXrBnXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSrUnGFXQJIkaVjmLz5l2FWQJmTPuiRJklQow7okSZJUKMO6JEmSVCjDuiRJklQow7okSZJUKK8GI0mSpFltplf1WbZk3x7VpPfsWZckSZIKZc+6NA2jvAcvSZLKYc+6JEmSVCh71iVJI8ujYJJmO3vWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCGdYlSZKkQhnWJUmSpEIZ1iVJkqRCeVOkKZjpzTUkDYZtVZI0KgzrkiRJGmu96OTp1x2PHQYjSZIkFcqwLkmSJBXKYTCSJGnW8hwVjTp71iVJkqRCFRHWI+LJEXFaRNwZEbdExKcj4uHDrpek3rKtS+PBti71ztCHwUTEpsDPgeuAVwCbA0cDjwT+fohVk9RDtnVpPNjWpd4aelgH3gzMAxZm5h8AIuI+4OsR8a+ZuXSotZPUK7Z1zTolX86tYLZ1qYdKGAazD/CzVoOufQdYXc+TNBps69J4sK1LPVRCz/oC4KvNCZm5OiKuruf1jGeMS0M1sLYuaahs61IPlRDW5wErO0xfAWzWbaGImAvMbUzaCmD58uVd/9B9q/7QdZ40SDfccEPXeY3P8LoDqczg2NY1lsawvQ+srYPtXeXoV1svIawDZIdp0WV6y8HA4e0Td911117VSeqbbb4wqWLzqU7QGiW2dY2dSbb3RzNa7d22rrHTr7ZeQlhfQbUX3m5TYKKTUI4Gjm28Xg/YBlgG3N+ryk3TFsCFwC7ATUOuS4ncPhPbCjgPuH7YFemxUWrrfoa7c9tMrH37rEv1433pMCvVY6PU1gfBNjOxUdk+027rJYT1pbSNYYuI9YFtaRvz1pSZq4BVbZOv6XntpiEiWv+8KTO7HxMZU26fiTW2zz3DrEcfjExb9zPcndtmYl22zyj1qMMItfVBsM1MbMS2z7TaeglXgzkV2DMiHtmY9jJg/XqepNFgW5fGg21d6qESwvoXqU5E+UFEvDAiXgt8Fvi612KVRoptXRoPtnWph4Ye1jNzJfB84E7gu1Rj1k4A3jTMes3QKuBDPPRwnipun4mN5PYZsbY+kv9HPeK2mdjIb58Ra+uDMPKfiRka++0TmROdmC1JkiRpWIbesy5JkiSpM8O6JEmSVCjDuiRJklQow7okSZJUKMP6DEXEuhHxnog4OyJujYgVEXFOROw5yeWzw2M236HrQXqwfR4WER+LiOUR8eeIODMint7veg9SRLwgIr4REVfX//+fm+Ry87t8fs7rd51HVX2ZudZndXVEXBMRR0fEJpNY9qCIuDIi7o6IyyLigEHUeVCmu20i4qwun9PtBlX3QYuIjSLihvp97jyJ8osiYln92bkwInYfQDU1YJP9rp9sLoiIJ0fEaRFxZ0TcEhGfjoiH9/+d9Mdktk/9nj8bEVfU7/u6iPhKRGzRVm73Ltvxm4N7R71Twh1MZ7uHA+8Fjgc+DtwLvB74aUS8JDN/NIl1fBb4RuP1KN25cqbb55PA64BDqG45/R7gjIjYITNHZafmRcBC4Gxgs2ks/17gzMbr23tRqTG1GXAu8CmqW6Y/Dfhg/bxXt4UiYn/gOGAJcDrwUuDEiPhTZp7e3yoPzLS2Te1XwKK2act6W72ifIBJ/r5GxCLgo1Tt+CKqyxv+OCJ2zcwp35ZcRZvKd/2EuSAiNgV+TnVHzFcAm1NdIvORwN/3qL6DNpntsxewG/Al4HfA1lTfQ7+uc8EdbeX/Abiy8foPvazwwGSmjxk8gHWBeW3TAvgtcOYklk9g0bDfR4nbB9gKuA94a2PaxlSNbcmw31sPt9E6jX8vAz43yeXm15+f/Yf9Hkb5QRWeEnjsBGWWAt9qm/YT4Lxh17+AbXMW8KNh13WA22Q74A7gzfW22XmCsutT3TzoqMa0dYErgG8O+7346PlnY1Lf9ZPJBcC/UF3H/lGNaX9XL7tg2O+1X9sHeBT1Zccb055ev++DGtN2X1v7m00Ph8HMUGben5kr2qYl1R7fY4dTq3LMcPvsRfXD9ZfDVpl5O3AysG+Pqzo0mfnAsOugCf2xfn5Yp5kR8QSqgHZC26xvALtGxKP6WLdhm3DbjKnPAMcAV02i7F8Dm9D47GTm/cCJwD4REX2poYaix9/1+wA/y8xmT/F3gNX1vFlnMtsnM/9QZ4imS4H7GeHMZVjvg4hYh+pLeLK3VV4cEfdGxMqIODEiHtfH6g3dFLbPAuDmzLytbfoVwFPq9Qi+EBH312MWvxwR0xlKo4b6XIsNIuIZwGHAyZl5XZfiC+rn9s/zFVRHkUZqbPYUt03LbvX40rvrce/PG0BVB64eDrUjcMQkF2l9dq5sm34F1VHErXpUNc0+a8sFC2j7zsnM1cDVrPlcjYtnU3XsdcoUp9a/jzdExMdn65h+x6z3xz8DT6E6DLo2XwN+BNxMNfbzA8AvI2LH9h7pETLZ7TOP6hBxuxVUPXkbMca3H6bqQfkC1XCLlcCzgPcBO9fjXe8dZuVmuetYE5ROA149Qdl59XP7Z7XVfkdt52kq2waq8adfA/6bqudrEfCziNgtM3/dt1oOWERsSDVm+NDMXDXJTvF5wOrMvKttevOzc0PvaqlZYjK5YKLfx1H7zukqIh5GdR7NVVTbrOVPwFHAOcBdwPOpvnsWAH874GrOmGG9g/rqBltOoui19Z5sc9ndqD4gn8jMc9a2gsw8qPHynIj4JWtOMjpq8rUenEFuH6oxZw+pwgTzhm4m22cqMnM58NbGpLMj4nKqL6yXAd+a7rrFPlQ7g0+l+qE8OSJeUA9R6Kb981j053QGprRtMvPw5uuI+BFweb3srDxc38X7qcLVcVNcbtZ9x6nSr+/6KeSCbp+dIj43A/ot/BzVDs3zMvO+1sTMvBi4uFHu5xGxHPhc3ZkXr0FqAAAK6klEQVR1wTT/3lAY1jt7GfAfkyi3E9XYawCiuqTgD4DvU538MWWZeUlEXAU8czrLD8igts8K1vRaNm1KdVWZOyexjmGY1vbpkVOpTm57Job1acvMS+p/nhsRFwG/ofp//XaH4s2erpsb0zdtmz8SprhtOi1/Z0ScAuzfpyoOXEQ8nuqKVS8D5ta96hvVszeKiI3yoVepgOqzsUFEbJCZdzemj+RnZwQN5Lu+Sy6Y6PdxskNw+62v2yciDgfeALw8M38ziUW+RRXunwnMqrDumN8OMvO4zIxJPJpBdFuq4QgXAa/tcALEVBR9UtEAt89SYPMOY7C3B64q9cTM6WyfHiv68zML/Y7q5KUndZnf+mFsHye6PVUPV/t45FGytm3Tzah9Rp8ArAecQhWiVlCdCA/VZVV/1mW5iT47twM39raa6qUBf9e3t5mltH1uImJ9YFsKCev93D4R8VaqSza+NTN/2PPKF8aw3gP1xfhPB24CXpqZ075OekQsBJ4MXNij6g3dDLbP6cADwCsb69oIeDHVj6Ie6m+BRzBCn58CtE5euqbTzMy8liqQH9g269XABW1Xaxg1E26bTiLiEVRXcxqlz+jvgD3aHu+q5/0TDx6u1nQu1djav3x2ImJdqu+8U2fY6aMR0SUXnArsGRGPbEx7GdXlQE8dYPUGLiJeRXUd+sMy80tTWPRV9fOs++5xGMwM1WcWn0Z1Q4KDge2bJxZl5nmNsr8HrsvMPevXi4AnUp2AdQvVuKv3AdcDxw7oLfTVTLZPZt4YEccAR0bEfVQntrVurPKpwbyD/qsPoe9Sv9wQ2La+qgSZ+e1GufuA4zPzDfXrT1DtzJxPdaLRrsChVMMSvj+wNzBCIuK7VNvvEqqTknakuhHXJdTbNCK+QnU93+b352FUN0G6GvgpsB/VpUf3Hlzt+2s62yYinkvVZr9H1X4fSzVcZAtgZO7wmpkrqa4n/xeN77nfZuZF9bQzgMdn5pPq5VZHxEeAj0bErVRHHt9I9bvwKjRSJvNdP4Vc8EWqizX8ICI+zJqbIn09M4voWZ+qSW6f3ahOwP0F1c0V/6qxilsz8+q63H8Cv6dqU3dTnWD6LuAHkxwyU5Ys4GLvs/nBmhvTdHy0lV0GnNV4/WLg18BtVGOw/xf4CrDlsN9XCdunnrYe1V0hb6IKCGcBOw77ffV4G71+ktsngeMar99AdXOpP9Wfn2VUd3ydO+z3NFsfwGKqk5JWUY39v4zqMnxzG2WOa/+/qacfRHVFgtVUJ1AeMOz3M+xtQzU85jRgOdUdGFdQHRXbddjvZwDba3fabspSf38taysXwLupdmbupur122PY9ffRl8/EWr/rp5ILqHrbf0J1/tatVNf4f/iw32eft88Hu5Vp+308tP6Our3+7rkKOBxYb9jvczqPqN+UJEmSpMI4Zl2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdHUXExhHxQES8om362yLi9oiIYdVNUpki4sKIeHvbtBdFREbE5l2WeX9E/HQwNZTGW0RcGhHHdZj+voi4JSI2HkK1tBaGdXWzEAjg/3WYfklm5uCrJKlUEfFy4PHAl9tmLQRuzMxbuiz6OeBZEfH8ftZPEgCXAQuaEyJiHrAI+Ghm3j6UWmlChnV1sxC4Hbi6w/T2AC9J7wS+kZl3tU1fCFzUbaHMXAl8D3hHH+smqfKQsA68h+r3/guDr44mw7Cubh7Sgx4R6wJPpQ7rEfHEiLh8SPWTxsp0h6ZFxDp1md9FxF0RcWtE/DAiNqvnrx8RH4qIayLizxHxm0693BGxZ0ScGRF/jIg76sPpr63nPRF4LvDtDlXYCbi0/hs3RsRtEXFsRGzYKHMSsE9EPHp6W0fSJF0ObBwRWwNExGOAfwY+mJmrh1ozdWVYVzedetC3AzZoTH8GcOEgKyWNsSkPTasD/H8CS6h6r/cD3g6sB9wREesDpwOvBY4AXkL1Y35qRDylsZ6/AU4FzgcOBA6gCtj31UX2BO6l7fsgIh4BbAscBGxW/50l9evDG0V/BcwBdp/ktpA0PZfVz63e9fcBNwDHdypcd9JpyOYMuwIqT0TMoepB/2LbrOcADwCX1q+fCayOiDOB+cB3MnPRoOopjZmJhqZd0GWZtwAvB56Tmb9tTD8BICKOAJ4GPD0zb6ynnQXsDbwGOKwu/xrgrMxc3FjHjxv/3hn4rw49cztSdQqdmJmH1NN+HhE71PX6F4DMXBER/wM8i2onQFJ/XAPcBSyIiKuANwOvycz7WwUi4vNUHQMLgAQ8n2TI7FlXJwuA9VnTa9YK8G8Crs7MO+vJzwA2BV4A7AC8KiLax8JJ6o3JDE17ckRcUv87gPcCX2wL6s1l3wp8oRXUATLzPuB6YKtG8T8Du0XE4ojYpkPdtgRu7TB9J+Buqt70psuA9iEvfwC26LAOST2SmQ8AS6l+5z9I1fn2nbZiC4F5wAsz06BeAHvW1clCqh70xRFxd/3vN1IF8usj4jmZ+SuqXrMd6x/3OyLiv4HNqb4IJPXWQuC8tmmdhqa1gvn2VIH7u13WtwPwSKDTZRO35sE954dT9bQdAnw0In4F/EtmnlvP3wDoNN51IfDbzGwP8o+hOvTedDfw8C51ldQ7lwN7UV29aZ+2DoB1qI62vSwz7xlS/dTGnnV1shD4L+BLwL8DnwbOAI6i6vnaISIeB9yVmTdDdZIa1Z76pR3XKGnaGkPT2serdxqa1rryypb18/Iuq231bN/U9rcWUIXpc1rTMvOOzDyYamd8d2ATHtwbdxvVUbZ2C4Hr2ta/LtXY+fadhHnAH7vUVVLvXAY8EfhFZp7eNu//ANdNcKlVDYFhXZ0sBC7OzE9k5qaZ+ZjM/FhmHpaZG2fmMVSh4NH1meRQnZx2YmbeNrRaS6NrKkPTWj3rrRD+1C7r/J/6edu26R8BrgJ+3r5AVs6h6q1vHpm9CnhCs2wdyp/WYf2vBbahur56q+w6wOPq9Ujqo8w8KjMjM/foMHshcPGg66SJOQxGnewIHLmWMs8APg98t74E3OnAu/tdMWlMTWVoWqv3fSnwO+AzEbEpsIwqOD8rM99EdfTsfOCTEbEB1bj0N1EdHn9e64SziPgScD9wJnALsAtVW29+R/wKOCwits7M1vCWBVTDYx4TEUcCpwF/DbwfeE9mNk+U3R54BPCLmW4oSTOyE4b14hjW9SD1tVcfyVoaa2Z+YDA1ksSaoWlfoRqathr4FHAu8C6qoWnLgZtbveyZeX9E7Ad8HPgYMBe4tl6ezMz6mu2fBY6t/87ZwC6ZeWXjby8FXlU/5tT1eGtmNi/1dhbVCaIvYs0dTBdSHQl4AfBVqpseXQO8OTO/1vb+9qEaLuOlYKUharvikwoR3jVeksoWEWdQBfG/m6DMAcBLMvO1g6vZg/7+vwE7TefqERFxEfD9zDyi9zWTpNnNMeuSVL4dWfuh6Wew5uTSYfg48KyI2GkqC0XEblT3afhMPyolSbOdw2AkqWBTGJp26GBq1PXv3xQRr+eh109fm7nA6zJzZe9rJUmzn8NgJEmSpEI5DEaSJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSqUYV2SJEkqlGFdkiRJKpRhXZIkSSrU/wecihhCtt6stQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 880x440 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(dpi=110, figsize=(8,4))\n",
    "\n",
    "#def plot_common():\n",
    "    \n",
    "plt.subplot(1, 3, 1)\n",
    "plt.hist(galactic_coords_w_pm.pm_b.to_value())\n",
    "plt.xlabel(r\"$\\mu_b$\")\n",
    "\n",
    "plt.subplot(1, 3, 2)\n",
    "plt.hist(galactic_coords_w_pm.pm_l_cosb.to_value())\n",
    "plt.xlabel(r\"$\\mu_l cos(b)$\")\n",
    "\n",
    "plt.subplot(1, 3, 3)\n",
    "plt.hist(galactic_coords_w_pm.radial_velocity.to_value())\n",
    "plt.xlabel(r\"$v_r$\")\n",
    "\n",
    "plt.subplots_adjust(wspace=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'$\\\\mu_l cos(b)$')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAFzCAYAAADG/Ph0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcZWV95/HPV5YWoavUuDQqKCJBxWATFsUoGBUSycTEZMZ9IsRoXBKNiYlpnZfgTJSYEDQDqAkhw/KK0WASNDoizkTcQJTFBQTHGHbo1hZTRQQahN/8cW/pTXVtd6m6VU9/3q/XeVWd5z7POb++NvLlOc85J1WFJEnSWne/cRcgSZI0CoYaSZLUBEONJElqgqFGkiQ1wVAjSZKaYKiRJElNMNRIkqQmGGokSVITDDWSJKkJhhpJktSENR1qkhyR5J+S3JKkkvzyrM+T5ITu53cmuTDJAUs47muTXJvkriSXJXnG8v0pJEnSKKzpUAPsDnwV+K15Pv8D4He7nx8KbAY+lWT9fAdM8kLgPcA7gIOAzwGfSLL3COuWJEkjllZeaJmkgOdX1Xnd/QC3AO+pqnd129YBW4A3V9VfzHOcS4DLq+o1PW1XA+dV1aZl/mNIkqQB7TzuApbRPsAG4IKZhqraluQzwNOA7UJNkl2Bg4E/nvXRBd0xc+qGpXWzmh8M3DZQ5ZIk7djWA7dUnzMvLYeaDd2fW2a1bwEePc+YhwA7zTNmw/bdf2QTcHy/BUqSpHk9Cri5nwEth5oZs1Ne5mgbdsyJwMk9++uBm2688UYmJiaWVKQkSYLp6Wn22msvgNv7HdtyqNnc/bkBuLWn/WFsPxMzYytwL9vPyiw0hqraBmyb2e8s54GJiQlDjSRJK2St3/20kGvpBJujZhq6a2aOBC6aa0BV3Q1c1jum66j5xkiSpNVhTc/UJNkDeFxP0z5JNgK3VdUNSd4DvCXJt4BvAW8B7gA+0HOM/wv8Y1Wd2m06GTgnyaXAxcCrgL2B9y/7H0iSJA1sTYca4BDg0z37M+tazgKOBf4E2A14L/Ag4BLg6KrqvU63L50FwgBU1YeS/ATwNmBP4ErgmKq6fpn+DJIkaQSaeU7NapJkApiamppyTY0kSX2Ynp5mcnISYLKqpvsZ2/KaGkmStAMx1EiSpCYYaiRJUhMMNZIkqQmGGkmS1ARDjSRJaoKhRpIkNcFQI0mSmmCokSRJTTDUSJKkJhhqJElSEww1kiSpCYYaSZLUBEONJElqgqFGkiQ1wVAjSZKaYKiRJElNMNRIkqQmGGokSVITDDWSJKkJhhpJktQEQ40kSWqCoUaSJDXBUCNJkppgqJEkSU0w1EiSpCYYaiRJUhOaDjVJrktSc2ynzdP/2Hn633+la5ckSf3ZedwFLLNDgZ169p8EfAo4d4Ex08D+vQ1VddfoS5MkSaPUdKipqu/27if5Q+DbwGcWHlab+zlPknXAup6m9f2MlyRJw2v68lOvJLsCLwP+uqpqga57JLk+yU1JPpbkoCUcfhMw1bPdNHzFkiSpHztMqAF+GXggcOYCfa4BjgWeB7wYuAv4QpL9Fjn2icBkz/aoIWuVJEl9ysKTFu1I8kng7qr6xT7G3A+4HPhsVb2+j3ETwNTU1BQTExP9FytJ0g5qenqayclJgMmqmu5nbNNramYkeTTwHOBX+hlXVfcl+TKw2EyNJEkasx3l8tNxwHeAj/czKEmAjcCty1GUJEkaneZnarqXkI4DzqqqH8767Gzg5qra1N0/Hvgi8C1gAng9nVDzuhUtWpIk9a35UEPnstPewF/P8dnewH09+w8E/hLYQOcupiuAI6rqS8tdpCRJGs4Os1B4JblQWJKkwQyzUHhHWVMjSZIaZ6iRJElNMNRIkqQmGGokSVITDDWSJKkJhhpJktQEQ40kSWqCoUaSJDXBUCNJkppgqJEkSU0w1EiSpCYYaiRJUhMMNZIkqQmGGkmS1ARDjSRJaoKhRpIkNcFQI0mSmmCokSRJTTDUSJKkJhhqJElSEww1kiSpCYYaSZLUBEONJElqgqFGkiQ1wVAjSZKaYKiRJElNMNRIkqQmNB1qkpyQpGZtmxcZc2SSy5LcleRfk7x6peqVJEmD23ncBayAq4Dn9OzfO1/HJPsA/xs4HXgZ8DPAe5N8t6r+flmrlCRJQ9kRQs0Pq2rB2ZkerwZuqKrf6e5fneQQ4E3AvKEmyTpgXU/T+oEqlSRJA2v68lPXfkluSXJtkg8meewCfQ8HLpjV9kngkCS7LDBuEzDVs900VMWSJKlvrYeaS4BfA34OeCWwAbgoyU/M038DsGVW2xY6M1oPWeA8JwKTPdujhqhZkiQNoOnLT1X1iZ7drye5GPg28HLg5PmGzdrPPO2959kGbPvRgGS+rpIkaZm0PlPzH1TVD4CvA/vN02UzndmaXg8Dfgh8bxlLkyRJQ9qhQk13Qe8TgFvn6XIxcNSstqOBS6vqnuWsTZIkDafpUJPkpO5zZ/ZJ8hTgw8AEcFb38xOTnN0z5P3Ao5OcnOQJSX4deAVw0ooXL0mS+tL0mho6C3b/ls4i3+8CXwSeWlXXdz/fE9h7pnNVXZvkGODdwOuAW4DX+4waSZJWv1TNu/5VA0oyAUxNTU0xMTEx7nIkSVozpqenmZycBJisqul+xjZ9+UmSJO04DDWSJKkJhhpJktQEQ40kSWqCoUaSJDXBUCNJkppgqJEkSU0w1EiSpCYYaiRJUhMMNZIkqQmGGkmS1ARDjSRJaoKhRpIkNcFQI0mSmmCokSRJTTDUSJKkJhhqJElSEww1kiSpCYYaSZLUBEONJElqgqFGkiQ1wVAjSZKaYKiRJElNMNRIkqQmGGokSVITDDWSJKkJTYeaJJuSfDnJ7Um+k+S8JPsvMubYJDXHdv+VqluSJPWv6VADHAmcBjwVOArYGbggye6LjJsG9uzdququ5SxUkiQNZ+dxF7Ccqurne/eTHAd8BzgY+OzCQ2vzctYmSZJGq/WZmtkmuz9vW6TfHkmuT3JTko8lOWihzknWJZmY2YD1I6lWkiQt2Q4TapIEOBn4fFVduUDXa4BjgecBLwbuAr6QZL8FxmwCpnq2m0ZRsyRJWrpU1bhrWBFJTgN+AXh6VS05dCS5H3A58Nmqev08fdYB63qa1gM3TU1NMTExMUTVkiTtWKanp5mcnASYrKrpfsY2vaZmRpJT6My8HNFPoAGoqvuSfBmYd6amqrYB23rON2ipkiRpQE1ffkrHqcCvAM+qqmsHOQawEbh11PVJkqTRaX2m5jTgJcAvAbcn2dBtn6qqOwGSnA3cXFWbuvvHA18EvgVMAK+nE2pet8K1S5KkPrQeal7T/XnhrPbjgDO7v+8N3Nfz2QOBvwQ20Fn0ewWdy1ZfWrYqJUnS0HaYhcIrqXtb95QLhSVJ6s8wC4WbXlMjSZJ2HIYaSZLUBEONJElqgqFGkiQ1wVAjSZKaMNQt3Ul2oXPr8wOA71bVYi+KlCRJWhZ9z9Qk2SPJbya5kM5zXK4DvgF8t/tm69OTHDraMiVJkhbWV6hJ8kY6IeaVwD/Tef3ARmB/4HDg7XRmfz6V5PxF3mwtSZI0Mv1efnoa8LNV9fV5Pv8S8NdJXg28AjiSzusGJEmSllVfoaaq/ssS+20D3jtQRZIkSQMY+u6nJA/reVGkJEnSWAwcapIcmOQq4Fbg5iQ3J/mjJLuPrjxJkqSlGWam5gxgC/B04CDgvwHPBS5N8qAR1CZJkrRkA7+lO8kPgIOr6pqetgDnAndV1ctGU+La41u6JUkazDBv6R7m4XuXAv9hRqaqKslbgMuGOK4kSVLf+go1ST4OfBX4CvB+4N1JfqmqtvR0mwS+P7oSJUmSFtfvTM3XgZ8GjgMe3m371yR/Ryfo7NT97I0jq1CSJGkJhllT83A6C4Q39myPA+4FvllVB46qyLXGNTWSJA1mLGtqupeczu9uACTZDXhyd5MkSVox/b77ae+FPq+qO6vqi1X1F93+jxymOEmSpKXq9zk1X+6+hfuw+TokmUzyyiRX0nnhpSRJ0rLr9/LTE4C3AOcnuYfObd23AHfRub37icAB3fbfr6pPjLBWSZKkeQ20UDjJ/YFjgGcAjwF2A7YCVwCfrKorR1jjmuNCYUmSBrPiC4Wr6i7gH7qbJEnS2A3zQsuzkhzXs//oJM9NMjma0iRJkpZumBda/hxwDUD3BZaX05m5+UaS/UdQmyRJ0pINE2omgZu6v7+AzoLhSeADwIlD1iVJktSXYULNjcA+3d9/BTirqu4GTgd+ZtjCRinJa5Ncm+SuJJclecYi/X81yTeSbOv+fP5K1SpJkgYzTKg5Ezg1yYnAs4Dzuu07AXsMWdfIJHkh8B7gHXRe6/A54BPzPUgwyeHAh4Bz6DwZ+Rzg75I8ZWUqliRJgxjm3U+h88yao4GPVtWfddt/DdhUVU8YWZVDSHIJcHlVvaan7WrgvKraNEf/DwETVfXcnrbzge9X1YuXeE5v6ZYkaQDD3NI98ExNdbyjqo6cCTRdD6ezrmbskuwKHAxcMOujC4CnzTPs8Dn6f3KB/iRZl2RiZgPWD1iyJEka0MAvtJxPVf3pqI85hIfQuRy2ZVb7FmDDPGM29NkfYBNw/CAFSpKk0RjmOTU7J3lrkouTXN59bs3RoyxuhGZfY8scbcP0P5HOnV8z26P6LVCSJA1nmIXCfwy8ls6lmQ/TmRH5aDfcZBTFjcBW4F62n2V5GNvPxszY3Gd/qmpbVU3PbMDtA9YrSZIGNEyoeQnwoqo6oareWVUvA36SzhqWN42kuiF1bzG/DDhq1kdHARfNM+ziOfofvUB/SZK0CgwTanYHbu5tqKobgNcDrxqmqBE7GfiNJL+e5AlJ3g3sDbwfIMnZ3dvSZ/w5cHSSNyd5fJI3A8+hc1u4JElapYYJNZ8HXj5H+7XAnkMcd6Sq6kPA7wBvA74CHAEcU1XXd7vsTU+9VXUR8CLgOOBrwLHAC6vqkhUsW5Ik9WmY59Q8CfgCcBZwSlV9K8kuwLuAn6mqHfZhdT6nRpKkwYzrOTVXAs+k81yXbya5C7gDeCnwhkGPK0mSNIiBn1OT5OnA7VV1aJLHA0+kc9fPJf0mK0mSpGENs6bmFDp3OlFV11TVP1TVp4ANSSZHUp0kSdISDRNq9gc+M0f7zwJ/O8RxJUmS+jZMqJkGHjxH++eAw4Y4riRJUt+GCTUfZe6H7N0H7DrEcSVJkvo2TKh5C/CMJP8nyYEASe4PvJnO810kSZJWzMB3P1XV1iSHA+8DvpJkW/d4U8Avjqg+SZKkJRk41AB0n8p7TJK9gIOAe+jc0n3bKIqTJElaqqFCDUCSk+lcbvoacFVVbRu6KkmSpD4NHWro3NZ9IPALwAFJ7gWuoht0qupjIziHJEnSgoYONVX1EeAjM/tJdgOeRCfoPBsw1EiSpGU3zGsS9qqqG2e3V9WdwJe7myRJ0ooYZqbm+iTfB77a3b7S/bkOeF1V/doI6pMkSVqSYULNY4GN3e0g4D8Dj+h+5gstJUnSihrmOTXXAdcB5820dZ9bcxadB/BJkiStmGGeKLydqroYeAPw30Z5XEmSpMUMHGqS7DLPR98CDhj0uJIkSYMYZk3ND5J8A7iCziLhK4BbgN8GLhhBbZIkSUs2TKh5FvDk7vZS4J3Abt3PLkjyDn78AL6rh6pSkiRpEcMsFP488PmZ/ST3A/bnx3dEHQz8OvAwYKfhypQkSVrYKF6TAEBV3Qdc3d3+dqY9ycNHdQ5JkqT5jPTup7lU1ZblPockSdKyhxpJkqSVYKiRJElNMNRIkqQmGGokSVITmg01SR6T5Iwk1ya5M8m3k7w9ya6LjLswSc3aPrhSdUuSpMGM7JbuVejxdELbbwL/AjwJOB3YHXjTImNPB97Ws3/nchQoSZJGp9lQU1XnA+f3NP1rkv2B17B4qLmjqjYv9VxJ1gHreprWL7lQSZI0Es1efprHJHDbEvq9NMnWJFclOSnJYiFlEzDVs900ZJ2SJKlPzc7UzJZkXzov2/y9Rbr+DXAtsJnOJasT6bzf6qgFxpwInNyzvx6DjSRJK2rNhZokJwDHL9Lt0Kq6tGfMI+hcijq3qv5qoYFVdXrP7pVJvgVcmuSnq+ryecZsA7b1nG+R8iRJ0qituVADnAosdjfSdTO/dAPNp4GLgVcNcL7LgXuA/bq/S5KkVWjNhZqq2gpsXUrfJI+kE2guA47rvnSzXwcAuwC3DjBWkiStkGYXCndnaC4EbqRzt9NDk2xIsqGnzyOTXJPksO7+vkneluSQ7nNujgHOBa4AvrDyfwpJkrRUa26mpg9HA4/rbrMX7c4setkF2B94QHf/buDZwBuAPegEoo8Db6+qe5e7YEmSNLhU1bhraE6SCWBqamqKiYmJcZcjSdKaMT09zeTkJMBkVU33M7bZy0+SJGnHYqiRJElNMNRIkqQmGGokSVITDDWSJKkJhhpJktQEQ40kSWqCoUaSJDXBUCNJkppgqJEkSU0w1EiSpCYYaiRJUhMMNZIkqQmGGkmS1ARDjSRJaoKhRpIkNcFQI0mSmmCokSRJTTDUSJKkJhhqJElSEww1kiSpCYYaSZLUBEONJElqgqFGkiQ1wVAjSZKasPO4C5C0Y/rhD3841PiddtppRJX0L8nYzi1pfs7USJKkJjQdapJcl6RmbX+8yJgkOSHJLUnuTHJhkgNWqmZJkjSYpkNN19uAPXu2P1qk/x8Avwv8FnAosBn4VJL1y1mkJEkazo6wpub2qtq8lI7pXCj/HeAdVfUP3baXA1uAlwB/Mc+4dcC6niYDkCRJK2xHmKl5c5LvJflKkrcm2XWBvvsAG4ALZhqqahvwGeBpC4zbBEz1bDcNX7YkSepH6zM1fw5cDnwfOAw4kU5w+Y15+m/o/twyq30L8OgFznMicHLP/noMNpIkrag1F2qSnAAcv0i3Q6vq0qp6d0/b15J8H/hwkjdX1fcWGF+zTztH2487d2ZztvXUuEh5kiRp1NZcqAFOBT64SJ/r5mn/Yvfn44C5Qs3M2psNwK097Q9j+9kbSZK0iqy5UFNVW4GtAw4/qPvz1nk+v5ZOsDkKuAKguwbnSODNA55TkiStgGYXCic5PMkbk2xMsk+SF9C5e+mjVXVDT79rkjwfoKoKeA/wliTPT/Ik4EzgDuADK/+nkCRJS7XmZmr6sA14IZ31N+uA64HTgT+Z1W9/YLJn/0+A3YD3Ag8CLgGOrqrbl7tgaRCdLD6YzZuX9LSDeX3oQx8aeOxJJ5001LlvvvnmgceeccYZQ537iCOOGHjsYx/72KHOfb/7NfvfotLQmg01VXU58NQl9Mus/QJO6G6SJGmNMPJLkqQmGGokSVITDDWSJKkJhhpJktQEQ40kSWqCoUaSJDXBUCNJkppgqJEkSU1o9uF70o5imCcK77bbbkOd+6CDDlq80zxe8YpXDHXuO+64Y+Cxj3jEI4Y690Mf+tCBx953331DndsnCkvz858OSZLUBEONJElqgqFGkiQ1wVAjSZKaYKiRJElNMNRIkqQmGGokSVITDDWSJKkJhhpJktSEDPM0Us0tyQQwNTU1xcTExLjLkSRpzZienmZychJgsqqm+xnrTI0kSWqCoUaSJDXBUCNJkppgqJEkSU0w1EiSpCYYaiRJUhMMNZIkqQnNhpokz0xS82yHLjDuwjn6f3Ala5ckSf3bedwFLKOLgD1ntf0P4DnApYuMPR14W8/+nSOsS5IkLYNmQ01V3Q1sntlPsgvwPODUWvwxyndU1eZF+kiSpFWk2ctPc3ge8BDgzCX0fWmSrUmuSnJSkvULdU6yLsnEzAYs2F+SJI1eszM1c3gF8MmqunGRfn8DXEtnludJwInAk4GjFhizCTh+FEVKkqTBrLkXWiY5gcUDxKFV9aN1M0keBVwPvKCq/r7P8x1MZw3OwVV1+Tx91gHreprWAzf5QktJkvozzAst1+JMzanAYncjXTdr/zjge8BHBzjf5cA9wH7d37dTVduAbTP7SQY4jSRJGsaaCzVVtRXYutT+6SSM44Czq+qeAU55ALALcOsAYyVJ0grZERYKPwvYBzhj9gdJHpnkmiSHdff3TfK2JIckeUySY4BzgSuAL6xo1ZIkqS9rbqZmAK8ALqqqq+f4bBdgf+AB3f27gWcDbwD2AG4EPg68varuXYFaJUnSgJoPNVX1kgU+uw5Iz/6NwJErUJYkSRqxHeHykyRJ2gEYaiRJUhMMNZIkqQmGGkmS1ARDjSRJaoKhRpIkNcFQI0mSmmCokSRJTTDUSJKkJhhqJElSEww1kiSpCYYaSZLUBEONJElqgqFGkiQ1wVAjSZKaYKiRJElNMNRIkqQmGGokSVITDDWSJKkJhhpJktQEQ40kSWqCoUaSJDXBUCNJkppgqJEkSU0w1EiSpCYYaiRJUhMMNZIkqQlrNtQkeWuSi5LckeTf5umzd5J/SvKDJFuT/M8kuy5y3AclOSfJVHc7J8kDl+dPIUmSRmXNhhpgV+Bc4H1zfZhkJ+DjwO7A04EXAb8K/Nkix/0AsBH4+e62EThnNCVLkqTlsvO4CxhUVR0PkOTYebocDTwR2Kuqbun2/T3gzCRvrarp2QOSPIFOkHlqVV3SbXslcHGS/avqm6P/k0iSpFFYs6FmCQ4HrpwJNF2fBNYBBwOfnmfM1EygAaiqLyaZAp4GzBlqkqzrHnfGeoDp6e1ykyRJWsAw/+5sOdRsALb0NlTV95Pc3f1svjHfmaP9OwuMAdgEHD+7ca+99lpapZIkabYHA30lnFUVapKcwBzhYJZDq+rSJR6y5jrNPO3DjDkROLlnfz1wE/Ao4PZFalSH39lg/N7653c2GL+3/vmdDWbme7ut34GrKtQApwIfXKTPdUs81mbgKb0NSR4E7MKsGZxZYx4+R/tDFxhDVW0DtvWcZ+bX2+dau6Pt+Z0Nxu+tf35ng/F765/f2WB6vre+rapQU1Vbga0jOtzFwFuT7FlVt3bbjqYTPi5bYMxkksOq6ksASZ4CTAIXjaguSZK0DNbsLd3dZ9BsBPYGdkqysbvt0e1yAfAN4JwkByV5NnAScPpMYk5yWJJrkjwSoKquBs4HTk/y1CRPBU4HPuadT5IkrW5rNtQA/x24Ang7sEf39yuAQwCq6l7gF4C7gC8AfwecB7yp5xgPAPanc0lqxkuBr9MJRRcAXwP+a5+1bevWtW2xjvoRv7PB+L31z+9sMH5v/fM7G8zA31uqFlr/KkmStDas5ZkaSZKkHzHUSJKkJhhqJElSEww1kiSpCYaaZZbko0luSHJXkluTnJPkEeOua7VK8pgkZyS5NsmdSb6d5O1Jdh13batdkrcmuSjJHUn+bdz1rFZJXtv9+3VXksuSPGPcNa1mSY5I8k9JbklSSX553DWtdkk2JflyktuTfCfJeUn2H3ddq1mS1yT5WpLp7nZxkuf2exxDzfL7NPACOreO/yqwL/DhsVa0uj2ezt/L3wQOAN4IvBp45ziLWiN2Bc4F3jfuQlarJC8E3gO8AzgI+BzwiSR7j7Ww1W134KvAb427kDXkSOA04KnAUXQedHtBkt3HWtXqdhPwh3Qey3II8M/AR5Ic0M9BvKV7hSV5Hp3n5ayrqnvGXc9akOT3gddU1WPHXctakORY4D1V9cBx17LaJLkEuLyqXtPTdjVwXlVtGl9la0OSAp5fVeeNu5a1JMlD6bwY+ciq+uy461krktwG/H5VnbHUMc7UrKAkD6bzcL+LDDR9mWSAF5tJvbqXMA+m81DNXhcAT1v5irQDmez+9P/HliDJTkleRGeW8OJ+xhpqVkCSdyX5AfA9Oq91+KUxl7RmJNkX+G3g/eOuRWveQ4Cd2P7ltFuADStfjnYE6byd8WTg81V15bjrWc2S/FSSf6fzJOH305kV/EY/xzDUDCDJCd0Fcwtth/QM+VM61++PBu4Fzs4wryFdgwb4zuguqD4fOLeq/mo8lY/XIN+bFjX7mnvmaJNG5VTgQODF4y5kDfgmsJHOWqT3AWcleWI/B1hVb+leQ04FPrhIn+tmful5+/j/616/v5HO/2h9TautcX19Z91A82k639Grlq+sVa+v700L2krnPypmz8o8jO1nb6ShJTkFeB5wRFXdNO56Vruquhv4l+7upUkOBd5A58aRJTHUDKAnpAxiZoZm3YjKWRP6+c66b03/NHAZcFxV3becta1mQ/5dU4+qujvJZXTuRvnHno+OAj4ynqrUou5M/CnA84FnVtW1Yy5prQp9/rvSULOMkhwGHAZ8Hvg+8Fg6bxf/NjvWLM2SdWdoLgRuoPNG9YfOXKmrqs3jq2z1696W/GA667Z2SrKx+9G/VNW/j6+yVeVk4Jwkl/LjWcC9cc3WvJLsATyup2mf7t+t26rqhjGVtdqdBryEzvrJ25PMzA5OVdWd4ytr9UryTuATdK5krAdeBDwT+Pm+juMt3csnyU8Bfw48mc4q7lvprBH5o6q6eZy1rVbd25H/11yfVdUOtQ6pX0nOBF4+x0c/W1UXrmw1q1eS1wJ/AOwJXAm80dts55fkmXRmTmc7q6qOXdlq1obure9zOa6qzlzJWtaKJGcAz6bzz+UU8DXgXVX1qb6OY6iRJEkt8O4nSZLUBEONJElqgqFGkiQ1wVAjSZKaYKiRJElNMNRIkqQmGGokSVITDDWSJKkJhhpJktQEQ40kSWqCoUZSs5I8Pck9Sdb1tO2TpJI8epy1SRo9Q42klm0Erq6qbbPa/q2qrh9TTZKWiaFGUsueDFwxq20j8NWZnST/KckpK1qVpGVhqJHUso3AV2a1HURPqAEOZPvgI2kNMtRIalKSnYAD2D6w/DT/MegcCDwhyWVJrkrykytVo6TRMtRIatX+wG7ALTMNSQ4HHsn2MzU3VNXBwCnA765kkZJGx1AjqVUbuz9/O8l+SZ4LnN1tWwfQvStqV+C0bvvXgIesaJWSRsZQI6lVG4FPAfsAVwLvBP4QmAZe1+1zAPCNqrqvu38QnWAjaQ3aedwFSNIyeTJwWVVtmtX+9z2/Hwjsm2QX4EHAscBzV6Y8SaNmqJHUqicDZy7S56fohJwv05m5fmNVbV3muiQtk1TVuGuQpJFKsgG4FXgi4KndAAAAUElEQVRSVV017nokrQxDjSRJaoILhSVJUhMMNZIkqQmGGkmS1ARDjSRJaoKhRpIkNcFQI0mSmmCokSRJTTDUSJKkJhhqJElSEww1kiSpCf8f8fjaGGc704kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.figure(dpi=100)\n",
    "\n",
    "ranges = ((-3, 3), (-10, 10))\n",
    "extent = sum(ranges, ())\n",
    "H, _, _ = np.histogram2d(galactic_coords_w_pm.pm_b.to_value(),\n",
    "                         galactic_coords_w_pm.pm_l_cosb.to_value(),\n",
    "                         bins=(30, 30),\n",
    "                         range=ranges)\n",
    "\n",
    "plt.imshow(H.T, cmap='gray_r', origin='lower', extent=extent, aspect='auto')\n",
    "plt.xlabel(r\"$\\mu_b$\")\n",
    "plt.ylabel(r\"$\\mu_l cos(b)$\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_xyz_to_pm_matrix(coord):\n",
    "\n",
    "    raw_rows = []\n",
    "    for vx, vy, vz in [(1,0,0), (0,1,0), (0,0,1), (0, 0, 0)]:\n",
    "        gc_coord = Galactocentric(x=coord.galactocentric.x,\n",
    "                                  y=coord.galactocentric.y,\n",
    "                                  z=coord.galactocentric.z,\n",
    "                                  v_x = vx * u.km/u.s,\n",
    "                                  v_y = vy * u.km/u.s,\n",
    "                                  v_z = vz * u.km/u.s,\n",
    "                                  # galcen_v_sun=coord.CartesianDifferential([0,220,0]*u.km/u.s)\n",
    "                                  )\n",
    "        \n",
    "        galactic_coords_w_pm = gc_coord.transform_to(Galactic)\n",
    "        raw_rows.append(np.array([galactic_coords_w_pm.pm_b.to_value(),\n",
    "                                  galactic_coords_w_pm.pm_l_cosb.to_value(),\n",
    "                                  galactic_coords_w_pm.radial_velocity.to_value()]))\n",
    "\n",
    "\n",
    "    center = raw_rows[-1]\n",
    "    matrix_rows = [row - center for row in raw_rows[:-1]]\n",
    "    return np.matrix(matrix_rows), center\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-7.26024968e-03, -1.60973685e-02, -5.47027551e-01],\n",
       "        [ 8.64080724e-03, -1.36333788e-02,  6.43839769e-01],\n",
       "        [ 1.78219598e-02,  5.23221446e-05, -5.35005804e-01]])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#coord = astropy.coordinates.SkyCoord(180*u.deg, 0*u.deg, 10*u.kpc, frame='galactic')\n",
    "ra, dec = 20, 30\n",
    "coord = astropy.coordinates.SkyCoord(ra*u.deg, dec*u.deg, 10*u.kpc, frame='icrs')\n",
    "m, center = make_xyz_to_pm_matrix(coord)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 300\n",
    "velocities = np.random.multivariate_normal([0, 0, 0], np.diag([80, 80, 80]), n_samples)\n",
    "converted_vels = np.matmul(np.matrix(velocities), m) + center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuEAAAGYCAYAAADshTwGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAQ6gAAEOoBgtMKmAAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X28XGV56P3fBQhIMRKkGgQ0FOshQUvgCBU5rW9HXitoBRE9ij2+YKmVFtEGilBjtdH2iVXsU4tooc8DiKkiIhRBCFCLChIQMYk9vATFhjcNRhGChOv8sdaYxWRmZ+/ZM/fs2fv3/XzmM3vuda8115o9a+5r7rnXvSIzkSRJklTOFsMOQJIkSZppTMIlSZKkwkzCJUmSpMJMwiVJkqTCTMIlSZKkwkzCJUmSpMJMwiVJkqTCTMIlSZKkwkzCJUmSpMJMwiVJkqTCTMIlSZKkwkzCJUmSpMK2GnYAmxMR2wAvBB4ANgw5HGmytgR+E/heZq4fdjBTice6phmP9TF4vGsa6flYn/JJONVBeuOwg5D6bD/gO8MOYorxWNd05LHemce7ppsJH+ujkIQ/AHDDDTew8847DzsWaVLWrFnD/vvvD/X7Wk/isa5pw2N9szzeNS1M5lgfhSR8A8DOO+/MrrvuOuxYpH7x59dNeaxrOvJY78zjXdPNhI91T8yUJEmSCjMJlyRJkgozCZckSZIKMwmXJEmSCjMJlyRJkgozCZckSZIKMwmXJEmSCjMJlyRJkgozCZckSZIKMwmXJEmSCjMJlyRJkgozCZckSZIKMwmXJEmSCttq2AFo6pi78NJJrb968eF9ikTSWDxWJY3HZD8rwM+LQbInXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqrOckPCLeFhHfjYhHI+L+iPhK2/LjImJVvfy2iDh68uFKkqRBsW2XyukpCY+IvwKWAOcBBwPHA2say48CzgEuAg4FrgIujIiDJheuJEkaBNt2qawJX7Y+IuYBpwGHZeYVjUUXNf7+ELA0M0+pHy+LiD2BRUBzHUlTWERcA7y0y+JjM/Pzdb3jgFOAucDtwAczc2mJGCVNnm27VF4vPeFvBe5sO0h/LSJ2B/YELmhbdD6wf0Ts1MNzShqOE4AD2m4XAo8DXwd7x6Rp4q3YtktF9ZKEvxj4XkR8oB4v9lhEXBsRC+rl8+r7lW3rrQCC6iCWNAIyc0Vmfqt5A/YHrsjMB+tqv+4dy8xlmXkicCVV75ik0WDbLhU24eEowBxgX2Av4F3AY8AZwJUR8dvA7LreQ23rra3vdxxr4xExC5jV9nySpoCIeAmwO9XP1s3esVPbqp4P/HNE7NRI1iVNXbbtUmG9JOFbANsDr8vM7wNExE3AXcA7gR/X9bJtvehS3u4kqgNf0tTzRuCXwMX14/H0jn2jTGiSJsG2XSqsl+EoPwXuax2kAJm5BlhF9Q269a14dtt6O9T3axnbEmC3xm2/HmKU1GcRsRVwNHBxZj5cF/fcOxYRsyJi19YNe8akYbJtlwrrpSd8JfDcDuUBPMHGHrF5VAdvy3yqb8qrGENmrgPW/XqjEWPUllTQq4BnUg01addL75g9Y9LUYdsuFdZLT/hXgWdFxAtaBRGxC9XPzt/NzLuoDsZj2tY7FrjB8aHSyHoj8BPga42yyfSO2TMmTR227VJhvfSEXwQsB74UEadRnbxxOvAA8Jm6zulUU5TdQTVLwpHAQcAhk45YUnER8VSq4/i8zPxVY1HPvWP2jElTim27VNiEe8IzcwPVXMA3AmdRXVnrXuCVrXGi9UU6/gg4iqrX7CDgmG7zj0qa8o4AnkbbUBR7x6TpwbZdKq+XnnAy837gTZupcy5wbi/blzTlvBH4IZ1nOrF3TJoGbNulsnpKwiXNHBExmyqh/vvM3OREy8xcGhHbUc0VfjLVZevtHZMkaQwm4ZLGlJlrgW02U8feMUmSJqCX2VEkSZIkTYI94dPI3IWXDjsESZIkjYM94ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYVsNOwBJ0uiZu/DSSa2/evHhfYpEkkaTPeGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhTlEoSZI0BU12KlBNbfaES5IkSYWZhEuSJEmFmYRLkiRJhZmES9qsiHhbRHw3Ih6NiPsj4itty4+LiFX18tsi4uhhxSpJ0igwCZc0poj4K2AJcB5wMHA8sKax/CjgHOAi4FDgKuDCiDiodKySJI0KZ0eR1FVEzANOAw7LzCsaiy5q/P0hYGlmnlI/XhYRewKLgOY6kiSpZk+4pLG8FbizLQH/tYjYHdgTuKBt0fnA/hGx02DDkyRpNJmESxrLi4HvRcQH6rHgj0XEtRGxoF4+r75f2bbeCiCoEnRJktTG4SiSxjIH2BfYC3gX8BhwBnBlRPw2MLuu91Dbemvr+x27bTgiZgGz2p5LkqQZYcI94RHx1ojIDrfFbfWcLUEafVsA2wOvy8wvZeZXgSOApwHvbNTLtvWiS3nTScCPGrcb+xKxpAmzbZfKm0xP+CHAzxqPf9z6ozFbwmKqE7NeQzVbws+6jS2VNCX9FLgvM7/fKsjMNRGxiqp3/La6eDZwX2O9Her7tXS3BDi78XgOJuLSsNm2S4VMJgm/KTMf7LLM2RKk6WEl8NwO5QE8wcax4POAVY3l86l6wVfRRWauA9b9eoMR3apKKse2XSqk7ydmOluCNK18FXhWRLygVRARu1Ad49/NzLuoEu1j2tY7FrhhjMZc0gixbZf6bzJJ+PcjYkNE3BkRp0TElnX5pGZLiIhZEbFr64Yna0nDdBGwHPhSRLw+Il4DXAI8AHymrnM68PqI+HBEvCwiPg4cVJdLGi227VIhvQxHWUM1O8K3qX5uPgL4a2AX4N1MYraE2kn19iUNWWZuiIhDgY8DZwFPAa4Fjs3Mh+s6SyNiO+BU4GTgduAYx4hKI8W2fQDmLrx02CFoCptwEp6ZXwO+1ii6IiIeAf48Ij7crNq26nhmSwBP1pKmlMy8H3jTZuqcC5xbJiJJ/WbbLpXXrzHhXwC2BBaw8Vvx7LY645ktgcxcl5n3tG7AvX2KUZIkjZ9tuzRA/UrCm9MaNGdLaNrsbAmSJGnKsG2XBqhfSfgxwAbgZmdLkCRpWrBtlwZowmPCI+JrwFVsvEjHEVRXzvtEZrZ+XjqdagL/O4ArgSOpZks4ZNIRS5KkvrJtl8rrZXaUVcDbgV2petL/E/gz4MxWBWdLkCRppNi2S4X1MjvKicCJ46jnbAmS1IHTlmmqsW2Xyuv7FTMlSZIkjc0kXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqbKthB6DpY+7CSye9jdWLD+9DJJIkSVObPeGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYZ6YOUX046RGSZIkjQZ7wiVJkqTCTMIlSZKkwkzCJUmSpMJMwiVJkqTCTMIlSZKkwkzCJUmSpMJMwiV1FRFvjYjscFvcVu+4iFgVEY9GxG0RcfSwYpYkaRQ4T7ik8TgE+Fnj8Y9bf0TEUcA5wGLgCuA1wIUR8bPMvKJkkJIkjQqTcEnjcVNmPthl2YeApZl5Sv14WUTsCSyiSsolSVIbh6NI6llE7A7sCVzQtuh8YP+I2Kl8VJIkTX0m4ZLG4/sRsSEi7oyIUyJiy7p8Xn2/sq3+CiCoEvSOImJWROzaugFz+h+2JElTk8NRJI1lDXAG8G0ggSOAvwZ2Ad4NzK7rPdS23tr6fscxtn1SvW1JkmYck3BJXWXm14CvNYquiIhHgD+PiA83q7atGl3Km5YAZzcezwFu7DVWSZJGicNRJE3UF4AtgQVs7PGe3VZnh/p+LV1k5rrMvKd1A+7te6SSJE1RJuGSJioaf7fGgs9rqzOfqhd8VZGIJEkaMSbhkibqGGADcHNm3kWVaB/TVudY4IYxpjWUJGlGc0y4pK4i4mvAVcBtddERwDuBT2Rma/jI6VQX57kDuBI4EjiI6gI/kiSpA5NwSWNZBbwd2JXql7P/BP4MOLNVITOXRsR2wKnAycDtwDFeLVOSpO5MwiV1lZknAieOo965wLmDj0iSpOlh0kl4RGxP1Vu2C7BfZn6nsew44BRgLlXv2Aczc+lkn1OSJA2G7bqa5i68dFLrr158eJ8imX760RP+gU7biYijgHOAxcAVwGuoxo3+zJ+pJWl4JtuoatqzXZcKmNTsKBGxJ/AndL7q3YeApZl5SmYuq3/WvhJYNJnnlCRJg2G7LpUz2SkKPwl8GvhBszAidgf2BC5oq38+sH9E7DTJ55UkSf1nuy4V0vNwlPpnqb2Bo4B92xa3Ltyxsq18BdWFPvYEvtFlu7OAWY2iOb3GKEmSxmdQ7Xq9bdt2qU1PSXg9HdkS4JTMXBcR7VVal7B+qK28dQnrHcfY/El0/hlMkiQNwIDbdbBtlzbR63CU04D7qE7QGEu2PY4u5U1LgN0at/16iE+SJI3fINt1sG2XNjHhnvCIeC7wXuC1wKz62/L29eLt66mNWt+MZ1Md1C071Pdr6SIz1wHrGs830RAlSdI4DbpdB9t2qZNehqPsDmwNdJrjahnwbeDY+vE8qrlGW+ZTfVtehSRJmgps16Uh6CUJvwV4eVvZAuDjwLuAGzPzrohYBRwDXNSodyxwQ2Y+2EuwkiSp72zXpSGYcBKemQ8B1zTLGj8r3ZSZy+u/T6eaxP8OqnlEjwQOAg7pNVhJktRftuvScPTjipkdZebS+mzrU4GTqS5ve4xX1ZIkafTYrkv91ZckPDOvYeMZ0s3yc4Fz+/EckiSpDNt1afAme8VMSZIkSRNkEi5JkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQVttWwA5A0GiJie2AVsAuwX2Z+p7HsOOAUYC5wO/DBzFw6jDglqV/mLrx02CFoGrMnXNJ4fYAOX9wj4ijgHOAi4FDgKuDCiDioaHSSJI0Qk3BJmxURewJ/ApzRYfGHgKWZeUpmLsvME4ErgUUlY5QkaZSYhEsaj08CnwZ+0CyMiN2BPYEL2uqfD+wfETuVCU+SpNHimHBJY6qHm+wNHAXs27Z4Xn2/sq18BRBUCfo3umx3FjCrUTRn0sFKkjQiTMIldRUR2wFLgFMyc11EtFeZXd8/1Fa+tr7fcYzNn0Tn4S2SJE17DkeRNJbTgPuoTrwcS7Y9ji7lTUuA3Rq3/XqIT5KkkWRPuKSOIuK5wHuB1wKz6l7w7evF29dTFrZ6vGdTJestO9T3a+kiM9cB6xrP15/AJUkaASbhkrrZHdga6DRR7jLg28Cx9eN5VHOIt8yn6gVfhSRJ2oRJuKRubgFe3la2APg48C7gxsy8KyJWAcdQzRPecixwQ2Y+WCRSSZJGjEm4pI4y8yHgmmZZY8jITZm5vP77dKqL89xBNT/4kcBBwCFlIpUkafSYhEualMxcWs+icipwMtVl64/JzCuGG5kkSVOXSbikccvMa9g480mz/Fzg3OIBSZI0opyiUJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSpswkl4RBwcEddGxAMRsT4i7oyIJRHx9LZ6h0XEzRHxaETcHhEn9C9sSZLUL7btUnm99ITvCFwPvBM4GFgCvAVY2qoQEQcAFwPLgUOBc4AzI+Ltk4xXkiT1n227VNiEr5iZmRcAFzSKromI9cBZEfHszPwv4HRgeWa+ra6zLCKeAyyKiM9l5hOTjlySJPWFbbtUXr/GhP+kvn9KRGwDvAL4fFud84CdgX369JySJGlwbNulAeo5CY+ILSNi24jYl+rb8SWZeTewB7A1sLJtlRX1/bxen1OSJA2ObbtUzoSHozTcDexS/305cGz99+z6/qG2+mvr+x3H2mhEzAJmNYrmTCJGSZI0frbtUiGTGY5yGHAg1UkcewGXRMSWjeXZZb1u5S0nAT9q3G6cRIySJGn8bNulQnruCc/MW+s/r4+I5cB3gNey8aep2W2rtB6vZWxLgLMbj+fgwSpJ0sDZtkvl9OvEzFuADcDzgDuAx9h0fNj8+r59PNmTZOa6zLyndQPu7VOMkiRp/GzbpQHqVxJ+ALAlcGdmrgeuBl7fVudYYA1wc5+eU5IkDY5tuzRAEx6OEhFfovp56lbgEWBv4P314y/X1RYB10XEZ6imLzoQeAdwvPOISpI0tdi2S+X1Mib8BuAYYCFVT/pq4Czg7zLzMYDM/GZEHAl8hOqKW/cA78nMsztucRqYu/DSYYcgSVKvbNulwnq5YuZiYPE46l0GXNZLUJIkqRzbdqm8fo0JlyRJkjROJuGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZI6ioiDI+LaiHggItZHxJ0RsSQint5W77CIuDkiHo2I2yPihGHFLEnSqDAJl9TNjsD1wDuBg4ElwFuApa0KEXEAcDGwHDgUOAc4MyLeXjpYSZJGyVbDDkDS1JSZFwAXNIquiYj1wFkR8ezM/C/gdGB5Zr6trrMsIp4DLIqIz2XmE4XDliRpJNgTLmkiflLfPyUitgFeAXy+rc55wM7APiUDkyRplNgTrill7sJLJ7X+6sWH9ykStUTElsBTgPlUPd+XZObdETEf2BpY2bbKivp+HnBTsUAlSRohJuGSNuduYJf678uBY+u/Z9f3D7XVX1vf7zjWRiNiFjCrUTRnEjFKkjRSHI4iaXMOAw6kOkFzL+CSune8Jbus16285STgR43bjZOMU5KkkWFPuKQxZeat9Z/XR8Ry4DvAa9k47GR22yqtx2sZ2xLg7MbjOZiIS5JmCHvCJU3ELcAG4HnAHcBjVGO/m+bX9+1jxZ8kM9dl5j2tG3Bvv4OVJGmqMgmXNBEHAFsCd2bmeuBq4PVtdY4F1gA3F45NkqSR4XAUSR1FxJeohp7cCjwC7A28v3785braIuC6iPgM1dSEBwLvAI53jnBJkrozCZfUzQ3AMcBCql/NVgNnAX+XmY8BZOY3I+JI4CNUV9O8B3hPZp7dcYuSJAkwCZfURWYuBhaPo95lwGWDj0iSpOnDMeGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhTlEoSZKkgZi78NJJrb968eF9imTqsSdckiRJKswkXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSpswkl4RBwdEV+OiB9FxMMRcWtE/HFEbNFW77CIuDkiHo2I2yPihP6FLUmS+sW2XSqvl57w9wLrgfcBfwB8Gfgk8NFWhYg4ALgYWA4cCpwDnBkRb59kvJIkqf9s26XCerls/asz84HG42URsT3w7og4LTPXA6cDyzPzbY06zwEWRcTnMvOJScYtSZL6x7ZdKmzCPeFtB2nLzcC2wI4RsQ3wCuDzbXXOA3YG9pnoc0qSpMGxbZfK69eJmb8H/BS4H9gD2BpY2VZnRX0/r0/PKUmSBse2XRqgXoajPElEvAj4I+CDmbkhImbXix5qq7q2vt9xM9ubBcxqFM2ZbIySJGn8bNulwZtUEh4Rc4AvAjfQOHmjll1W61bechJwxmTikiRJvbFtl8roeThKRDwd+Dfgl8ARmfmrelHrW/HstlVmty3vZgmwW+O2X68xSpKk8bNtl8rpqSc8IrYFvgI8CzggM3/SWHwH8BjV+LDLG+Xz6/v28WRPkpnrgHWN5+olREnSFDZ34aWTWn/14sP7FIlabNulsnq5WM9WwBeAvYFDMvPu5vJ6GqOrgde3rXossIbqbGtJkjRF2LZL5fXSE/4PwKuB9wPbRcSLG8tW1N92FwHXRcRnqKYvOhB4B3C884hKkjTl2LZLhfWShB9c33+sw7KXA9dk5jcj4kjgI8BbgHuA92Tm2b2FKUmSBsi2XSpswkl4Zs4dZ73LgMsmun1JklSWbbtUXr8u1iNJkiRpnEzCJUmSpMJMwiV1FBFHR8SXI+JHEfFwRNwaEX8cEVu01TssIm6OiEcj4vaIOGFYMUuSNCpMwiV1815gPfA+4A+ALwOfpHEFvYg4ALgYWA4cCpwDnBkRby8drCRJo2RSl62XNK29OjMfaDxeFhHbA++OiNPqeYNPB5Zn5tsadZ4DLIqIzzltmSRJndkTLqmjtgS85WZgW2DHiNgGeAXw+bY65wE7A/sMNkJJkkaXSbikifg94KfA/cAewNZsernqFfX9vIJxSZI0UhyOImlcIuJFwB8BH8zMDRExu170UFvVtfX9jpvZ3ixgVqNoTl8ClSRpBJiES9qsiJgDfBG4gcaJmbXsslq38paTgDMmGVpxcxdeOuwQJEnTgMNRJI0pIp4O/BvwS+CIzPxVvajV4z27bZXZbcu7WQLs1rjtN/loJUkaDfaES+oqIrYFvgI8CzggM3/SWHwH8BjV2O/LG+Xz6/v2seJPkpnrgHWN5+pHyJIkjQR7wiV1FBFbAV8A9gYOycy7m8vrKQqvBl7ftuqxwBqqmVQkSVIH9oRL6uYfgFcD7we2i4gXN5atqHuyFwHXRcRnqKYmPBB4B3C8c4RLktSdSbikbg6u7z/WYdnLgWsy85sRcSTwEeAtwD3AezLz7EIxSpI0kkzCJXWUmXPHWe8y4LLBRiNJ0vTimHBJkiSpMJNwSZIkqTCTcEmSJKkwx4RLkqRpx6vbaqqzJ1ySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqzBMzJUmSNCX14wTb1YsP70Mk/WdPuCRJklSYSbgkSZJUmEm4JEmSVJhJuCRJklSYJ2ZqWpnsCRxT9eQNSZI0vdgTLkmSJBVmEi5JkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQV1lMSHhHPi4hPR8QtEfF4RNzWpd5xEbEqIh6NiNsi4ujJhStJkvrNdl0qr9ee8L2Aw4HbgRWdKkTEUcA5wEXAocBVwIURcVCPzylJkgbDdl0qrNfL1l+SmRcDRMQ5wIs61PkQsDQzT6kfL4uIPYFFwBU9Pq8kSeo/23WpsJ56wjPzibGWR8TuwJ7ABW2Lzgf2j4idenleSZLUf7brUnmDOjFzXn2/sq18BRBUB7IkSRoNtutSn/U6HGVzZtf3D7WVr63vd+y2YkTMAmY1iub0MS5JkjRxPbfrYNsudTKoJLwl2x5Hl/Kmk4AzBhNOd3MXXlr6KSVJGjW9tOswpLZdmsoGNRyl9c14dlv5Dm3LO1kC7Na47dff0CRJ0gRNpl0H23ZpE4NKwltjxua1lc+n+ra8qtuKmbkuM+9p3YB7BxSjpM1w7mBJtZ7bdbBtlzoZSBKemXdRHZDHtC06FrghMx8cxPNK6jvnDpZkuy4NQE9jwiNiO+Cw+uFzgVl1QwxwbWY+AJxO1RDfAVwJHAkcBBwyuZAlFeTcwdIMYLsuldfriZnPBJa2lbUevxy4JjOX1gf1qcDJVD1px2SmjbI0IiYwd/CpbYvOB/45Inayh0waCbbrUmE9JeGZuZqNZ0SPVe9c4NxenkPSSBjP3MHf6LSiU5ZJU4ftulTeoKcolDS9TWbuYKcskyTNWIOaHUXSzNLL3MFOWSZJmrHsCZc0Gc25g+9rlG927uDMXAesaz2O2Owv4ZIkTRv2hEuajEnNHSxJ0kxlEi6pZ84dLElSbxyOIqkr5w7WVDV34aWTWn/14sP7FIkk9cYkXNJYnDtYkqQBMAmX1JVzB0uSNBiOCZckSZIKMwmXJEmSCjMJlyRJkgozCZckSZIKMwmXJEmSCjMJlyRJkgozCZckSZIKMwmXJEmSCjMJlyRJkgozCZckSZIKmxaXrZ+78NJhh6Bpoh/vpdWLD+9DJJIkaTqzJ1ySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSrMJFySJEkqzCRckiRJKswkXJIkSSpsWlwxU5IkTS9eDVv9Mtn30qCuhG1PuCRJklSYPeGSZhR71yRJU4E94ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhJuGSJElSYSbhkiRJUmEm4ZIkSVJhA03CI+L5EXF5RDwcEfdHxCci4qmDfE5J5XmsSzODx7rUPwO7WE9E7ABcDdwNvA54JrAEeAbwvwb1vJLKKn2se7Ed9UM/3keDupT1VGW7LvXXIK+YeTwwG1iQmQ8CRMTjwHkR8eHMXDnA55ZUjse6NDN4rEt9NMjhKIcBX28dqLUvAuvrZZKmB491aWbwWJf6aJA94fOAzzULMnN9RNxRL+soImYBsxpFuwCsWbOm6xM9vu7Brsuk0u65556uyxrv4y2LBFNGsWMdPN41dXisj+9YB9t2jbZBHeuDTMJnAw91KF8L7DjGeicBZ7QX7r///n0KSxqs3f5xXNV+k2pc5XTgsa4ZyWP91zZ3rIPHu0bYoI71QSbhANmhLLqUtywBzm483hrYDVgNbJjg888BbgT2A+6d4Lp6Ml/L/tiS6kD93rAD6bNhH+vD4DGxka/Fk7Vej8PxWG8pebzP5PfjTN33Ye53z+36IJPwtVTfmtvtAHQ9eSMz1wHr2orv7CWAiGj9eW9mdv8tQZvla9lX06VXrGXox/oweExs5GvxZI3X49bMXD/MWPqsp2Mdyh7vM/n9OFP3fQrsd0/t+iBPzFxJ2xixiNgG2IPNHKySRorHujQzeKxLfTTIJPwy4JUR8YxG2WuBbeplkqaVO8pUAAARK0lEQVQHj3VpZvBYl/pokEn4P1GdwHFxRBwcEW8GzgTOKziX6Drgg2z6E5gmztdS3UyFY30YPCY28rV4sun6eozKsT5dX//xmKn7PpL7HZmbO5diEhuPeD7VAfo/gF8CFwB/kZmPDOxJJRXnsS7NDB7rUv8MNAmXJEmStKlBDkeRJEmS1IFJuCRJklSYSbgkSZJU2LRKwiNiy4h4f0RcGxEPRMTaiLguIl45zvWzw20mXXEK6Mvr+JSI+JuIWBMRv4yIZRHxO4OOW+qHetaH1nt/fUTcGRFLIuLp41j3uIhYFRGPRsRtEXF0iZgHqdfXIyKu6fKZumep2ActIraPiHvq/XrROOqfHBGr6/fHjRHxsgJhTlsR8aqIOD8i7qj/B58ao+5eEfGViPhZRPwiIr4TES9pq/P8iLg8Ih6OiPsj4hMR8dTB78nETGS/G+t8olvdUdlvGN++1/tzZkSsqPfp7oj4bETM6VB3TkRcGBHrIuKhiPiXiNixzN5MsyQceCpwKnAL8EfAG4AfA1dGxB+McxtnAgc0bocNIM6pbrKv48eBPwFOB44EHgeu6nQASFPQjsD1wDuBg6kut/0WYOlYK0XEUcA5wEXAocBVwIURcdAggy2gp9ej9h88+fP0AKrLlE8XH2CcV56OiJOBjwCfompXbgf+LSJeOLjwpr1DgQXAtVRTJ3ZUdwJdD/yCqj17LdX7d7tGnR2Aq4GnAa8DTgbeBHxmQLFPxrj2u6V+j/1vOkzfN2L7DePb94OAlwJnAYcDf1k//mZEbN+qFBFbAZcDLwTeDLydatafiyM2XoJzoDJz2tyALYHZbWUB3AQsG8f6CZw87P0Y9m0yryOwC1XSfUKj7GnAg8DiYe+bN2+93IB31J8Pzx6jzkrgC21lXwO+Nez4h/R6XAN8ddixDvA12JMqqTu+fi1eNEbdbagSho81yrYEVgCfH/a+jOoN2KLx92rgU13qXQ+cv5lt/QXwMLBTo+yN9f923rD3tZf9btS5lmoO7U3qjtJ+j3ffgZ2oZ/9rlP1OvU/HNcqOqcv2apS9pC47pMT+TKue8MzckJlr28qSqkf32cOJavRM8nU8iKpx+Xxj3Z8Dl1B9I5VG0U/q+6d0WhgRu1MlZRe0LTof2D8idhpgbMMw5usxQ3wS+DTwg3HUfQnwdBrvj8zcAFwIHFas122aycwnNlcnIuZR/QJz5maqHgZ8PTMfbJR9EVjPFPtFfDz73RIRbwJ2Bz7apcrI7DeMb98z88E6Z2n6HrCBJ+cwhwG3Zub3G+teT5XcF8lXplUS3klEbEH1ATjeq3ktjIhf1WODLoyI5wwwvJExgddxHnBfZv60rXwF8N/q7UhTXlTnRmwbEftSDa26JDPv7lJ9Xn3ffnysoPoVaeTHQU/w9Wh5aT0m89F6XPnvFwh14OqhR3sDi8a5Suv9saqtfAXVL4W79Ck0berF9f3TI+KWiHi8Hpf/p2315tF2/GbmeuAONv7/RkpEPA34W+B9mfnLLtWm3X53cQBVB2FzXzfZ99oKCu37TEiI/hT4b1TjGDfnX4A/Bl5JNSb694FvRMTswYU3Msb7Os6m8zittVS9Ztt3WCZNRXcDj1ANw1oDHDtG3dZnRPt7v/WLUrETfQZoIq8HVD+BnwgcAhxHNf726xFxwCCDHLSI2I7qc/CUzBzvJbJnA+tz06tKTqf3x1TVOhfpPKpfHl5Fdd7GJ+te4pax2q5R/f/8FXB7Zl44Rp3puN9PEhFPAf6e6lerrzYWDX3fx3VCyTDVZ+DvPI6qd9Xf3prrvhT4GPB3mXnd5jaQmcc1Hl4XEd8AllONf/zY+KOeekq+jlTjqTYJYYxl0lR0GNWXxr2oTsC7JCJeVQ8j6Kb9/T2d3vcTej0y84zm44j4KvD9et0p9zP3BJwG3Ed1Eu5E+Lm4GZNpp8bQ6mz8bGb+Tf33sojYg+qEvfMadbv9jwb6/xnEfkfEfKoJEl68uboMab9hYP/zdp8CXgD8fmY+3rZsaPsOI5CEU53F/M/jqLcP1Zhl4NdnQ18MfJnqxIMJy8xbI+IHwH/vZf0pptTruJaNvYJNOwC/ojoBRJryMvPW+s/rI2I58B2q4+hfO1Rv9WjOpkrQWnZoWz6yJvh6dFr/4Yi4FDhqQCEOXEQ8F3gv1X7Pqodyt37d2z4its/MX3RYdS2wbURsm5mPNsqnzfujT3pqpzajNTTy6rbyq6nG4z8lM3/F2G3XeIez9moQ+72EagaY1fUMKFB9Idm6fryuHl89zP2Gwez7r0XEGcDbgD/MzO+0LR5r34sck1N+OEpmnpOZMY5bM3Hcg2pWguXAmzsM0J+IaXHCTMHXcSXwzA7zbM4HfjCRE0qkKeQWqpN6ntdleauxah9HOJ+qR6V9LPCo29zr0c2of57uDmwNXErVSK+lOukcYBnw9S7rjfX++DnVFLAzXi/t1Dh0SyQDeIKNPZ4rafv/RMQ2wB5jbKMvBrTfewL/i43v07XAblS/7K8Fnl/XG9p+w8D2HYCIOIFqSM4JmfmVDlU22ffafMp8AZn6SfhE1XNRXwHcC7wmMx+bxLYWUL1Rb+xTeCNjEq/jFVQfbK9vbGt74NVUDZc0ilon9dzZaWFm3kWVaB/TtuhY4Ia2mQemgzFfj04i4jeoZhwY5c/TW4CXt93+vF72LuCELutdD/yMxvsjIrak+py8bJIdRRrb9VRJ5/9sK38lsKIxPOEy4JUR8YxGnddSTS952cCj7L83sOl79T6qX7VfDvywrjfd9huAiHgD1Yw4p2fmWV2qXQa8sJ5Bp7Xei4G5FMpXRmE4yrjVV3i6HHgmcBIwvznzU2Z+q1H3duDuzHxl/fhk4LeoTia6n2r80F8CPwLOLrQLU8JkXsfM/HFEfBr4aEQ8TnUy18l19b8vswdS7yLiS1RDLW6lOhFxb+D99eMv13U+SzXfbPMz9HSqi/PcAVxJdaGqg6hOTBxZvbweEfF7VMf9RVSfAc+mGsYxBxjZq4hm5kNU85//WuOz8abMXF6XXQU8NzOfV6+3PiL+GvhIRDxA9evi26nanDeUiX76qYcH7Vc/3A7Yo565hsz81/r+sYhYBHwsIh4Cvk3VKXQ4VbLZ8k9UExBcHBEfomr/lgDnZWaRXtHxGud+f6vDeo8CP87MaxrFI7PfML59r89j+xfg36kuMtgcF/9AZt5R//1Fqs+xf42IU6hy4r8FvkE1CmDwcgpMvt6vG9W3l+x2a6u7Grim8fjVwDepxo/9Cvgv4LPAzsPer1F6HeuyrYHFVL3oj1A1WnsPe7+8eRvPDVgI3Ex1dblfALdRTUU3q1HnnPZjoS4/juoM/PVUJyEePez9GcbrQTVM5XKqWVQeo+qJvBTYf9j7M4DX52W0Xayn/sxb3VYvgPdRfSl5lOoXgZcPO/5RvgFvHU87Vdc9keqXm8eohhq8pUOd51MlXw8DD1DNBf/UYe/nZPa7bb3VdL64zUjs93j3nWoISrcc5py27e0MfKH+fPsZ8P8Bzyi1P1EHIUmSJKmQaTcmXJIkSZrqTMIlSZKkwkzCJUmSpMJMwiVJkqTCTMIlSZKkwkzCJUmSpMJMwiVJkqTCTMIlSZKkwkzCJUmSpMJMwiVJkqTCTMIlSZJGWER8LyLO6VD+lxFxf0Q8bQhhaTNMwmeYiHhaRDwREa9rK393RPw8ImJYsUmamiLixoh4T1vZoRGREfHMLuucFhFXlolQmvFuA+Y1CyJiNnAy8JHM/PlQotKYTMJnngVAAN/tUH5rZmb5kCRNVRHxh8Bzgc+0LVoA/Dgz7++y6qeA342IVwwyPklAhyQceD/wc+Afy4ej8TAJn3kWUB2Ud3Qob0/MJenPgPMz85G28gXA8m4rZeZDwEXAiQOMTVLl+8DTImJXgIh4FvCnwF9l5vqhRqauTMJnnk16vCNiS2Av6iQ8In4rIr4/pPikGaXXIWIRsUVd55aIeCQiHoiIr0TEjvXybSLigxFxZ0T8MiK+06lXOiJeGRHLIuInEfGLemzpm+tlvwX8HvCvHULYB/he/Rw/joifRsTZEbFdo85S4LCI+M3eXh1J43Rbfd/qDf9L4B7g3E6V63ZfQ7bVsANQcQuAb7WV7Qlsy8ae8H2BG0sGJc1gEx4iVifm/z9wBPC3wDeBZwDHAb+IiG2AK4DdgEVUjfGbgcsiYu/M/EG9nf8BXAZ8HPgw8BRgP+Dx+qleCfyKts+DiPgNYI/6+S6qt/2iehs/Af6irvofVO3My6gSckmDcSfwCDAvIn4AHA+8KTM3tCpExD9QfdbMAxJwqNiQmYTPIBGxFVWP9z+1LToQeAL4Xv34vwPrI2IZMBf4YmaeXCpOaYYZa4jYDV3W+WPgD4EDM/OmRvkFABGxCHgB8DuZ+eO67BrgEOBNwOl1/TcB12TmwsY2/q3x94uA/+zwc/beVL+kXpiZ763Lro6IF9Zx/QVAZq6NiB8Cv4tJuDQwmflERKykSrAXULXnX2yrtgD4IXBwZj5WOER14HCUmWUesA0be7laifk7gDsy8+G6eF9gB+BVwAuBN0RE+wkfkvpjPEPEnh8Rt9Z/B3Aq8E9tCXhz3ROAf2wl4ACZ+TjwI2CXRvVfAi+NiIURsVuH2HYGHuhQvg/wKLC4rfw2oH3oyYPAnA7bkNRf3wcOAt4CnNr2mbIF1RfzE03Apw6T8JllAVWP98KIeGNEvAG4nCrRjog4sK63N/CezHw8M38B/B+g4zRkkiat00nRnYaItRLu+VSJ9Je6bO+FVENTOk0PuCvwX43HZwD/L/Be4O6I+PeIeElj+bZAp5O6FgA3ZWZ7gv4sqqEvTY8CT+0Sq6T+uQ34LeDfM/OKtmW/Ddw9xmxGGgKT8JllAfCfwFlUDe8ngKuAj1H1VL0wIp4DPJKZ90F1chdVD/r3Om5RUs8aQ8Tak/BOQ8RaM5HsXN+v6bLZVk/0vW3PNY8qSb6uVZaZv8jMk6i+ZL8MeDpP/gn7p1S/irVbANzdtv0tgSPZNPmfTTVOXNIAZebHMjMy8+UdFi8Abi4dk8ZmEj6zLABuzsy/y8wdMvNZmfk3mXl6Zj4tMz9N1dj/Zj29EVQndV2YmT8dWtTS9DWRIWKtnvBWcr1Xl23+sL7fo638r4EfAFe3r5CV66h615vnCv0A2L1Zt062X9Bh+2+mOhH0U426WwDPqbcjaXj2wSR8yvHEzJllb+Cjm6mzL/APwJfqqc6uAN436MCkGao5ROzR+u+3Uw0p+VFEHJiZ/0F17LZ6y1cCtwCfjIgdgNVUCfHvZuY7qH7t+jbw8YjYlmrc9zuoxor+fmu2hIg4C9gALAPup5oV5X08+TPiP4DTI2LXzGwNM5lHNUzlWRHxUaohbS8BTgPen5nNE0znA78B/PtkXyhJvWs7+VpThEn4DFFP4P8MNvNNODM/UCYiSWwcIvZZqiFi64G/B64H/pxqiNga4L5Wr3hmboiII6mmJvwbYBZwV70+mZn1nONnAmfXz3MtsF9mrmo890rgDfVtqzqOEzKzOa/wNVQnVh7KxitmLqDquX8V8Dmqi/HcCRyfmf/Stn+HUQ1bccpTSWoTXqVckoYjIq6iSrDfOEado4EjMvPN5SJ70vP/P8A+mTnhOYUjYjnw5cxc1P/IJGm0OSZckoZnbzY/TnNfxrg8fAF/C/xuROwzkZUi4qVU1xn45CCCkqRR53AUSRqCCQwRO6VMRF2f/96IeCubzv+9ObOAt2TmQ/2PSpJGn8NRJEmSpMIcjiJJkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQVZhIuSZIkFWYSLkmSJBVmEi5JkiQVZhIuSZIkFfZ/AY1GpUxCRc67AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 880x440 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(dpi=110, figsize=(8,4))\n",
    "\n",
    "plt.subplot(1, 3, 1)\n",
    "plt.hist(converted_vels[:,0])\n",
    "plt.xlabel(r\"$\\mu_b$\")\n",
    "\n",
    "plt.subplot(1, 3, 2)\n",
    "plt.hist(converted_vels[:,1])\n",
    "plt.xlabel(r\"$\\mu_l cos(b)$\")\n",
    "\n",
    "plt.subplot(1, 3, 3)\n",
    "plt.hist(converted_vels[:,2])\n",
    "plt.xlabel(r\"$v_r$\")\n",
    "\n",
    "plt.subplots_adjust(wspace=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3-SMJ Spark",
   "language": "python",
   "name": "spark-smj"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}