minimoons overlap.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt \n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lsst.sims.movingObjects as mo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>objId</th>\n",
       "      <th>a</th>\n",
       "      <th>e</th>\n",
       "      <th>inc</th>\n",
       "      <th>Omega</th>\n",
       "      <th>argPeri</th>\n",
       "      <th>meanAnomaly</th>\n",
       "      <th>H</th>\n",
       "      <th>epoch</th>\n",
       "      <th>g</th>\n",
       "      <th>sed_filename</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NESC0000001a</td>\n",
       "      <td>0.997277</td>\n",
       "      <td>0.032099</td>\n",
       "      <td>0.323426</td>\n",
       "      <td>283.584339</td>\n",
       "      <td>158.333342</td>\n",
       "      <td>306.015266</td>\n",
       "      <td>37.008092</td>\n",
       "      <td>60603.262969</td>\n",
       "      <td>0.15</td>\n",
       "      <td>S.dat</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NESC000000B3</td>\n",
       "      <td>1.011793</td>\n",
       "      <td>0.006768</td>\n",
       "      <td>0.040201</td>\n",
       "      <td>208.493888</td>\n",
       "      <td>119.170625</td>\n",
       "      <td>271.591057</td>\n",
       "      <td>37.021498</td>\n",
       "      <td>54602.716701</td>\n",
       "      <td>0.15</td>\n",
       "      <td>S.dat</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NESC000000C8</td>\n",
       "      <td>0.985949</td>\n",
       "      <td>0.027715</td>\n",
       "      <td>0.159922</td>\n",
       "      <td>155.755741</td>\n",
       "      <td>326.535489</td>\n",
       "      <td>64.935801</td>\n",
       "      <td>36.509016</td>\n",
       "      <td>55288.209681</td>\n",
       "      <td>0.15</td>\n",
       "      <td>S.dat</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>NESC000000GM</td>\n",
       "      <td>1.023415</td>\n",
       "      <td>0.019560</td>\n",
       "      <td>0.234426</td>\n",
       "      <td>97.532517</td>\n",
       "      <td>308.932078</td>\n",
       "      <td>187.279126</td>\n",
       "      <td>37.022154</td>\n",
       "      <td>58616.217523</td>\n",
       "      <td>0.15</td>\n",
       "      <td>S.dat</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NESC000000J9</td>\n",
       "      <td>1.016520</td>\n",
       "      <td>0.012195</td>\n",
       "      <td>0.059523</td>\n",
       "      <td>357.285344</td>\n",
       "      <td>68.994836</td>\n",
       "      <td>115.973768</td>\n",
       "      <td>36.829602</td>\n",
       "      <td>59660.286367</td>\n",
       "      <td>0.15</td>\n",
       "      <td>S.dat</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          objId         a         e       inc       Omega     argPeri  \\\n",
       "0  NESC0000001a  0.997277  0.032099  0.323426  283.584339  158.333342   \n",
       "1  NESC000000B3  1.011793  0.006768  0.040201  208.493888  119.170625   \n",
       "2  NESC000000C8  0.985949  0.027715  0.159922  155.755741  326.535489   \n",
       "3  NESC000000GM  1.023415  0.019560  0.234426   97.532517  308.932078   \n",
       "4  NESC000000J9  1.016520  0.012195  0.059523  357.285344   68.994836   \n",
       "\n",
       "   meanAnomaly          H         epoch     g sed_filename  \n",
       "0   306.015266  37.008092  60603.262969  0.15        S.dat  \n",
       "1   271.591057  37.021498  54602.716701  0.15        S.dat  \n",
       "2    64.935801  36.509016  55288.209681  0.15        S.dat  \n",
       "3   187.279126  37.022154  58616.217523  0.15        S.dat  \n",
       "4   115.973768  36.829602  59660.286367  0.15        S.dat  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orbitfile = 'NESCv9reintv1.TCO.withH.kep.des'\n",
    "mini = mo.Orbits()\n",
    "mini.readOrbits(orbitfile)\n",
    "mini.orbits[0:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyo = mo.PyOrbEphemerides()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lsst.sims.maf.db as db"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "opsdb = db.OpsimDatabase('kraken_2026.db')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({1: 'GalacticPlane',\n",
       "  2: 'SouthCelestialPole',\n",
       "  3: 'WideFastDeep',\n",
       "  4: 'NorthEclipticSpur',\n",
       "  5: 'DeepDrillingCosmology1'},\n",
       " {'WFD': [3], 'DD': [5], 'NES': [4]})"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "opsdb.fetchPropInfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "nightlimit = 365 * 100\n",
    "simData = opsdb.fetchMetricData(['fieldRA', 'fieldDec', 'observationStartMJD', 'rotSkyPos', \n",
    "                                 'fiveSigmaDepth', 'filter', 'seeingFwhmGeom'], \n",
    "                                sqlconstraint='night < %d and proposalId = 3' % (nightlimit))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lsst.sims.photUtils as photUtils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsst = {}\n",
    "for f in ('u', 'g' , 'r', 'i', 'z', 'y'):\n",
    "    lsst[f] = photUtils.Bandpass()\n",
    "    lsst[f].readThroughput(os.path.join(os.getenv('LSST_THROUGHPUTS_DEFAULT'), 'total_' + f + '.dat'))\n",
    "V = photUtils.Bandpass()\n",
    "#os.listdir(os.path.join(os.getenv('SIMS_MOVINGOBJECTS_DIR'), 'data'))\n",
    "V.readThroughput(os.path.join(os.getenv('SIMS_MOVINGOBJECTS_DIR'), 'data', 'harris_V.dat'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u -1.863\n",
      "g -0.369\n",
      "r 0.263\n",
      "i 0.457\n",
      "z 0.400\n",
      "y 0.407\n"
     ]
    }
   ],
   "source": [
    "colors = {}\n",
    "ctype = photUtils.Sed()\n",
    "ctype.readSED_flambda(os.path.join(os.getenv('SIMS_MOVINGOBJECTS_DIR'), 'data', 'S.dat'))\n",
    "for f in ('u', 'g' , 'r', 'i', 'z', 'y'):\n",
    "    colors[f] = ctype.calcMag(V) - ctype.calcMag(lsst[f])\n",
    "    print(f, '%.3f' %colors[f])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "v_equiv = np.zeros(len(simData), float)\n",
    "for f in ('u', 'g' , 'r', 'i', 'z', 'y'):\n",
    "    match = np.where(simData['filter'] == f)\n",
    "    v_equiv[match] = simData['fiveSigmaDepth'][match] + colors[f]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'Equivalent V band depth per filter (color corrected for S type)')"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAF3CAYAAABXB2nBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu4XFWd4P3vTyAEb1wjIokkL0Zt2ge5RMFLO1wcQF8F7NYBx7fJKL4ZbbqxR2kvw7Qw0s7o20zzSjNq04pgN4K39pH20ohCXrXllpAYMIgEEsgBDCERhJaoyO/9Y6+CyqHOteqcXbXr+3me85yqVat2/fauOudXa+2114rMRJIkNcPT6g5AkiT1joldkqQGMbFLktQgJnZJkhrExC5JUoOY2CVJahATuyRJDWJilySpQUzskiQ1iIldkqQG2bHuAKZrr732yoULF9YdhiRJs2LlypUPZOa8ieoNbGJfuHAhK1asqDsMSZJmRUTcNZl6dsVLktQgJnZJkhrExC5JUoMM7Dn2Tn77298yMjLCtm3b6g5lTHPnzmX+/PnstNNOdYciSWqgRiX2kZERnvWsZ7Fw4UIiou5wniIz2bJlCyMjIyxatKjucCRJDdSorvht27ax55579mVSB4gI9txzz77uUZAkDbZGJXagb5N6S7/HJ0kabI1L7JIkDbNGnWMf7btrN/V0e689YO+ebk+SpF6zxd5jGzZs4CUveckT988991zOPvvs+gKSJA0VE7skSQ1iYpckqUFM7D2244478vjjjz9x30vbJEmzycTeY3vvvTf3338/W7Zs4de//jXf+MY36g5JQ+KOlddzx8rr6w5DUs0mTOwRcVFE3B8Rt3R47IyIyIjYq9yPiDg/ItZFxJqIOKSt7tKIuL38LG0rPzQibi7POT8G/ELvnXbaiQ9/+MMcdthhvOENb+DFL35x3SFJkobIZC53uxi4APh8e2FELAD+PXB3W/HrgMXl5zDgU8BhEbEHcBawBEhgZURckZm/KHWWAdcB3wKOA749/V16Ul2Xp51++umcfvrptby2JGm4Tdhiz8zvA1s7PHQe8H6qRN1yAvD5rFwH7BYR+wDHAldl5taSzK8CjiuPPTszr83MpPrycGJ3uyQNt1aXvN3y0nCa1jn2iDgeuCczfzzqoX2BjW33R0rZeOUjHcolSdI0TDmxR8TTgTOBD3d6uENZTqN8rNdeFhErImLF5s2bJxOuJKmPLd+4vO4QGmc6Lfb9gUXAjyNiAzAfuCkinkvV4l7QVnc+cO8E5fM7lHeUmRdm5pLMXDJv3rxphC41j93uktpNObFn5s2Z+ZzMXJiZC6mS8yGZ+XPgCuCUMjr+cOChzLwPuBI4JiJ2j4jdgWOAK8tjD0fE4WU0/CnA13u0b9LQM+lLw2cyl7tdBlwLvCgiRiLi1HGqfwu4E1gH/D3wJwCZuRU4B7ix/HyklAG8G/hMec4d9GhEvCRJw2jCy90y860TPL6w7XYCp41R7yLgog7lK4CXPPUZkqTZ8PDV1/Cso46sOwz1SKOXbeW2Hjf+X/S63m5P6pLd7JJGa3Zir8E555zDpZdeyoIFC9hrr7049NBDOeOMM+oOS5LG9fDV1wDYcm8AE3sPrVixgq9+9ausWrWKxx57jEMOOYRDDz207rAkSUPExN5DP/zhDznhhBPYZZddAHjjG99Yc0SS1J+8fn3mmNh7qBo7KPWf1rn4/Q89rOZI1E9a3e91ayX5IxYcUWscTeGyrT306le/mn/+539m27ZtPPLII3zzm9+sOyRJ0pAxsffQy172Mo4//nhe+tKX8od/+IcsWbKEXXfdte6wJGkg2D3fG83uiq/h8rQzzjiDs88+m1/96le85jWv4X3ve9+sxyBJ/cgu99nR7MReg2XLlrF27Vq2bdvG0qVLOeSQQ+oOSZL6ii3zmWVi77EvfOELdYcgSRPql4Fzo9mq756JXRpAzjgnaSwmdkkaIhO11J2BbvA5Kl4aIi7jKjWfLXZJ0lP0csU3B8vNLlvskqS+s3zjcr8QTFOjW+y9/lBMdZRmZpKZPO1pfn+SJM0OM06Pbdiwgd/7vd/jT/7kTzjkkEPYuHFj3SGpQTxHrm706yVu6i0T+wy47bbbOOWUU1i1ahX77bdf3eFIkoaIiX0G7Lfffhx++OF1hyFJXXn46mts5Q8gE/sMeMYznlF3CJL0hDoTdLdjnRxAN3UmdkmSGsTELklSgzT6crc6FhFYuHAht9xyy6y/riRJYItdktTnnKxmakzs0gDw+nXVyZHxg6XRXfGSpHrYwq6PLXZpCNn6l5rLxC5JUoOY2CVJahATe4+98pWvrDsESdIQa/TguV6P5HzWUUdOWOdHP/pRT19TkrpRx4h2B87VyxZ7jz3zmc+sOwRJ0hBrdItdkoaV154PL1vskqQJ9cMSrnbxT46JXZKkBpkwsUfERRFxf0Tc0lb21xHx04hYExFfi4jd2h77UESsi4jbIuLYtvLjStm6iPhgW/miiLg+Im6PiC9GxJxe7qA06JxMZnzr1zxQdwhSX5lMi/1i4LhRZVcBL8nMA4GfAR8CiIgDgJOB3y/P+WRE7BAROwD/G3gdcADw1lIX4OPAeZm5GPgFcGpXeyRp6Kxf84AJXiomHDyXmd+PiIWjyr7Tdvc64M3l9gnA5Zn5a2B9RKwDXl4eW5eZdwJExOXACRFxK3AU8B9LnUuAs4FPTWdnRpvM5Wm99sgjj8z6a0qS1NKLc+zvAL5dbu8LbGx7bKSUjVW+J/BgZj42qlzSDGviinG23KUuE3tEnAk8BlzaKupQLadRPtbrLYuIFRGxYvPmzVMNV5Kkxpt2Yo+IpcAbgLdlZisZjwAL2qrNB+4dp/wBYLeI2HFUeUeZeWFmLsnMJfPmzZtu6JIaYDKtc1vvGkbTSuwRcRzwAeD4zPxV20NXACdHxM4RsQhYDNwA3AgsLiPg51ANsLuifCG4hifP0S8Fvj69XZEkSRMOnouIy4AjgL0iYgQ4i2oU/M7AVREBcF1mviszfxIRXwLWUnXRn5aZvyvb+VPgSmAH4KLM/El5iQ8Al0fEXwGrgM/2cP8kDSFb6hpmkxkV/9YOxWMm38z8KPDRDuXfAr7VofxOnhw5L0mSuuDMc5KkSat7WllNrNGLwPS6O27RgXv1dHuSps5u9omZfIebLfYe+vSnP81BBx3EQQcdxKJFizjyyNmfIEeSNNxM7D30rne9i9WrV3PjjTcyf/583vve99YdkiRpyDS6K74u73nPezjqqKN44xvfWHcoGmCzNStc63X2P/SwWXk9NZfLqvYHE3uPXXzxxdx1111ccMEFdYciSRpCJvYeWrlyJeeeey4/+MEPeNrTPMshSZp9Zp8euuCCC9i6dStHHnkkBx10EO985zvrDkmSGmX5xuV2+U+g0S322b487XOf+9ysvp4kSaPZYpckTcnDV1/jtfJ9zMQuSVKDmNglNdpklneVmqTR59glSTPPwWz9xcQuSQ3heW+BXfGShoRd8hoWttgl9T0TsjR5jU7svZ5r27m0NRtma454Sc1kV3wP/eVf/iWf+MQnnrh/5plncv7559cYkSRp2JjYe+jUU0/lkksuAeDxxx/n8ssv521ve1vNUUlS8zgSf2yN7oqfbQsXLmTPPfdk1apVbNq0iYMPPpg999yz7rAktWmdr5/tKael2WJi77F3vvOdXHzxxfz85z/nHe94R93hSNKMsuXcf+yK77E3velN/Mu//As33ngjxx57bN3hSJKGjC32HpszZw5HHnkku+22GzvssEPd4UiShkyjE3sdl6c9/vjjXHfddXz5y1+e9deWNLycdU4tdsX30Nq1a3nBC17A0UcfzeLFi+sORxp4zhYnTV2jW+yz7YADDuDOO++sOwxpWloT4zgRkzTYbLFLktQgjWuxZyYRUXcYY8rMukNQn3IqWUm90KgW+9y5c9myZUvfJs/MZMuWLcydO7fuUCSpaw7Y60+NarHPnz+fkZERNm/eXHcoY5o7dy7z58+vOwxp6K1f84Czz6mRGpXYd9ppJxYtWlR3GJK65Ej4/nfzAzdXNxYfWG8geopGdcVLkjTsTOySJDVIo7riJWkqmrDSmwPYNNqELfaIuCgi7o+IW9rK9oiIqyLi9vJ791IeEXF+RKyLiDURcUjbc5aW+rdHxNK28kMj4ubynPOjn69VkySpz02mK/5i4LhRZR8EvpeZi4HvlfsArwMWl59lwKeg+iIAnAUcBrwcOKv1ZaDUWdb2vNGvJUmSJmnCxJ6Z3we2jio+Abik3L4EOLGt/PNZuQ7YLSL2AY4FrsrMrZn5C+Aq4Ljy2LMz89qsLj7/fNu2JEka0/KNy10PvoPpDp7bOzPvAyi/n1PK9wU2ttUbKWXjlY90KJckSdPQ61Hxnc6P5zTKO288YllErIiIFf08CY0kSXWZbmLfVLrRKb/vL+UjwIK2evOBeycon9+hvKPMvDAzl2Tmknnz5k0zdElSN56YnEZ9abqJ/QqgNbJ9KfD1tvJTyuj4w4GHSlf9lcAxEbF7GTR3DHBleezhiDi8jIY/pW1bkqQxPHz1NV7qpo4mvI49Ii4DjgD2iogRqtHtHwO+FBGnAncDbynVvwW8HlgH/Ap4O0Bmbo2Ic4AbS72PZGZrQN67qUbe7wJ8u/xIGlJOJyt1Z8LEnplvHeOhozvUTeC0MbZzEXBRh/IVwEsmikPS7GgtH7v/oYfVHImk6XBKWakPuBZ7vdavecCegmmac+0a5ly7pu4w1MYpZSX1BROr1Bu22CVJahATuyRJDWJil6TCc+1qAhO7JEkNYmKXJKlBTOySpIHmKm/bM7FLktQgJnZJkhrECWqkGjnjnKReM7FLkibF5VoHg4ldUq28blzqLc+xS5LUICZ2SRowD199Td0hqI+Z2CV1NMwD+zw9MHUu3do/TOySJDWIiV2SpAYxsUuS1CAmdkmSGsTELklSg5jYJUlqEBO7JEkNYmKXJDWCa7JXTOyS1MH6NQ84UY0GkovASNKAcCpZTYaJXarBME/XqsHkkq2Dw654SZIaxMQuSVKD2BUvaUytUwb7H3pYz7ftwDRpZthilyT1xJxr17h8ax8wsUuS1CAmdkmSGsTELklSg5jYJUlqkK4Se0T8l4j4SUTcEhGXRcTciFgUEddHxO0R8cWImFPq7lzuryuPL2zbzodK+W0RcWx3uyRJ0vCadmKPiH2B04ElmfkSYAfgZODjwHmZuRj4BXBqecqpwC8y8wXAeaUeEXFAed7vA8cBn4yIHaYbl6T+56VuU+d0spqsbrvidwR2iYgdgacD9wFHAV8pj18CnFhun1DuUx4/OiKilF+emb/OzPXAOuDlXcYlSeqBmx+42elkB8y0E3tm3gOcC9xNldAfAlYCD2bmY6XaCLBvub0vsLE897FSf8/28g7PkSRJU9BNV/zuVK3tRcDzgGcAr+tQNVtPGeOxsco7veayiFgRESs2b9489aAlSWq4brriXwusz8zNmflb4J+AVwK7la55gPnAveX2CLAAoDy+K7C1vbzDc7aTmRdm5pLMXDJv3rwuQpckNdHyjctZvnF53WHUqpvEfjdweEQ8vZwrPxpYC1wDvLnUWQp8vdy+otynPH51ZmYpP7mMml8ELAZu6CIuSeqZ9WsecLCfBsq0F4HJzOsj4ivATcBjwCrgQuCbwOUR8Vel7LPlKZ8F/iEi1lG11E8u2/lJRHyJ6kvBY8Bpmfm76cYlSdIw62p1t8w8CzhrVPGddBjVnpnbgLeMsZ2PAh/tJhZpELRWS5OkmeLMc5IkNYiJXZKkBjGxS5rQHSuv9zSCNCC6OscuSVPh6PLB0c1sc3OuXcNvXnFgD6PRVNhilySpQUzskiQ1iIldkqQG8Ry7JPUxl2vVVNlilySpQUzskiQ1iIldkqQG8Ry7pBnXhOvXW/uw6MC9ao5EGp8tdmmWOHObpNlgYpckNc7yjcvrDqE2JnZJ0na6mU5W9TOxS5LUICZ2SZIaxMQuSVKDmNglSWoQE7skTUETrslXs5nYJUlqEBO7JKnn5ly7hjnXrqk7jKFkYpekPuWSrZoOE7ukSbtj5fVOjSv1OReBkTRjHGgmzT5b7JIkNYgtdmmG2XUtaTaZ2CVJgIu/NIVd8ZIkNYgtdknqM17mpm7YYpckqUFM7JIkNYiJXZKkBjGxS5IaafnG5SzfuLzuMGadg+ekGeL1683VmlFv0YF71RxJ/2stBPObVxxYcyTDo6sWe0TsFhFfiYifRsStEfGKiNgjIq6KiNvL791L3YiI8yNiXUSsiYhD2raztNS/PSKWdrtTkiQNq2674j8B/Etmvhh4KXAr8EHge5m5GPheuQ/wOmBx+VkGfAogIvYAzgIOA14OnNX6MiBpcDlPvFSPaSf2iHg28BrgswCZ+ZvMfBA4AbikVLsEOLHcPgH4fFauA3aLiH2AY4GrMnNrZv4CuAo4brpxSZI0zLppsf8fwGbgcxGxKiI+ExHPAPbOzPsAyu/nlPr7Ahvbnj9SysYqlyTNkilPJ7tl3fa32++rVt0k9h2BQ4BPZebBwL/xZLd7J9GhLMcpf+oGIpZFxIqIWLF58+apxitJ6iUTel/qJrGPACOZ2Rr6+xWqRL+pdLFTft/fVn9B2/PnA/eOU/4UmXlhZi7JzCXz5s3rInRJkppp2ok9M38ObIyIF5Wio4G1wBVAa2T7UuDr5fYVwClldPzhwEOlq/5K4JiI2L0MmjumlEnqU17KN8Rspfe9bq9j/zPg0oiYA9wJvJ3qy8KXIuJU4G7gLaXut4DXA+uAX5W6ZObWiDgHuLHU+0hmbu0yLklSt0zgA6mrxJ6Zq4ElHR46ukPdBE4bYzsXARd1E4skSXLmOUk95vXrUr2cK16S9FTT6Ya3674vmNglSTOuNWe8Zp6JXeqxO1Ze76jxIbF+zQM9P/Xw8NXX9HR7YuhWeDOxS5LUICZ2SZIaxMQuSVKDeLmbJA2xKS/+MpHWyPg9X9Db7WrSbLFLktQgttgl9YQT0ww4W9qNYWKXJD3JSWYGnl3xkiQ1iC12SVLvtbf87d6fVSZ2SdPSml1v/0MPqzmSZnDGOfWKXfGSJDWIiV2SpAYxsUvSsJulkfBzrl1T2ypvyzcuH5rFYEzsUg+5qttw8hp+9RMHz0nqiklNE3Lym1lli12ShlTP54lXX7DFLknDylnmGskWuyTVzGvY1UsmdkmSGsTELklSg3iOXX3tu2s3AfDaA/auORKpQe5ZWXcEmkG22DUQWgle/WXT+l/WHYKkUWyxqy+ZyDVoWtfzLzpwr5ojmUAftNbnXLuG37ziwLrDaCxb7Oob3127yYQuNdmWdV5iNwtssUs94FSyaunblnsftNQ1O2yxq+/Yah8s9/1sFff9bFXdYUgqTOySJDWIXfGSVJO6Zpy7+ZG7a3ldzQ4Tu2phd7s0xLasAxwVP1PsitfAcNS8pG4t37ic5RuX1x3GjOo6sUfEDhGxKiK+Ue4viojrI+L2iPhiRMwp5TuX++vK4wvbtvGhUn5bRBzbbUzqbybnwbdp/S+dnEbqU71osb8HuLXt/seB8zJzMfAL4NRSfirwi8x8AXBeqUdEHACcDPw+cBzwyYjYoQdxSTPujpXXe6mb+ts9K73Ubch0ldgjYj7wfwKfKfcDOAr4SqlyCXBiuX1CuU95/OhS/wTg8sz8dWauB9YBL+8mLknqd3UMnLv5kbsdODcEum2x/7/A+4HHy/09gQcz87FyfwTYt9zeF9gIUB5/qNR/orzDcyRJ0hRMO7FHxBuA+zOzvY8nOlTNCR4b7zmjX3NZRKyIiBWbN2+eUrySJA2Dbi53exVwfES8HpgLPJuqBb9bROxYWuXzgXtL/RFgATASETsCuwJb28pb2p+zncy8ELgQYMmSJR2Tv/qXg+akWdbH59bnXLsGwMVgZsC0W+yZ+aHMnJ+ZC6kGv12dmW8DrgHeXKotBb5ebl9R7lMevzozs5SfXEbNLwIWAzdMNy41n18QJGlsMzFBzQeAyyPir4BVwGdL+WeBf4iIdVQt9ZMBMvMnEfElYC3wGHBaZv5uBuJSg7SS+2sP2LvmSNTSmi9+nxceXHMk6ueWumZeTxJ7Zi4Hlpfbd9JhVHtmbgPeMsbzPwp8tBexaLiY4GeX165PXt+u8tYvnli+1a74XnPmOUmSGsTELkkaOk2eVtbErkZwQJ00PiemGR6u7qYZZ9LVMFu/5oHtzrPXtVRr32oN9Nv30HrjaBATuzQNzg+vvuRoeGFXvGaYrXVJml0mdkmT4qVumlH2NvSMiV2SpAYxsUuS1CAOnpOmyIFzmg5Hw/ef1rXsRyw4otY4es0Wuxrju2s3OVhP0tAzsUvqqdZiMJplAzr4bM6qu+oOoXHsipekQTagCb0jJ6vpCVvsahy75NVv1q95gI33Zt1haEjYYteMMLE2h9evDzbniB8+ttglSWoQE7t6zta6pEHStCVc7YpXY3137SZee8DePdue169ruh5ds6buEDRETOxqtFbvQS8TvNQXmjQaXj1lV7wkSQ1iYpekQTPJ1roj4oeTXfHqGQfNSeNrXcu+4HlRcyR9zolqumJil9SR169Lg8mueA0FexNm130/W+Wc8Zq0Oavucs74HrLFLkkzxMvcVAcTu7pma1jSjLhnpefZp8HErqEx3WvanZhG0zEjrXVHw2sSPMcuSbNs473pam+aMSZ2dcVueEnqLyZ2DR3Xa5f6U50j45dvXN6YxWA8xy5pO16/Lg02W+ySJDWILXZpDI6GV99wJTdNgYldknqo7klpvNRNJnZNi4PPNBmtaWX3eeHBNUfSnzbemy4Io56b9jn2iFgQEddExK0R8ZOIeE8p3yMiroqI28vv3Ut5RMT5EbEuItZExCFt21pa6t8eEUu73y1pYn45kQbAPSs9FTFF3Qyeewx4X2b+HnA4cFpEHAB8EPheZi4GvlfuA7wOWFx+lgGfguqLAHAWcBjwcuCs1pcBSRp6JjVN0bQTe2bel5k3ldsPA7cC+wInAJeUapcAJ5bbJwCfz8p1wG4RsQ9wLHBVZm7NzF8AVwHHTTcuaSrGuqbdgXPScGrCtew9OcceEQuBg4Hrgb0z8z6okn9EPKdU2xfY2Pa0kVI2VrkkDa8pttSbMmiuNUnNbw7er+ZIBlfXiT0ingl8FfjzzPxlxJgDQTo9kOOUd3qtZVTd+Dz/+c+ferDqmuelm83JabpT94j4Rmt90XG1twl1NUFNROxEldQvzcx/KsWbShc75ff9pXwEWND29PnAveOUP0VmXpiZSzJzybx587oJXZKkRupmVHwAnwVuzcy/aXvoCqA1sn0p8PW28lPK6PjDgYdKl/2VwDERsXsZNHdMKZM0Szat/6Wt9Zq40pt6rZuu+FcBfwzcHBGrS9l/BT4GfCkiTgXuBt5SHvsW8HpgHfAr4O0Ambk1Is4Bbiz1PpKZW7uIS5o2B81JGnTTTuyZ+UM6nx8HOLpD/QROG2NbFwEXTTcWqVvfXbuJ1x6wd91hSFLXnHlOk+bAOU1X02egq3PQXFNGw6t3TOyakAldUt+4Z6Uj4ydgYpeAh25bzR2P7lZ3GBpAPW+tO9OcumRil6Q+0BoZv2CCepp5rdnnjlhwRK1xTFdX17FLTbL67gdZffeDdYchTVqTz6+3ZqDT1Nlil4aY165LzWOLXeNy4JwkDRYTuzSKXfKSBpld8RpqD922euJKDVRXF/x9P1vVmGvZZ+Ta9S13wB6936yGiy12SZIaxBa7OhqGc+sTtdZX3/0gBz3fa9s1C7bcMaXqTR4NPyku4TouW+yS1Ec2bt2VjVt3rTuMvjBn1V21Xva2fOPyJ65pHyS22LWdYWipDzsvcRtsQ99a14RM7NI4WqPj7ZLvnaYvCDNTTOiaLBO7hs6wjoRX7/RsRPwUz61rFM+1d2Ri1xPshh9bE1rudsF3r87lWaXJcvCcJNVhgta6A+g0XbbYZUtdUl9rjYz/zcH7da7gGu3bscWuofHQbau7Pr/uVLPq2pY7vG59wAzaJW+22CVpAnWdWzehazpM7ENsWLrgez0KftAG0jlork84Ar5rc1bdNX53PNglj13xQ2tYkrokDRsTuzRNnm8fDj3php9ia/2uRzdz16Obeeih5/LQQ8/t/vWHyT0rn2y9Dym74ofMsLXUh3kymn7vgh+KGejsfm+M1gC6IxYcUWsck2Fil6QOnIxGg8rErkYa5pY69H9rvd/VmdTvenRzba89CCa8pl0mdqkb/TZC3oTeJ6bRBW9C77EhnrTGxD5EhuH8+rC31AfRUJxr71JrAN2uu/685kgGzJBeAueo+CHw3bWbhiKp12n13Q86Sr5LrQRfp0fXrKmlG97W+uAYhFnobLE3mMl8eNgF30ccCT8rxp2sZrQh65Y3sasRhrkL3qTeG1231GfhvLpd8poME3vDDGMrvZ+S+uq7H+ybgXSDaLbPt/es232WW+kPPfRck7vGZGKXBpQt9RrZ3T54Rs9G10XXfL9PVmNib4BhbKVDf7XU2830JXDDkNBnuuXuALlm8Jr2zkzsA8yEPlyGIaHPhjrOpc8Ez7d3qcED6vomsUfEccAngB2Az2Tmx2oOqa8NY1IftITey/Ptw5rUe9Vyr2Mhl3a21mfWlEbIt2vode59kdgjYgfgfwP/HhgBboyIKzJzbb2R9YdhTOKjDVpSbxl9bftkE/2wJvJemXYi73FrfDYSui33Slfd8u3n36eQ5Pv1XHtfJHbg5cC6zLwTICIuB04Ahiaxm7wHN3lPxXiteJP52Nonr2m13vthRHsrce+3y7zaW+Xty7sOe5Ifdv2S2PcFNrbdHwEOqymWaWkl5tcesPe4jw+rYUjanWz7+aNPKbtuVNmiec+YrXAG1q83bHji9oZye+/n7v9khRk47z2VRF13Uh9t9Bruw5ToWy13qFrvU27JT6N7vn02un5ovfdLYo8OZfmUShHLgGXl7iMRsQV4YCYDGwJ74THsBY9j9zyGveFx7F6/HsNJfTvpl8Q+Aixouz8fuHd0pcy8ELiwdT8iVmTmkpkPr7k8hr3hceyex7A3PI7dG/Rj2C+LwNwILI4XKLcPAAANtElEQVSIRRExBzgZuKLmmCRJGjh90WLPzMci4k+BK6kud7soM39Sc1iSJA2cvkjsAJn5LeBbU3zahRNX0QQ8hr3hceyex7A3PI7dG+hjGJlPGaMmSZIGVL+cY5ckST0wMIk9IhZExDURcWtE/CQi3lPK94iIqyLi9vJ797pj7VfjHMO3lPuPR8TAjgSdDeMcw7+OiJ9GxJqI+FpEuHbrOMY5jueUY7g6Ir4TEc+rO9Z+NdYxbHv8jIjIiNirrhj73Tifw7Mj4p7yOVwdEa+vO9apGJiu+IjYB9gnM2+KiGcBK4ETgf8EbM3Mj0XEB4HdM/MDNYbat8Y5hgk8DvwdcEZmrqgxzL42zjGcD1xdBoJ+HMDP4djGOY4jmfnLUud04IDMfFeNofatsY5hZq6NiAXAZ4AXA4dmZj9ek127cT6H/wF4JDPPrTXAaRqYFntm3peZN5XbDwO3Us1YdwJwSal2CdWbog7GOoaZeWtm3lZvdINhnGP4ncx8rFS7jirRawzjHMf2eXWfQYeJqlQZ538iwHnA+/H4jWuCYziwBiaxt4uIhcDBwPXA3pl5H1RvEvCc+iIbHKOOoaZhnGP4DuDbsx3PoBp9HCPioxGxEXgb8OH6Ihsc7ccwIo4H7snMH9ca1IDp8Pf8p+W00EWDdop34BJ7RDwT+Crw56O+3WuSPIbdG+sYRsSZwGPApXXFNkg6HcfMPDMzF1Adwz+tM75B0H4MqT57Z+IXoinp8Dn8FLA/cBBwH/C/agxvygYqsUfETlQH/9LM/KdSvKmcJ2mdL7m/rvgGwRjHUFMw1jGMiKXAG4C35aAMXqnRJD6LXwD+aHajGiwdjuH+wCLgxxGxgeqU0E0R8dyxtzLcOn0OM3NTZv4uMx8H/p5qBdKBMTCJPSIC+Cxwa2b+TdtDVwBLy+2lwNdnO7ZBMc4x1CSNdQwj4jjgA8DxmfmruuIbFOMcx8Vt1Y4HfjrbsQ2KTscwM2/OzOdk5sLMXEi1DschmTk8y7tNwTifw33aqr0JuGW2Y+vGII2KfzXwA+BmqhHcAP+V6nzIl4DnA3cDb8nMrbUE2efGOYY7A38LzAMeBFZn5rG1BNnnxjmG51Mdxy2l7DpHc49tnON4KvCiUnYX8K7MvKeWIPvcWMewzOLZqrMBWOKo+M7G+Ry+laobPoENwH9ujeUaBAOT2CVJ0sQGpitekiRNzMQuSVKDmNglSWoQE7skSQ1iYpckqUFM7A0TEb9rW5FodVkYZ7rb+tE0n7cwIqZ93WdE/HlEPL1D+dkR8T9HlR0UEbd2qLthJla1KjGcMYl6j3TxGtvtfzfb6iKGPyirXa2OiH0j4iul/IiI+Ebb7Vf24LX2aW1zGs9dHgO6ImH57E551bCx9nnUe7bLNGN6R0TcXKZSvSUiTuhQ58SIOGA625/gtd8QEf+919sdRib25nk0Mw9q+/nYdDeUmV3/056mPweektiBy4CTRpWdTDVDWZOMtf8zIiqj/xe8DTi3fIbuycw3d3jqEcCUPiMRsWOH4vdSze414yJih5nYzhjHcCIHAb1cDrT9PXt0osod9mE+1XS0r87MA4HDgTUdnnoi0PPEDnwTOL7Tl3pNUWb606AfqqUGO5UfRzWL1w+pJlP5Rik/m2qp1la9W4CF7dsCvgi8vq3OxVRTfS6kmtzhpvLzyvL4QuCWcnsH4K+BG6n+SfznUn4EsBz4SonrUiCA04HfUE0YcU2H/bgJOKzt/p3A4g71NgAfB24oPy8o5W+kmtRoFfBdqkWEWsfhohLTncDpbds6E7it1L+s/Xi11VkEXFv285z29wH4i7b9/+9tx+inVCsSrinH4emd9h94BPgo8GOqleP27vD6ZwP/AFwN3A7835N4/VuBT5ZjsV9b/XcCW4H15X1pfz+PAL5Ryn4O3AOsBv6AaoKjr5bXuhF4VVtsFwLfAb7QIfY7gZ3bPi/nlv1fA/xZKT+6xHlzeZ9a9ZdTTcAC1aQiN1N9hj/e/jcBfKS8768e9dovKO/rj6k+W/tTfQ7/umznZuCktn2/huqL5NpOxxA4pnwObgK+DDyzPPdlwI/K69wA7Eo1odbmcvxOolrN7qJy7FYBJ5Tn7gJcXo7HF8t+LBm1H6Pfs0ntw6htHFJi2WGc/y+vbHud1eV43dT2+GJg5QR/gx0/J+Wx84D/UPf/0UH/qT0Af3r8hsLvyh9c6+ckYC6wsfzRBdVMfVNJ7G8CLim355Rt7UKViOaW8sXAinJ7IU8mgmXAfyu3dwZWUCXBI4CHqOayflr5Z/jqUm8DsNcY+/cXwHnl9uHAjWPU2wCcWW6f0ra/u/PkxEzvBP5X23H4UYlxL6oZ5HYCDi3/GJ8OPBtYR+fEfgVwSrl9WtuxO4YqqUXZz28ArynHKHky+V3U2u7o/S/13lhu/z+t4znq9c+mShq7lPg3As+b4PUfBw4f4/hdDLy5w/t5BGN/dr7Q9h4+n2qazla9lcAuHV5nESURlPvvpvqnv2O5vwdPfn5fWMo+T7VYB5TEXvb1bqqksSPVF5wT245fx2RBlSTfVG7PLe/zHwFXUX3J2Ltsd5+y7/8GLGo7Lk8cw3Lcvw88o9z/ANViLHOovry8rJQ/u8T4n4AL2mL5H8D/VW7vBvyMKtm/F7iolB9ItdDLkg770v6eTWofRj1/B+DKUvdzlM/ceK9T7l8DHNS2D60vYxvo/DfY8XNS7r8N+NvZ/r/ZtB+74ptndFf8F4EXA+sz8/as/nr+cYrb/DZwVETsDLwO+H5WXX07AX8fETdTtU46dc8dA5wSEaup/onuSfUlAOCGzBzJaqGF1VT/KCdyOfDm0u15MlULeiyXtf1+Rbk9H7iyxPwXwO+31f9mZv46q+k376f6h/gHwNcy81dZrfp0xRiv9aq21/uHtvJjys8qqlbci3ly/zdm5r+W2/8IvHqMbf+GKiFDlSAXjlHv65n5aIn/GqqFK8Z7/bsy87oxtjUdrwUuKO/1FcCzI+JZ5bErsnP38D5Urdb2bXw6y9r2WU0P/SKqz+/PSp1LqL6ctHsZsDwzN5fnXtpW53dUXxa2U2LbNzO/Vl5rW1bz/L8auCyrRUA2Af9f2T5Un9n1bZtpP4aHU/0N/Gs5BkupWvEvAu7LzBvL6/yytX+jHAN8sDx3OdUXjeeX/fjH8tw1dO4eH20q+0DZ9u+oevbeTPWl4ryIOHsSr/UZ4O2la/8ktj811ulvcLzPyf1UX9LUhU7nu9RMY80d/Bjbj7WY+5QnZm6LiOXAsVR/uK0/1v8CbAJeWraxrcP2g+ob/JXbFUYcAfy6reh3TOLzmJkby/zX/46qVfKK8ap3uP23wN9k5hUlhrPb6owVz2TnXe5UL4D/mZl/t11htfbz6Ppjvc5vyxey0XFN9Po5wev/2xjbma6nAa8YncCrdTbGfK1H2f4zFzx1P2ISrz1enW0laU32OeNta/R+tN8P4KrMfOt2G4s4kMl9hgL4o8y8bdTzmeTzR29rLGO+7+VzdgNwQ0RcRdVyP3uC1/oqcBZVL8nKzNzS9linv8GOn5NiLtVnQl2wxT4cfgosioj9y/32fzwbqM6tERGHUHWNdnI58HaqFmwrSe9K1RJ5HPhjqq680a4E3l2WRiQiXhgRz5gg3oeBZ43z+GVU5+LuyMyRceqd1Pb72raYW4uKLJ0gDqi6Vt8UEbuUVsUbx6j3r1Q9CFB1J7ZcCbwjqvWeKSPMn1Mee35EtL6YvJVq/ANMvP9jOSEi5kbEnlRdrjdO8PrdGh3nd2hbPz0iDprENn7G9j0Q3wHe1RpkFxF7UH1+F0bEC0qdP6Zqgba7Hvh3EbFXaTm+tUOd7ZQemJGIOLG81s5l4Nb3gZMiYoeImEfVYr5hEvtyHfCqVpwR8fSIeGGJ/3kR8bJS/qyyf6OP35XAn0XJ5BFxcCn/PuUzFREvoeqOn8iU9yEinlf+B7QcRLUQz2jbxZ2Z20rsn6L6ItCu09/geJ+TFzJgK6n1IxN78+wS21/u9rHyh7cM+GZE/JDt/1i/CuxRusXeTfWPtpPvUP1z+G5m/qaUfRJYGhHXUf1BdmoJfIZqoNFNUV0C93dM3DK/EPh2RFwzxuNfpupCv3yC7ewcEdcD76HqXYCq9fHliPgBMOGKV5l5E9WApdVUx+oHY1R9D3BaRNxI9eWh9fzvUHVNXlu6/7/Ck/8Ub6U6fmuoziV/qpRPtP9juYFqZPF1wDmZee8Er9+tf6b60rM6Iv6AauDfknKp1FpgwtXtMvPfgDvakvZnqM7xromIHwP/sXx+3071vrVW4fr0qO3cB3yI6hTEj6kGdE1mCec/Bk4v78GPgOcCX6Pq7v4xVSv0/TmJZU8zczPVefPLyvauA15c/l5OAv627NNVVC3Ta4ADyvE7iWrQ5U5l328p96H6XDyzbPP9TO5LxnT2YSfg3Ij4afl/cBLV53q0y4G/iIhVbY2FS6la5N8ZVbfT3+B4n5MjqT7D6oKruw2h0gV9Rma+oe5YhlXpCv9GZr6kR9s7m2rA3rm92N5siog3AYdm5n+rOxZNT1RzO+yamX/ZVraBKSwZGxF7U101cfTMRDk8PMcuqVaZ+bVy+kADKCK+RnXZ21Fdbur5wPu6j0i22CVJahDPsUuS1CAmdkmSGsTELklSg5jYJUlqEBO7JEkNYmKXJKlB/n8qyr88baDBOwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "for f in ('u', 'g' , 'r', 'i', 'z', 'y'):\n",
    "    match = np.where(simData['filter'] == f)\n",
    "    n, b, p = plt.hist(v_equiv[match], 100, alpha=0.3, label=f)\n",
    "plt.legend(loc='upper left', fancybox=True, numpoints=1)\n",
    "plt.xlabel(\"Equivalent V band depth per filter (color corrected for S type)\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'Five Sigma Depth per filter')"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAF3CAYAAABXB2nBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu0XHWV4PHv5hEi8giESNNJ4KYxLdKM8kgDPsbhoTwcJbhaZ8E4Q0ZxMrTYYLdO2zbTwoj0tNNMMzK2dNOCQLeKijoggohCljoCJiEQIIgECORKhEACBBAksuePOgXFTd1XVd06Vae+n7XuulW/c+rUrpMLu/bv/H6/E5mJJEmqhq3KDkCSJHWOiV2SpAoxsUuSVCEmdkmSKsTELklShZjYJUmqEBO7JEkVYmKXJKlCTOySJFWIiV2SpArZpuwAWrXbbrvl0NBQ2WFIktQVy5cvfywzZ423X98m9qGhIZYtW1Z2GJIkdUVEPDiR/eyKlySpQkzskiRViIldkqQK6dtr7M288MILDA8P89xzz5UdyqimT5/OnDlz2HbbbcsORZJUQZVK7MPDw+y4444MDQ0REWWHs4XM5PHHH2d4eJh58+aVHY4kqYIq1RX/3HPPMXPmzJ5M6gARwcyZM3u6R0GS1N8qldiBnk3qdb0enySpv1UusUuSNMgqdY19pB+seqSjx3v7vrt39HiSJHWaFXuHrVmzhv322++l5+eeey5nnXVWeQFJkgaKiV2SpAoxsUuSVCEm9g7bZpttePHFF1967tQ2SVI3mdg7bPfdd+fRRx/l8ccf5/nnn+fqq68uOyQB9y2/pewQJKkrxk3sEXFxRDwaEXc22fbxiMiI2K14HhFxfkSsjoiVEXFgw76LIuLe4mdRQ/tBEXFH8Zrzo88nem+77bZ86lOf4pBDDuFd73oX++yzT9khSZIGyESmu10CfB64rLExIuYC7wAeamg+Fphf/BwCXAAcEhG7AmcCC4AElkfEVZm5sdhnMXAzcA1wDHBt6x/pZWVNTzvttNM47bTTSnlvSdJgGzexZ+aPImKoyabzgD8HrmxoWwhclpkJ3BwRMyJiD+Aw4PrM3AAQEdcDx0TEEmCnzLypaL8MOJ4OJXap0cju+L0POqSkSCRp6rR0jT0ijgN+mZm3j9g0G1jb8Hy4aBurfbhJuyRJasGkV56LiO2BM4Cjmm1u0pYttI/23oupdduz5557jhurJA26JWuXAHDY3MNKjUPd00rFvjcwD7g9ItYAc4BbI+J3qFXccxv2nQM8PE77nCbtTWXmhZm5IDMXzJo1q4XQNWjuW36LI+IlDZRJJ/bMvCMzX5OZQ5k5RC05H5iZvwKuAk4qRscfCjyZmeuA64CjImKXiNiFWrV/XbFtU0QcWoyGP4lXXrOXpoxJX1IVTWS621eBm4DXRcRwRJw8xu7XAPcDq4F/Aj4MUAyaOxtYWvx8uj6QDvhj4IvFa+7DgXOSJLVsIqPiTxxn+1DD4wROHWW/i4GLm7QvA/bb8hWSVL5NN9wIwI5HHF5yJNLEVPq2rdzT4eL/dcd29niaMnaxSxpULinbYWeffTb77LMP73jHOzjxxBM599xzyw5JUgfUK3ep11W7Yu+yZcuW8c1vfpMVK1awefNmDjzwQA466KCyw5IkDRATewf95Cc/YeHChbzqVa8C4N3vfnfJEUkaZPU57BosdsV3UG3soPqN097UzKYbbqxU97tJfnCY2DvorW99K9/5znd47rnnePrpp/nud79bdkiSpAFjV3wH/eEf/iHHHXccb3zjG9lrr71YsGABO++8c9lhSRLg8rKDotqJvYTpaR//+Mc566yzePbZZ3nb297Gxz72sa7HIGkwmbgFVU/sJVi8eDGrVq3iueeeY9GiRRx44IFlhyRpwHg9fbCZ2DvsK1/5StkhSGpTlQbNNWNlX20mdlWKo9slDToTuyQVxqvUXTde/cDpblLBal9SFZjYJWmSevEavAPmVGdil6QB5ZeBaqr0NfZO/9FOdgRpZpKZbLWV358kSd1hxumwNWvW8PrXv54Pf/jDHHjggaxdu7bskAaC672rXb3YvS61wsQ+Be655x5OOukkVqxYwV577VV2OJKkAWJinwJ77bUXhx56aNlhSJpCVbv7m6qj0tfYy/LqV7+67BAkTVC/J+d2xxK5Cl31WLFLklQhJnZJkiqk0l3xZXQtDQ0Nceedd3b9fSVJAit2SRK1a+0uWFMNJnb1Peevq0yOjlevqXRXvCRVldW1RmPFLjVwBTtJ/c7ELklShZjYJUmqEBN7h735zW8uOwRJ0gCr9OC5To9U3fGIw8fd56c//WlH31PS1OjnkewOnNNYrNg7bIcddig7BEnSADOxSxo4/VytS+MxsUtSB1Tly4Ld/P3PxC5JUoWMm9gj4uKIeDQi7mxo+9uI+HlErIyIb0fEjIZtn4yI1RFxT0Qc3dB+TNG2OiL+oqF9XkTcEhH3RsTXImJaJz+gqsvFZCbugZWP8cDKx8oOQ1IXTKRivwQ4ZkTb9cB+mfkG4BfAJwEiYl/gBOAPitd8ISK2joitgb8HjgX2BU4s9gX4LHBeZs4HNgInt/WJJI3K5C5V37jT3TLzRxExNKLt+w1PbwbeWzxeCFyemc8DD0TEauDgYtvqzLwfICIuBxZGxN3AEcC/L/a5FDgLuKCVDzPSRKanddrTTz/d9feUJKmuE9fYPwhcWzyeDaxt2DZctI3WPhN4IjM3j2iXSlXlbn675aVqayuxR8QZwGbgy/WmJrtlC+2jvd/iiFgWEcvWr18/2XAlSaq8lleei4hFwLuAIzOznoyHgbkNu80BHi4eN2t/DJgREdsUVXvj/lvIzAuBCwEWLFgw6hcASTUTqczr+8x7w25THY6kLmipYo+IY4BPAMdl5rMNm64CToiI7SJiHjAf+BmwFJhfjICfRm2A3VXFF4Ibefka/SLgytY+iiRJGrdij4ivAocBu0XEMHAmtVHw2wHXRwTAzZl5SmbeFRFfB1ZR66I/NTN/WxznI8B1wNbAxZl5V/EWnwAuj4jPACuAizr4+SSNwuvsUjVNZFT8iU2aR02+mXkOcE6T9muAa5q038/LI+clSVIbXHlOkjpk0w03VmZpWfWvSt+2tdNdjQ4uUr+wm705k64GgRV7B/3DP/wD+++/P/vvvz/z5s3j8MO7v0COJGmwmdg76JRTTuG2225j6dKlzJkzhz/7sz8rOyRJ0oCpdFd8WU4//XSOOOII3v3ud5cdSiV1c0W4+nvtfdAhXXtPaTTeUlUTYWLvsEsuuYQHH3yQz3/+82WHIkkaQCb2Dlq+fDnnnnsuP/7xj9lqK69ySJK6z+zTQZ///OfZsGEDhx9+OPvvvz8f+tCHyg5JkiZtydoldvv3sUpX7N2envalL32pq+8nSdJIVuyS1GHOl1eZTOySJFWIiV0SUFutzhXrpP5nYpekPuBgNk2UiV3SQPC6twaFiV3SK9gdL/W3Sk93kwaJCVkSVDyxd3pNcdcLL1c314iXpH5lV3wH/dVf/RWf+9znXnp+xhlncP7555cYkSRp0JjYO+jkk0/m0ksvBeDFF1/k8ssv5/3vf3/JUUlSa1xatj9Vuiu+24aGhpg5cyYrVqzgkUce4YADDmDmzJllhyVNWv16fbeXZZbUPhN7h33oQx/ikksu4Ve/+hUf/OAHyw5HUp+zYtZk2RXfYe95z3v43ve+x9KlSzn66KPLDkeSNGCs2Dts2rRpHH744cyYMYOtt9667HAkSQOm0om9jOlpL774IjfffDPf+MY3uv7ekrbkinMaNHbFd9CqVat47Wtfy5FHHsn8+fPLDkcDwpu3SGpU6Yq92/bdd1/uv//+ssPQFKgvjuMiRZJ6nRW7JEkVUrmKPTOJiLLDGFVmlh1C33EpWUmauEpV7NOnT+fxxx/v2eSZmTz++ONMnz697FAkTbFNN9zowD2VolIV+5w5cxgeHmb9+vVlhzKq6dOnM2fOnLLDkCbEFeik/lOpxL7tttsyb968ssOQusKR8NXminNqVaW64iVJGnQmdkmSKqRSXfGSpkY/Xmt34JoG1bgVe0RcHBGPRsSdDW27RsT1EXFv8XuXoj0i4vyIWB0RKyPiwIbXLCr2vzciFjW0HxQRdxSvOT96ea6aJEk9biJd8ZcAx4xo+wvgh5k5H/hh8RzgWGB+8bMYuABqXwSAM4FDgIOBM+tfBop9Fje8buR7SZKkCRo3sWfmj4ANI5oXApcWjy8Fjm9ovyxrbgZmRMQewNHA9Zm5ITM3AtcDxxTbdsrMm7I2+fyyhmNJknqAI/T7S6uD53bPzHUAxe/XFO2zgbUN+w0XbWO1DzdplyRJLej0qPhm18ezhfbmB49YHBHLImJZLy9CI0lSWVpN7I8U3egUvx8t2oeBuQ37zQEeHqd9TpP2pjLzwsxckJkLZs2a1WLoktTb7PpWO1pN7FcB9ZHti4ArG9pPKkbHHwo8WXTVXwccFRG7FIPmjgKuK7ZtiohDi9HwJzUcS5Ja4lQ3DbJx57FHxFeBw4DdImKY2uj2vwG+HhEnAw8B7yt2vwZ4J7AaeBb4AEBmboiIs4GlxX6fzsz6gLw/pjby/lXAtcWPpFG4lKyksYyb2DPzxFE2Hdlk3wROHeU4FwMXN2lfBuw3XhxSL7hv+S3sfdAhZYchSaNy5Tn1LO/D3nv6cQW6stUvC+x4xOElR6JBYWKX+oRd8JImwpvASJJUISZ2SZIqxMQuadK8LCD1LhO7JEkVYmKXJKlCTOySpHEtWbvEpW77hIldkqQKMbFLklQhJnb1JFedk6TWuPKcJPUIr2GrE6zYpT7gvHFJE2VilySpQuyKl1QZ9TupSYPMil2apPuW3+LgPmqXB7xEMHF+6VC3mNglSaoQE7skSRViYpckqUJM7JIkVYiJXZKkCjGxS5JUISZ2SZIqxMQuSZow17PvfSZ2SW1xoRqpt5jYJVWCK7tJNa4Vr57iUq0aRHZvq5Os2CVJqhATuyRJFWJXvNSi+mWDvQ86ZMrew0FpkibLil2SumTTDTc6yE9TzsQuSVKFmNglSaoQE7skSRViYpckqULaSuwR8acRcVdE3BkRX42I6RExLyJuiYh7I+JrETGt2He74vnqYvtQw3E+WbTfExFHt/eRJEkaXC0n9oiYDZwGLMjM/YCtgROAzwLnZeZ8YCNwcvGSk4GNmfla4LxiPyJi3+J1fwAcA3whIrZuNS6pClx/feIcaS69Urtd8dsAr4qIbYDtgXXAEcAVxfZLgeOLxwuL5xTbj4yIKNovz8znM/MBYDVwcJtxSVJfcDlZdVrLiT0zfwmcCzxELaE/CSwHnsjMzcVuw8Ds4vFsYG3x2s3F/jMb25u8RpIkTUI7XfG7UKu25wG/C7waOLbJrll/ySjbRmtv9p6LI2JZRCxbv3795IOWJKni2umKfzvwQGauz8wXgG8BbwZmFF3zAHOAh4vHw8BcgGL7zsCGxvYmr3mFzLwwMxdk5oJZs2a1EbokqVVL1i7xEkIPayexPwQcGhHbF9fKjwRWATcC7y32WQRcWTy+qnhOsf2GzMyi/YRi1Pw8YD7wszbiklQCB/xJvaHlm8Bk5i0RcQVwK7AZWAFcCHwXuDwiPlO0XVS85CLgnyNiNbVK/YTiOHdFxNepfSnYDJyamb9tNS5JkgZZW3d3y8wzgTNHNN9Pk1Htmfkc8L5RjnMOcE47sai/1e+UJklqjyvPSZJUISZ2SZIqxMQutem+5bd4KUFSz2jrGrukznJU+WDYdMONLJ/fbAkPqX1W7JIkVYiJXZKkCjGxS5JUISZ2SX3L27VKWzKxS5JUISZ2SZIqxMQuSVKFOI9d6gFVmr9e/yzz3rBbyZFIg8nErlK5YpskdZZd8ZKklixZu4Qla5eUHYZGsGKXpC6647E7ag/mv6HcQFRZVuySJFWIiV2SpAoxsUuSVCEmdkmSKsTELmlKVGluvtRPTOySJFWIiV2SSjDtppVMu2ll2WGogpzHLqnveLtWaXRW7FKH3Lf8FpfIlVQ6E7tUMgeZSeokE7skSRXiNXaVwi5rSZoaVuyS1CUv3QBGmkImdkmSKsSueEl9w2lu0vis2CVJqhATuyRJFWJilySpQkzskqS2LFm7pOwQ1MDErq5zDvvgeGDlY66sNw5vBqNOayuxR8SMiLgiIn4eEXdHxJsiYteIuD4i7i1+71LsGxFxfkSsjoiVEXFgw3EWFfvfGxGL2v1QkiQNqnYr9s8B38vMfYA3AncDfwH8MDPnAz8sngMcC8wvfhYDFwBExK7AmcAhwMHAmfUvA1KVWc1KmgotJ/aI2Al4G3ARQGb+JjOfABYClxa7XQocXzxeCFyWNTcDMyJiD+Bo4PrM3JCZG4HrgWNajUuSpEHWzgI1vwesB74UEW8ElgOnA7tn5jqAzFwXEa8p9p8NrG14/XDRNlq7JFWCS8kW7rm29vt1x77ysTqqna74bYADgQsy8wDgGV7udm8mmrTlGO1bHiBicUQsi4hl69evn2y8kqReUE/qmhLtJPZhYDgz60Ocr6CW6B8putgpfj/asP/chtfPAR4eo30LmXlhZi7IzAWzZs1qI3RJkqqp5cSemb8C1kbE64qmI4FVwFVAfWT7IuDK4vFVwEnF6PhDgSeLLvvrgKMiYpdi0NxRRZvUl+5bfotT+qRG91xrld5F7d4E5k+AL0fENOB+4APUvix8PSJOBh4C3lfsew3wTmA18GyxL5m5ISLOBpYW+306Mze0GZckqQwm8NK1ldgz8zZgQZNNRzbZN4FTRznOxcDF7cQiSZK8bavUdc5dlzSVTOySpPa12gXvtLeOc614SeoBrhevTrFiV9c4Unxw1S8/zHvDbi0fY9MNN3YqHE2BJRtXAXDYLvuWHIms2CVJqhATuyRJFWJilySpQkzskjSFvAHMBLmwTceY2CVJqhBHxUtd4sI0qpxfLq/9nn1QuXHoFazYJUntqSd49QQTuyRJFWJXvCSpN4wcQOcysy0xsUtT5L7lt7D3QYeUHUbfc8U5aXLsipckqUJM7JIkVYhd8ZLUI+p3ePvNm95QciTjGGMUvDeDKZ+JXVPOu7qprhN3eZM0NhO7NMVcmEZqUX2UvKPjJ8XELklTwDXiVRYTuyRpYlxhri84Kl5Sz3IOuzR5JnZJkirExC5JUoWY2CX1hdvXPlF2CIPN6+t9w8QuTSGnunXW7WufMMFL43BUvKSum8xCNYOYyKfdtLJ3Vp+zUu87VuySes4PVj1SdgjqJSNv56oxWbFryriUrMbzwMrHXF62V1mp9y0TuzSF1v1iBQB7/P4BJUfSf36w6hGmD2A3vNQuu+IlSaoQK3ZJPWf60p+WHULLXCNeZTOxSypVs4Fy00uIQz3OO71NmF3xkvqO89l735KNq8oOYWC1ndgjYuuIWBERVxfP50XELRFxb0R8LSKmFe3bFc9XF9uHGo7xyaL9nog4ut2YpLI98sBTPPLAU2WH0dNuf+gJbn/I5Cx1Wicq9tOBuxuefxY4LzPnAxuBk4v2k4GNmfla4LxiPyJiX+AE4A+AY4AvRMTWHYhLJblv+S1OdZP6mVPd+lpbiT0i5gD/Fvhi8TyAI4Aril0uBY4vHi8snlNsP7LYfyFweWY+n5kPAKuBg9uJS1J/mr70pw6ck9rUbsX+v4E/B14sns8EnsjMzcXzYWB28Xg2sBag2P5ksf9L7U1eI0mSJqHlxB4R7wIezczGPptosmuOs22s14x8z8URsSwilq1fv35S8UqSNAjame72FuC4iHgntdkpO1Gr4GdExDZFVT4HeLjYfxiYCwxHxDbAzsCGhva6xte8QmZeCFwIsGDBgqbJX1Jvc8DcxEy7aSVAd28G47X1Smi5Ys/MT2bmnMwcojb47YbMfD9wI/DeYrdFwJXF46uK5xTbb8jMLNpPKEbNzwPmAz9rNS5Jg8Mpb9KWpmKBmk8Al0fEZ4AVwEVF+0XAP0fEamqV+gkAmXlXRHwdWAVsBk7NzN9OQVxSaVwzfurUk/sb584oOZI+ZqVeKR1J7Jm5BFhSPL6fJqPaM/M54H2jvP4c4JxOxCKVybnrk/f0/ZsA2OH3dmzrOCb4AeEKdONy5TlJkirExC5JmhJLNq5yadkSeBMYSZVy+9onut4d78I06iVW7JK6Zqypbs9f9zOe3GCtoQmqX2vXFvyvSB3lGvFSH3E0fCWZ2CVNORelkbrHxC51gNPcpBI49a0pr7FLklQhJnZJkirErnh1hIPm1Kpt77277BA0xZZsXMVhu+xbdhgDw4pdUuXcvvYJbxCjgWVil7po3S9WvHRDGGki6rdv7TinulWWXfGSNEiqmNAdHf8KVuySesqTG7bp2Ap0dslrEFmxS21yDvvoqr4wjWvEqxdZsUuSVCEmdklTourVuibHW7h2j13xaovz19WObs1hL+NWrlJZTOySBkJ9EN3AJvgqjoZXU3bFS5JUISZ2SWpBX42It1ofKHbFS+qoTg2aq89l33nXzR05ngaAC9UAJnapZc5fl9SL7IqXSuCa8eXpx9Xopt20curWjFflWLFL6jpv1SpNHSt2SR1x+0NPuCiNekP9WvuAsmJXS1yYRq3opUq91XntjoZXr7Nil9TTOnWnN2lQmNgltc0ueKl3mNglDbR+GiFfhdHx3gxm6tnHJU2S89cl9TIrdkmSKsSKXZPiaHipDzgafqBZsUvqil6a6taqvprqpoFlYpdK1O/LynZrUZonN2zjtDdpglpO7BExNyJujIi7I+KuiDi9aN81Iq6PiHuL37sU7RER50fE6ohYGREHNhxrUbH/vRGxqP2PJUkT14/rx2sc91w7sCvQtVOxbwY+lpmvBw4FTo2IfYG/AH6YmfOBHxbPAY4F5hc/i4ELoPZFADgTOAQ4GDiz/mVAkjQJv1zu9XW1ntgzc11m3lo83gTcDcwGFgKXFrtdChxfPF4IXJY1NwMzImIP4Gjg+szckJkbgeuBY1qNS1PjvuW3OHBOlWfl3j3OZZ86HbnGHhFDwAHALcDumbkOaskfeE2x22xgbcPLhou20dolqSfc8dgdPTVwrukiNVbqKrQ9GiUidgC+CXw0M5+KiFF3bdKWY7Q3e6/F1Lrx2XPPPScfrNQGF6Z52WQGzFVhNLz6WP06++uOLTeOLmqrYo+Ibakl9S9n5reK5keKLnaK348W7cPA3IaXzwEeHqN9C5l5YWYuyMwFs2bNaid0SZIqqZ1R8QFcBNydmX/XsOkqoD6yfRFwZUP7ScXo+EOBJ4uu+uuAoyJil2LQ3FFFm9QTHnngKav1HuG0N2l87fwX8hbgPwJ3RMRtRdtfAn8DfD0iTgYeAt5XbLsGeCewGngW+ABAZm6IiLOBpcV+n87MDW3EpQ5z0Jwk9Y+WE3tm/oTm18cBjmyyfwKnjnKsi4GLW41Fkjrp9rVP8Ma5M8oOQ2qJfVpSyeqrz+3x+weUHMnEDNLAuV4aCS9NlEvKSpqQbi0fK02JAVqFzopdo/LaulrV75V6X3H+ukawYpfUdxwdXw1LNq5yBbop4H8ZktTEVT+/GYB5u7265Eiaq68+9xvX6tIIVuzSGJy/LqnfWLFLGpMD5qT+YsUuSWN44LFneOCxZ8oOQ5owK3ZtwdHw5XTB99t89mYcDd9Fj6+u/d5zr3LjUM+xYpckqUKs2PUSK3WNNIjX1x985udN2x947JneGSFfr9Y1OQNyC1crdkl9y/nsMG3Fg2WH0Dbns3fWYP8XIY3g9LYaK3Wpf1mxSz1m3S9WvDSQTt0x0aTuCHn1Ayt2eW1dHeGI+C7w2npnVPxau4ldwi74ula64Acxoder9p4ZTCc1MLFL6ntPbtiGnXfdXHYYU8dKXZPgNXZJwGAOmKuKaSserMTo+K6759pK3qfdxD7gvL6uQfbgMz9vazT8lA+ke3z1QFXrTnnrDBO7pJb18/V1p7epqrzGPqCs1GscNGcXfLumbCBdC5V6vTv+Nwe4fvyk3HNtpUbIm9ilAWVCl6rJxC5p0nqxC76+tOx4o+Onqgu+Y5X7AF1T7ykVmttuYh8wdsG/rNe74afyNq5W6z3IhA68PIDusF32LTmS/mVilzRhvVipS3olE/uAsFJ/Wa9X6r2q35P6II2Cn7biQQfQDTATuwaGCX1wdTOpt3Tfdrvhe0cFrrWb2CvOSr3/rfvFio5dZx+Ea+sTHUQnjamPp8CZ2KU+0O5AunYSej93wZfV/e5NYtq3ZOMqB9C1yMSuyhvkLvhBqND72hR2wbtYTQf0abe8ib3CBr0bfpATeqf0a7W+/vlh1q8DmMG2O/XglxuvqWsKmdgraNATepVNpEt+ULvdR/PCUzMASknwW3TJm9DVBSZ2VZLVunrF9pvWQGxXdhh9qWcWqxl5a9ce75o3sVeIlfrgJPSRlXsnrqX3a7W+/vnhskMY1yNPPc/uO5WT3J3TPnhM7BVgQh+chD7SICf0yXrhqe5eb99+05quvZe6rMenwvVMYo+IY4DPAVsDX8zMvyk5pJ5nQh/chP7A+meKR3cC8Oq5+03q9VVJ5r1arTdL6o889TxAKZW7I+SnQA+PmO+JxB4RWwN/D7wDGAaWRsRVmbmq3Mh606An9EFL5i8n8fZUJZlD6wm9kwPpWq3I6wkeup/k+zXB98y19mYar7/3SJLvicQOHAyszsz7ASLicmAhMLCJfdCTd6NBSeTtJPBn1t7J1uvWArDjLkMdiqg89cQ9a7s5U1KV1xN83UQS/VR0rZdZxau6eiWxzwbWNjwfBg4pKZYpZ9KuXrKuJ+V5s5qvNNapqruunsSb2bRxDdCbCX6ySbrTSX27Xz/afMOvX344ffuHO/qeE9HtKr5euUOteu+XSr5euUOPV+8lV+69ktijSVtusVPEYmBx8fTpiLhnSqMqz27AY2UHUSGez87yfHae57Szqno+J/TNq1cS+zAwt+H5HGCLr82ZeSFwYbeCKktELMvMBWXHURWez87yfHae57SzBv18blV2AIWlwPyImBcR04ATgKtKjkmSpL7TExV7Zm6OiI8A11Gb7nZxZt5VcliSJPWdnkjsAJl5DXBN2XH0iMpfbugyz2eH4ZlGAAAH9klEQVRneT47z3PaWQN9PiNzizFqkiSpT/XKNXZJktQBJvYSRcTciLgxIu6OiLsi4vSi/X3F8xcjYmBHdrZijHP6txHx84hYGRHfjogZ4x1LY57Ps4tzeVtEfD8ifrfsWPvBaOezYfvHIyIjYreyYuwnY/x9nhURvyz+Pm+LiHeWHWs32RVfoojYA9gjM2+NiB2B5cDx1Obwvwj8I/DxzFxWYph9ZYxzOge4oRio+VmAzPxEiaH2hTHO53BmPlXscxqwb2aeUmKofWG085mZqyJiLvBFYB/goMys4jzsjhrj7/PfAU9n5rmlBlgSK/YSZea6zLy1eLwJuBuYnZl3Z2ZVF9+ZUmOc0+9n5uZit5upJXqNY4zz2bh04KtpsqCUtjTa+Sw2nwf8OZ7LCRvnfA4sE3uPiIgh4ADA9WY7ZIxz+kHg2pH7a2wjz2dEnBMRa4H3A58qL7L+1Hg+I+I44JeZeXupQfWxJv+9f6S4XHRxROxSWmAlMLH3gIjYAfgm8NERlZBaNNo5jYgzgM3Al8uKrR81O5+ZeUZmzqV2Lj9SZnz9pvF8Uvt7PAO/HLWsyd/nBcDewP7AOuB/lRhe15nYSxYR21L7g/xyZn6r7HiqYLRzGhGLgHcB708Hl0zYBP5GvwL8UXej6l9NzufewDzg9ohYQ+0y0a0R8TvlRdk/mv19ZuYjmfnbzHwR+CdqdxAdGD2zQM0giogALgLuzsy/KzueKhjtnEbEMcAngH+Tmc+WFV+/GeN8zs/Me4unxwE/LyO+ftPsfGbmHcBrGvZZAyxw8Nz4xvj73CMz1xVP3wPcWUZ8ZXFUfIki4q3Aj4E7qI2CB/hLYDvg/wCzgCeA2zLz6FKC7DNjnNPzqZ3Xx4u2mx3FPb4xzufJwOuKtgeBUzLzl6UE2UdGO5/Fypv1fdZgYp+QMf4+T6TWDZ/AGuC/NCT6yjOxS5JUIV5jlySpQkzskiRViIldkqQKMbFLklQhJnZJkirExC5NUET8tuFuUbdFxFBELIiI8ztw7O0j4ssRcUdE3BkRPylW0yIiftp+9JOKZU0Rxx0RsSoiPhMR27VxvI9GxPYNz5/uTKSTiuFfF3f/ui0iZkfEFUX7YRFxdcPjN3c7NqnTXKBGmrhfZ+b+I9rWAJ24+97pwCOZ+a8AIuJ1wAsAmVlGsjk8Mx8rvlxcWPwsavFYHwX+BejKwkDFoiVRrDpW937g3Mz8UvH8vU1eehjwNDDhL1IRsU3DzYWknmDFLrWhXvFFxFZFpTujYdvqiNg9ImZFxDcjYmnx85Ymh9oDeGmBl8y8JzOfL47zdPF7q4j4QlF5Xh0R10TEe4ttayLiryPipohYFhEHRsR1EXFfRJxS7LNDRPwwIm4tqvGF432+zHwaOAU4PiJ2LY7zX4vPsTIi/nvRNhS1+91fWrRfUfRCnAb8LnBjRNzYcG7OiYjbI+LmiNi9yXk9KyL+OSJuiIh7I+I/N2wb7f3vjogvALcCcxv2/xC123h+qugVGYqIO0e831DxOf+0qOr/9Wj/bkVsF0bE94HLxjuHUtdlpj/++DOBH+C3wG3Fz7eLtsOAq4vHnwM+UDw+BPhB8fgrwFuLx3tSW/5y5LH3Bx4FbgI+A8xv2PZ08fu9wDXUvpD/DrAReG+xbQ3wx8Xj84CVwI7UVi98tGjfBtipeLwbsJpikaoRsawBdhvRdlvxmY6iVr1HEcfVwNuAIWqrfL2l2P9i4OPNjlfs9+7i8f8E/luTGM4CbgdeVcS6ltoXhLHe/0Xg0FH+7S5pOFdDwJ1N/v3Oqsc81r9bsd9y4FVl/03640+zH7vipYlr1hXf6GvU7tD1JeCE4jnA24F9az3EAOwUETtm7f7RAGTmbRHxe9QS19uBpRHxpsy8u+H4bwW+kbUu5l81VsCFq4rfdwA7FMffFBHPFT0JzwB/HRFvo5YEZwO7A7+awGevB39U8bOieL4DMB94CFibmf+vaP8X4DTg3CbH+g21hAy1BPmOUd7zysz8NfDr4rMeTO0cjPb+D2bmzRP4LBPV9N+teHxVEZvUc0zsUufcBLw2ImYBx1OrvKFWWb5pvESQtW7vbwHfiogXgXcCjYk9mr7wZc8Xv19seFx/vg2168yzgIMy84WorUk+fbwPVSSzIeAXRQz/IzP/ccQ+Q9Qq8Vd8pFEO+UJm1rf9ltH/P9TseGO9/zOjHKdVTf/dikTf6feSOsZr7FKHFMnq28DfUeu2rd9w5vs03K88Irao+iPiLRGxS/F4GrAvtZurNPoJ8EfFtfbdqXUjT8bO1LrlX4iIw4G9xntBMXjuC8D/zcyNwHXABxtG7M+OiPqdyfaMiDcVj08s4gXYRO2ywGQtjIjpETGT2mddOs77t2tknOP+u0m9yIpd6qyvUUtA/6mh7TTg7yNiJbX/5n5EbaBWo72BC4oR3VsB36V2j+lG3wSOpHYLyl8AtwBPTiK2LwPfiYhl1K6Zj3Wr1RsbYvk2cDZAZn4/Il4P3FRUrk8D/4Fa5X03sCgi/hG4F7igONaFwLURsS4zD59EvD+jdh72BM7OzIeBh8d4/3Z9B7iiGFT4J0zs303qOd7dTeojEbFDZj5dVLE/ozZYbSLXyKc6riFqg9D269DxzqI2aLDZNXpJY7Bil/rL1cVAuGnUqtjSk7qk3mLFLklShTh4TpKkCjGxS5JUISZ2SZIqxMQuSVKFmNglSaoQE7skSRXy/wHuyX4V3jrhpwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "for f in ('u', 'g' , 'r', 'i', 'z', 'y'):\n",
    "    match = np.where(simData['filter'] == f)\n",
    "    n, b, p = plt.hist(simData['fiveSigmaDepth'][match], 100, alpha=0.3, label=f)\n",
    "plt.legend(loc='upper left', fancybox=True, numpoints=1)\n",
    "plt.xlabel(\"Five Sigma Depth per filter\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "sortIdx = np.argsort(simData['observationStartMJD'])\n",
    "simData = simData[sortIdx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3649"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tstep = 1.0\n",
    "times = np.arange(int(simData['observationStartMJD'][0]) + 0.56, \n",
    "                  int(simData['observationStartMJD'][-1]) + 0.56 + tstep, \n",
    "                  tstep)\n",
    "len(times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Fraction of population/dates')"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg4AAAF3CAYAAADNfsOtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcFdW97/3Pl0FxREXOjYIGVIzBCbHF8cQxESdIHE4wJgGHS46PxtzcmCgxMSec63nCNdEcoybHEwf0iYJDoh0l1+gxPmYEW0EiGGKrqMQMRJQ4BBT83T9qNW6bPVRDV/fevb/v12u/umrVqlW/XQ29f3vVqlWKCMzMzMzy6NfbAZiZmVnjcOJgZmZmuTlxMDMzs9ycOJiZmVluThzMzMwsNycOZmZmlpsTBzMzM8vNiYOZmZnlVmjiIGm8pCWS2iVdXGb7ppJmp+1zJY0o2TYtlS+RdGwX2vyOpNfzHMPMzMy6prDEQVJ/4BrgOGA0cLqk0Z2qnQ28EhG7AVcCM9K+o4FJwJ7AeOBaSf1rtSmpBdgmzzHMzMys6wYU2PY4oD0ingWQNAuYCCwuqTMR+Je0fCdwtSSl8lkRsRp4TlJ7ao9Kbaak4nLgE8DHah0jqsy1vf3228eIESM25D2bmZk1pMcee+yvETG0Vr0iE4dhwIsl68uAAyvViYg1klYCQ1L5bzrtOywtV2rzfKA1Iv6Y5R41j/HXSoGPGDGCtra2Wu/PzMysz5D0fJ56RSYOKlPW+Vt+pTqVystdWglJOwKnAUdsYBxImgpMBdh5553L7GJmZmZFDo5cBuxUsj4ceKlSHUkDgMHAiir7VirfD9gNaJe0FNg8Xd6odoz3iIjrIqIlIlqGDq3ZU2NmZtaUikwcHgVGSRopaROywY6tneq0ApPT8qnAQ2nsQSswKd0RMRIYBcyr1GZE3BcR74uIERExAngzDYasdgwzMzProsIuVaTxBOcD9wP9gRsiYpGk6UBbRLQC1wO3pN6BFWSJAKne7WQDKdcA50XEWoBybdYIpewxuurtt99m2bJlrFq1akN2tzo0aNAghg8fzsCBA3s7FDOzhiF/+V5fS0tLdB4c+dxzz7HVVlsxZMgQOg2+tAYUEbz88su89tprjBw5srfDMTPrdZIei4iWWvU8c2ROq1atctLQh0hiyJAh7kEyM+siJw5d4KShb/Hv08ys65w4NBBJfOpTn1q3vmbNGoYOHcqJJ54IQGtrK9/4xjeqtvHSSy9x6qmnFhpnHkuXLmWzzTZjzJgx614333xzt7X/8MMPrzsvZmbWfYqcx6FPu3XuC93a3icOrD13xBZbbMGTTz7J3//+dzbbbDMeeOABhg0btm77hAkTmDBhQtU2dtxxR+68886Njrc77LrrrixYsKC3wzAzsy5wj0ODOe6447jvvvsAuO222zj99NPXbbvppps4//zzAZgyZQoXXHABhxxyCLvsssu6ZGHp0qXstdde6+p/9KMf5aSTTmLkyJFcffXVXHHFFey3334cdNBBrFiRTXexYMECDjroIPbZZx8+9rGP8corr1QtP+KII7jooosYN24cu+++Oz//+c+79B633HJLvvCFLzB27FiOPvpoli9fXvV47e3tHHPMMey7776MHTuWZ555BoDXX3+dU089lT322IMzzjgDDwQ2M9t4ThwazKRJk5g1axarVq1i4cKFHHhg51m83/XHP/6RX/ziF9x7771cfPF6DxIF4Mknn+TWW29l3rx5XHLJJWy++ebMnz+fgw8+eN2lg09/+tPMmDGDhQsXsvfee/P1r3+9ajlkl1HmzZvHt7/97feUl3rmmWfec6miI8F44403GDt2LI8//jiHH354zeOdccYZnHfeeTzxxBP86le/YocddgBg/vz5fPvb32bx4sU8++yz/PKXv+zKqTYzszJ8qaLB7LPPPixdupTbbruN448/vmrdj370o/Tr14/Ro0fz5z//uWydI488kq222oqtttqKwYMHc9JJJwGw9957s3DhQlauXMmrr77K4YcfDsDkyZM57bTTKpZ3OPnkkwHYf//9Wbp0adljV7pU0a9fPz7+8Y8D8MlPfpKTTz654vFee+01/vCHP/Cxj2XPNRs0aNC6dsaNG8fw4cMBGDNmDEuXLuWwww6res7MzKw6Jw4NaMKECVx44YU8/PDDvPzyyxXrbbrppuuWK3XTl9bp16/fuvV+/fqxZs2aDY6xo53+/ftvVDtQ/e6HapcfSt9bd8RhZma+VNGQzjrrLC699FL23nvvwo81ePBgtt1223WXEW655RYOP/zwiuXd4Z133lk3JuPWW2/lsMMOq3i8rbfemuHDh3P33XcDsHr1at58881uicPM+pC2G7OXbTT3ODSg4cOH87nPfa7Hjjdz5kz++Z//mTfffJNddtmFG2+8sWp5Xh1jHDqcddZZXHDBBWyxxRYsWrSI/fffn8GDBzN79uyqx7vlllv4zGc+w6WXXsrAgQO54447uumdm3VR243QcmZvR2FWKE85XUa5KaefeuopPvjBD/ZSRM1lyy235PXXX++RY/n3at3KiUP96uht8O+nIk85bWZmZt3OiYPVnZ7qbTAzs65z4mBmZma5OXEwMzOz3Jw4mJmZWW5OHMzMzCw3Jw4NZNmyZUycOJFRo0ax66678rnPfY633nprg9o655xzWLx48Qbt2/FI7P32248PfvCDjBs3jpkzZ9bcb8GCBcyZM2eDjmlmZvXBE0BtqO6egazGvcURwcknn8y5557LPffcw9q1a5k6dSqXXHIJl19+eZcP9/3vf39DIwWy50zMnz8fgGeffZaTTz6Zd955hzPPrPw+FixYQFtbW81nbJiZWf1yj0ODeOihhxg0aNC6D+b+/ftz5ZVXcsMNN/Dmm29y0003MXHiRMaPH88HPvCBdU+OfOONNzjhhBPYd9992WuvvdbNwnjEEUfQMcnVlltuyUUXXcT+++/PMcccw7x58zjiiCPYZZddaG1trRnbLrvswhVXXMFVV10FwLx58zjkkEPYb7/9OOSQQ1iyZAlvvfUWl156KbNnz2bMmDHMnj2bN954g7POOosDDjiA/fbbj3vuuaeIU2dmZt3IPQ4NomMK5lJbb701O++8M+3t7UD2gf3kk0+y+eabc8ABB3DCCSfw/PPPs+OOO3LfffcBsHLlyvXafuONNzjiiCOYMWMGH/vYx/jKV77CAw88wOLFi5k8eTITJkyoGd/YsWP53e9+B8Aee+zBI488woABA3jwwQf58pe/zF133cX06dNpa2vj6quvBuDLX/4yRx11FDfccAOvvvoq48aN45hjjmGLLbbYqHNlZmbFceLQICKi7FMiS8s//OEPM2TIECB7rPUvfvELjj/+eC688EIuuugiTjzxRP7xH/9xvTY22WQTxo8fD2SP0950000ZOHAge++9d8VHYpeLo8PKlSuZPHkyTz/9NJJ4++23y+7z05/+lNbWVr75zW8CsGrVKl544QVPAW1mVsd8qaJB7LnnnnR+fsbf/vY3XnzxRXbddVdg/cdPS2L33XfnscceY++992batGlMnz59vbYHDhy4bt8NfbT2/Pnz133gf/WrX+XII4/kySef5Mc//jGrVq0qu09EcNddd7FgwQIWLFjgpMHMrAE4cWgQRx99NG+++SY333wzAGvXruULX/gCU6ZMYfPNNwfggQceYMWKFfz973/n7rvv5tBDD+Wll15i880355Of/CQXXnghjz/+eLfHtnTpUi688EI++9nPAlmPw7BhwwC46aab1tXbaquteO2119atH3vssXznO99Z11vRMdjSzKzp1fFjwJ04NAhJ/OhHP+KOO+5g1KhR7L777gwaNIh/+7d/W1fnsMMO41Of+hRjxozhlFNOoaWlhd/+9reMGzeOMWPGcNlll/GVr3ylW+J55pln1t2O+U//9E989rOfXTdw80tf+hLTpk3j0EMPZe3atev2OfLII1m8ePG6wZFf/epXefvtt9lnn33Ya6+9+OpXv9otsZlZL6rjDzzrHn6sdhmN+Fjtm2666T0DDy2fev+9WoPxY7Xr9/HV9RpXJb0Qrx+rbWZmZt3Od1X0EVOmTGHKlCm9HYaZmfVx7nEwMzOz3Nzj0AWV5lKwxuTxPWZNwAM1u517HHIaNGgQL7/8sj9s+oiI4OWXX2bQoEG9HYqZWUNxj0NOw4cPZ9myZSxfvry3Q7FuMmjQIIYPH97bYZiZNZRCEwdJ44F/B/oD34+Ib3TavilwM7A/8DLw8YhYmrZNA84G1gIXRMT91dqUdD3QAgj4PTAlIl6XNAW4HPhDOuzVEdHlR0MOHDiQkSNHdnU3MzOzPqWwSxWS+gPXAMcBo4HTJY3uVO1s4JWI2A24EpiR9h0NTAL2BMYD10rqX6PNz0fEvhGxD/ACcH7JcWZHxJj02rjnSZuZmTWxIsc4jAPaI+LZiHgLmAVM7FRnIjAzLd8JHK1s9OFEYFZErI6I54D21F7FNiPibwBp/80AD0YwMzPrZkUmDsOAF0vWl6WysnUiYg2wEhhSZd+qbUq6EfgTsAfwnZJ6p0haKOlOSTttxHsyMzNrakUmDuXuW+zcC1CpTlfLs4WIM4EdgaeAj6fiHwMj0iWMB3m3h+O9gUhTJbVJavMASDMzs/KKTByWAaXf7ocDL1WqI2kAMBhYUWXfmm1GxFpgNnBKWn85Ilanzf9JNhBzPRFxXUS0RETL0KFDc75FMzOz5lJk4vAoMErSSEmbkA12bO1UpxWYnJZPBR6KbKKEVmCSpE0ljQRGAfMqtanMbrBujMNJwO/S+g4lx5tA1hthZmZF85My+6TCbseMiDWSzgfuJ7t18oaIWCRpOtAWEa3A9cAtktrJehompX0XSbodWAysAc5LPQlUaLMfMFPS1mSXM54Azk2hXCBpQmpnBTClqPdsZmbW1xU6j0NEzAHmdCq7tGR5FXBahX0vAy7L2eY7wKEV2pkGTOtq7GZmvWpjHqvcaI+QtobiKafNzMwsNycOZmb2rr4wJsFjKwrlxMHMrC/zh6h1MycOZmZmlpsTBzMzM8vNiYOZmZnl5sTBzMzMcnPiYGZmZrk5cTAzM7PcnDiYmZlZbk4czMzMLLdCn1VhZmY9qLsmevKzLqwK9ziYmZlZbk4czMzMLDcnDmZmZpabEwczMyvWxjxoq6ce0OWHgeXmxMHMzMxyc+JgZmZmuTlxMDMzs9ycOJiZmVluThzMzMwsNycOZmZmlpsTBzMzM8vNiYOZ2Ybyvf/WhJw4mJnZxnEC1VScOJiZmVluThzMzMzqVR325DhxMDMzs9ycOJiZmVluThzMzMyqqcPLBb3JiYOZWSPrC3c09IX30BPq5DwVmjhIGi9piaR2SReX2b6ppNlp+1xJI0q2TUvlSyQdW6tNSddLekLSQkl3Stqy1jHMzMysawpLHCT1B64BjgNGA6dLGt2p2tnAKxGxG3AlMCPtOxqYBOwJjAeuldS/Rpufj4h9I2If4AXg/GrHMDMzs64rssdhHNAeEc9GxFvALGBipzoTgZlp+U7gaElK5bMiYnVEPAe0p/YqthkRfwNI+28GRI1jmJlZX1EHXfjNosjEYRjwYsn6slRWtk5ErAFWAkOq7Fu1TUk3An8C9gC+U+MYZmZm1kVFJg7lvtVHzjpdLc8WIs4EdgSeAj7ehTiQNFVSm6S25cuXl9nFzMzqgnsXelWRicMyYKeS9eHAS5XqSBoADAZWVNm3ZpsRsRaYDZxS4xh02u+6iGiJiJahQ4fmfpNmZtZH1cldDPWmyMThUWCUpJGSNiEb7NjaqU4rMDktnwo8FBGRyielOyJGAqOAeZXaVGY3WDfG4STgdzWOYWZmZl00oKiGI2KNpPOB+4H+wA0RsUjSdKAtIlqB64FbJLWT9QJMSvsuknQ7sBhYA5yXehKo0GY/YKakrckuTTwBnJtCKXsMMzMz67rCEgeAiJgDzOlUdmnJ8irgtAr7XgZclrPNd4BDK7RT8RhmZmbr8eWJqjxzpJmZWT1okISl0B4HMzNrEA3yoWW9zz0OZmbWN1VKhurhbol6iGEDOXEwM2sEHR80Dfph03T68O/JiYOZmZnl1qXEQVK/dMujmZmZNaGaiYOkWyVtLWkLsnkVlkj6YvGhmZmZNZkGuByVp8dhdHry5EfJ5k/YGfhUoVGZmZlZXcqTOAyUNJAscbgnIt6mzEOizMzMrO/LM4/DfwBLyaZxfkTS+4G/FRmUmVlDqfOu5YbRcR5bzuzdOKyqmolDRFwFXFVS9LykI4sLyczM+iwnBw0vz+DI/ybpekk/Seujefdpk2Zm1te5R8VK5BnjcBPZ0yh3TOu/B/5HUQGZmdWlBhjtbn1Unf3by5M4bB8RtwPvQPa4bGBtoVGZmZlZXcqTOLwhaQjpTgpJBwErC43KzMwaWx19Q7buleeuiv8JtAK7SvolMBQ4rdCozMysvnhQoyV5EodFwOHABwABS/AzLszMGlPbjf7wt42SJwH4dUSsiYhFEfFkmgDq10UHZmbWFMp16bubv/HU2QDGIlXscZD0PmAYsJmk/ch6GwC2BjbvgdjMzKyeNckHpb1XtUsVxwJTgOHAFSXlrwFfLjAmM7PG1xfGBDgxsDIqJg4RMROYKemUiLirB2MyM7NqPE7BelGeKafvknQCsCcwqKR8epGBmZmZWf3JM+X094CPA58lG+dwGvD+guMyMzPrW/rIpZ88d1UcEhGfBl6JiK8DBwM7FRuWmZmZ1aM8icPf0883Je0IvA2MLC4kM7MG0Be+PTbRLYTWffIkDvdK2ga4HHgcWArMKjIoMzPrBo2QFHgei4aTZ3Dkv6bFuyTdCwyKCD+rwszMGk+tpKQ7k5a8t+Q2WKJUbQKok6tsIyJ+WExIZma2ToN9qDQcn98uq9bjcFL6+Q/AIcBDaf1I4GHAiYOZmVlP6+V5PKpNAHUmQLo8MToi/pjWdwCu6ZnwzMz6KH/T7Tqfs7qQZ3DkiI6kIfkzsHtB8ZiZmVkdy5M4PCzpfklTJE0G7gN+VnBcZmbWU/xNvuua+JzVTBwi4nzgP4B9gTHAdRHx2TyNSxovaYmkdkkXl9m+qaTZaftcSSNKtk1L5UskHVurTUk/SOVPSrpB0sBUfoSklZIWpNeleWI3MzOz9VW7q+J+4P8AP0l3UHRpMKSk/mRjIT4MLAMeldQaEYtLqp1NNiPlbpImATOAj0saDUwiez7GjsCDkjouj1Rq8wfAJ1OdW4FzgO+m9Z9HxIldid/MzKxwDdhzUa3HYTLwCvAvkh6X9F1JEyVtmbPtcUB7RDwbEW+RTRo1sVOdicDMtHwncLQkpfJZEbE6Ip4D2lN7FduMiDmRAPPIHgduZmZm3ahi4hARf4qImyJiEtAC3AzsD9wv6UFJX6rR9jDgxZL1ZamsbJ2IWAOsBIZU2bdmm+kSxafIeks6HCzpCUk/kbRnjbjNzMysgpozRwJExDvAr9PrUknbA8dW3wuVaypnnUrl5RKdzm1eCzwSET9P648D74+I1yUdD9wNjFovWGkqMBVg5513LnMYMzMzq5k4SBoK/HdgRGn9iDirxq7LeO9TNIcDL1Wos0zSAGAwsKLGvhXblPQ1YCjwmZI4/1ayPEfStZK2j4i/lgYSEdcB1wG0tLR0TkbMzAzyT6Oct53e1F0x1MN76UF5ehzuAX4OPAis7ULbjwKjJI0E/kA22PETneq0ko2l+DVwKvBQRISkVuBWSVeQDY4cRTZuQZXalHQOWS/I0amHhFT+PuDPqd1xZL0WL3fhfZiZFau7PoyteE2WJJSTJ3HYPCIu6mrDEbFG0vnA/UB/4IaIWCRpOtAWEa3A9cAtktrJehompX0XSbodWAysAc6LiLUA5dpMh/we8Dzw62x8JT+MiOlkCcm5ktaQPSJ8UhpAaWZmG8ofoE0rT+Jwr6TjI2JOVxtP+8zpVHZpyfIq4LQK+14GXJanzVRe9r1ExNXA1V0K3MzMzMrKkzh8DviypLeAt1NZRMTWxYVlZmZV+Ru/9ZKaiUNEbNUTgZiZNQV/4HcPjwvpNblux5Q0AfhQWn04Iu4tLiQzMyuUkxfbCDWfVSHpG2SXKxan1+dSmZmZmTWZPD0OxwNjOm5xlDQTmA+s99AqMzMz69tyXaoAtiG7XRKySZrMzKyR+PKEdZM8icP/C8yX9DOyCZg+BEwrNCozM+t7nLz0CXnuqrhN0sPAAWSJw0UR8aeiAzMzsybmJKNuVRwcKWmP9HMssAPZ8yNeBHZMZWZmZtZkqvU4/E+yp0V+q8y2AI4qJCIzM7O+pI/1nlRMHCJialo8Lk0NvY6kQYVGZWZmZnWp5jwOwK9ylpmZmVkfV7HHIT2OehiwmaT9yAZGAmwNbN4DsZmZ1Z8+1u3c8Pz76HHVxjgcC0wBhgNXlJS/Bny5wJjMzMy6rhGSiEaIsYZqYxxmAjMlnRIRd/VgTGZmZlan8szjcJekE4A9gUEl5dOLDMzMzMzqT83EQdL3yMY0HAl8HzgVmFdwXGZmfUcf6J4265DnropDIuLTwCsR8XXgYGCnYsMyMzOzepQncfh7+vmmpB2Bt4GRxYVkZmZm9SrPQ67ulbQNcDnwONmskd8vNCozMzOrS3kGR/5rWrxL0r3AoIhYWWxYZmbWJR5HUSyf33WqTQB1cpVtRMQPiwnJzMzM6lW1HoeTqmwLwImDmZm9y9/Km0K1CaDO7MlAzMzMrP7lmcfh0nLlngDKzKyb+Ru7NYA8d1W8UbI8CDgReKqYcMzMzKye5bmr4lul65K+CbQWFpGZmZnVrTwTQHW2ObBLdwdiZmZm9S/PGIffkt1FAdAfGAp4fIOZmVkTyjPG4cSS5TXAnyNiTUHxmNWVW+e+AMAnDty5lyMxM6sPecY4PC9pLHAYWc/DL4D5RQdmZtaQfGeE9XE1xzik2zFnAkOA7YGbJH0lT+OSxktaIqld0sVltm8qaXbaPlfSiJJt01L5EknH1mpT0g9S+ZOSbpA0MJVL0lWp/sKUBJmZmXW/thv7fPKYZ3Dk6cABEfG1iPgacBBwRq2dJPUHrgGOA0YDp0sa3ana2WSP694NuBKYkfYdDUwC9gTGA9dK6l+jzR8AewB7A5sB56Ty44BR6TUV+G6O92xm3azjso+ZNbY8icNSsvkbOmwKPJNjv3FAe0Q8GxFvAbOAiZ3qTCTrzQC4EzhaklL5rIhYHRHPAe2pvYptRsScSIB5wPCSY9ycNv0G2EbSDjniNzPrVnOfW8Hc51b0dhhmGyVP4rAaWCTpJkk3Ak8Cr6fu/6uq7DcMeLFkfVkqK1snDbhcSXZJpNK+NdtMlyg+BfyfLsRhZmZmOeS5q+JH6dXh4Zxtq0xZ5KxTqbxcotO5zWuBRyLi512IA0lTyS5lsPPOHkFvZs1p7nMrOHDkdr0dhtWxPHdVzJS0CbB7KloSEW/naHsZsFPJ+nDgpQp1lkkaAAwGVtTYt2Kbkr5GNs/EZ7oYBxFxHXAdQEtLy3qJhZmZmeW7q+II4GmyQYnXAr+X9KEcbT8KjJI0MiUek1h/qupWYHJaPhV4KI1RaAUmpbsuRpINbJxXrU1J5wDHAqdHxDudjvHpdHfFQcDKiPhjjvjN+oxb577gwYlm1i3yXKr4FvCRiFgCIGl34DZg/2o7RcQaSecD95PNOHlDRCySNB1oi4hW4HrgFkntZD0Nk9K+iyTdDiwmm3TqvIhYm46/XpvpkN8Dngd+nY2v5IfpCZ5zgOPJBli+Cfhx4WbWc/r4rXnWfPIkDgM7kgaAiPh9xxwJtUTEHLIP7tKyS0uWVwGnVdj3MuCyPG2m8rLvJfVgnJcnXjMzM6suT+LQJul64Ja0fgbwWHEhmZmZWb3KkzicS/aN/QKyOxQeIRvrYGZmZk0mz10VqyVdDfwX8A7ZXRVvFR6ZmZmZ1Z08d1WcQDZT5L8DVwPtko4rOjAzay6+68OsMeS9q+LIiGgHkLQrcB/wkyIDMzMzs/qTZ8rpv3QkDcmzwF8KiseSvnLffV94D0XrK79rM2sOeXocFkmaA9xONlXzacCjkk4GiIgfFhifmZmZ1ZE8PQ6DgD8DhwNHAMuB7YCTgBMLi8zMrIuK7Lnxky3NMnnuqvBMi2ZmZgbku1RhZmZmnTXpdOJ5LlWYmZmZAVUSB0mfSz8P7blwzKwe9YW7Pnz3iln3qNbj0DG24Ts9EYiZmZnVv2qJw1OSlgIfkLSw5PVbSQt7KD6rwN+czMysN1QcHBkRp0t6H3A/MKHnQrJmcuvcF/jEgTv3dhhmuc19bgUHjtyut8Mw6zVVB0dGxJ8iYl/gj8BW6fVSRDzfE8GZ1Tv3/Fij8VwUtrFq3o4p6XDgZmAp2WO1d5I0OSIeKTg2K5i/7fd9PZHYdBzD/5bMmkOeeRyuAD4SEUsAJO0O3AbsX2RgZmZmVn/yzOMwsCNpAIiI3wMDiwvJzLrCl0vMrCflSRzaJF0v6Yj0+k/gsaIDs43TFz5M6um++3qJo9n1xu+hnv4dmtWDPJcqzgXOAy4gG+PwCHBtkUGZNRJf4zezZpLnIVerycY5XFF8OGZmZlbP/KwKq6hRu2fzdi03yvtrlDjNrDn46Zi2Qdw9n5/PlXWnjnkYemMSqs5zQHgirObkHgczMzPLLc8EULsDXwTeX1o/Io4qMC4zMzOrQ3kuVdwBfA/4T2BtseGYmTWv3rwMYZZXnsRhTUR8t/BIzKxP8uBOs74lzxiHH0v6fyTtIGm7jlfhkVmvKb0roREmv6n3+BqBz6GZ5ZWnx2Fy+vnFkrIAdun+cKzR+I6BjeMPbOtuRV3u8FM1rUPNHoeIGFnm5aTB+hR/gNevjfndNEKPmVmjyXNXxUCyaac/lIoeBv4jIt4uMC7r4/zH3MysMeUZ4/BdskdoX5te+6eymiSNl7REUruki8ts31TS7LR9rqQRJdumpfIlko6t1aak81NZSNq+pPwISSslLUivS/PEbmZmZuvLM8bhgIjYt2T9IUlP1NpJUn/gGuDDwDLgUUmtEbG4pNrZwCsRsZukScAM4OOSRgOTgD2BHYEH03wSVGnzl8C9ZD0inf08Ik7M8V7rzq1zX/D4gQZWqWelt3qzunfmAAASoElEQVRcOv97KjqOeulZqoc4PEbA+oo8PQ5rJe3asSJpF/LN5zAOaI+IZyPiLWAWMLFTnYnAzLR8J3C0JKXyWRGxOiKeA9pTexXbjIj5EbE0R1x9Rm9dv62HP8LWXHry3/qGHseJgTWLPD0OXwR+JulZssdqvx84M8d+w4AXS9aXAQdWqhMRayStBIak8t902ndYWq7VZjkHp16Sl4ALI2JR5wqSpgJTAXbe2d/we1t3PKSqdHsz9to4wes7emtiKCdDVk6ex2r/l6RRwAfIEoffpUdt16JyzeWsU6m8XA9J5zY7exx4f0S8Lul44G5g1HqNRFwHXAfQ0tJSq82G1uwfKD3x/rvjGM3+e8rL58msZ1VMHCQdFREPSTq506ZdJRERP6zR9jJgp5L14WTf+MvVWSZpADAYWFFj31ptvkdE/K1keY6kayVtHxF/rRG/WZ9W7QPX83OYWSXVxjgcnn6eVOaVZ6Dho8AoSSMlbUI22LG1U51W3p1g6lTgoYiIVD4p3XUxkqyHYF7ONt9D0vvSuAkkjUvv+eUc8RdqQ67Z+puVFWljxxFs6L6NPFanXrry5z63om5isb6vYo9DRHwtLU5PAxTXSR/mVaUxC+cD9wP9gRsiYpGk6UBbRLQC1wO3SGon62mYlPZdJOl2YDGwBjgvItamY6/XZiq/APgS8D5goaQ5EXEOWUJyrqQ1wN+BSSk5sQI1y2j9ZuJeCDODfIMj7wLGdiq7k2w+h6oiYg4wp1PZpSXLq4DTKux7GXBZnjZT+VXAVWXKrwaurhWrmZmZ1VZtjMMeZPMoDO40zmFrYFDRgdnGy/MNsa99i/Tln43nc7Lh6vWx2HOfW1F3MVnjqtbj8AGysQzbkI1r6PAa8N+LDMqsKM3wodgM77G7NetEa/Wa6Fh9qzbG4R7gHkkHR8SvezCmPqs7vt3X84fChr6/eu31yHPXQXcfr97OQT3r7t9BPf/fyssDJK0n5Bnj8M+SnoqIVwEkbQt8KyLOKjY0606NkrT0hT/e3cHnoXFsyIf1hl46KDIxcNJheeWZcnqfjqQBICJeAfYrLiSz9/Kjkc3W150f9E4arCvyJA79Ui8DAJK2I19PheXgb/FmxaiXf/e15ljwh7Y1mjwJwLeAX0m6M62fRpnbJM3MzKzvy/OsipslPQYcSfYMiZM7PRrbuqhevgmVszGz/1lx6nUA6cbwv5n65V4QqybXJYc0k+Ny0vwNknaOCP+vL0Bf/mPaFz/8usJ3TXSPvvx/xKwR1EwcJE0gu1yxI/AXssdqP0U2OZRZj/KYEOuL/A3fGkmeHod/BQ4CHoyI/SQdCZxebFhWlEb+UNyYeRXyvO9GPje2Yfw7N+u6PInD2xHxsqR+kvpFxM8kzSg8Musz/Me5b6n3Sy7+9/Yu92RYEfIkDq9K2hJ4BPiBpL+QPbHSrKxm/8Nd5Pvv6XEizf677C1+toTVszyJw0Syx1F/HjgDGAxMLzIo23D+Q2+NrNkH0Jo1gqqJg6T+wD0RcQzwDjCzR6IyMzOzulQ1cYiItZLelDQ4Ilb2VFBmjagZe3ua8T33FD+50upVnksVq4DfSnoAeKOjMCIuKCwq6zL/AW9u9dLFXy9xVOP/K93HYzGaU57E4b70MrM6Ue93Nlj98R0W1l0qJg4ds0NGhMc1WLfyN77i9NS53dipyZ30FM+JghWlWo/D3cBYAEl3RcQpPROSdeYP2vrgDz0zs+qJg0qWdyk6kL6sET74GyHGZuDfQ6a3L8X492BWWbXEISosW0H8x6ox1MvvqV7iMLPmUi1x2FfS38h6HjZLy6T1iIitC4/OzMzM6krFxCEi+vdkIGZmtj4PcrR6k+d2TDMz2wD1+qFfr3FV1HZjb0dgJZw4mJn1gob78DZLnDjUAQ9yy/g89IyeOM/N/rt0UmB9Wb/eDsA2XrP/kbbG43+zZo3LPQ5mVpecXJjVJ/c4mJmZWW5OHCwXf/szMzMoOHGQNF7SEkntki4us31TSbPT9rmSRpRsm5bKl0g6tlabks5PZSFp+5JySboqbVsoaWxx79isulvnvuAkzLrN3OdWeCCm9bjCEgdJ/YFrgOOA0cDpkkZ3qnY28EpE7AZcCcxI+44GJgF7AuOBayX1r9HmL4FjgOc7HeM4YFR6TQW+253v08zMrJkUOThyHNAeEc8CSJoFTAQWl9SZCPxLWr4TuFqSUvmsiFgNPCepPbVHpTYjYn4q6xzHRODmiAjgN5K2kbRDRPyxW9+tmTUNf8u3ZlbkpYphwIsl68tSWdk6EbEGWAkMqbJvnjY3JA4zMzPLocjEYb2v/qz/lM1KdbpavrFxIGmqpDZJbcuXL6/RpJmZWXMqMnFYBuxUsj4ceKlSHUkDgMHAiir75mlzQ+IgIq6LiJaIaBk6dGiNJs3MzJpTkYnDo8AoSSMlbUI22LG1U51WYHJaPhV4KI1FaAUmpbsuRpINbJyXs83OWoFPp7srDgJWenyDmZnZhilscGRErJF0PnA/0B+4ISIWSZoOtEVEK3A9cEsa/LiCLBEg1budbCDlGuC8iFgL2W2XndtM5RcAXwLeByyUNCcizgHmAMcD7cCbwJlFvWczM7O+TtkXfCvV0tISbW1t3dae79s3a3y7vnBHb4dQlw4cuV1vh9CcWrr/O7CkxyKipVY9zxxpZmZmuTlxMDMzs9ycOJiZmVluThzMzMwsNycOZmZmlpsTBzMzM8vNiYOZmZnl5sTBzMzMcnPiYGZmZrk5cTAzM7PcnDiYmZlZbk4czMzMLDcnDmZmZpabEwczMzPLzYmDmZmZ5ebEwczMzHJz4mBmZma5OXEwMzOz3Jw4mJmZWW5OHMzMzCw3Jw5mZmaWmxMHMzMzy82Jg5mZmeXmxMHMzMxyc+JgZmZmuTlxMDMzs9ycOJiZmVluThzMzMwsNycOZmZmlpsTBzMzM8vNiYOZmZnlVmjiIGm8pCWS2iVdXGb7ppJmp+1zJY0o2TYtlS+RdGytNiWNTG08ndrcJJVPkbRc0oL0OqfI92xmZtaXFZY4SOoPXAMcB4wGTpc0ulO1s4FXImI34EpgRtp3NDAJ2BMYD1wrqX+NNmcAV0bEKOCV1HaH2RExJr2+X8DbNTMzawpF9jiMA9oj4tmIeAuYBUzsVGciMDMt3wkcLUmpfFZErI6I54D21F7ZNtM+R6U2SG1+tMD3ZmZm1pSKTByGAS+WrC9LZWXrRMQaYCUwpMq+lcqHAK+mNsod6xRJCyXdKWmnjXlTZmZmzazIxEFlyiJnne4qB/gxMCIi9gEe5N0ejvcGIk2V1Capbfny5eWqmJmZNb0iE4dlQOm3++HAS5XqSBoADAZWVNm3UvlfgW1SG+85VkS8HBGrU/l/AvuXCzYirouIlohoGTp0aBfeppmZWfMoMnF4FBiV7nbYhGywY2unOq3A5LR8KvBQREQqn5TuuhgJjALmVWoz7fOz1AapzXsAJO1QcrwJwFPd/D7NzMyaxoDaVTZMRKyRdD5wP9AfuCEiFkmaDrRFRCtwPXCLpHaynoZJad9Fkm4HFgNrgPMiYi1AuTbTIS8CZkn6X8D81DbABZImpHZWAFOKes9mZmZ9nbIv61aqpaUl2trauq29W+e+0G1tmVnv2PWFO3o7hLp04MjtejuE5tRyZrc3KemxiGipVc8zR5qZmVluThzMzMwsNycOZmZmlpsTBzMzM8vNiYOZmZnl5sTBzMzMcnPiYGZmZrk5cTAzM7PcnDiYmZlZbk4czMzMLDcnDmZmZpabEwczMzPLzYmDmZmZ5ebEwczMzHJz4mBmZma5OXEwMzOz3Jw4mJmZWW5OHMzMzCw3Jw5mZmaWmxMHMzMzy82Jg5mZmeXmxMHMzMxyc+JgZmZmuTlxMDMzs9ycOJiZmVluThzMzMwsNycOZmZmlpsTBzMzM8vNiYOZmZnl5sTBzMzMcnPiYGZmZrk5cTAzM7PcCk0cJI2XtERSu6SLy2zfVNLstH2upBEl26al8iWSjq3VpqSRqY2nU5ub1DqGmZmZdU1hiYOk/sA1wHHAaOB0SaM7VTsbeCUidgOuBGakfUcDk4A9gfHAtZL612hzBnBlRIwCXkltVzyGmZmZdV2RPQ7jgPaIeDYi3gJmARM71ZkIzEzLdwJHS1IqnxURqyPiOaA9tVe2zbTPUakNUpsfrXEMMzMz66IiE4dhwIsl68tSWdk6EbEGWAkMqbJvpfIhwKupjc7HqnQMMzMz66IBBbZd7lt95KxTqbxcolOtft44kDQVmJpWX5e0pMx+1n22B/7a20E0AZ/nnuHz3DN8ntc5q4hG35+nUpGJwzJgp5L14cBLFeoskzQAGAysqLFvufK/AttIGpB6FUrrVzrGe0TEdcB1XXyPtoEktUVES2/H0df5PPcMn+ee4fNcH4q8VPEoMCrd7bAJ2WDH1k51WoHJaflU4KGIiFQ+Kd0RMRIYBcyr1Gba52epDVKb99Q4hpmZmXVRYT0OEbFG0vnA/UB/4IaIWCRpOtAWEa3A9cAtktrJegEmpX0XSbodWAysAc6LiLUA5dpMh7wImCXpfwHzU9tUOoaZmZl1nfzl23qDpKnp8pAVyOe5Z/g89wyf5/rgxMHMzMxy85TTZmZmlpsTB9sokpZK+q2kBZLaSso/m6YGXyTpf5eUd8tU4s2m3HmWNEbSbzrKJI1L5ZJ0VTqXCyWNLWlncjqXT0uaXFK+f2q/Pe3blJOkSdpG0p2SfifpKUkHS9pO0gPpnD0gadtU1+d5A1U4z5en9YWSfiRpm5L6/rtRTyLCL782+AUsBbbvVHYk8CCwaVr/h/RzNPAEsCkwEniGbJBr/7S8C7BJqjM67XM7MCktfw84t7ffcx2d558Cx6Xl44GHS5Z/QjaHyUHA3FS+HfBs+rltWt42bZsHHJz2+UlHu832Iptl9py0vAmwDfC/gYtT2cXADJ/nQs7zR4ABqWxGyXn23406e7nHwYpwLvCNiFgNEBF/SeXdOZW4ZROZbZ2WB/Pu3CUTgZsj8xuyOU52AI4FHoiIFRHxCvAAMD5t2zoifh3ZX9qbacLzLGlr4EOkO7Ii4q2IeJX3TlvfeTp7n+cuqnSeI+Kn8e7sv78hm48H/Hej7jhxsI0VwE8lPaZs9k2A3YF/TF2F/7+kA1J5d04l3mzKnef/AVwu6UXgm8C0VN7V8zwsLXcubza7AMuBGyXNl/R9SVsA/y0i/giQfv5Dqu/zvGEqnedSZ5H1yID/btQdJw62sQ6NiLFkTyw9T9KHyOYH2Zas+/aLwO3pW0BXpxjPNV14kyh3ns8FPh8ROwGf5925S3yeN8wAYCzw3YjYD3iD7NJEJT7PG6bqeZZ0Cdn8PT/oKCrThs9zL3LiYBslIl5KP/8C/Iis+3AZ8MPUhTsPeIdsjvlKU4lXKl83lXin8qZT4TxPBn6YqtyRyqDr53kZ73YLl5Y3m2XAsoiYm9bvJPuA+3O6zED6+ZeS+j7PXVfpPJMGkp4InJEu53TU99+NOuLEwTaYpC0kbdWxTDa46UngbrJrjEjanWzg0l/p3qnEm0aV8/wScHiqdhTwdFpuBT6dRv0fBKxMXez3Ax+RtG26M+AjwP1p22uSDko9Q5+mCc9zRPwJeFHSB1LR0WSz15ZOW995Onuf5y6qdJ4ljSebAXhCRLxZsov/btSb3h6d6VfjvsiuVT6RXouAS1L5JsD/R/bh9jhwVMk+l5CNhF5CyYhyshHqv0/bLul0jHlkA6LuIN2p0UyvKuf5MOCxVD4X2D+VC7gmncvfAi0lbZ2VzmU7cGZJeUv6fT0DXE2aHK7ZXsAYoA1YSJYAb0t2zfy/yBKz/wK283ku5Dy3k41ZWJBe3yup778bdfTyzJFmZmaWmy9VmJmZWW5OHMzMzCw3Jw5mZmaWmxMHMzMzy82Jg5mZmeXmxMHMeoykkHRLyfoAScsl3ZvWp0i6Oi3/i6Q/KHv659OSfihpdG/FbmYZJw5m1pPeAPaStFla/zDwhyr1r4yIMRExCpgNPCRpaNFBmlllThzMrKf9BDghLZ8O3JZnp4iYTfYo8U8UFJeZ5eDEwcx62iyyKYQHAfuQzXqZ1+PAHoVEZWa5OHEwsx4VEQuBEWS9DXO6uHu5Jx+aWQ9y4mBmvaEV+CY5L1OU2A94qvvDMbO8BtSuYmbW7W4ge5rkbyUdkWcHSaeQPWnyC0UGZmbVOXEwsx4XEcuAfy+zaQCwumT985I+CWxB9lTJoyJieQ+EaGYV+OmYZlY3JF0JPB0R1/Z2LGZWnhMHM6sLkn4CbAKcHBErezseMyvPiYOZmZnl5rsqzMzMLDcnDmZmZpabEwczMzPLzYmDmZmZ5ebEwczMzHJz4mBmZma5/V++TCkhyS+gmAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "n,b,p = plt.hist(mini.orbits.epoch, 200, density=True, alpha=0.4, label='Minimoon Epoch')\n",
    "n,b,p = plt.hist(simData['observationStartMJD'], 100, density=True, alpha=0.4, label='Opsim Date')\n",
    "plt.legend(loc='upper left', fancybox=True, numpoints=1)\n",
    "plt.xlabel(\"MJD\")\n",
    "plt.ylabel('Fraction of population/dates')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "59853.016794 61405.5500411 4.25351574545\n"
     ]
    }
   ],
   "source": [
    "print(simData['observationStartMJD'].min(), mini.orbits.epoch.max(), \n",
    "      (mini.orbits.epoch.max() - simData['observationStartMJD'].min())/365)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5232"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = simData['observationStartMJD'].min() - 260\n",
    "late = mini.orbits.query('epoch > @x')\n",
    "len(late)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "35.3 ms ± 130 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "ephs = pyo.generateEphemerides(times, timeScale='TAI', obscode='I11', ephMode='2body')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "ephs = pyo.generateEphemerides(times, timeScale='TAI', obscode='I11', ephMode='2body')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.7325801556433014"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.diff(ephs['ra']*np.cos(np.radians(ephs['dec']))).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.8829690540305677"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.diff(ephs['dec']).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Looks like objects can travel about 8 degrees per day (is there higher variation on shorter timescales?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "raTol = 10\n",
    "decTol = 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "eph_idxs = np.searchsorted(ephs[0]['time'], simData['observationStartMJD'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(623453, 623453, 1094)"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(eph_idxs), len(simData), len(ephs[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "59853.56 59853.016794\n",
      "59853.56 59853.0206366\n",
      "59853.56 59853.021088\n",
      "59853.56 59853.0215278\n",
      "59853.56 59853.0219792\n",
      "60946.56 60946.3950347\n"
     ]
    }
   ],
   "source": [
    "for i in (0, 1, 2, 3, 4, -1):\n",
    "    print(ephs[0]['time'][eph_idxs[i]], simData['observationStartMJD'][i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "dts = ephs[0]['time'][eph_idxs] - simData['observationStartMJD']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xb6b0c3a20>]"
      ]
     },
     "execution_count": 79,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD8CAYAAACcjGjIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X9w3Hd95/HnW6sfdmQZWYqMTGys5MZJ7btrYqzK7NK4FG6kQGdiOFKGVjpcfpwE9NcNN5XcaXud3s2EkVuGg+HX7oXJJSMxQJ00ca+AlHKUXLqLhHKYQOwkdoKNRWT8I3ZjWziypff9sd/VbSzZWsva/e5qX4+Zz+xXH333u++vJO/L3+/38/2suTsiIiLZKsIuQEREio/CQURE5lA4iIjIHAoHERGZQ+EgIiJzKBxERGSOnMLBzOrNbK+ZPWdmB80sama/bWbPmtmMmbVesf6fmtlhM3vezDqy+u8J+g6b2e6s/lvNbMTMDpnZ182seul2UURErleuRw6fBb7t7r8C3AkcBH4C/HvgyewVzWwL8AHgXwP3AF80s4iZRYAvAO8CtgC/E6wL0A98xt03AWeAj9zQXomIyA1ZMBzMbDWwA/gKgLtPuftZdz/o7s/P85SdwNfc/TV3/ylwGGgL2mF3f8ndp4CvATvNzIB3AHuD5z8EvOdGd0xERBavMod1bgNOAg+a2Z3A08Afu/uFq6x/C/D9rK/Hgz6AY1f0bwcagbPufnme9a/q5ptv9paWlhzKFxGRjKeffvqUuzcttF4u4VAJvAX4Q3cfMbPPAruBv7jK+jZPnzP/UYpfY/25GzbrBroB3vzmNzM2NrZA6SIiks3MjuayXi7XHMaBcXcfCb7eSzosrrX+hqyv1wMvX6P/FFBvZpVX9M/h7gl3b3X31qamBYNPREQWacFwcPfjwDEzuyPoeidw4BpP2Qd8wMxqzOxWYBMwCvwA2BSMTKomfdF6n6dn/vsucF/w/F3A44vaGxERWRK5jlb6Q2DQzJ4B7gLuN7P3mtk4EAX+wcyGANz9WeAbpAPk28Dvu/t0cE3hD4Ah0qOdvhGsC9AHfNLMDpO+BvGVpdk9ERFZDCvVKbtbW1td1xxERK6PmT3t7q0Lrac7pEVEZA6Fg4iIzKFwKJBUKsXWrVuprKzEzGZbJBKhsrKSSCQyu1xZWUlVVRV1dXXcfvvtNDY2zq6TvW5VVRVVVVWzz8n0r1ixgpaWFhKJBACJRII3velNVFRUvO61r2yRSIRVq1bR19eX95/Fxz/+cbZu3cpNN91EJBKhurqarq4uALq6uli5ciWRSGROjQvtQ+ZnMt9z59vfysrKBbd55evP97vK/B6yf2+33norjY2NrFmzhg0bNlBXV0dNTQ01NTWzz1mxYsXsY2NjI1u2bKGjo4O6urrrqmsp25V/o2o31tauXUsqlcrrv6m8cPeSbNu2bfNSEI/Hvaqqyknfu6Gmplamrb29Pey3I3d3B8Y8h/dYHTkssb6+PlavXj37P9eenh4uXboUdlkiErLh4eHZo+NSoHBYAl1dXbOnAPbs2cO5c+eYmZkJuywRKTKDg4O0lMi0PwqHG9TV1cXg4CBeokOCRaSwjh49SmNjY9hlLCiXuZXkCqlUiocffhiARx55JORqRKTUvPLKK3R0dDA0NBR2KVelcMhBX18fX/3qV3F3Tp48ydTUVNgliUiJGx4exszYvHkzBw5ca0aicCgcriGVSnHvvfdy6tSpsEsRkWXq4MGDtLS0cOTIkbBLeR1dc5hHV1cXkUiEWCymYBCRvDt69Cg1NTVFdT+EwiGQSqVYu3YtZsbg4KBGG4lIQU1NTRGLxYpmuKvCgfQdxLFYjJMnT4ZdioiUucHBwbzPUpCLsg6Hvr6+2RvVRESKxZ49e0IPiLKdsrulpYWjR3P6tDwRkVBs3LhxyS9Um6bsvrqOjg4Fg4gUvaNHj9LR0RHKa5dlOAwPD4ddgohIToaHh6mtrS3465ZdOITxQxYRuRGTk5NUVhb2trSyCoft27czOTkZdhkiItdtenq6oKeYyiocRkdHwy5BRGTRhoeH2b59e0Feq6zCQUSk1I2OjrJly5a8v47CQUSkxBw8eDDv90GUVTg0NDSEXYKIyJL49Kc/ndftl1U4fOpTnwq7BBGRJTE9PZ3X7ZdVOHR3d9Pc3Bx2GSIiRa+swgFgYmKCioqy220RWWZqamryuv2yfJecnp7O+w9WRCSf8v2JlDmFg5nVm9leM3vOzA6aWdTMGszsCTM7FDyuCdZ9g5n9vZn9yMyeNbMPZW1nV7D+ITPbldW/zcx+bGaHzexzZmZLv6uvd/HiRQWEiJSslStX5nX7uR45fBb4trv/CnAncBDYDXzH3TcB3wm+Bvh94IC73wm8Hfi0mVWbWQPwl8B2oA34y0ygAF8CuoFNQbvnRncsFxcvXqSurq4QLyUisqQ+85nP5HX7C4aDma0GdgBfAXD3KXc/C+wEHgpWewh4T7DsQF3wv/9VwCvAZaADeMLdX3H3M8ATwD1mtg5Y7e4pT88f/nDWtvLu1VdfZePGjYV6ORGRG9bZ2Ul3d3deXyOXI4fbgJPAg2b2QzN7wMxqgTe6+wRA8Lg2WP/zwGbgZeDHwB+7+wxwC3Asa7vjQd8twfKV/QVz5MgRjWISkZLQ3t7OwMBA3l8nl3CoBN4CfMndtwIX+P+nkObTAewH3gTcBXw+OPqY7zqCX6N/DjPrNrMxMxtb6o/0nJiYYPPmzUu6TRGRpdTb28vQ0FBBXiuXcBgHxt19JPh6L+mw+EVwSojg8UTw/Q8Bj3raYeCnwK8E29mQtd31pI8uxoPlK/vncPeEu7e6e2tTU1Mu+3ddDhw4QGdn55JvV0TkRrg77k5/f3/BXnPBcHD348AxM7sj6HoncADYB2RGHO0CHg+Wfxasg5m9EbgDeAkYAtrNbE1wIbodGApOSZ0zs7cG1yk+mLWtghsYGMDddRQhIqGrqqoitI9yziTStRrp00NjwDPAY8AaoJH0KKVDwWNDsO6bgGHS1xt+AnRlbefDwOGgfSirvzVY90XS1yxsoZq2bdvm+dbe3u6kT3GpqampFbS1t7fn5X0NGPNc3vdzWakYWyHCwd09mUyG/keipqZWHq2mpsY7Ozvz+p5GjuFQ2M+dK0HRaJRkMkksFgu7FBFZppqbm5mYmAi7jNcpy+kzrlc0GtV1CBFZcrW1tSSTyaILBlA4XJcDBw6QTCbDLkNESpSZEY/HZ0/dnD9/nmg0GnZZ81I4XKfMaaZVq1aFXYqIlIiamhri8TgzMzN5v7N5qSgcFiEajXLu3Dl6e3vDLkVEipSZ0dnZibtz8eLFkgmFDIXDDejv78fdaWtrC7sUESkSmUCYmZkpyDQX+aJwWAIjIyOz5xDvv//+sMsRkTxauXLlNYeAlnIgZFM4LLG3v/3t+qQ5kWXs7rvvDruEgtC72BKLRqM89dRT1NfXh12KiCyx9vb2gk18FzaFQx5Eo1HOnDlDPB5nxYoVYZcjIjegqamJZDKJu5dNMIDCIa+6u7v55S9/ibvT0NAQdjkich3a29txd06cOFG09yLkk8KhQE6fPq37I0SKXG9v7+yF5XI6SpiPwqGAMvdHJJNJXZMQKRK1tbWzoVDIz0sodpp4LwSZaxIAHR0dDA8Ph1yRSPmpqqpiamoq7DKKlo4cQjY0NIS7k0wmufPOO8MuR6QsdHZ2KhgWoCOHIhGNRtm/fz+pVIqdO3ey1J+RLSJpHtYnq5UYHTkUmWg0yokTJ3B3ent7dUOdyBJSMORO7zxFrL+/n+npadxd90yILIF169aFXULJUDiUiMw9E5oJVmTxjh8/Tm1tbdhllASFQ4nJzATb2dkZdikiJWlycpJIJBJ2GUVP4VCiBgYGcHc2btwYdikiJWdmZgYzw8zo6OgIu5yipHAocUeOHNGpJpEbMDw8jJlRV1dHKpUKu5yioXBYBnRXp8iNO3/+PLFYTCERUDgsAxrFJLJ0MiFR7qecFA4lLJFIYGa89tprYZcisiwNDw9TU1NTlkcSCocS1dHRQU9PT9hliCx7U1NTxGIxEolE2KUUlMKhxKRSKbZu3arJ+kQKrKenZ3aE05YtW8IuJ+8UDiWgr69v9o8yFouxf//+sEsSKWsHDx6c/Te5XI8ocgoHM6s3s71m9pyZHTSzqJk1mNkTZnYoeFyTtf7bzWy/mT1rZt/L6r/HzJ43s8Nmtjur/1YzGwm29XUzq17a3Sw9XV1ds398e/bsCbscEbmKzBFFX19f2KUsrcynHl2rAQ8BHw2Wq4F6YA+wO+jbDfQHy/XAAeDNwddrg8cI8CJwW7CNHwFbgu99A/hAsPxl4OML1bRt2zZfjtra2hxQU1MrwVZXVxf2W8iCgDHP4X1/wSMHM1sN7AC+AuDuU+5+FthJOjQIHt8TLP8u8Ki7/yxY/0TQ3wYcdveX3H0K+Bqw08wMeAewd55tlYXMdQQzY3R0NOxyRGSRzp07t2xONeVyWuk24CTwoJn90MweMLNa4I3uPgEQPK4N1r8dWGNm/2RmT5vZB4P+W4BjWdsdD/oagbPufvmK/jnMrNvMxsxsrJQ/7yCVSrF+/XpdRxBZpnp6ekp+gr9cwqESeAvwJXffClwgfRrpWutvA34L6AD+wsxuB2yedf0a/XM73RPu3ururU1NTTmUXjwSiQTr1q2bDYOf//znYZckInk0OTmJmdHY2Bh2KYuSSziMA+PuPhJ8vZd0WPzCzNYBBI8nstb/trtfcPdTwJPAnUH/hqztrgdeBk4B9WZWeUV/yUulUjQ1NWFm9PT0cPz48bBLEpECe+WVVzAztm/fHnYp12XBcHD348AxM7sj6Hon6QvO+4BdQd8u4PFg+XHgbjOrNLObgO3AQeAHwKZgZFI18AFgX3CB5LvAffNsq+QkEgmqq6tnjxBOnToVdkkiUgRGR0dL61RTLletgbuAMeAZ4DFgDelrBd8BDgWPDVnr/wnpAPkJ8J+y+t8NvEB61NKfZfXfBowCh4G/BWoWqqnYRivF43GPRCKhj5ZQU1Mr7tbW1hbqexU5jlYyL9HPVG1tbfWxsbHQXj+VSvGJT3yCZ599lsuXL+uzaUUkZ2G+X5jZ0+7eutB6lQutIHN1dHRo+goRWTQzo7m5mYmJibBLuSpNn5GDRCJBY2Pj7NBTBYOI3Kjjx48X9WgmHTksoK+vT9NXiEjeZEYzFdupaYXDVWzfvl13K4tIwdTW1nLhwoWwy5il00pZurq6WLFihaaxEJGCm5ycLKpPnlM4BLq6uhgcHNSnqolIaIaHh4tmdteyD4eOjg4qKioYHBwMuxQREfbs2YOZ0dXVFWodZRsOmQ/QGR4eLroLQSIig4ODrFixIrTXL8sL0plTSCIixey1115j9erVvPrqqwV/7bI7ckgkEgoGESkZ586dY8OGDaRSqYK+blmFQyqVoqenJ+wyRESuy/j4OLFYrKABUVbh8O53vzvsEkREFu1tb3tbwV6rrMLh7NmzYZcgIrJo7k76k5Xzr6zCQURkOSjE50IoHERESszk5GTeb5Yrq3Bobm4OuwQRkSXx6KOP5nX7ZRUOf/VXfxV2CSIiS+L8+fN53X5ZhUN3dzednZ1hlyEicsOOHz+e1+2XVTgADAwMUF9fH3YZIiJFrezCAeCb3/xm2CWIiBS1sgyHaDRKMplk5cqVYZciIrIokUgkr9svy3CAdEBMTk7m/QcsIpIPX/ziF/O6/bINh4zLly9riKuIlJR4PE53d3deX6Msp+y+0sTEBB0dHQwPD4ddiojIVTU3NzMxMVGQ1yr7I4eMoaEhOjs7CzZviYjI9WhoaChYMIDC4XUGBgaYmZnRvRAiUlTi8TinT58u6GsqHOYxMDDA5s2bwy5DRMpcJBIhmUzm/frCfHIKBzOrN7O9ZvacmR00s6iZNZjZE2Z2KHhcc8Vzfs3Mps3svqy+XcH6h8xsV1b/NjP7sZkdNrPPWRGc2zlw4ADt7e1hlyEiZaq3t5fLly8TjUbDKcDdF2zAQ8BHg+VqoB7YA+wO+nYD/VnrR4D/DXwTuC/oawBeCh7XBMtrgu+NAlHAgG8B71qopm3btnmhNDc3O6CmpqaWt7Zp0yZPJpN5fz8DxjyH9/0FjxzMbDWwA/gKgLtPuftZYCfp0CB4fE/W0/4QeAQ4kdXXATzh7q+4+xngCeAeM1sHrHb3VFD4w1dsK3QTExM6ihCRvInH47zwwgvhHSXMI5fTSrcBJ4EHzeyHZvaAmdUCb3T3CYDgcS2Amd0CvBf48hXbuQU4lvX1eNB3S7B8Zf8cZtZtZmNmNnby5MkcSl86mdFMFRW6TCMiS6O5ubkg9ywsRi7vdJXAW4AvuftW4ALp00hX89+BPnefvqJ/vusIfo3+uZ3uCXdvdffWpqamhStfYgMDAzz11FMFf10RWR7a29tfd+pmYmKiKIMBcguHcWDc3UeCr/eSDotfBKeECB4zp5Baga+Z2RHgPuCLZvaeYDsbsra7Hng56F8/T39RyszLVF1dHXYpIlIiKioq6OzsZGhoKOxScrZgOLj7ceCYmd0RdL0TOADsA3YFfbuAx4P1b3X3FndvIR0kn3D3x4AhoN3M1gQjm9qBoeCU1Dkze2swSumDmW0Vq2g0ymuvvUY8Hg+7FBEpYm1tbbg709PTDAwMhF3O9cnlqjVwFzAGPAM8Rnq0USPwHeBQ8Ngwz/P+J8FopeDrDwOHg/ahrP5W4CfAi8DnAVuopkKOVlpIPB73ysrK0Ec7qKmphdsikYh3dnaG/ZZ0TeQ4WimncCjGVkzhkBGPx72hoSH0P1A1NbXCtvXr1xdkGOpSYKmGskruuru7OX36tKbfEFnG4vH4nDfSY8eOFdUw1KWgcMiDgYEBBYTIMlNRUVG0w07zQeGQJwMDA7i7QkKkxFVVVRGPx5meni6bYACFQ95lQiKZTIZdiohch8xIo6mpqbIKhQyFQ4Fk7o9YtWpV2KWIyDzMjGQyOXsdYWRkZOEnLWMKhwKKRqOcO3cOd6etrS3sckQkEI/HmZmZWXYXlW+EwiEkIyMjuDu9vb2sXLky7HJEyk5lZeXskUI5njZaiMIhZP39/UxOThKPx6mtrQ27HJGy0NnZyaVLl3SkcA0KhyLR3d3N+fPndV1CJI/Wr19PMpksvaksQqBwKDLZ1yU0RbjI0lqON6vli955itjAwADT09P09vZSX18fdjkiJa+rqyvsEkqGwqEE9Pf3c+bMmdkL2CKyOF/96lfDLqFkKBxKTH9/v0JCZJHcHTNj5cqV9PX1hV1OUVM4lKhMSOizrUWu38WLF9mzZw9mplNNV6FwKHFDQ0OamkPkBgwODmJmrF27llQqFXY5RUPhsAxo9IXIjTt58iSxWIyOjo6wSykKCodloKWlJewSRJaN4eFhzIxVq1aV9XUJhUMJS6VS1NXVcfTo0bBLEVl2Lly4MHtdohyPJhQOJaqjo4NYLMb58+fDLkVk2cscTZRTSCgcSkhXVxdmhpkxPDwcdjkiZScTEtXV1SQSibDLySuFQwnIhMLg4GDYpYgIcOnSJXp6epb1UFiFQ5FKJBK84Q1vUCiIFLnBwUHWrVsXdhlLTuFQZBKJBC0tLfT09PDqq6+GXY6I5OD48eM0NjaGXcaSUjiELJFIcOutt1JXV0dlZSU9PT0afSRSgl555RXMjLq6umVxM11l2AWUo1Qqxa5duzh06FDYpYjIEjt//jyxWAyAtra2kv0sah05FFBfXx+RSIRYLKZgECkDo6OjmFlJHkkoHPIolUrR1NQ0O/x0z549zMzMhF2WiBRYLBYruYDIKRzMrN7M9prZc2Z20MyiZtZgZk+Y2aHgcU2wbqeZPRO0pJndmbWde8zseTM7bGa7s/pvNbORYFtfN7Pqpd/Vwunr66OqqopYLMapU6fCLkdEisD73//+sEu4Pu6+YAMeAj4aLFcD9cAeYHfQtxvoD5ZjwJpg+V3ASLAcAV4Ebgu28SNgS/C9bwAfCJa/DHx8oZq2bdvmxaa3t9dra2sdUFNTU5vT6uvrPZlMhvo+BYx5Lu/7C64Aq4GfAnZF//PAumB5HfD8PM9dA/w8WI4CQ1nf+9OgGXAKqJxvvau1sMMhmUz6jh07vLKyMvQ/ODU1tdJqmzZtCi0kyDEccjmtdBtwEnjQzH5oZg+YWS3wRnefAAge187z3I8A3wqWbwGOZX1vPOhrBM66++Ur+ucws24zGzOzsZMnT+ZQen6kUil+/dd/nSeffJLLly8v/AQRkSyHDh0iFosV9RQcuYRDJfAW4EvuvhW4QPo00jWZ2W+SDofMnLc2z2p+jf65ne4Jd29199ampqYcSl86fX19sxeWY7GYLiyLyA3r6ekp2oDIJRzGgXF3zwzW3Us6LH5hZusAgscTmSeY2a8CDwA73f101nY2ZG13PfAy6VNK9WZWeUV/6LZv3/66kUYiIkutp6enKEcyLRgO7n4cOGZmdwRd7wQOAPuAXUHfLuBxADN7M/Ao8B/c/YWsTf0A2BSMTKoGPgDsC86BfRe478pthWn79u2Mjo6GXYaIlIGdO3eGXcJcuVyYAO4CxoBngMdIX2huBL4DHAoeG4J1HwDOAPuDNpa1nXcDL5AetfRnWf23AaPAYeBvgZqFasrHBen29vbQL1SpqamVb7vppps8Ho8v+XtbNpZqtFKxtqUOBwWDmppasbR8BgRLOFppWctMfKcPzxGRYvHII4+EXUL5hkMqlWLNmjX09PRw5MiRsMsREZmV+cS5MD9IqCzDoauri1gsxtmzZ8MuRUTkqgYHB6muDmc2obKbsrujo0OnkESkZFy6dAkzywzeKZiyOnJIJBIKBhEpSWbz3S+cP2UVDp/97GfDLkFEZNEKGRBlFQ4vvPDCwiuJiBSx2tragrxOWYWDJskTkVI3OTlZkNcpq3AQEVkO+vr6Fl7pBpVVOBR6JlcRkXx49NFH8/4aZRUOjz8e+nx+IiI37OWX8z9xdVmFQzQaJZlMsmrVqrBLERFZtMnJSdatW5fX1yircIB0QJw7d46VK1eGXYqIyKIdP348r9svu3DIuPvuu8MuQUSkaJVtOAwNDXHTTTeFXYaIyKJ0dnbmdftlGw4AFy5cYPPmzWGXISJyXRoaGhgYGMjra5R1OAAcOHCA9vb2sMsQEcnZpz71qby/RtmHA6RPMekIQkRKQW9vL93d3Xl/HYVDQEcQIlLs4vE4/f39BXkthUOWoaEhent7wy5DROR1IpEIyWSyIEcMGQqHK/T39+PuRCKRsEsREaGtrY3Lly8TjUYL+roKh6u4fPkybW1tYZchImUqmUzi7oyMjITy+gqHaxgZGSEej4ddhoiUkR07dpBMJgt+pHAlhcMCuru7SSaTBf+IPhEpL52dnbg73/ve90IPBoDKsAsoBdFolJmZGVasWMFrr70Wdjkiskxs3ryZAwcOhF3GvHTkcB0uXrzIxo0bwy5DREpc5iihWIMBFA7X7ciRIxruKiLXLRMI7p73qS+WQk7hYGb1ZrbXzJ4zs4NmFjWzBjN7wswOBY9rgnXNzD5nZofN7Bkze0vWdnYF6x8ys11Z/dvM7MfBcz5nRX6CPzPcVTfNici1tLe3l1QgvE6m8Gs14CHgo8FyNVAP7AF2B327gf5g+d3AtwAD3gqMBP0NwEvB45pgeU3wvVEgGjznW8C7Fqpp27ZtXiwaGhocUFNTK/NWUVHhO3bs8GQyGfbb0lUBY57D+/6CRw5mthrYAXwFwN2n3P0ssJN0aBA8vidY3gk8HNTxfaDezNYBHcAT7v6Ku58BngDuCb632t1TQeEPZ22rJJw+fZre3l6NaBIpUxs3bsTdmZ6eLprRRjcql9NKtwEngQfN7Idm9oCZ1QJvdPcJgOBxbbD+LcCxrOePB33X6h+fp38OM+s2szEzGzt58mQOpRdOf38/MzMzui9CpAxkblDLtCNHjoRd0pLLJRwqgbcAX3L3rcAF0qeRrma+/z77Ivrndron3L3V3VubmpquXXVIuru7FRAiy9TNN99cFDeoFUIu4TAOjLt75h7uvaTD4hfBKSGCxxNZ62/Iev564OUF+tfP01+yuru7dcFaZJno7e2dPUI4efJkWQQD5BAO7n4cOGZmdwRd7wQOAPuAXUHfLuDxYHkf8MFg1NJbgX8JTjsNAe1mtiYY2dQODAXfO2dmbw1GKX0wa1slbWhoCHcnHo9TV1cXdjkikqOKiorZUCjUFNlFJ5er1sBdwBjwDPAY6dFGjcB3gEPBY0OwrgFfAF4Efgy0Zm3nw8DhoH0oq78V+EnwnM8DtlBNxTRaKVfJZNKbmppCH1GhpqY2f1u5cmVRjzRaCuQ4WsnS65ae1tZWHxsbC7uMRVu9ejXnzp0LuwwRIX0tYd++fWVxysjMnnb31oXW0x3SIXn11VeJx+OsWLEi7FJEylbmruVyupaQK4VDiLq7u/nlL385exhXU1MTdkkiZSETCiV313IBKRyKyMWLF3U0IZJHZkYymVQo5EDhUGSyjyZ6e3upqqoKuySRZWNmZkanj3KkcChi/f39TE1NkUwmqa+vD7sckZLX1dUVdgklQ+FQAqLRKGfOnCGZTHLnnXeGXY5IyRocHGTLli1hl1ESFA4lJBqNsn//ftydtra2sMsRKUkHDx7EzNi+fXvYpRQ1hUOJGhkZIZlMUqxzTIkUu9HRUaqrq8Muo2gpHEpYNBrlxIkT+mQ6kUW6dOkSLS0tYZdRlBQOy0DZzv0isgSOHj2Kmc22DRs2kEqlwi4rdAqHZUDnTkWWzvj4OLFYDDMr69FNCocSlkqliEQijI6Ohl2KyLI0ODhYtiGhcChBqVSKpqYmYrEYMzMzYZcjsuxlQqKurq5sTjkpHEpIKpXive99L7FYjFOnToVdjkjZOX/+PLFYjFWrVtHX1xd2OXmlcCgBXV1dmBmxWIzHHnss7HJEyt6FCxfYs2fPsj6SUDgUqUQiwbp16zAzBgcHwy5HROaROZJYjsNhFQ5FJJVK8Ru/8RtUVFTQ09PD8ePHwy5JRHKQGQ5bUVFBIpEIu5wloXAIWebicua00ZNPPkmpfjqfSLlzd3p6eqisrCz5000Kh5C7Kk1YAAAGvElEQVSkUilqamp0cVlkGZqeniYWi5V0QCgcCiyRSFBbW0ssFmNqairsckQkj9rb28MuYdEqwy6gnCQSCXp6esIuQ0QK5Pz582GXsGg6ciigRx55JOwSRKSAmpubwy5h0RQOBfS+970v7BJEpECam5uZmJgIu4xFUzgUUHd3N/F4nLq6ugXXjUQiVFRUUFFRMWe5srKSmpqa2f751q2qqqKzsxN3p7Ozk+rq6tnvZ1pFRQVmVoA9X/y+Zh6bm5uJx+O4O+5OPB6nubl5dp3KykoqKytn962mpmZ2/92dZDLJxz72sVA/bjW7zux9zP59XPkzqKmpoa6ujpUrV855XubxrrvuIplMzu7n/fffT2dnJ/X19df1WtnLV/48M9+7XpnhnYt57mJEIpGCvM58/242b948+/fm7iUdDABWqsMmW1tbfWxsLOwyRERKipk97e6tC62nIwcREZlD4SAiInPkFA5mdsTMfmxm+81sLOi708xSQf/fm9nqoL/KzB4K+g+a2Z9mbeceM3vezA6b2e6s/lvNbMTMDpnZ181MH+wqIhKi6zly+E13vyvrXNUDwG53/7fA3wF/EvT/NlAT9G8DesysxcwiwBeAdwFbgN8xsy3Bc/qBz7j7JuAM8JEb2isREbkhN3Ja6Q7gyWD5CSAzTtOBWjOrBFYCU8CrQBtw2N1fcvcp4GvATktf9n8HsDd4/kPAe26gLhERuUG5hoMDw2b2tJl1B30/Ae4Nln8b2BAs7wUuABPAz4C/cfdXgFuAY1nbHA/6GoGz7n75in4REQlJrtNnvM3dXzaztcATZvYc8GHgc2b2X4B9pI8QIH2EMA28CVgD/B8z+0dgvgH1fo3+OYJgyoTTeTN7/ir13gyU22x22ufyoH0uD/nc5425rJRTOLj7y8HjCTP7O6DN3f8GaAcws9uB3wpW/13g2+5+CThhZv8MtJI+atiQtdn1wMukfwD1ZlYZHD1k+uerIwEsOFm6mY3lMo53OdE+lwftc3kohn1e8LSSmdWaWV1mmXQg/CQ4isDMKoA/B74cPOVnwDssrRZ4K/Ac8ANgUzAyqRr4ALDP03fhfRe4L3j+LuDxpdpBERG5frlcc3gj8JSZ/QgYBf7B3b9NerTRC6Tf+F8GHgzW/wKwivQ1iR8AD7r7M8FRwR8AQ8BB4Bvu/mzwnD7gk2Z2mPQ1iK8syd6JiMiilOz0GddiZt3BKaiyoX0uD9rn8lAM+7wsw0FERG6Mps8QEZE5SjYcrjYVR9b3d5jZ/zWzy2Z233zbKDU57PMnzeyAmT1jZt8xs5yGrBW7HPb7Y1nTuzyVded9yVpon7PWu8/M3MxKfjRPDr/n3zOzk8Hveb+ZfTSMOpdSLr9nM3t/8O/6WTP7asGKy55/vFQaEAFeBG4DqoEfAVuuWKcF+FXgYeC+sGsu0D7/JnBTsPxx4Oth112g/V6dtXwv6aHUodeez30O1qsjPUvB94HWsOsuwO/594DPh11rgfd5E/BDYE3w9dpC1VeqRw7zTsWRvYK7H3H3Z4CZMArMg1z2+bvuPhl8+X3S94yUulz2+9WsL2u5yk2UJWTBfQ78N2APcLGQxeVJrvu8nOSyz/8R+IK7n4H0vWaFKq5Uw+FqU3EsZ9e7zx8BvpXXigojp/02s983sxdJv1n+UYFqy5cF99nMtgIb3P1/FbKwPMr17/t9wWnTvWa2YZ7vl5Jc9vl24HYz+2cz+76Z3VOo4ko1HHKecmMZuZ5pRrpI35X+13mtqDBy2m93/4K7/yvS98z8ed6ryq9r7nNw4+lngP9csIryL5ff898DLe7+q8A/kp6ks5Tlss+VpE8tvR34HeABMyvIZ92WajiMM/9UHMtZTvtsZv8O+DPgXnd/rUC15dP1/q6/RunP6rvQPtcB/wb4JzM7QnoWgn0lflF6wd+zu5/O+pv+H6Q/EqCU5fK3PQ487u6X3P2nwPOkwyLvSjUc5p2KI+Sa8m3BfQ5ONcRJB0PBzk3mWS77nf2P5beAQwWsLx+uuc/u/i/ufrO7t7h7C+nrS/e6eyl/qHouv+d1WV/eS3qmhVKWy/vYY6QHmmBmN5M+zfRSIYoryXDwq0zFYWb/1czuBTCzXzOzcdLTicfN7Nmrb7H45bLPpE8jrQL+NhjqV/KBmeN+/0EwzG8/8EnS83OVrBz3eVnJcZ//KPg9/4j0daXfC6fapZHjPg8Bp83sAOk56P7E3U8Xoj7dIS0iInOU5JGDiIjkl8JBRETmUDiIiMgcCgcREZlD4SAiInMoHEREZA6Fg4iIzKFwEBGROf4fFObz41IO8FoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(dts, simData['observationStartMJD'], 'k.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "dras = (ephs[0]['ra'][eph_idxs] - simData['fieldRA']) * np.cos(np.radians(simData['fieldDec']))\n",
    "ddecs = (ephs[0]['dec'][eph_idxs] - simData['fieldDec'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "matchIdx = np.where((dras < raTol) & (ddecs < decTol))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(118064, 623453)"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(matchIdx[0]), len(simData)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "118064"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "time_acc = simData['observationStartMJD'][matchIdx]\n",
    "len(time_acc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12.1 s ± 378 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "ephs = pyo.generateEphemerides(time_acc, timeScale='TAI', obscode='I11', ephMode='nbody')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "ephs_acc = pyo.generateEphemerides(time_acc, timeScale='TAI', obscode='I11', ephMode='nbody')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "68.62777777777778"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(35.3/100 + 12) * 20000 / 60. / 60. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Then check what is really in fov? \n",
    "from lsst.sims.utils import angularSeparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [],
   "source": [
    "sep = angularSeparation(ephs_acc['ra'], ephs_acc['dec'], simData['fieldRA'][matchIdx], simData['fieldDec'][matchIdx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(179.79016681290727, 0.10970210680756547)"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sep.max(), sep.min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "accMatchIdx = np.where(sep < 2.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "93"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(accMatchIdx[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "fephs = ephs_acc[0][accMatchIdx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 37.33056666,  37.32565026,  31.91883082,  31.95051175,\n",
       "        31.95146478,  32.01076165,  32.01078091,  32.01124576,\n",
       "        32.01126605,  32.04238469,  32.0424044 ,  32.04290912,\n",
       "        32.04339749,  32.13677374,  32.13751594,  32.13757527,\n",
       "        32.13817401,  32.16907173,  32.16966973,  32.20070793,\n",
       "        32.20138792,  32.26330566,  32.26411595,  32.35517469,\n",
       "        32.35530286,  32.35610341,  32.35623843,  32.41596286,\n",
       "        32.41691657,  32.50271069,  32.50321685,  32.6744164 ,\n",
       "        32.69683777,  32.74993981,  32.74995599,  32.75061583,\n",
       "        32.75063265,  32.77539006,  32.77540518,  32.77552908,\n",
       "        32.77597443,  32.77599056,  32.77611934,  32.82693251,\n",
       "        32.87756934,  32.87759982,  32.87834405,  32.87837642,\n",
       "        32.90250044,  32.95204466,  32.97664037,  32.97666996,\n",
       "        33.00107158,  33.00110157,  33.04895552,  33.04955804,\n",
       "        33.0956765 ,  33.09576164,  33.09635497,  33.09644411,\n",
       "        33.14146075,  33.14210597,  33.2071849 ,  33.20755433,\n",
       "        33.2082808 ,  33.22856254,  33.24932799,  33.35205639,\n",
       "        33.38632439,  33.3863369 ,  33.40466871,  33.40480457,\n",
       "        33.40481709,  33.42284105,  33.4231641 ,  33.47608877,\n",
       "        33.47640882,  33.56186091,  33.56225779,  33.81697022,\n",
       "        33.81721572,  33.8400277 ,  33.84059139,  33.85165623,\n",
       "        33.85166689,  33.85194027,  33.85195129,  34.07565295,\n",
       "        34.07567176,  34.07569082,  34.08231542,  34.08261173,  34.08262141])"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fephs['magV']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "28.098116774930414"
      ]
     },
     "execution_count": 123,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ephs['magV'].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "bright_objs = []\n",
    "tstep = 2.0\n",
    "times = np.arange(54466, 61406 + tstep/2., tstep)\n",
    "pyo = mo.PyOrbEphemerides()\n",
    "for i, obj in enumerate(mini):\n",
    "    pyo.setOrbits(obj)\n",
    "    ephs = pyo.generateEphemerides(times, timeScale='TAI', obscode='I11', ephMode='2body')\n",
    "    if (ephs[0]['magV'].min() < 26.5):\n",
    "        bright_objs.append(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20265, 1342)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(mini), len(bright_objs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "br = mini[bright_objs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Number of objects')"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF9ZJREFUeJzt3X2QZXV95/H3BxTxAQYV3BBkHGQoI4spdBvQ7K7Gh2wwOGIIKrNb2dUiTjCLSSVbWcfo+hBjSULilo9ljU+YDUIoVllGR3FjBHVDIRKUB4GIBMIsJkSN48PGB+C7f9wzzJ2uc7tPz/Tpc2/3+1XV1X3PPffcT9+uvt97fr/z+/1SVUiSNN8BQweQJE0nC4QkqZUFQpLUygIhSWplgZAktbJASJJaWSAkSa0sEJKkVhYISVKrhwwdYH8cfvjhtWHDhqFjSNJMue66675ZVUcstt9MF4gNGzbwpS99aegYkjRTktzVZT+bmCRJrSwQkqRWFghJUisLhCSplQVCktTKAiFJajWTBSLJpiTbdu3aNXQUSVq1ZrJAVNX2qtqybt26oaNI0qo10wPlJGm127D1Ew/+fOd5p63oc8/kGYQkqX+eQUjSFBk/YxiaZxCSpFYWCElSKwuEJKmVBUKS1MoCIUlqZYGQJLXyMldJWiGTBr1N06Wt4ywQktSjaX3z78ImJklSq6kpEEmenOS9SS5N8sqh80jSWtdrgUjywST3Jrlp3vZTk9yW5PYkWwGq6paqOgd4CTDXZy5J0uL6PoO4ADh1fEOSA4F3A88Hjgc2Jzm+ue+FwBeAz/ScS5K0iF4LRFV9Dvj2vM0nA7dX1R1V9WPgYuD0Zv/Lq+rngP/QZy5J0uKGuIrpKODusds7gVOS/DxwBvAwYMekByfZAmwBWL9+fX8pJWmNG6JApGVbVdWVwJWLPbiqtgHbAObm5mpZk0nSEgy5mM9KGKJA7ASOHrv9eOCeAXJI0pJ1Gdcwy2Mfxg1xmeu1wHFJjklyEHAWcPlSDpBkU5Jtu3bt6iWgJKn/y1wvAq4GnpRkZ5Kzq+o+4FzgCuAW4JKqunkpx62q7VW1Zd26dcsfWpIE9NzEVFWbJ2zfwQId0ZKk4U3NSOqlsIlJkvo3kwXCJiZJ6t9MFghJUv+c7luSGrO2XkPfZrJAJNkEbNq4cePQUSTNuLX65t/FTDYx2QchSf2byQIhSeqfBUKS1GomC4TjICSpfzNZIOyDkKT+zWSBkCT1byYvc5Wk/eGlrd3MZIFwHISkSSa9+a/GBX36NpMFoqq2A9vn5uZeMXQWSVopK72CnX0QkqRWFghJUisLhCSplQVCktRqJguEI6klqX9exSRp5jmuoR8zWSAkrR2OaxiOBUKSBjALZz0z2QchSeqfBUKS1MoCIUlqZYGQJLWayQLhOAhJ6t9MXsXkOAhp9s2/imf8stVZuMJnLZjJAiFpNq3mN/7V+LvNZBOTJKl/FghJUisLhCSp1ZIKRJJHJ/nZvsJIkqbHogUiyZVJDk3yGOArwIeSvK3/aJKkIXU5g1hXVd8FzgA+VFX/Cnhev7EkSUPrcpnrQ5IcCbwEeG3PeSTNqPHLPB3TsDp0KRBvAq4AvlBV1yZ5IvC1fmNJmgW++a9uXQrEN6rqwY7pqrpj6D6IJJuATRs3bhwyhqRVzOLXrQ/inR23rZiq2l5VW9atWzdkDEla1SaeQSR5BvBzwBFJfmfsrkOBA/sOJkka1kJNTAcBj2r2OWRs+3eBM/sMJUka3sQCUVVXAVcluaCq7lrBTJKkKdClk/r9SV5cVd+B0Whq4OKq+sV+o0maRrPaeTuruYfUpUAcvrs4AFTVPyV5XI+ZJE2BSeMatHZ0KRAPJFlfVX8HkOQJQPUbS9IQ/JStcV0KxGuBLyS5qrn9TGBLf5EkSdNg0QJRVZ9K8jTg6UCA366qb/aeTJI0qC6zuQY4FXhasxb0I5Kc3HsySdKguoykfg/wDGBzc/t7wLt7SyRJmgpd+iBOqaqnJbkeHryK6aCec0mSBtalQPwkyYE0Vy4lOQJ4oNdUkrQIr7jqX5cC8Q7gY8DjkryF0TQbr+sjTJIXAacBjwPeXVWf7uN5JO3hG60m6XIV04VJrgOey+gqphdV1S1dnyDJB4EXAPdW1Qlj208F3s5o4r/3V9V5VXUZcFkzWvuPAQuEtB8mvfk78E1dTOykTnJo8/0xwL3ARcBHgH9I8uim2amLCxhdBTV+7AMZdXQ/Hzge2Jzk+LFdXocd4ZI0qIXOID7C6JP/dYz6HzLv/kcleV9V/d5CT1BVn0uyYd7mk4Hbq+oOgCQXA6cnuQU4D/hkVf11599CkrTsFprN9QXN92Pa7m/OAm4CFiwQExwF3D12eydwCvAq4HnAuiQbq+q9Lc+7hWYk9/r16/fhqSVJXXTppCbJGcC/YXQm8fmquqyq7geevI/PO/9sBKCq6h2MOsUnqqptwDaAubk554SSpJ50GUn9HuAc4EZGZwznJNnf/oGdwNFjtx8P3LOfx5QkLaMuZxDPAk6oqt3jID7MqFjsj2uB45IcA/xf4Czg33d9cJJNwKaNGzfuZwxp9fGyVS2XLlNt3AaMN/YfDdzQ9QmSXARcDTwpyc4kZ1fVfcC5wBXALcAlVXVz12NW1faq2rJu3bquD5EkLdHEM4gk2xn1OawDbknyxeauk4G/6voEVbV5wvYdwI7uUSWt5nENnvlMn4WamP54xVIskU1MktS/iU1MVXXV7i/gVuCQ5uuWZttgbGKSpP51uYrpJcAXgRcDLwGuSXJm38EkScPquuToSVV1Lzw4m+tfAJf2GWwhNjFJUv+6XMV0wO7i0PhWx8f1xiYmSepflzOITyW5gtFkfQAvxauPJGnV6zLd9++OTbURYFtVfaz3ZJKkQXWai6mqPgp8tOcsndkHodViNY9r0OzrVCCmTVVtB7bPzc29YugsUt8cQKahDNrZLEmaXgutKPeZ5vsfrlwcSdK0WKiJ6cgkzwJe2Kz4ttcaDq74JmkSm8VWh4UKxOuBrYzWanjbvPsKeE5foRZjJ7Uk9W+hJUcvBS5N8t+q6s0rmGlRdlJrlvnpWrOiyziINyd5IfDMZtOVVfXxfmNJkoa2aIFI8lZGa0Bc2Gz6rST/uqpe02syacY4pkGrTZdxEKcBJ1bVA/DgkqPXAxYISVrFug6UOwz4dvOzM+RJM+LOg/cs9b7hhx8ZMIlmUZcC8Vbg+iSfZXSp6zMZ+OzBq5gkqX9dOqkvSnIlcBKjAvHqqvr7voMtksmrmKQZtNbPaGbt9+86Wd83gMt7ziJpYLP2BqZ+zeRkfZJWlmM3Vs54kYZhC7WT9UmSWi14BpHkAOCGqjphhfJIa4pNOkuz1Ner7/1XuwULRFU9kOQrSdZX1d+tVChpGow3q0wa7LbSTS++gWkldemDOBK4OckXgR/s3lhVL+wt1SK8zFVaur2Ky9Z9/zQ9v428bZ+1YrUX7C4F4k29p1giL3PVSutyNrG/bxZ7v/HuWvLjpeXWZRzEVUmeABxXVX+R5BHAgf1Hk1anSZ/ANYxp+HsslGHIs5Quk/W9AtgCPAY4FjgKeC/w3H6jSdJsWi1NT12amP4zo9lcrwGoqq8leVyvqaSBeL3/8puGT+jaN10KxI+q6sfJaMXRJA9htKKcpCnkG/J0meW/R5cCcVWS3wMenuQXgN8AtvcbS9JKmeU3MPWrS4HYCpwN3Aj8OrADeH+foaQ+uKBPv1ZLu7v26HIV0wPNIkHXMGpauq2qbGKSpFWuy1VMpzG6aunrjKb7PibJr1fVJ/sOJ2mySQPf7jx4iDTTwyaz5dOlielPgGdX1e0ASY4FPgFYIKSe7DUwzzc8DaRLgbh3d3Fo3AHc21OeTpxqQ10Nddmqo6L3z1KL4koW0bVUsCcWiCRnND/enGQHcAmjPogXA9euQLaJnGpDkvq30BnEprGf/wF4VvPzPwKP7i2RJGkqTCwQVfXylQwi7Y8uk+lJWpouVzEdA7wK2DC+/5DTfUvQf/+C1/VrrevSSX0Z8AFGo6cf6DeOJGladCkQP6yqd/SeRJI0VboUiLcneQPwaeBHuzdW1V/3lkqSNLguBeIpwK8Cz2FPE1M1t6U1wXENa9taGvswrkuB+GXgiVX1477DSMutj45mRzkvj/mvnRcCTJ8DOuzzFeCwvoNIkqZLlzOIfwHcmuRa9u6D8DJXrbiVvLRVWuu6FIg39J5CkjR1uqwHcdVKBJEkTZcuI6m/x541qA8CHgr8oKoO7TOY1p6+p8twZLR266MpcTU2T3Y5gzhk/HaSFwEnL3eQJE8EXgusq6ozl/v4kqSl6dIHsZequizJ1i77Jvkg8AJGa0qcMLb9VODtwIHA+6vqvKq6Azg7yaVLzSQt1UqOa5jfsb7WV3zT7OjSxHTG2M0DgDn2NDkt5gLgXcCfjh3vQODdwC8AO4Frk1xeVV/teExJ0grocgYxvi7EfcCdwOldDl5Vn0uyYd7mk4HbmzMGklzcHK9TgUiyBdgCsH79+i4P0Qxa69N3r8b2bM2eLn0Qy70uxFHA3WO3dwKnJHks8BbgqUleU1VvnZBnG7ANYG5uruuZjCRpiRZacvT1CzyuqurN+/icmXC8bwHn7OMxJUnLbKEziB+0bHskcDbwWGBfC8RO4Oix248H7lnKAZJsAjZt3LhxHyNoJXRpJlqukdFOpqdZMwvNiBPnYqqqP9n9xahJ5+HAy4GLgSfux3NeCxyX5JgkBwFnAZcv5QBVtb2qtqxbt24/YkiSFrLgZH1JHpPkD4AbGJ1tPK2qXl1V93Y5eJKLgKuBJyXZmeTsqroPOBe4ArgFuKSqbt6v30KStOwW6oM4HziD0dnDU6rq+0s9eFVtnrB9B7Bjqccby2YT0xo1C6fl881i5tXIv8PSLXQG8V+AnwZeB9yT5LvN1/eSfHdl4rWziUmS+jfxDKKquqwVIUlapSwCkqRWS56LaRrYBzGsrqOcvfR0eq2G9vi+f4dpfI1W+n9qJs8g7IOQpP7NZIGQJPXPAiFJamUfhCaaNA3GpHbQBdc9eOOe5sDx1dwmrfK21/at49sXjb1ke/WpdDj+9LdNS8tjJs8g7IOQpP7NZIGQJPXPAiFJamWBkCS1spN6jZrYAb0Gl/fU8rPTfHWYyTMIO6klqX8zWSAkSf2zQEiSWlkgJEmt7KRejcZGLfPGfZ/xcVo6GtfirJ1rkX+H6TOTZxB2UktS/2ayQEiS+meBkCS1skBIklpZICRJrSwQkqRWFghJUivHQcyCCeMaJq/41n6YvVdOm3TNef/jJibt53Xw0nSZyTMIx0FIUv9mskBIkvpngZAktbJASJJaWSAkSa0sEJKkVhYISVIrC4QkqZUFQpLUygIhSWrlVBsr5Y0TRn2PLwk6ts+GH37kwZ8nTZ0xPjXF+P6T9DGVxd7Tdyz74ffJNGaadU6DsjbN5BmEU21IUv9mskBIkvpngZAktbJASJJaWSAkSa0sEJKkVhYISVIrC4QkqZUFQpLUygIhSWplgZAktbJASJJaWSAkSa2mZjbXJI8E3gP8GLiyqi4cOJIkrWm9nkEk+WCSe5PcNG/7qUluS3J7kq3N5jOAS6vqFcAL+8wlSVpc301MFwCnjm9IciDwbuD5wPHA5iTHA48H7m52u7/nXJKkRfRaIKrqc8C3520+Gbi9qu6oqh8DFwOnAzsZFYnec0mSFjdEH8RR7DlTgFFhOAV4B/CuJKcB2yc9OMkWYAvA+vXr9z3FEld4637c/Xx8Y6kreE3cfz8yzCpXP5OWxxAFIi3bqqp+ALx8sQdX1TZgG8Dc3FwtczZJUmOIppydwNFjtx8P3DNADknSAoYoENcCxyU5JslBwFnA5Us5QJJNSbbt2rVr8Z0lSfuk78tcLwKuBp6UZGeSs6vqPuBc4ArgFuCSqrp5Kcetqu1VtWXdurXXvi5JK6XXPoiq2jxh+w5gR5/PLUnaPzN5OalNTJLUv5ksEDYxSVL/ZrJASJL6N5MFwiYmSepfqmZ3rFmSfwTuGujpDwe+OdBzd2XG/Tft+WD6M057Ppj+jMud7wlVdcRiO810gRhSki9V1dzQORZixv037flg+jNOez6Y/oxD5ZvJJiZJUv8sEJKkVhaIfbdt6AAdmHH/TXs+mP6M054Ppj/jIPnsg5AktfIMQpLUygLRQZKDk3wxyVeS3JzkTc32Y5Jck+RrSf68mZ12mvKd26z7XUkOHyJbh4wXNuuT39SsYf7QKcz4gWbbDUkuTfKoaco3dv87k3x/iGxjGSa9hhck+dskX26+TpyyfEnyliR/k+SWJL85RL5FMn5+7PW7J8llvYepKr8W+WK0yNGjmp8fClwDPB24BDir2f5e4JVTlu+pwAbgTuDwKX0Nf6m5L8BFQ72Gi2Q8dGyftwFbpylfc3sO+B/A96f073wBcOaQ2RbJ93LgT4EDmvseN20Z5+3zP4H/2HcWzyA6qJHdn8we2nwV8Bzg0mb7h4EXDRBvYr6qur6q7hwi03wLZNzR3FfAF9mzLvk0ZfwujD5lAg9n9LefmnxJDgTOB/7rELnGLfC/MhUWyPdK4Per6oFmv3sHirjoa5jkEEbvPb2fQVggOkpyYJIvA/cC/xv4OvCdGq1vAaOV8o6alnxVdc1QWSZZKGPTtPSrwKeGytfkaM2Y5EPA3wM/A7xzyvKdC1xeVd8YKte4Bf7Ob2ma6f57kodNWb5jgZcm+VKSTyY5bqh8C2Tc7ZeBz+z+4NInC0RHVXV/VZ3I6BPuycCT23Zb2VRjTzwvX5IThsoyySIZ3wN8rqo+P0y6kUkZq+rlwE8zWuTqpVOU75nAixmwaM034TV8DaPiehLwGODVU5bvYcAPazRa+X3AB4fKt0DG3TYzao7tnQViiarqO8CVjNotD0uye9GlqVhbeyzfqQNHmWh+xiRvAI4AfmfAWHtpex2r6n7gz4FfGSjWg8byPRvYCNye5E7gEUluHzDag8Zfw6r6RtN08iPgQ4w+ZA1q3t94J6N2fYCPAT87UKy9tPyvPJbRa/eJlXh+C0QHSY5Icljz88OB5zH6JPlZ4Mxmt/8E/K8pynfrEFkmmZQxya8Bvwhs3t3+O2UZb0uysdkWYBMDvbYT8l1XVT9VVRuqagPw/6pq4xD5Fsh4a5Ijm21h1Fd30zTlY9Se/5xmt2cBfzNEvibXQv/PLwY+XlU/XIksvS45uoocCXy46Qw8gNE62h9P8lXg4iR/AFwPfGDK8v0mo47LnwJuSLKjqn5tyjLex2hG3qtH7x18tKp+f1oyMvqk9vkkhzK6uuQrjDo0pyJfVX18oCyTTPo7/2WSIxi9hl8GzpmyfF8ALkzy28D3gaH+TyZmbO47CzhvpYI4klqS1MomJklSKwuEJKmVBUKS1MoCIUlqZYGQJLWyQGjVyryZTZO8LMm7hsozLsmOJIc1X7+xD4//+STTdomrVhkLhDSAqvqlZpTsYcCSC4S0EiwQWvOSvDHJh5N8OsmdSc5I8kdJbkzyqWYiQZK8Psm1Ga1dsa0ZFUySk5pJ6K5Ocn6Sm5rtL0vy0eYYX0vyR2PPeWdGa3ScBxyb0Rz/588/M0jyriQva34+NcmtzaCuM8b2eWRGa2lcm+T6JKevxOum1c8CodXs4dmzwMqXgYVGaB8LnAacDvwZ8Nmqegrwz812gHdV1UlVdQKjab9f0Gz/EHBOVT0DuH/ecU9kNLnfUxjNFnr0vPu3Al+vqhOr6ncnhUtyMKNJ5DYB/5bR6PjdXgv8ZVWdxGhupvOTPHKB31XqxAKh1eyfmzfeE5uZMV+/wL6frKqfADcCB7Jn2vEbGS26BPDsjFYQvJHRvD3/spkz55Cq+qtmn4/MO+5nqmpXM3fOV4En7OPv8jPA31bV15q1M/5s7L5/B2xtiuCVwMHA+n18HulBzsUkjfwIoKoeSPKT2jMHzQPAQ5pP8O8B5qrq7iRvZPRGnC7HbdzP4v9z97H3B7eDx36eNC9OgF+pqtsWOba0JJ5BSN3sfqP+ZkZrUp8JUFX/BHwvydOb+89a4nG/Bxwydvsu4PgkD0uyDnhus/1W4Jgkxza3N4895grgVWN9Ik9dYgaplQVC6qC54uh9jJqcLgOuHbv7bGBbkqsZfZrftYTjfgv4P03H9/lVdTejWWRvAC5kNEswTRPVFuATTSf1XWOHeTOjZSlvaDrI37xvv6W0N2dzlfZTkkftXkM4yVbgyKr6rYFjSfvNPghp/52W5DWM/p/uAl42bBxpeXgGIUlqZR+EJKmVBUKS1MoCIUlqZYGQJLWyQEiSWlkgJEmt/j9KSeX6uNFogwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "n, b, p = plt.hist(mini.orbits.H, bins=100)\n",
    "n, b, p = plt.hist(br.orbits.H, bins=b)\n",
    "plt.semilogy()\n",
    "plt.xlabel('H magnitude')\n",
    "plt.ylabel('Number of objects')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Number of objects')"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFvxJREFUeJzt3X20JHV95/H3F0bEKMyAiGERHGCMG4wJwoXIJsGnjRnFESX4QPbsEXWZoNHDml1PhuAiCSfrAzHnLIsmGeMDJhFEDmEZBRHdgHElPIzIk0gYEJYJLoNGB+KJoM53/6jfZXpuqu+t+1Bd1X3fr3P63Orq6upPV9ftb//qVw+RmUiSNNNuXQeQJPWTBUKSVMsCIUmqZYGQJNWyQEiSalkgJEm1LBCSpFoWCElSLQuEJKnWiq4DLMZ+++2Xq1ev7jqGJI2VzZs3fzcznzHXdGNdIFavXs1NN93UdQxJGisRcX+T6dzEJEmqZYGQJNWyQEiSao1lgYiIdRGxcfv27V1HkaSJNZYFIjM3Zeb6lStXdh1FkibWWBYISVL7LBCSpFoWCElSrbE+UE4aB6s3fP6J4fvef3yHSaT5sQUhSaplC0JSb9n66pYFQr3X5EtiMV8kfglJ9SwQDfklokm1XNbtPr/PvmazQKi3K6cmg+vX+LJALNIoVv62N7Fo/LkZbjJ1vdwtECPQ9ENue2XoemWbqW95JtHgMh7la/l5TgYLhBqZ7z9/X74s+pJjLsNyzvyC7/N7mGTjsh4tNQuEFmWpNm2M4ouwb5tYRvnrvi19eQ99/gJfzI+rrlkgxljbK9J85z+KFXuUXwRN3k/TX/5L/bp9t5h1Zyk/18V8hn3TRc6xLBARsQ5Yt2bNmiWZ37isIJIqfSmik/7dMZYFIjM3AZumpqZO7TrLfE3ar2yNVh8+266+nPtSFIYZ501Jw4xlgVBzfVsJ+5ZHo7VcPv9JeZ8WiAUY9uGP+piIroy6FTRJxul9jVPWNiz39w/LuED44Y+/NjqCXS/UpnFbv5ZtgWiiyz1RutpDadxW4Dp927tL/TPJ6/9SskDM4AqihRr1umMLqhvLaRl5wSBJUi1bENKEGfdfuOOef5LYgpAk1bIFIWmi2AJZOrYgJEm1LBCSpFoWCElSLQuEJKmWBUKSVMsCIUmqZYGQJNWyQEiSalkgJEm1LBCSpFoWCElSLQuEJKmWBUKSVKs3BSIiXhMRH42I/xURL+86jyQtd60WiIj4eERsi4jbZ4xfGxF3RcSWiNgAkJmXZeapwCnAG9rMJUmaW9stiE8CawdHRMTuwIeBVwCHAydHxOEDk7ynPC5J6lCrBSIzvwL804zRxwBbMvPezHwcuAg4ISofAK7MzK+3mUuSNLcu+iAOBB4YuL+1jHsn8O+BkyLitGFPjoj1EXFTRNz08MMPt5tUkpaxLi45GjXjMjPPA86b68mZuRHYCDA1NZVLnE2SVHTRgtgKHDRw/1nAgx3kkCTNoosCcSPwnIg4JCL2AN4IXN5BDknSLNrezfVC4DrguRGxNSLempk/Ad4BXAXcCVycmXfMc77rImLj9u3blz60JAlouQ8iM08eMv4K4IpFzHcTsGlqaurUhc5DkjS73hxJLUnqFwuEJKmWBUKSVGssC4Sd1JLUvrEsEJm5KTPXr1y5susokjSxxrJASJLaZ4GQJNWyQEiSalkgJEm15lUgImKfiPjFtsLMI4d7MUlSy+YsEBFxTUTsHRH7ArcAn4iIP2k/2nDuxSRJ7WvSgliZmY8AJwKfyMyjqC7sI0maYE0KxIqIOAB4PfC5lvNIknqiSYH4A6pTc2/JzBsj4lDg7nZjSZK61uR039/JzCc6pjPz3q77ICRJ7WvSgvifDceNjHsxSVL7hrYgIuJY4N8Bz4iI3x14aG9g97aDzcYLBklS+2bbxLQH8LQyzV4D4x8BTmozlCSpe0MLRGZeC1wbEZ/MzPtHmEmS1ANN+iD+IiJWTd8pR1Nf1WImSVIPNCkQ+2XmD6bvZOb3gf3biyRJ6oMmBWJHRBw8fScing1ke5EkSX3Q5DiIM4GvRsS15f5xwPr2Is0tItYB69asWdNlDEmaaHO2IDLzC8CRwGeAi4GjMrPTPghP1idJ7WtyNtcA1gJHluMPfiYijmk9mSSpU036ID4CHAucXO4/Cny4tUSSpF5o0gfxy5l5ZETcDNVeTBGxR8u5JEkda9KC+HFE7E7ZcykingHsaDWVJKlzTQrEecDfAPtHxB8BXwX+e6upJEmdm3MTU2b+dURsBl4GBPCazLyz9WSSpE7NdjbXvTPzkXIt6m3AhQOP7QM8kpk/HUFGSVIHZmtBfBp4FbCZqv8hZjz+tIj4aGb+flvhhvFAOUlq39A+iMx8Vfl7SGYeWv4+cQN+FnjtqILOyOaBcpLUsia7uRIRJwK/StWS+LvMvKxsXvr5NsNJkrrT5EjqjwCnAbcBtwOnRYQHyknShGvSgngR8AuZOX0cxAVUxUKSNMGaHAdxF3DwwP2DgFvbiSNJ6ovZdnPdRNXnsBK4MyJuKA8dA3xtBNkkSR2abRPTH48shSSpd4YWiMycvkAQEfFM4Ohy94bM3NZ2MElSt5rsxfR64AbgdcDrgesj4qS2g0mSutX0kqNHT7caytlcvwRc0mYwSVK3muzFtNuMTUrfa/g8SdIYa9KC+EJEXMXOk/W9AbiivUhz81xMktS+OVsCmflu4M+BXwR+CdiYmb/XdrA5MnkuJklqWaNzMWXmpcClLWeRJPWIfQmSpFoWCElSraEFIiK+XP5+YHRxJEl9MVsfxAER8SLg1RFxETOuKJeZX281maSRuG/P33piePWPPt1hEvXNbAXiLGAD8CzgT2Y8lsBL2wolLTeT8CU9Ce9Bu5rtXEyXAJdExH/LzHNGmEmS1ANz7uaamedExKuB48qoazLzc+3GkiR1rcnJ+t4HnA58s9xOL+MkSROsyYFyxwNHZOYOeOKSozcDZ7QZTJLUrabHQawaGPb8FpK0DDRpQbwPuDki/pZqV9fjsPUgSRO/51aTTuoLI+IaqivKBfB7mfn/2g4mSepW05P1fQe4vOUsy9Kk/wKRhnHd779GBaJvvB6ENFksFv00lifr83oQktS+WQtEROwWEbePKowkqT9m3cSUmTsi4paIODgz/++oQmm0JqF5P/geBo3r+5H6oEkfxAHAHRFxA/DD6ZGZ+erWUkmSOtekQPxB6ykEDP8lPwm/8MfFKJf1sFaPn3e/DfvcJlGT4yCujYhnA8/JzC9FxM8Au7cfTZLUpTkLREScCqwH9gUOAw4E/gx4WbvRtNTa+OWzVL92Z2YblxZUn7NNApdvt5rs5vo7wK8AjwBk5t3A/m2GkiR1r0mBeCwzH5++ExErqK4oJ0maYE06qa+NiN8HnhIRvw68HdjUbqzxtNybw+P0/scp6zB9eA/LqcN2OWrSgtgAPAzcBvw2cAXwnjZDSZK612Qvph3lIkHXU21auisz3cTUssX8MlvML8s2Op37/gu9z533Gq22P7dxO6CzyV5Mx1PttXQP1em+D4mI387MK9sOJ0nqTpM+iA8BL8nMLQARcRjwecACIUkTrEmB2DZdHIp7gW0t5emtvjcNJ62zcNLeT9uW46bBYZbqjARtHeMzToYWiIg4sQzeERFXABdT9UG8DrhxBNkkSR2arQWxbmD4IeBFZfhhYJ/WEkmSemFogcjMN48yyKi13XychKb6Ulnuy6LL999kc8ti5tkXfcvUtzwL1WQvpkOAdwKrB6f3dN+SNNmadFJfBnyM6ujpHe3GkST1RZMC8aPMPK/1JD00Kc3EpdD2wWRdPH8cdPke+/aZt7Wprg/rUV83wzYpEP8jIt4LfBF4bHpkZn69tVSSpM41KRDPB/4j8FJ2bmLKcl+SNKGaFIjXAocOnvJbkjT5mhSIW4BVtHz0dEQcCpwJrMzMk9p8rZn6uv2vzjhl7YM+bF8eteXwnu2/Go0mp/t+JvCtiLgqIi6fvjWZeUR8PCK2RcTtM8avjYi7ImJLRGwAyMx7M/Ot838LkqQ2NGlBvHcR8/8kcD7wqekREbE78GHg14GtwI0RcXlmfnMRryNJWmJNrgdx7UJnnplfiYjVM0YfA2zJzHsBIuIi4ATAAiFJPdLkSOpH2XkN6j2AJwE/zMy9F/iaBwIPDNzfCvxyRDwd+CPgBRFxRma+b0ie9cB6gIMPPniBEUbLfoOd3Pa7PPRlnR/l+jaJ63aTFsReg/cj4jVUrYCFivqXye8BpzXIsxHYCDA1NeWV7SSpJU06qXeRmZexuGMgtgIHDdx/FvDgIuYnSWpBk01MJw7c3Q2YYucmp4W4EXhOOQngPwJvBCavbSZJY67JXkyD14X4CXAfVafynCLiQuDFwH4RsRV4b2Z+LCLeAVwF7A58PDPvmE/oiFgHrFuzZs18ntYLS7mdsu1tnpO4TVWjN9/1yPWuP5r0QSz4uhCZefKQ8VcAVyxivpuATVNTU6cudB6SpNnNdsnRs2Z5XmbmOS3kkST1xGwtiB/WjHsq8Fbg6YAFQpIm2GyXHP3Q9HBE7AWcDrwZuAj40LDnSZImw6x9EBGxL/C7wH8ALgCOzMzvjyLYbMa5k3rUJq3Db9zfz7jnV/v6cpAhzHIcREScS7VL6qPA8zPz7D4UB6g6qTNz/cqVK7uOIkkTa7YD5f4L8G+A9wAPRsQj5fZoRDwymniSpK7M1gcx76OsJUmTwyIgSarV5Ejq3mmzk9pORE2bb2dh39edvufTv9Z1h/VYtiDspJak9o1lgZAktc8CIUmqZYGQJNWyQEiSarkXkzRP7g20/CzXz3wsWxDuxSRJ7RvLAiFJap8FQpJUywIhSaplgZAk1bJASJJqWSAkSbXGskBExLqI2Lh9+/auo0jSxBrLAuFxEJLUvrEsEJKk9lkgJEm1LBCSpFoWCElSLQuEJKmWBUKSVMsCIUmq5QWDpAaW6wVjtLyNZQvCA+UkqX1jWSAkSe2zQEiSalkgJEm1LBCSpFoWCElSLQuEJKmWBUKSVMsCIUmqZYGQJNUaywLhNaklqX1jWSA81YYktW8sC4QkqX0WCElSLQuEJKmWBUKSVMsCIUmqZYGQJNWyQEiSalkgJEm1LBCSpFoWCElSLQuEJKmWBUKSVMsCIUmqZYGQJNWyQEiSaq3oOsBCRMQ6YN2aNWu6jiJJI3Hfnr81cG80F0sbyxaEFwySpPaNZYGQJLXPAiFJqmWBkCTVskBIkmpZICRJtSwQkqRaFghJUi0LhCSpVmRm1xkWLCIeBu7vOkeN/YDvdh1iHszbLvO2y7zz9+zMfMZcE411geiriLgpM6e6ztGUedtl3naZtz1uYpIk1bJASJJqWSDasbHrAPNk3naZt13mbYl9EJKkWrYgJEm1LBBziIj7IuK2iPhGRNxUxp0dEf9Yxn0jIl45MP0ZEbElIu6KiN8YGL+2jNsSERsGxh8SEddHxN0R8ZmI2GMps5bx7yyvfUdEfLAPWYflLfOdXq73RcQ3ep73iIj4++lxEXFMGR8RcV7JdGtEHDkwnzeVTHdHxJsGxh9V5r+lPDdayPtLEXFdGb8pIvYemL7r5bsqIi6JiG9FxJ0RcWxE7BsRV5fXuDoi9inT9mH51uV9Xfk/2xERUzOm73T5LkhmepvlBtwH7Ddj3NnAf62Z9nDgFuDJwCHAPcDu5XYPcCiwR5nm8PKci4E3luE/A962xFlfAnwJeHK5v38fsg7LO+PxDwFn9Tkv8EXgFWX4lcA1A8NXAgG8ELi+jN8XuLf83acM71MeuwE4tjznyun5LnHeG4EXleG3AOf0aPleAPynMrwHsAr4ILChjNsAfKBHy7cu788DzwWuAaYGpu18+S7kZgtiaZ0AXJSZj2Xmt4EtwDHltiUz783Mx4GLgBPKL5iXApeU518AvGaJM70NeH9mPgaQmdt6nPUJ5fVeD1zY87wJTP8KXwk8OJD3U1n5e2BVRBwA/AZwdWb+U2Z+H7gaWFse2zszr8vqG+FTLeV9LvCVMnw18JsDeTtbvqUlcxzwMYDMfDwzf1ByXVDzGp0u32F5M/POzLyr5il9XX9nZYGYWwJfjIjNEbF+YPw7StP249PNXuBA4IGBabaWccPGPx34QWb+ZMb4pcz6c8CvlabqtRFxdE+yDss77deAhzLz7p7n/c/AuRHxAPDHwBkLzHtgGZ45fqnz3g68ugy/DjhogXmXevkeCjwMfCIibo6Iv4iIpwLPzMzvAJS/+y8w71Iv32F5h+l6+S6IBWJuv5KZRwKvAH4nIo4D/hQ4DDgC+A7VphComq4z5QLGL2XWFVRN7RcC7wYuLr9Ous46LO+0k9nZemABuUaV923AuzLzIOBdlF+UPc77ljK8GdgLeLwneVcARwJ/mpkvAH5ItUlpGPOOgAViDpn5YPm7Dfgb4JjMfCgzf5qZO4CPUjUToaryBw08/VlUmxyGjf8uVdN4xYzxS5a1vPalpSl+A7CD6lwwnWadJS/lNU4EPjMweV/zvgm4tEzyWRa+LmwtwzPHL2nezPxWZr48M4+iKsD3LDDvUi/frcDWzLy+3L+E6gv4obJ5iPJ328D0XS7fYXlnm77T9XdBRt3pMU434KnAXgPDXwPWAgcMTPMuqm2LAM9j146oe6k6oVaU4UPY2RH1vPKcz7JrR9TblzjracAflvE/R9WcjS6zzpa33F8LXDtj+l7mBe4EXlzGvwzYXIaPZ9dO1BvK+H2Bb1O16vYpw/uWx24s0053or6yhbzTOynsRrUd/i19WL5lHn8HPLcMnw2cW26DndQf7MPyHZZ34LFr2LWTuvPlu6D3OOoXHKcb1XbGW8rtDuDMMv4vgduAW4HL2bVgnEn1q+wuBvaSoNrr4h/KY2fOeI0bqDqtPkvZ22gJs+4B/BXVtuevAy/tOutsectjnwROq3lO7/ICvwpsLuOvB44q4wP4cMl024wvi7eUTFuANw+Mnyqf0z3A+ZQDWZc47+llWf0D8P7B1+hy+Zb5HQHcRPV/dRnVF/zTgS8Dd5e/01/2nS7fWfK+lqpV8BjwEHBVX5bvQm4eSS1JqmUfhCSplgVCklTLAiFJqmWBkCTVskBIkmpZIDSxIuKfZ9w/JSLO7yrPoIi4opwNdFVEvH0Bz39xRHyujWzSNAuE1IHMfGVWJ6NbBcy7QEijYIHQshfV9T0uiIgvRnUNhRMj4oPl2gFfiIgnlenOiogbI+L2iNg4fT2BiDi6nLjxuog4NyJuL+NPiYhLyzzujl2vxXFfROxHdbDaYVFds+HcmS2DiDg/Ik4pw2vLtQe+SnUqkulpnlpOGnljOXHcCaNYbpp8FghNsqfEzosPfQP4w1mmPYzq9A0nUB15/reZ+XzgX8p4gPMz8+jM/AXgKcCryvhPUB35fSzw0xnzPQJ4A/B84A0RcdCMxzcA92TmEZn57mHhImJPqvN+raM60+3PDjx8JvC/M/Noqut/nDvHmUWlRiwQmmT/Ur54j8jMI4CzZpn2ysz8MdVpG3YHvlDG3wasLsMvKadNv43qXP3Pi4hVVOc8+lqZ5tMz5vvlzNyemT8Cvgk8e4Hv5d8C387Mu7M6/cFfDTz2cmBDKYLXAHsCBy/wdaQnrJh7EmlZmL6g0o6I+HHuPAfNDmBF+QX/Eapz/jwQEWdTfRHPddnKxwaGf8rc/3M/YdcfbnsODA87L04Av5n1F6qRFswWhNTM9Bf1dyPiacBJAFldtezRiHhhefyN85zvo1TXZZh2P3B4RDw5IlZSnSEW4FvAIRFxWLl/8sBzrgLeOdAn8oJ5ZpBqWSCkBsoeRx+l2uR0GdWpo6e9FdgYEddR/ZrfPo/5fg/4P6Xj+9zMfIDqWsS3An8N3Fym+xGwHvh86aS+f2A25wBPAm4tHeTnLOxdSrvybK7SIkXE0zLzn8vwBqrTv5/ecSxp0eyDkBbv+Ig4g+r/6X7glG7jSEvDFoQkqZZ9EJKkWhYISVItC4QkqZYFQpJUywIhSaplgZAk1fr/8jTqyvtpR6oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "n, b, p = plt.hist(mini.orbits.epoch, bins=100)\n",
    "n, b, p = plt.hist(br.orbits.epoch, bins=b)\n",
    "plt.semilogy()\n",
    "plt.xlabel('H magnitude')\n",
    "plt.ylabel('Number of objects')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "br.orbits.to_csv('allbright_NESCv9reintv1.TCO.withH.kep.des', index=False, sep=' ')"
   ]
  },
  {
   "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}