hipparcos proper motion uncertainties

hipparcos_uncertainties.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "85b21bd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import Table\n",
    "from astropy.coordinates import SkyCoord, Distance\n",
    "from astropy.time import Time\n",
    "import astropy.units as u\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4bb527ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "hipparcos = Table.read('hipparcos.fits')\n",
    "now = Time.now()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4146e674",
   "metadata": {},
   "outputs": [],
   "source": [
    "hipparcos = hipparcos[(hipparcos['Vmag'] < 6) & (hipparcos['Plx'] > 0)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4cfd6e42",
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_coord(catalog):\n",
    "    coord = SkyCoord(\n",
    "        ra=catalog['RAICRS'].quantity,\n",
    "        dec=catalog['DEICRS'].quantity,\n",
    "        distance=Distance(parallax=catalog['Plx'].quantity),\n",
    "        pm_ra_cosdec=catalog['pmRA'].quantity,\n",
    "        pm_dec=catalog['pmDE'].quantity,\n",
    "        obstime=\"J1991.25\",\n",
    "    )\n",
    "    return coord"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "68731fcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_position = to_coord(hipparcos)\n",
    "position_today = catalog_position.apply_space_motion(now)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c3e5b8db",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAATqUlEQVR4nO3dcbSkdX3f8ffHXcEABtRdEwKsF1203dgUkxurkTSohC7qQpt4ErbaQyJlDzlHTdpjI9TUpM05OZBa29po6KqEnkihhBjC6iaYqJQYQVkQEESSlaBsSArUk01rTQH99o/nuez0Onede2funbk/3q9z7tmZ3zPPM987s/O5z/ye5/n9UlVIktrytGkXIEmaPMNdkhpkuEtSgwx3SWqQ4S5JDdo47QIANm3aVHNzc9MuQ5LWldtuu+3Rqto8bNlMhPvc3Bz79u2bdhmStK4k+fJSy6baLZNkR5LdBw8enGYZktScqYZ7Ve2pql3HHnvsNMuQpOZ4QFWSGmS4S1KDDHdJapDhLkkNMtwlqUGeCilJDZrqRUxVtQfYMz8/f8FKtzF30UefvP3AJa+dRFmStO7ZLSNJDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkFeoSlKDnKxDkhpkt4wkNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGjTxcE/yt5NcluTaJD8z6e1Lkr69kcI9yeVJHk5y96L27UnuS7I/yUUAVXVvVV0I/AQwP/mSJUnfzqh77lcA2wcbkmwA3gucBWwDdibZ1i87G/gU8PGJVSpJGtlI4V5VNwFfXdT8UmB/Vd1fVY8BVwPn9I+/vqp+CHjDJIuVJI1m4xjrngA8OHD/APD3kpwO/BhwJLB3qZWT7AJ2AWzZsmWMMg6Zu+ijT95+4JLXTmSbkrQejRPuGdJWVXUjcOO3W7mqdgO7Aebn52uMOiRJi4xztswB4KSB+ycCDy1nA86hKkmrY5xwvxU4JcnJSY4AzgWuX84GnENVklbHqKdCXgXcDLwoyYEk51fVE8CbgRuAe4Frquqe1StVkjSqkfrcq2rnEu17OcxB028nyQ5gx9atW1e6CUnSEFMdfsBuGUlaHY4tI0kNmmq4e7aMJK0Ou2UkqUF2y0hSgwx3SWqQfe6S1CD73CWpQXbLSFKDDHdJapDhLkkN8oCqJDVonMk6xlZVe4A98/PzF0x6287KJOmpzG4ZSWqQ4S5JDTLcJalBHlCVpAZ5haokNchuGUlqkOEuSQ0y3CWpQYa7JDXIcJekBnkqpCQ1yFMhJalBdstIUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktSgqc6hulacT1XSU4177pLUIIcfkKQGOfyAJDXIbhlJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBqxLuSf5hkvcn+d0kZ67Gc0iSljbyeO5JLgdeBzxcVS8eaN8O/EdgA/CBqrqkqq4DrkvyLOBdwMcmWvUYHNtd0lPBcvbcrwC2DzYk2QC8FzgL2AbsTLJt4CG/0C+XJK2hkcO9qm4Cvrqo+aXA/qq6v6oeA64GzknnUuD3qur2yZUrSRrFuH3uJwAPDtw/0Le9BTgDeH2SC4etmGRXkn1J9j3yyCNjliFJGjTuHKoZ0lZV9R7gPYdbsap2A7sB5ufna8w6JEkDxt1zPwCcNHD/ROChUVd2mj1JWh3jhvutwClJTk5yBHAucP2oKzvNniStjpHDPclVwM3Ai5IcSHJ+VT0BvBm4AbgXuKaq7lmdUiVJoxq5z72qdi7RvhfYu5InT7ID2LF169aVrC5JWsJUhx+wW0aSVodjy0hSgwx3SWrQuOe5j2Xafe6D48yAY81Iaod97pLUILtlJKlBUw13r1CVpNVht4wkNchuGUlq0FTPlpk1ztIkqRXuuUtSgzygKkkN8oCqJDXIbhlJapDhLkkNMtwlqUEeUJWkBnlAVZIaZLeMJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDnuUtSgzzPXZIaZLeMJDXIcJekBhnuktQg51BdgvOpSlrP3HOXpAa55z4C9+IlrTfuuUtSgwx3SWqQV6hKUoOm2udeVXuAPfPz8xdMs46Vsi9e0qyyW0aSGmS4S1KDDHdJapDnuS/TYD+7JM0q99wlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDZp4uCd5fpIPJrl20tuWJI1mpHBPcnmSh5Pcvah9e5L7kuxPchFAVd1fVeevRrGSpNGMuud+BbB9sCHJBuC9wFnANmBnkm0TrU6StCIjjS1TVTclmVvU/FJgf1XdD5DkauAc4AujbDPJLmAXwJYtW0atd11z/HdJa2WcPvcTgAcH7h8ATkjynCSXAS9JcvFSK1fV7qqar6r5zZs3j1GGJGmxcUaFzJC2qqr/CVw40gaSHcCOrVu3jlGGJGmxcfbcDwAnDdw/EXhoORuoqj1VtevYY48dowxJ0mLjhPutwClJTk5yBHAucP1kypIkjWOkbpkkVwGnA5uSHAB+sao+mOTNwA3ABuDyqrpnOU/eUreMB0slzZJRz5bZuUT7XmDvSp+8qvYAe+bn5y9Y6TYkSd/K4QckqUGGuyQ1aKoTZLfU5z5olP53++glraap7rl7KqQkrQ67ZSSpQXbLrLLB7hdJWit2y0hSg+yWkaQGGe6S1CDDXZIa5AHVRnkevfTU5gFVSWqQ3TKS1CDDXZIaZLhLUoMMd0lqkGfLzDDPeJG0Up4tI0kNsltGkhpkuEtSgwx3SWqQ4S5JDTLcJalBngq5Tiw1o9PgKZKTmvVp8Xac5FtafzwVUpIaZLeMJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUFeoToD1vJKz/V0Velya11Pv9ss83Vsg1eoSlKD7JaRpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMmPipkkqOB9wGPATdW1ZWTfg5J0uGNtOee5PIkDye5e1H79iT3Jdmf5KK++ceAa6vqAuDsCdcrSRrBqN0yVwDbBxuSbADeC5wFbAN2JtkGnAg82D/sG5MpU5K0HCN1y1TVTUnmFjW/FNhfVfcDJLkaOAc4QBfwd3CYPx5JdgG7ALZs2bLcutUbnFhhuY8ZZd3Fj1tq8oalHjNO+3Jra8GsTZQxyXpa/t1m8bnHOaB6Aof20KEL9ROADwM/nuTXgT1LrVxVu6tqvqrmN2/ePEYZkqTFxjmgmiFtVVVfA356jO1KksY0zp77AeCkgfsnAg8tZwNJdiTZffDgwTHKkCQtNk643wqckuTkJEcA5wLXL2cDzqEqSatj1FMhrwJuBl6U5ECS86vqCeDNwA3AvcA1VXXPcp7cPXdJWh2jni2zc4n2vcDelT55Ve0B9szPz1+w0m1Ikr6Vww9IUoMMd0lq0FTD3T53SVodqapp10CSR4Avr3D1TcCjEyxnkma1NutavlmtzbqWb1ZrW0ldz6uqoVeBzkS4jyPJvqqan3Ydw8xqbda1fLNam3Ut36zWNum67HOXpAYZ7pLUoBbCffe0CziMWa3NupZvVmuzruWb1domWte673OXJH2rFvbcJUmLGO6S1KB1He5LzOE6jTpOSvLJJPcmuSfJz/btz07yB0n+tP/3WVOqb0OSzyX5yIzVdVySa5N8sX/tXj4LtSX5Z/37eHeSq5I8Yxp1DZu7+HB1JLm4/yzcl+QfTKG2f9u/l3cl+Z0kx611bUvN99wve1uSSrJpres6XG1J3tI//z1JfnVitVXVuvwBNgBfAp4PHAHcCWybUi3HA9/f334m8Cd088r+KnBR334RcOmU6vvnwH8FPtLfn5W6/gvwT/vbRwDHTbs2utnE/gz4jv7+NcBPTaMu4O8D3w/cPdA2tI7+/9udwJHAyf1nY8Ma13YmsLG/fek0ahtWV99+Et0Itl8GNs3Qa/ZK4A+BI/v7z51UbWv2oVmFF+rlwA0D9y8GLp52XX0tvwv8KHAfcHzfdjxw3xRqORH4OPCqgXCfhbq+sw/RLGqfam0cmj7y2XSjpn6kD62p1AXMLQqDoXUs/v/fB9nL17K2Rcv+EXDlNGobVhdwLfB3gQcGwn3qrxndzsMZQx43dm3ruVtmqTlcp6qfSPwlwGeA76qqvwDo/33uFEr6D8DPA98caJuFup4PPAL8Rt9l9IEkR0+7tqr6c+BdwFeAvwAOVtXHpl3XgKXqmLXPw5uA3+tvT7W2JGcDf15Vdy5aNAuv2QuBH07ymST/PckPTqq29RzuQ+dwXfMqBiQ5Bvht4Oeq6q+nWUtfz+uAh6vqtmnXMsRGuq+ov15VLwG+RtfNMFV9H/Y5dF+Fvwc4Oskbp1vVSGbm85DkHcATwJULTUMetia1JTkKeAfwzmGLh7St9Wu2EXgW8DLgXwDXJAkTqG09h/vYc7hOUpKn0wX7lVX14b75fyQ5vl9+PPDwGpf1CuDsJA8AVwOvSvKhGagLuvfvQFV9pr9/LV3YT7u2M4A/q6pHqupx4MPAD81AXQuWqmMmPg9JzgNeB7yh+v6EKdf2Aro/1Hf2n4MTgduTfPeU61pwAPhwdT5L9w170yRqW8/hPvYcrpPS/6X9IHBvVb17YNH1wHn97fPo+uLXTFVdXFUnVtUc3evziap647Tr6mv7S+DBJC/qm14NfGEGavsK8LIkR/Xv66vpppGcdl0LlqrjeuDcJEcmORk4BfjsWhaWZDvwduDsqvo/A4umVltVfb6qnltVc/3n4ADdyQ9/Oc26BlxHdzyMJC+kO7Hg0YnUtpoHD1b7B3gN3ZkpXwLeMcU6TqP7ynQXcEf/8xrgOXQHM/+0//fZU6zxdA4dUJ2JuoBTgX3963Yd3dfTqdcG/Gvgi8DdwG/SnbGw5nUBV9H1+z9OF0rnH64Ouu6HL9EddD1rCrXtp+snXvgMXLbWtQ2ra9HyB+gPqM7Ia3YE8KH+/9rtwKsmVZvDD0hSg9Zzt4wkaQmGuyQ1yHCXpAYZ7pLUIMNdkhpkuGtZknwjyR39iIm/1V8BuJz1vyfJtf3tU5O8ZmDZ2ZnQ6J5JvjvJ1Um+lOQLSfYmeWGS09OPjjkL0o2E+f5p16H2GO5arq9X1alV9WLgMeDC5axcVQ9V1ev7u6fSXQ+wsOz6qrpk3AL7i49+B7ixql5QVduAfwl817jbXgXbgd9f7krp+PnVkvzPoXH8EbC1H2P8un4c71uSfB9Akh/p9/Lv6AcHe2aSuX6v/wjg3wA/2S//ySQ/leTX+nWfl+Tj/TY/nmRL335Fkvck+XSS+5O8fkhdrwQer6rLFhqq6o6q+qP+7jE5NI78lf0fA5K8M8mtfX27B9pvTHJpks8m+ZMkP9y3H5Xkmr7G/9YP/jTfLzszyc1Jbu+/4RyzxGv4arohX5+U5Jj+d749yeeTnNO3z6Ub9/59dBe8nJTk5/vH3Jnkkv5xb+2/rdyV5Oq+7eh044nf2r8XC9vckORd/TbuSvKW0d9+zbTVvsrOn7Z+gP/d/7uR7tL3nwH+E/CLffurgDv623uAV/S3j+nXmaMf8pRunPRfG9j2k/f7dc/rb78JuK6/fQXwW3Q7JtuA/UNqfCvw75eo/3TgIN1YHU8DbgZO65cNXu35m8CO/vaNwL/rb78G+MP+9tuA/9zffjHdYFnzdGOD3AQc3S97O/DOIbVsAj45pH0j8J0Dj9lPN5DUHN3YIy/rl50FfBo4arB+ujFIFsYHP67/91eANy600V3ZfXT//v02h8Zhn9pV1P5M9sc9dy3XdyS5g27YgK/QjalzGl0YUlWfAJ6T5Fjgj4F3J3krXcg8sYzneTndBCP02z5tYNl1VfXNqvoCK+tq+WxVHaiqb9JdJj/Xt7+y3/v+PN0fqe8dWGdhMLjbBh5/Gt2AbFTV3XTDKEA3wt824I/71+o84HlD6jgT+NiQ9gC/kuQuur36Ezj0e365qm7pb58B/Eb147hU1Vf79ruAK9ONZrnwmp8JXNTXcyPwDGBLv43LFt6bgW1onds47QK07ny9qk4dbFjovlikquqSJB+l29u9JckZwN+s8HkHx8n4v4NPP+Sx9wDDumuGrf8NYGOSZwDvA+ar6sEkv0QXgIvX+QaHPjfDnnuh/Q+qaudhaoBuz/vdQ9rfAGwGfqCqHk83muFCLV9b9DzDxg95Ld2sP2cD/yrJ9/aP/fGquu//K7R77xyDpEHuuWsSbqILJJKcDjxaVX+d5AXVjcp3Kd2e/t9atN7/opuWcJhP041kSb/tTy2jnk8ARya5YKEhyQ8m+ZHDrLMQno/2/eOH++Ow4FPAT/Tb3wb8nb79FuAVSbb2y45KN+Lfk/pQ/T66bw6LHUs3Dv/jSV7J8L1+6Pb635T+jKX+2MfTgJOq6pN0k7QcR9cldgPwloHjCC8Z2MaFSTYubGOE31vrgOGuSfglYL7vRriEQ0PS/lx/cPJO4OscmplnwSeBbQsHVBcteyvw0/02/wnws6MWU1VFN83bj6Y7FfKevsYlx8Ouqr8C3g98nm6EyltHeKr3AZv7Gt9O1x1ysKoeoTt+cFW/7Ba+9Q/bDwCf62td7Eq613Mf3R+2Ly5R8+/TDQ27r+9ueRvd3MIf6ruWPkd37OGvgF8Gng7clW6C5l/uN/MBuu61u/r36R+P8HtrHXBUSGmFkmwAnl5Vf5PkBXRD8L6wqh4bYd1foDsYfPVq16mnJsNdWqEkz6T79vF0uj7tt1fV4m8n0lQY7pLUIPvcJalBhrskNchwl6QGGe6S1CDDXZIa9P8AZYDSmB8lls4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(position_today.separation(catalog_position).to_value(u.arcsec), bins=100)\n",
    "plt.yscale('log')\n",
    "plt.xlabel('Position Change / arcsec')\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ff466b19",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "01d555b9eda7461caa7048b72860a994",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/1000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXEElEQVR4nO3df7RldXnf8ffHkWuMPzA6JDED4xCH0kxsFsotRE0sWhsHdcQaV2SaajSUKa6FMSsrXZLWVm1titWaH0pCJkqnJgYkinSmjsVIxDEFDQMFZaDYCdFyCyuIuiYabRR8+sfZF4+He+7d955z5py75/1aa9bs8z17f/ez7773ud/7nL2/O1WFJKlbHjHtACRJ42dyl6QOMrlLUgeZ3CWpg0zuktRBj5x2AAAbN26sLVu2TDsMSVpXbrrppvur6oSl3ptqck+yA9ixdetWDh48OM1QJGndSfLFYe9NtSxTVfuqatfxxx8/zTAkqXOsuUtSB001uSfZkWT3kSNHphmGJHWOZRlJ6iDLMpLUQZZlJKmDLMtIUgdZlpGkDpqZm5jWastFH3lo+QsXv2gMUUnS+mdZRpI6yLKMJHWQyV2SOsjkLkkdZHKXpA7yJiZJ6iCvlpGkDrIsI0kdNBOP2RsXb2iSpB5H7pLUQSZ3SeqgsZdlkjwC+HfA44GDVfVfxr0PSdLyWo3ck1yW5L4ktw20b09yZ5LDSS5qms8BNgHfBhbGG64kqY22ZZk9wPb+hiQbgEuAs4FtwM4k24BTgRuq6leA144vVElSW62Se1UdAL4y0HwGcLiq7qqqbwFX0Bu1LwBfbdZ5cFifSXYlOZjk4Je+9KXVRy5JGmqUD1Q3AXf3vV5o2q4CXpDkXcCBYRtX1W7gLcDNc3NzI4QhSRo0ygeqWaKtquobwHltOqiqfcC++fn580eIQ5I0YJSR+wJwUt/rE4F7VtOBc8tI0mSMktxvBE5JcnKSOeBcYO9qOnBuGUmajLaXQl4O3ACcmmQhyXlV9QBwIXANcAdwZVUdWs3OHblL0mS0qrlX1c4h7fuB/WvduTV3SZoM53OXpA5yPndJ6iBH7pLUQY7cJamDnPJXkjrIsowkdZBlGUnqIMsyktRBJndJ6iBr7pLUQdbcJamDxv6A7Fmx5aKPPLT8hYtfNMVIJOnos+YuSR1kcpekDvIDVUnqID9QlaQOsiwjSR1kcpekDjK5S1IHmdwlqYPGntyTnJXkU0kuTXLWuPuXJK2sVXJPclmS+5LcNtC+PcmdSQ4nuahpLuDrwPcBC+MNV5LURtuR+x5ge39Dkg3AJcDZwDZgZ5JtwKeq6mzgDcBbxheqJKmtVsm9qg4AXxloPgM4XFV3VdW3gCuAc6rqO837XwUeNbZIJUmtjTJx2Cbg7r7XC8CZSV4GvAB4AvDuYRsn2QXsAti8efMIYUiSBo2S3LNEW1XVVcBVK21cVbuT3AvsmJubO32EOCRJA0a5WmYBOKnv9YnAPaOFI0kah1GS+43AKUlOTjIHnAvsXU0Hzi0jSZPR9lLIy4EbgFOTLCQ5r6oeAC4ErgHuAK6sqkOr2bmzQkrSZLSquVfVziHt+4H9a915Ve0D9s3Pz5+/1j4kSQ/nfO6S1EFTfYbq0Rq5+zxVSccaJw6TpA6yLCNJHeRj9iSpgyzLSFIHWZaRpA6yLCNJHWRZRpI6yLKMJHWQZRlJ6iDLMpLUQSZ3Seogk7skdZDJXZI6yKtlJKmDjokpf/s5/a+kY4FlGUnqIJO7JHWQyV2SOsjkLkkdNJHknuQxSW5K8uJJ9C9JWl6r5J7ksiT3JbltoH17kjuTHE5yUd9bbwCuHGegkqT22o7c9wDb+xuSbAAuAc4GtgE7k2xL8nzgduCvxhinJGkVWl3nXlUHkmwZaD4DOFxVdwEkuQI4B3gs8Bh6Cf+bSfZX1XcG+0yyC9gFsHnz5jUfgCTp4Ua5iWkTcHff6wXgzKq6ECDJq4H7l0rsAFW1O8m9wI65ubnTR4hDkjRglOSeJdrqoYWqPSt1MI07VPv1360K3rEqqTtGuVpmATip7/WJwD2r6cC5ZSRpMkZJ7jcCpyQ5OckccC6wdzUd+CQmSZqMtpdCXg7cAJyaZCHJeVX1AHAhcA1wB3BlVR1azc4duUvSZLS9WmbnkPb9wP617nzaNXdJ6irnc5ekDppqcrfmLkmT4chdkjrIkbskddBUH7M3a3wEn6SusCwjSR1kWUaSOsgnMUlSB5ncJamDrLlLUgdN9WqZWZ5+wCtnJK1nlmUkqYNM7pLUQSZ3SeogP1CVpA7yJiZJ6iDnlmnBK2ckrTfW3CWpg0zuktRBJndJ6qCx19yT/BjwemAjcG1V/e649zFN1t8lrQetRu5JLktyX5LbBtq3J7kzyeEkFwFU1R1VdQHwc8D8+EOWJK2kbVlmD7C9vyHJBuAS4GxgG7AzybbmvZcAfwZcO7ZIJUmttSrLVNWBJFsGms8ADlfVXQBJrgDOAW6vqr3A3iQfAf5oqT6T7AJ2AWzevHlt0U+ZJRpJs2qUmvsm4O6+1wvAmUnOAl4GPArYP2zjqtoN7AaYn5+vEeKQJA0YJblnibaqquuA61p1kOwAdmzdunWEMCRJg0ZJ7gvASX2vTwTuGS2c9csSjaRZMsp17jcCpyQ5OckccC6wdzUdOLeMJE1G20shLwduAE5NspDkvKp6ALgQuAa4A7iyqg6tZufOCilJk5Gq6X+WOT8/XwcPHlzTtv3lkFlkiUbSpCS5qaqWvJ/I+dwlqYOcz12SOsj53CfMq2gkTcNUk/uxdp27iV7S0TLV5F5V+4B98/Pz508zjmkw0UuaJOdzl6QOsiwzAxzFSxo3r5aRpA7yapkZ0+amLEf3klZiWWYdsowjaSVeLdMhJn1JiyzLrHOzPreOpOnwUkhJ6iBH7h3lB7PSsc3kfgyzRi91l1fLCFh+pG/il9Yfb2KSpA6yLKMVDRvVO6KXZpdXy0hSBzly15o5opdml8ldY+dVONL0TSS5J3kp8CLgB4FLqupjk9iPZp+JXpqO1sk9yWXAi4H7quppfe3bgd8CNgDvqaqLq+pq4OokPwC8AzC5a9X8xSCt3WpG7nuAdwPvW2xIsgG4BPhHwAJwY5K9VXV7s8obm/elocnaJC6NX+vkXlUHkmwZaD4DOFxVdwEkuQI4J8kdwMXAR6vq5qX6S7IL2AWwefPmNYSu9WzYh7FOhCaNx6g1903A3X2vF4AzgdcBzweOT7K1qi4d3LCqdie5F9gxNzd3+ohxSJL6jHqde5Zoq6r67ao6vaouWCqx963oHaqSNAGjjtwXgJP6Xp8I3NN2Y+eWUVteUy+tzqgj9xuBU5KcnGQOOBfY23ZjR+6SNBmruRTycuAsYGOSBeBNVfXeJBcC19C7FPKyqjq0ij4duWskjuilpa3mapmdQ9r3A/vXsnOfoSpJkzHVicOS7Eiy+8iRI9MMQ5I6Z6pzyzhy16R4w5SOdT6JSZ3njVE6Fjly1zGrzSjekb7WK6f8lTCJq3v8QFWSOsiyjNTSYO3eEb5mmWUZaYAfwKoLTO7SGnl3rGaZNXdJ6iBr7tIEtSnxONLXJFiWkcZstTV7L8PUJEy1LCNJmgxH7tI640hfbfiBqiR10FSTu09ikqTJsOYuSR1kzV2aUd4pq1E4cpekDnLkLs0QR+sal7GP3JP8aJL3JvnguPuWJLXTKrknuSzJfUluG2jfnuTOJIeTXARQVXdV1XmTCFaS1E7bkfseYHt/Q5INwCXA2cA2YGeSbWONTpK0Jq1q7lV1IMmWgeYzgMNVdRdAkiuAc4Db2/SZZBewC2Dz5s1t45XUZ5S7Vb3TtdtGqblvAu7ue70AbErypCSXAk9P8mvDNq6q3VU1X1XzJ5xwwghhSJIGjXK1TJZoq6r6MnBBqw6SHcCOrVu3jhCGJGnQKMl9ATip7/WJwD2jhSNprdo8GcpLLY8do5RlbgROSXJykjngXGDvajpwbhlJmoxWI/cklwNnARuTLABvqqr3JrkQuAbYAFxWVYdWs3PLMtJsGNcHs2vZXpPR9mqZnUPa9wP717pzH7MnSZPhfO6S1EHO5y5JHeSskJLUQVOdFdIPVKXZ452r3WBZRpI6yLKMJHWQZRmp40a5K3VSd7Ra+pk8yzKS1EGWZSSpgyzLSJqYNmWdNhOeafUsy0hSB1mWkaQOMrlLUgeZ3CWpg0zuktRBXi0jaawmcePTaq+oaXOT1LB1Rtl22DrLrTcpXi0jSR1kWUaSOsjkLkkdZHKXpA4yuUtSB439apkkjwF+B/gWcF1VvX/c+5AkLa/VyD3JZUnuS3LbQPv2JHcmOZzkoqb5ZcAHq+p84CVjjleS1ELbssweYHt/Q5INwCXA2cA2YGeSbcCJwN3Nag+OJ0xJ0mq0KstU1YEkWwaazwAOV9VdAEmuAM4BFugl+FtY5pdHkl3ALoDNmzevNm5JWvUNU8PWX237sHWWu1HpaE9tPMoHqpv47ggdekl9E3AV8LNJfhfYN2zjqtoNvAW4eW5uboQwJEmDRvlANUu0VVX9DfCaNh1U1T5g3/z8/PkjxCFJGjDKyH0BOKnv9YnAPavpIMmOJLuPHDkyQhiSpEGjJPcbgVOSnJxkDjgX2LuaDpxbRpImo+2lkJcDNwCnJllIcl5VPQBcCFwD3AFcWVWHVrNzR+6SNBltr5bZOaR9P7B/rTu35i5JkzHV6QccuUvSZDifuyR1kCN3SeqgVNW0YyDJl4AvrnHzjcD9YwxnFh0LxwjHxnF6jN0wK8f4lKo6Yak3ZiK5jyLJwaqan3Yck3QsHCMcG8fpMXbDejhG53OXpA4yuUtSB3Uhue+edgBHwbFwjHBsHKfH2A0zf4zrvuYuSXq4LozcJUkDTO6S1EHrJrkPeV5r//tJ8tvN+59N8oxpxDmKFsd4VpIjSW5p/v2bacQ5imHP4+17vwvncaVj7MJ5PCnJJ5LckeRQktcvsc66Ppctj3F2z2VVzfw/YAPwF8CPAnPArcC2gXVeCHyU3kNEfhL4zLTjnsAxngX8t2nHOuJxPgd4BnDbkPfX9XlseYxdOI9PBp7RLD8O+HwHfybbHOPMnsv1MnJ/6HmtVfUtYPF5rf3OAd5XPZ8GnpDkyUc70BG0OcZ1r6oOAF9ZZpX1fh7bHOO6V1X3VtXNzfLX6E37vWlgtXV9Llse48xaL8l92PNaV7vOLGsb/zOT3Jrko0l+/OiEdlSt9/PYVmfOY5ItwNOBzwy81Zlzucwxwoyey1GeoXo0Lfm81jWsM8vaxH8zvbkkvp7khcDVwCmTDuwoW+/nsY3OnMckjwU+BPxyVf314NtLbLLuzuUKxziz53K9jNzbPK915Ge6TtmK8VfVX1fV15vl/cBxSTYevRCPivV+HlfUlfOY5Dh6Se/9VXXVEqus+3O50jHO8rlcL8m9zfNa9wKvaj6h/0ngSFXde7QDHcGKx5jkh5OkWT6D3vn78lGPdLLW+3lcURfOYxP/e4E7quqdQ1Zb1+eyzTHO8rlcF2WZqnogyeLzWjcAl1XVoSQXNO9fSu9xfy8EDgPfAF4zrXjXouUxvhx4bZIHgG8C51bzkf16kd7zeM8CNiZZAN4EHAfdOI/Q6hjX/XkEng28Evhckluatn8JbIbOnMs2xziz59LpBySpg9ZLWUaStAomd0nqIJO7JHWQyV2SOsjkLkkdZHKfYUm2DM4smOTNSX51wvt9aZJtLda7IMmrVljntObOvXHE9agkH29m33vFOPoc6P9HknywWf6euJO8JEvM1DnGfe9J8vJJ9b9aSX4vybOnHYfWzuSu75HkkcBLgRWTe1VdWlXvW2G10+hd6zwOTweOq6rTquoDY+rzIVV1T1UtJtjT6Iu7qvZW1cXj3ucMOxP49Go3ar5/NANM7utYkuuSvC3Jnyf5fJKfbto3JHlHks8182i/rmk/Pcknk9yU5JrFGfqafn49ySeBNwAvAd7ejJCfmuT8JDc2kyN9KMn3N9s99FfEUrE0d9r+W+AVi6PtJP87yQnNNo9Ib67vjQPH9cQkVzexfzrJTyT5QeAPgdMW41ria/GbSa5Pcltzt+CSfTXt/yDfnYP7fyZ53OJfSkPifnWSdzfbPiXJtU2f1ybZ3LTvSW/+8uuT3DVsJJ7kVc22tyb5g763njO4bZLHNvu4uTmf5zTtW9KbZ/z305tr/GNJHt289/eb/m9I8vY0f/013xdvb87lZ5P88yHx/Rjw+ap6cKB9R5LPNF+vjyf5ob7vg91JPga8L8kPJflwc3y3JnlWksck+Ujz+rY0f3kt8z25tdnHrc2xP3UwTq1g2nMO+2/4P2ALA3OCA28GfrVZvg74T83yC4GPN8uvpTcfxiOb10+kd4fk9cAJTdsr6N0Fu9jP7/TtYw/w8r7XT+pbfivwulXE8mrg3X3bv4neBEwAPwN8aInjfhfwpmb5ecAtzfJZDJk7u9n/7zfLz1n8ui3T1z7g2c3yY+ndrb2lb7vBuB963Wz7C83yLwJX933d/pjeoGkbvSmcB+P8ceBOYOPiuVlu2yauxzfLG+nd7Zkm1geA05r3rgT+abN8G/CsZvnivmPaBbyxWX4UcBA4eYkYfwX4xSXaf4Dv3vj4z/rO95uBm4BHN68/0HeONwDHAz+7eH6a9uNZ/nvyM8A/bpa/D/j+af88rrd//gk124bdPtzfvjiZ0U30fuABng9cWlUPAFTVV5I8DXga8CfpTYWxAeif52O5MsfTkrwVeAK9RHjNkPWWimXQZcB/BX6TXmL8z0us81P0kgFV9adJnpTk+GXiW3R5s82BJI9P8oRl+vofwDuTvB+4qqoWmq9LG88EXtYs/wHwH/veu7qqvgPcvjiyHfA84INVdX8T01dW2DbAryd5DvAdelPmLr73l1V1S7N8E7ClOebHVdX1TfsfAS9uln8G+Im+vyiOpzeD4V8OxPgClp4q4ETgA83oem5gu71V9c2+Y3xVc3wPAkeSfA54R5K30fsF/alh35NJHgdsqqoPN338vyVi0QpM7rPty/RGS/2eyPf+UP1t8/+DfPd8hqWnRD5UVc8csq+/WSaOPcBLq+rWJK+mN4JeylKxfI+qujvJXyV5Hr267s8vsdpap4odXKeG9VVVFyf5CL2/Mj6d5PnAWpNI/37/tm95qX0vdW6W2/bngROA06vq20m+QG8kO7j+g8Cjh+yzv8/XVdWwX86kV3J7QlUtNXvju4B3VtXeJGfRG7EvWu77h6r6fJLT6X29/0NTwvkwS3xPJnn8cn2pHWvuM6x6U4nem+QfQq9+DGwH/myFTT8GXJDmw61muzuBE5I8s2k7LsMfLPA1eo8VW/S4Jo7jWDoZL2ewL4D30KufX1kDdd3GgcX9NEnk/nr4PNpLWazj/hS9GQiPDOsryVOr6nNV9TZ65Ym/2yLuRdfTm7WTpu+Vzke/a4GfS/KkJqYnrrD+8cB9TWJ/LvCU5Vauqq8CX0tvFkb64oTeX1yvbc4jSf5OkscMdPFc4BPLxPJ/m+VfWCaMa+mVBhfr/I9P8iPAN6rqD4F30HsM4ZLfk825Xkjy0qb9Uc0vHa2CyX32vQp4Y3qz0v0p8Jaq+osVtnkP8H+Azya5Ffgn1Xt038uBtzVttwDPGrL9FcC/aD44eyrwr+nVQP8E+F+rjP8TwLZ87+WLe+mVd5YqyUBvRDif5LP0asbLJZJ+X01yPXApcN4Kff1y88HerfRm8/toi7gX/RLwmqbPVwIPe3DyMFV1CPj3wCebfQ+bLnfR+5v4D9L7RdLm638esDvJDfRG60ea9vcAtwM3Nx+y/h4P/wvrbOC/D+n3zcAfJ/kUcP8y+3898NymFHMTvc8Z/h7w58338b8C3rrC9+QrgV9qvsbXAz+88mGrn7NC6qhLMg/8RlX99Bj7vI7eh7sHx9XnepXksc1ffaR3bf6Tq6rVL6AkNwNnVtW3JxmjJs+au46qJtm8ltWXd9Tei5L8Gr2f7y/Su9Knlap6xqSC0tHlyF2SOsiauyR1kMldkjrI5C5JHWRyl6QOMrlLUgf9fzyjZsjiWiDxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rng = np.random.default_rng(42)\n",
    "\n",
    "seps = []\n",
    "for i in tqdm(range(1000)):\n",
    "    random_pm = hipparcos.copy()\n",
    "   \n",
    "    unit = hipparcos['pmRA'].unit\n",
    "    random_pm['pmRA'] = rng.normal(\n",
    "        hipparcos['pmRA'].quantity.to_value(unit),\n",
    "        hipparcos['e_pmRA'].quantity.to_value(unit)\n",
    "    ) * unit\n",
    "    random_pm['pmDE'] = rng.normal(\n",
    "        hipparcos['pmDE'].quantity.to_value(unit),\n",
    "        hipparcos['e_pmDE'].quantity.to_value(unit)\n",
    "    ) * unit\n",
    "\n",
    "    random_coord = to_coord(random_pm)\n",
    "    random_coord_today = random_coord.apply_space_motion(now)\n",
    "    seps.append(random_coord_today.separation(position_today))\n",
    "\n",
    "plt.hist(np.concatenate(seps).to_value(u.arcsec), bins=100)\n",
    "plt.xlabel('Uncertainty of position change / arcsec')\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b479fdd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}