{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# AsteroGaP: A Guide \n", "\n", "## Motivation\n", "\n", "The motivation for this code stems from a drive to incorporate new computational methods with new astronomical problems. \n", "\n", "Asteroids are leftover remnents from the creation of the solar system. As they orbit the Sun, they also rotate along an axis. This rotation makes it so that the light reflected and measured off an asteroid will vary with time, meaning if we monitor the asteriod for long enough, we will be able to measure a complete rotation and figure out what the asteroid's rotational period is (an example lightcurve is plotted in the cell below). This rotational period tells us a long about that the asteroid might be made off and how we can expect it to behave in the future.\n", "\n", "Traditionally, these rotational periods were measured by training a telescope on an asteroid for several hours, and observing as frequently as possible so that you could see the lightcurve." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, 'Phaethon 3200')" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOy9d5hkV3nn/3krx66O05PzKIxyAAkLSeRowOBE8GKz2DILmPX+lvU6go0DXsP+1ssaY4PNYgwGY4wBI4EQSCCB4gjFkTQ59nQO1ZXj2T/uvdU1o+7qCvdW6D6f5+lnum/d6jp9597zPW847ytKKTQajUajuRBXpweg0Wg0mu5EC4RGo9FolkULhEaj0WiWRQuERqPRaJZFC4RGo9FolkULhEaj0WiWRQuEZk0hIj8QkV9t02f9oYh8oR2fpdF0Ai0Qmp5DRE6KSEZEkiIyKSKfE5GIw5/5EhE56+Rn1PjsYRH5sYjMisiCiDwgIjdVvf7LIvKoiCyKyFkR+QsR8VS9Pigi/yYiKRE5JSJvv+D3v908nhKRr4vIYDv/Pk33ogVC06u8QSkVAa4Frgd+v8PjcZIk8B+BEWAA+B/Av1eJQAj4TWAYuAF4OfDBqvd/EsgDo8A7gE+JyGUA5r9/C/wH8/U08NcO/z2aHkELhKanUUqNAd8GLq86vMNccSdE5LsiMmy9ICL/IiITIhIXkXutidJ8zS8iHxeR06Zl8jciEhSRsPkZm02rJSkim823+UTk8+ZnHRSR66t+36Wmy2vBfO2NVa99TkQ+KSK3m+99SET2rPA3ZpVSh5RSZUCAEoZQDJqvf0opdZ9SKm9ejy8CN5mfEwZ+FvgDpVRSKfUj4JsYggCGYPy7UupepVQS+APgLSISbfT/QrP20AKh6WlEZBvwOuCxqsNvB94FbAB8nL+a/jawz3ztJxiTqcWfAxcBVwN7gS3Ah5RSKeC1wDmlVMT8Ome+543Al4F+jIn3r8xxeYF/B75rftZvAF8UkYurPu+twB9hTPZHgT9d5W99Esian/N3SqmpFU69BThofn8RUFRKHa56/QnAEsbLzJ8BUEodw7A2Lqo1Fs36QAuEplf5uogsAD8Cfgj8WdVr/1cpdVgplQG+gjHhA6CU+qxSKqGUygF/CFwlIjEREeA24L8opeaUUgnzd751lXH8SCl1h1KqBPwjcJV5/EYgAvy5ubK/G/gW8Laq9/6bUuphpVQRQ6iupgZKqSuBPgwB/NFy54jIf8RwuX3cPBQBFi84LQ5Eq16P13hds47xrH6KRtOV/IxS6nsrvDZR9X0aYxJERNwYq/Sfx/Dnl81zhgE/hi//UUMrAMOd415lHBd+VsCMDWwGzphuIYtTGFZJzXHWQimVBb4kIs+KyONKqcrqX0R+Bvgo8Aql1Ix5OIkhKtX0AYk6X9esY7QFoVlPvB14E/AKIAbsNI8LMANkgMuUUv3mV8wMhAM0Wvb4HLBNRKqfse3AWLODvwAvsNv6QUReA3wGI3j/VNV5hwGPiOyrOnYVSy6ogyxZPYjIbgyxrHZJadYpWiA064kokANmMayFilvKXOl/BvhfIrIBQES2iMirzVMmgSERidX5WQ9hWAW/JSJeEXkJ8AaMeEVDiMiNIvJiEfGZQfP/jpFx9JD5+sswXFQ/q5R6uPq9Zvzka8BHRCRspse+CcMdhvm+N4jIzWZA+yPA10wXm2adowVCs574PIabZwx4Bnjwgtf/O0aw+EERWQS+B1wMoJR6DvgScNzMStpMDZRSeQxBeC2GdfLXwDvN39MofoxU1Vlz7K8DXl8VKP8DDIvojqosq29Xvf+9QBCYMv+G/6SUOmiO8yDwHgyhmMIQ0fc2MUbNGkR0wyCNRqPRLIe2IDQajUazLFogNBqNRrMsWiA0Go1GsyxaIDQajUazLGtmo9zw8LDauXNnp4eh0Wg0PcWjjz46o5QaWe61NSMQO3fu5MCBA50ehkaj0fQUInJqpde0i0mj0Wg0y6IFQqPRaDTLogVCo9FoNMuiBUKj0Wg0y6IFQqPRaDTLogVCo9FoNMuiBUKj0Wg0y6IFQqPRrAkmF7N87scnyBfLq5+sqYs1s1FOo9Gsb377X5/knkPThPwefuH6bZ0ezppAWxAajWZNcHgyCcDDJ+Y6PJK1gxYIjUbT82TyJcYWMgAcm052eDRrBy0QGo2m5zk+Y4hCf8jL8ekUulOmPWiBWGfE0wV+fHRGP0BtYjFboFTW19ppjk+nAHjlpaPEMwXm04UOj2htoAVinfFfvvI47/i7h/jx0dlOD2XNM5XIctNH7+Y3vvSTTg9lzTO5mAXgpr3DAJyY0W4mO9ACsY4olxU/OjIDwH1Hpjs8mrXPfYdnSOSK3PHUBLliqdPDWdNMJ3P43C4uGo0CMLmY6/CI1gZaINYRs6k8+ZKRI64Dec5zdj5T+f7ETKqDI1n7zCTyDEV8jPb5AZgyLQpNa2iBWEdMxI2HxuMSjk3rCctpJqomKSsFU+MMM8kcwxE/AyEfHpcwndQWhB1ogVhHWBPWT+0d5vRcWrs9HGYinmHfhggugaOTiU4PZ01jCIQPl0sYjviZ0i4mW9ACsY6oCMSeIUplxZm5zCrv0LTCxGKO7YMhtgwEOTWX7vRw1jSWBQEwEvVrC8ImtECsIybjWdwu4YotMeNn7ad1lMnFLBtjATb2BSruPY39lMuK2WSe4aghEBuifqYTWiDswDGBEJHPisiUiDy9wusiIp8QkaMi8qSIXHvB630iclZE/sqpMa43xuNZNkT9jPYFACMNU+MM2UKJuVSejX0BRvsCWowdJJ4pUCyrigUxHNECYRdOWhCfA15T4/XXAvvMr9uAT13w+h8D9zoysnXK5GKW0b4AGyqZHvohcgrr2o6aFsTkYk5vTnQApRQzpjtpOOIDjN3U8YzeKGcHjlVzVUrdKyI7a5zyJuDzynhqHhSRfhHZpJQaF5HrgFHgO8D1To1xvTGxmGXfhghRv4eg182UXmU5hhXv2dgXYDFTIFMosZgtEgt6OzyytYNSijd98sc8O74IULGMYyEvuWKZbKFEwOvu5BB7nk7GILYAZ6p+PgtsEREX8D+BD672C0TkNhE5ICIHpqf1xq/VmIgbFoSIsKHPrwXCQcbjRgLApligMnFpN5O9nJxN8+TZOIWSYZltGwwBVERYWxGt041B6vcCdyilzq52olLq00qp65VS14+MjLRhaL1LMlckmSuyMWZMVhuifr2ZyEEsMRiNBRgxg6czWpBt5bHT85Xvwz43m817uz9ouJoWdD2mlulkw6AxoLqrx1bz2IuAm0XkvUAE8IlIUin12x0Y45phssrlATAU9lcqYGrsZyKeI+RzE/V79IrWIU7PpRGB23/jZoYiPkQEMGIQoK+3HXRSIL4JvF9EvgzcAMSVUuPAO6wTRORXgOu1OLSOJRBWgLov6NEPkINYKa4iUpmwFvT1tpXTc2k29gXYv7nvvOOWIC+k850Y1prCMYEQkS8BLwGGReQs8GHAC6CU+hvgDuB1wFEgDbzLqbGsZ8plxVs//SAHz8WBqkBe0MtiptjJoa1Jjkwm+Nt7j/Pk2ALbBrRP3EnOzmfYOhB83nF9ve3DySymt63yugLet8o5n8NIl9U0yeGpBA+fXGrBaAlEX8BLplAiXyzj83RjKKo3+dQPjvG1x8YAuHHXEABBrxuf26V94jYzm8xx8cbo845rF5N96JlhjXPg5FIgbzjiI+I31gQx8yFazOqHyE6eNi01gBt2GwIhIvQFdW6+3cynCwyEfM87HvF7cLtEC7INdDIGoWkDp+fS+DwuvvqeF7Glf8kc7wssrbKsHaia1iiVFSdn07zjhu1ct2OAN1y1ufJaLOghntE+cbsolRXz6TxD4ecLhIgQ04JsC1og1jhn5tJsHQhy5db+845rP639nFvIkC+WuXxLjLdcu/W81/pDPn2tbSSeKaAUDCwjEGDc3zopoHW0i2mNYwTyQs873hc01gaL+iGyDWv39Ob+5QOnWiDsYy5l7CkZrCUQbc5iUkpxaCJBvlhu6+c6iRaINc5MMsdo9PkuJG1B2M98ypiQBpfxixsTlr7WdjGbNK71UHh592h/yNv2xc8/PXyaV//lvfz6Px5o6+c6iRaINc58Or+sGW7FILQFYR+WAFhZNNX0BTwksjqt2C7mTDEeCC9f26rdLqapxSwfveM5AO45NF1ZLPQ6WiDWMNlCiWyhvGyBuL6glcWkJy27mE9bk9Yyghz0ksgWKJd1RVc7mE2tYkG02aV358EJkrkif/TGywB4/MxC2z7bSbRArGGsFe1yqYB+jwuvW/Sq1kbm0nl8bhdh3/MriPYFvJQVpPL6ettBPRZEPNM+QT5wap5NsQA/f/1WAl4X331msi2f6zRaIIBcsURpDa7sKivaZVweIkI0YKxqNfawkCrQH/JWagJVYyUFaEG2h7lUnqjfg9+zfDnvWMiHUu273mPzGbYPhgj5PLzi0lF+cGiKbKHE1x8b6+mg9boXiHMLGV76sR/w0o//YM2VY7YEIraMQIDhF2+3i+meQ1Pc8Gff4+9/dKKtn9sO5tP5Za01gGig/RsTlVL8t395gts+f2DNJSPMpvIMRpa/1mC4mKB9SRjj8Wwle23/5j7G41k+fe9xfvOfH+cz9x1vyxicYN0LRH/Iy8ZYgNNzab715Hinh2Mr8RouJqDtFsRitsCHv3GQycUcf/m9w2vOH7+QLiwboIbqpID2CfLtT43zL4+e5bvPTPK1n6xaPb+nmE/lV0xxhaqCfW3YnFgqKyYWs2zuN8rY7B2JAEbZFYDbe3heWfcCEfJ5+Np7b2JD1M/TY/HV39BDzK8qEO3NrPnao2c5PZfm567bSiJb5PhMqm2f3Q7m0itPWu3ed6KU4i++c4jtgyGGwj4eOj63+pt6iNnU8ruoLdpZj2k6kaNUVmyKGRbEjqEwAJlCCb/HxTPji4wtZBwfhxOse4GwePHeYe56ZpJUbu34iC0XU61VbTstiMNTSfpDXv7TS/YA8N1nJtr22e1gIZ2nfxUXUyLXnus9lchxei7Nu27aya0Xj/Dwybk11RN7LpWrz4Jow94Ta/K3LAjrX4Dff/2lADxwbNbxcTiBFgiTn71uK8lckUdOzpHIFiiUejewZBHPFAh4XSv25Y0GPG11eYzNZ9jSH2TPSIRrtvdzl5npsRZcTUopFtKFZRMCwIj3QPtcTMemjWZQ+zZEuXH3EHOpPEemkiilel4olFLMpfIMrpDiCktxt3bshVhqL2tYENZiAODnrtvGUNjH/UdnHB+HE2iBMLl8cwyAg+cWedX/upd3/v3DHR5R68ynVg6aQvtjEGMLmUrBwCu2xDgymeRbT55j/4e/w4GTve0CSeSKFMtq9SB1m1xMx6YN993ukXCl7PhDx2f5D3//MO/5wqNtGYNTJHJFCiVV08VkWRDtuN7jCyuXWAn63Ny4Z4j7j832pDDrYn0msZCXobCPf3tsjPF41vzKVFYFvch8urCiywMMv3gqb6T4ul3PT820E6UUY/MZbt43DMBFo1GSuSJ/dfdRsoUy//zIGa7fOejoGJxkPrXyJjkAn8dFwOsi0SYX5tn5ND63i419AURgtM/PVx89yxNnjThbL9/bc2aZjVouJr/HTdDrbks9prGFDGGfu2IlAtz3Wy/Fb/ZZefHeYW5/cpzDk8ll+1d0M9qCqGL3SJijU0t9mu9+bqqDo2mdeCZfSfdbDiuzph1WxHy6QKZQqlgQu4aNQN5zEwkAHjjemz5ai6WEgNrXu10WxEQ8y2jMj8sliAh7RiIVcQC470hvujxgaRd1rTRXMGJv7YhBjMczbOoPnrf/ZdtgiA1mc66XX7IBgO8923ub57RAVGG1iLxp7xDbBoPc0+MCMZ8urLjTFJZ2oc634SE6ZwbyrBaR1eb4Ky4d5ex8htOzacfH4RRLCQG1LDZv2/ZBTMSzbOpbusY7TUHe0h9kJOrn3sPTbRmHE8zVKIpYTbsq6I7Hs2yKBVZ8fUNfgEs2RnsyUK0FogqrPtG12we4ac8wD52Y6+kd1rWyamAp/dUqnewkZ+cNgdjSb4hw9QP1gZfvBeD+Y727ql2osWvdoi/gadvGrYnFLKNV13iXmXr5oj1D3LxvmB8dnenZe3u1Ut8W7SrYd24he14zruW4cfcQB07N9dyuai0QVVjZPns3RLhh9yCJbJFnxxc7PKrmsLJqarmYrAdsLtW+VMAtpgVRnVl1xZYYI1F/T7uZrGtYKylgIORjvg3XWillWBBVAmFZEC/aPcQt+0ZYSBc4eK439/1UCvXV4WJy2qWXK5aYSeZWjedcv3OAbKHM4cmEo+OxGx2kruL9L9vLSNTPT1+5mZlkDhG465lJLt8S6/TQGia5SlYNLAlEO0oTj81nCHrd562wv/ObNzMc8SMi3LxvmO8enCSRLZyXJtgrLKTzuGTJCl2O/pCvLQuOhXSBXLHMaN+SQNy8b5j/+sqLeN0VmyoFA+87MvO8ToO9wFwyT8DrIuSrPX21owfHKdMtum2wtkBcPGoEpw9PJnpqPtEWRBURv4d3v3gXbpcw2hfgxl1D3HmwNzdzWSvVlbJqoMqCaEumR5otA+cH8i7Z2Ffph/3WF2wnmSvyox4Nns6n88SC3prZYAMhb1viPeNxI+2y2oIIeN38xsv3EfS5GY74uWxzX8/GIebS+RXLfFfTH/I5XmrjKTPwv9qkv3M4jNctHJ5M1jyv29ACUYMbdw9xaDLR1gJrdjGdNCaJ4RpmeNDrxu9xVYJ+TlK9B2I5rtneT9jn5sc9GoeYTxdqWmtgiHWmUCJbKDk6Fqvo5MYagdOb943w6Kl5kj1YOWBulTpMFrGgl2yh7Oj1fmosTtDrZo9Zf2klvG4Xe0YiHOkxF5MWiBpct2MApeDx073X/GM6YUz6I8u0G7UQEQbDvvYIxHymEn9YDq/bxQt3DXJ/D2Z6AMwkcqv6xC0BmXfYYrMsiI19KwvELRcNUywrHuzB692IQIBzm+WyhRJ3PTPJ1dv669pHdNFotJLW3StogajBVdtiuAQePTXf6aE0zHTSyPSoJRBgBU6dnbDS+SLz6UJdmR7Hp1PMJp3PqrKbmWSujmttphU7HKieiGcQqf1/f92OAYJeN/ce6T0302yydqE+C6cL9n3z8XOMLWR470v31HX+xRujjC1keqoHixaIGkQDXi7e2MdPTvegQCRyuGTllowWg2Gf4zEIK8V1aw0LAuCSTX0APeenBeN6W/GUlbBSjp3e3Xt0Osn2wRBe98qPt9/j5qptMZ7qwQrG9VoQ1j6JaYcWHI+emmco7OPFe4frOv8KM07RS+W/tUCswhVb+nh2vLfMQjBWtINh36qm72DYeQvCKqN+yca+muddNGr4cY9M9db1zhVLLGaLqwpEJWvM4UD1c+MJLl3lWgPsGYlwzCzg1ytk8iUyhdKqu6iByj4QpxqBjS9mn5d4UYub9w1z7fZ+/uaHxxwZjxNogViF3SMRZpK5nuvIVc+KFoxJa9ZhgXjizAJhn5u9G2oH8jb2BYgGPD2XKz6bXD3eA0suJicttky+xInZFJdsWr3mz56RCIvZouP//3Yya26Sq8fFZMVgJuLOWBBTi9nzUolXQ0R41WUbOTmb5sxcb1QN0AKxCrvNDUbHp3vL7TGdWN0nDkYMIpEtOlbevFRW/PjYLJdvia1qzYgIF41Ge87FNJ0wJqC6XUwOTsiPnppHKbhs8+q59rtHjHv72FTvXG8roWK1jDGAsN9D1O9xzIKYWMwy2rf6M1bN1duMfSc3/8U9PRFr0wKxCnvMVa9VPrlXmE7kGKnLgrACp85MWl9/bIyjU0neceOOus7ftyHSUxMWGO48qJ1SDEZF14jf4+iK/X9//zDDET8/tWdo1XOt1Mxe6uxX7y5qi5Go35EYRLZQYiFdqJkpthxXbl0S7v/ylSe63r3nmECIyGdFZEpEnl7hdRGRT4jIURF5UkSuNY9fLSIPiMhB8/gvOjXGetg+GGIg5OUbj491/X+mhVKqrqwaWHKLTCWcWc388PA0G/sCvOHKTXWdv3dDhNlUnifO9E5q8UydGWNgbF4751D7yXi6wIFT87zjhu2E/asXSdjSHyTkc/PIid7pxbFU6ru+lXtf0JlyG5ZV0oiLCYwWx8985NV84GV7uffwNP/5y4/bPjY7cdKC+BzwmhqvvxbYZ37dBnzKPJ4G3qmUusx8/1+KSMfqAXjdLn7tlt3cd2SmZ3KYE7kiuWK5zgnLyCyaiDtjho8tZNg1HK47kPfGqzYTDXj4P3cfcWQ8TjBjTlr1xHy2DgQd6098cDyOUkbdn3pwuYRffME2vvHEOU7N9oYVUankWkcMApyr6Dq5aCwKGhUIMETiPS/Zw9aBIN984lxX96t2TCCUUvcCtZYmbwI+rwweBPpFZJNS6rBS6oj5O84BU8CIU+OshzdfswWgZ8pAWD7xele0sNQ20W4mG/TTbugL8PorNvHIyfmesdimEzmifs+KrV2r2TIQdCxAaf2/N+L2eM+teygrxb89NubImOxmLp3H45LzmvPUwimBmKhjt3otQj4Pn3z7tcBSll830skYxBbgTNXPZ81jFUTkhYAPWDYvTERuE5EDInJgetq5DT/WKvtP73i2J0pSz9QZNLXO8biksvvWTpRSTC3mzis7XQ9Xb+snninwjr97qCdEYrpOdx4Yk/dituhI+YdGLBmL0b4AWweC/OX3jvBYD+z3WUjnGQj76rZIHbMg4s25mKqpJAl0cQJM1wapRWQT8I/Au5RSy6bYKKU+rZS6Xil1/ciIs0bG600f+ge/8kRbqp+2Qr27qMFwMwyEfY6Uf5hPF8iXyoxGG3uIrtthuEjuPzbL7U91/6aimTpTimGpeKITVUZnkjk8LqmUmKiXl11sdDx79z8c6PpdvnOpfM2eGxcSM2MQZZt7X0wsZgl4XXVbMssRDXjZEPVzvIsTYDopEGPAtqqft5rHEJE+4Hbg90z3U8f5nz9/FX/19ms4F89y3Z/c1dUparN19Oytxjk/bXNm+L7RKH/zS9dx0WiE9//TY3z+gZO2j81OppM5hqP1XetBB+sxWfWgXA32F/+d113Ke1+yh7lUntf85X1k8s4WE2yF1fqsX0gs6KWsIJm3tyjh5GLW7PfdWi/37YMhzs53756ITgrEN4F3mtlMNwJxpdS4iPiAf8OIT3y1g+M7j4DXzU9fuZnffd0llBV85cDZTg9pRSrtL+tcSTrtp200VxzgNZdv5FO/dB0Af/ytZ7p60pqpM6UYliwIJ6zQmWT9lkw1Aa+b33rNJdx2y27GFjJ88aFTto/NLuZT+VVbjVZjWVNxmy22yQY3ya3EtsEQZ+bWYZBaRL4EPABcLCJnReTdIvIeEXmPecodwHHgKPAZ4L3m8V8AbgF+RUQeN7+udmqcjXLbLXvYv6mPHx/t3ljEfMroTeCpUYunmn6HGqtMmQKxoUEXk8WekQif/g/XUSgpnpvozs5+9ZbZsHCyB8dMMt+UQFj87usu5ZKNUe7u4l7sq/VZvxCrgZPdC6DJxRwbbBCIrQNBxuMZxzaqtopjHeWUUm9b5XUFvG+Z418AvuDUuOxg32ikqyu8Gr0JGvPTHnKgvIWVCrihCQvCYp/ZievIVJJrtteXvtlOKoHhOoPUVoVRJ+oxzSRzXDS6eomNWlyxJcYPu7SRkNFGN1/XLmqLvqAxxdm9F2I2Wb/VWIuNsQBlZbiFm82IcpKuDVJ3MzuGwpxbyJArdqfbY97M9KiXPgddTINhH37P6umfK7F9MITP7eraTA8rY6zeyaLPbKdqdzBYKcVsMl93LGQldgyFmErkutKll6ijje6FWNd7MWtfDCJbKJHKl+rezV0Ly+Kb6dKYphaIJtg5FKKslspYdxtGpkf9N29/yEsiW6Rkc6ZHo8XMlsPtErYMBBnr0mtdKbNRpwXh97jwuISEjRMWwGKmSL5UbnlVu2PISL083YXF5Ky4TUOLHwcEuVLuo4FxrIQlEE6VJG8VLRBNYD1E3br7dKGO9pfVONV5q5liZsuxud+58hStYglEvZOFiBANeEjaLBDTyfr3vtRix1AIgJNdeG9bbrlG3KdRMw3VTkG2MhiHbHAxWYI+41Cpm1bRAtEEVmc0JzaX2YHRUKWxGATAggOBvEb3QCzH5phz5SlapTJpNbCajAa8truYZmwSiG0DhkB0YznqZiyIiCMC0VjBwFpYLkErltVtaIFoguGID5GlIGw3kS0YDVUazRUHey2Iclkxm8y1FKC22NwfZCqRI1/svkyPhXQBr1sI++qPs0T8HpI5ey2IJVdXa5NWf8iLz+OqlO3oJqz07UasY6/bRdDrdsTFNFxnwcBahHwe/B6X410Gm0ULRBN43C6GI/5KGmc3YT1E9W6SA2NFC/aushYyBcqqsXGsxJb+IEo5V1CwFRbSefpD9Zd+AMPtYWfQFBorr1ILEWFD1O9Ydd9WqBTqa0AgwLjeTriY6ulqVw8DIWcqGdiBFogmGe3zVzaCdRNLDVWa8dPat8qaS9nnp7WskOlk913v+XRjpR8AR2IQM8k8Lmlsdb0So30Bx5rstMJCuoBLlu7XeokGPCRydt7befweV0NWYy36Q17mUt1Z4kQLRJMMR/wVX2Q3sVAJ5DViQdjvp7V8qnZkegyZpnw3Xu9GSz+AGYOwccICqwe5f9WuffXQrRbEvGmtNVpKpC/otfXeNmJ8jVmNtRgI+bSLaa0xGPJVVuvdxFwTgbxoJVfc3lUW2ONiskz5brzeC01aEHanuRplNuxxeQyGfV1ZkLIZaw2M+9tOl54lEHYxEPZqF9Naw6kKqK2y0EQgL+J3MhXQDgvC+B1OtupslvkGU4rBDFJni7aWMp9O5usuOb4ag+a9bXcF1FaZTzV+rcESZBsXP2l7BaI/5HOk1I0daIFoksGwj3S+5Ehd/1awfJn9Day03C4h4rd3VWu5mOzwiQe8bkI+d9dZEFbph2ZcTMWyIluwLyurkZLjqzEQ8lFW9lqUdtBohQCLPpsttkY3oq7GQMjLQqbQlb1PtEA0yYCDZZtbYT6dJxrw4K2zULgvSjYAACAASURBVJ+F7asss2Bgo+NYicFw97n0UvkShZJqSIzB/qQAqwe5nS4m6D6XXksuJhtTuG13MYV8lMrK9sw2O9AC0STWRrRufIiauXnt9ovPpfK2uJcshsK+rnMxzTeRMQZLFUbtWqEnzR7ktlkQ4e5b/CilDBdTM/e230OuWLZlH02+WCaRLdruYgK6MlCtBaJJKhZEl6WnzaUad3mAFcizczNRruF89VoYFkR3ZdZYBQ4bvd5WFzK7VozNtBqthfX/1k2pl+l8iXyp3HQMAuyx2CoxPlstCOcq/LaKFogmcbKufysspAsMNmWG22tBzKcKtq6yBsN+5roszbWZnb1QlTVmk9uj0YKBq2H1W+imTKZmN8mBvRtBrXI0zbi6VqK/S93VoAWiaZzsDNYKzQbQ7K4PNGu3iyliuJi6KZDXTPE4gFjQZguisova5hhEF01YlS6JTS5+wB6BsKxGq0qsHVhZet1YsE8LRJNY7Ty7MQbRlJ/WRguiXFZNx0JWYjDsI1csk+6iPgULlUmrUReTvSWoLQvCjgY2AEGvG7/H1VWLH6s2VDOpvFEbr7dl9cXqbOdbD1algG7cnKgFokk8bhexYHdtcMkWSqTzpSYzPewTiMVsgVKDjV1Woxsza+abSCmGaheTPdd7OplHxJ5NiWDUY+q2rDGrMGYz/UWiNsZ8KhaEjQIR8nmIBjxdWSBRC0QLdNtDtNBE6WmLvoCXfKlsy74O65rYncUE3bVZbj6dJ+pvPKU44HXhdYttSQEzyRwDIV/dPcjrodsKyE0ljNpQzQTiYzZmjTlhQYBR3uRz95/kUz84ZuvvbRUtEC0w2GW7qZsNmoK9ftpWxrESAxULontWWQvpPP0N9N2wMJoG2RfzMTbJ2XetofsWP1OJHINhHz5P41NWLGRfUkDctPr6GiwYuBo5MwX3f3znOW5/ctzW390KWiBaYCDk66pUQMtn3Fogz46HyP5VVsWC6KJMprl0oelU3r6AxzYX02wqXyloaBcD3SYQizk2NJmlFfF5cAm29F1fzBYI+9y2WmsAH37DZZXv3/dPP+mafvdaIFpgMOztrkBeC8HKqN++VEBr4rNTILoxBjGXyjXt9zcqjNqzuIhnCpXUVLsYDHm76lpPJ7JsaLK/ucsl9AW9tghEPFOw3b0E8Mr9o5z46Ot4xw3bATgymbT9M5pBC0QLDIR9zKW7J/VyygzkNfMgOZEKaOeDFPF78LldXTVpzSXzDDa5crezaVA8U7A17RKMe3sxW6RQ6o4ufpMtWBBg3It2CYSdAepqRIT/+OJdADw7vujIZzSKFogWGI0GyBfL3PjR73dFD9+pRBa/x9WUf9TOVEAnMj2szJq/vfc4X330rG2/t1mUUi3t9eizsT7QogOrWssy+r1/e4qUze1RG6VcNmpNtSoQdlRMXXRQIAB2DoXxeVwcndIWRM9z+ZYYYKxu3vzX93fckphczDHaF2iqkYndFkTI57atUJ+FtZr94L88wbHpzj5A6XyJXLHcvIvJptIm2YIxDrsnre2DIQC+cuAsf3L7s7b+7kaZS+cpllXXWBBOuJgs3C5hS3+QswsZxz6jEbRAtMD1Owb44KsuYtdwmJlkju8/O9XR8Uwlsk0/RHYWkHNiRQvwW6+5uPL9T3/iR7b//kawNjU1f709tgVNjd9n7/W+ZvsAW/qDAHzp4dPcd2Ta1t/fCK24Ti36gvZYbIls0XZ33oVsHQhydl4LRM/jcgnvf9k+7vjAzQDcfajDAmFaEM0Q9XsQsSsV0H6fOMAvvmA7R/70tVyyMUqmUOroxqKJuJGXv7HJ690f8pEtlFvOVlmslH6wN+0yFvRy72+9lAO//wrcLuEbj5+z9fc3wqS5B2K0b+1bEABb+oOMzXfeZQ1aIGwh6HPzwl2DHJ5IdHQcU4lc013FXC4x/LRd/hB53S4+9Ib9ADzTwUDe5KIxaTW7qrVW/K1OWnEHMsYs3C5hOOLnRbuHONJBn/i0ZUFEm7cgLIFoxQ1cLJVJ5or0Be0V4wsZjviZS3VHRz8tEDZx0WiEQ5OJjsUhUrkiyVyxaQsCjPpStgTyskVHA3mXbuwD6KggT5gCsTHW3PWu7O5tUSAWHUgIuJBdw2FOTCc7dm9bu6hbaakaCxpd/Fqp5WXF55y2IAbCRkc/u/uWN4MWCJvYORQmkS3atvmpUVr1iYPZG9cOv3im4Ogqqz/kJer3cLaDZvhEPEvE76n0824Ua5JpVZCtGISTk9au4TCL2WLHypxMJXL0BTwEvO6mf0fMBovNifTt5ag0I+uCKg2OCYSIfFZEpkTk6RVeFxH5hIgcFZEnReTaqtd+WUSOmF+/7NQY7WRTzAjojXUo+2Bq0fLTtmBBhLy2dLVy2k8rImwdDHGmg4G8ycVsyz5xsMPF1AaBGAkDcHw65dhn1GJqMddSgBrsud6VhACHg9T9oe7ZFFqXQIjI/mWOvWSVt30OeE2N118L7DO/bgM+Zf7eQeDDwA3AC4EPi8hAPePsJJv7jRt4PN6ZSWvSsiBamLTscDFZflqnV1nbBoIdtSAmF7NNu5fAPoFYClI7d72tbCYr7tJu5tP5lrsT2mpB2NgsaDkGu6gFab0WxFdE5L+bq/6giPwf4KO13qCUuheYq3HKm4DPK4MHgX4R2QS8GrhLKTWnlJoH7qK20HQFm82HaDzemYeoYkG0EMjrD/lavinb5acdjvo7usKabCFjDJb6idhhQQS97qaK2NWLVUHV6jvRbuKZQlP1xaqxxYJwMCGgmoFesyAwVvPbgPuBR4BzwE0tfvYW4EzVz2fNYysd72qWelR35j91OpHD53G15PuPBb0sZouUWsiecKLj1nIMhnzMpwsdyfQol5VhQbSYlw92WBDOZ9X0B724XdIxgVhId4dAtOvetupq2ZEw0ir1CkQByABBIACcUEp1vEiLiNwmIgdE5MD0dOc28gD4PC7CPnfHGo9bPvFmdlFbWA9hLwTyBsI+SmXVkUyP2ZSxs7cVC8LtEqJ+T8uTQDvy8l0uo8xJpyrpzqfzDXftu5A+G7LG2nVvR/wePC7pqSD1IxgC8QLgZuBtIvIvLX72GIZVYrHVPLbS8eehlPq0Uup6pdT1IyMjLQ6ndexw0TTLVCLXUp442CMQlawah/20Vte8TjxElunfTPOaauzY3evUpsQLGY74O2JBWKVEWp2UrY2grd7bXrcQ8Dqb/Cki9Ifsq9XVCvX+pe9WSn1IKVVQSo0rpd4EfLPFz/4m8E4zrnEjEFdKjQN3Aq8SkQEzOP0q81jXMxC2Z6NZM8yl8pWeCc3SW2Z45/y0Sw2RWnd7tOxiyjpvQYB9e2QaxfrMVl1M1kbQVu/tWNDbkpVeLxG/h2SHiyQC1Ou8nBKR7Rcc+2GtN4jIl4CXAMMichYjM8kLoJT6G+AO4HXAUSANvMt8bU5E/hjDagH4iFKqVrC7a+hkm0Y7A3mtWEFtyxXvYMzHuj6tWkn9odYFIp4pcNFotKXfUQ+xoLcjBRIXMvZ1J2xV5CbjWUZatNLrJRLwkDTdp+Wy4lw8w9aBUFs+u5p6BeJ2QAGCEYPYBRwCLlvpDUqpt9X6hcrYlvm+FV77LPDZOsfWNcSC3o4V2TICea1aEMb7eyEGUWkg1AFBtuJMrU5asaC35bLO7YhBgH21jBqlYkHY8DfGWtwIei6eZXMLqc2NEPZ5SJgWxN/ce4y/+M4hfu66rXz8569qy+db1OViUkpdoZS60vx3H8b+hAecHVrv0SkLIlsokSmUWp4o7Cj/MJvME/K5Cfqa3/VaD5aLqRMWhF09t1utfWUF6Z0ss2HRF/TYUum3Ueyy1sAQmXgLz+e5hUwlnd1polUWxPeemQTgq4+ebfven6aiLUqpn2CkvmqqGDBdBq2kiTaDNaHb52JqfiIwguX29kdejrDPjUvsKU/eKAvpAn6Pq2URbHVVnmhDmQ2LWNBrS/XZRlmKQdjgYgo1L8ipXJF4psCm/ja5mMwYRLmseGZ8kVsvMpJw7j0805bPt6jLxSQi/1/Vjy7gWoy9EJoqYiEfShkPrh03dL3MV8zw1j7T53ER9LpbmrSmFrMtFVWrFxExHqIOpLkupPO2+MT7gl7yxTLZQqmpOkPt2rhV/RnxTIENUWetw2qsCd0OF1N/sPke8laFhM2x9lgQkYAhEDOpHNlCmZdfuoGfnJrnuYn2VjCu14KIVn35MWISb3JqUL3KkovGeHBnkzmyBedXXJYZ3qoFYf2OVgRiOtl6um29RAPeip8WaMu1BkOQ7brW0HzMp13xHrBnH0EzzKfy+NwuQja4LPtDvqY3gk7EjRTfVva+NELE7yWZLXJuwaiQsDkWZO9ohMOT7a1gXJcFoZT6I6cHshaoXmVNLma59WP3UCwpfvzbL3P0xlqwcaJo1S8+vZjjln3OWxBwvp/20/ce48/ueI7ff/2l/OrNux39XLssiOr7pZn7YyZlTFqtptvWQ/VYS2XFrR+7h4V0gYd+9+WEm6xoWw8Ti1lGY61tALWwBHkxU6jEsOpl1rzWw032IG+UaMBDvlTm1KxRIHFzf5A9I5G2d/araUGIyL+LyDdX+mrXIHsFq6tXPFPg/mMzZAtlimXF5+4/6ejnxm3KFQfTL95kDCKTL5HIFdviYgLDT2vtpP6H+08B8Ce3P0vK4fzx+XShUg6hFVqN+ZyeNQKWVv9oJ6kWiMfPLHB2PkMyV+Svf3DU0c+diGfZ1GePW8d6PppZAFmuqWZ7kDeKVUbeshg29wfY0h9kKpEjXywzncjxD/efdNxqXk36P+7op68xrEyLxWyBZ8cT+NwuLt4Y5cBJZ7dxWLnidgXyTs40lykxbUNPikaIBjzMJPMUS8YDc9XWGE+cjXP/sVleuX/Usc9dsKH0A7S+MfHcQgafx9UWQa6uHXVsyljVDoZ9PHBs1tHPnVjMctXWflt+lxWjM1yy4YbeO5fKI2LPM1YPlkAcnUri97iIBb1sGQiilBEP+ccHTvF3PzrBiZkUf/jGFXcbtMxqMYgTSqkfrvTl2Kh6lOoH/shkgj0bItywa5AnzsYdzf5YSBfwuISwHX7aoK8iOI1idf5qtXZ/vUQCXpK5ImfnM+RLZX72uq2IwNNjccc+UynFQrpgi1unVYEYj2fZFAu0ZWdv1LSOk9kip+ZSxIJefu66rTw9tnjevf3dgxN88h57rAqlVOVvtINYCxbEbCpfKVrYDiLm9T4ylWS0z/g/3mqm2I7NZ3jwhCHMn7v/JI+dnndsHKsJxNetb0TkXx0bxRrh/BhEjs2xAFdu6ydfLHNixrlmK/Np+0oAGE2DmpuwLAuiXX5ay8V0as6weC7d1MfOoTCHHGxFmsgVKZZVyxlj0LpATMRbqyjbCFbplMVskVOzabYPhrhmWz/5Upnnxpeu923/+Cgfu/OQLXWb5tMF8sVyS303qhlooc/CfDrfNvcSLGVtHZ9OVf6PtwwYAvHkWJyD5xb59Vt3E/K5+cqBs46NYzWBqJ5xnI38rQGCXjcel7CYKTCVyDES9bPN/E89O+fcDuvphH2ppbGQl5yZetko1q7mVovY1Us04CGZK5zn2tozEnFUjBdS9sV7ogFvSwXkxhcztk2eq+H3uPC6hUS2yNhChm2DwUqnOUuggcpO428+3noWvJVaapcF0d9CzGc22V6BqP6sUfPv3xQL4ve4+PNvP4dS8DNXb+FV+0e5/clzjsUiVhMItcL3mmUQMQqCzafzzKaMDWPbzADiGXMH5P3HZirBRbuw0wxf8tM2/hDNJe2rm1MPUb+HbKFc6XQ2FPGzdSDI2EIGpRTZQomHT8xhVHWxByubZcgGK8kq+d3s7t65ZL5tYlzZd5IrmIUh/WwzawNZu3vPzKWJmpbGXebu31aYMJtvbbRp70FfCwIxl2qzBVH1DI2aiz+fx8W7X7yrcvySjVF++srNLGaLPHnWGbfqagJxlYgsikgCuNL8flFEEiLS3h0bPYJR1CyFUjDSF2Ao7CPodXN2PsN0IsfbP/MQt3zsHltXuRPxrG0P0VKmR+OT1lw6T9TvcbS7WTWWn/bETIqA1+jHsbk/QDJXZDFb5P+/6zC/8LcP8L/uOmzbZ04uWtaKfX7xZiyIXLFEKl9qS4qrRTRguB8X0gUGwz7Cfg+DYR9n5jIksgVe/Zf3csjMujnSYo0pWOrOaNfix+0S+gKepq53PFOwxa1YL9UWarWV+DPXLPVOExGu3BYDnIu71XySlVJupVSfUiqqlPKY31s/9zkyoh4nGvRWUtM2WsGlgSBn5tJ8/9mlVdX//p49k1a2UGI2le8KM3w+lW84v7wVrNXqiZkUQ2EjV35Lv7GqHZvP8MNDRs74J+4+Wsknb5XpSiDeJpdek+U27CxBUS/RgIczpjvJsqCs3uCPnJwjnV9yc8wkc02nS1tMxLO4XWKrldTfZL20ZK5YCdS3A697aWquTvrYtyHCO27Yzv/9lRcYr0UDDEf8PDvuzHq9PUu9dUSsqqSwNWlvGwxxbDrJXc9MsqU/yM9dt5XvPzdFvth6U74pc0Vrly+6kunRjJ+2zWa4tXo+PJGoTNhWIO/RU3McmkzwSzcaVeq//fSELZ85uZjD7RKGwp0VCLsKBjZCNOCpxBusz906GOLMXJq7n5si5HPzrd94MX/25isAODrdWrLAeDzLhqjf1syhZq53qaxI50sVi7Xd7B5eSskVEf70zVfw0ks2VI7tGQlz3KG4mxYIm6nezWxVfhwK+zg2neL7z03xyv2jvPqyjSSyRR460XoOue2BvJBV8rsHMj1MgUjkimwxr/WOwRAugT/+1rMAvOumXVyxJcaXHj5tS//q8XiWkYh9k5aRVtyMtWaWHLdhw169RPxLix+rOdW2gRAnZ9N88aHT3LxvmMu3xLhx9yAAp1qMtU04EISPNdHFz2rcE3Fwx/hyWAX6Lttc21nz3pfu5TdetteRMWiBsJmY2UDe73FVVri/ctPOyus37Brk5n3DBL1u7jzY+qp2YtFeP21rLqZCW1e01e4Vq5nKQNjHW67dSr5URsRYff3Sjds5NZu2ZZV1ajbFjiH7di4323a0ExZEX9UK2kq+uOWiYQCUghftHgKoZNS1muo6YWPyhUUzFoQlEO10MQF85p3X88xHXr1q+vqtF43wkos31DynWbRA2IyVL169gemyzTGCZrXOSzf1EfC6efG+YX50pPXSveM2Z3qEfG58bldTq9q5VL6tQdOB8wRi6e//9VuMjOzNsSAiwlXbjJ24dgTyTs6m2TnU2C7cWlgTVqOZVp1yMQF43VKxjn9qz3BlZb1/sxEwjfg9+D0uZpPN916wNslttKnMhkVf0EM801gpFqveV8TfvnsbjKylkK8zbi0LLRA2Y7mYLgzW/tOv3cA7btheqZtz6cYop+fSLe+wnohnifo9tpm/IkKsic1yxVKZTKHUluY1FtUr2mqB2Lshwu+89hK+8KtGy5K9IxG8buG5FjfQZfIlZpI5tttoQcSCXgolRabBPHa7ejU3gnVPb4wFznOxfenXbuQVl27gyq2GQIgYgeXpFiyIRK5IOl+y3YKwLLZGBDmZM651p2IQnWT9/cUOY1W2vPiCPsHXbB/gmu0DlZ/3bIhQVoaftpWewuNx+/20/UFvwzGIVM6Y4Nrpp/VUZXrsGYlUvhcRfv3WPeedt3UgVMnAaRbLnWfn7uXq3dSNrBbjmQIBr6upPhLNsqUSUzs/QH/F1hh/98svOO/YcMTXkgUxZrbudSIGkS+VyRXLdV+7RLYzMYhuYP39xQ7zxqs388jJOT7w8n01z7MsiTNzrQmEsQfCZoFowoJIWKusNj9Ev/yiHcylCxWf+EpsHwxxaq61GMTSxi37BWIhXWBTA27CZK7Y9mt9y0UjDEd8/OdX1L63wdi0aG1gbIb/+d1DeFzC1dvsKdRnUS3I9QqEtfhpdwyiG1h/f7HD9AW8/O+3XrPqeVZud6uBvPF4los3Ni8wy9EX8FZiG/VSyfRo80P0R2+6vK7ztg4EeeLsQkufZRUjtLO3R7NNg1K5oqN9GJZjtC/Agd9/ZV3nDkd8PHOuudz8p8fifO/ZKd7/0r2rCn+jNNODI9mhxU83oGMQHWIp06N5M7xQKjOdzNkWoLboC3orFkG9JLvcDB+O+FlIFyiWmt97YlkQozZtkoPmC/Yls+23IBphKOJnNpVrqszJT8zqpO8w97DYSaXoYAPXu+JiWocWhBaIDhHwuon4PZVCc80wlcihlH0prhbRwFIjnnpJdMiCqBdr5+98C7t7JxazhH3uyg5uO2haIDpgQTTCUNhHoaTOawlbL8+OJ+gPeR2pVNvM9bas43CHM4o6gRaIDjIc8THbZBN1gPEFZwJ5lkA0svpLdWgzUb1YG/isYnvNcGbO/oSAZns9p/LdbUFY6bfzTdzfz00scsnGqCN9LqwFTLIB4crkSwS8rrb1gugmtEB0kKGIn5kmLQilFH933wlcYtRnsZNowEup3FjqZbe7mKzMm7kmXXqnZ9M8eHyWq7cNrH5yA0T9HkQa35iYypW62oKwBHmuQYEolxWHJhJcstGZUm/W/dmQQBRKlX1M6w0tEB1kOOJrOkj9+JkFvnNwgv/0kj2VXcR2seSnrf8h6lSQul4sF1OzFtvHvnuIYrnMr92ya/WTG8DlkqZ390b83TtpWXsmGi2MN7aQIZ0v2Z54YWEJRCN9yzN5LRCaDjAc8TctEA8cN+o43XbznlXObBwrnS+RbTyQ161+2oqLqYnrXS4rvv/sJG++ZqsjK9tmBCLVgTTXRhgMWde7MYE4Z7pNqzc+2knI50ZkyeKth2wDeybWGlogOshg2Ec8U6DURBG5QxMJNscCleqrdmIJxGIDD1EyVyTkc3etn3Yg5EOkcZcHLK1qrZ3CdtOoQFjVRbvZxTQa8+N1S8P1r5Z6QDgjECJCxOdpKHhuxCC0QGjazGDYR1k1HqAEQyCcMsOtLJ1GLIh0vruzatwuYSDUXFKA1d/D7liPRaMCkcp3d7wHwO9xs3dDlO88PdFQFd1zZnXizf3OtVIN+z0NuZiyhRJBnxYITZupBPIa9NMWSmWOTSe52KFAnlWRthELohf8tEPh5so/WN3R9rWw470WfY0KhJV22cUCAfCayzZyYibFR771TN3vGV/IEgt6HS1SFwl4dJC6TrRAdBArFbBRt8eJmRSFkuKSLrIgeuEhGgz7mnIxHZ5MMNrnP6/Xh530r1GB+MDL9/L6Kzfx+QdO1h37GY9nbN/XcyFhv4dkrv4MPSvNdT3i6F8tIq8RkUMiclREfnuZ13eIyPdF5EkR+YGIbK167S9E5KCIPCsinxAnkqI7jGVBNJrqetpsxLJz2L6y09UsBakbWWWVCXS5GT4c8TPTxD6IY9Op84oB2k2jJb8zeWM3eKjLBVlE+KUbdlBW8FSdpdbHHegBcSFRv4dkA4ufbEHHIGxHRNzAJ4HXAvuBt4nI/gtO+zjweaXUlcBHgI+a7/0p4CbgSuBy4AXArU6NtVNYjWcaDeQ57acNeo1gcyMWRDZfItjlq6yhiK+pfScT8UylkqkTxILGvpNUvr5VbdYsEd8Lk9Z+sxvawTrrMo3Hs2xy8FoDhP3uSgG+esj2gHXsFE4+0S8Ejiqljiul8sCXgTddcM5+4G7z+3uqXldAAPABfsALTDo41o4QDXjZHAvw4PHGWo+eW8jic7sYtqkv8oWISMPlNrLF7l9l7RgKs5gtNlT2u1gqM53IObqqtfaO1Bs4zZhCEvR1tyCDIX6bYgGOTScrx1K5Ip+85yhPXlA8MVsoMZfKs9lhCyLi9zYeg+hy69gpnLzDtgBnqn4+ax6r5gngLeb3bwaiIjKklHoAQzDGza87lVLPXvgBInKbiBwQkQPT09O2/wHt4KWXbOC+IzP87Q+P1f2ecwsZNvUHcDmYUtqoQPRCkPpV+0cBuPkv7mGhzsSA6WSOsrKvY99yhMzJp26BKPSOBQGwazjMCdNKvvu5SS778J187M5D/OY/P37eeXZ3R1yJiN+95uJrTtHpJcgHgVtF5DEMF9IYUBKRvcClwFYMUXmZiNx84ZuVUp9WSl2vlLp+ZGSkneO2jT/46f1sigX4xPeP1N1d7txChs0OP0RhX2OpgL3wEG0bDPHfXn0xAF986HRd71matJyx1oBKxk66XhdTDwvE3//oBEGvm1svGuH4dOq8nhGTDjRkWo5IwEMqX6or5lMuK7KFMv4eudZ246RAjAHbqn7eah6roJQ6p5R6i1LqGuD3zGMLGNbEg0qppFIqCXwbeJGDY+0YAa+b33/9flL5Eocnkqu/ActP63ymR70TFpiBvB4ww9/30r3s2xDhwMm5us6vNAmyuTdyNeEmBaLbBdli13CYhXSBu56Z5MdHZ/m567ZWhPqBY0vu1YpAOCjGYNzbJXPiX41c0TinV6613TgpEI8A+0Rkl4j4gLcC36w+QUSGRcQaw+8AnzW/P41hWXhExIthXTzPxbRWuGKLsUO3nkyPclkxuZh1fJUV8rkrG7LqIVso98xDdNW2fp4aqy9oOlHZ2evc9bb82/VebysG0SsWxO4RI9vuQ994GoB3v3gX+zf10R/ynicQU4tGAsEGh+/taAMF+5bEuNPOls7g2F+tlCoC7wfuxJjcv6KUOigiHxGRN5qnvQQ4JCKHgVHgT83jXwWOAU9hxCmeUEr9u1Nj7TRbB4IEvC6OVwXyiqUyX374NEenzrcqZlN5imVle9npCwn7PKTrzPRQyqj82iu54rtHwswkc+dNEIcmEtx5cOJ5504sZvF5XJXOb04QNovuZeq0IDKF3lrV7ho2UoTH41lef8Umdg6HcbmEK7bEODi+tCiaXMwS9LorE7hThBsQCCves16D1I7+Tyil7gDuuODYh6q+/yqGGFz4vhLw606OrZtwuYSdQ2FOqTrQrQAAGrlJREFUzhp+2u8/O8n7/uknZAtlrt7Wz9ffd1PlXMsMt7Pt5XKE/PVbEIWSolRWvTNhDRkr2pMzKS7b3McffONpvvCgEZP44q/ewE17hyvnTsQNa83JbTiWi6nemI+1qvV7ekOQqwvvXbppaXPnJRuj/MMDpyiWynjcLiYTOUb7/I5ea1hKCqhHkHstIcBueuMOWwfsHApX9kP81T1HEYR9GyI8NRZnsSrjol0CEfbVH4PotYdohyUQsykeOTnPFx48zdXb+gG478jMeee2w51nrU4biUH4PS5Hs9jsxOtemmZeuGuo8v1V2/rJF8v87b3HUUoxtZhlQ9TZaw1L92k9/U56zZ1nN1oguoSdw2HOzKV56Pgsj51e4N0v3sWf/MzllMrqPD/tRJsyPUJ+d90r2lyPmeE7h40Niqdm03z63uO4XcI/vOuFXL2tn0dPnR+8nlzMssHGHtTL0WiQuhfz8m/eN8xQ2McLdi41XHrpxRsYDPv42J2H+PbTE0wlco5fa1hyzdVjQfRaQoDdaIHoEnYNhyiUFH9gBvJ+8QXbuGb7ACGf+/xMj3gWlxjNhpwk7POQK5YpllbP9KhYEJ7eeIhCPg+jfX7uPTzN956d5Fd+aiexkJfrdwzwxNk4eTNzRSnF5GLOcTEOeF2IGBVx66EXd/Z+7l0v5MHfffl57qOw38Odv3kLYZ+bu56ZZHIx67hlDEtpxfVYEFamU68Jsl1ogegSdppuj8OTSX7h+q1sGwzh87i4ZGOUZ8eXMm4mFrOMRP143M7+11U2bzXgp+2lh2jHUJiHThjWws9ea5QAu27HAPlimYPnjMBpIlckUyg5PmmJiLnvpP4gda+5PNwuOc/VZDES9fOiPcP8+OgM6XyJ0XZYEOYO9HoEudcWP3ajBaJL2FVVeG//pqUy3hdv7OO5iUSlpv5EG1a0sJTpUddDlO89M3yb2aY15HNz0aiRZXOF2RDIqhs0Zbrz2uL28LnrtiDWWgObq7bGmDJrZLXDggiaFkS2nhhEJb62PqfK9flXdyEj0aVJ6OrtS37a63cMEM8U+MJDpwBj0nI6Txyqyz+szUwPa6V6+eZYxRrb0h+kP+Tlnx85QzpfZCLevkkr7HM3FKReS3n5V5oJAkBbgtSNxCAypmiHury0ulOsnbusx6n2zV5V1dryFftHiQY8fOgbB/nh4Wkm2pBVA9WB09VXtTnTT9tLq6xLTSvtXTftrBwTEd71U7t4aizOu/7vI20r/QCGX7yRGEQvifFqXLll6X63EgicxFr8pOuwICzR7vbS6k6xPmWxS3n4916Ox+U6TyxiQS/f/s838+L/cQ+3P3mOhXTB8U1yYGQxQWMWRC/FIF5/xSau3THwvDLeH3j5Xk7PpfnXn5ytiEg7XEwhX/0lqDOFEn0ONS/qBAPhpYSL0TZYEH6PkRSQrcOCSOd77962k95Z8q0DNkQDlSZC1WwdCHHDrkG+9+wU0C6Xx9qOQbhcsmyPBxHhl27cDsAdT40TDXgcbX9pEfJ76lrRQm+mua7GjbsHGQh527K3Q0QIeutz6aXzRVzSO5sS7UZbED3ClVtjlayb9gSpm8hi6iGBqMVlm2P4PC6mEjn2bXCuk1w1YZ+b8YVMXedm8qU15/L44q/eWHdHPTsIet11pbmm8yXCPo/ju7u7lfUpiz3IFVuXAnlOV7uEqhLUDRQ0WyslkX0eF5ea/b7bYa2BFYNYvxaE2yWOp25XE/C66wxSr71r3QhaIHqE6kBeO11M9VgQlUDeGnqQ9m82rnc74g9gXLtG0lzXirXWKUK++i2ItXRfN4oWiB7B6l8NRqtSp1kKUq8+aaXyRXwe17IboXqVXWY2TaRN6Y1GccTVJ6xSWZErltf1qtYOgg0JxPr1xK/fv7zHEBGu3BqjXW5ar9uF1y31BfJyJcJrbMJ62SWjfOPxc7z9hu1t+bywz0O+WKZQKtcU2vVeG8guGglSr2cLQgtED/H19960+kk2EvJ5KhuFapHKF9fcKmvvhgi3f+B5XW4dI1RV0TUWXFkgejGluBsJ+tzMpVbvS57Ol4gG1ta93QhrxyewDnC5pK0lno2ucnVaEH49YbVCpYDcKte7F1OKu5FQnTvXE9mCFgiNZjlCvvoyPdaiBdFultKKa1ts2oKwh3qzmBLZIlH/2tmU2ChaIDQrUm/5h0xeWxCtspRWrC2IdlDvPohEtqgtCI1mOep1MaXWeaaHHSzFIGoLspVVFl6nxePsoh7ruFAqkymU2pI12K1ogdCsSL0upnS+uOaymNpNdZC6FglTINqVfrtWsSyIWru3k1njWmsLQqNZhpDPs6pPHIyCfuu1HLJdhCobE2tfbz1p2cNST4iVOyYm9LXWAqFZGW1BtA/LgljtelsCol1MrWH106gVh1jMFgAtEBrNstSTClguq3W/29QOrKyk1QKn1qpWu5haI1RHteKFtCEQ/SFn+793M1ogNCsS8q+exWRNaOt5t6kd1NvlLJkr4nXLui0/bRcB836t1XZ0NmV0FBxapgT/ekHfZZoVCXndFEqKQml1P21kHZvhdmAJxGoWWypXJOJfv+Wn7SJUx/WeN3daL9ejZb2gBUKzIsE6Mmvm08ZDNLCOzXA7cLkMq6DWihYMQdbxh9ap596eS+UR0S4mjWZZrImolpvJEoj+0PrNFbeLYB0xn4V0XouxDdQT85lN5ekPenG3sbxNt6EFQrMi9eTmW4E8PWm1TqiO3b1z6YIWYxuoJ2tsLpVf1+4l0AKhqUE95R+0i8k+AnWkFS+k9aRlByGvZR3XFoihcHsaRnUrWiA0K1JP+YelVEC9qm2VerqczaW0i8kOKi6mGvf2dDLHUGR9X2tHBUJEXiMih0TkqIj89jKv7xCR74vIkyLyAxHZWvXadhH5rog8KyLPiMhOJ8eqeT4Vgagxac2n8oR8bgK6eFzLBFepMFoolUlki1ogbGC1GES2UOLUbJo9I5F2DqvrcEwgRMQNfBJ4LbAfeJuI7L/gtI8Dn1dKXQl8BPho1WufBz6mlLoUeCEw5dRYNctTj4tpTgdNbSPo89QU40q8J6yttVZZLa346FSSUllxyaZoO4fVdThpQbwQOKqUOq6UygNfBt50wTn7gbvN7++xXjeFxKOUugtAKZVUSqUdHKtmGepxMc2n8tq9ZBNBr4tszYQAHe+xC7eZVrySxXZ8JgUYnQXXM04KxBbgTNXPZ81j1TwBvMX8/s1AVESGgIuABRH5mog8JiIfMy2S8xCR20TkgIgcmJ6eduBPWN+E6kgFHI9n2RQLtmtIa5qQz0O6sLIYWy0ytUDYQ61SMmfmjPXo9sFQO4fUdXQ6SP1B4FYReQy4FRgDShi9sm82X38BsBv4lQvfrJT6tFLqeqXU9SMjI20b9HqhUmF0BReTUoqz8xm2DmiBsAOjy9nKu9bPzGcA2BgLtGtIa5pa1YpPzaYYjvjXfY0xJwViDNhW9fNW81gFpdQ5pdRblFLXAL9nHlvAsDYeN91TReDrwLUOjlWzDAGvC5GVMz3imQLJXFELhE0Y1XNXtiCOTCbwuV3sGFrfq1q76At6Wcwsf71Pz6X1dcZZgXgE2Cciu0TEB7wV+Gb1CSIyLCLWGH4H+GzVe/tFxDILXgY84+BYNcsgIoS8K3eVO3huEdB+WrtYrYnNM+OLXLQxgtfdacN/bRALeljMFJZ97fRset27l8BBgTBX/u8H7gSeBb6ilDooIh8RkTeap70EOCQih4FR4E/N95Yw3EvfF5GnAAE+49RYNSsT9HlqZnoA7N/c184hrVmCPjdlBbni8m6mZ8cTXLpRX2u76At4iS8jELliifHFLNu0QOCog00pdQdwxwXHPlT1/VeBr67w3ruAK50cn2Z1wv6V3R5jCxl8HhfD63y3qV1YqZfZQul5+0rS+SIzyRw7h8OdGNqaJBb0VpoCVTM2n0Ep2KEFouNBak2XE6zhYjo7n2ZrfxDXOi5mZie1al+dNQPUOt5jH7Hg8hbEaSuDSccgtEBoahP2e1bMFR+bz7BFT1i2UWt379l5Y9LSbg/76At6SedLzyuxPh7PArC5X9/bWiA0NTFyxVd2MW3RD5FtVHb3LpNWfGbOsCC2DWiBsIuLRo1d0k+cWTjv+OSiIRAbotp1qgVCU5Ogd/nNRNlCiZlkXguEjUQDxo70xDJ+8bPzafweF8PrvHicnbxozxAAD5+YO+/45GKO4YhPZ4uhBUKzCmH/8llMYwvGila7mOwjarZtXcw+32Ibj2fZ3B/UrUZtJBb0snUgyKHJxHnHpxazbIjqzYigBUKzCsEVXExjZtBUWxD20WdaEMtl1kzEs4z2aZeH3Vw8GuXQxPkCMbGor7WFFghNTaJ+D4vZIsXS+bn52oKwn42xAD63q7K/pBpd88oZrts5wJGpJOfM+xkMF5MuZ2KgBUJTkyu2xsgXyzw1Fj/v+Nh8BrdL2NinHyS78Hlc7N/cx+MXBE3LZcVUIsuovta287JLNgDw0IlZwOi5MZvKaReTiRYITU1u2GUE8g6cnD/v+NhCho19ATw6kGcrV2/r5+mx+HkW22wqT6Gk2KRXtbazeziC2yUcmzLKe88kcyiFFmMT/XRrajIS9bOxL8DBc+dbEGfn0zr+4ADXbO8nnS9xeHLJzTRh5uXrSct+fB4X2wdDHJ8xrre1B2JjTMcgQAuEpg4u39JXKcxnoTfJOcPV2/oBznMznYvrhAAn2TMSrlgQVixiS7/ebwJaIDR1cPmWGMemk8wkcwAUS2UmFrN6wnKA7YMhBsM+Hj+z5NKzJq3N/dqCcILdIxFOzKYolRX3PGc0HtukrzWgBUJTBzfvG6Gs4PHTxqp2KpGjrHQpAicQES7bfL7FNh7P4ve4GAzrTXJOsGckTL5Y5g+/eZB//clZYCnleL2jBUKzKla/h2PT5/tpddDUGfZv6uPIZJKCGageW8joTXIOsnvEuL//8cFTgOFS1RhogdCsSizoZTjirwjERCWQpwXCCS7d1Ee+VF4S5IWMFmMHuWRjtPL9b75iH1++7UUdHE13oQVCUxd7RsIcnzYCeROL2oJwkks2GROWtcPXKrOhcYZowMtn3nk9N+wa5O0v3E7Ev777UFejr4SmLvZsiPDtp8aJpwt86gfH8LqFWFD7aZ1g13AYlxgd+8YWMoZAaDF2lFfuH+WV+0c7PYyuQwuEpi52D4eZTxf4r//yRCWbSfvEncHvcbNjKMzDJ+b4P3cfBWD0/7V357FWlGccx78/uSggawVlU24VEVFZFBVcCi2N4h6rTbGKYmsstq5Va9VU0zZNbTBWrVbiFheI0SoqVSwgiDStKNuVRdTgUiruG4oaI/L0j/c99w7jnMVyFu89zye5Yc7MnLnPy7kzz7zvnHnGE4SrAR9iciXZLV6ofnzNWwBcMn5wLcNp8wbu2JmnE2Wod/LSD64GPEG4kuw3oEfz9M0T9+OssbvVMJq2b/SuocRJh/bbcP/k0Yzbc8caR+TqkQ8xuZJ07dCe+ReO4eGm1xm7hx+sKm3SQY1s/HwTw3fuzsjGb9U6HFenZGa1jqEsRo4caUuWLKl1GM4516pIWmpmI7OW+RCTc865TJ4gnHPOZfIE4ZxzLpMnCOecc5k8QTjnnMvkCcI551wmTxDOOecyeYJwzjmXqc3cKCfpHeA/Zd5sT+DdMm/zm87bXB+8zfWhlDYPMLNeWQvaTIKoBElL8t1h2FZ5m+uDt7k+bG2bfYjJOedcJk8QzjnnMnmCKOzmWgdQA97m+uBtrg9b1Wa/BuGccy6T9yCcc85l8gThnHMuU90nCEl7SGpK/Hwk6fzUOmMlbUisc0Wt4i2HUtoc1xsbl6+W9GQtYi2XEj/nixPLV0n6UlKrfZxbiW3uJunvkp6Nn/PptYq3HEpscw9JD0paIekZSXvXKt5ykXRB/PxWSbpHUofU8u0k3StpraSnJTWWtGEz85/4A7QD3iTcOJKcPxZ4pNbxVbnN3YHngF3i6x1rHWul25xa5xhgfq1jrcLnfBnwpzjdC3gf2LbW8Va4zVOAK+P0YGBerWPdynb2A14BOsbX9wGTUuv8HJgapycA95ay7brvQaSMA14ys3Lfkf1Nlq/NPwZmmNk6ADN7u+qRVU4pn/NJwD1Viqca8rXZgC6SBHQmJIhN1Q6uQvK1eQgwH8DMngcaJe1U7eDKrAHoKKkB6AS8nlp+HHBnnL4fGBc/84I8QWxpAvkPCqNjN/wxSXtVM6gKy9fmQUAPSQskLZV0apXjqqRCnzOSOgHjgQeqFlHl5WvzDcCehAPKSuA8M9tczcAqKF+bnwV+ACDpAGAA0L+KcZWVma0HrgbWAW8AG8xsTmq1fsB/4/qbgA3ADsW27QkikrQtcCzwt4zFywjd1GHAX4CHqhlbpRRpcwOwH3AUcDjwG0mDqhheRRRpc84xwL/M7P3qRFVZRdp8ONAE9AWGAzdI6lrF8CqiSJuvArpLagLOAZYDX1YxvLKS1IPQQ/g24XPcXtIp5di2J4gWRwDLzOyt9AIz+8jMNsbpWUB7ST2rHWAF5G0z8Bow28w+MbN3gYXAsKpGVxmF2pxTsIfRChVq8+mEoUQzs7WEsezBVY2uMortz6eb2XDgVMK1l5erHWAZfR94xczeMbMvgBnAQal11gM7A8RhqG7Ae8U27AmiRd4xZ0m9c+N1sUu6DSX857YChcbZHwYOkdQQh1wOBNZULbLKKXhtQVI3YAyh/W1FoTavI4zVE8fh96B1HyxzCu3P3WMPA+AMYKGZfVS1yMpvHTBKUqd4nBrHV/fVmcBpcfpEwhcwit4l7XdSA5K2J/wn72pmG+K8yQBmNlXS2cBZhIt3nwG/NLN/1yrecijW5vj6YsIZ5mbgVjO7tkbhlkWJbZ4EjDezCbWKs5xK+NvuC9wB9AEEXGVm02oUblmU0ObRhAu2BqwGfmpmH9Qq3nKQ9FvgR4Rj1HJC4rscWGJmM+PXXu8GRhC+iDDBzIqeCHiCcM45l8mHmJxzzmXyBOGccy6TJwjnnHOZPEE455zL5AnCOedKIGmKpOdjkb8HJXUvsG47ScslPZKx7HpJGxOvB0iaF7e7QFL/xLJdJM2RtEbSc7kie5KmS3ohFue7XVL7EuIfKumpWNRvZbqgXxZPEK7VkbRDolrnm5LWJ15X5OvHkkZIui1OD4472ueSLkqss7OkJ+KOvFrSeXm2pXiQWBsPCvtmrNMoaVUl2pKPpMfjXbl1T6GS8R2p2XOBvc1sKPAicGmBTZxHxn1DkkYC6f/jq4G74nZ/B/wxsewuYIqZ7QkcAORqok0n3NC4D9CR8LXWQu1pAKYBk81sL0IB0i8KvQc8QbhWyMzeM7Ph8U7YqcCfc6/NLH0HablcBlwfp98HziXs2EmbgAvNbAgwCviFpCEZ2zoC2D3+nAncVJGIM8QDRT53E6p+ugxmNifWMQJYRJ76TbEHcBRwa2p+O0Il2V+l3tJcPBB4glA2g/i302Bmc+Pv32hmn8bpWfHudwOeycUiafvYo3gm9mCOi9s9DFhhZs/G979nZkXLi3iCcG1KrusezwCflPSwpJclXSXp5LjjrJS0W1yvl6QHJC2OPwdnbLMLMDSxc71tZotJnYGZ2RtmtixOf0w4g+yXEeZxhDNGM7NFhLpAfTLWayfpltgbmSOpY4xnuKRFiaGOHnH+gniGiqSekl6N05MkzZQ0H5gnqY+khWp57sWh8ffNJNyB7Ir7CfBYnmXXEpJAuujh2cBMM3sjNb+5eCBwPKG67g6EgpkfSpoRD/ZTYpJpFoeWJgL/iLMuJ9wlfQDwXWBKvHFwEGCSZktaJimdpDJ5gnBt2TBgMqFa6URgUNxxbiUUaQO4jtAD2R84gdRZXzQS+FrDPXGseATwdMbi5sqa0WtkJ5LdgRvjkMCHMT4Iww6XxCGJlcCVJYS0L3CimY0hlHKfHXtgwwjF+oh3E28XD051SeFhOk2Ev4NjE0OXhyfWuZzQW5ye8f6jgbfNbGlqfl/gh4Rin2kXAWMkLSeUeVlPKB7YABwal+8P7ApMSr33r4RSIf+Mrw8Dfh3bsADoAOwSt3UIcHL893hJ44r9fxTqbjrX2i3Ona1JegnIlUBeSTi7glDobIhaSuN3ldQ5V5wx6gO8U+ovldSZUCr8/K2s8fOKmTXF6aWE5xZ0A7qbWe4Jf3dSuDJtztxEddrFQO7C5kOJ3wFhjLsvbaPW2NdmZgdC6IESHrozKblcoRTL0cC4PLWMDiYkliMJB+eukqYR6kINBNbGv7VOktaa2UAze52W8uOdgRPM7ENJrwFNuZIYkh4iDF3mroVdSSg0+LNkiPH9L6TiHkZIJO/G17MIJw3zCv1/eA/CtWWfJ6Y3J15vpuXkaBtgVOIaRr9UcoBQf6voNz6gucv/ADDdzGbkWa25smbUP84rFH/ujLKQTbTs0+l4P8lNmNlC4Dvxd96hLZ/10YHQXpciaTxh6OjY3LWANDO71Mz6m1kjoSrwfDM7xcweNbPeZtYYl31qZgPjdntKyn1ulwK3x+nFhOHHXvH19whPeUTSGYRS7SfZls/vmA2cIzUXFx2RmL+PQkG/BkJP5blibfYE4erdHFqGm5A0PGOdNYSzv4LiTnkbsMbMrimw6kzg1PhtplGEB7ykx6UzxeJzHySuG0wEcr2JVwnP8IBQsTNfnAOAt8zsFsJQyr6J+HvH7bivugHoAsyNw065Ao994xn5/2ss8IKkF4GdgD8AxIvIFxGuG60k9A5uie+ZGtd9KsZyRZz/e6A9sELS6vg6N3x4DSHpNBFKoT9aLDAfYnL17lzgRkkrCPvDQsJ1i2Zm9rykbpK6mNnHknoDS4CuwGZJ5xO+iTKUcMBeGceAAS4zs1nasmrsLOBIYC3wKaFi7tdxGjBVoQz7y4n3Xw3cJ+lMoNDOPxa4WNIXwEbCMxEgJJdFiW/q1C0zW0AYw0/OyzxJiENER5ayjcSyzonp+wmPAc1aby7h7yo9P/PYbWafseWQU3LZNMJXXUvm1VydK4GkC4CPzSzrInabIOk6wrdsCo5Lu/rhQ0zOleYmtrwm0Bat8uTgkrwH4ZxzLpP3IJxzzmXyBOGccy6TJwjnnHOZPEE455zL5AnCOedcpv8Ba/oymHS0MIAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# import packages to run\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import numpy as np\n", "\n", "# find a way to include or fetch this data\n", "filename = \"../../../CometGPs/data/paper_plots/3200/3200_lc_49627_to_49787.txt\"\n", "\n", "data = pd.read_csv(filename, delim_whitespace=True, header=None)\n", "\n", "time = data[0]\n", "flux = data[1]\n", "\n", "x = 1440 # 12 hours, every index is 30 seconds\n", "\n", "plt.plot(time[0: x], flux[0: x])\n", "plt.xlabel(\"Time (%.1f hours)\"%(x*0.5/60))\n", "plt.ylabel(\"Flux\")\n", "plt.title(\"Phaethon 3200\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Currently (as of 2020) we are seeing a slew of new survey telescopes being developed (e.g. LSST), intended for all-sky coverage. While these telescopes are great for developing time-series data, their expansive field makes it so that the same objects are only photographed every couple of nights. Rather than looking like the asteroid above, it would instead look sparse, like the following plot, if it was only observed every 12 hours for a month." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAezklEQVR4nO3df7RcZX3v8feHkGCUQICkFE6AgEDKb4JHfjSikVoDVCHEtgu0Cv5YXC/irb2LVLjWK2K5oNBra7V4QSNGvSBKSFHwBkpESivICQFCwGBEkJygCWKASBQJ3/vHfoZMhj3nzMk5e/bsmc9rrbMys589s58n+5z5zvNjf7ciAjMzs0bblV0BMzPrTA4QZmaWywHCzMxyOUCYmVkuBwgzM8vlAGFmZrkcIKyrSLpd0gfadKwLJX29HccyK4MDhFWOpMckbZK0UdIvJV0taceCjzlb0poijzHEsadI+g9Jv5K0QdIPJc2qKz9T0jJJz0paI+kzkravK99V0g2SfiPpcUnvbHj/d6btv5G0WNKu7WyfdS4HCKuqt0fEjsBRQD/wdyXXp0gbgfcBU4FdgE8D36kLAq8GPgJMAY4B/gQ4r+71XwBeAHYH3gVcIekQgPTv/wHencqfB/6l4PZYRThAWKVFxCDwPeDQus37pG/cz0m6RdKUWoGkb0n6haRnJN1R+6BMZTtIulzSz1PP5IuSJkp6TTrGnqnXslHSnullEyQtTMdaKam/7v0OSkNeG1LZKXVlV0v6gqSb0mvvlvTaJm38bUSsioiXAAGbyQLFrqn8ioj494h4If1/fAOYlY7zGuAdwMcjYmNE3AncSBYQIAsY34mIOyJiI/BxYJ6kSSM9F9Z9HCCs0iTtBZwMLK/b/E7gvcAfABPY+tv094ADUtm9ZB+mNZcCBwJHAvsDfcD/jIjfACcBayNix/SzNr3mFOBaYDLZB+/nU73GA98BbknH+jDwDUkz6o53OvBJsg/71cDFw7T1AeC36Thfioh1TXZ9I7AyPT4QeDEiHqkrvx+oBcZD0nMAIuKnZL2NA4eqi/UGBwirqsWSNgB3Aj8A/ldd2Vci4pGI2ARcR/aBD0BELIiI5yLid8CFwBGSdpYk4GzgbyLi6Yh4Lr3n6cPU486IuDkiNgNfA45I248FdgQuTd/slwLfBc6oe+0NEfGjiHiRLFAdyRAi4nBgJ7IAeGfePpLeRzbkdnnatCPwbMNuzwCT6sqfGaLcetj2w+9i1pHmRsS/NSn7Rd3j58k+BJE0juxb+l+Qjee/lPaZAuxANpa/LIsVQDacM26YejQe61VpbmBP4Ik0LFTzOFmvZMh6DiUifgtcI+lhSfdFxMvf/iXNBS4B3hIRT6XNG8mCSr2dgOdaLLce5h6E9ZJ3AqcCbwF2Bqan7QKeAjYBh0TE5PSzc5oIBxhp2uO1wF6S6v/G9gYGt7XyDcYD+9WeSDoRuIps8n5F3X6PANtLOqBu2xFsGYJayZZeD5L2IwuW9UNS1qMcIKyXTAJ+B/yKrLfw8rBU+qZ/FfBZSX8AIKlP0py0yy+B3STt3OKx7ibrFfytpPGSZgNvJ5uvGBFJx0p6g6QJadL8o2Qrju5O5SeQDVG9IyJ+VP/aNH+yCLhI0mvS8thTyYbDSK97u6Tj04T2RcCiNMRmPc4BwnrJQrJhnkHgIeCuhvKPkk0W3yXpWeDfgBkAEfFj4Brg0bQqaU+GEBEvkAWEk8h6J/8CvCe9z0jtQLZU9Vep7icDf1Y3Uf5xsh7RzXWrrL5X9/pzgInAutSG/xoRK1M9VwIfJAsU68iC6DnbUEfrQvINg8zMLI97EGZmlssBwszMcjlAmJlZLgcIMzPL1TUXyk2ZMiWmT59edjXMzCpl2bJlT0XE1LyyrgkQ06dPZ2BgoOxqmJlViqTHm5V5iMnMzHI5QJiZWS4HCDMzy+UAYWZmuRwgzMwsV9esYira4uWDXLZkFWs3bGLPyROZP2cGc2f2Df9CM7OKcoBoweLlg1ywaAWbfr8ZgMENm7hgUZZy30HCzLqVh5hacNmSVS8Hh5pNv9/MZUtWlVQjM7PiOUC0YO2GTSPabmbWDRwgWrDn5Ikj2m5m1g0cIFowf84MJo7f+t71E8ePY/6cGSXVyMyseJ6kbkFtItqrmMyslzhAtGjuzD4HBDPrKR5iMjOzXA4QZmaWywHCzMxyOUCYmVkuT1JXhHNBmVm7FdaDkLRA0jpJDzYpl6TPSVot6QFJRzWU7yRpjaTPF1XHqqjlghrcsIlgSy6oxcsHy66amXWxIoeYrgZOHKL8JOCA9HM2cEVD+aeAOwqpWcWUnQtq8fJBZl26lH3Pv4lZly51YDLrEYUFiIi4A3h6iF1OBRZG5i5gsqQ9ACS9DtgduKWo+lVJmbmg3Hsx611lTlL3AU/UPV8D9EnaDvgH4Lzh3kDS2ZIGJA2sX7++oGqWr8xcUGX3XsysPJ24iukc4OaIWDPcjhFxZUT0R0T/1KlT21C1cpSZC8qZbM16V5mrmAaBveqeT0vbjgOOl3QOsCMwQdLGiDi/hDp2hDJzQe05eSKDOcHAmWzNul+ZAeJG4FxJ1wLHAM9ExJPAu2o7SDoL6O/l4FBTVi6o+XNmbHU3PXAmW7NeUViAkHQNMBuYImkN8AlgPEBEfBG4GTgZWA08D7y3qLrYtnMmW7PepYgouw5jor+/PwYGBsquhplZpUhaFhH9eWWdOEltZmYdwAHCzMxyOUCYmVkuBwgzM8vlAGFmZrkcIMzMLJcDhJmZ5er5Gwb5Rjydz+fIrBw9HSBqqaxraSRqqawBfwB1CJ8js/L09BCTU1l3Pp8js/L0dIBwKuvO53NkVp6eDhBl3ojHWuNzZFaeng4QZd6Ix1rjc2RWnp6epHYq687nc2RWHqf7NjPrYU73bWZmI9bTQ0xl8wVgZtbJHCBK4gvAzKzTeYipJL4AzMw6nQNESXwBmJl1OgeIkvgCMDPrdA4QJfEFYGbW6TxJXRJfAGZmnc4BokRzZ/aVGhC8zNbMhuIA0aO8zNbMhuM5iB7lZbZmNhwHiB7lZbZmNhwHiB7lZbZmNhwHiB7lZbZm1bd4+SCzLl3KvuffxKxLl7J4+eCYvr8nqXuUl9maVVs7FpoUFiAkLQDeBqyLiENzygX8E3Ay8DxwVkTcK+lI4ApgJ2AzcHFEfLOoevayspfZmtm2G2qhyVj9XRc5xHQ1cOIQ5ScBB6Sfs8mCAmTB4j0RcUh6/T9KmlxgPc3MKqcdC00KCxARcQfw9BC7nAosjMxdwGRJe0TEIxHxk/Qea4F1wNSi6mlmVkXtWGhS5iR1H/BE3fM1advLJB0NTAB+mvcGks6WNCBpYP369YVV1HpX0ZOAZtuqHQtNOnYVk6Q9gK8B742Il/L2iYgrI6I/IvqnTnUnw8ZWbRJwcMMmgi2TgA4S1gnmzuzjknmH0Td5IgL6Jk/kknmHjem8YpmrmAaBveqeT0vbkLQTcBPwsTT8ZNZ27ZgENBuNohealNmDuBF4jzLHAs9ExJOSJgA3kM1PfLvE+lmP89Xm1uuKXOZ6DTAbmCJpDfAJYDxARHwRuJlsietqspVL700v/UvgjcBuks5K286KiPuKqqtZnj0nT2QwJxj4anOr181ZkQsLEBFxxjDlAXwoZ/vXga8XVS+zVs2fM2OrC5HAV5vb1ro9K3LHTlKbla0dk4BWbd2eFdmpNsyG4KvNbSjdPk/lHoSZ2Tbq9qzIDhBmZtuo27Mie4jJStHNKz+sd3R7VmQHCGu7bl/5Yb2lm+epPMRkbdftKz/MuoUDhLVdt6/8MOsWHmKytvMVyjaWPJ9VHPcgrO26feWHtY8z7hbLAcLazlco21jxfFaxPMRkpejmlR/WPp7PKpZ7EGZWWd1+JXPZHCDMrLI8n1UsDzGZWWV1+5XMZXOAMLNK83xWcTzEZGZmuRwgzMwslwOEmZnlcoAwM7NcDhBmZpbLAcLMzHJ5mauZjZozqnYnBwgzGxXfIbB7eYjJzEbFGVW7lwOEmY2KM6p2Lw8xWU/ymPnY8R0Cu5d7ENZzqnQXssXLB5l16VL2Pf8mZl26tCPr6Iyq3csBwnpOVcbMqxLIfIfA7uUhJus5VRkzHyqQddqHrzOqdqfCehCSFkhaJ+nBJuWS9DlJqyU9IOmourIzJf0k/ZxZVB2tN1XlLmRVCWTWvYocYroaOHGI8pOAA9LP2cAVAJJ2BT4BHAMcDXxC0i4F1tN6TFXGzKsSyKx7FRYgIuIO4OkhdjkVWBiZu4DJkvYA5gC3RsTTEfFr4FaGDjRmI1KVMfOqBDLrXi3NQUg6OCIeatg2OyJuH8Wx+4An6p6vSduabc+r19lkvQ/23nvvUVTFek0Vxsx9O00rW6uT1NdJ+hrwGeBV6d9+4LiiKtaKiLgSuBKgv78/yqyLWRGqEMise7U6xHQMsBfwn8A9wFpg1iiPPZjes2Za2tZsu9mQqnDNgFmVtBogfg9sAiaS9SB+FhEvjfLYNwLvSauZjgWeiYgngSXAWyXtkian35q2mTVVlWsGzKqk1QBxD1mAeD1wPHCGpG8N9QJJ1wA/BGZIWiPp/ZI+KOmDaZebgUeB1cBVwDkAEfE08Kl0zHuAi9I2s6aqcvGbWZW0Ogfx/ogYSI+fBE6V9O6hXhARZwxTHsCHmpQtABa0WDczXzNgVoBWA8Q6SY3LhH4w1pUx21ZOGGc29loNEDcBAYhsDmJfYBVwSEH1MhuR+XNmbHXTGmj/NQPOEGvdpqUAERGH1T9PaTHOKaRGZtug7GsGfFc160bblKwvIu6VdMxYV8ZsNMq8ZqBKifXMWtXqldT/ve7pdsBRZNdCWIfxMEc5PElu3ajVHsSkuscvks1JXD/21bHR8DBHeTxJbt2o1TmITxZdERs9D3OUpxMmyc3G2pABQtJ3yFYv5YqIU8a8RrbNOmGYo1eHuMqeJDcrwnA9iMvbUosuU9aHZNnDHL0+xOXEetZthgsQP4uIn7elJl2izA/Jsoc5PMRVnl7tuVmxhsvFtLj2QJInpVtQZk6gsm+E0wlDXL3IiQqtKMP1IFT3eL8iK9Ityv6QLHOYo+whrl5VVM/NvRIbrgcRTR5bE718H2HfIrMcRXwp6YReie/vUb7hAsQRkp6V9BxweHr8rKTnJD3bjgpWTS9/SJY9xNWrivhSUnb69E4IUDbMEFNEjBuq3F6p15c7eiVP+xWxOKHsoVIveOgM25SLyYbmD0lrpyK+lJQ9n1R2gLKMA4RZFxjrLyVlL5kuO0BZptVbjppZDyl7PqmX5/I6iXsQZm1WleWjZQ6V9vpcXqdwgDBro15PRzISnssrn4eYzNqo7OWjZiPhAGHWRl6dY1XiAGHWRr18pb1VjwOEWRt5dY5ViSepzdrIq3OsShwgzNrMq3OsKjzEZGZmuRwgzMwslwOEmZnlcoAwM7NchQYISSdKWiVptaTzc8r3kXSbpAck3S5pWl3ZZyStlPSwpM9JUuPrzcysOIWtYpI0DvgC8KfAGuAeSTdGxEN1u10OLIyIr0o6AbgEeLekPwZmAYen/e4E3gTcXlR9zcxqqpJQsWhF9iCOBlZHxKMR8QJwLXBqwz4HA0vT4+/XlQfwKmACsAMwHvhlgXU1MwN8u9N6RQaIPuCJuudr0rZ69wPz0uPTgEmSdouIH5IFjCfTz5KIeLjxAJLOljQgaWD9+vVj3gAz6z1OqLhF2ZPU5wFvkrScbAhpENgsaX/gIGAaWVA5QdLxjS+OiCsjoj8i+qdOndrOeptZl3JCxS2KDBCDwF51z6elbS+LiLURMS8iZgIfS9s2kPUm7oqIjRGxEfgecFyBdTUzA5xQsV6RAeIe4ABJ+0qaAJwO3Fi/g6Qpkmp1uABYkB7/nKxnsb2k8WS9i1cMMZmZjTUnVNyisAARES8C5wJLyD7cr4uIlZIuknRK2m02sErSI8DuwMVp+7eBnwIryOYp7o+I7xRVVzOzmrLvx91JFBFl12FM9Pf3x8DAQNnVMDOrFEnLIqI/r6zsSWozM+tQTvdtZtZhOuVCPQcIM7MOUrtQr3YtRu1CPaDtQcJDTGZmHaSTLtRzgDAz6yCddKGeA4SZWQfppAv1HCDMzDpIJ12o50lqM7MOUpuI9iomMzN7hbkz+zriym0PMZmZWS4HCDMzy+UAYWZmuRwgzMwslwOEmZnlcoAwM7NcDhBmZpbLAcLMzHI5QJiZWS4HCDMzy+UAYWZmuRwgzMwslwOEmZnlcjZXsx6zePlgR6SSts7nAGHWQxYvH+SCRStevufx4IZNXLBoBYCDhL2Ch5jMeshlS1a9HBxqNv1+M5ctWVVSjayTOUCY9ZBmN75vtt16mwOEWQ9pduP7ZtuttzlAmPWQ+XNmMHH8uK22TRw/jvlzZpRUI+tknqQ26yG1iWivYrJWOECY9Zi5M/scEKwlhQ4xSTpR0ipJqyWdn1O+j6TbJD0g6XZJ0+rK9pZ0i6SHJT0kaXqRdTUzs60VFiAkjQO+AJwEHAycIenght0uBxZGxOHARcAldWULgcsi4iDgaGBdUXU1M7NXKrIHcTSwOiIejYgXgGuBUxv2ORhYmh5/v1aeAsn2EXErQERsjIjnC6yrmZk1KDJA9AFP1D1fk7bVux+Ylx6fBkyStBtwILBB0iJJyyVdlnokZmbWJmUvcz0PeJOk5cCbgEFgM9nk+fGp/PXAfsBZjS+WdLakAUkD69evb1ulzcx6QZEBYhDYq+75tLTtZRGxNiLmRcRM4GNp2way3sZ9aXjqRWAxcFTjASLiyojoj4j+qVOnFtUOM7OeVGSAuAc4QNK+kiYApwM31u8gaYqkWh0uABbUvXaypNqn/gnAQwXW1czMGhR2HUREvCjpXGAJMA5YEBErJV0EDETEjcBs4BJJAdwBfCi9drOk84DbJAlYBlxVVF3NrD2carxaFBFl12FM9Pf3x8DAQNnVMLMmGlONQ5bm45J5hzlIlEjSsojozysre5LazHqEU41XjwOEmbWFU41XjwOEmbWFU41XjwOEmbWFU41Xj7O5mllbONV49ThAmFnbONV4tXiIyczMcjlAmJlZLgcIMzPL5QBhZma5HCDMzCyXA4SZmeVygDAzs1wOEGZmlssBwszMcjlAmJlZLgcIMzPL5QBhZma5nKzPzHqG74k9Mg4QZtYTGu+JPbhhExcsWgHgINGEA4SZ9YSh7ondjgBRxd6LA4SZ9YQy74ld1d6LJ6nNrCeUeU/soXovncwBwsx6Qpn3xC6z9zIaDhBm1hPmzuzjknmH0Td5IgL6Jk/kknmHtWWIp8zey2h4DsLMekZZ98SeP2fGVnMQ0L7ey2g4QJiZFawWlLyKyczMXqGs3stoeA7CzMxyOUCYmVkuBwgzM8vlAGFmZrkcIMzMLJciouw6jAlJ64HHy65HC6YAT5VdiTHUbe2B7mtTt7UHuq9NZbZnn4iYmlfQNQGiKiQNRER/2fUYK93WHui+NnVbe6D72tSp7fEQk5mZ5XKAMDOzXA4Q7Xdl2RUYY93WHui+NnVbe6D72tSR7fEchJmZ5XIPwszMcjlAmJlZLgeIgkh6laQfSbpf0kpJn0zb95V0t6TVkr4paULZdW3VEG26WtLPJN2Xfo4su64jIWmcpOWSvpueV/Yc1eS0qbLnSNJjklakeg+kbbtKulXST9K/u5Rdz5Fo0qYLJQ3WnaOTy66nA0RxfgecEBFHAEcCJ0o6Fvg08NmI2B/4NfD+Eus4Us3aBDA/Io5MP/eVV8Vt8tfAw3XPq3yOahrbBNU+R29O9a5dK3A+cFtEHADclp5XTWObIPu9q52jm0urWeIAUZDIbExPx6efAE4Avp22fxWYW0L1tskQbaosSdOAPwO+lJ6LCp8jeGWbutSpZOcGKniOqsIBokCpm38fsA64FfgpsCEiXky7rAEqdQeRxjZFxN2p6GJJD0j6rKQdSqziSP0j8LfAS+n5blT8HPHKNtVU9RwFcIukZZLOTtt2j4gn0+NfALuXU7VtltcmgHPTOVrQCcNmDhAFiojNEXEkMA04Gvijkqs0ao1tknQocAFZ214P7Ap8tMQqtkzS24B1EbGs7LqMlSHaVMlzlLwhIo4CTgI+JOmN9YWRrdWvWk82r01XAK8lG759EviHEusHOEC0RURsAL4PHAdMllS71es0YLC0io1CXZtOjIgn0/DT74CvkAXDKpgFnCLpMeBasqGlf6La5+gVbZL09QqfIyJiMP27DriBrO6/lLQHQPp3XXk1HLm8NkXEL9MXsJeAq+iAc+QAURBJUyVNTo8nAn9KNmn4feDP025nAv9aTg1Hrkmbflz3hyqyseAHy6tl6yLigoiYFhHTgdOBpRHxLip8jpq06a+qeo4kvUbSpNpj4K1kdb+R7NxAxc5RszbVzlFyGh1wjrYffhfbRnsAX5U0jiwQXxcR35X0EHCtpL8HlgNfLrOSI9SsTUslTQUE3Ad8sMxKjoGPUt1z1Mw3KnqOdgduyOIa2wP/NyL+n6R7gOskvZ8szf9flljHkWrWpq+l5ccBPAb8l/KqmHGqDTMzy+UhJjMzy+UAYWZmuRwgzMwslwOEmZnlcoAwMxslSZdJ+nG6CvqG2nLwJvtulUixoexzkjbWPd9H0m3pfW9PaVRqZXtLukXSw5IekjS9hXp+ONVzpaTPDLe/A4RVgqTd6rJc/qIh6+V/FnTMmZK+nB7Przveg5I2p4yie0n6fvoDXSnpr+te31LGUUlnpn1+IunMJvs8JmlKEe1scrzLJZ3QruNViaTZkq5u2HwrcGhEHA48QnblejN5iRSR1A80/o5cDixM73sRcEld2ULgsog4iOyiuiEvFpT0ZrIcVkdExCHpvYcWEf7xT6V+gAuB89pwnG+lP6bG7W8nuwANsmtDjkqPJ5F9OBycnn8GOD89Ph/4dM577Qo8mv7dJT3eJWe/x4ApY9w+Ads1KdsHuKXsc92JP8Bs4Oohyk8DvtGkbBpZ9tkTgO/WbR9HdoHmHsDGuu0rgb3qztez6fHBwJ1NjvE64AfAMmAJsEfafh3wlpG01T0Iq7xalzx9s/uBpH+V9KikSyW9S9k9LFZIem3ab6qk6yXdk35m5bznJODwiLg/55BnANcARJbC4t70+Dmyb4a15H6tZBydQ5b08OmI+DXZN9ETmzT1w5LuTW35o1TPXSUtTkMQd0k6PG2/UNJ5de15UNL09LNK0kKyK3X3UnaviAfT+/5NasvjwG6S/rBJXay59wHfa1LWLJHiucCNsSUBYc39wLz0+DRgkqTdgAOBDZIWpeGqy9LQ1Xjgn4E/j4jXAQuAi9PrDwSOV3avkx9Iev1wDfGV1NZtjgAOAp4m+zb+pYg4Og39fBj4CFm+pc9GxJ2S9ib7lnVQw/v0k5PqQNKryT7Az80pmw7MBGoZblvJONoHPFH3fKjssU9FxFGSzgHOAz4AfBJYHhFz05DQQrJkb0M5ADgzIu6S9DqgLyIOTW2oHzu/lyy30/XDvF9PkHQ3sAOwI7CrsqzGAB+NiCVpn48BLwLfyHn9y4kUJc2u274n8BdkPZNG5wGfl3QWcAdZXrDNZJ/dx5P9vv0c+CZwFtnv3qHArelK7XFkif9Ir9kVOJYsaeN1kvaL1L3I4wBh3eae2oeypJ8Ct6TtK4A3p8dvAQ5Of0AAO0naMbbc6wKyrv76nPd/O/AfEfF0/UZJO5J9kH4kIp5tfFFEhKTRpi1YlP5dxpZvlW8A3pGOsTTN1ew0zPs8HhF3pcePAvtJ+mfgJrb8f0E2pr3nKOvcNSLiGMh6qsBZEXFWfXn6EH8b8CdNPnRriRRPBl5F9nv3dbLe6P7A6vQ7+WpJqyNi/4hYSzrX6XfsHRGxQdIa4L6IeDSVLSb74P8RsDIijss5/hpgUarbjyS9BEwh//cc8CS1dZ/f1T1+qe75S2z5QrQdcGxsuXNXX0NwANhE9kfc6HTS8FJN6tZfTzbuvKiuqJWMo4PAXnXPh8oeW2tL7RvkUF5k67/v+rb8pvYgDWsdAdxOlp/pSw2v2TTMcQyQdCLZ0NEpEfF83j7RJJFiRNwUEX8YEdNT2fOR3c0QSVMk1c7jBWRDRgD3kGUdnpqenwA8BKwCpko6Lr1+vKRD0j6LSV+SJB0ITACeGqpdDhDWi24hG24CQPn3Z36Y7FsddfvtDLyJusyhyr7yfRl4OCL+d8N7tJJxdAnwVkm7pFVOb03bWvXvwLtSXWaTDUM9SzapfVTafhSwb96L08qo7SLieuDvaq9JDqQDMopWxOfJFincmla6fRGy4SNJo7l16GxglaRHyIYoL4bsvixkw0+3SVpBNoF9VUS8QJaJ+NOS7idLzPjH6b0WkPUWHyRLBX/mUMNL4CEm603/DfiCpAfI/gbuoCG7aUT8WNLOkialyWfIJglviYjf1O06C3g3sKJuTPp/RHY/4UvJyTialjN+MCI+EBFPS/oU2TdCgIsah6+GcSGwILXlebYEpOuB90haSTYu/UiT1/cBX2n4llrrFe0PDIygLj0hIm4n63HVb9u/yb5rgZNbeY+6sh3rHn+bLbe/bdzvVuDwnO33AW/M2f4C8Fd579WMs7maNZFW9DwXEd18b+dckk4jW7778bLrYuXxEJNZc1ew9ZxGL9meDrjlpZXLPQgzM8vlHoSZmeVygDAzs1wOEGZmlssBwszMcjlAmJlZrv8PyirB5288nmUAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = 1440*2*30 # 30 days\n", "\n", "plt.plot(time[0: x: 2880], flux[0: x: 2880], \"o\")\n", "plt.xlabel(\"Time (%.1f hours)\"%(x*0.5/60))\n", "plt.ylabel(\"Flux\")\n", "plt.title(\"Phaethon 3200\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visually, it's a nearly impossible to guess what the rotational period should be, and even Lomb-Scargle Periodograms sometimes have a hard time determining what the period might be if the lightcurve is sparse enough. \n", "\n", "It was with these sparse incoming datasets in mind that we decided to develop this method, named AsteroGaP for **Astero**id **Ga**ussian **P**rocesses." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Gaussian Processes\n", "\n", "A Gaussian Process (also known as Gaussian Process Regression) is a generative kernel-based framework for solving regression problems [(Roberts et al. 2012)](https://ui.adsabs.harvard.edu/abs/2012RSPTA.37110550R/abstract). Rather than modeling the data points directly as in standard linear regression, a Gaussian process models the *covariance* between data points.\n", "\n", "Run the following cells to view an interactive visualization of how GP models are fit." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from ipywidgets import interact\n", "import george\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def create_data(gamma = 10, period = 1, amp = 1):\n", " \"\"\"\n", " Generate some periodic data using a periodic kernel.\n", " \"\"\"\n", " \n", " # generate 10 random numbers for x\n", " x = 10 * np.sort(np.random.rand(10))\n", " \n", " # determine the y error and make it an array as long as x\n", " # have the error scale with the amplitude\n", " yerr = 0.2 * amp * np.ones_like(x)\n", " x_possible = np.linspace(0,10, 1000)\n", "\n", " # create a sinusoidal plot based on inputs\n", " # establish the kernel type we are using\n", " kernel = amp * george.kernels.ExpSine2Kernel(gamma, period)\n", " gp = george.GP(kernel)\n", " \n", " # generate a ExpSine2 function using our x-values (0-10) and our y error\n", " gp.compute(x, yerr)\n", "\n", " # sample the y-values from the function we made\n", " y_possible = gp.sample(x_possible) # a subset of possible x-values\n", " \n", " #return our simulated data with the original function\n", " return(x_possible, y_possible, gp)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def sample_data(x_possible, y_possible, yerr_amp, random_n=10):\n", " \"\"\"\n", " Samples a n-sized random array of points from the sample data generated in\n", " the create_data function or another type of light curve sample.\n", " \"\"\"\n", " \n", " idx = y_possible.size\n", " \n", " # randomly select n points from 1-500\n", " idx_choice = np.random.choice(idx, random_n, replace= False)\n", " idx_choice.sort()\n", " \n", " # filter out the randomly chosen indicies from earlier\n", " x = x_possible[idx_choice]\n", " \n", " # get the predicted y values corresponding to x values\n", " y_base = y_possible[idx_choice]\n", " \n", " # generate some uncertainty\n", " yerr = yerr_amp * 0.01 * np.ones_like(x) #turn 0.2 into an input (1%)\n", " \n", " # add noise to initial y values\n", " y = y_base + yerr * np.random.randn(len(x))\n", "\n", " return(x, y, yerr)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def plotting_gp_ex(x, y, yerr,pre_y, x_possible, gp):\n", " \"\"\"\"\"\"\n", " plt.figure(figsize=(15,10))\n", " pred, pred_var = gp.predict(y, x_possible, return_var=True)\n", " plt.fill_between(x_possible, pred - np.sqrt(pred_var), pred + np.sqrt(pred_var),\n", " color=\"red\", alpha=0.4)\n", " plt.plot(x_possible, pred, \"red\", lw=1.5, alpha=0.7)\n", " plt.errorbar(x, y, yerr = yerr, fmt=\"ok\", capsize=0, )\n", " plt.plot(x_possible,pre_y)\n", " plt.xlim(x_possible.min(), x_possible.max())\n", " plt.ylim(pre_y.min()-0.5, pre_y.max()+0.5)\n", " plt.xlabel(\"x\")\n", " plt.ylabel(\"y\");\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def fit_data(x, y, yerr, gamma, period, gp):\n", " \n", " gp.compute(x, yerr)\n", " x_possible = np.linspace(0,10,500)\n", " pred, pred_var = gp.predict(y, x_possible, return_var=True)\n", " ln_likelihood = gp.log_likelihood(y)\n", "\n", " #print(ln_likelihood)\n", "\n", " return gp, ln_likelihood" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def data_comparison(gamma, period, amp):\n", " \n", " x_pos, py, gp1 = create_data(gamma, period, amp)\n", " \n", " def _do_work(n = 10):\n", " x, y, yerr = sample_data(x_pos, py, 1, n)\n", " \n", " gp2, ln_likelihood = fit_data(x, y, yerr, gamma, period, gp1)\n", " \n", " plotting_gp_ex(x,y,yerr,py,x_pos,gp2)\n", "\n", "\n", " return _do_work\n", "\n", "#n will only work if m is set equal to 0" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "5642754004ba43fd8cdf33899b271b35", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(IntSlider(value=10, description='n', max=20, min=2, step=2), Output()), _dom_classes=('w…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "._do_work(n=10)>" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vary_nm = data_comparison(gamma=10, period=np.log(3), amp=1)\n", "interact(vary_nm, n=(2, 20, 2), continuous_update=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The main thing to focus on is this interactive plot down here at the bottom. For the last cell, there is a periodic lightcurve in **blue** set up to repeat every 3 units (the kernel takes the logarithm of the period, hence the log(3) notation), with a gamma of 10, and an amplitude (amp) of 1. Feel free to test out different values for gamma, the period, and the amplitude to get a feel for how the lightcurve changes.\n", "\n", "The number **n** in the bar above is the number of **black** observations/samples of the blue lightcurve we are taking. Once we have n number of randomly selected observations (time/x, flux/y, and flux error), we can then feed those samples into a Gaussian Process and see how the model would generate a fit given what it knowns.\n", "\n", "As you can hopefully see, with few points, the GP model is only certain with its fit in specific areas (small red shading), and is very uncertain about what happens in other parts of the period (large red shading). The red line is what the predicted fit is, and the shading is the uncertainty. \n", "\n", "If we increase the number of points, the GP model has more information and can make better estimates about what is going in with the lightcurve, reducing the uncertainty. \n", "\n", "Even with just a few observations (~10-20), if the lightcurve has been broadly and evenly sampled, the GP model is able to produce a really great fit for the lightcurve data. \n", "\n", "## Kernels\n", "\n", "The main component driving this Gaussian Process model is the kernel, which is a matrix telling the model how the different data points should relate to one another. \n", "\n", "We have specified an Exponential Sine Squared kernel, which is significant because this kernel expects data to repeat after a certain amount of time, like a period. " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "%matplotlib notebook\n", "import matplotlib.pyplot as plt\n", "import matplotlib.gridspec as gridspec\n", "\n", "import seaborn as sns\n", "sns.set_style(\"white\")\n", "sns.set_palette('viridis')\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "import george\n", "from george import kernels" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": false }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
');\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('