{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set 1, Problem 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster Color Magnitude Diagram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this problem is to make a color magnitude diagram for M67, a nearby open cluster, using\n",
    "Gaia data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (a) Determine the location of M67 using public tools."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that you can also just resolve \"M67\" in the Gaia archive search. \n",
    "From the SIMBAD database (http://simbad.u-strasbg.fr/simbad/sim-fid):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Coordinates (ICRS)\n",
    "ra_m67 = '08 51 18.0'\n",
    "dec_m67 = '+11 48 00'\n",
    "\n",
    "# Angular size (arcmin)\n",
    "d_m67 = 25.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (b) Download a list of stars near M67."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download from the Gaia archive at https://gea.esac.esa.int/archive/. Make sure you download all the necessary data for the steps below. You'll need the following: ra, dec, parallax, pmra, pmdec, phot_g_mean_mag, bp_rp, a_g_val, e_bp_min_rp_val (the last two are for extinction corrections)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I've saved all stars within a radius of 12.5 arcmin of M67's coordinates to a file called 'gaia_result.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the file with data\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "data = pd.read_csv('gaia_result.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's clean up the data a little before we plot anything."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\ndata = data[~np.isnan(data['a_g_val'])]\\ndata = data[~np.isnan(data['e_bp_min_rp_val'])]\\n\""
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get rid of any stars that don't have all the measurements we need\n",
    "data = data[~np.isnan(data['parallax'])]\n",
    "data = data[~np.isnan(data['phot_g_mean_mag'])]\n",
    "\n",
    "# In case we want to consider reddening/extinction\n",
    "'''\n",
    "data = data[~np.isnan(data['a_g_val'])]\n",
    "data = data[~np.isnan(data['e_bp_min_rp_val'])]\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (c) Plot a histogram to determine a rough estimate of M67’s distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average: 1.16 +/- 0.52\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'N')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Compute average and standard dev of a quantity (weight by the inverse variance)\n",
    "def compute_avg(values, error):\n",
    "    avg = np.average(np.asarray(values), weights=1./np.power(error,2.))\n",
    "    std = np.sqrt(np.cov(np.asarray(values), aweights=1./np.power(error,2.)))\n",
    "    print('Average: '+\"{0:.2f}\".format(avg)+' +/- '+\"{0:.2f}\".format(std))\n",
    "    return avg, std\n",
    "    \n",
    "parallax_avg, parallax_std = compute_avg(data['parallax'], data['parallax_error'])\n",
    "\n",
    "# Make plot\n",
    "plt.hist(data['parallax'], bins=25)\n",
    "plt.axvline(parallax_avg, color='r', linestyle='--')\n",
    "plt.axvspan(parallax_avg - 2.*parallax_std, parallax_avg + 2.*parallax_std, color='r', alpha=0.2)\n",
    "plt.xlabel('Parallax (milliarcsec)')\n",
    "plt.ylabel('N')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's estimate the distance to M67 based on its parallax, using d = (1 AU)/parallax. (Note that in actual science cases, this is NOT the best way to get distances from Gaia! But it's fine for our purposes here.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Distance: 860.70 pc\n"
     ]
    }
   ],
   "source": [
    "dist_pc = 1./(parallax_avg/1000.)\n",
    "print('Distance: '+\"{0:.2f}\".format(dist_pc)+' pc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's decide how many stars to throw out. Start with a simple +/- 2 sigma (95%) cut:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data1 = data[data['parallax'] > (parallax_avg - 2.*parallax_std)]\n",
    "data1 = data1[data1['parallax'] < (parallax_avg + 2.*parallax_std)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (d) Further select M67 stars by plotting the proper motions of the selected stars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average: -9.91 +/- 3.53\n",
      "Average: -3.16 +/- 3.04\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEOCAYAAACuOOGFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcVOWZ9//PVdXdIAraogZkleAKZqFR2yRO1OBEHUccTaLRR81klCQ/jfFJMonLM5gwj/nFmUzGZHQmIcTXRH8IbrhMRo1LNMvERmniAhJji4DNIiAtoEAvVdfvj3Oq+lR1dXd1Uyv9fb9eRVWdc6rO1dXNfdW9nPs2d0dERGRvxcodgIiI7BuUUEREpCCUUEREpCCUUEREpCCUUEREpCCUUEREpCCUUEREpCCUUEREpCCUUEREpCBqyh1AKR1yyCE+efLkcodRfXbsgK6unts3bgzux47N/bqaGhg1qnhxiUhJNDc3b3X3Q/s7bkgllMmTJ7Ns2bJyh1F9Fi+GMWN6br/22uD+1ltzv27TJrjoouLFJSIlYWZr8zlOTV4iIlIQVVFDMbMJwJ3AGCAJzHf3H5nZwcA9wGRgDfA5d28rV5xDzqWXljsCEakg1VJD6QK+4e7HAo3AVWZ2HHAd8LS7Hwk8HT6XUmloCG4iIlRJQnH3je6+PHy8E1gFjANmA78ID/sFcF55IhyiWlqCm4gIVZJQosxsMvBRYCnwAXffCEHSAQ4rX2RD0G23BTcREaosoZjZAcADwLXuviPP18wxs2VmtmzLli3FDVBEZAirmoRiZrUEyWShuy8JN79tZmPD/WOBzdmvc/f57j7T3Wceemi/w6j3Oc1r27j9mRaa12qsgogUV7WM8jLg58Aqd/9hZNcjwOXA98P7h8sQXsVqXtvGJQua6OhKUlcTY+EVjTRMqi93WCKyj6qWGsrHgUuB083sxfB2NkEiOcPMXgfOCJ9LqGn1O3R0JUk6dHYlaVr9TrlDEpF9WFXUUNz994D1svtTpYylmjROGU1dTYzOriS1NTEap4wu7AmuuKKw7yciVa0qEooMTsOkehZe0UjT6ndonDK68M1d06cX9v1EpKopoezjGibVF6/fZMWK4F6JRUSonj4UqUQLFgQ3ERGUUEREpECUUEREpCCUUPYhuohRRMpJnfJl0ry2raCjr3QRo4iUmxJKGRSj8M91EWPRE8rVVxf3/UWkqqjJqwyKcQV76iLGuFGcixhzmTo1uImIoBpKWRTjCvaiX8SYS3NzePLMRbaa2xI0bUvQSIxiL79V6KZDERk8JZQyKFbhn30RY9EL27vuCk/cnTaa2xJc8sIeOpJQZyNYuLataAW9+o1EKosSSpkU9Qp2ylfYNm1L0JGEJNDpZPTlFDrBlaXfSER6pYSyjypXYdt4cJy6WCedSYgDG97dTfPaNl7btJO5D68gkXSG1RYmwRV98ksRGRAllH1UuQrbhvo4C08YzgPrO7l/fReLnl/HfcveIuGQSDoAHZ2FSXBl6TcSkV4poeyDUk1Lc8+ZRtuujpIXtg31cZq2JejysOkr4XhkfyxmBUtwxW46FJH8KaHsY4rRd9L8foymNzpoPDgOBP0kjQfHafj613t9TePBceoMOoF4zMCMrkSSmBnzZk8v2MWcqp2IVA4llH1MoftOmte2cUnLCDq8kxrrBIOuJNTFOll4wjga6oMkkx4qfHCchvrgtnDqLpqmzEjXRnor/AeTGDTCS6TyKKHsYwrZd9K8to1bn/oz7Q5OMGqL1OMkNL34Jg3DNtN87EndQ4VjnSw8YXiQVPZP0nBa94WPuQr8wSYGjfASqTxKKPuY/jqq860NRAt6J5hSocYAg0QSamPQ+JuH4d23aPrqzO6hwsmgSSxVc8n1vtHzRxNDR1eSW5/6M9fOOqrf5NA4ZTQ18SBxxuM9E6eaw0RKTwllCBlIbSBa0MeA40cZ00fFmTYqRlunB01bj74FZA4Vro2R7muJnrdp9TvUj6hj3i9XZpy/ccpoamJGR8JJOvzu9a28sGZbfjUV98z7QfycIlI4Sij7mL4K0/6aiaLf6qNNZ3Fg1U7nlR1d1MVIN2kBNB80gaZtCeYeU9edaCK1k2g8MTMSyWDEV2dXkgeWt7J1ZzudicyE0N6ZZMny1owaRnaN44HlrenRY4mkZ/wsag4TKQ8llH1MX4VpX/0ruRJRquls/bJXWPxOXY8mreaDJnDJCV+k4/VOamLwmcPjBJcz5o4HnHjMcA/u71v2Vo9kEhwF9y17i66kU1cTY+450zJqNnPPmcb9za3pocjxrGHI/fUj9dYcNtBmMjWriWRSQtnH9FWYNkyqZ+4503hsxUbOmj42oxDMTkRLlrdy+EH7Ba9f3cmStrrg6neDDbs9GNV18BF0WJwk0JGERa0JlmyI1FaI9YgndW3Mhnd3s3Dpupw/w/iDhrNh+550LI+t2JgR22MrNtLZlUwf/9mZEzJ+lux+JIDbn2lJP85VgxtoM1l2zWve7OlcfNLEQf3ORPYVSij7mL465ZvXtqW/6b+wZhtHjxmZs/YSj8cyaggLJ9N99fuGBItau3hgQxdzP3Uqda2x9CgwJ0gsc1d1kHCIM4J503em46kfUZe+0PK1TTt7/RnO+dDh/Odza+joTGJmTBs7ihfWbEvHtqczkXGh5MhhNemEkfp5UvdLlrdm/CwXzBifswY30Gay6PFJd+Y+vCLj8xQZiqo+oZjZmcCPCNpaFrj798scUtn1dvV4X4VmNBGtf3c3i59f133cezVc9cHw6vdkIt301Tb8ABaeGGfJ+i7uW99FV1jKp++BGx98hTOO+wCnHn1YOpnFzDjtmMN6jX/11veZe8609NxfP/v9m8yYeBD1I+p49rXNLFvTvcSxAT/93WrcIWYwc1I9Uz8wkumHH8jcR1bQFWlSa+9Msnlne84a3ECHWzdOGU3MjGQ4ICC7H6cU1OQmlaaqE4qZxYHbgTOAVuAFM3vE3V8tb2SVqb8+lGgT0ZLlrd3HHdAVvD57NNebL9Kw9j346CfZ3J7k11tTfSUphgNPvPo2T/9pc3our6Q7T696u9c4n3r1bQ4dOYykd3e6v7CmjRjB0OQoT/8DSYfn17TxfCThZB/7zGubmXfu9B5T0gx0uHXDpHqu+MQR/OS3q9PvXT+irtefqdA0kk0qUVUnFOBEoMXdVwOY2WJgNqCEkkNvhWZvHfJLlocd39tWBa8PJ35s2pagvtZ44PGd/LRuJM927qHTIXf3ugHdE0OmJHse3L0P2LyzHTPLGBKcnUwGoyvhPPvaZuZfNrPHvt5qds1r2/j8z5rSCXbRlUHhPXK/WmIW/Cwxg7ZdHQWIMD8aySaVqNoTyjjgrcjzVuCkMsVSFXItwjXvv1aypzMorlOFU+OU0elRWPcwgulv7+bkg2OMrDXqa43vrOqgY0L4UfeRHFLiMctIKqmCuDe//tNmhsWNXX0dNEhPvvo2dy9dl3cn+pLlrXSEgwA6wgELDZPqyzp9vqbul0pk7oX/D1sqZvZZ4NPufkX4/FLgRHf/auSYOcAcgImHHdaw9kc/KkuslaT5/RhN79VQH0/yndbhRL9Xx4HpIxLUAS/sigNhLcG6j6nBSWB4emP0b8i6nzscOzzBjAOSTNsvwbz1w+kIL5ScMSLB86n375GR+tsWCabHtlyxRPd1v89xwxPUxuDCgzu5+JAuenPDW8O4+53a9GsvHt3J9ya0A92fZeMBXTTs37MO1d/+vVHM9xaJss9/vtnde1brs1R7DaUVmBB5Ph7YED3A3ecD8wFmTpnijBlTuugqSGryxvpaY94bHXQkg1pCdjGaAF7alfVnYal/ggI5/Zrwy0jcYEysk/XJWjIKeHP+3FHDhFicFcDcY+O0dTr1tcY9rZ0Zx/ZMB9GkMZBtfT23jPtX9wQ/50u7alhXU8PIWutxYSbABcMS3L8taNarNeOCqaMgPKYhvOXS3JbgkpdTc5wNy7ggtBD6OrdIOVR7QnkBONLMjgDWAxcBF5c3pMoTXec9ZpDwzM7s/kW/6fd8UQLCZBI9Njg+4fDElgQAcRJcObmGf3ujk47IF+p4eHQ5v2P/ZE0XBtTGOvnM4XEuGFebLvwb6uMsOnF4xmzK+chYDrmfOc5E9gWxcgewN9y9C7ga+BWwCrjX3VeWN6rKEy3Ykh7UKOJAXQzOG5NvARf9dh/ezCIVgOw6Rs/kkwB+uqaL9jCWGHDK6BhXTq4hlqvyUWKp62gWtQYJuLktkd7XUB/nqg/WDSghBKPigs861xxnIvuaqu5DGaiZ++3ny44+OnPjqafCeefBnj1w3XU9X3TmmcFt+3a46aae+889F04/HTZvhu99r+f+z30OPvYxWLcOfvjDnvsvvRQaGqClBW67ref+K66A6dNhxQpYsKDn/quvhqlTobkZ7rqr5/6vf53mkeO4ZOkuOpNOrSeYu+q/aasbQeO2N2n4fy7l7j313LPqHV7p2o9kLCz0LN8SPtpHkdoU+ZtKvY978DjcF/MkNe6cOnYYz25J0JFMhm9lGcf1iCPrfXrdn/04O7bUe0T3R/clk5yxu5WfXXBssG3uXNixI/O9ZsyAyy4LHn/729Denrn/5JNp/svPBDWbu/6Nhnffytw/BP72mDgR/vAHuPfenvtvuAEOOwx+/Wt45JGe+7/7XTjwQHj88eCW7fvfh+HD4aGH4Nlne+6/9dbg/p574LnnMvcNGwa33BI8vvNOWL48c/+oUTBvXvD4Zz+DlVnfUw89FG68MXh8223BZxg1fjx885vB4x/8AFpbM/dPnRp8fgA33wxbtmTunzYNrrwyeDzIvz0uvDB4fO219DDAvz176aUh0YcieWioj7NwxBs0vbg6SCJZBdvFE2u5uOVVmp/7I1/68MVsHT4y2JGrQO4hq0AOjz9093a27Hdg5qHhfvMkU97fwpv7H8oTm8NagMXo7sz3ngkj9d65EkFGOJZ5bPb5cz0On8dwSCaDpGrGkyMmcPe6Ti6eWMtgpRYbIzuZiOyDhlYNZcoUX3bHHeUOo2I1tyX4/PN76BjUn0RmU9e4WCcbkrV42OGe/Za5LlKMiofDikv11xkD/u9xddzT2sVLO7ojO2V0jLtO2K8g58he1VKkWthpp+VVQ6nqPhTZe81tCW5/oyNd2HX1WoJ71n32vsw+lPXJ2u7ZgMmegzgzmWT/ERphJSWP+KNGxGDcsAG+KBLPyh1JLhyfWWmfNjKW/nyiop9bPlIDI/7l9c4e/TMi+wo1eVW4vf1W29frUzWSYDhsJ985to4aI0cNJZpMIokj4xKQ6Ciw7nsDZtTHeKEtd32kxuC0Q2I8uaV7vxN0YncmB5ZU2pOwvr33/amo4gZ//YE4D2/KnGTS8XTz1j2tnQyLG3esDeYoiy5tHB01F93eF434kqFACaWCZQ737WTesXUDas/vr+B7YH1nOnl0eFCInnpovLtfI7QfCXaTfRFi+NiduIFjJLNqMTGcmpjR3JbMmRjiwBWTatjZBdkNYGOHBdd8tOzKP6X0950/lQrjwP41mf0rceCCccFne/TIGK+957SHi4FBZhIYTHLob1XLXPr7MqDmM6k0Qyuh1NTApk3ljiJvTW/X0ZGsI0kwq+3cV9s5urMt76uio6/vTDpN67bT0N59XbztHgZ0XwH+0g6nbkeSGroveKw16EhapBaSEj6OBZ3pY9nDeoYT1EqcWdbGR9q3sn7sZBZFrjLvFiSg+Ws6s5JNcNza3dHz9XYdTO7rYnoKjkuGjzvd2bxjN7XU0EnQ5PaP4/fQ0L4TNnV/bh6+zsLPoZF3YVOS+t01xBie/nxS2/vSACz8YOTK9vBcvWl+P8YlLSPocKgzWDh1V/r33tc+kXIaWgll1Ci46KJyR5G3xrVtxH76HMlkUJgmgaYpM2g4bWrer69bkJrUME7jeadBZB6v89e2cd/PmsJ5qsJJHA0uPDGY4ypVXC9auja934AzDo3zdruzckcy3bC1nu6O6zjGl088nIb2WppPPpUlC5poD9c2OffDh/P4yk10dCZJhrMR99TXFfH9XRnfm+7jHOOZ9+uwOJBw4jUxjv7rT6U/m+jnFo/H+UzDeC6YMT69ENe8BU0kLZiGf+7s6TScNDGvqeQHcmV70zMtdLS8FtSCyPy997VPpCg+//m8DhtaCaXKNEyqZ97s6cx9eAXJpFNXO7BJAPubkr1hUj2LrmxML0KVSDq14SJU0ZmI71+6Nj3fV43BIcPgS1OCqdqbtiV44u1ExsioaaNiQTPMpu5VIuc+vIKkO4+v3JReNfL3r2/dq1Fc+9fF6Uw6XYlg1t186ysAiUR3c1ZXove1YbI/t+gsv4bTtquj19mI90Zfkz9qYkipVEooFe7ikyZy9JiRg15Iqbcp2bP3nz9jfM5zNEyqZ9GRu3hg9yi2dsCzWxIsDpf6Xbj7ea7y7dTP+CteerW7Ke3kg4NxW83vx2h6poX17+4m6U7Sg9l623Z1cO2so3hhzTbaO3P3r0BwGcmXTpnCjvYutu5sZ/WW92jZ8n56/56uJP84O1jbpH5EHSs2bGfR0nVZqznG2dne3btSGzeSScdill58K+mwc3dnXp9brsK8t9mI90ZfSa2/Lwoi5aKEUgX6SwrFPkfD/kkaPjic29/o4KnN3Ss2Nq1to2H1b7n4vPNYtyvJ/DVdOPCf67qYOCLGvJYRdLS8Rjxm6anqkx4sRJUqFOf910peat3e45zxmPGPs6dz9JiR6bVaauKxjGnvk0nnsRUbuXbWUenYn1/9TkbSea8jEUwUYzDnlCmcMW0MTavf4aW33uWJV7sX+Vrw+zeZOHp/VmzYjgHnR2ppKalmrbnnTMtYoOuB5ZlXQRfq2pno7yTXAl9KJFJplFAkbz1GKm17M71vZG3QR+EEyeaxt7vo8HD+sMgyvDG6F6JqmFTP9HEH9kgopxx5SDpJ3P5MS8ZCUseMGclrb+9MX6fyu9e3svTNbelmpi9+Ygo3PPhK95uFx8WAkfvVpgvhf3nitYxzdiWdf3joFVKh3tfcyqIrG4Ggmat+RF16CePsFRIvmDGe+8O1Y2rjxvTDD+yxxv1AZCcPrc4o1UIJRfIWXbGx8eA4DY92TyeSnWzO+kANS7clMq5piRnUZbX5nz9jPPeGhTFAXdy4dtZRANz+TAv1I+qoicfo6AqaxlZt2gmEs6qkhjxHmplSi2Y9tmIj08aO4j+fW9Ojr2HJ8tYei3sZmQt+dYbvmVpkzCI1o+wVEhsm1bNozsn9Jp585EoeWp1RqoUSigA9vxX3Jj03VY7tGcmmPs4zre/z5I7u62aOH3cgc/96Wvpbd+p8i+ecnF5u+IIZ4wEyCtVPHnUoT736duag5eypuCKPjx4zMt0klWriiv5cuZqkZn/kcB59ZSMdYWKrrYmxeWd7+nn0fPF4z47wVBNUdo0qV+Hf12edK3moE16qhRKKpL8Vt3cmiceMebOn5708blR2sjmsNrPonj7uQBom1XP30nXpUV+pb+E3/83x6eOyC+XDRg5jWG0sMroqMynUxLoTUa5v+FdlDalNNVF1JIJrTL70F1O47uxjufTkyTywvDXdh7Ikq28EgnN/pqFn/0pKf4V/f81XuV6vTnipFkooQtPqd9KjrbqSztyHV3D0mJH9F1yp6cF7cf7Bndy3fXi6cDx/xnia17Yx9+EVdIXtRx05vsVnF6rnzxifHoVWP6KOx1Zs5H9atqaTy4UnTEy//oHlremfpbcaQrSJKlpA5+rovq+5lc6wuS0G1NXG0skrl/4K//6ar3p7vTrhpRoooQwhvTW1NE4ZTTxm6UI+6V6QdvqG/ZMsujKzcLz9mRYSkc6KmBmNU0b3iK23QrV5bRsrNmynJh4jkehOOKmf7/7m1u5JKWPWa/NQPgV06jqdVCKLjuzq73WDrcHken2+zZEi5aaEMkT01dSScQFl2AyVVzv9PfcE96mFfPLQOGV00HzVmSQWNq8BOWPL1feQHkIcMy46cWLG8N6m1e/QlUhNsAKfnTmhINeDFLIQz7f5KpVE9raTX6SUlFCGiP6aWgZ1AWVqFbxeEkrz+7GciSK7QM2nIzv7Z0gkncMP2q/XprJ4PIYTFMyFLoD3tsbQXw0kmjhjZiTCSSo1wksqnRLKEJFd2K5/d3ePwrbQ38ab3qvJmSiyz5PvKKb+jkslq9Rw38XPr2PJ8taCfqsv9DUh/Q0TBiceM9xdI7yk4imhVKHBfEMuRWGbEWNbgg0dRk3M0nOE9dWfkU8zUL7Hrdu2i66kD+q6jf4+20JfE5LPMOHsK/NFKpUSSpXZm2/IDZPqg36GQRa2eceYXoellpoa48ITJ2RMONlbbHsbR3T4sxNcSDmQb/X5fLaFviZkIMOEm9e27dUV+CLFpoRSZfb2G3K+BWLGN/Xe3mxY7vV2uxegMhKJJOOy+joGq78CP/XZpIb4fnzqIRnzfGW/V/RiylSy7e2zjX4ehbwmJN9hwpp+RaqBEkqV2dtvyPk0G/UovCbHcieVW27JHWN6GhYP1mEpULt/f8k0+7PpK5l8Pr0ODNy/7C0WzTm518+2v4slC91JP5ifXaQSKKHkqVDXAhSi8Nnbb8j9FWA9Cq/3avJeGAoi07Cs207jeacVrODLt1O+v8+mafU7dHZ1r9/SmQiuu7nqtKk5X99fzaUUNQdNvyLVoOITipn9M/DXQAfwBvC37v5uuO964O8IlhO/xt1/VYwYClVoFOp9in3VdI/C64Cu3AfeeWdwf9llPWOsjwfLDZf4Go58PpvGKaOprYmlayi18e4LIHO9vq/CvFQ1B02/ItWg4hMK8CRwvbt3mdktwPXAt83sOOAiYBpwOPCUmR3l7ok+3mtQClVoVFOzxfkzxnevC/JcL3l6+fLgPkdCKZZ8k2lfNcHUFfDZfSh9nbO3wryUNQdNvyKVruITirs/EXnaBHwmfDwbWOzu7cCbZtYCnAg8V+gYClVoVEOzRXYt6vw+5q2qVPnUBAdaOPd2fDFqDppqRapVxSeULF8Ewvk+GEeQYFJaw20ZzGwOMAdg4sSBz6ALhSs0qqHZImctqtxBDVCpa4KFrDloNJdUs4pIKGb2FDAmx64b3f3h8JgbgS5gYeplOY7vsdSFu88H5gPMnDlz0KuzFqrQqPRmi5y1qLf7f10lqYaaYG+qqVlUJFtFJBR3n9XXfjO7HDgH+JR7eqmjVmBC5LDxwIbiRDh05KxF9daIOGpUSWPLVzXUBHtTzclQxDx76bsKY2ZnAj8EPunuWyLbpwF3E/SbHA48DRzZV6f8zJkzfdmyZUWOeB+0eDGMyVWB7MemTXDRReoT6EOuz0afVyZ9HuVnZs3uPrO/4wZcQzGz/YE9xRhN1YvbgGHAk2YG0OTuX3b3lWZ2L/AqQVPYVSWMaUhqbktkLPGb12uGcJ9AfwVhb59NpTeLltJQ/vupRv0mFDOLEQzPvQQ4AWgHhpnZFuBRYL67v16sAN19ah/7bgZuLta5pVv3/FxQF+tk4QnDabj/jmDnlVf2+rqh2ieQT0E4VD+bgdBnVF1ieRzzDPBBgus/xrj7BHc/DDiFYJTV983sfxUxRqkA3fNzQWcyeM7KlcGtD6k+gbiRMW1+qaQmVCzlOSF3QZgt+tmovyQ3fUbVJZ8mr1nu3mlmVwPrgHcB3H0b8ADwgJnVFjFGqQDd83NBbSx4no9ST5sfVc7mknyX+q3WwQOlos+ouvSbUNy9M3w4BlhmZsuBO4BfpUZcRY6RfVR6fq4B9qFA6abNB7h76ToeW7GRs6aPpW1XR9maSwayxosKyb7pM6oeeXfKu/v/MbN/AP4S+FvgtrBT/Ofu/kaxApTK0VA/sEQSVYrhsHcvXccND74CwO9e38qX/2JKWYfgqiCUoWZAo7zc3c1sE7CJYGRVPXC/mT3p7t8qRoBSwQ49NO9DS9F08diKjRnPV27coeYSkRLKO6GY2TXA5cBWYAHw92HfSgx4HVBCGWpuvHFAhxf7G/tZ08fyu9e3ZjxXLUGkdAZSQzkEON/d10Y3unvSzM4pbFgiA3fxScFcbak+lNRzXRgnUhr5XIdyMsHFhHN7O8bdVxU0KqkOt90W3F99dXnjiLj4pInpRAK6ME6klPK5DuVyoNnMFpvZF8xsEHNwyD6ppSW4VbB8rgcRkcLIZ9jwlwHM7BjgLOA/zexAggseHwf+R1OeSCEVsolKky2KlM5Ahg3/CfgT8K9mth9wGvBZgokb+500TCQfhW6iGuoXxqn/SEppIKO8Pgs87u47gW8AM4D/6+5fLVZwUlqVUPgUY+6moTrSS/1HUmoDGeX1D+5+n5l9Avg08APgP4CTihKZlNSgCp/xhV8eWE1UhaOJFaXUBpJQUv0kfwX8h7s/bGbfKXxIUg6DKny++c2Cx1GIJqpKqGlVAiVnKbWBJJT1ZvZTYBZwi5kNI79RYlIFKqnw2ZsmKjXzdBvq/UdSegNJKJ8DzgR+4O7vmtlY4O+LE5aU2qAKnx/8ILgvQk1lsNTMk2mo9h9JeQxklNcuM3sGONLM/iLcvKc4YUk5DLjwaW0tXjCDVEk1LZGhZiCjvK4AvgaMB14EGoHngNOLE5rIwKmZR6R8BtLk9TWCJYCb3P208ELH7xYnLJHBUzOPSHkMpFN9j7vvATCzYeGFjkcXJywREak2A6mhtJrZQcBDwJNm1gZsKE5YUhWmTi13BCJSQQbSKf834cPvhJ3zBxLM5SUVpmTXYVTQLMMiUn4DWrExxd1/U+hApDB0HYaIlEvefShmNtPMHjSz5Wb2cupWzOCyzv9NM3MzOyR8bmb2YzNrCWOZUapYKllJp2u/+ebgJiLCwGooCwkuZHwFSBYnnNzMbAJwBrAusvks4MjwdhKaVwwo8XUYW7YU771FpOoMJKFscfdHihZJ3/6VYM36hyPbZgN3ursDTWZ2kJmNdfeNZYmwQug6DBEpl4EklJvMbAHwNNCe2ujuSwoeVYSZnQusd/eXzCy6axzwVuR5a7htSCcU0HUYIlIeA0kofwscA9TS3eTlwF4nFDN7Csi1tPCNwA3AX+Z6WY5tnuO95wBzACZOnNjjBSIiUhgDSShxT43LAAASpUlEQVQfdvfjixGEu8/Ktd3MjgeOAFK1k/HAcjM7kaBGMiFy+HhyXBfj7vOB+QAzZ87skXBkL0ybVu4IRKSCDCShNJnZce7+atGiyeLurwCHpZ6b2RpgprtvNbNHgKvNbDFBZ/z2od5/UnJXXlnuCESkggwkoXwCuNzM3iToQzHA3f1DRYmsf48CZwMtwC6CJjkRESmTgSSUM4sWRZ7cfXLksQNXlS8aYe7c4H7evPLGISIVod+EYmbmgbX9HVPY0KTi7dhR7ghEpILkc6X8M2b2VTPLGCJlZnVmdrqZ/QK4vDjhiYhItcinyetM4IvAIjM7AngXGA7EgSeAf3X3F4sXooiIVIN+E0q4Bsq/A/9uZrXAIcBud3+32MGJiEj1GNBsw+7eia5El5QZmo+zkpRs2QKRXgxq+noRAC67rNwRSEjLFkglGMgSwCJSoUq6bIFILwaVUMzs6sjjgwoXjlSVb387uEnZpZYtiBvFX7ZApBeDbfKaFHl8PaBSZShqb+//GCkJLVsglWCwCSVmZqcA/wPoq5BIBdCyBVJug+1D+RbwYeBnZC56JSIiQ9RgE8qPCWb9/Tsii22JiMjQNdgmrw7g7fDx6QRXzMtQc/LJ5Y5ARCrIYBPKLuDA8Mp5LYM4VF14YbkjEJEKMtgmr5uAN4DbgbsLF46IiFSrwdZQEu7+7wWNRKrPtdcG97feWt44RKQiDDihmNkXgEvN7H1gHfBtd3+/0IGJiEh1GUwN5VR3/xSAmX2IoPnrWwWNSkREqs5g+lDSy/S5+8togkkREWFwyaDRzH4MNIe3usKGJCIi1SifNeVPB15x9y0A7n6imY0HGoDPAZOLGqFUrlNPLXcEIlJB8qmhPAVsNrMksAJ4GXglvH/c3XWl/FB13nnljkBEKkg+CeUagjXl7wX+ABxNUDv5AnAsMKZYwUmF27MnuB8+vLxxiEhF6LdT3t1vAz4OOHAr0Al8zd1Pc3clk6HsuuuCm4gIeY7ycvfd7n4LcCowFXjezE4qZmBRZvZVM3vNzFaa2T9Ftl9vZi3hvk+XKh4REekpn075Uwiato4J7w8DdlKidVDM7DRgNvAhd283s8PC7ccBFwHTgMOBp8zsKHdPlCIuERHJlE8fym+Al4BFwI/dfU1RI+rpK8D3U53/7r453D4bWBxuf9PMWoATgedKHJ+IiJBfk9dXCFZm/CuCpq5XzeweM/s/ZlaKYT5HAaeY2VIz+42ZnRBuHwe8FTmuNdwmIiJl0G8Nxd1/Gn0eXoPyIeB44ALgob0NwsyeIvdosRvDGOuBRuAE4F4zmwJYrnBzvPccYA7AxImaab+gzjyz3BGISAXJpw/lMuCHBLWZ/wKudvdHgUcLFYS7z+rj/F8Blri7E9SQksAhBDWSCZFDxwMbcrz3fGA+wMyZM3skHNkLSigiEpFPk9dc4AyCTvl1wPeKGlFPDxGsComZHUUw1ctW4BHgIjMbZmZHAEcCz5c4tqFt+/bgJiJCfp3yO9z9j+HjfzCzpcUMKIc7gDvMbAXB0sOXh7WVlWZ2L/Aq0AVcpRFeJXbTTcG91kMREfJLKGPDfohVwJ+A2uKGlMndO4D/1cu+m4GbSxmPiIjklk9CuYmgE/4Sgo74A8zsUYKhxC+7+6IixiciIlUin1Fe86PPs0Z5nU1wfYqIiAxxA14Pxd1bCUZYFWyUl4iIVD+ttiiDd+655Y5ARCqIEooM3umnlzsCEakgg1lTXiSweXNwExFBNRTZG98Lr3HVdSgigmooIiJSIEooIiJSEEooIiJSEEooIiJSEOqUl8H73OfKHYGIVBAlFBm8j32s3BGISAVRk5cM3rp1wU1EBNVQZG/88IfBva5DERFUQxERkQJRQhERkYJQQhERkYJQQhERkYJQp7wM3qWXljsCEakgSigyeA0N5Y5ARCqImrxk8FpagpuICEoosjduuy24iYighCIiIgVS8QnFzD5iZk1m9qKZLTOzE8PtZmY/NrMWM3vZzGaUO1YRkaGs4hMK8E/Ad939I8Dc8DnAWcCR4W0O8B/lCU9ERKA6EooDo8LHBwIbwsezgTs90AQcZGZjyxGgiIhUx7Dha4FfmdkPCBJgas70ccBbkeNaw20boy82szkENRgmTpxY9GCHlCuuKHcEIlJBKiKhmNlTwJgcu24EPgX8b3d/wMw+B/wcmAVYjuO9xwb3+cB8gJkzZ/bYL3th+vRyRyAiFaQiEoq7z+ptn5ndCXwtfHofsCB83ApMiBw6nu7mMCmFFSuCeyUWEaE6+lA2AJ8MH58OvB4+fgS4LBzt1Qhsd/eNud5AimTBguAmIkKF1FD6cSXwIzOrAfYQ9ocAjwJnAy3ALuBvyxOeiIhAFSQUd/890GPSKHd34KrSRyQiIrlUQ5OXiIhUASUUEREpiIpv8pIKdvXV5Y5ARCqIEooM3tSp5Y5ARCqImrxk8Jqbg5uICKqhyN64667gXis3igiqoYiISIEooYiISEEooYiISEEooYiISEGoU14G7+tfL3cEIlJBlFBk8LRgmYhEqMlLBu8PfwhuIiKohiJ74957g/uPfazv40RkSFANRURECkIJRURECkIJRURECkIJRURECkKd8jJ4N9xQ7ghEpIIoocjgHXZYuSMQkQqiJi8ZvF//OriJiKAaiuyNRx4J7k8/vbxxiEhFUA1FREQKoiISipl91sxWmlnSzGZm7bvezFrM7DUz+3Rk+5nhthYzu670UYuISFRFJBRgBXA+8NvoRjM7DrgImAacCfy7mcXNLA7cDpwFHAd8PjxWRETKpCL6UNx9FYCZZe+aDSx293bgTTNrAU4M97W4++rwdYvDY18tTcQiIpKtIhJKH8YBTZHnreE2gLeytp9UqqAk9N3vljsCEakgJUsoZvYUMCbHrhvd/eHeXpZjm5O7qc57Oe8cYA7ARK3fUVgHHljuCESkgpQsobj7rEG8rBWYEHk+HtgQPu5te/Z55wPzAWbOnJkz6cggPf54cH/mmeWNQ0QqQqV0yvfmEeAiMxtmZkcARwLPAy8AR5rZEWZWR9Bx/0gZ4xyaHn+8O6mIyJBXEQnFzP7GzFqBk4H/NrNfAbj7SuBegs72x4Gr3D3h7l3A1cCvgFXAveGxIiKSpXltG7c/00Lz2rainqciOuXd/UHgwV723QzcnGP7o8CjRQ5NRKSqNa9t45IFTXR0JamribHwikYaJtUX5VwVUUMREZHiaFr9Dh1dSZIOnV1Jmla/U7RzKaGIiOzDGqeMpq4mRtygtiZG45TRRTtXRTR5SZX6/vfLHYGI9KNhUj0Lr2ikafU7NE4ZXbTmLlBCkb0xfHi5IxCRPDRMqi9qIklRk5cM3kMPBTcREZRQZG88+2xwExFBCUVERApECUVERApCCUVERApCCUVERArC3IfOBLxmtgVY28chhwBbSxTOQFRqXFC5sSmugVFcA1epsRUjrknufmh/Bw2phNIfM1vm7jP7P7K0KjUuqNzYFNfAKK6Bq9TYyhmXmrxERKQglFBERKQglFAyzS93AL2o1LigcmNTXAOjuAauUmMrW1zqQxERkYJQDUVERApCCQUws380s5fN7EUze8LMDg+3m5n92Mxawv0zShzXP5vZn8JzP2hmB0X2XR/G9ZqZfbrEcX3WzFaaWdLMZmbtK1tc4fnPDM/dYmbXlfr8WbHcYWabzWxFZNvBZvakmb0e3hd/CtiecU0ws2fMbFX4e/xaJcRmZsPN7HkzeymM67vh9iPMbGkY1z1mVlfKuCLxxc3sj2b2y0qJy8zWmNkrYdm1LNxWvt+juw/5GzAq8vga4Cfh47OBxwADGoGlJY7rL4Ga8PEtwC3h4+OAl4BhwBHAG0C8hHEdCxwNPAvMjGwvd1zx8JxTgLowluPK+Hf1F8AMYEVk2z8B14WPr0v9Tksc11hgRvh4JPDn8HdX1tjC/2cHhI9rgaXh/7t7gYvC7T8BvlKm3+fXgbuBX4bPyx4XsAY4JGtb2X6PqqEA7r4j8nR/INWxNBu40wNNwEFmNraEcT3h7l3h0yZgfCSuxe7e7u5vAi3AiSWMa5W7v5ZjV1njCs/V4u6r3b0DWBzGVBbu/ltgW9bm2cAvwse/AM4raVCAu2909+Xh453AKmBcuWML/5+9Fz6tDW8OnA7cX664AMxsPPBXwILwuVVCXL0o2+9RCSVkZjeb2VvAJcDccPM44K3IYa3htnL4IkFtCSorrqhyx1Xu8+fjA+6+EYKCHTisnMGY2WTgowS1gbLHFjYrvQhsBp4kqHG+G/liVa7f6a3At4Bk+Hx0hcTlwBNm1mxmc8JtZfs9DpkVG83sKWBMjl03uvvD7n4jcKOZXQ9cDdxEUAXPVtBhcf3FFR5zI9AFLEy9rBLiyvWyHNtKOYyw3OevKmZ2APAAcK277wi+dJeXuyeAj4T9hQ8SNK/2OKyUMZnZOcBmd282s1NTm3McWo6/tY+7+wYzOwx40sz+VIYY0oZMQnH3WXkeejfw3wQJpRWYENk3HthQyrjM7HLgHOBTHjaKVkJcvSh6XBV+/ny8bWZj3X1j2Hy6uRxBmFktQTJZ6O5LKik2AHd/18yeJehDOcjMasLaQDl+px8HzjWzs4HhwCiCGku548LdN4T3m83sQYJm37L9HtXkBZjZkZGn5wKpLP8IcFk42qsR2J6qSpYorjOBbwPnuvuuyK5HgIvMbJiZHQEcCTxfqrj6UO64XgCODEff1AEXhTFVkkeAy8PHlwO91faKJmz//zmwyt1/WCmxmdmhqZGMZrYfMIugf+cZ4DPlisvdr3f38e4+meBv6tfufkm54zKz/c1sZOoxwSCeFZTz91jqUQmVeCP4prYCeBn4L2BcuN2A2wnacV8hMqKpRHG1EPQJvBjefhLZd2MY12vAWSWO628IagPtwNvAryohrvD8ZxOMWnqDoHmunH9Xi4CNQGf4ef0dQdv708Dr4f3BZYjrEwTNMy9H/rbOLndswIeAP4ZxrQDmhtunEHwxaQHuA4aV8Xd6Kt2jvMoaV3j+l8LbytTfezl/j7pSXkRECkJNXiIiUhBKKCIiUhBKKCIiUhBKKCIiUhBKKCIiUhBKKCIiUhBKKCIiUhBKKCIhM/uSmW0K15ZYbWZfiGx3M/tk5Nirw22DmaIGM9vPzH5jZvEChd/beX5qZh/vY3+dmf3WzIbMNExSPEooIt0+BHzH3T9CMKXGv0S2v0w4UaGZjSC46n0LwQwKg/FFYIkHkyEW00kESx/k5ME0/08DFxY5DhkClFBEuh1PMHcUBFOlxCPbFwHHhM+vIZhqI+nubw/yXJcQzrFkZpMtWJlzgZmtMLOFZjbLzP4nXHUvvaaMmT0UTlW+MjVdeTin03+HKx2uMLMLw+3HAn9294QFq5J+LfI+N5vZNeHTh8J4RPaKEopIt+OBP4WTJ14D/DLcfizB6nzHmNmBBN/m/0Aw39SAhRNXTnH3NZHNU4EfEdSGjgEuJphz65vADZHjvujuDcBM4BozGw2cCWxw9w+7+3Tg8fDYsyKPf044YaCZxQgmOUwth7ACOGEwP4tIlBKKCME668ABwK8IJvyrB64Kt7/j7qsJFir6FvBvwFEEzWCY2Qtm9h9m9lS4xkjqPe8xs2/kON0hwLtZ295091fcPUkw0d/THky09wowOXLcNWb2EkEz1gSCGZ1fAWaZ2S1mdoq7bw+P/TRhQgmT1ztm9lGCWWn/6O7vhPsSQEdq5lqRwVJHnEjgQwSF+JnRjWb2Cbr7SXYS1AZOJFgPY3mYcJ5z92vM7C7gUOA9M5tNUMPJ1Wm/m2Bdjaj2yONk5HmS8P9puLjTLOBkd98Vrhcy3N3/bGYNBDMG/79m9gTwA+AgD9fLCC0AvkCwcNodWecfBuzJEatI3lRDEQkcTzANeK7tqYTyz8DV4Tf64wlqKA3AUWb2JPCau79pZsOBz7r7XcCB2W/o7m1APDxuIA4E2sJkcgzB4lOY2eHALnf//wgSyQzgNIL1OqIeJEiIJxDUxAhfPxrY4u6dA4xHJINqKCKB44FHe9n+AIC7/zKy/TjgVeA84H8TrHtyW7jv74EDzOwnwDQz28/dd2e97xMEfSRPDSDGx4Evm9nLBOvNpEZvHQ/8s5klCdZe+QpBTeT+6IvdvcPMniFYCz06uuw0cv/sIgOi9VBE9oKZLSGojSTM7H7gGwQLQ/1duP8m4HF3X5r1uo8CX3f3S4sU13LgpGitI+yMXx7G+3rWz3C9u79WjFhk6FBCESkTM/si8IsSXIuCmR1H0KfzoLt/I7K9DrjI3e8sdgyy71NCERGRglCnvIiIFIQSioiIFIQSioiIFIQSioiIFIQSioiIFIQSioiIFIQSioiIFIQSioiIFMT/DzyBfGXukt4YAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute average and standard dev of RA pm (weight by the inverse variance)\n",
    "pmra_avg, pmra_std   = compute_avg(data1['pmra'], data1['pmra_error'])\n",
    "pmdec_avg, pmdec_std = compute_avg(data1['pmdec'], data1['pmdec_error'])\n",
    "\n",
    "'''\n",
    "plt.hist(data1['pmra'], bins=25)\n",
    "plt.axvline(pmra_avg, color='r', linestyle='--')\n",
    "plt.axvspan(pmra_avg - 0.5*pmra_std, pmra_avg + 0.5*pmra_std, color='r', alpha=0.25)\n",
    "plt.xlabel(r'$PM_{RA}$ (mas/y)')\n",
    "plt.ylabel('N')\n",
    "plt.show()\n",
    "'''\n",
    "\n",
    "# Plot proper motions\n",
    "plt.errorbar(data1['pmra'], data1['pmdec'], marker='.', linestyle='None')\n",
    "plt.axvline(pmra_avg, color='r', linestyle='--')\n",
    "plt.axvspan(pmra_avg - 0.5*pmra_std, pmra_avg + 0.5*pmra_std, color='r', alpha=0.25)\n",
    "plt.axhline(pmdec_avg, color='r', linestyle='--')\n",
    "plt.axhspan(pmdec_avg - 1.*pmdec_std, pmdec_avg + 1.*pmdec_std, color='r', alpha=0.25)\n",
    "plt.xlabel(r'$PM_{RA}$ (mas/y)')\n",
    "plt.ylabel(r'$PM_{Dec}$ (mas/y)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's decide how many stars to throw out. Start with a simple +/- 0.5 sigma cut on pmra and a +/- 1 sigma cut on pmdec (a slightly stricter cut than on parallax, in an attempt to avoid contamination):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "data2 = data1[data1['pmra'] > (pmra_avg - 0.5*pmra_std)]\n",
    "data2 = data2[data2['pmra'] < (pmra_avg + 0.5*pmra_std)]\n",
    "\n",
    "data2 = data2[data2['pmdec'] > (pmdec_avg - 1.*pmdec_std)]\n",
    "data2 = data2[data2['pmdec'] < (pmdec_avg + 1.*pmdec_std)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (e) Plot a color-magnitude diagram."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that we have to convert apparent magnitude to absolute magnitude, and we have to account for extinction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The history saving thread hit an unexpected error (OperationalError('unable to open database file')).History will not be written to the database.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11d36e780>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvX98VPWd7/98zyThh0ZI+VHAEJSiqMG2kqjxtrW6Vbd2sVrQauXWWots79f73XrX3VurLfVLv7u3fbR7223rXUttr9t9IGUFBNurW7Vaf+waJKFaAoggmhBAwBADyo8kM+/7x5lzcmZy5lcyk/mR9/PxQDMzZ868TwY+r/N5/xRVxTAMwzASCRXaAMMwDKM4MYEwDMMwAjGBMAzDMAIxgTAMwzACMYEwDMMwAjGBMAzDMAIpCoEQke+LyGsi8icReVREJhbaJsMwjNFOUQgE8BQwT1U/DLwOfKPA9hiGYYx6ikIgVPVJVe2PPWwGagtpj2EYhgEVhTYggNuA1cleFJGlwFKAU045peGcc84ZKbsMwzDKgtbW1ndUdUq642SkWm2IyNPAtICX7lXVDbFj7gUagYWagWGNjY3a0tKSW0MNwzDKHBFpVdXGdMeN2A5CVa9I9bqIfAlYAHwqE3EwDMMw8ktRuJhE5NPA14FPquqxQttjGIZhFEmQGvgpUA08JSKviMgDhTbIMAxjtFMUOwhVnVNoGwzDKF36+vro7OzkxIkThTalqBg7diy1tbVUVlYO6f1FIRCGYRjDobOzk+rqas444wxEpNDmFAWqSldXF52dnZx55plDOkexuJgMwzCGzIkTJ5g0aZKJgw8RYdKkScPaVZlAGIZRFpg4DGa4vxMTCMMwDCMQEwjDMIwcICLcdddd3uMf/OAH3HfffRm/v6Ojg6uuuopzzz2X8847j7feegtwYgn33nsvZ599Nueeey4//vGP4963adMmwuEwa9asycVlxGFBasMwjBwwZswY1q1bxze+8Q0mT56c9ftvueUW7r33Xq688kree+89QiHn/v2hhx5iz549vPbaa4RCIQ4ePOi9JxKJ8PWvf50///M/z9l1+LEdhGEYo5LW9m7uf3YXre3dOTlfRUUFS5cu5Yc//GHW7922bRv9/f1ceeWVAJx66qmMHz8egH/6p39i2bJlnmBMnTrVe99PfvITFi1aFPdcLjGBMAxj1NHa3s3iB5v5hyd3sPjB5pyJxB133MHKlSvp6emJe37lypV89KMfHfTn+uuvB+D1119n4sSJLFy4kAsuuIC//du/JRKJAPDGG2+wevVqGhsbufrqq9m5cycAe/fu5dFHH+WrX/1qTmwPwlxMhmGMOpp3d9HbHyWq0NcfpXl3Fw2zaoZ93tNOO41bbrmFH//4x4wbN857fvHixSxevDjp+/r7+3nhhRf44x//SF1dHTfeeCMPPfQQX/nKVzh58iRjx46lpaWFdevWcdttt/HCCy9w55138r3vfY9wODxsu5NhAmEYxqijafYkqipC9PVHqawI0TR7Us7OfeeddzJ//ny+/OUve8+tXLmS73//+4OOnTNnDmvWrKG2tpYLLriA2bNnA3DdddfR3NzMV77yFWpra1m0aBEAn/vc57zztrS0cNNNNwHwzjvv8Pjjj1NRUcF1112Xs2sxgTAMY9TRMKuGlUuaaN7dRdPsSTnZPbh84AMf4POf/zy/+MUvuO2224D0O4gLL7yQ7u5uDh06xJQpU3jmmWdobHS6cV933XU888wz3HbbbTz33HOcffbZALz55pve+2+99VYWLFiQU3EAEwjDMEYpDbNqcioMfu666y5++tOfZnx8OBzmBz/4AZ/61KdQVRoaGrj99tsBuPvuu1m8eDE//OEPOfXUU3nwwQfzYnMQIzYwKB/YwCDDMAC2b9/OueeeW2gzipKg302mA4Msi8kwDMMIxATCMAzDCMQEwjCMsqCU3eX5Yri/ExMIwzBKnrFjx9LV1WUi4cOdBzF27Nghn8OymAzDKHlqa2vp7Ozk0KFDhTalqHAnyg0VEwjDKDFa27u9/H0gL7n8pUZlZeWQp6YZyTGBMIwSwu0h1NsfpSIcAlX6o0pVRYiVS5pGtUgYucdiEIZRQiT2EOqLaFw/IcPIJbaDMIwSwt9DKBzbQUSimvN+QoYBRSYQIvJp4B+BMPCgqn63wCYZRlGR2EMILAZh5I+iEQgRCQP3A1cCncAmEXlMVbcV1jLDKC4SewiZMBj5ophiEBcBu1R1t6r2Ar8Gri2wTYZhGKOWYhKI04E9vsedsefiEJGlItIiIi2W82wYhpE/ikkgJOC5QWWRqrpCVRtVtXHKlCkjYJZhjAy5npFsGMOlaGIQODuGmb7HtcC+AtliGCOKv77BahqMYqGYdhCbgLNE5EwRqQJuAh7L94faXZtRDATNSDaMQlM0OwhV7ReR/wr8DifN9ZequjWfn2l3bUaxkM8ZyYYxVIpGIABU9XHg8ZH6vKC7NhMII1v8vZGG+vcnnzOSDWOoFJVAjDTlcNeWi8XJGDq53IXmc0ayYQyFUS0QpX7XZi6ywuPfhZ7oi7L8N1tZdk29fQ9GWVBMQeqC0DCrhjsun1OS/6AtsFl4asZXEfUlY7/a2cPnH/gPHt7YUTijDCNHjHqBKGVcF1lYKFkXWTExlIy27mO9g56LKCzb0GaZcUbJM6pdTKVOqbvIiomhuuuaZk8iLI4o+ImqWtKDUfLYDqLEKWUXWTExVHddw6wavnPd+YR9fQAEqLIdnVEG2A7CMBheRtvNF9cxd1o1zbu7qBlfRfexXtvRGWWBCYRhMHx3nT9F1U09dp83jFLFBMIwYuSiDsFSj41ywmIQJYD1iyodLPXYKCdsB5FnhlvpbHekpUU5VOcbhosJRB7JxeJu/aJKA/+NgKUeG+WCCUQeycXibnekxU/QjcAdl88ptFmGMWxMIPJILhZ3K4bLH7lqdGi7PKNcMYHII7la3Iuhy2e5dY3NZWzHdnlGuWICkWeKYXEfLuUYKM/VXX9rezfrNnfyibOmMLV6DAvn15b878YwXEwgSoRC3sGXowtluHf9re3drN3cySMte+iLNWKqCgsL59fmw1zDKAgmECPEcBb4Qt/Bl6MLZTjuP/f7ONkXxd+jry9iDfqM8sIEYgQY7gJf6Dv4kQqUj/QuaajuP/f7SGjgSmVYykI8DcPFBGIEGO4CXwx38PmOpRR6l5QN/u8jHBIumzuVydVjWGTxB6PMMIEYAYa7wI+GVNdC75KyoWFWDcsW1PNE236unjedmy+uK7RJhpEXTCBGgFws8OWQDZWKYtglZUprezfLf7uV3v4om946zNxp1WX93Rijl6IRCBGZCfwKmAZEgRWq+o+FtSp3lPsCny2J8YZS2iWV0m7HMIZD0QgE0A/cpaqbRaQaaBWRp1R1W6ENM3JLsnhDqYhoKe12DGM4FI1AqOp+YH/s56Mish04HTCBKDNK/Q68lHY7hjEcikYg/IjIGcAFwMaA15YCSwHq6iw4WIqUwx14qex2DGM4iGpiNndhEZFTgeeAv1PVdamObWxs1JaWlpExzMgp5dbbyTBKCRFpVdXGdMcV1Q5CRCqBtcDKdOJglDaldgdugmaMRopGIEREgF8A21X1fxbaHsNwKaUiPsPIJUUjEMDHgC8CW0Tkldhz96jq4wW0qagp57vaYrg2t1Nr294eL6h+si/Kus2dZff7NowgikYgVPVFQAptR6lQqne1/oUfCBSBYri21vZuvvBzxwY/CjzSssfaehujgqIRCCM7SjFV1L/wV4QEROiPDBaBYri25t1d9CWIg0skal1bjdFBqNAGGEPDTRUNCyWTKhq38EeUvgQRcEl1ba3t3dz/7C5a27sDH+eKptmTqKwY/M9DKJ3ft2EMF9tBlCilWKxVM76KkAiq6vgSBUQHL7jJri3R9bRsQb3XEynXrqiGWTWsur2Jnz33Br/ffgAFKsIhrm+ota6txqjBBKKEKaVUUbfBXX/UqbvR2H9CArdecoa3g3Cvx+9uch8nup6eaNufV1dUw6waVtzSWBQBc8MoBCYQxoiwdnMnJ/sG+/SjCite2A0QtwsIClQnVmDXTz+Nl97oAtWcuH2SCUEpCbFh5BITCGMQQ71jTnyf+7hmfBVrWjsHTWBziW0q4nYB/t2Cm1q6cH4tC+fXIkD9jAnc9xtnRxIOCcsW1Afamum1xAXQzZVkGIAJhJHAUFNM3bRQ9+7+vmsG4gMhESLR9C1dRPB2AU2zJ1ERDnmjPVe37OGRlj30R5WqihA7Dxz1UlAjUWXrvp5hXYtfkHr7o6za2MG6zZ0lkz5sGPnAspiMOIJSTDNh3eZObzHv7Y+yelOHd55oVJPuHvz0R2HH20e9u/5Pnj3FK4yJRJS+iDoLeF+UTW/FZy0FnT+ba3HdV+7nKdldv2GUI7aDMOIYaqfVg0dPxj2eetpYqg4c5WRfFBWCV/AAVm/qYEdsd1ARDlEZdnYf4XAIVIlEFRHxgt0u82ZMyOhakhXqAXzirClOxlLs1OGQUDO+ivuf3WUBamNUYgJhxDGU9NnW9m6ee/2Q97giBF/95Ic4pSrM+lf2ZSwOAB88bSxb9vYQVYhEotx0UR0zJo6LW9BrxlexbEObJxICdB/rTXstQFycAVX6o+oV7fXFdkDuOS+bOzVvabSGUQqYQBiDyDZrp3l3F/0RJx4gwI0XOnM6Hnt1X8bnEOAvL53NlfXTeH7nIe+uP7Glhf/nZRvaiKoTk0i20/Ffy/3P7nJ2NDhuMHdj0xdRQOPEYUxliCnVYwpe0W0YhcQEwsiKoKygRFfOwvm1NO/uIoO4tIcCR072Z7yDufniOuZOq85qp1MzvipuMxMOQTTqBMclJGgsI+qGxpksnF8LOOm5pTzYyDCGgwmEkTGpZkkHLepjK0NO7YPAOR+s5rW3j6b0NrkB4kx3MOmOS4w3PNG239s1hAT+7JwP8sxrB4mqEhbh+otmDkptLbVqdcPIJSYQRsakaqKXuFi7orFucyerN3Ww/e2jKc8dDol3154LEusaUCcLyhWHqooQk6vHEFX14h2nTxw3SASsSM4YzViaq5Ex2TYIbJhVg+Kkr6bjzMmnsCONiGRDoph54gB8bM5kVi5pYtH82pJreGgYI4ntIIyMCcoKclNAIXi2Q6YDPnYdfI97Ht0COPGF4eKPi/hTZCsrQtx5xdmejeZCMozkiGoWkcQio7GxUVtaWgptxqgkyIXjVjn700Fb27v5/M9eyqiSGuATZ03mX75ycc5sTDecKFdYQz+jlBCRVlVtTHec7SCMIZHowoH46mN/bOLsqaemjUG4XD1vetzjbBfexOOTpcjmkmKYgGcY+cAEwhgSqVw4iQN+MhWHcAjmTquOe282C2+hFupimIBnGPnABMIYEkHxiKA7/Wx6GUWj8L0ntjOmMszV86bTfaw3q4W3UAv1UNuTGEaxYwJhDJlMXDhNsyeRaSsmBV6ONeF7Yec7fPXS2VktvIVaqEtxup9hZIIFqY28c/uvWnhq24Gs3/eR2gksu6Z+WDEIwzAGk7cgtYicApxQ1ciQLEt//jDQAuxV1QX5+AxjZJlaPWZI72vbd4S1mzuzGtxjhW2GkTvSFsqJSEhEbhaR/yMiB4HXgP0islVEvi8iZ+XYpq8B23N8TqOALJxfS2U404qIASJR5eGNHXz+Zy9x+69aaG3vTv8mwzByRiaV1M8CHwK+AUxT1ZmqOhX4BNAMfFdE/nMujBGRWuAvgAdzcT6jOGiYVcP/99l5KY+ZOL4y6WuRqPLUtgN8YcVLJhKGMYJkIhBXqOp3VPVPquo1TVDVw6q6VlUXAatzZM+PgP8OZNCcwSglguY1+Lns7Clxj4XBfzl7I8razZ25NcwwjKSkFQhV7cvFMekQkQXAQVVtTXPcUhFpEZGWQ4cOpTrUKCKaZk9ibGXyv26b2rs5pSrsPVag8YwawqF419TqTXt4eGNHvsw0DMNHxllMIvLXAU/3AK2q+sqwDRH5H8AXgX5gLHAasE5Vk7qvLIuptGht7+b2f97E4WOZ3U8svriO+hkT+PHvX+ftIwMjTUMC//9158f1bLLsJcPInEyzmLLp5toIfBU4PfZnKXAZ8HMR+e9DMdKPqn5DVWtV9QzgJuCZVOJglB4Ns2r4mz8/J6NjK8JC/YwJLP/tVg4ciZ93HVVnmpwbj3ArqP/hyR0sfrA54zhFa3s39z+7y+IahpGEbARiEjBfVe9S1btwBGMKcClwax5sM8qQmy+u46uXzk573PyZE71K6qA9biSqXpV2UAV1OoYqKoYxmshGIOoAf6SxD5ilqseBk8FvGRqq+gergShf7v7Mufz9585PeUzH4WPsPOD0cBKgIiEWIeKMEIXs51TA0ETFMEYb2RTKPQw0i8iG2ONrgFWxwrltObfMKGu6j/WmbMHx9pGTrH9ln/f4P31oEv/xRhf9sbbhqrD8t1uZO62ahlk1LFtQzxNt+71usO6cioZZNRnN0bb+SYYxmIwFQlW/IyKPAx/Huan7qqq6EeLF+TDOKF+aZk+iMiz0RjJLknhh1ztO6qs4MQgFevsG7vyX/3Yrvf1RNu7uAhH6I05H12UL6r3XMpmjbRjGANm22tgNhHGyjMaLyKWq+nzuzTLKnYZZNaxaeglrN3ey68BRDr/veC93HXo/8HjVwbuNKI6bye8ucgTHOfJEX5Rfvrg74znahmHEk7FAiMgSnDYYtcArQBPwEvBn+THNKHcSF+iHN3Z4Y0czQXBcVU2zJ1ERDtEbMPx616H3CcWODYfEXEmGkQXZBKm/BlwItKvq5cAFgFWqGTkjXbX1IATPPXR9Q23S+ddRnD2FlecbRnZkIxAnVPUEgIiMUdXXgLn5McsYjaSrth6Ewo7YtLp5MyYQDgkhgYokp+hPaNVhdRCGkZpsYhCdIjIRWA88JSLdwL407zGMjHEDx2s3d2bUTkNxCuY6ut7nwRffpD+qVISEJR8/k80d3d7wIT/uLsPmSDtYBbqRiiENDBKRTwITgH9T1Sz9ArnDWm2UL63t3XzloZd593h/2mPdzCYYiDVEVRERItGBv9/hkPCda+fRfayXfe8eZ9XLHUQVwgJ/fdVc7rh8Tp6upjgxkRy95G1gEICqPjeU9xlGpjTMquEXt17EF37eTF9/lIqwEI0qQVmxPg1ABKKqznOqhEOCqhISZ2fhprxWhISKcIhIZPTWQRRqhrdROmSTxdQI3AvM8r9PVT+cB7uMUY7r+rjvmnovU+mprW/zwPO7U77vsx+ZweNb9nv1FarK+adP4MYL62jb18PJPqd1RySq3HjRTE6fOC5lQV05Y8WCRjqy2UGsBP4W2IIlhBh5JJnrI5N2GG37jvDRmRO9+ENU4dXOHrbvb2OgQsJxN/lHmY5Gd4sVCxrpyEYgDqnqY3mzxDBiJHN9uFlO7i4giF0H3wt8vs/nmxLghsaZwEBLjtHqbrFiQSMV2QjEt0XkQeD3+Jrzqeq6nFtljGqSuT78d7w146u499EtSYUiEZFY8DqqVFaEqJ8xIW7HsGxBfUp3y2hzP5Uy9l3ljmwE4svAOUAlAy4mBUwgjJySyvXhv+OdO62abz66hdcPHA0MXvuJKogqN11Ux8L5tYN2DN3HepN+ZqL7admCgbiILUDFxWh0FeaTbATiI6qaukezYeSITFwfDbNqeOLOS2lt7+aGB/4jLpspiEgUDh09ScOsGna8fZSQCKh6O4ZknxnX66kvyrfWbyESS4/9TsJkO6OwjFZXYb7IppK6WUTOy5slhjFEGmbVsPQT6YcQATy17QDffXw7y3+7lagqoZCwbEF9ykXEP29CBG+3ElH41votGVViW9X2yDCU2SBGcrLZQXwc+JKIvIkTgxBALc3VKAbu/sy51E06he8+sZ0jJ5IX1ynEpcoKmrYHlH/exMm+SFyFdlQZdJea6ANP5vYwX3nuscys3JKNQHw6b1YYRg64+eI6/rDjIE9uO5Dxe5SByXSJuAt4zfiquAK7cGigQrsyLNSMr+L+Z3dRM76Krft6eKRlD/1R9cQg2fQ685XnB8vMyh1pBUJERB3a0x2TW9MMI3PcxTxbogrf2tDG6k0dXDJ7EtXjKqkZX8Ufdhzk968d9KqwI1H1FdjVedPw5s2Y4IlHYgzEFYOgrCzzlRulQCY7iGdFZC2wQVW9DmoiUkXM7QQ8CzyUFwsNIw1+F47TQkOIuJXUGbw/ElVe7ezh1c6ewNfV17KjsiIUV2B376NbAusyBOKC30FuD6tiLn3K3U2YiUB8GrgNZ/70mcC7OBPlwsCTwA9V9ZX8mWgYqfHfjbt3+G4LjR1vHw2MHWRDOCQsjzX58y8Ere3dPNKyJ1CEPlw7gWXXDAS/E90e2fjKy30RGiqF/r2MhpTatAIRmwHxv4D/JSKVwGTguKq+m2/jDCMTEl04/jv8hlk13HxxnfePua/fudtPlxLrEhZYfu08L5XVzUZy3UT9sRNJ7I/ixCVccRjuIjYaFqGhUAy/l9HgJsyqm6uq9gH782QLsXkTDwLzcP6t3aaqL+Xr84zyIJO78cQqbLeWIRUfOKWSv7nqnDhx8C9Kt15yhlNLgQYW0KXKXlq3uXNQMDsoE2rvu8fLfhEaCsWwOI+GZodDavedR/4RZ8bE9bEYx/hCG2SUBpkW1vmP+eb6LSl3Et3v97FsQxvgZEjFFcz1R3nwxTeJRJ1aik+cNYW506rjzr9uc6cXn0jMXvLHLRIXuIc3drBsQxtRVSrCISpiWVPluggNhWJYnEdDSm0mWUwfAMaqal6nx4nIacClwK0AsUFEBRtGZJQ3N19cx9xp1fzsuTfYfeg9dh16f9AxCvRHlW+t38LWfT3Uz5jgLUqSkNn09LYDPLfjIDc0zmTh/FqAuPiEhCQue8l7nviCrtb2bpZtaPNcV/2RKF+4qI4ZvrbkRvEszuWeUpvJDuIHwE7gfwCIyH8AncBm4F9UdW+ObJkNHAL+t4h8BGgFvqaqg//lGkYOaJhVw4pbGrn/2V38w5M7ku4mIgorN3ZQFRbu+6wTrHZrI9ydgAK9EeXhjR2s3dzJwvm1cR1k3STwptmTqAiHvCFINzTOjIuZNO/uipuCFxJhYex1f/yjnBelTCn3xbkYyKTVRgPwXd/jauAXOMHqb+TQlgpgPvBPqnoB8D5wd+JBIrJURFpEpOXQoUM5/HhjtOK6K0ICFSHhqvM+6M2u9tMbUdr29XDH5XOYO62aRfNrufK8D1IVFu94BU70Rfm3trcR30lUdaBOQwcC235xcG0ZUxkihGPL8mvnxcUz/uHJHSx+sNladhgjQiY7iJMJRXDPqOrvRORJIJcB5E6gU1U3xh6vIUAgVHUFsAKcmdQ5/HxjlOJvpXH1vOnMnVbNM68d9Nw8foTBwer7PjuPtn09rN7UQSTW5/jw+453NBQTiSpfgVxfxHFN9UcGRMPvKglynSQGZddt7iy4e8UofzIRiBMiMsutpFbVr8X+r7G015ygqm+LyB4RmauqO4BPAdtydX7DSEZre7dXDb3prcOsXNLE8mvnsWxDmxdncGMFyVqFL5pfy+pNe0gszZs9+RSmTxzH1fOme11k3SOiwNHjfYGZTomLvj8oGw6HUmZAGUauyEQg/g5YLyJfUNXX3CdFZHqG78+G/xdYGctg2o0zg8Iw8kpQyqTrRnLTYhOL5KoqQvTGgtU146to3t1FNGDH8WbX+7xx6H02vnmYudOq6T7WS0icOoyQwNb9RzJK1/TvLPa+e5xfv9yR9j2FLiQzSp9MCuV+F8swelZEXgHaYi8tBL6ZS2NiFdmNuTynYaQj1QS7ZIv1sgX1XrbRsg1tLPjw9MB2G67Lqbc/yu2/auFDU06hIhwiEnE+6+p509n01mHvs93Gf4kV2+5Cf8flc7w6Cvc9R4/38cVfbOTqedOT1mwE7TJMQIx0ZLQDUNVHROT/AJ8B6oHjwEJVfTWfxhnGSDCUlMnuY71etlF/VHns1cFZ4ImCcfj9Xg6/30tY8CbbASycX4sA9b7Gf+6iDgNdXyvCIa5vqGXR/FovZjLplCqvffkLO98BBtdsnOxzYhappuSZm8oIImMXkaoewwkcr8mfOYZRGLJNmWyaPYlwSLxAtsZcRoleJrf9hp+IQsfhY+x4+2icIACBbcH9xXmrNnawprUTVAOD6E+07efmi+s8+6KxgPgjLXu8dFmAtQFFfCYQRiLZTJQzDCNGw6wall87j4qQEBIYUxli9uRTAo+98IzBC++/73qHb63fwom+gVGmbXt7nHkTvmloNeOr4kTHXdB7I0pUBwvS1fOmez/7YyL90YGMqdb2bta0dnrCFY4V8RlGIsXWasMwSga3Gtt1Te14+yj3PLol7hgFLps7lYa6Gta/spdwSNjfc2LQwh4F/tTZQ2VFiBsvGiiea97dFbcLEaAiLHFFeNd9dAZd7/fGxSCad3fFfUZIBkSgeXcX/bHgiAA3NM4smt2DxUWKCxMIwxgGftdUw6waOrrejxtpCk4q60MvveXFESpC4tVC+FEcN5LEzgWOK6syljEFjhvr8rlTeXr7AS8T6qwPVvOjy+fEncstuOvtixLyFdy5r/mD8m4spNBYXKT4MIEwjBxy92fO5cjJflZt7EAZnMoaiUS56aI6FFjT2um1H/ez6uUO6mdM4OaL62iYVcP1DbXe+VThwJETcZlQQe6hVIH3TIPy/rt5IPD4XN7xF0OHViMeEwjDyDGL5td6i39FeHAqqxssXjS/1mv73etzGbljUMHJlpo3Y4K3G4gCW2KxCjcTKlX66h0JOwuXdEH5xCl9iNAfGdy2PJM7/kxFpBg6tBrxmEAYRj5wu9OoMndadeAdu7tIL5xfy/LfbI0beRqJ1VdEdWDWxBNt+/n3Xe94k/NmTBzHjreP8qOnX/fiD8N10wTOoYg47QgTM54yuePPxp5i6dBqDGACYRg5xp0057YCdyuzU91d33hhHVv3bSEWanBSVFXj2nncecXZbHrrML19TgX3zgNHWf+KU3/h1kA8u+MgJ/qck5zsi7I2of4h6LODBhz551CEYzuIRJdWJnf82bqNSq1Da7kH1U2DvCC6AAAbN0lEQVQgDCPHZOoqSby7Xn7t+bTt64krmvOfw1/BHdXBxXk/+N1rHD7W5z124xyLYkHoVGJwfYNTrJcYK5kxcRw146vYuq8HhUHjXNPd8Zez22g0BNVNIAwjx2TqKglq+vf3nzvfe93fC8qtYeg+1uvtLBLxi4NLJBLlZ8+9wTOvHfTcVa5tiQV4lWGJC3672U3+RXBRQsZTujv+cnYbjYagugmEYeSBTFwl6e6u3ff7F+hlC+qpCA+kvaZCgFBIeGrbAS9Tqje2kLmf7R94FIkqN140k9N90+vuf3bXsBfBTALipSgg5bw7cjGBMIwCkcndddAuw5/2GsQHxlcyZ+qpAGx6qzvuOBFh77vHAVi5pMnLonJnXvvdUZD/RbCU3TTlvDtyMYEwjAKS7u462QK9ztdLyY8AR0/209LeHeiGUlVWbexg3eZOVi5p4u8+dz71MyZ4w5KAQQt2PhfBUnfTlFpQPVtMIAyjiPHfpbqxiKbZk1i5pIm1mztZ09rptc1AnZYd/jYcfsTXTLC3b6AZoNswcOPuLs6dflrgbIyhpMpmIijZBPTL+U69WDGBMIwiJygWsXJJE3//ufNZFJtw5w4RcrcUQV1k/YODQyFn0NGPnn7d24n0RtSrxRBAQsKTW9+mZnyV1+PJT9Ci7bqMTvZFCcdafAS9139tiTuUVOm3mbih3PcHDXoyssMEwjBKgGSuGPePf4hQOCRcNncq3cd62dzxLhpVKsID1dAhEZZ8/Exv5xC033BnZr/a2cOrnU4DQv9Cn2zRbt7d5QmOO0xp7rTqtMOKEgXGf95s3FB+gXJbnZRabKOYMIEwjBIgk4wnvyvKXfwrQsINFw8MJ1q32WnzffRkv7fohsSZnf3GO+/H7TL8PNG2Py7t9om2/YHzJBLnZERV4xZ0V8iSzdQOEoNs3FA/evr1ONFLJiq2y8gMEwjDKAEyyZhx78T9qaluSw53l/FIrEdU2FfzEA4JHd3HB/ukfJzsi/CFFS/RH3VqMFwXVggGjWldfu28uDYh7muJd/cwePEOEoNMrj3x3J59wiBRCbKjqiLEqtttl5GICYRhlAiZZsykynxy6yf6I8pV503lIzMnsu/d46x6uSNuYfUjDE6XdRffj82ZzJ1XnB1nl39Ohr/Iz90d+GdbJC7eycQg3bX7zx0CPnbWZK6eN92rAE92rEtv/+CxrIYJhGGUHckW2cSF8sCRE97ivNaNX4RDRKNRrydUSJz/u9lP/p1DVUVokDj4bYD4wPqtl5wRl3p75Xkf5C8/+aFB7x9K6miiKN55xdnAQIbWv27a4wXM3WPdnlUuKTZQoxYTCMMoQxIDv827u5g3YwJVYfFai2/Z28PiB5sH1TqAIxiJPaHCsZ5N82ZMyMhvnxhP2Lr/SJzAfGTmxJxVWAeJot/VFtX4gLmbJvyvLXuIRJTKsAxqI2KYQBhGWZOYFXTfZ+fFtQ1PVuvg/3nutGovuL0oYP5EMhLv6hPnYqSqyh5KhXXizqNp9iRC4nTFBWdGtz/7CwZ2RNlc12iiqARCRP4bsATnO9sCfFlVTxTWKsMoXYJadbhtw7Npn7E2Fr9wK7CDFtOHN3Z4FdnuNLzEncmi+bUZLchB2Uzu85lmHDXMqmHJx89kxQu7UYWqysEB82RNCA2HohEIETkd+CvgPFU9LiL/CtwEPFRQwwyjhBlqVpCfTBbrhzd2cM+jTr3ECzvfoaPrfarHVXpT7fwLcjgkbN3bw40X1iUtoku0u2Z8FV9Y8RK9ESUs8J3rzk9ZgAeOCDz00luAM19j2YL6lOm0toMYTNEIRIwKYJyI9AHjgX1pjjeMUUmm/vlMs4JSnS9osU50/zzRtj/uPSte2A0QNw3PiwekKMBLZvfazZ1e7CSi8M1HtwQW4Pnxi4CgdB/rTXpN5diJNRcUjUCo6l4R+QHQARwHnlTVJxOPE5GlwFKAurrUdxCGUY5k65/PpN12qvMlLtZBd99Xz5vuTbWD+J5Pyza0EYlN2EtkxfNvJF3o/Xav29wZ91oUUk7Lg3gRCIdD7H33OK3t3d55y70Tay4IFdoAFxGpAa4FzgRmAKeIyH9OPE5VV6hqo6o2TpkyZaTNNIyCk8zlk8/zNcyq8QLZ7sIbFryFd+60av7+c+fzkdoJzojSGE6DwGBxAGjvOsbiB5tpbe9Oal9rezeHjp4c9LwEHJto88olTdx0UR2o8uuXO+I+y39NRjBFIxDAFcCbqnpIVfuAdcB/KrBNhlF0+BfoXLhHgs7X2t7N/c/uSrpwf+KsKZw55VSivoV37rRqrqqfhsayhgQ4c8qpVISdc1dVhKgKi9MIMHYexSlS+9HTr3uf1drezT2PbuHeR7fw8Ebn3E9tOxD3+eGQeO1DUtEwq4YZE8d5FeDJBDDd9Y5WisbFhONaahKR8Tgupk8BLYU1yTCKj1y7R4KyjZK5nFrbu71gsZ/EvknujuSNg+8RDsG80ydw44XxFdZuEVtU4cWd77DprcMsW1DPfY+1eecPh4Sozz0lseeWXzsvo+tube9m77vH40apJgpqKQ8tyjdFIxCqulFE1gCbgX7gj8CKwlplGMVJptXG2QSz3ddTjRlt3t01aN6Ev2WGKzY/evp1Xtz5TqyrK/yps4cdB7ayckkTd1w+B3DqK/zH9fVHeaJtf9z5o1ElHBJU1SvUc1Nk012bf+GvCAmfOveDTK4eM+g4y2hKTtEIBICqfhv4dqHtMIxyYKh3xqkyfJpmT6LSV41dERZubJzJQl9dQ8OsGu684mw2vnnY6/2U2PXVPe7qedN56Y0uorGRp1fPm87G3V3e+SsrQtx3Tf2gyu1Mrs2/8PdHlGdeO0hUdVAth2U0JaeoBMIwjNwx1DvjVC6shlk1rFp6ideKY2GKgrdoQu/wcHhwV9Xlv91KVJVQrE7BbfT3wHNvcPDIiaS1Ev5rO9kX5WfPvcFHZk6MGzK0793jVISESFSRWEV10O/CMpqSYwJhGGXKcO6MU7mwMnFvNe/uIpLgiqqrGceOt48GpsuCsnpTB1v39VA/YwIv7DxEb3+UHQe2AgzaQTTNnkRFOOR1ZX1y2wGe3n7Aq7vw5mGEQ9x40Uzm+XpK+QPxQUOLXGzMqQmEYZQthbwzbpo9icpYsNpl16H3uefRLQgwptJZyP0Bbad4rseLOUTVyXD61oY2zwXldzdd31DLqo0dg4YD+YvyIpEop08cF9eCPFkgHsgoUD+aMIEwjDJmKK2zc/W5q25vYvlvtnpzrl3ceET3sV4voO0vsotElYqQICiI4yKCAbHQ2CCiZQvqGVMZorcvSpSB4UDJmgKmCsSv29zp9Zuqqghx6VlTAifmjTZMIAyjTCg2l0jDrBqWXVPPF37eHLeT8E+hcwPaL73xDr5DWPLxM6keV8kre96Nq4FwxcIvMEGjQ/27haDfRaL7za3HiKpT/f371w56O5NwSEZt4NoEwjDKgGLN5Xd3Em67cP8sCXDu5JtmT+LGC+tYubEDcHYC1eMquePyOTy8sSNOICrD4rmbksUO3M9Nd/0L59d6gXZw2nn09Ue9gDY4Kbw3NM4sit9lITCBMIwyoJhz+ZMFgP2C5sYj+mKBZVdAuo/1EhK8OdiXzZ3KRxOyldwdRNu+nrSZVUGf7R7v3434A9qZVGyXKyYQhlEGlFou/7rNnXE+/q37esBNi1X1sp1qxldREXLqLhR4dsdBpsaK3fwLvX+U6SOtnay6PfkOKpmY+oUsnYtqtGACYRhlQCnl8re2d/NIy54BH384FKu4dkSgP+KMB43GgtGXzZ3KU9sOeK89vLGDtZs7WTi/dpA4QPodVCZiWqjgfrExKgWi2IJ5hpELSmVRa97dRX90wMfvts9IjAG4d/hTqscwpjLk7TjcXYfgNAB0n3cJh1MHlQslpqW47ow6gSjWYJ5hjBYS7+AXZRADqJ8xgW+t34JbexcOh6ifMQGAtr09/KmzB8URnBszCConE9N8LeKluu6MOoEo5mCeYYwGMplylxgDaN7dFdfR9bKzpwxUS4eEyoqBbq1DDSrnaxFvbe/mR0+/XpJ1FaNOIEotmGcY5Ug6d1ji64n/bidXjxmolo4qN140k9MnjhvWnX8+bh79oqPE14CUAqNOIEopmGcY5chQ3DhBMyvcmIXfTTUcmmZPoiIk9EU0Z8VxftEJCXxszmTuvOLskll3Rp1AQOkE84zRRSkGMbNlOG6cxH+3ub7R2/H2USLqBMGRdANN0xM0rKiUxAFGqUAYRrFRqkHMbBmKGyeZcObyRq+1vZtlG9q8Vh79AbZlI+CJw4puuqgubQFfMWICYRhFwGhJnsg2BjgU4QxayNMt7s27u+LmV4QSXEzZ2uH/PiNRZcbEcSX5fZpAGEYRMFqSJ7KNAWYrnEELOaRv3e2fpR2SwTOvs7WjXL5PEwjDKAJGU/JENq6hZAttsh2BfyHvjbXx7jh8LG2Kabrff7YLfrl8n6IJYwFLicbGRm1paSm0GYZh5JFEMYjz74dDXiV2w6wavvv4dh54frf33opY99eoOimmVZVDj++UUxKBiLSqamO642wHYRhGUeJfkO+4fI73fOIuYdXGDtZt7mTZgnoefPHNuHNEYk3+/CmmMNBmPNOFvpzEIRtMIAzDKDpSBYVdd09ib6Yn2vbHB5oFRMSbQOeKg/+8t15yBlv3H+HqedO5+eK6lLac7IsSDjnxiWTHlhuhkf5AEfmliBwUkTbfcx8QkadEZGfs/6NHog3DGERQUNjF9e/ffHEdVWEh7Bs1WlURIiTOFLhQyHEvhURYtqDea9nhnvdkX5QHnt/NCzvf4Z5Ht/BwbGBRkC2uGPVHnU6zre3dI/SbKCyF2EE8BPwU+JXvubuB36vqd0Xk7tjjrxfANsMwCoR/+M++d49TEXLmUQcFhd1A98L5tXHV1e6UOAV+/XIHCkSjyhNt+5k7rTou2KwMjKAAWL2pI3BnUDO+CpGBY6OqZZuGnMiIC4SqPi8iZyQ8fS1wWeznfwb+gAmEYZQlyeoU/MN/BGdXcONFM5k3Y4K3gwjKPkoMXPsn1J3sixIFXtz5DpveOszKJU1edtHOA0dZ/8o+71zb9h+htb17UHHc8t9ujROSkAg146vy9vspJoolBvFBVd0PoKr7RWRqsgNFZCmwFKCubnT4AQ2jXEgWW/C7fsDZAUQizswHt2uru/C7M61T1Sl0H+tl2YJ6vrl+Cxprn9Ebc1W5Ae+fPLMzzrb+iLJuc2fged1W4iLODmL5b7cyd1r1iLYMLwTFIhAZo6orgBXgpLkW2BzDMLIgcSFfu7nTcyu5hWpuSmplRchb2N2MJf+kuaDAtb9OoXl316A7f9cV5drhR4FHWvbEtcTwn1fEcXmlqqcot5YpxSIQB0Rkemz3MB04WGiDDMPIPf4FNxwS1rR20h9x+hVdNncqU6rHUD9jgrdLAOImzSVboJMVpo2pDNHbFyUUiq+O9tsBEI3tMiJRTXrexEFGQcVy5dYypVgE4jHgS8B3Y//fUFhzDMPIB/4Fd9+7x1n1coezO4goT207wJjK0KCmdpku0IkV2qmqmbNZ+FMNMkqkXFpsuIx4JbWIrMIJSE8GDgDfBtYD/wrUAR3ADap6ON25rJLaMEoXf32BuwqFBf76qrlxhXGJ78mmo2o+jk1HKcQgMq2ktlYbhmEUjIc3drB6Uwfb9h8hEqtZyEUh2kjGAkpBEBKxVhuGYRQ1bgqp00EVL8aQKkMoU5IV2mXbBjyTayinoHQiJhCGYRQE/yLuODJSZwhlQ2IsoGZ81aA6ibZ9PV6QfKiLe7qgdCnuLvyYQBiGURASM5oQ8UZzDje4mxigTmzw506Pcx3s6UQp2UKfKihdDrsLEwjDMApC4iIOg11Awz2//zwVIaEvogjEiYNbtZ1MlFI160uVKZXY92ltQhFeKWACYRhGwQhKTc0b4nRpCoWEMI5IhEPCDY0zU86LDmrW54+RJBuA1DR7EhUhoTfWcnxNa6c3t6JUMIEwDKOsaW3v5kdPv05/JOo177vpojpmTBxHzfgquo/1Jn2fWycRDgn9sT4g2TTrO3f6abza2QM4rUNKrXDOBMIwjLIlsdYiJBAOOy083AK5oD5PED83YsnHz+TBF9/02nyki5H4P5fY55Zi4ZwJhGEYZYu/2V4IOP/0CWzff4Rfv9xBSISoamCfp4Xza+Oyk6rHVbL6Ly/JOEaS+LnuNLtS2j2ACYRhGGVEYrZRYpZR/ekT2LK3x+kaq048QtBBfZ7eOXqSUGwIhHvnnyzWEETi55aiOIAJhGEYZUJrezdf+Hmztyivur0pMFPKbf5X6XMr+fsxhcMh/rDjINGYgLjT6LIhVXZTKWECYRhGWbBuc6fXwru3P+rNdki880+2cLuN+Pa+e5xfx5oICpo0iJ2ObHYc2TJSBXgmEIZheJRy5W9iV7lkXeaSLdz+6XT+XUbT7ElF9XsZyQI8EwjDMIDiq/zNdlFeNL+WNS176IsolWFh0fzaIX1ukFuqmH4vIzlzwgTCMAyguIbdDEWsGmbVsGpp5plG6c7lvv/+Z3cVze8FRnbmhAmEYRhAcQ27GapY5cPvX0y/FxjZALgJhGEYQHFl3hTTolxMvxe/TSNhhw0MMgyjKCmmwHC5YQODDGOUUi4L60jdJRvJMYEwjDKi2DKREikX8RotmEAYRhlRDJlIyUSg2MVrqJSz6JlAGEYZUejgbioRKAbxyjXlKnouJhCGUUYUOuMmlQgUWrzyQTmKnp8RFwgR+SWwADioqvNiz30fuAboBd4Avqyq7460bYZRDhQyuJtKBAotXvmgHEXPz4inuYrIpcB7wK98AnEV8Iyq9ovI9wBU9evpzmVproZRfJSzTz6IUrzeok1zVdXnReSMhOee9D1sBq4fSZsMw8gdoy09tZyvtxhjELcBq5O9KCJLgaWxh++JyI4RscphMvDOCH7eSFKu12bXVVrYdY0MszI5qCCV1LEdxG9dF5Pv+XuBRmChFmGJt4i0ZLItK0XK9drsukoLu67iomh2ECLyJZzg9aeKURwMwzBGG0UhECLyaeDrwCdV9Vih7TEMwzAgNNIfKCKrgJeAuSLSKSJfAX4KVANPicgrIvLASNuVISsKbUAeKddrs+sqLey6ioiS7uZqGIZh5I8R30EYhmEYpYEJhGEYhhGICUQAIvJpEdkhIrtE5O6A18eIyOrY6xsTC/+KlQyu61YRORSLA70iIksKYWe2iMgvReSgiLQleV1E5Mex6/6TiMwfaRuHQgbXdZmI9Pi+r2UjbeNQEJGZIvKsiGwXka0i8rWAY0ruO8vwukrrO1NV++P7A4Rx+kHNBqqAV4HzEo75f4AHYj/fBKwutN05uq5bgZ8W2tYhXNulwHygLcnrnwGeAARoAjYW2uYcXddlOPVEBbc1y+uaDsyP/VwNvB7wd7HkvrMMr6ukvjPbQQzmImCXqu5W1V7g18C1CcdcC/xz7Oc1wKdEREbQxqGQyXWVJKr6PHA4xSHX4vT+UlVtBiaKyPSRsW7oZHBdJYmq7lfVzbGfjwLbgdMTDiu57yzD6yopTCAGczqwx/e4k8FfsneMqvYDPUCxt3HM5LoAFsW29GtEZObImJZ3Mr32UuQSEXlVRJ4QkfpCG5MtMffsBcDGhJdK+jtLcV1QQt+ZCcRggnYCibnAmRxTbGRi82+AM1T1w8DTDOySSp1S/L4yYTMwS1U/AvwEWF9ge7JCRE4F1gJ3quqRxJcD3lIS31ma6yqp78wEYjCdgP/OuRbYl+wYEakAJlD8roC016WqXap6Mvbw50DDCNmWbzL5TksOVT2iqu/Ffn4cqBSRyQU2KyNEpBJnEV2pqusCDinJ7yzddZXad2YCMZhNwFkicqaIVOEEoR9LOOYx4Euxn6/HmWVR7Hc3aa8rwcf7WRwfajnwGHBLLDOmCehR1f2FNmq4iMg0N/YlIhfh/HvuKqxV6YnZ/Atgu6r+zySHldx3lsl1ldp3VhS9mIoJdYYW/VfgdziZP79U1a0ishxoUdXHcP4S/IuI7MLZOdxUOIszI8Pr+isR+SzQj3NdtxbM4CyItW+5DJgsIp3At4FKAFV9AHgcJytmF3AM+HJhLM2ODK7reuC/iEg/cBy4qQRuVAA+BnwR2CIir8Seuweog5L+zjK5rpL6zqzVhmEYhhGIuZgMwzCMQEwgDMMwjEBMIAzDMIxATCAMwzCMQEwgDMMwjEBMIAzDMIxATCAMwzCMQEwgjLJERMaJyHMiEhaRvxSR/bH++7tEZH2smnwo550jIlsSnhsjIm/GZjcEFp/67RnK52ZhX5WIPJ/MDsPIBhMIo1y5DVinqhHgw8A9qvpR4GxgXuy5obAbmCki/n87S4HngPuBGzOwJ2/EWrn/PoUdhpExJhBGubIY2BD7+Xzgj7Gf5+B0Cn19KCdV1SjQAZwBzs4AuAu4D6cz5+J09ojIGSLymog8KCJtIrJSRK4QkX8XkZ2xHj3Ejl0vIq2xCWVLfc9/K3aOp0RklYj8je+zUtlhGBlj21Cj7Ii5j2ar6luxp+qBX8U6bdYCfxHQhjkbtgPn4Owm7gAeU9W3Yu6jCzOwBxyhugFn97EJuBn4OE6TxHuA62LH3aaqh2NCtElE1gJnAotw5g1U4LSQbvWduy3IDsPIFhMIoxyZDLwLzpxg4GBsxgUicgvwLeBK/xtE5GlgWsC57lXVDQnPbQfmisjzOALRBKCqERHpFZHq2ESxQfb4eFNVt8Q+eyvwe1XVWHzjDN9xfyUin4v9PBM4K/Z5G1T1eOz9v/GfOIUdhpEVJhBGOXIcGBv7+cPANt9rr+K4hOJQ1SuyOP924M+Ar+H0/T/ge20McCKFPS4nfT9HfY+jxP5dishlwBXAJap6TET+EDtPJuNtg+wwjKywGIRRdqhqNxAWkbE48Yft4PXr/xLOtLzhsB1nxvdtwPfdJ0VkEnBIVftS2JMNE4DumDicQ2ynArwIXCMiY2PTy/7C/6ZkdhhGtphAGOXKkzg+/fOBW0Xkjzh++rE4LqbhsCN23hWq2uN7/nKcOQap7MmGfwMqRORPwHeAZgBV3YQzUOdVYB3QgjMXPRM7DCNjbB6EUZaIyAXAX6vqF0fwM9cB31DVHfm2R0ROVdX3RGQ88DywVFU3p7PDMLLBYhBGWaKqfxSRZ0UknO/aA/AyldYnW5TzYM8KETkPZ0f0zz5xSGmHYWSD7SAMwzCMQCwGYRiGYQRiAmEYhmEEYgJhGIZhBGICYRiGYQRiAmEYhmEEYgJhGIZhBGICYRiGYQTyfwGFRO7U00D9SAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Convert G magnitude to absolute magnitude (accounting for extinction)\n",
    "g = data2['phot_g_mean_mag'] + 5. - 5.*np.log10(1000./data2['parallax']) #- data2['a_g_val']\n",
    "'''\n",
    "# Correct B-R for reddening\n",
    "br = data2['bp_rp'] - data2['e_bp_min_rp_val']\n",
    "'''\n",
    "#g = data2['phot_g_mean_mag']\n",
    "br = data2['bp_rp']\n",
    "\n",
    "plt.plot(br, g, marker='.', linestyle='None', label='N='+str(len(data2['bp_rp'])))\n",
    "plt.xlabel(r'$(B-V)$ (mag)')\n",
    "plt.ylabel(r'$G$ (mag)')\n",
    "plt.ylim(12,-2)\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}