{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Light Curve Analysis Module\n",
"\n",
"**Lecturer:** Melissa Hayes-Gehrke \n",
"**Jupyter Notebook Authors:** Melissa Hayes-Gehrke & Cameron Hummels\n",
"\n",
"This is a Jupyter notebook lesson taken from the GROWTH Summer School 2020. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-school-2020.html\n",
"\n",
"## Objective\n",
"Identify periodic behavior in noisy observational data to generate light curves\n",
"\n",
"## Key steps\n",
"- Extract a periodic signal from noisy data when we know the period of the oscillation\n",
"- Use the Lomb Scargle algorithm to extract a light curve from data when we do not know the period\n",
"- Apply the Lomb Scargle to real Cepheid observational data\n",
"\n",
"## Required dependencies\n",
"\n",
"See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with `pip install `. The external astromatic packages are easiest installed using package managers (e.g., `rpm`, `apt-get`).\n",
"\n",
"### Python modules\n",
"* python 3\n",
"* astropy\n",
"* numpy\n",
"* matplotlib\n",
"\n",
"### External packages\n",
"None"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from astropy.timeseries import LombScargle\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import os"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"cwd = os.getcwd()\n",
"data_dir = os.path.join(cwd, 'data')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Generate a Synthetic Lightcurve\n",
"\n",
"First, we need to generate a lightcurve to work with. The function below generates a series of timestamps t, and a sinusoidal signal y with noise dy added to it. We will use this lightcurve to test our functions. Note that our y values are generated with a sin function with a frequency of 2pi, or equiaveletly, a period of 1."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'Flux')"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAASZUlEQVR4nO3dfawld13H8feHLQVUkIeuUtquW7Q+oMilXspjyFUhKWi6PhGKT5VAVhIqihotkvCUGCEhPqANWqFSjaE0RWHVRpSHVRMC7lZuhLZWlqL02krLQ0EEgdKvf5xZenI5d/Y+nHNmzpn3K7nZM3Nm73znzjnzmd/vN2dOqgpJkrZyn64LkCT1m0EhSWplUEiSWhkUkqRWBoUkqdVpXRcwbWeccUYdPHiw6zIkaaFcf/31n6iq/ZOeW7qgOHjwIMePH++6DElaKEn+c6vn7HqSJLUyKCRJrQwKSVIrg0KS1MqgkCS1MigkSa0MCklSK4NCktTKoJB6bm1tjbW1ta7L0IAZFJK0A0MMboNC0twN8WC7yAwKSVIrg0KS1MqgkCS1MijUC/ZZS/1lUOirPFhLmsSgkCS1Mig0V7ZapMVjUEhTYghqWRkUPTbNA48Hsf6b1T5y32+Pf6etGRSSlp4hsDcGhaTe8IDeTwaFtAUPWtKIQSFJamVQSANli0nbZVBI2jFDZlg6DYokVya5I8mHtng+SV6X5ESSf01y/rxr3A3fRJKWSdctijcBF7Y8/wzgvObnMPD6OdQkDY4nN/O3SH/zToOiqv4R+FTLIoeAP62R9wEPTnLmfKpbHmtra6yvr3ddhqQF1XWL4lTOAm4dm95o5kmS5qTvQZEJ8+prFkoOJzme5Pidd945h7IkaTj6HhQbwDlj02cDt21eqKquqKrVqlrdv3//3IqTpCE4resCTuEIcGmSq4HHA5+pqts7rkmaq6NHj3ZdQm+dHAye599oiPuj06BI8mZgDTgjyQbwcuC+AFX1h8B1wDOBE8Dnged2U6kkDVenQVFVzznF8wW8cE7lSJIm6PsYhSSpY30fo5C0hIbYz7/IDAp9lW/evfHvp2Vl15MkqZUtCklLz9be3tiikCS1MigkSa3selIv2DUg9ZdBMQBHjx6d633vu7itwjKY1d/L/dBPi7RfDAppoMYPVIt00NL8GRSSemOngeUXcs2HQaG58sxVWjwGRY95UF0esxi3cSxI82JQSNqxvoTTysrK1H5XX7apjwyKgfBNIGm3/MCd5mZtbW2ul+lKmg6DYsY8OKovfC1qt+x6krZgd500YotCktTKFoVmyks4NUt+4G4+bFHskP280+ffVOo3WxQz4NmzpGVii6KFZ7qSZFBIkk7BridNnV1vmiUvkJg/WxSSpFYGhSSplUEhaaq8CGT5GBQLyDeipHkyKKQOGfpq05fXR6dBkeTCJDcnOZHksgnP/1ySO5OsNz/P76LORdSXF5ikxdfZ5bFJ9gGXA08HNoBjSY5U1Y2bFn1LVV069wI1cwaZtBi6/BzFBcCJqroFIMnVwCFgc1BowLxmfmv+TTQvXQbFWcCtY9MbwOMnLPfjSZ4K/Dvw4qq6dfMCSQ4DhwEOHDgwg1Il9dE0vzNbW+tyjCIT5tWm6b8CDlbV9wLvBK6a9Iuq6oqqWq2q1f3790+5zOXgmIWk3eoyKDaAc8amzwZuG1+gqj5ZVV9sJv8Y+L451aYFZihK09Vl19Mx4Lwk5wL/BVwM/OT4AknOrKrbm8mLgJvmW+LeDf2LVexHlxZfZ0FRVXcnuRR4B7APuLKqbkjyKuB4VR0BXpTkIuBu4FPAz3VV727Zhypp0XV699iqug64btO8l409fgnwknnXJUm6l7cZl+bALjgtMm/hoc6sr68PfgxHWgS2KKSBsFWj3bJFIUlqZVBIklrZ9aS5setDWky2KCRJrQwKSVIrg0LqEe9TpT5yjGKH7GefL//ei8d9tnwMCnXOA4vUbwaFpIUyfmLhScZ8OEYhSWpli6KFZysCv7db3enLa86gUGf8rg5pMdj1JElqZYtC6lBfuhakNgbFAvLgImme7HqSJLUyKCRJrQwKSZqyWd+za973BDMoJEmtHMyesa4Gnh3wljQtBsVALFpw+GloqT8MCvXS+vr6rv+v4SJN17bGKJI8asK8talXsyT88hlJy2S7g9nXJPn1jDwgye8DvzXLwiRJ/bDdrqfHA68B3gs8EPhz4MmzKkrDYBeRtBi2GxRfBr4APAC4P/DRqrpnZlVJA2NXpfpsu11PxxgFxeOApwDPSXLtzKqSJPXGdoPieVX1sqr6clX9d1UdAt6+15UnuTDJzUlOJLlswvP3S/KW5vn3Jzm413VKknZmu11PdyQ5sGneP+xlxUn2AZcDTwc2gGNJjlTVjWOLPQ/4dFV9W5KLGY2TPHsv65Uk7cx2g+JvgALCaIziXOBm4Lv3sO4LgBNVdQtAkquBQ8B4UBwCXtE8vhb4gySpqtrDeiX1kB+y7K9tBUVVPXp8Osn5wM/vcd1nAbeOTW8wurpq4jJVdXeSzwAPAz6xqZ7DwGGAAwc2N3ykvfHA1W/jAWPYzMaubgpYVf/CaGB7LzLpV+9iGarqiqpararV/fv377GsfvBDe5L6YlstiiS/PDZ5H+B84M49rnsDOGds+mzgti2W2UhyGvCNwKf2uF5J0g5st0XxwLGf+zEaszi0x3UfA85Lcm6S04GLgSObljkCXNI8/gng3Y5PTI+tlvny761Ftd0xildOe8XNmMOlwDuAfcCVVXVDklcBx6vqCPBG4M+SnGDUkrh42nVI0rQt2xhJa1Ak+SsmjAmcVFUX7WXlVXUdcN2meS8be/x/wLP2sg5J0t6cqkXx2rlUIQlYvjNRLYdTBcVHq+pjc6mkx7zkTlouvqd35lSD2W87+SDJW2dcy8JxcFLSEJyqRTH+OYZHzrKQPpjWWYZnKZKWyalaFLXFY0nSQJyqRfGYJJ9l1LJ4QPOYZrqq6kEzrU5fw75VaTEt8nu3NSiqat+8CpGkoetrmGz37rGSpJ5YX1+f6/p2dVNASdJwGBSSpFZ2PamX7rrrrq5LkNQwKCRpwaysrMx1fXY9SZJaGRSSpFZ2PUk90Lfr5qVxBoWkXthtWBqys2fXkySplUEhSWpl15OkpWE31GwYFD3lC15SX9j1JElqZVBI6gW/Wri/7HraA7uHJA2BQbENfQ+E3XzZydraGuvr63O/Z4ykxWNQSNKC6KprzqCQ5qTvLVNpKwaFpMExtHfGoBizCC+eada4srKyENssqVsGhST1RF9P3Dr5HEWShyb5+yQfbv59yBbLfSXJevNzZN51ShoGP8PRrqsP3F0GvKuqzgPe1UxP8oWqWml+Lppfed1atBftotUraWe6CopDwFXN46uAH+moDknSKXQ1RvHNVXU7QFXdnuSbtlju/kmOA3cDr66qt01aKMlh4DDAgQMHZlFvr62vr3ddgqQp280HaWdlZkGR5J3Awyc89dId/JoDVXVbkkcC707ywar6yOaFquoK4AqA1dXV2lXBkqSJZhYUVfW0rZ5L8vEkZzatiTOBO7b4Hbc1/96S5CjwWOBrgkKSNDtdjVEcAS5pHl8CvH3zAkkekuR+zeMzgCcDN86tQkmaoqNHj/aiG2k3uhqjeDVwTZLnAR8DngWQZBV4QVU9H/gu4I+S3MMo0F5dVQaFpMGbd+B0EhRV9UngByfMPw48v3n8XuDRcy5NkrSJX1zUIT9/IGkRGBSSpFYGhSSplUGhhWb3nTR7BoUk7dDQTlAMij0Y2otF0jAZFFoa8wzuLk4SprVOT3C0UwbFFPkGlLSMDApJg+HJ3O74VagDtaj3nJE0f7YoBsyzK0nbYYuihzzbn74+fQmMtGgMiiWwsrLSdQnSthjYi8mgWAK+6STNMoQNCkm94AlPfxkUc7ZV6i9yk3wRa5YWUVfvNa96kqQpWdYrCW1RSAOxjAcwzYdBIU2wyF2BWg59eu0ZFJIGr08H5T4yKCRph4YWLA5mS5Ja2aLYg6GdVUgaJlsUkuZqfX3dK7AWjC0KLbT19fVO1mtrUkNii0KS1MoWxRR5ltkt//7SbBgUHfLANj/+raXds+tJktTKoJAkteokKJI8K8kNSe5Jstqy3IVJbk5yIsll86xRGrplvRPqsjp69OjMuli7GqP4EPBjwB9ttUCSfcDlwNOBDeBYkiNVdeN8StQi8Gtgt2dtbY319XVWVlYcr9GOdRIUVXUTQJK2xS4ATlTVLc2yVwOHAINCg+QBXl3p81VPZwG3jk1vAI+ftGCSw8BhgAMHDsy+MkmaYFnDfGZBkeSdwMMnPPXSqnr7dn7FhHk1acGqugK4AmB1dXXiMn3XxQtsWV/UkqZrZkFRVU/b46/YAM4Zmz4buG2Pv1PSgHlytDt9vjz2GHBeknOTnA5cDBzpuCZJGpyuLo/90SQbwBOBv0nyjmb+I5JcB1BVdwOXAu8AbgKuqaobuqhXkoasq6ue/hL4ywnzbwOeOTZ9HXDdHEubOZu+GrKjR4/62YwF1OeuJ0lSDxgU0kD4YTvtVp8/RyFpCRlWi8cWhSSplS0KaQLPeqV72aKQJLWyRaGF5pm/NHu2KCRJrQwKSVIrg0KS1MoxCkkTOf6jk2xRSJJaGRSSpFYGhSSplUEhSWplUEiSWnnVkzQAXsGkvbBFIUlqZVBIkloZFJKkVgaFJKmVQSFJamVQSJJaGRSSpFYGhSSplUEhSWqVquq6hqlKcifwn7v4r2cAn5hyOYvA7R4Wt3tYdrLd31JV+yc9sXRBsVtJjlfVatd1zJvbPSxu97BMa7vtepIktTIoJEmtDIp7XdF1AR1xu4fF7R6WqWy3YxSSpFa2KCRJrQwKSVKrwQdFkguT3JzkRJLLuq5nVpKck+Q9SW5KckOSX2zmPzTJ3yf5cPPvQ7qudRaS7EvygSR/3Uyfm+T9zXa/JcnpXdc4bUkenOTaJP/W7PcnDmF/J3lx8xr/UJI3J7n/su7vJFcmuSPJh8bmTdzHGXldc6z71yTnb3c9gw6KJPuAy4FnAI8CnpPkUd1WNTN3A79SVd8FPAF4YbOtlwHvqqrzgHc108voF4GbxqZfA/xOs92fBp7XSVWz9XvA31bVdwKPYbT9S72/k5wFvAhYrarvAfYBF7O8+/tNwIWb5m21j58BnNf8HAZev92VDDoogAuAE1V1S1V9CbgaONRxTTNRVbdX1b80j/+H0UHjLEbbe1Wz2FXAj3RT4ewkORv4IeANzXSAHwCubRZZuu1O8iDgqcAbAarqS1V1FwPY38BpwAOSnAZ8HXA7S7q/q+ofgU9tmr3VPj4E/GmNvA94cJIzt7OeoQfFWcCtY9MbzbylluQg8Fjg/cA3V9XtMAoT4Ju6q2xmfhf4NeCeZvphwF1VdXczvYz7/ZHAncCfNF1ub0jy9Sz5/q6q/wJeC3yMUUB8Brie5d/f47bax7s+3g09KDJh3lJfL5zkG4C3Ar9UVZ/tup5ZS/LDwB1Vdf347AmLLtt+Pw04H3h9VT0W+F+WrJtpkqY//hBwLvAI4OsZdblstmz7ezt2/bofelBsAOeMTZ8N3NZRLTOX5L6MQuLPq+ovmtkfP9n8bP69o6v6ZuTJwEVJ/oNR1+IPMGphPLjpmoDl3O8bwEZVvb+ZvpZRcCz7/n4a8NGqurOqvgz8BfAkln9/j9tqH+/6eDf0oDgGnNdcEXE6o0GvIx3XNBNNv/wbgZuq6rfHnjoCXNI8vgR4+7xrm6WqeklVnV1VBxnt33dX1U8B7wF+ollsGbf7v4Fbk3xHM+sHgRtZ8v3NqMvpCUm+rnnNn9zupd7fm2y1j48AP9tc/fQE4DMnu6hOZfCfzE7yTEZnmPuAK6vqNzsuaSaSPAX4J+CD3NtX/xuMximuAQ4wepM9q6o2D44thSRrwK9W1Q8neSSjFsZDgQ8AP11VX+yyvmlLssJoAP904BbguYxODpd6fyd5JfBsRlf6fQB4PqO++KXb30neDKwxup34x4GXA29jwj5ugvMPGF0l9XnguVV1fFvrGXpQSJLaDb3rSZJ0CgaFJKmVQSFJamVQSJJaGRSSpFannXoRSZMkeRijm64BPBz4CqPbZgB8vqqe1Elh0pR5eaw0BUleAXyuql7bdS3StNn1JM1Aks81/64l+Yck1yT59ySvTvJTSf45yQeTfGuz3P4kb01yrPl5crdbIN3LoJBm7zGMvg/j0cDPAN9eVRcw+tT0LzTL/B6j70t4HPDjzXNSLzhGIc3esZP31EnyEeDvmvkfBL6/efw04FGjuywA8KAkD2y+O0TqlEEhzd74PYXuGZu+h3vfg/cBnlhVX5hnYdJ22PUk9cPfAZeenGhu6Cf1gkEh9cOLgNXmS+9vBF7QdUHSSV4eK0lqZYtCktTKoJAktTIoJEmtDApJUiuDQpLUyqCQJLUyKCRJrf4fZ/WEW1DIGgkAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# \"Generates a sample lightcurve with times t, magnitudes y, and errors dy\n",
"\n",
"rand = np.random.RandomState(42)\n",
"t = 100 * rand.rand(100)\n",
"y = np.sin(2 * np.pi * t) + 0.1 * rand.randn(100)\n",
"dy = 0.1 * (1 + rand.rand(100))\n",
"\n",
"plt.errorbar(t,y,dy,ls='none',c='k')\n",
"plt.xlabel('Time')\n",
"plt.ylabel('Flux')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can see, this lightcurve is not very interesting when we plot it. This is where we need to use period finding analysis to learn more. Below, we try plotting the function in a new way--phase folding it. When we do this, we assume a period, and compute what phase each timestamps corresponds to between 0 and 1, assuming that period. Below, we try folding our lightcurve with the period equal to 1."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAYOElEQVR4nO3dfbBjdX3H8c+nQRENbcDFSoHtQksfpGqkF22VadOqFemU1VqnWDuKRVdsse040xGHKbX+I32YaceWqa7KqB1HoPi01nUQxdRuGSiLE56LrovKdqmsINiMFCV++0dO8HBJzk1ucp6S92vmzibnnNx879nvyff8fr+T33FECACASX6k7AAAANVGoQAAZKJQAAAyUSgAAJkoFACATIeVHcCibdmyJbZt21Z2GABQKzfeeOO3IuKYceuWrlBs27ZNe/fuLTsMAKgV21+ftI6uJwBAJgoFACAThQIAkIlCAQDIRKEAAGSiUAAAMlEoAACZKBQAgEwUihrpdDrqdDplh4ElRo5hHArFEuDgRhHIs9VFoVgiHMgoAnm2eigUAIBMFAoAQCYKBQAgE4WiYPTvIm/kGBaNQlFRHOzIGzmGaVEoAACZKBQVwhke8kaOYTMoFBXT6/U4kJErcgyzolBUwKSzvD179qjVaj1m2bQHOWeOSCPHMA8KRYX0ej31+/2yw8CS6/f76vV6ZYeBGqFQlKzT6XDQIne9Xo88w6ZRKEpEkUAROp0OLVXMhUIBAMh0WNkBrCpaE8gbA81YFFoUFdfv9zngkTuuYEKWUguF7Utt32v71gnrbftdtvfZvtn2qUXHWHd8AKAI5NlyK7tF8QFJZ2Ssf6mkk5OfHZL+qYCYCpHXpbB8mQppeV3txFVUq6XUQhERX5R0f8Ym2yV9KIauk9SyfWwx0eWr3+9rMBiUHQaWGN/LwaKU3aLYyHGS7k49P5AsAwAUpOqFwmOWxeM2snfY3mt776FDhwoICwBWR9ULxQFJJ6SeHy/p4PqNImJnRKxFxNoxxxxTWHAAsAqqXih2SXpNcvXTL0l6MCLuKTuoPA0Gg4mDhN1uV+12e+zybrebc2RYFoPBYOLYxaQcG60jz1ZTqV+4s/0RSR1JW2wfkPQXkp4gSRHxbkm7JZ0paZ+k70p6XTmRAsDqKrVQRMSrNlgfkv6ooHAAAGNUvesJAFAyCsUSyepfBhaFPFs9FIqKaTQaarfb6na7ajabZYeDJdRoNNRsNskxTI1CAQDIRKGokHa7zRkectdsNuk6wkwoFACATBQKAEAm7nBXsNE3Ww87jF2PfIxyrNVqlRsIlgafVjkZ3RNi0pQHp59++obz+bfb7cdtwxQKSMvKs/Q4xJ49e8a+ftKUMEAaXU8lyetadObjwUieucB3KVYLhSIHnU5n6rt/tdttDjhsyix5xpVOmAeFokI4S0PeyDFsBmMUOZj1FpRFdhVtNHaC+pglz0bf9i8CObZ8KBQVkz64ONCQB3IMs6JQLKnRB8Do7A5YNIrM6mCMYsE6nc5M3U7AZpBnKBItigWa5SqUqqA/uV7q2EIkx+qPQlEiDhwUgTzDvOh6AgBkolDMqdVqMacOctXpdNRqtWrZ7YTlQKFYgH6/v1QHcafTWaq/Zxn0+/3ajX9lIcfqhTGKJUf/NIpAni03WhQz4CwIeSPHUEUUCgBAJrqeNmmjs75ms0lzHHOZpmVBnqEItCgAAJkoFACATHQ94VEMoiJv5Fg90aKYU7/f12AwKDsMLLFer0eOoVQUik2o4+R/0+j1ekv5d9UROYYqKbVQ2D7D9p2299m+YMz6c2wfst1Lfl5fRpzS8h64qA5yDFVV2hiF7YakSyS9WNIBSTfY3hURt6/b9PKIOL/wAFcQ9zdA3sixeiqzRfFcSfsiYn9EfE/SZZK2lxjPzCb1HTebTW5gj4WZ9OFKnqEoZRaK4yTdnXp+IFm23its32z7StsnjPtFtnfY3mt776FDh/KIdSZF3sgeq4s8Q1HKLBQesyzWPf+UpG0R8SxJn5P0wXG/KCJ2RsRaRKwdc8wxCw4TacxFhLyRY9VTZqE4ICndQjhe0sH0BhFxX0Q8nDx9r6RfLCg2iAMWxSDPqq/ML9zdIOlk2ydK+m9JZ0v6vfQGto+NiHuSp2dJuqPYEDfWaDTKDiE3XIFTDc1mUw8++GDZYeSGPKu+0gpFRDxi+3xJV0lqSLo0Im6z/Q5JeyNil6Q/tn2WpEck3S/pnLLiXVb0cSNv5Fj9lTqFR0TslrR73bKLUo/fJultRccFAPgh5npaMZzdIW/k2PJhCg9kWrb7gaOayLNqo1AsULfb5QtQyFW32+WMHYWjUAAAMlEoAACZGMxeMLoFUATyDEWiRQEAyEShAABkolBgatxYB3kjx6qJMYo5tNtt9Xq9pbkZS7PZlPTD/u9Wq1ViNBhpt9vas2dP2WEsxPocQz3QoliAZrNJ4iNX3KQIZaJQbNLopjHtdpsDGLlI35iImxShTBQKAEAmxijWGc03M+7sjbM6LAI5hrqhUOBRdKEhb+RYPdH1BADIRIsCj6LLA3kjx+qJQjEnEh95I8dQNrqeAACZKBQAgEwUCsyk3+8zF88mdDodbvU5JXJsc/LMMQoFACATg9lTWsUBxdGkhyjGKuaYRJ7VAS0KzGTc5HR0q2CRJk2ASJ6Vh0KBiUaTHk6LAxmbQZ5V31RdT7afERG3r1vWiYhuLlGhFjhYUQTyrHzTtiiusP1WDx1h+x8kvTPPwAAA1TDtYPbzJP2VpGslHSnpw5JekFdQqI5VHWBFsUZ5RuuhmqZtUXxf0kOSjpD0JEl3RcQPcouqJNyvF3kjx1BH0xaKGzQsFKdJOl3Sq2xfmVtUAIDKmLZQnBsRF0XE9yPifyJiu6RPzvvmts+wfaftfbYvGLP+cNuXJ+uvt71t3vcEAMxm2jGKe21vXbfs3+Z5Y9sNSZdIerGkA5JusL1r3dVV50r6dkT8tO2zNRwn+d153hcAMJtpC8WnJYUkazhGcaKkOyWdMsd7P1fSvojYL0m2L5O0XVK6UGyX9Pbk8ZWS/tG2IyLmeF8gV1m3OgUWoegcm6pQRMQz089tnyrpjXO+93GS7k49P6Dh1VVjt4mIR2w/KOmpkr61Lp4dknZI0tat6xs+s+GexbPp9Xrq9/sTv02LxyPHZrdnzx5J0umnn15yJKtpU3M9RcSXbJ8253t73K/exDaKiJ2SdkrS2toarQ2UiquakLeic2zab2a/JfX0RySdKunQnO99QNIJqefHSzo4YZsDtg+T9GOS7p/zfQEAM5j2qqcjUz+HazhmsX3O975B0sm2T7T9RElnS9q1bptdkl6bPP4dSdcwPlE9g8FA/X6/7DAqpd/v8+WxBRoMBhoMBmWHUSlF5ti0YxR/ueg3TsYczpd0laSGpEsj4jbb75C0NyJ2SXq/pH+2vU/DlsTZi44Ds6F/HXkjx6ons1DY/pTGjAmMRMRZ87x5ROyWtHvdsotSj/9P0ivneQ8AwHw2alH8bSFRoBY4y0PeyLFq2qhQ3BUR3ygkEiyNVqslSXrggQdKjgTLqtfrqdVq0U1VkI0Gsz8xemD7oznHghobDAZcForckWPl2KhQpL/HcFKegaB+su5MxlU/WJRut6tGozF2Xa/XI88KsFGhiAmPAQArYqMximfb/o6GLYsjksdKnkdE/Giu0QE1wBkt8pae26nT6RT+vaXMQhER49t7AB6H/nMUodfraTAYTOyOy8O038wGMKV+v8+3iJG7InOMQgEAyEShAABkolBgbo1GQ41Ggy8/IXejHOPeJ8Xa1P0oAMxmVEC5Qgp5yfMkjRYFACAThQIAkImupxT615E3cgx1RIsCc+l2u2o2m2WHUSnNZrPQL0OtAvbp4/GFOwBAZVAoAACZKBSYW7vdpvsJuSLHykWhABaMDzXkrdFoFJpjFAoAQCYKBTCDTqfDt6uRqyrmGIUCC8E8T8gbOVYeCgUWbqM++k6no1arVbmzJtTHNBMDtlottVqtgiJabhQK5KqKzWgsl16vR47ljCk8MDe6A5A3cqxctCiwcNwvAEUgz4pDoUAuut3uyp0FttvtRz+4GHgtxioWi9EYYJE5RqEAFmQViyOKV0ZxLKVQ2D7a9tW2v5L8e9SE7Qa2e8nPrqLjBBiMR97qkGNlDWZfIOnzEXGx7QuS528ds91DEbFa7col1+l01Ov1yg5janWKFUNV/9AdZ3TlVlVbpGV1PW2X9MHk8QclvaykOAAAGyirRfHjEXGPJEXEPbafNmG7J9neK+kRSRdHxCfGbWR7h6QdkrR169Y84sWC9fv9Wp2t1/EsddX1+/2yQ5hJlb8PkluhsP05SU8fs+rCGX7N1og4aPskSdfYviUivrp+o4jYKWmnJK2trcWmAgYAjJVboYiIF01aZ/ubto9NWhPHSrp3wu84mPy733ZX0nMkPa5QAADyU1bX0y5Jr5V0cfLvJ9dvkFwJ9d2IeNj2FkkvkPTXhUYJTKGqA5BYHmXnWFmF4mJJV9g+V9I3JL1SkmyvSTovIl4v6eclvcf2DzQcdL84Im4vKV5gamUf1Fh+RedYKYUiIu6T9MIxy/dKen3y+FpJzyw4NADAOnwzG0urDl9kQv2tQp5RKAAAmSgUAIBMFAogsQpdCChfHfOMQoFS1fGgQf2QZ/OhUKByOKiRN3JsNtwKFaUZzfc0mlt/0oE7Wj7tteN8AGBkMBg8Lh9GMxin7+kwa46lX7MKaFFgJXFGibwtU45RKADV7z4ZqKe65hmFAoUb3fMXyFOz2VSj0Sg7jKVAoQAm6HQ6tbunAeqn1+tVPs8YzEahRoOFrVar3ECm0Ov1NBgMyg4DM6pTjtXlZIQWBUqxme6nZRocRP7IscWhUKAU3W73MZcnTqPX69VyIBDl2EyOSdW+JWlZ6HoCZsC9JpC3KuYYhQK5qlLST/OlqtHZ5GibRqOxqbNSFKduOTb6omlas9msdJ7R9YRao08ZRVj1PKNQAGPU5WoU1Ftdrqyj6wlLZ9SsX2RTvkrdG6iGRedZlXOMFgWQodlsVvoARv3VYRyMFgVK1W63+SBG7six+dCiAABkokWByuHsD3mjJTsbWhRYOXWd6hn1skwzCdCiQGmmPaNbloMNxZul1UCeTUaLAgCQiRYFKqPIPmP6qFdXUf/vVb/kdRYUCiyteT8Qqj7/Dqphnjxrt9u16PKi6wkAkIlCAQDIVEqhsP1K27fZ/oHttYztzrB9p+19ti8oMkZUT7vdnrorqN/vr/Rsn9i8acev+v3+ykwcWdYYxa2SflvSeyZtYLsh6RJJL5Z0QNINtndFxO3FhIi6mtTvu9EkboxHYBbjbrM6zXhDHce+SmlRRMQdEXHnBps9V9K+iNgfEd+TdJmk7flHh6qZpSUBbMZmb5u6Kqp81dNxku5OPT8g6XnjNrS9Q9IOSdq6dWv+kaEyuMQVRVj1PMutUNj+nKSnj1l1YUR8cppfMWZZjNswInZK2ilJa2trY7cBZj1j5LsWmNVmWiV1yLPcCkVEvGjOX3FA0gmp58dLOjjn7wTGGh2oDIAjT91ut5Y5VuXLY2+QdLLtE20/UdLZknaVHBMArJyyLo99ue0Dkn5Z0qdtX5Us/wnbuyUpIh6RdL6kqyTdIemKiLitjHgBYJWVMpgdER+X9PExyw9KOjP1fLek3QWGhgoroh+XLqjVVsR4QdXHI8apctcTAKACqnx5LFCaOp71oX7qkmcUCqycuhycqLdlyjO6ngAAmWhRYOnU9Vp11MsqTflBiwIAkIkWBSpvmfp6UV3k2WS0KAAAmSgUAIBMFAoAQCbGKLCU6G9G3lYpx2hRAAAyUSgAAJkoFACATBQKAEAmCgUAIBNXPWFlrNJVKijHsuYYLQoAQCYKBQAgE4UCAJCJQgEAyEShAABkolAAADJRKAAAmSgUAIBMFAoAQCZHRNkxLJTtQ5K+vsmXb5H0rQWGsyjENbuqxkZcs6tqbMsW109GxDHjVixdoZiH7b0RsVZ2HOsR1+yqGhtxza6qsa1SXHQ9AQAyUSgAAJkoFI+1s+wAJiCu2VU1NuKaXVVjW5m4GKMAAGSiRQEAyEShAABkWolCYfsM23fa3mf7gjHrD7d9ebL+etvbUuveliy/0/ZLCo7rLbZvt32z7c/b/snUuoHtXvKza5FxTRnbObYPpWJ4fWrda21/Jfl5bcFx/V0qpi/bfiC1Lrd9ZvtS2/favnXCett+VxL3zbZPTa3Lc39tFNerk3hutn2t7Wen1n3N9i3J/tq7yLimjK1j+8HU/9lFqXWZeZBzXH+WiunWJK+OTtblts9sn2D7C7bvsH2b7T8Zs00+eRYRS/0jqSHpq5JOkvRESTdJesa6bf5Q0ruTx2dLujx5/Ixk+8MlnZj8nkaBcf2apCcnj980iit53i95n50j6R/HvPZoSfuTf49KHh9VVFzrtn+zpEsL2me/IulUSbdOWH+mpM9IsqRfknR93vtryrieP3o/SS8dxZU8/5qkLSXus46kf503DxYd17ptf0vSNUXsM0nHSjo1eXykpC+POS5zybNVaFE8V9K+iNgfEd+TdJmk7eu22S7pg8njKyW90LaT5ZdFxMMRcZekfcnvKySuiPhCRHw3eXqdpOMX9N5zx5bhJZKujoj7I+Lbkq6WdEZJcb1K0kcW9N6ZIuKLku7P2GS7pA/F0HWSWraPVb77a8O4IuLa5H2lYnNsmn02yTz5uei4isyxeyLiS8nj/5V0h6Tj1m2WS56tQqE4TtLdqecH9Pid++g2EfGIpAclPXXK1+YZV9q5Gp4pjDzJ9l7b19l+2YJimjW2VyTN2yttnzDja/OMS0k33YmSrkktznOfbWRS7Hnur1mtz7GQ9FnbN9reUVJMv2z7JtufsX1KsqwS+8z2kzX8sP1oanEh+8zD7vHnSLp+3apc8uywzQRZMx6zbP01wZO2mea1mzX177b9+5LWJP1qavHWiDho+yRJ19i+JSK+WmBsn5L0kYh42PZ5GrbIfn3K1+YZ18jZkq6MiEFqWZ77bCNl5NjUbP+ahoXi9NTiFyT762mSrrb9X8nZdlG+pOH8Q33bZ0r6hKSTVZF9pmG3039ERLr1kfs+s93UsDj9aUR8Z/3qMS+ZO89WoUVxQNIJqefHSzo4aRvbh0n6MQ2bntO8Ns+4ZPtFki6UdFZEPDxaHhEHk3/3S+pqeHaxKBvGFhH3peJ5r6RfnPa1ecaVcrbWdQnkvM82Min2PPfXVGw/S9L7JG2PiPtGy1P7615JH9fiul2nEhHfiYh+8ni3pCfY3qIK7LNEVo7lss9sP0HDIvHhiPjYmE3yybM8Bl2q9KNhq2m/ht0Qo4GvU9Zt80d67GD2FcnjU/TYwez9Wtxg9jRxPUfDQbuT1y0/StLhyeMtkr6ixQ7mTRPbsanHL5d0Xfxw0OyuJMajksdHFxVXst3Pajio6KL2WfJ7t2nywOxv6rGDjP+Z9/6aMq6tGo69PX/d8qdIOjL1+FpJZywyrilie/ro/1DDD9xvJPtvqjzIK65k/ehk8ilF7bPkb/+QpL/P2CaXPFvof3pVfzS8EuDLGn7oXpgse4eGZ+mS9CRJ/5IcMP8p6aTUay9MXnenpJcWHNfnJH1TUi/52ZUsf76kW5ID5BZJ55awz94p6bYkhi9I+rnUa/8g2Zf7JL2uyLiS52+XdPG61+W6zzQ8s7xH0vc1PHs7V9J5ks5L1lvSJUnct0haK2h/bRTX+yR9O5Vje5PlJyX76qbk//nCHHJso9jOT+XYdUoVs3F5UFRcyTbnaHihS/p1ue4zDbsFQ9LNqf+vM4vIM6bwAABkWoUxCgDAHCgUAIBMFAoAQCYKBQAgE4UCAJCJQgFsQmom2ltt/4vtJ9veNmnGUaDOKBTA5jwUEe2I+AVJ39PwWnZgKVEogPn9u6SfTh43bL83uV/AZ20fIUm232D7hmSCu48mE8rJ9iuTVslNtr+YLGvY/ptk+5ttv7GcPwsYolAAc0jmBnupht+ClYaT1l0SEadIekDSK5LlH4uI0yLi2RpOD31usvwiSS9Jlp+VLDtX0oMRcZqk0yS9wfaJ+f81wHgUCmBzjrDdk7RXwzmI3p8svysiesnjGzWcM0iSfsH2v9u+RdKrNZxHTJL+Q9IHbL9BwxvySNJvSHpN8vuv13DK+5Pz/GOALKswzTiQh4ciop1eMLzXlR5OLRpIOiJ5/AFJL4uIm2yfo+Hd2xQR59l+noaTufVstzWcr+fNEXFVnn8AMC1aFEAxjpR0TzJN9KtHC23/VERcHxEXSfqWhlNBXyXpTcm2sv0ztp9SRtCARIsCKMqfa9iN9HUNxzOOTJb/je3RzXg+r+HMozdr2GX1peSWvIckFX1HPuBRzB4LAMhE1xMAIBOFAgCQiUIBAMhEoQAAZKJQAAAyUSgAAJkoFACATP8PWviHiV/xYkIAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"rand = np.random.RandomState(42)\n",
"t = 100 * rand.rand(100)\n",
"y = np.sin(2 * np.pi * t) + 0.1 * rand.randn(100)\n",
"dy = 0.1 * (1 + rand.rand(100))\n",
"\n",
"\n",
"\n",
"# \"this function takes times t, mags y, and errors dy, and a period and phase folds the lightcurve at this period\n",
"\n",
"def phase_fold(t,y,dy,period):\n",
" phases=np.remainder(t,period)/period\n",
" phases=np.concatenate((phases,phases+1))\n",
" y=np.concatenate((y,y))\n",
" dy=np.concatenate((dy,dy))\n",
" plt.errorbar(phases,y,dy,ls='none',c='k')\n",
" plt.xlabel('Phase')\n",
" plt.ylabel('Flux')\n",
" \n",
"phase_fold(t,y,dy,1)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Phase-folding can be used as a rudimentary method of period determination. One may make educated guesses about the period of the phenomenon and phase-fold the lightcurve to see if a repeating pattern emerges. Try this for this observations of the asteroid 1856 Ruzena, observed by non-Astronomy majors from the University of Maryland in April 2018 along with collaborators in Malta. Asteroids typically have rotation periods between 2 and 12 hours.\n",
"\n",
"The code below plots the raw, or unphased lightcurve, for all four nights of observations of the asteroid."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'Differential Magnitude')"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2de5gddX3/X+8sB7MBZKPk8QkLMWgl1BCTlSBIWn+C2ig3Vy5FKirqT6qtChS3Df3RJiiWPEbFPvrrBRCxwg+5dksINdpyq1TADZsAkaReEMgGJZSsXLIlm83n98eZszl7MjNndvfMObMzn9fznCfnzPnOzGcnc+bz/X6uMjMcx3Ecp5ZprRbAcRzHySauIBzHcZxQXEE4juM4obiCcBzHcUJxBeE4juOEsk+rBWgkBx10kM2dO7fVYjiO40wZ1q1b95yZzQr7LlcKYu7cufT19bVaDMdxnCmDpCejvnMTk+M4jhOKKwjHcRwnFFcQjuM4TiiuIBzHcZxQXEE4juM4oeQqislxHKcI9PYPsGrtZrYODnFwRzs9S+fR3dXZ8PO4gnAcx5lC9PYPcPFtjzI0PALAwOAQF9/2KEDDlYSbmBzHcaYQq9ZuHlUOFYaGR1i1dnPDz+UKwnEcZwqxdXBoXNsng5uYHMfJJM2ys081Du5oZyBEGRzc0d7wc6W2gpB0jaRnJT1Ws/2zkjZL2ijpyyH7HSrpbkmPB2POT0tGx3GyScXOPjA4hLHHzt7bP9Bq0VpOz9J5tJfaxmxrL7XRs3Rew8+VponpWuC91RskHQ+8H3iLmc0HvhKy3y7gIjP7XeBY4E8lvTlFOR3HyRjNtLNPNbq7Orn8tAV0drQjoLOjnctPWzC1opjM7D5Jc2s2fxpYaWavBGOeDdnvGeCZ4P2Lkh4HOoGfpiWr4zjZIsyEAunY2aci3V2dTTG3NdtJfTjw+5IelHSvpKPjBgcKpgt4MGbMeZL6JPVt27atocI6jtN84sxIB7aXmiiJ02wFsQ8wk7LpqAe4SZLCBkraH7gVuMDMXog6oJldaWaLzWzxrFmhJc0dx5lCxJmRwp8WTlo0W0FsAW6zMg8Bu4GDagdJKlFWDteb2W1NltFxnBYSZ0Ya3DHcREmcZiuIXuAEAEmHA/sCz1UPCFYU3wIeN7OvNVk+x3FaTMeMaDNS3HdO40nNSS3pBuCdwEGStgDLgWuAa4LQ153AR83MJB0MXG1mJwJLgA8Dj0paHxzuL83szrRkdZxG4vH7k8NsYt85jSfNKKazI746J2TsVuDE4P2PALc0OlOSZtbJySu/HYo2I8V95zQeL7XhOA3E4/cnT1xGcBrZwk40riAcp4E0s07OVKa3f4AlK+/isGVrWLLyrjGhrccfER6N2DZNqWQLO9F4LSbHaSDjqZNTVF9FPTPc3ZvC85kOeNU+hbg+kJ17w1cQjtNAktbJKXKtoXpmuKjVVlH8D1m6N1xBOE4DSVonp8i+inplNKL8DEXJos7SveEmJsdpMEnq5BTVV9HbP4CAsGjVimLoWTqPnps3MLx77KiXd+6it38g92amLN0bvoJwnBYQNUvOe5TOqrWbQ5WDYNQM193Vyb777P1oGh6xQqywopIBW3FvuIJwnBSIi9KB5tb0zxJRs2BjT55Ib/8AL+8cCR1XhBXWS/+za6/t00RL7g03MTlOg0mSLFf5NwuRKs0kKsqrs2p2HLdKKMIKq9a0BhCyqSmRTr6CcJwGk8TJmJUwxmYTleNQvT1ulVDUFRaw1/3TjEgnVxCO02DqORmzFMbYbKJyHKq3R60SOtpLuVeicSuk6vuqWZFOriAcp8HUc0BnKYyx2SSJ0Inyz6w4dX6qsmWBua9NVmYkKlQ4avtEcQXhOA2mngM6S2GMzSYql2F6ac+jqJk9l7NEb/8A9//i+cjvq81rbRGdk6K2TxR3UjtOg6nngD6wvcRgSFZwERLBop5fQ8O7x+Q4NKvncpaot4Ksvh4jEXXPo7ZPFFcQjpMCcQ+4qIdkEdppxnWEW7V2c+GUQjVJV5BxvqpGryDcxOQ4TWZ7xEMyanueiOsIVwQTWxxxDuqOqtVl3Eqj0SsIVxCO02TiZnl5jmSKSgKrUP2ArJdomEeiQoABTl44e/R9nCLtbHCeiCsIx2kycbO8PEcyRSWBwVgnflHDgKNCgAFufOjp0b8/zlfV6DwRVxCO02RmFtTMEve3nX5U5xjnfhHDgOOuz/BuY8XtG4FoX9WM0rSG+3DcSe04Daa3f4AVt28cjVSaOaPE8lPmj/5448zEcTb6qU5UmQ0YO3suahhw3PUBRu+nKF/VjuHdDZfJVxCO00B6+wfouXnDmDDW7TuG6bllw6iJIK7xTYN9jJkizvxR/fAvaqXbpOahKB9WGkFwriAcp4FE2dmrS1XHPejC8iPyQndXZ6R5rfqaFLXSbT3zUOXaRfmwjMYHObiCcJwGEmcGqXwX96BrdBx71jjpLbP3munWPvyLmkkN8f//y08plxqJi1TyWkyOk2HiVgeV7+IedI2OY88Svf0D3LpuYEzDIDHWQV2hu6uT+5edwBMrT+L+ZScUQjkAnH3MoaHbl7zxNaPXIKmprhG4gnCcBtKzdB7TQiaBpTaN+WFHzQIbHceeJcKik4y9wzuLmANR4bLuBZxz7JzRlUSbxDnHzuH6T759dExSU10j8Cgmx2kwbRK7a1YCZx196JhZcM/SeWOaCkH+7exRs9uBwaHROkxJmi3lncu6F3BZ94LYMctPmd+U+8cVhOM0kCgnde0suYgd5eLCOHtu2QBE50BcdFP5+zxfn/HQrPtHliOb5+LFi62vr6/VYjgF5rBlawj7RQl4YuVJzRYnU/T2D3DhjetDrw+UzWtbg+zpMNpLbYVxVjcTSevMbHHYd+6DcJwGUtQY/iR0d3VGPvyB0ZlwFEXIpq6QFT+MKwjHaSBFjeFPSpwTvmImqb1+1eQ9mxqyVYvKFYTjNJAix/AnoWfpPEohYV6VKK/K9YvKByjCSixLtajcSe04DaaI3dCSUrkucbWqKv8WLcqrQr1aVL39A00LbnAFEUIz/wMcp2gkUaBFjPKqEBXtdXBHe2gY8IU3rqfvyefrhsZOhMQKQtJ+ZvZywyXIGB6H7TjZoKgrsbgcmahkw+sfeIrFr39Nw69XXR+EpOMk/RR4PPi8UNLfNVSKDJEl+5/jOMUjzo8VZX4y0mk2lcRJfQWwFPhvADPbALyj3k6SrpH0rKTHarZ/VtJmSRslfTli3/cGY34uaVkCGRtGUWvRO04zyUoYZ1bp7uqkZ+k8Dg5yQ1at3Uxv/0Cskz6NZ1SiKCYze7pm00jowLFcC7y3eoOk44H3A28xs/nAV2p3ktQG/F/gfcCbgbMlvTmJnI3A49idRuAPwGiyFMaZVaKu0fFHzIrs+5DGMyqJgnha0nGASdpX0ucJzE1xmNl9wPM1mz8NrDSzV4Ixz4bs+jbg52b2SzPbCXyPslJpCh7H7kyWJA/AIiuQembcIl+bClHX6O5N2/jQsXPqlkxvFEkUxKeAPwU6gS3AouDzRDgc+H1JD0q6V9LRIWM6geoVy5ZgWyiSzpPUJ6lv27bopt9J8Th2Z7IkeQAWeQYdZ8Yt+rWpKMeomlVbB4e4rHsBV5y1aDTpsE0avb8afZ3qRjGZ2XPAhxp4vpnAscDRwE2S3mBjC0KFraAiM/TN7ErgSijXYmqEkEWNnnAaQz0/VpwCKcJ9FxfGWeRrUxtBGUZtT5G0Iy4jFYSkbxD/YP7cBM63BbgtUAgPSdoNHARsqxlT3TXjEGDrBM41KTwXwpkoUQ/AA9vLNfyLHghx/BGzuO6Bp0K3Xx+yHYpxbcKUYzW1ZqRmKNM4E1MfsA6YDrwV+FnwWkQyJ3UYvcAJAJIOB/YFnqsZ8xPgTZIOk7Qv8EHg9gmeb2JCFnyZ60yOqHISL+/cFRuJUpRAiNrS59Xbi3xt4pRgmKm7GRONSAVhZt8xs+8AbwKON7NvmNk3gHdRVhKxSLoB+DEwT9IWSZ8ArgHeEIS+fg/4qJmZpIMl3RmcdxfwGWAtZWf4TWa2cXJ/5vjwXAhnMnR3dbL/9L0X58Mjxqq1myML0u0IFEjeiXuwFTlIJEoJdna0h7ZdbYYyTeKkPhg4oOrz/sG2WMzsbDObbWYlMzvEzL5lZjvN7BwzO9LM3mpmdwVjt5rZiVX73mlmh5vZG83sS+P9oyZL0U0AzuQZ3DEcun3r4NBoIERH+9i2kdt3DBdipRr3YKsNEuloLzG9NI0Lb1yf+4im8SrHZijTJApiJdAv6VpJ1wIPA3/TMAkySEeT+r06+aXe7K67q5P9XrX3KqMIK9V6D7burk7uX3YCV5y1iFd27Wb7juFCmHrHG0HZjIjLJFFM35b0r8AxwaZlZvbrhkmQMXr7B3jpf3bttb226bzjxJGk53RRV6pJC/EVMaJpvBGUaUdc1lUQkiplNbYH/x4u6fAgES53XLp6Y2hPYcxye1M6jSfuIViJkIsKESzCSjXJg62oCrQe1RGWB7aXkMomzTSiLZNUc+2pej+dcqbzOoJopDzR2z/A9gjb8fBuuKT30VRK6jr5JOwhWC/WvSgO2STE5UsUldr7p9JTA9LJg6jrgzCzU6pe7wGOBH7TkLNnjHq23xserC1J5TjjIy7W3bP2x1LkiKYo6uVKNNqHNZGGQVsoK4ncUW/pOmINSdR2CkzUPSbg/mW5W5RPiiI3DYoiiXmtkSa4JD6I6ozqaZRzIDY0TIIMEbWkrRDVJ9dxklJ0s8l4KxR42Zux1HtGVcY0iiQriL6q97uAG8zs/oZJkCF6ls6j5+YN4U5q4OxjDg3d7jhhhD0Mk0Q35ZWk3Rq9zE00UWVKKrQiD6KjklVtZteb2f2Szm+YBBmiu6uTVWcupL009rIIOOfYOe6gdhITVa4FKGy14CQVCrzMTTxRZUognXspiYL4aMi2cxsmQSYZa0qaXmpj8etf0yJZnKmIl2vZmyRhq37d4qnnw2paT2pJZ0taDRwm6faq190E7UfziN+gTiOI+iFXZsRFnCFH2cYNRstoeO5DPM0uZhi3gvhP4KvApuDfyusialqJ5gm/QZ1GEPWDrTR3qaYoE5CoIoWwR1EWvcxNvW56zQ79javm+qSZ3WNmbzeze6teDwcVV3NJkcsNO40j6occFSo9MDiU+2J01bWDwhgaHsGMwuY+JPG/NLvjZZyJ6UfBvy9KeqHq9aKkF1KRJgOE1fIvTfM6TM74iPohRz0coRjmpkohvqiA8d8ODbsTv4qw1WXlGj6x8qRU/A7VRIa5mtnvBf8eEDUmj/Q9+fzeYa6e/uBMgKgY/rhSG3kvRlchLh+kqLkPWTRvJ4liQlJb0NRnTuWVtmCtoLd/ILTlYaXRi+NMlnpmFiiGv8vLaOxNFs3bdRWEpM9Srr30Q2BN8LojZblaQlyFzSL8aJ3mUDERRCmJIvi7mm1LnwpkUWkmyaQ+H5hnZrkNba0QpwSK8KN1mkuRs6rBy2jUksXaU0kUxNPAb9MWJAtE2UUFhfnROs0jiw8Ep7VkTWkmURC/BO6RtAZ4pbLRzL6WmlQtImxGJ+BDx87J1H+akx+y9kBwnGqSKIingte+wSu3+IzOcRxnD7Ic9ThYvHix9fX11R/oOCnhlUgnj1/DeBp9fSStM7PFYd8l6QexGvYK7vkt5TLg/2hm/zNhyRwnRyQtZ+2UCXvQAX4NY2j2PZYkD+KXwEvAVcHrBcphr4cHnx3HwQs9joeoshKXrt7o1zCGZt9jSXwQXWb2jqrPqyXdZ2bvkLQxFakcZwqSxUzYrBL1oIvKMPdrWKbZ91iSFcSs6szp4P1BwcedqUjlOFOQLGbCZpXxPtD8GpbJUrnvChcBP5J0t6R7gP8AeiTtB3wnFakcZwqSxUzYrBL1QOtoL/k1jKHZ91hdE5OZ3SnpTcARlNMCNlU5pr+eilSOMwXxMOn6VBzTA4NDiLHRL+2lNlacOh/waxhFs++xRGGuko4E3gxMr2wzs39KRaJJ4GGujpNdaiNwgFEl0emKoGVMNsx1OfBOygriTuB9wI+AzCkIx5lqFCnmP8wxXVEO9y87oTVCZYgs3gtJopjOABYC/Wb2MUmvA65OVyzHyT9Fy5vwKK9oJnovpK1Ukjiph8xsN7BL0quBZ4E3NEwCxykoRcub8CivaCZyL/T2D9Bz84YxuSQ9N29oaEfCJAqiT1IH5aS4dcDDwEMNk8BxckS9pvPVFG1G7VFe0UzkXlhx+8a9ul8O7zZW3N649LQkUUx/Erz9B0nfB15tZo80TALHyQnjNRPEtd3MIx7lFc1E7oXBoeFxbZ8IkQpC0lvjvjOzhxsmhePkgDgzQdhDsIgNg7y8eThZvRfiVhB9wEZgW/BZVd8ZEBt2IOka4GTgWTM7Mti2Avhk1TH/0szuDNn3QuB/B+d5FPiYFwV0ss54zQQ+o96bS3of5YYHn2bEjDaJs485lMu6F7RarNTp7uqk78nnx/ztpx8VrkwrjukoZs4oNUyuOAVxEXA6MAR8D/hnM3tpHMe+Fvgme4fDXmFmX4naSVIn8DngzWY2JOkm4IPB8Rwns0zETOAz6j1c0vso1z3w1OjnEbPRz3lXEr39A9y6boCRIC9txIwbH3qaNY88w+CO4chqt7WU2sTyU+Y3TK5IJ7WZXWFmvwd8BjgU+HdJN0lalOTAZnYf8PwE5doHaJe0DzAD2DrB4zhO00jqhB2PI7tI3PDg0+PanifCzJPDu43tO4brVrut0NnRzqozFjZ0wpHESf2EpH8B2oEPUy7zvX4S5/yMpI9QNmFdZGbba843IOkrlLvYDQE/MLMfRB1M0nnAeQBz5syJGpaYLCarOFODJCaj3v4Bem7ZwPBIeaY4MDhEzy0bxuxfVEYiqjpEbc8TSSLX4qrdClJJNowstSHpDZRNO+8HnqZsZrpjPL4ASXODfSo+iNcBz1H2LXwRmG1mH6/ZZyZwK3AWMAjcDNxiZtfVO99kS22ElQJoL7Vx+WkLCv/jdcZP2GTj0tUb2b5j7yiTmTNK9P/1H7RAytZSfY2i1ECbxC8uP7GpcjWbJSvvCjVPJmUy2ehxpTbi8iB+Dvwh8H3gx8Ac4E8k/ZmkP5uIIGb2GzMbCRLvrgLeFjLs3cATZrbNzIaB24DjJnK+8VK0xCUnPaIa4oQpB4DtO4YLZ26qvUZRnH3MoU2TqVWEmSfDaHa12zgT0xfYU2xx/0acTNJsM3sm+PgB4LGQYU8Bx0qaQdnE9C7K5qjUKVrikpMeUZONOPJeaqOWsGtUTdGimGCPefLA9hIv79w1aoqE1lS7jVQQZrZiMgeWdAPlIn8HSdoCLAfeGTi5DfgV8MfB2IOBq83sRDN7UNItlDO2dwH9wJWTkSUpRUtcctJjopOKuLyJvBF1jQQ8sfKk5gqTAWoj2uL8oc26P5IU65sQZnZ2yOZvRYzdCpxY9Xk5ZYXSVHqWzhvjQIRy2Firk1WcqUH1D3qaFOpc7Wgv8fIru/YqkVBNUVasPiGLJwsh0ElqMRWL2t9t/gMonAZQa0+Pirw5eeFsVp25kM6Yh2BRHpBemyn7uIKoYtXazaHFr9xJ7dSjnj29wt2bttHd1cn9y07g62ctKvQDsrurk8tPW0BnRzuiHInjEYNjaXXOTFwtpthIJTP7WuPFaS3upHYmStJ7pHpc5UG44vaNowXWppeKNWfLghklq4QVf+y5eQOXrt44Jrs6zesXdzceUOeVO7xevTNRkt4jYeNe2bV79P32HcNcfNujhQp3dcJJml2d5r0SF8V0aWpnzShZrajoZJ+we6eWsHtpvBVg84RXLYgnaXZ1mvdKkp7U04FPAPOB6ZXttRnQecCrazoTJezeOf6IWdyx4ZlY81GRzJrVCqE2zr9oOSBJiIryqiXNeyVJmOt3gU3AUsrJcx8CHk9NohbjNlFnooTFsd+6bs/yv2I+qoyF4oR61trTw5raFGXllJQkq1JI915JoiB+x8zOlPR+M/uOpP8HrE1NIsfJCVHmoxW3bxwzky61aa+M2byZNZNGeeVx5ZSUMJPb5actqJtdnea9kkRBVFT9oKQjgV8Dc1OTyHFyQtTDbnBoeHQGPTg0TGmamDmj1LTIlFaQ9MGft5VTUqLa1V5+2oIxRfia7bdJoiCuDCqs/hVwO+W6TH+dmkSOkxOS2pArkSmdOVUOkOxa5HHllJSkwQrNNoEn6QdxdfD2XuAN6YrjOFOX2tnd8UfM4tZ1A4lMK5BvR22YPb00Tew/fZ9cr5ySktVghbhEuXPM7LqohLk8Jso5zkQJMxHcum6A04/q5O5N20aVxo6duyJLfkN+HbUeIRhPVoMV4lYQ+wX/hiXFeYUix6kiykRw96Zte9mQ60WmtHrWmBYeIRhNVnOw4hLl/jF4+29mdn/1d5KWpCqV40wxkpoIqmfSUTb5Vs8a0yKJg7WoyXNZXWElcVJ/A3hrgm2OU1iSmgiqH4AdLQhbbBVRUTqw5+GYZEyeSbLCykwUk6S3U271OavGD/FqoH5vPMcpEFEmguOPmMWSlXeNxrG/+MouRoKKwYNDw0wTuQ9xhWRROkUuOxJFq7PP41YQ+1IOad2HsX6IF4AzUpHGcaYoUaU2qqOYwrKHK9Xl895BLcoENzA4NKpAoxybefXJ1CML2edxPoh7gXslXWtmT6ZydsfJEbUmgiUr70oU4hoX1ZQXokxwgrr5EXn1ydQjC9nnSYrPv0rSlZJ+IOmuyis1iRwnJxR15htGWPc4UT8cMq8+mSRkIfs8iZP6ZuAfgKuBZBk/juMkzqTuaC81QZrWEmaCi7s2glz7ZJKQ5P4ptanltZh2mdnfpyaB4+SU44+YxXUPPBU7pjRNrDh1fpMkai1hJriwB2BnR/uY3JGiMKFM/JQz0pKYmFZL+hNJsyW9pvJKVyzHmfrcvWlb6PY2abQH86ozFxZ2hhxmdiqqSanikB4InPXVmfiVnt1t0l77De82Vq3dnJpcSVYQHw3+7anaZnhdJseJJcqGvNss91FLSchqclgrSJKJf9iyNaH7trRhkJkdltrZHSfHRNmQp0kctmxNoR+IFbz8RpkkmfitqNdU18QkaYakSyRdGXx+k6STU5PIcXJCmAkFYMSsaU3nnalB1EO+ensrTHJJfBDfBnZSzqoG2AJclppEjjNF6e0fYMnKuzhs2RqWrCxHgl9+2oJYG3Il0ckpNlEP/0om/mHL1rBq7WbeOufA0fuoTeL0o9Jdgcks3g0uqc/MFkvqN7OuYNsGM1uYmlQTZPHixdbX19dqMZwCElaltb3UxuWnLRj9AR+2bE1k0ElnR3vh7fBQ3GJ9MLEoptp7bCJIWmdmi8O+S+Kk3impnSCgStIbgVcmLI3j5JAkdYSSZBMXrUBdNV6sb/yZ+GmX2khiYloOfB84VNL1wL8Df56KNI4zRUniZEyaTZxXs1OtCa7W9xKnZItI0uiklkUxSRKwCTgNOJby/Xy+mT2XmkSOMwVJEmEynmzivJXpSLI6yGrbzVaRNBO/ZVFMVnZQ9JrZf5vZGjO7w5WD4+xN0giT7q5O7l92Ak+sPIn7l51AZ4LolTyQZHWQJJKnSBx/xKy6Y7IQxfSApKNTk8BxckB3V+eYiKXOjvZEzsOiZBNP1ASXx2uRlCSZ+JN1UNcjiZP6eOBTkn4FvExgNjWzt6QmleNMQSaS9FWUbOIoc8mB7aXRfhAHd7Rz+lGd3L1pW66vRVKykImfREG8L3UpHCcHTKbnct4fgmEd90rTxMs7d402wqnUH0p7VjxVaEXmdC11TUxBs6BDgROC9zuS7Oc4RSKs2FptlnSSMXklzAS3//R9xvTjhmJHLdWSBZNbklIby4G/AC4ONpWA6xLsd42kZyU9VrVthaQBSeuD14kR+3ZIukXSJkmPB/2xHSezJHHCFj2Ms9ZBPxjRSa+oUUu1TNSv1UiSmJg+AHQBDwOY2VZJB8TvAsC1wDeBf6rZfoWZfaXOvn8LfN/MzpC0LzAjwfkcp2UkccJ6GOdYsmBCyTqtNj8myqQ2M5NUyaTeL8mBzew+SXPHK5CkVwPvAM4NjrOTci0ox8ksSZyw0yRGQkrbFPWB2LN0Hj03b2B4955rUpqWboe0rJLVEiNJfAk3SfpHoEPSJ4F/A66axDk/I+mRwAQ1M+T7NwDbgG9L6pd0dZxSknSepD5Jfdu2hYeFOU7ahNmLK07Yis8hTDkUOYwTKMdEVjG827jgxvXMXbaGri/8oBD+mSz7piKL9Ul6lZm9Erx/D/AHlP8715rZDxMdvLyCuMPMjgw+vw54jnJ1gS8Cs83s4zX7LAYeAJaY2YOS/hZ4wcz+qt75vFif00pqZ4E7du5ie4idvU1it1mmZorNIOn1qabUJladke+ue1GtVzvaS6xf/gepn3+ixfp+DLxV0nfN7MNAIqUQh5n9pkqoq4A7QoZtAbaY2YPB51uAZZM9t+OkTa29OKoDWBE7yoWV2kjC8IilWowuC0T5oAaHhuntH8isD2JfSR8FjpN0Wu2XZnbbeE8mabaZPRN8/ADwWO0YM/u1pKclzTOzzcC7gJ+O91yO02q8o9wewiK4kpJ3J35czaVKhNuqtZsZGByiLfBjdTbp3onzQXyKcoG+DuCUmlfdjnKSbqC8CpknaYukTwBflvSopEcoZ2hfGIw9WNKdVbt/Frg+GLcI+Jtx/2WO02K8o9weJvOQz7sTP84HVblHKgqk4sdq1r0Tt4KYbWafDhoFXTneA5vZ2SGbvxUxditwYtXn9UCoTcxxpgq1ZTTCopjSruefFZJWJq2l1Jb/qKburk4uXb0x0l8VtfJqxr0Tt4KoJMZ9KrWzO07OqU4O2x0REJJ3EwpE98KIY+aMUu4d1BWWnzI/NGs6LPKtmrTvnbgVxH9Luhs4TNLttV+a2anpieU4+aPIiWFhRQnDWmo2ooXmVCSqaGPF9xBF2vdOnII4CXgr8F3gq6lK4Tg5oF6yU1jBuqLkQURdm8Wvf00mE8RaQVTWdO09U6EZ5rdIBRFkMD8g6Tgz8ww0x4khrmMa7JkZHtheYnppGoM7hgvzQKzXTa7y91eUyIU3ri/MtalHd+VyR8sAABDsSURBVFcnfU8+z/UPPLVXa9q9NzSeSAUh6etmdgFwTaXMRjVuYnKcPUQV4ltx+0Ze2bV79LvBoWHaS21ccdaiwjz84ooUViuHei1Ji0hv/wC3rhsI1QXDu9PPEYkzMX03+LdeYT3HKTxxyU61FCVyqUKSIoVJlEgRqZc/0jIntZmtC/69V9Ks4L2bmhwnhPGGcRYhcqlCEue8V7oNp97fn7aTOjLMVWVWSHoO2AT8l6Rtkv46VYkcZwoS1dxl5oxS6PgiRC5VqNf4prd/gGkKD3ot0nUKI+7vb0aAQ1wexAXAEuBoM3utmc0EjgGWSLowVakcZ4oR1dwlKr69CJFLFeIa31R8D2Hx/m0FLf1dTVQ2/swZpaaEA8f5ID4CvMfMnqtsMLNfSjoH+AFwRaqSOc4UI665i4dyhhNnYx/ZbfQ9+XwhrlVcr/K+J5/nhgefHqNEZ+ybpJXP5Ik7S6laOVQws22SwtfNjlNQ4nIgWt0VrNXERSjVs7Hf8ODTXNa9IHUZW0m9EOlb1w3stcJqVpRXnIKI6+LmHd4cJ6BeiGZWu4U1i7gIpXrO/XqlJvJAvV7lWa3FtFDSCyGvF4F8q3THGQdxP/De/gF6btkwplvYBTeuL0y3NIiPUIqysVdoi3Be54Xe/oFIBbl1cKjuCquVYa7R/2s5pegzPWdixD0AL129keGRvWfB23cMFyYRLC7MtfK3X3zbIwwN795rzNnHHJq6fK2isvKMomNGCbPwXJrqMWmSpCd1IchyX1gn20SFIh7c0R7bUrPajJBn6oW5dnd18vgX38c5x84ZXTG0SZxz7Jxc+x/qJcGZQb0FVNoWuOa4wqcAnsnpTJS4InwX3Lg+dt+J9EiYakRVKq39XV3WvSDXCqGWeuah38asHMYzZjIUXkFUzEpxdkDHiSPuAbji9o2xJoK829grFD2SK4x6DvrKyjSr5b5zT230SRhFz+R0khH1AFxx6nx6bt7A8O5wW0ARonSccMJWnhWqTXBJxqRFoRVEPRtg0TJencZTURoX3bQhVBl0+gQEgEt6Hx1T0nq/fdv40gfy3TioeuU5MDhEW9CStjPEBFdr5WiTOP2o9FdlhXZSx5mPqssBOM5k6O7q5Kt/uLDwJTei+NBVP+a6mn4HL+8c4aKbN3iQCOX7p9bRP2LGresGUr8+hVYQUeajzo527l92gisHp2HE1SMqMr39A9z/i+dDvxsJ+h3klerISdhjbgyLoLx09cbYZLq0KLSCqBd+5zhOuly6emPs93kOEokzcVc//C/pfTQyXLpliXJFIGn4neNMFu+YFk5cngjkO0gkSZZ0b/8A1z/wVOQYj2JKGQ+/c9Kmt38g1EnteTb1yfNqPkmY66q1m2NbT3sUU8pUl9c4sL2ERKEayjvpUqnFFBXOmmcTShIEsQ/APP/+koS5XhiTaCnSvz6FVhC1y/7qhCY3AThJqFe/K6oWU4U8m1CSEKcc8p5EmCTMNS6JtxkZNIVWEPXyINwE4MQR5le44Mb1XLp6I8tPmU93V2ddG/vxR8xqhqiZpTPGzFKEJMJ6Ju64ci0d7em35Sl0FFOS5X3RTQBONFETjEql1iQx6ndv2paGaFOGnqXziFoneBJhWYHst294Ye1mLLAKvYKo5ySqjHGcMOLuncrqs6O9FFuLqegTkEpLzetrEuWKFG7e2z8wpmbXzBml0RUowI6d4VaOwTqr00ZQ6BVEvWYlRbpJnfFTz0Y+MDjEilPnx47xCUi5iusVZy0qZBLhJb2PcsGN68dMIrbvGKbnlj1Z5HHl5NOm0AqiNru1o73EzBmlwt2kzsSoZyNvk+ju6mRmRFMXke8wTiee3v4BrovIcRge2ZNFHuWnaob/qtAmJvA8CGfi1DMfVRTI8lPCK7p+6Ng5fu9R3CTCemUyKubHNY88E/r9mkeeSb1/RqFXEI4zGeo5CaudrLXNNNumicWvf03jhZqCxDXryjP1/E8VE1JUJFy9CLlGUHgF0ds/wJKVd3HYsjUsWXmXV490ElPPSVgxH126eiMjNauHkd1Wtw5RUYjr6Z1n6vkQsmB+LLSC8D7UzmSo9wOvmEdaOQOcCrTSCdtK6imAyv0Tle8wpfMgJF0j6VlJj1VtWyFpQNL64HVizP5tkvol3ZGWjEVd2jqNod4P3CcayShqVeXurk6mRZgpqyPkVpw6n1LNwNI01Y2QawRpriCuBd4bsv0KM1sUvO6M2f984PFUJAso6tLWaQzdXZ3s2xbtiKhMNKJG5LuQRHKK3Cvjj46ZE7r97GMOHX3f3dXJqjMXjrk+q85c2JTrk1oUk5ndJ2nuRPaVdAhwEvAl4M8aKNYYohLlpkn09g8U4gZ1JsfOmDpLlYlG1Ij8F5JITlGjCStRSDc8+DQjZrRJnH3MoXtFJ7Xq+rQizPUzkj4C9AEXmdn2kDFfB/4cOKDewSSdB5wHMGdOuDaOIqqa4ohZIcLsnHSp2NCjwmGbYUN2ss9l3QtSD1edKM12Uv898EZgEfAM8NXaAZJOBp41s3VJDmhmV5rZYjNbPGvW+BJHKkvbsIxY90U4SYhKgoM9PoqocNicFyt1ckBTFYSZ/cbMRsxsN3AV8LaQYUuAUyX9CvgecIKk69KSqburk91eq9+ZIMtPmU8pxA9xTlUSXFQ4bDNq6TjOZGiqgpA0u+rjB4DHaseY2cVmdoiZzQU+CNxlZuekKdeBEUv9qO2OU6G7q5NVZ4x1IH79rEVjTAZFDeN0pj6p+SAk3QC8EzhI0hZgOfBOSYso++d+BfxxMPZg4Goziwx7TZOonhBuAnCSkKSmf62vqwhhnM7UJ80oprNDNn8rYuxWYC/lYGb3APc0VLAaLul9lFd21RZCKOMmAKcRVHcOi+o85zhZpPDF+q6PqKYIbgJwklGv7SgUN4zTmdoUXkHExaIXvR2kU5+iViJ1ikGhazHVo+jtIJ36eLkWJ88UWkHUq5XjYa5OPbxci5NnCq0g6s3y3Afh1MNDWJ08U2gFUW+W5z4Ipx5FrUTqFINCK4h6szz3QTj1KHIlUif/FDqKKapYXwW3IztJ8BBWJ68UegURV6wP3I7sOE6xKbSCgLKSqG7OUY37IBzHKTKFVxAQ7WtwH4TjOEXGFQQey+44jhOGKwg8lt1xHCcMVxB4LLvjOE4YhQ5zreDlmB3HcfbGFUSAx7I7juOMxU1MjuM4TiiuIBzHcZxQXEE4juM4obiCcBzHcUJxBeE4juOEIrO4rsxTC0nbgCdbLUcGOAh4rtVCZBy/RvH49YknT9fn9WYWWnguVwrCKSOpz8wWt1qOLOPXKB6/PvEU5fq4iclxHMcJxRWE4ziOE4oriHxyZasFmAL4NYrHr088hbg+7oNwHMdxQvEVhOM4jhOKKwjHcRwnFFcQOUPSeyVtlvRzSctaLU+WkHSNpGclPdZqWbKKpEMl3S3pcUkbJZ3fapmyhKTpkh6StCG4Ppe2WqY0cR9EjpDUBvwX8B5gC/AT4Gwz+2lLBcsIkt4BvAT8k5kd2Wp5soik2cBsM3tY0gHAOqDb76EykgTsZ2YvSSoBPwLON7MHWixaKvgKIl+8Dfi5mf3SzHYC3wPe32KZMoOZ3Qc832o5soyZPWNmDwfvXwQeB7xRSoCVeSn4WApeuZ1lu4LIF53A01Wft+A/bmeCSJoLdAEPtlaSbCGpTdJ64Fngh2aW2+vjCiJfKGRbbmc3TnpI2h+4FbjAzF5otTxZwsxGzGwRcAjwNkm5NVe6gsgXW4BDqz4fAmxtkSzOFCWwrd8KXG9mt7VanqxiZoPAPcB7WyxKariCyBc/Ad4k6TBJ+wIfBG5vsUzOFCJwwn4LeNzMvtZqebKGpFmSOoL37cC7gU2tlSo9XEHkCDPbBXwGWEvZuXiTmW1srVTZQdINwI+BeZK2SPpEq2XKIEuADwMnSFofvE5stVAZYjZwt6RHKE/Ifmhmd7RYptTwMFfHcRwnFF9BOI7jOKG4gnAcx3FCcQXhOI7jhOIKwnEcxwnFFYTjOE4CJK2StEnSI5L+uRLuGjG2TVK/pDuqtl0r6Ymq6LBFwfYDJa2uKgD4sWD7Ikk/DrY9IumsqmN9Kxj/iKRbgsTGevK/MzjvRkn3JvmbXUE4hUbSS3W+v0fS4uD9nXEPhXGc81xJ24IHyM8krZV0XIL9uiW9ebLnd+oTPEyvrdn8Q+BIM3sL5aKYF8cc4nzKoea19JjZouC1Ptj2p8BPzWwh8E7gq0Ee0w7gI2Y2n3Iy3ter7r8LzWxhIMtTlMPb4/6eDuDvgFOD450ZN76CKwjHSYiZnRhkzzaCG82sy8zeBKwEbpP0u3X26QZcQbQIM/tBkGsE8ADlSgV7IekQ4CTg6qSHBg4IkhT3p1xQcpeZ/ZeZ/Sw491bKtZ9mBZ9fCM4loD04RiWR71ZJPwleS4Jz/BFwm5k9Fez/bBLBXEE4hSeYLVabAr4p6dyQcb+SdFDwvlfSumC5fl7VmJckfSlY/j8g6XX1zm9md1PucXxecIxPBj/uDcGPfUawwjgVWBWYCd4YvL4fyPEfko6Y9MVwkvJx4F8jvvs68OfA7pDvvhSYha6Q9Kpg2zeB36VcFudRyuXDx+wr6W3AvsAvqrZ9G/g1cATwjWDz3wJXmNnRwOnsUVKHAzODFfE6SR9J8ke6gnCcifFxMzsKWAx8TtJrg+37AQ8E5oL7gE8mPN7DlH/oUJ7pHR0c43HgE2b2n5TLplRMFL+grFQ+G8jxecomBGcSSHowqNR6NXBqlb9gadWY/wPsAq4P2f9k4FkzWxdy+Isp/x8fDbwG+Itg+1JgPXAwsAj4pqRXVx1zNvBd4GPVisPMPhbs8zhQ8U+8O9h/PeX75dUq9/XYBziK8spmKfBXkg6vdz32qTfAcZxQPifpA8H7Q4E3Af8N7AQqq5F1lJs3JaG6Eu+Rki4DOiibHNbuNbjslDwOuLlsZQDgVbXjnPFhZsdAeVUJnGtm51Z/L+mjwMnAuyy8DMUSyorlRGA65Qf0dWZ2jpk9E4x5JZj9fz74/DFgZXC8n0t6grIieShQFGuAS8KaEpnZiKQbgR7g25Qn/W83s6EaubcAz5nZy8DLku4DFlL2pUTiKwjHKc8Gq38L0+MGBw+Pd1P+IS4E+qv2Ga56cIyQfBLWxR6n5rXAZ8xsAXBphDzTgMEqh+ciM6vnw3AmgaT3Up71n2pmO8LGmNnFZnaImc2lXCzzLjM7J9h/dvCvKPuTKq1vnwLeFXz3OmAe8MvAUf3PlDsg3lwlhyT9TtWxTmFPwcAfUOWwrkRKAf8C/L6kfSTNAI4h3Ik+BlcQjgNPAm+W9CpJBxL8WGM4ENhuZjsCu/+xkzm5pP9F2f9wVbDpAOAZlctuf6hq6IvBdxUn5ROSzgyOIUkLJyOHU5dvUr7+PwzMTv8AIOlgSXcm2P96SY9S9jMcBFwWbP8icFzw3b8Df2FmzwF/CLwDOFdjQ2MFfKfqWLOBLwTH+hywOPBz/BT4FICZPQ58H3gEeAi42szq9mb3Yn1OYZG0D/AbM3utpC9Tbs/6M8pmotvN7FpJ9wCfN7M+Sb+i7HN4Eeil3K1vM+XIkhVmdo+kl8xs/+D4ZwAnh5gpzgVWAQPADOAJ4Atmdn/w/acpOzmfpPwAOMDMzg0iUq4CXgHOoOwE/XvKD4gS8D0z+wKO0yBcQTiFJZhxX2Vmb2u1LI6TRdzE5BQSSZ8CbgAuabUsjpNVfAXhOI7jhOIrCMdxHCcUVxCO4zhOKK4gHMdxnFBcQTiO4zihuIJwHMdxQvn/nUhmzhms04EAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"data=np.loadtxt(os.path.join(data_dir,'asteroid1.csv'),delimiter=',',skiprows=1)\n",
"t=data[:,0]\n",
"y=data[:,1]\n",
"plt.plot(t,y,ls='none',marker=\"o\")\n",
"plt.xlabel('Julian Date')\n",
"plt.ylabel('Differential Magnitude')\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Why does the plot above show four vertical strings of points? Each vertical \"string\" of points is one night of observations. Compared to the entire period of time along the x-axis, each night is a very small length of time. If you alter the code in the cell and change the limits on the x axis, you, can view each night of data individually. (This will help you get an idea of the length of the asteroid's rotation period, so this is an especially useful exercise.)\n",
"\n",
"It's not unusual for astronomers to have large gaps in lightcurves due to daytime and longer periods of time when observations were not able to be made (i.e., weather, length of observing cycles). Phase-folding, or phasing, the lightcurve to trial periods allows one to remove these gaps and visualize the data more easily.\n",
"\n",
"Below, you will try phasing the data for 1856 Ruzena to trial periods to see if you can determine the asteroid's rotation period by the trial-and error method. Don't spend more than 5 minutes trying different test periods - you could do this all day and never get lucky with the correct period, since you'll need several decimal places of precision. (The students from the class published the rotation periods to 4 decimal places of precision using only these data.)\n",
"\n",
"Also, consider carefully: when you have correctly determined the asteroid's rotation period and have a phased lightcurve, how many maxima and minima should it have in one complete phase?"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEGCAYAAACQO2mwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2df5hcVZnnv293KpMKMOkg2R1SJCSw0JEQkoZGsmYcDTOPzYBACwhmRcdZRmbc1VVke40ua8KsDhlbxN2dWWdRWR4XjAHDUwPCTNxnEgcHTaQznRACZB9FCKlkJE7SqEmTVHe/+0fd27lddc79UXV/VHV9P8/DQ/etU9Unp773vue8533fI6oKQgghxERH1h0ghBDSvNBIEEIIsUIjQQghxAqNBCGEECs0EoQQQqzMyLoDcXLWWWfpokWLsu4GIYS0DDt37vyFqs6zvT6tjMSiRYswNDSUdTcIIaRlEJFX/V6nu4kQQogVGglCCCFWaCQIIYRYoZEghBBihUaCEEKIlWkV3UTah+JwCYNb9uHgyCjmd+Ux0NeN/p5C1t0iJFGy0D2NBGk5isMlfOaxPRgtjwMASiOj+MxjewCAhoJMW7LSPd1NpOUY3LJv8kZxGS2PY3DLvox6REjyZKV7GgnSchwcGY10nZDpQFa6p7upQegbT5/5XXmUDDfG/K58Br1pT6j79MlK94mtJETkARF5XUSer7r+cRHZJyJ7ReSLhvctEJFtIvKi0+YTSfWxUVwfYWlkFIpTPsLicCnrrk1rBvq6kc91TrmWz3VioK87ox61F9R9NmSl+yTdTQ8CuMp7QURWA7gewCWquhTAlwzvGwNwp6q+FcBKAP9eRC5KsJ91Q994NvT3FHDPDctQ6MpDABS68rjnhmWcyaYEdZ8NWek+MXeTqj4tIouqLn8UwAZVPeG0ed3wvkMADjk//0pEXgRQAPBCUn2tF5sv0LQkJPHS31OgUcgIm765J5Q8Weg+7Y3rCwG8Q0R2iMjfi8jlfo0dI9MDYIdPm9tFZEhEhg4fPhxrZ4Pomp2zvsalN5mO+Ol6Tt5+P5DWJW0jMQPAXFTcSAMAHhERMTUUkdMBbAbwSVX9pe0DVfV+Ve1V1d5586wl0RNB1f4al95kOuKna/OdTFqdtI3EAQCPaYUfA5gAcFZ1IxHJoWIgHlbVx1LuY2hGRsvW17j0JtMRP12PHLffD6R1SdtIFAFcCQAiciGAmQB+4W3grCy+AeBFVf1yyv0LTZA7ieGYZDri52L1e420LoltXIvIRgDvAnCWiBwAsA7AAwAecMJiTwL4A1VVEZkP4OuqejWAVQA+CGCPiOxyPu6zqvpUUn2thyB3EsMxo8PY++bHz8Xq9xqx0+y6TzK6aY3lpVsNbQ8CuNr5+R8ANL13M8id1ExfcivAekytwRs+Lla/14iZVtA9y3LUiZ87qUBXU2QYe98a+OmeLtbotILuaSSqKA6XsGrDVixe+yRWbdhq3XtYvcQcSdXZIXQ11QHrMWULdZ8NraB71m7yEGXpt+0lc07GGb8xo2mWia1EvXVpmt2f2wpQ99kxJ58zRkk2k+65kvAQZelns/T0y9aHbYZquw6whlBcUPfZUBwu4djJsZrruYBVWdq6p5HwEKXcgM3Sd4jwIVUHthmq7TrQGv7cVoC6z4bBLftQHq8NCTt9lv+qLG3d00g4FIdL1pAq041hqsgIAOOqnM3WQT2+2Vbw5zY71H122HQalJSYtu5pJBwGt+yDKcxbYM556O8p4MbLzNaes9no2BKx6ommYZRNeKj77KhXv2nrnkbCwWaFFeZ45eJwCZt32mdNnM2Gpzhcwq/frPXNAv57EjxXonGo++ywafv1X476rsjS1j2NhIPNCttyHkx+wTCfR2oZ3LIP5Qlzuu6mH79mvGHc6I7R8jg6ncpyPFciOtR9dtj228oTwMCju5tG9zQSDlGja/xmTJzNRsNvLMsTWuPC8EZ3ABV/uKDyXdFARIO6z45W0T2NhEPU6BrbjKlThLPZiATNPqtvJtNsVgE8tH0/N04jQt1nR6vonkbCIeopcza/4L03L+eNEpFFb4m2Uec3A1v/+N5Y+tQuRI2Uoe7joThcwuu/9N+/aRbd00g4+Fn1u4p7aq7xnOV4KA6X8MxPj1hfNyUW+X1Xfmd8kFpsp8nNypkfDdR94xSHSxh4dDfKE/Y2zaR7luVwGOjrxh2bdhnDAR/evh+9555ZcyPwnOXGCQqZHHxf7Qx1oK8bn9y0y/IOEgXbaXKj5QkUh0tGfVP3jeEXqDHZxqD7RW8xl65JGq4kHPp7CkYDAVT8foz/TgaWXM8Wv8Qtaj4ZgjRf6MrX6L44XMIPfVbcHQkerkAj4WGuz8lajP9OBr8ldJfFFcKHV3z4nSZHzSeDn+ZtSYy2pEeXgIVJQ9BIOPgldAHmLzZseWVixy9Z7j3LzzZeD3p48XsIBzWfDQN93dYH78wZ5lfCGGwW+EsYPz+hKf6bFUjjwa+Any2RzrbZ6sKVRjj8NJ/rrN04pebjob+ngDmWFdyJsQljIl2YJEUW+EsYP0t942W1G3WsQBoPQQlFptA+v4zfoM8kp/AdJ4PtoObjw28vyKT7MEmKLPCXMH6W2jTbZQXSeLCFWrpUh/YVh0s4MeYTOwh/Pzs5hZ/mTRm/1Hx8RNV9GFjgL2H8LHWUuvqsXRONoAc+MNXXGmbWqglu4k0nbGW/Xap1T83HRxK6Z4G/hOnvKVijm8LW1WftmuiEicrw3iBhZq1MqAuHmxjXaUmWqNY9NR8fSeg+KWgkHIrDJeMM1HYTMPM0HmwPKC/eGyTMrDXMZ5JTnDGrNqfWpHtqPj6S0H1Se0PMuEbtQfAuc2fnsO7apdabgJmnjbPmigV4aPt+3zbeG2Sgr9v4XXkZp78pFPXonpqPhyR0n9Rqg0YC9hr5s2faz5p167ofHBnF/K48Bvq6efPUwef7lwEAvrVjv3EJXh2K6Y7xnY/sthoD21kIZCrUfXYkoXtuXCeIXwVYU8IQ48Xj5fP9y/Dlm1cgZ6gtcMvlC4w1s+69eTlynbXtTYXRiBnqPltaRfc0EvC3wKYbwRYv/slNu5iFWie2xC5bsl1/TwGDNy2fEmzQlc8ZC6MRM9R99rSC7ulugn8FWOBUwpD7Jfj5/tybC2BxuijUE4NP/3hjUPfZ0wq650oC/hVgXaJEGjALNTq2MQ0qwUHqh7rPnlbIPaGRcAja7KyONPBLQgKYhRqVgb5uo2/22MkxujEShLrPlqhnjGcBjYSD3w1QHTPujRe30UwzgVagv6eA0w3x+uXx2vIQJD6o+2yJesZ4FnBPwsH18Q1u2YfSyCg6RTCuioIlzM/1C5pizZmFWh+2omfVs1OGYcYHdZ8tYfckstQ8jYSHejaEvDcZH1qNMb/LfDyjd3Z66nzgije9NDKKgUd3A+CGab1Q99kRVvNeg1waGcUdm3Zh6NUjk/kWSUIjEQOMsomH1UvmGbNQvf7Z9Y/vrQkZdEsr8ztIF+q+ccJo3hR6rAAe3r4fveeemfh3wD0J0jSE8c/aivexqB9pRcJo3uaSUqRzwFZiRkJEHhCR10Xk+arrHxeRfSKyV0S+aHnvVU6bn4jI2qT6aINHNGYDzyvIFuo+ffyy3l38ggHSuDeSXEk8COAq7wURWQ3gegCXqOpSAF+qfpOIdAL4SwC/D+AiAGtE5KIE+zkFlh7IDltOhPe6rZy77ToJB3WfDX4GwB37gb5u2GrGphFNlpiRUNWnARypuvxRABtU9YTT5nXDW98G4Ceq+rKqngTwbVQMSyqEPaKRs674KY+bD2LxXl937VJj7ZprLjk7sX61A9R9NvhFg939ROUI0/6eAj6wcmGNoRCkk0+R9sb1hQDeISJfAPAmgP+oqs9WtSkAeM3z+wEAV9g+UERuB3A7ACxcuLDhDoZxeZiiDViSoD68oX227N9jJ089vPp7Chh69Qge3r5/SvuHnc2/NKI9piPUfbp4dW/jqCck3NW1V/fu5rX39SRIe+N6BoC5AFYCGADwiEjN6RumlZW1eoCq3q+qvaraO29e41Y1TJo8D4SPh2oXR1i2vXS4pr17w3BmWx/UfXq0mu7TNhIHADymFX4MYALAWYY2Czy/nwPgYEr9C5Umzw3WeLCdZ1BN9ZnxWUd7TEeo+/RoNd2nbSSKAK4EABG5EMBMAL+oavMsgAtEZLGIzATwfgCPp9XBMCFprVCUqxUI+3ApTwB3FfdM/p51tMd0hLpPj1bTfZIhsBsB/AhAt4gcEJHbADwA4DwnLPbbAP5AVVVE5ovIUwCgqmMAPgZgC4AXATyiqnuT6mc1YWZLPBA+HqI8XDbuOLVN5Rft0SFCl1MdUPfpkYTuuxKM7ksyummNqp6tqjlVPUdVv6GqJ1X1VlW9WFUvVdWtTtuDqnq1571PqeqFqnq+qn4hqT6aCDNbqj4Qviufw6xcB+7g4SuRCFNV1MV7ZKMt2sNtx9DN6FD36dGo7k38+s3kqiUz47qKsLOl/p4Cnll7Je67ZQVOjE3g6PEy48sjUv3Qcf9vorMqvuHz/ctw3y0raq4D3EytB+o+PUy6t2HSfZchn6g8kVy1ZBqJKkxf4D03LLOG+DHiozHch87PNlyDZ9Zeibeff6ax3crz5hrfO2E5FJ57E9Gg7tPFq/uBvm7rg3jNFQtqrr1hKUGTlOZZ4M9AlMJljPiIj+JwCT98uTr/ssLeg7+abOMta20LIeRmanSo+2y4+4m9MKWRduBU/oNX9zaS0jyNRIOEKfVLgnFLgFsWBhgZLeOu4h5s3lmanMGOWxrnOoSbqQlD3cdDcbg0JWnOywQq0U29555Zc3aHiaSyr+luahBGfMTD4JZ9NSXAq9m447VQ8eWnz5rBDOCEoe7jIcg9t3HHa6HzKpI6zY4riQbh4SvxEMZNYVs5VGM74Y7EB3UfD0G6H1cN7cLjnkTCNHI8IA9faRyb+8KLe7RmmM8i4aDusyVI950i+K05swLvDfezkoDuJtRXJpnVMONloK8buQ5bACzQIeYIp2ro8ggPdZ89Qbpfc8WCUHsNSeqeRgLRw/lYez9++nsKGHzfcuSrC9Y4dHYIXjj0K9/PCArbJFOh7rPH1b1J9p0dgt5zzwzca0ha93Q3IXo4n9/NxQdUo5hnVeVxtUaBAMBXblnBsY8Idd88jBliYMedBDm/vYZXNlyTYK8qcCUBuy+vQ8S4rGaMeDKEjeIwwRltdGy6V8DoSqLuk2Fwyz5rvo+7V2RCgFQ0TyMBey2VcVXjsprVMOPD6+MO2pzryuesNW+Y7RsdvxpCJlcSdR8fYXXvBhPYDtlJQ/M0EqgtSRBUD4gx4vEQ5fCVfK4T669bintusJ/AVRoZ5WZqBLy6N1FteKn7eAirewEmo81s7Uojo4kHEdBIOHhrqQTVA4pa54aYCetemjs7Nzm+/T0F34Jo3EyNhqt7W3yN15VE3cdDWN1/YOXCybH103zSQQTcuDYQpuQAY8QbJ4wvO9cpWHft0iljPdDX7VumgJup0QlbZoO6b5ywuu8991SxyyDNA8npnisJA1xWp0MYX3Z5vLYEcpCbBOBmalSo+fSoR/fVqzgbSeieRsIAl9XpEPbwFZPwXTeJzVBwMzUa1Hx61Kt7r0s8Td3T3WSBy+rk8db/qbcEsmkZzhlwfVDz6dBquudKgmRK0IrAjfDwez9nwKTVaCXdi4asrNkK9Pb26tDQUNbdIHXghgV6Z0aCSoSHe/AKIdONZtC9iOxU1V7b63Q3kaaApadJO9IKuqeRqINGyiuTWqrH8z6fOkwc++zg2MdLq+ieRsKC7UupXh66SSwAeMPUQZTx5NgnD3WfDq2ke25cG/AriRy1vDLxJ8p4cuyThbpPj1bSPY2EAb8vhZUw4yXKeHLsk4W6T49W0j2NhAG/L4WVMOMlynhy7JOFuk+PVtI9jYQBvy+F5QviJcp4cuyThbpPj1bSPTeuHbwbdnPyOeQ6BeXxUzkk7pfSCiFrrUSU8eTYx0v1JvXqJfOweWfJmMXLsY+XVtI9k+lgTmjJdQhOnzUDI8fLvCHItMOk+XyuEzdeVsC2lw7TELQRTKYLgWnDrjyhmD1zBoY/9+6MetUexBH/zfj96Ng2qbe9dBjPrL0yo161B3HpNS3d00gg++iBdiWO+O+7invw8Pb9kyd3MX4/HNR8NsSV85Cm7rlxjeyjB9qVRuO/i8MlPOS5Uer5jHaFms+GOHIe0tY9jQSyjx5oV2yzVr/yyV7WP7438meTCtR8NvhpPuzRo2nrnkYCLDedFbZZqwChbpiR0XLkzyYVqPls8NNl2DOq09Y99yQceOBK+gz0deOOTbtqls0KNHxWL2fEwVDz6eN3VnUcZ1QnofvEjISIPADgPQBeV9WLnWvrAXwEwGGn2WdV9SnDe+8A8EeoPC/2APhDVX0zqb5G5a7iHmzc8RrGVdEpgjVXLOCZB3ViC8C2LZu9ER0igCmC+7SZnXz4JQB13zj9PQUMvXoED23fb3zdpPvqKKbTZnbi2MlaI5OU7kMZCRG5SFVfqLr2LlX9vs/bHgTwFwC+WXX9PlX9ks/fKgD4DwAuUtVREXkEwPudz8ucu4p7pnzB46qTv/OGCY8b5WGjQwSL1z7pW4nUZGFynYIvvJffQ9xQ9/FQHC5h8067S6la9wBqoqFyHYLODsH4xKkbIEndh92TeEREPi0V8iLyPwDc4/cGVX0awJE6+zUDQF5EZgCYDeBgnZ8TOxt3vBbpOjFjivLwMq4aqhIpAHQIJv3qgzct5yoiAaj7eIiq+7uf2GvM4ZqYUMydnUtF92HdTVcA+HMAPwRwBoCHAayq829+TEQ+BGAIwJ2qetT7oqqWRORLAPYDGAXwPVX9nu3DROR2ALcDwMKFC+vsUnjGLRnqtuvETJQojKBKpBMKfMXnwBbSONR9PETVvc2gKIBfnxjzPagoLsKuJMqoPLDzAGYB+JmqTtTx974K4HwAKwAcAnBvdQMRmQvgegCLAcwHcJqI3Gr7QFW9X1V7VbV33rx5dXTJTnG4hFUbtmLx2iexasNWFIdL6BQxtrVdJ2aiRmGUfCqRAsDdT9jDAkk0qPvkiDP6qDyuqeg+rJF4FhUjcTmA3wawRkS+E/WPqerPVXXcMTBfA/A2Q7PfQ8UIHVbVMoDHALw96t9qFNsBLCvPm2tsv+aKBel2sMUxxen7IQBWL7FPAo4eL9c82Eh0qPtkiap7oLLfYOPo8XLimg9rJG5T1c+pallV/0lVrwfw11H/mIic7fn1vQCeNzTbD2CliMwWEQHwuwBejPq3GsWWGfnKP4/i1pULJ2dQnSK4deVCbt5FpDpOPwgFsO2lw75tTCeqkWhQ98ni1X1YTpvpvyuQtOZDVYEVEaOzX1XNcVyV92wE8C4AZwH4OYB1zu8rUPl3vQLgj1X1kIjMB/B1Vb3aee/dAG4BMAZgGMAfqeqJoH7WWwXWxOK1TxrDMwXAzzZcE8vfIKdYtWFrYKa1AJiTz/kmE3kpdOVZrC4i1H16hNE8EE339Wg+qAps2JXEkwC+6/z/7wC8DOBv/N6gqmtU9WxVzanqOar6DVX9oKouU9VLVPU6VT3ktD3oGgjn93WqukRVL3beE2gg4oa1bdKjOFzCyPGTge3md+Wx/rqlyHWE84OzNEd0qPv08HOfeomi+8zKcnge7MtU9QJU9hL+IfbeNBGsbZMOxeESBr6z25gcVM3qJfPQ31PA4PuWTyknMXd2ztieD7boUPfpUBwuYdOz4cKHTbq3BQw0TVkOVf1HEbk87s40E1mfBtUuDG7ZN+UEQD/cPYnqchK2A3T4YIsOdZ8Ojeo+Tc2Hzbj+lOfXDgCX4lRpjWkLa9skT5TlcWlkFKs2bK15eLmlDrwlI268jN9dvVD3yROH7oFKRVh3r2JWLpl6rWFXEmd4fh5DZW9ic/zdIe3G/K586NLgglNlxL2HrADA5p2lycSucVVs3llC77ln8mFHmpK4dH9i7FS62tHj5UQOHuIZ1z7wWMzkcfckwi69q3FDCU03HKObokPNp0Mz6b6hM65F5AnYC3VCVa8L3ZMmxnRjALWFtXgsZvy4Y3n3E3tx9Hi40FYvfst2Rjf5U6371UvmYfPOEjWfAq2ke9+VhIi80+/Nqvr3sfamQepZSdg2gGblOoxfHmenyRJ1hsWVRH2YdC8wzwg5jsmTpe4bWkmgUh7DmjA3HbBlmNoKa3F2Gh8214Z3M86L6SF27MQY3rP87CkzYIDRTUGYdB/1bA9SH62m+6Dt8KL7g4hMy43qqDcAY+/jwVYjqDhcwhuWzFIFanIiRkbL2LyzhBsvK/AozghE0T01Hx+tqPsgI+HN2Dgv1r/cJNhugK58jklFCWJbwQ1u2Wf9Tgpdecw21LEZLY/ju7sPJdLP6Yrf+eJeqPl4aUXdBxkJtfw8bbBlmK6/bikPik8Q20y2NDKK1UvmWQ207X0jo2UW+IuATfcfWLmQmk8QP937ZbtnqfugPYnlIvJLVCYYeednOL+rqv5mrL3JgKAMU94gyeAXJ+4uo7e9dLjmOxncsi9UfHkch8pPZ5hZnQ023bsruHtuWGb8TrLUPfMkSCaYomu82CI0gt7nhZVLSbNRHC7hjk27IkeRJan7RqOb2oKoCURMOGocd7w+uWmX8XXb8to0Az5+cswYrswNV3+i6Jiaj4f+nkJkzbvvA7LRfdsbiWoLHZRAFLU9seO3jDYJ3X1QlUZG0SkyORu75hKGwEYlio6p+XgpWFxOtod7cbg0JemuK58zJvwCyei+7d1NtoM/bEu/qO2JP6ZldK5DcPqsGTh6vIxOEYyroiufw7GTY8Zko3yu07qHQcxE0TE1Hy+2BN4bLyvgyecOTTEG71l+NjY9+1qN7jsE+PLNKwA0vq9Ed1MAtiXewZFR4xLbrz2JTvUyeo5jDNwbxS3a53cq12h5HNteOswHVgSi6N62YUrN14fJdbR6ybwaYzAyWsZD2825zBNaKekx/Ll3Jz4ZansjYbsJ5uRzxiV21+wc/d8x4y1NvWrD1tDHk3rhAysaUXRvK9dBzddPdTn2VRu2Ri72V0/Np3pIpgB5C2GLTRaBMelFFUyyS5B6H/Z8YEUjiu4VTLJLmmae5LS9kejvKRiT5kYsVvqN0TKT7BKknoe9AHxgRSSq7tVpQ80nQz2678qbj+2Nm7bfuLbBzbrkCVOiPQy3rlyIz/cvS6qbbQV1nzw23UepAtsB4Mu3rIjFUAdtXLf9SsIGD4RPFluhM6CSdWo76L2aubNzNBAxQt0ni5/uB29ajo5wssec2bnUVnI0EhZsy3EusePBr9BZf08BEyFXuDb3CKkP6j5ZgnQf1rGTpu7bPrrJDx4InxxBocRhzwDump2OX7adoO6ToxV1z5UEyQTbRp173eT2MPHrN8dY7ZW0DK2oexoJkgmmm0EArF4yb/L3WblT8sznzFItTygGt+xLpI+ExE0r6p7upoiw0Fk89PcUMPTqETy8ff9kopaiUiYczv+9vtvR8oT1s2zZ8fxe4oPjGw+tqHuGwEbAVnOFG3v1YQu3dOs1haUrn8OJsQl+LwlB3cdLs+meIbAhKA6XsGrDVixe+yRWbdhq9fX5RSaQ6Ng28aLcKH7Z8fxe7ITVPEDdx02r6b7tjYTfweTVsLhfvNg28cLmSHSK+GYJ83sxE0XzAHUfN62m+7Y3ElFmSUGRCSQatsStlefNDXxvPteJe29ejv6eAr+XiERdGXB846XVdN/2RiLKLInZqPFiS9x65Z/N30mniDHBi99LNKKuDDi+8dJqum/7jWu/WjUDfd3GGiuM8kiWxWufNJam9ju7l9E34Qmqz2SrLcTxTZasdB+0cd32RsLvlCjTkZiM6IgXk8htR5raDDe/j2j4RSsB5iMxqft4iaL7ubNzmD1zRmKaZ3RTALal37aXDjOiI2FsG6irl8wzLqNXL5kXacOVmPGrz8RIpuSJovtcp+DXb45lqvnEkulE5AEA7wHwuqpe7FxbD+AjAA47zT6rqk8Z3tsF4OsALkYl1+TfquqPkuqrqVbNHZt2GdsyoiM+bA+kbS8dxj03LDPOtPyKo5Hw2OozMZIpeaLo/tiJsZqTGtPWfJIZ1w8C+AsA36y6fp+qfingvf8NwN+q6k0iMhPA7AT654ut0BYjOuLD74FEw50N1H3yRNH94rVPRvqMJEjMSKjq0yKyKOr7ROQ3AfwOgA87n3MSwMk4+xaGgb7uGt9sdY0VUj/F4RI6LBmm7gOp2m/L88WTZ6CvGwOP7kZ54tT3kusQRjLFRFTdB7VNgyxqN31MRD4EYAjAnap6tOr181BxR/1vEVkOYCeAT6jqMdOHicjtAG4HgIULF8bWSVuNlW9t349v7diPCa2Epq25YgEPvYmI65M1id8N4aveXC2NjCLXIch1ypTTuxiKmQDVOV0CPDq0H3c+shvjqtR9ndSje7+2aZH2xvVXAZwPYAWAQwDuNbSZAeBSAF9V1R4AxwCstX2gqt6vqr2q2jtvXvRZvl95gm0vHa4JSZsA4E6yxlXx0Pb9uKu4J/LfbWdMPlmg8mzy20AtTyhOmzmDB+LEgE33g1v21RyhWR5XPPPTI5MPLOq+Pmy6dzOobbp322Sl+VRXEqr6c/dnEfkagO8amh0AcEBVdzi/fwc+RqIRTLNV9yjB/p5CaL/fxh2vcVYVAdu4aog2b4yWsWvduxPoVfvgp/sovm7qPhq2w4TGVScf+rbxn1C15kokTaorCRE52/PrewE8X91GVf8JwGsi4q6nfhfAC0n0Jyjcb04+3OlPUQpzEX9/qjv2tjYKYNHaJ7Hi7u8x9LVO/HQfxddN3UfDVpvJe932zMlS94kZCRHZCOBHALpF5ICI3AbgiyKyR0SeA7AawB1O2/ki4g2F/TiAh512KwD8WRJ9DAr3C1lvK3RhLlLBz5/qjn1QgMDIaBkDj+6moagDP90P9HXXbEnYoO6jYTOq7vXicAnHTo75fkYWuk/MSKjqGlU9W1VzqnqOqn5DVT+oqstU9RJVvU5VDzltD6rq1Z737nL2GS5R1X7D5nYsBBXICnvY+JorFsTWp0nVKysAABEISURBVHagv6eAuZYzet2x3/bSYePrXngqXX346b6/p2AsDWGCuo9GwTLu7nXTfpCJtHXf1hnXQQWygpbenSK4deVC+mXrYN21S33HPqxvnDkS0QnSve1h5kLd10fQuEfR8rTIk2gF3M0iWy0gU64E69jEQ9DY25K6qmGORHSo+2yIS/Nu27Ro65VEEH41bkjj9PcUMNDXjfldeRwcGcXgln2TvlbTrKsaJnklA3WfHI1qHkhf9229kggKgXX/77053PhyViFtHL/xB4DfmNFhjBl3ueVtCzj2dRBV924G8B2bdlHzDRI09kOvHplM1rWRtu7b2khELRgX5uYi4bGN//rH99Yc8G4izOY2qSWK7qn5eAkKu9/07Gu+BgJIX/dt7W6KWvGSZZTjxTbOI6PlQAPh937iTxTdU/Px4jf2YaOb0tZ9WxuJqGfEsoxyvDS6+cZN6/qIontqPl78xj7smKat+7Y2EmHOiPXWuOmwJA/xYVUftvG35VBUt+OmdX1E0b1tXkvN14ff2IcZ0yx039Z7EkEhac1YkXE6YRt/oPYITS9zZ+ew7tql9InXSVTdV9PJqLK6CQw//s5uq8upkFHQQFsbiSBsFRm9XLpwDh9WdVB9VsR9t6yoGcf1j++dciqXCKAKzJ5J2SZJkO7HJxRDrx6h7uvAdLa1N6IMAD772HM4Xp6Y8r6sDATQ5u4m21mzbtxyGB/h9pcTqRgyrTGN+x2bdtWUnj4xNvVGcRdyPNu6MeLQ/cYdryXcy+lHWN2roXpWlppvayMRFLkRxkfISpjRMY27Anh4+/4p5xr4zWYZYVM/1H02NKr7rDTf1kYiKHIjTAYkK2GGx90MtZUeUJwqFR5mNhu2hAGZCnWfLnHqvjQyOn1KhbcCQaGA3vIENlgJMxzepbYf7k0Q5iwPcT6XRIO6T48kdJ+226mtjUSYUMD+ngKeWXslXtlwDW5duXByBsVKmNEIEwTg8pnH9qA8PhHYzjsDI+Gh7tMjCd2n7XZq6zCRoHC0aj7fv4w3R51ESb4Ke1MBdDnVA3WfHtNB921tJIDaAn4kGaKUQY4CfeP1Qd2nw3TQfVu7m0h6hC2D7NKVzyHXGXwjMMqGNDPTQfdtv5IIQ3G4VJPYxazfaHhdHEEzq3yuE+uvWwoAuPuJvTjqc4xs0ClqpH6o+8Zxx+nOR3YHPtibVfc0EgEUh0v41CO7asr3Hj1exsB3dgNgyeSwuOP0yU27rG1MmaW2EhEsiZIc1H189PcUfDUvgHFfqFl0TyMRwN1P7LXWdy+Pq/XsCVJLcbiEgUd3W18vdOXxzNorp1y7+4m9xhulU4SnpSUIdR8ffuGqXfkcdq17d831ZtI99yQC8FvyASyZHIXBLftQtjx5BKiZHd1V3GMdf+5FJAt1Hx9+4aqm/ecg3adtnGkkGoQlk8Pj92BRTHVfFIdLeHj7ft/PY/2m7KDuw+On+5EqY1AcLuEhH91nkUBKI9Eg9ImHx+/BUr0RN7hln/UsAxfWb0qGMA8h6j48frqvfm3943t9PyuLBFIaiQCCgtHolw3PQF83ch21I9ohwLETY1i89kms2rAVxeFSaHcG3R7xE+YhRN2HJ4ruvZFkNnh8aZPhN5tlIlc0+nsKeNviuTXXJ7RyrrW3bHVXiNPpALo9kiDoIcSw42hE0X0YeHxpk+F3Q3DzNBrF4RJ++NMjge1Gy+MIM7S5Tp6QlgRBDyGOeTSi6D7MvHP1knkx9Co8NBIBDPR1W11OnFFFI8w+g8sbIZbdp82cQbdHAvhpviuf45hHJIruw0yOtr10uKH+RIVGIoD+ngI+sHJhzU3DRK7oRPGlzu/KBxrhMIaERMdP825GMAlPlNpNYSae3JNoQj7fvwz33bICha48BJUvkolc0QnrS3UNcJAR5n5EclDz6ePqfm7Afhz3JMi0ZdFbwonbfRj19xR8BcqVXLL09xQw0NeN+V15HBwZxeCWfcxLiUiU8XJ1v+5a/9Va2rpnWY4QuKdLuWny3kgEzqzCs/3lo6HauWNaHC7BdgTLb8zo4NgnDHXfOFFyGsKMaRa650oiBEEHx5NwhIkG8/pk/cb3xNgEZ7UJQ903Ttj9g2bWPY1ECIIOjifxUB0MEDS+fFglC3XfOGH2D5pd9zQSIQg6OJ7Ew6UL50xZSgeNLx9WyULdN06Y/YMbLys0te4TMxIi8oCIvC4iz3uurReRkojscv672uf9nSIyLCLfTaqPYQlzcDwJJii875mfHpmylGZ0U7ZQ943T31MIjFbavPPAlN+bTfdJriQeBHCV4fp9qrrC+e8pn/d/AsCLifQsIv09BdxzwzKGAzbIQF83DCVspuBdSvf3FHxrZ/FhlSzUfTysu3apr+5Hy1PDM5pN94lFN6nq0yKyqJ73isg5AK4B8AUAn4qxW3XDg+MbJ8zJdNWJR35b3fw+koe6b5wwui8Ol6aMczPpPos9iY+JyHOOO6q26lWFrwD4T4A1AnISEbldRIZEZOjw4XTT1Ul0ggReXS/f5qLKoq4+IfUSpPvqs1GaSfdpG4mvAjgfwAoAhwDcW91ARN4D4HVV3RnmA1X1flXtVdXeefPSLXxF6sPPR1tdL99WRyiLuvqENMLsnP1xWx1a3Ey6T9VIqOrPVXVcVScAfA3A2wzNVgG4TkReAfBtAFeKyEMpdpMkzLprlyLXafe6el1O/T0F69Kb0U2klfizGy7x3Zs42KS6T9VIiMjZnl/fC+D56jaq+hlVPUdVFwF4P4CtqnprSl0kKdDfU8DgTcut53GEdTkxuom0Ev09BXz55hXWcuBz8lNX2M2i+yRDYDcC+BGAbhE5ICK3AfiiiOwRkecArAZwh9N2voj4RTqRacjMGea7xeRyYigmmTZYlgjVxqNZdJ9kdNMaw+VvWNoeBFCTM6Gq3wfw/Vg7RjKnOFzCwKO7UZ6wx3BUL72BiuE4ODKK+V15DPR1M+qGtBSu7m2qHzk+tfR9s+ieBf5I6gxu2edrIIDaJTVDMUmrE6R7kxupGXTPshwkdYIOYeGxpGQ6EqT7tI8lDQuNBEkd24a1C48lJe1I2seShoVGgqROUMlwHktK2pFmDemmkSCpE7SSYGgraUeaVfc0EiR1glYSzeqbJSQpmnkfjkaCpE5QyfBm9c0S0gh+um/mfTgaCZI6piQhL83qmyWkEfxWCs28D0cjQVLHPafAtjfRrL5ZQhrB7wCi6pIczQSNBMmE/p4C1lyxwPga9yTIdGXdtUuRM1T5O3ZyrGlL39NIkMyw7T1wT4JMV/p7Cjh9Vm2hi/K4Nm3pexoJkhm2vQfuSZDpTHWNJpdm1T2NBMkM294D9yTIdKbVdE8jQTKjWUohE5ImraZ7VoElmdEspZAJSZNW071oQPZrK9Hb26tDQ0NZd4MQQloGEdmpqr221+luIoQQYoVGghBCiBUaCUIIIVZoJAghhFihkSCEEGJlWkU3ichhAK8m+CfOAvCLBD+/XtivaDRjv5qxTwD7FZVW7Ne5qmotmDatjETSiMiQX6hYVrBf0WjGfjVjnwD2KyrTsV90NxFCCLFCI0EIIcQKjUQ07s+6AxbYr2g0Y7+asU8A+xWVadcv7kkQQgixwpUEIYQQKzQShBBCrNBIVCEiV4nIPhH5iYisNbz+KRF5QUSeE5G/E5Fzm6FfnnY3iYiKSCpheGH6JSI3O2O2V0S+1Qz9EpGFIrJNRIad7/LqlPr1gIi8LiLPW14XEfnvTr+fE5FLm6RfH3D685yI/FBEljdDvzztLheRcRG5qRn6JCLvEpFdjub/Puk+hemXiMwRkSdEZLfTrz8M9cGqyv+c/wB0AvgpgPMAzASwG8BFVW1WA5jt/PxRAJuaoV9OuzMAPA1gO4DeZugXgAsADAOY6/z+L5qkX/cD+Kjz80UAXklJY78D4FIAz1tevxrA3wAQACsB7GiSfr3d8x3+frP0y/N9bwXwFICbsu4TgC4ALwBY6PyeuOZD9uuzAP7c+XkegCMAZgZ9LlcSU3kbgJ+o6suqehLAtwFc722gqttU9bjz63YA5zRDvxz+K4AvAngzhT6F7ddHAPylqh4FAFV9vUn6pQB+0/l5DoCDKfQLqvo0KjenjesBfFMrbAfQJSJnZ90vVf2h+x0iPd2HGS8A+DiAzQDS0FaYPv0bAI+p6n6nfbP0SwGcISIC4HSn7VjQ59JITKUA4DXP7wecazZuQ2XWlzSB/RKRHgALVPW7KfQndL8AXAjgQhF5RkS2i8hVTdKv9QBuFZEDqMxAP55Cv8IQVYNZkJbuAxGRAoD3AvirrPvi4UIAc0Xk+yKyU0Q+lHWHHP4CwFtRmRDtAfAJVZ0IehOPL52KGK4ZY4RF5FYAvQDemWiPnD9nuDbZLxHpAHAfgA+n0BcvYcZrBioup3ehMvv8gYhcrKojGfdrDYAHVfVeEfnXAP6P06/AmyZhQmswC0RkNSpG4rez7ovDVwB8WlXHKxPkpmAGgMsA/C6APIAfich2Vf1/2XYLfQB2AbgSwPkA/q+I/EBVf+n3JhqJqRwAsMDz+zkwuCFE5PcA/GcA71TVE03QrzMAXAzg+86N8lsAHheR61Q1yfNcw4zXAQDbVbUM4Gcisg8Vo/Fsxv26DcBVAKCqPxKRWagUQUvFNeBDKA1mgYhcAuDrAH5fVf856/449AL4tqP7swBcLSJjqlrMsE8HAPxCVY8BOCYiTwNYDiBrI/GHADZoZVPiJyLyMwBLAPzY7010N03lWQAXiMhiEZkJ4P0AHvc2cNw6/wvAdWn5GoP6papvqOpZqrpIVReh4jNO2kAE9suhiMpmP0TkLFSW4i83Qb/2ozLTg4i8FcAsAIcT7lcYHgfwISfKaSWAN1T1UNadEpGFAB4D8MEmmBFPoqqLPbr/DoB/l7GBAIC/BvAOEZkhIrMBXAHgxYz7BEzV/L8E0I0Q9yJXEh5UdUxEPgZgCyoREw+o6l4R+VMAQ6r6OIBBVDZ9HnVmL/tV9bom6FfqhOzXFgDvFpEXAIwDGEh6FhqyX3cC+JqI3IGKO+fDzgwrUURkIyqut7Oc/ZB1AHJOv/8Klf2RqwH8BMBxVGZ/iROiX58D8BYA/9PR/ZimUO00RL9SJ6hPqvqiiPwtgOcATAD4uqr6hvCm0S9UAlseFJE9qLg1P62qgWXNWZaDEEKIFbqbCCGEWKGRIIQQYoVGghBCiBUaCUIIIVZoJAghhFihkSAkAk6l0V0i8ryIPCois0VkUVCVUkJaFRoJQqIxqqorVPViACcB/EnWHSIkSWgkCKmfHwD4V87PnSLyNadO//dEJA8AIvIREXnWqeG/2cnAhYi8z1mN7HbKNkBEOkVk0Gn/nIj8cTb/LEJOQSNBSB2IyAxUzlXY41y6AJWS6EsBjAC40bn+mKperqrLUSnNcJtz/XMA+pzrbsb+baiU4bgcwOUAPiIii5P/1xBih0aCkGjkRWQXgCFUauF8w7n+M1Xd5fy8E8Ai5+eLReQHTimEDwBY6lx/BpUSCR9BpXQIALwblbpNuwDsQKUMxgVJ/mMICYK1mwiJxqiqrvBecGoZeasBj6NSIhoAHgTQr6q7ReTDqNTWgar+iYhcAeAaALtEZAUq9XQ+rqpbkvwHEBIFriQISZYzABwSkRwqKwkAgIicr6o7VPVzAH6BSnnwLQA+6rSFiFwoIqdl0WlCXLiSICRZ/gsqrqNXUdm/OMO5PigiF6Cyevg7VM7hfg4VN9U/OkdMHgbQn3aHCfHCKrCEEEKs0N1ECCHECo0EIYQQKzQShBBCrNBIEEIIsUIjQQghxAqNBCGEECs0EoQQQqz8f5r5cI+wvo31AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"data=np.loadtxt(os.path.join(data_dir,'asteroid1.csv'),delimiter=',',skiprows=1)\n",
"t=data[:,0]\n",
"y=data[:,1]\n",
"period=2\n",
"#period in units of days\n",
"\n",
"def phase_fold(t,y,period):\n",
" phases=np.remainder(t,period)/period\n",
" phases=np.concatenate((phases,phases+1))\n",
" y=np.concatenate((y,y))\n",
" plt.plot(phases,y,ls='none',marker=\"o\")\n",
" plt.xlabel('Phase')\n",
" plt.ylabel('Flux')\n",
" \n",
"phase_fold(t,y,period)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From your experience with the trial-and-error method above, you can see that a better method of finding test periods is needed. This is where the Lomb-Scargle Periodogram method is used. This method essentially measures how much power is present at each frequency within the data: more power at a certain frequency should indicate a repeating signal within the data. However, that signal could be from a physical phenomenon, such as a pulsating star or rotating asteroid, or it could be an artifact, such as the day/night cycle or other artificial frequencies caused by the observing process. The astronomer must interpret the data to make a final decision.\n",
"\n",
"The code below applies the Lomb-Scargle Periodogram to the synthetic sine curve constructed at the beginning of this notebook. In this case, we know that the period is 1 day; usually one does not know this ahead of time, and so we will use the LS Periodogram to tell us this.\n",
"\n",
"The function below returns the best period and the power spectrum (which shows us the relative strength of signals at different frequencies (frequency = 1/period). You should see a large spike at period = 1. Try changing the value of the period and watch the periodogram (the power spectrum graph) change."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0003718852975687"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxcVfk/8M+TvWm60n0vUKSlsrYsoli1sqkUBQV+KqIioqCCiqJfQTZFFkFQQIsCiihFEGiltLTQBVraJt23pE3TLW3S7NtkmczM8/vjLrmzJJnJ5GY6mc/79eqrs9zcOXfm3vuc85xzzxVVBRERpa60RBeAiIgSi4GAiCjFMRAQEaU4BgIiohTHQEBElOIyEl2AWI0YMUKnTJmS6GIQESWVjRs3VqnqyEjvJV0gmDJlCgoKChJdDCKipCIiBzt7j6khIqIU51ogEJHnRKRCRHZ08r6IyJMiUiwi20TkbLfKQkREnXOzRfACgEu7eP8yANPMfzcBeMbFshARUSdcCwSquhpATReLzAPwDzWsAzBURMa6VR4iIooskX0E4wEcdjwvNV8LIyI3iUiBiBRUVlb2SeGIiFJFIgOBRHgt4gx4qjpfVWep6qyRIyOOfiIioh5KZCAoBTDR8XwCgKMJKgsRUcpKZCBYCOB6c/TQ+QDqVbUsgeWhJLZ4exmqm9oSXQyipOTaBWUi8m8AcwCMEJFSAL8GkAkAqvpnAIsBXA6gGEAzgG+6VRbq32o8Xnz/pU04a9JQvP79CxNdHKKk41ogUNXrunlfAdzi1udT6vAFAgCAwzUtCS4JUXLilcWU9NLEGHcQ4N32iHqEgYCSXroZCPwBBgKinmAgoH4jwEBA1CMMBJT0rNO/n6khoh5hIKCkZ/UNsI+AqGcYCCjpWed/ZoaIeoaBgJKemskh9hEQ9QwDASU/8/zPPgKinmEgoKRnnf4ZB4h6hoGAkh47iYniw0BASY9xgCg+DASU9BgHiOLDQEBJT9kkIIoLAwElPcYBovgwEFDSYyAgig8DASU9ZS8BUVwYCCjpsUVAFB8GAkp6jANE8WEgoKTHC8qI4sNAQEmPcYAoPgwE1A8wEhDFg4GAkh5bBETxYSCgpMc4QBQfBgJKeuwsJooPAwElPcYBovgwEFDSYyAgig8DASU9TjFBFB8GAkp6bBEQxYeBgJIeAwFRfBgIKOkxNUQUHwYCSnpsERDFx9VAICKXikiRiBSLyJ0R3p8kIitEZLOIbBORy90sD/VPjANE8XEtEIhIOoCnAFwGYAaA60RkRshivwLwiqqeBeBaAE+7VR7qv3hBGVF83GwRnAugWFVLVNUL4GUA80KWUQCDzcdDABx1sTzUTzEOEMUnw8V1jwdw2PG8FMB5IcvcA+AdEfkBgIEA5rpYHuq3GAmI4uFmi0AivBZ6xF4H4AVVnQDgcgAvikhYmUTkJhEpEJGCyspKF4pKyYwtAqL4uBkISgFMdDyfgPDUz7cBvAIAqvohgBwAI0JXpKrzVXWWqs4aOXKkS8WlZMU4QBQfNwNBPoBpIjJVRLJgdAYvDFnmEIDPAICITIcRCFjlp5gEAgwFRPFwLRCoqg/ArQCWAtgNY3TQThG5T0SuMBf7CYDviMhWAP8GcIMqG/oUG+4wRPFxs7MYqroYwOKQ1+52PN4F4EI3y0D9H6sORPHhlcWU9DjFBFF8GAgo+TEOEMWFgYCSHvuKieLDQEBJj6khovgwEFDSY2cxUXwYCCjpMQ4QxYeBgJIeZx8lig8DASU/xgGiuDAQUNJjZzFRfBgIKOkxM0QUHwYCSnoMBETxYSCgpMfOYqL4MBBQ0mMYIIoPAwElPTYIiOLDQED9ACMBUTwYCCjpsUVAFB8GAkp6nH2UKD4MBJT0eEEZUXwYCCjpMTVEFB8GAkp6jANE8WEgoKSnbBIQxYWBgJIe4wBRfBgIKOmxs5goPgwElPTYIiCKDwMBJT0GAqL4MBBQ0uPso0TxYSCgpMcwQBQfBgJKfowERHFhIKCkx1FDRPFhIKCkx0nniOLDQEBJj33FRPFhIKCkx9QQUXxcDQQicqmIFIlIsYjc2ckyXxGRXSKyU0T+5WZ5qH9ii4AoPhlurVhE0gE8BeCzAEoB5IvIQlXd5VhmGoBfALhQVWtFZJRb5aH+i3GAKD5utgjOBVCsqiWq6gXwMoB5Ict8B8BTqloLAKpa4WJ5qJ/i7KNE8XEzEIwHcNjxvNR8zekUAKeIyBoRWScil0ZakYjcJCIFIlJQWVnpUnEpWTEOEMXHzUAgEV4LPWQzAEwDMAfAdQD+KiJDw/5Idb6qzlLVWSNHjuz1glJyY4uAKD5uBoJSABMdzycAOBphmTdVtV1V9wMoghEYiKLGMEAUHzcDQT6AaSIyVUSyAFwLYGHIMm8A+BQAiMgIGKmiEhfLRP0QGwRE8XEtEKiqD8CtAJYC2A3gFVXdKSL3icgV5mJLAVSLyC4AKwDcoarVbpWJ+ifOPkoUH9eGjwKAqi4GsDjktbsdjxXAj81/RESUALyymJIeGwRE8WEgoKTHKSaI4sNAQEmPs48SxafbQCAiaSLylb4oDFFPMDVEFJ9uA4GqBmCM/iE6LjlTQ7y4jCh20aaGlonIT0VkoogMt/65WjKiKDnP/YwDRLGLdvjot8z/b3G8pgBO7N3iEMWHcYAodlEFAlWd6nZBiHoqEAhNDUWa5oqIOhNVakhEckXkVyIy33w+TUQ+727RiKKjnTwmouhE20fwPAAvgI+Zz0sBPOBKiYhixD4CovhEGwhOUtWHAbQDgKq2gO1vOk4EjRpim4AoZtEGAq+IDIDZ8haRkwC0uVYqohiwRUAUn2hHDd0DYAmAiSLyEoALAdzgUpmIYuK8doCBgCh20Y4aekdENgI4H0ZK6EeqWuVqyYiiFNxZzEhAFKuoAoGIvAhgNYD3VbXQ3SIRxYapIaL4xDJqaCyAP4rIPhF5TUR+5GK5iKIW3FlMRLGKNjX0noisAjAbxq0lbwZwGoAnXCwbUVQCQS0ChgKiWEWbGnoXwEAAHwJ4H8BsVa1ws2BE0QpKDSWuGERJK9rU0DYYF5TNBHA6gJnmcFKihAuefTSBBSFKUtGmhm4HABHJA/BNGH0GYwBku1c0oigxNUQUl2hTQ7cC+ASAcwAcBPAcjBQRUcIFDR9lHCCKWbQXlA0A8BiAjarqc7E8RDELmn00geUgSlbRpoYeEZEzANwsIoBxPcFWV0tGFKXgFgFDAVGsop2G+ocAXgIwyvz3TxH5gZsFI4oWRw0RxSfa1NCNAM5TVQ8AiMhDMIaS/tGtghFFi6OGiOIT7fBRAeB3PPeD01DTcSK4RcBIQBSraFsEzwNYLyKvm8+vBPA3d4pEFBtlbogoLtF2Fj8mIisBfBxGS+CbqrrZzYIRRct57g8wEBDFrMtAICI5MOYVOhnAdgBPc/goHW+YGiKKT3d9BH8HMAtGELgMwKOul4goRuwsJopPd6mhGar6UQAQkb8B2OB+kYhiE2AXAVFcumsRtFsPepISEpFLRaRIRIpF5M4ulrtaRFREZsX6GUTKuYaI4tJdi+AMEWkwHwuAAeZzAaCqOrizPxSRdABPAfgsgFIA+SKyUFV3hSw3CMAPAazv4TZQymNqiCgeXbYIVDVdVQeb/wapaobjcadBwHQugGJVLVFVL4CXAcyLsNz9AB4G0NqjLaCUx5M/UXyivaCsJ8YDOOx4Xmq+ZhORswBMVNX/dbUiEblJRApEpKCysrL3S0pJjfcsJoqPm4Eg0pXH9mEqImkAHgfwk+5WpKrzVXWWqs4aOXJkLxaR+oOA4+wfYCQgipmbgaAUwETH8wkAjjqeD4Jxx7OVInIAwPkAFrLDmGKlnTwmoui4GQjyAUwTkakikgXgWgALrTdVtV5VR6jqFFWdAmAdgCtUtcDFMlE/xFFDRPFxLRCYw01vBbAUwG4Ar6jqThG5T0SucOtzKfUEXVCWwHIQJatoJ53rEVVdDGBxyGt3d7LsHDfLQv0YO4uJ4uJmaoioTwR3EDMSEMWKgYCSHm9eTxQfBgJKerwdAVF8GAgo6QXfj4ChgChWDASU9Jwnf8YBotgxEFDy46ghorgwEFDSC76OgJGAKFYMBJT0OOkcUXwYCCjp8eRPFB8GAkp67Cwmig8DASW94OuKGQmIYsVAQEnP2QoIMA4QxYyBgPoBZ2qIkYAoVgwElPQCnGKCKC4MBJT0lJ3FRHFhIKCkp108I6LuMRBQ0uMFZUTxYSCgpMfb0hDFh4GAkh77CIjiw0BASS/4OgJGAqJYMRBQ0guafZRxgChmDASU9IJvVclIQBQrBgJKesreYqK4MBBQ0guoQsR4zDhAFDsGAkp6CiDNjATsIyCKHQMBJT8F0uwWASMBUawYCCjpKZQtAqI4MBBQ0gsokG42CXgdAVHsGAgo6ak6WgQJLgtRMmIgoKSngD1qiJGAKHYMBJT0VB2jhhgJiGLmaiAQkUtFpEhEikXkzgjv/1hEdonINhF5V0Qmu1ke6p+M4aPmY8YBopi5FghEJB3AUwAuAzADwHUiMiNksc0AZqnq6QBeBfCwW+Wh/ktV7c5iBgKi2LnZIjgXQLGqlqiqF8DLAOY5F1DVFarabD5dB2CCi+WhfirAzmKiuLgZCMYDOOx4Xmq+1plvA3jbxfJQPxUIABl2i4ChgChWGS6uWyK8FvEoFZGvAZgF4JOdvH8TgJsAYNKkSb1VPuonAqpIs68jSHBhiJKQmy2CUgATHc8nADgaupCIzAXwfwCuUNW2SCtS1fmqOktVZ40cOdKVwlLyUscFZUwOEcXOzUCQD2CaiEwVkSwA1wJY6FxARM4C8BcYQaDCxbK4rr6lHU++uxd+Vkn7XEAV6ZxigqjHXAsEquoDcCuApQB2A3hFVXeKyH0icoW52CMA8gD8R0S2iMjCTlZ33Lt30U48tmwPVu+pTHRRUo4zNcQ4QBQ7N/sIoKqLASwOee1ux+O5bn5+X6psNLJaEqlnhFyl6uwsTnBhiJIQryzuJS1eP4COK1yp7wS3CBgJiGLFQNBL/GZVtNkMCNR3jNlHjcdsERDFjoGgl2SaZ6KWdl+CS5J6AqpITzO+f8YBotgxEPSSLDMQsEXQ91SBdHuuIYYColgxEPQSq2ug3RdIbEFSUIBzDRHFhYGgl/l4HUGfC6hCOA01UY8xEPSydj9PRH3NSA2xRZDMLnp4BZ5eWZzoYqQsBoJe5vMzNdTXnPcsZiBITodqmvHwkqJEFyNlMRD0Euum6e1MDfU55ZXFSS3AYybhGAh6iZUSYoug7xlzDRmPj9dRQweqPFhTXJXoYhyXvDxmEo6BoJdYAYCdxX3PSA0d39cRXPKH1fjqX9cnuhjHpWQJBC9+eAAL8g8luhiucHWuoVRiBYD2BO7Ub245gsrGNtz4iRMTVoZEMIaPGo+PtxbB82v248UPD6LNHFZc39yOIbmZCS7V8cWbJEOu73pzJwDgmtn9754oKdsiCAQU60uqe2191s7sS+CooR+9vAUPvLU7YZ+fKM77ETz3wYHjairwexftQkmVx35e2RTxlhsprS1JAkGiqCqW7ChHXbPXtc9I2UDw+PI9uGb+Omw+VNsr67NaBL4Ad+q+5rxncdGxRry2qTTBJepcstR++xK/k67lH6jFzf/ciCffdW94bcoGgvcKjfvgtPTSlBBWH8HxcB3B8ZYecZvzymIA2FfZlMDSdK3Nd/xNQVJc0YhjDa0J+3wGgs41e334yl8+BADUtbBF0OusKSE8vRQIjqdRQ01t3U9898bmI8f1CTMWAccFZQAS2mP8/Jr9eHzZHvt56Kzkx2MaZO5jq3Hh795L2Oc7A8HxWolxlqsvj/HNh+rsx7lZ6a59TsoGggxzlElTW3uvrM/qJE7UdQTOsdit7V3vqAerPbhtwRbc/OJGt4vVY15fAFVR5tOd1xEAHdd0JMK9i3bhiXf32s9Di3I8BgIgsaPdvP6Oytjx0KKOxPm7tfbhb1jj6WgFtPvc+25SNhBYs4U2tvbOtNF2H0EPagvNXl/cHZzOHbW79ENJpdF5ubeiKaaLedr9gT4bFfXz17Zh1gPLo/peQlsEx0NfcWc12+MtDXI8pKpi2Xf7yq3/2oQvPr3Gfu6cVbi30snRsFr3uVnpaOylSmskKRsIzAYBPG29lBoyd+bOTlx1zV5c/sT7WFlUgT+9txdvbSuz35tx91J8/6X4auet7R3b0V2ts6y+Ix+8+XD0neWX/GE1PvP7VbEXrgf+t+0oAKA2ipESAVWkpx8fLQJLZ9ORHy8nOkt9s3snl2g5g2M0gfKWlzbhqmfWulkk/G9bWVBaxuNItzqPNbc1mRXV0YNz0NZNSz8eKRsIrHNFSy/9qO2BrjuL1++vwa6yBryw9gAefWcPbvnXJgAdEX/pzmNxfX6r4wTT3Q7jTLkcqmmO+jNKKj0xLR8rVbVrWzmZRj7Uuhd0V3+joS2CBDUJnK3BTgOBiwdzT/TW/h+LQEDxu7cLUVzRCCD45B9N6uyt7WXYeLC2T4YJW4Hb+Xt+4uEVOFTt3nHg1GieH4bmZrqaVkzZQGB9qS3eXkoN+bsePmrVKPY7xpQDQEUvjdZw9gt0V+us8XhhpdSPNcQ+rt2t9NDCrUcx/e4lKK5oslN3tR4vjtS1YF0n13xYAd05asjbi3nmWLbVmTvuLH0QzcHs8wdw49/z8Z+CwwCMWvuPF2xBdRzXIDR7fVi7z5jiIhBQfP1v67FqT2W3/Unx6Gzce1lDK/68ah+ufMqo1TuvLI7lZNdb4+qrmtrw81e32ZUyZ0Wi1mO0mDwh54kVRRVRrz8QUOwua+hRBaWp1Ye87AwMyEx3tTWZsoHAat71xh3FVNVxZXHkH9sKBOX1wSd+54k4ntEIsaSGajxeTByei7zsjB4NG+ytfpVQK8whvR+WVNvpnWavH9fNX4dr56+LeCBYy2U6UkPRHjCqigX5h1DfEjk98sHeKkz7v7ex9XBdxPdDOX+D5k5uWRpN2fZXebB8dwV+93YhAGDZ7mP47+YjeNB83hM3/WMj/t+z61FjBtb391bhlpc2xZzmOFjtwf3/29VtgDxY7cGZ9y3Dcx/sD3vPCpLWidcbQx+Bs+/F2ZEaj3sW7sSCgsNYvafSKJ/jO7HK2BySQo7lpPyvDYdw2RPvY/GOsu4XDtHU1o687AxkZ6SxReCGjhaBFRB8qG5qw9Mri1FwoCamdTlP/p2dzJvMHSn0x6xo7DgRx9Nf0RZD87rG48XwgVkYNTgbFT1oETS2upNXHpRjTL1wuKbZ7vD1eH12OmpfhSfsb6zlrHtGA9GnX3aXNeLnr23HrWaaLtQH5iRxK4sqo1pfSxQditHkwBvMQFttnuiyMoxtK60NT0ccrPZge2l9t+u0tqWisRWHze8zOyMtKBBEU2N98t1i/O2D/Vi2q+tUppVf/+/m4Iv71pVU47eLg69+j6WPwDncu7qXAoFVEbAqc87+AOtxaIugxRv9SbmwvAEAsL8yfP/tTlObD3k5GcjOSHc1rZiycw2FtgjOum8ZhuZm2jX0A7/7XNTrcqaDIg3DO1zTHLHTU1WDTsQVja0YmJ2OjPTY47PzgO6ullft8WLckBzkZKSjPMoWgTMf+8lHVuKFb87GnI+MirmcXakzD8jy+la75tfi9SNNjBN+RWMrZmBw0N90tAgcgSDK2toxMwjvq4h8PYW1yjSJ+HYY5/fe4vVHrBREU6sLDbQN5vdS5+jYVVW0+QL45vP5KKnyoPD+S+1+la5UNXqDWkCtIRWIARHGqtc3t2NfVRN+8dp2HK1vAQAcrTP+b/cH8PrmI/jiWeODfoMj5vu5mcGnmOuf2xB2so8lNVTrOPn3Vosgw/yBa8zUmzPYWIGgOSQQdNbii6S01vzO6nvW+s7LzkB2ZhpTQ26wA0F7R0091ny5quJHL28OqjGGpob8AcUnHl6B+atLwv6+zRcIahF89vHV+G4Px/ZHkxr67eLd+NHLm1HraBFEO1Y/tFPxzte296icXbFyvscaWu3cv8frtwNjRYSOYyvwOi+2ibYJbXVEO0ccOVn58+ZOAmuNx4vfv1Nkf/fO76il3R+xIzaagzn0gsCG1vBAsCD/ME69a4k9j5G1LaqKjQc7HwlW2+y115eeJkEtl87Kdv6D7+JLT69F0bFGOy1ojTx7fdMR/OzVbfb+bXTeq12ehpCglhkhqgalhrqp9TpbAdUeL66bvw7f+UdBl39jeWL53qCL/SzWbU7rW8zav+P7b7QDgfHdfOGMccbzNj8WbT2Kqb94q8u+ClXFjiNGi80KjrFoavNhUA5TQ66xDtIWr6/HY7sbWnx4c8tRfP+ljtRCaC2wq+GPja2+sODzbmEFNh4MTk1tPVyHZ1bu67LpHtRZ3MmJa/7qEry55aiRGsrLwsi87G5H5QDAK/mHMfPXS4Nei2ZYZyhVxXXz1+HNLUeCanYWq4ZX0dgGv9VH0Oaza2xWWQMBtTtOrWG7VvoEiD41ZAeC0Mt/TVbNuaGTPoTn1+zHH98rxoJ8o1PXeVJt9kYOBDuONHQ74iS0D6bBPEFVNbXZ+9frm49ELOt/Nx3BVc+sxeLtZWjz+cM6NRtbffay6WkSdPJ/Y/ORiK3JSNthtQgOVBuByBr0cPebO3Hl02vtCkaFI0C9urE04pX8zhOc8+KySIJaBE1efFhS3W2ayvL48j1BF/tZQtM/kVJDVh/Bg1/6KMYPHYBmrx+/eWs3VIHtRzpPzZVUeVDVZJT5SITUXneszuIsBoLel3+gxj5xtrT7I3YWRnPRSKSTYWhqqKvmq6fNhyN1LfaJznLVMx9i1Z6OVsa8p9bgoSWF2ObY4Xz+QFAKwXlAt/oCYRc0OZ97/QEMz83CyEHZaPb6g3b8SJ5ZtS/stTQRBAIaUxCtaGzDhyXV+NHLW3DW/cvCaqDWQV5e32rv9E1tPvvkZJ1snnh3L855YDlqPF572G5PUkNWIOhsmhGrBt5Z53hReWPQeoJGDbX70Rohj7xqTyUuemRFl+Vqcnxeuz9g16p9AbVr4ukh+8wzK/ehxetHmZm62XSwFo8t24NvPp+PfEefV0Nru72/B1SDTvz3LNqFB97aFbTe0MENFivdYb1/pK4Ff32/BC+uO4ith+uwxexgr/F44fUFsLW0Hj/9z9aw9VgpLktbewD1ze3YXdaA+xbtwvNr9gdVgIJbBB2VGGs7VBWvby4N2wecqc3QNI/1+zZF6A8IDRIDMtORm5WOlnYfBg8w0l6Hazqv6f997QGIAJ87fSz2VXqiboFbmtrM1FBGeqcVvN6QMoGgvL4Vb20rQ2NrO7785w/t15u9kQOBcyeLZPWeSjy0JHgUR05mWthoiuqmzgPBL1/fjo0HazF97OCw9/ab8wA5A5Jz9MrPX9uOj97zjn2QOA/ou97YgXsX7cKWw3XYebTefD+4XMMHGoEA6H6sfqQLtNIEuOvNHfjY796LerTTwZCa8O6yxqDnteaJt6Xdbx+4ZfWtdoewVbv81wbj5iCbD9XaqbisoEAQW4ug1uONeCWw1RJYuPUonloRPPNjfXM73t9rdMBa/SzO36qsrhU3vdh5ysIfUPgDGrHC0egIzHXN7UEtEus7DJ3u/K3tZfjD8o60R0CBUvME5ewDaWxttwOup80ftl+E/iZ/jlAJAIxa8O6yBjswLd9dETQFuhUoAKMlc6Aqckepx+sPu47gF69vw2VPvI/n1uzHvYt24Yz73sFjZkrHKvuIvGwccOxP1m+5oqgCty/Yij8sD675Oytkoa1wK9D+a/0h3PSPgqBBG01tfry+uRTPfbAfOZlpSE8T5Galw9Pmt2e8daZ3azzeoM86WteKU8cMxvknngAAuH3BlojfQ2c6OovZIugVL60/iFv+tQn3Lgqu8bR4/aiPMKtfdx1R1z+3AW/vKA96LS87M+wA7Wo9a/cZY+NnTxke9p41csQ5UqToWCP2HGvEo0uL7KmWrRxx6DDYF9YewJVPrcHnnvwAAMKC3Ql5jkBgpVn8Afziv9vw7u7gpnZ7hB0wTQQvrT+Eqqa2qC8yC02xOFMRrWZO/SOjBwUt41y3FQisuvAdr26zb/+YETR8NLoDxgr2voBGrAw4X3tkaZGdDgGMIZ0t7UZHtjUE1xmM//pBCQrLg0+qThWNrfjBvzfh3N8st//uSF0LPtxXHdTSq2v2oqHVh1Hmb2V9H5E6K8vqW+00REu7H9lmusz5HTa0+Ox90uP1haV9MtIE7xUesys0kUYqTR0xENkZafjTiuIuBxtMGDYAgHGSdp4snXYdbQgKNrct2ILF24OPq8ZWH558dy/a/QH8ZXUJMtIEU0fkorCswV7GWr8VEJ5ZuQ+/eWuXHeCdNfHQVo5zv3xn1zEcrO4IWp42H25fsBUNrT4MzDJaAAOy0tHi9dstCGf66+z7lwVNTVHX7MWw3Ex8+ZwJABC0Dz2zch/+vvZAxO/FWl9Tmw+DzBaBL6CuTXiXMoHgtrmnYNbkYXh1Y8dwthF5WWhoae+kRWAcLK/kH8YNz2+I6irGIQMygq7wNdbTfVPw2nMnhr1m7dDWQZyRJthT3oiHlxThT47a6b7KJlw7/8Mub0hT3dQWNoXtsNzwFsGG/TX494bDuNu8ExNg5OMj3UzFWWstq2/F798pwoNvd31TnNC5Ul5YewDPmp2MVprtI2MiB4LRg7Ptg91Kv9V4vPjZq9sABKeGQmvZz64uwZNmbrjZ68MX/vgBlu86FtRaq4rQcgvdL5y12jXFVRiRl4W500fbJxbnSbXbkVtNXizeXo7GNh8KDhiduz95ZQuue3Ydih01+FqzRXDK6EHISk/DwRoPVBVldeEn1vKGVvuEV9fstSdAfHplx4m2sbXdDgSqsK9VsKzfX4NvvVCA1zcZfRAHqpsxc3xwi/WS08bgspljsL6kxk5FRXLO5GEAjBNlZ61Oa4rlaDy6tMjoJwkoJp8wMGjwQEWD8RnbHENpn31/v93KdAYC57UzgYCGdc4/+o7R+hiYlR6UirImNszNykBzu89OHVoj/6xlD1Y32wGottmLobmZyMlMx9XnTMC+Sg++98+N8AcUDy0pxK8X7kRheZTHBOwAABUuSURBVEPEa1WavX6owmgRZBr7t1u39UyZQJCeJvjsjNFBr40bOgAerx/l9eE76Tefz8fGg7X42WvbsLKoEvNXl9jDGkNzjJYTBho5902HavHg4t3wBzToZDM0wi0KV/x0DiYOyw17vbKxDccaWvE/c06iS2aOwe6yBrtzzrJo61GsK6mBP6BhUx5bisobw+aUOWFgNkbmdQSC3WUN9j11ncGrorGt2xkhj9S24I/vFeMvq0rCAma7PwB/QHG0rgWrIozH31fZhLrmjub0qWODA4H1+vSxg1HR0IZAQCN23jovKHOOVGnx+vGbxbvx2LI92HOsEVsO1WH7kXrc+I8C1Hi8mDpiIACgpLIp7MQfFgiqm/H65lLMfWwV1pVU4+xJwzBu6ACU17ei3R9A/v6OXHx339nW0o4D3xrjbwUEK+UEdIzyGZqbiQnDBuAvq0rwpWfWRhxDf6S2xT7h1TZ7I16J3NDqQ00UHf27yxvgDygOVTdj1uSOFus/v30ebv/sNJx/4gmoampDa3sAI8z9yPL18ycDgH28GS2C+O/M9tZ241i4b95p9u9mOVLXgs89+T5eWh98T+GOzuzIgaC+pb3TSQoH5WSixDFVu9VfkJuVjrrmdjuAHK5pxj0Ld+INRwd+vWPI79DcLACwW3Vv7yjH/qqO9V76h/cx76k1qG9uxzMr99lXgFvrz8vOtFt3bl1LkFLXEUwbnRf0fMbYwdhWWo+i8oaIyzsntnpoSSHeKzyGK84YZ9+71DJkQCbqW9oxftgAbDhQgy89bfzd508fF5QaGjM4J2gIIABMHp4b8QRe2dSGm17caNcUPn7yCLy1rSyotgjADhRA+JTHlpIqj137twwbmImBWRlITxNUNrYFjS5p8wVQUtmEE0fmdZn2GZSdgcY2n32AAsD+qiaU1rZg9pThaG33GxPVnTraqPVEuPDpjS1H8caWo/bzUx0tgqyMNDt/PGPsYKwsqkRheWPEazWcLYJmrx+PLC1EdZMXXz1vsv366j2VdicmYNTgZk8Zjv1VHtxkDtstvP9SLNp6FOdNPSEsbXKw2oPn1x6wy3TFmeMwdEAWGtt8uPvNnfiveSIYOyQnaGK/SNYWG2nBrPQ0/OPDA1i9pzJou04cORAllR68u/sYSio9Ro7Z3E+ck6EBwCemjYAqsHZflR0Q15WEXxRptYBrPV77t+vMziMNKKtvgdcfCGqlfXzaCADAeWbOGwDOmDAE7xZ27D93f2EGvvGxKZh8grFvVzS2RjU6rTultS04bdxgXH/BFCzeHnyVbmct4h1H6nHauMH20M2MNAlKZ3V2Udpr3/sYfvbq1qDbjFrp19ysdLsPJCNNUHSsEUXHgtOAR+taMWRAJupa2jHMrAA6j0FnsLfcs2inPRps/4OX253YA7NjHxodK1dbBCJyqYgUiUixiNwZ4f1sEVlgvr9eRKa4WZ6TRnYEghU/nYMrzjTGBK/fH3zQDHPU3Ac6xqfnH6gNCwIA8Oz1s/DDz0zDaeOCm9BbSuuCAsHE4bkYP3QAnrj2TGRnpGH62MFISxN7HLPTxoO1Qc3FE0NqQE4XmAdlZifj4X+7eDceWVoU9FpedgbS0gQj8rJQ2diGHUfqMWl4Lh6++nSoAp/+/Sos2VFu7+Dfm3MSZo4fjJU/nYObP3kSgI4d2znC6fFle3HD8/n42avb8EpBKaqavFhQcDgsCJwxcSjGDckJK+sER+tounkCSk8TnGL2HVz+5Pudfg9OT63Yh5fzD9t3ogOMk4UzcALhqajfvV2IO17dhmvnGymLbMew1F1lDUEdm5OHD8SYIcZ38Io5LxCAsKBrefBLH8UPP30yAGD9fiMQXH/BZDR7/dhVFlwZGW7WIl8pMFKZg3MyMXPcEPt9535529xTMO/McQgogjpQnU4cORDnTB6GA9XNqGtpx8Th4a1Qpw0HanDPQqM/bcoJ4fvelBM6/v6SmWOQmS546KqPYvs9FyMzPQ0nj8pDZnoaxgzOQVF5Y9yBYOJwo79hklnu0BYBgKCWySfMgHX3mzvxkV8twWPL9mDUoGxMOiEXz685gF/8dzv+vGoflpjTPgwZENxaP2V0HvJyMiOW26rhA8D355wUsbxffHoNHlpSBH9AMcxuEXTs7//ecCjsb5xDgj8orsLn/2js69Z1BADCKoK9xbVAICLpAJ4CcBmAGQCuE5EZIYt9G0Ctqp4M4HEAD7lVHiD4JDN1xECcNnYIMtIEheWN9o4GACt/+ikMyjEaSw98cWbQOgZlhzeizp06HD/+7ClhV2W+vOFQUG15cE4m1tz5acw7czw2/HIuXv/+x+z3Hr7qdPz4s6dELPc3LpiMqSM7dvxLTxsT9P7nzxiL/9x8Ad77yZywvx2RZ6SrQncgK/iMHJSNo/Ut2FZaj4+OH4LJjhPEzf/ciLve2IHcrHTccfFH8L8ffAJTRgzEyaOMgDowOwOPXH160Hqt7X1rexkeWlJoN4edbvnUSXjpxvNw0qi8sPdOGNhxkM0wT3xjBudg1ODIJ1dLZ53yT68sRnZGGi52pAXvvOxU+/GpIYHgBbPzzroK1Hm1bmgtbsoJuRg92Di4nSmxkXmRy/rFs8bjdvM3rmryYtSgbFw6M/i3vPmTJ2Hu9FG454rTgl4fOSgbv7x8Oq6ZZfQnOW+ANGRAJqZ0UlH47kUnouS3l+O9n8zB7CnDUdXUBlVjOGNnrArFcnPQwJQRuVh2+0X48BeftpcREdzyqZNw8qg8XHnmeBTefxmumT3JnibEcvakYXh7Rzn2RnECs1JKkZxsVuImm0EpUnD613fOw4kjB+K2udPw4rfPw29Cjt2PjBlkB/J/bziE371daPcHTAoJjINyMjEgM/Lp0blPfPX8yfj5paciJ2TZNl/A7gS3KgbOix73HAv+Pr7zialBz294Pt8e0TVkQCYmm4F3V1n304n0hJstgnMBFKtqiap6AbwMYF7IMvMA/N18/CqAz0ik6nEvCR17PSQ3E9/42BQAwOkThmJQTgZGD87GkNxMuwZ+5sRh+M/NF2DpbRfhle9egIK75uL6C4wdNis9DXOnd5xgrBPA9+echLnTR2PnUaOWZ7Uw/t95k4I+27lDfWX2RHzPrF3c+PGpmDt9FJ7+6tnY88BluOeK04JOLo98+XQUPXCpXSv67PTRmD1leFAtb9Udc/DSjedh8Y8+br921+dD4zDwsZNG4P29VSitbcEZE4fglNGDMDQ3E9NG5dlTK0wanht0BzDroBk1KBtfnjURf/n6OfjGBZPx8ZONWtiVZ47DxOEDMHZIDpbcdhE+ecpI5GVn4DOnjsK3LpyKOy45FXnZGWG5ZcAY1mrVfi4+zfhuPzN9VND2D8hMxwNXzkROZpq97JghOSj41Vxsu+di3GD+poBxQJ40Mg/XnTsJg3Iy8PJN5+PmT55kB4AzJg61g7t1UdrtczsC8g8/Mw0AcNNFJwIw0j6WGeMGY+yQjgqEZeZ4I4D94ZozseeBy+zglpOZDhGx94erzpmAMyYOxVmThtp/e9XZ4/HXb8zGzPFDMGvyMIzIy8LsKcNw8YzRGDMkB9+5yDhhjBmcg8s/OsYu05mOFtbDV59u7+tjhuTYv91XZncMSpg7fXTQiR3oGIJ72cyxuG9eRyAaPSgH00YPCtvWOy45Fct//ElkZaSFHVuW++adhgtPNlqsoUH3wO8+h7V3GmWYOHwA7r9ypv1d/Plr5+CXlxsB++Mnj7ADwHSzD2lAVjpu/PhU3HHJR/DkdWfha+dPwrRReXjvJ3Nwm/n7ffW8yfjS2eMxdkgOsjPScNXZE+yWwl2fn4Hb556Cq86egD9ccyY+dpJRxjQBvnTWePOxsU3ZGWk4Y+JQPHHtmQCMitmtnzoZG381F6MH5+B7c07CG7dcaH5mxzFusSqgF50yEvdfORMnjgwPYtfMNv7ulNF5mHxCblDF4sQReZg1ZTj2PHAZbroocgskbtYl4b39D8DVAP7qeP51AH8KWWYHgAmO5/sAjIiwrpsAFAAomDRpksZj2c5yXb6r3H7u8wf01YLDeqyhRY/Vt2hTa7uqqpZUNukr+YcirqPd59cjtc3a4vWpzx+wXw8EArrpYI36/AEtKm/QX72+Xe9duFN3Ha3XQCAQcV2xeHd3ub6xudR+fqjao6uKKoKWKSxr0M2HaoNee3v7UV1TXGmvY++xBvu9umav3rNwh96zcIc2tHjt7VBV9fr8+vulhbrpYE3Q+nz+gP5+aaGW1bUEvX6gqknvW7RTjzW0aH2LN2h9LV5f2PbUNLXpo0sL9YU1+/X1TaW6wPy+91U0alF5gwYCAd17rFH9/oAGAgF9dvU+La1tttdr2VPeEPS8srFVX1izX4vKG/TXb+7QDfurVdX43SzH6lt02+E6VVVdX1KtC7cc0cKyBt10sEa9Pr8+vqxIn1qx1/4umtt8+vulhbq9tE73lDfoyxsO2tv2l1XF+uzqfXqo2qN7yhu0rd2vK4sq1G/uG2V1LZpvlkHV2LdCv9O9xxr0tY2Hg17z+vza2h78vfn9AX10aaEeqGrSOo9XD9d47Pc8be26Zm+lBgIB9fkDmr+/Wr2ObVZV3XGkThdsOBT0/GCVR5ta29XT1q4PLt6tlY2tqqq6YX91ULl7qs7j1fsW7dQ6j1fL61v0ULXH/gxVY5+t9bSpavhxV3CgRisbW/VQtUefXL5H29r9YeuPhnWctrb7wvZb6/VX8g/Zv5mqsT8/vqxID1V7wpaPZE95g3ra2vXt7Ud1/qp9+szKYn14ye6w36DW06br9lXpmr2VumjrEX2v8JiqGuemfRWNuvVwrT6ypFA/2FupT63Y26PtjQRAgXZyvhZ16W5OIvJlAJeo6o3m868DOFdVf+BYZqe5TKn5fJ+5TOTJ5wHMmjVLCwqim1uEiIgMIrJRVWdFes/N1FApAOcA+QkAjna2jIhkABgCILY5oImIKC5uBoJ8ANNEZKqIZAG4FsDCkGUWAviG+fhqAO+pW00UIiKKyLXrCFTVJyK3AlgKIB3Ac6q6U0Tug5GrWgjgbwBeFJFiGC2Ba90qDxERRebqBWWquhjA4pDX7nY8bgXwZTfLQEREXUuZKSaIiCgyBgIiohTHQEBElOIYCIiIUpxrF5S5RUQqARzs4Z+PABA+7V//xe3t31Jpe1NpWwF3tneyqo6M9EbSBYJ4iEhBZ1fW9Ufc3v4tlbY3lbYV6PvtZWqIiCjFMRAQEaW4VAsE8xNdgD7G7e3fUml7U2lbgT7e3pTqIyAionCp1iIgIqIQDARERCmuXwYCEblURIpEpFhE7ozwfraILDDfXy8iU/q+lL0niu29QUQqRWSL+e/GRJSzN4jIcyJSISI7OnlfRORJ87vYJiJn93UZe1MU2ztHROodv+3dkZZLBiIyUURWiMhuEdkpIj+KsEy/+X2j3N6++X07u3VZsv6DMeX1PgAnAsgCsBXAjJBlvg/gz+bjawEsSHS5Xd7eGxBym9Bk/QfgIgBnA9jRyfuXA3gbgAA4H8D6RJfZ5e2dA+B/iS5nL23rWABnm48HAdgTYV/uN79vlNvbJ79vf2wRnAugWFVLVNUL4GUA80KWmQfg7+bjVwF8RkQi3337+BfN9vYbqroaXd/Fbh6Af6hhHYChIjK2b0rX+6LY3n5DVctUdZP5uBHAbgDjQxbrN79vlNvbJ/pjIBgP4LDjeSnCv1x7GVX1AagHcEKflK73RbO9AHCV2ZR+VUQmRni/v4j2++hPLhCRrSLytoiclujC9AYzXXsWgPUhb/XL37eL7QX64Pftj4EgUs0+dIxsNMski2i2ZRGAKap6OoDl6GgN9Uf96beNxiYYc8icAeCPAN5IcHniJiJ5AF4DcJuqNoS+HeFPkvr37WZ7++T37Y+BoBSAs8Y7AcDRzpYRkQwAQ5C8ze9ut1dVq1W1zXz6LIBz+qhsiRDN799vqGqDqjaZjxcDyBSREQkuVo+JSCaMk+JLqvrfCIv0q9+3u+3tq9+3PwaCfADTRGSqiGTB6AxeGLLMQgDfMB9fDeA9NXtmklC32xuSQ70CRi6yv1oI4HpzdMn5AOpVtSzRhXKLiIyx+rdE5FwYx3R1YkvVM+Z2/A3AblV9rJPF+s3vG8329tXv6+o9ixNBVX0iciuApTBG1DynqjtF5D4ABaq6EMaX/6KIFMNoCVybuBLHJ8rt/aGIXAHAB2N7b0hYgeMkIv+GMZJihIiUAvg1gEwAUNU/w7hH9uUAigE0A/hmYkraO6LY3qsBfE9EfABaAFybxJWaCwF8HcB2EdlivvZLAJOAfvn7RrO9ffL7cooJIqIU1x9TQ0REFAMGAiKiFMdAQESU4hgIiIhSHAMBEVGK63fDR4kiERE/gO2Ol65U1QMJKg7RcYXDRykliEiTquZ18X6GOe8UUcphaohSlnmfhv+IyCIA75iv3SEi+eYEffc6lv0/854Py0Xk3yLyU/P1lSIyy3w8QkQOmI/TReQRx7q+a74+x/ybV0WkUEReclw5OltE1poTjG0QkUEi8r6InOkoxxoROb2vviNKDUwNUaoY4Lh6c7+qftF8fAGA01W1RkQuBjANxtTeAmChiFwEwAPj6vOzYBwzmwBs7Obzvg1j+oPZIpINYI2IvGO+dxaA02DMkbMGwIUisgHAAgDXqGq+iAyGcSXpX2FcCX6biJwCIFtVt8X1TRCFYCCgVNGiqmdGeH2ZqloTDl5s/ttsPs+DERgGAXhdVZsBQERC566K5GIAp4vI1ebzIea6vAA2qGqpua4tAKbAmAq9TFXzAWOyMfP9/wC4S0TuAPAtAC9Eu8FE0WIgoFTncTwWAA+q6l+cC4jIbeh8qmMfOlKsOSHr+oGqLg1Z1xwAbY6X/DCOQ4n0GaraLCLLYNyQ5SsAZnWzPUQxYx8BUYelAL5lzg8PERkvIqMArAbwRREZICKDAHzB8TcH0DGt99Uh6/qeOc0wROQUERnYxWcXAhgnIrPN5QeZU6QDRnroSQD5jtYLUa9hi4DIpKrviMh0AB+a/bdNAL6mqptEZAGALQAOAnjf8WePAnhFRL4O4D3H63+FkfLZZHYGVwK4sovP9orINQD+KCIDYPQPzAXQpKobRaQBwPO9tKlEQTh8lChGInIPjBP0o330eeMArARwqqoG+uIzKbUwNUR0HBOR62Hcx/b/GATILWwREBGlOLYIiIhSHAMBEVGKYyAgIkpxDARERCmOgYCIKMX9f83bVsD2S7JZAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"period=1\n",
"\n",
"\n",
"rand = np.random.RandomState(42)\n",
"t = 100 * rand.rand(100)\n",
"y = np.sin(2 * np.pi /period* t) + 0.1 * rand.randn(100)\n",
"dy = 0.1 * (1 + rand.rand(100))\n",
"\n",
"def lomb_scargle(t,y,dy):\n",
"\n",
" frequency, power = LombScargle(t, y, dy).autopower()\n",
" plt.plot(frequency, power)\n",
" plt.show\n",
" plt.xlabel('Frequency')\n",
" plt.ylabel('Power')\n",
" return 1/frequency[np.argmax(power)]\n",
"lomb_scargle(t,y,dy)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use the frequency with the most power to find the best period to represent the signal in the data. In the code below, we put the LS algorithm together with the phase-folding code to plot the phase-folded data with the best period."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"period=1\n",
"\n",
"rand = np.random.RandomState(42)\n",
"t = 100 * rand.rand(100)\n",
"y = np.sin(2 * np.pi /period* t) + 0.1 * rand.randn(100)\n",
"dy = 0.1 * (1 + rand.rand(100))\n",
"\n",
"def lomb_scargle(t,y,dy):\n",
"\n",
" frequency, power = LombScargle(t, y, dy).autopower()\n",
" return 1/frequency[np.argmax(power)]\n",
"\n",
"def plot_best_period(t,y,dy):\n",
" phase_fold(t,y,lomb_scargle(t,y,dy))\n",
" \n",
" print(period)\n",
"\n",
"plot_best_period(t,y,dy)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us apply the Lomb-Scargle Periodogram to real data. First, let's take a look at unphased Cepheid data using the code below. (These data are from the American Association of Variable Star Observers.) Please complete the code block below to plot the unphased data. Notice how these data span a very long time period. You may alter the code block in order to examine shorter time periods. Do you have an educated guess about the period of this Cepheid?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data=np.loadtxt(os.path.join(data_dir,'cepheid.csv'),delimiter=',',skiprows=1)\n",
"t=data[:,0]\n",
"y=data[:,1]\n",
"#complete plotting the raw lightcurve; please check data file to label axes appropriately"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we apply the LS algorithm to the Cepheid data. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data=np.loadtxt(os.path.join(data_dir,'cepheid.csv'),delimiter=',',skiprows=1)\n",
"t=data[:,0]\n",
"y=data[:,1]\n",
"#uncomment code below to continue\n",
"#def lomb_scargle(t,y):\n",
"\n",
" # frequency, power = LombScargle(t, y).autopower()\n",
" # return 1/frequency[np.argmax(power)]\n",
"\n",
"#def phase_fold(t,y,period):\n",
"# phases=np.remainder(t,period)/period\n",
"# phases=np.concatenate((phases,phases+1))\n",
"# y=np.concatenate((y,y))\n",
"# plt.scatter(phases,y,c='k')\n",
"# plt.xlabel('Phase')\n",
"# plt.ylabel('Flux')\n",
"# print(period)\n",
"\n",
"#def plot_best_period(t,y):\n",
"# phase_fold(t,y,lomb_scargle(t,y))\n",
"\n",
"#plot_best_period(t,y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, let's apply the LS algorithm to the 1856 Ruzena data in order to see which frequencies the algorithm says have the most power. This should tell us the rotation period of the asteroid. Please complete the code in the block below in order to plot the periodogram for the 1856 Ruzena data. It would be most useful in interpreting the periodogram if the \"best\" period was outputted in hours, since astronomers are used to discussing asteroid rotation period in hours.\n",
"\n",
"Once you get the outputted best period (in hours), when you put it into the phase-folding code that we used for the asteroid above, does it yield a good phased lightcurve for the asteroid? BE CAREFUL HERE: should an asteroid's lightcurve have one or two maxima (and minima) per rotation? And does the LS algorithm output take that into account? Human interpretation is always necessary in science."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"data=np.loadtxt(os.path.join(data_dir,'asteroid1.csv'),delimiter=',',skiprows=1)\n",
"t=data[:,0]\n",
"y=data[:,1]\n",
"\n",
"#complete the Lomb-Scargle code for the asteroid data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Challenge: The above data sets were relatively straightfoward to work with using the Lomb-Scargle Periodogram because the data were rather \"clean\" and did not exhibit much power at alternative frequencies. You can find two additional data sets in the data director for this module. One is \"asteroid2\", which is for the asteroid 4404 Enirac observed by student at the University of Maryland in Spring 2017. The other is \"delta-scuti\", which is for a Delta Scuti-type variable star. You will find that both of these objects have very \"messy\" periodograms, in which is its unclear what the \"best\" period is. You can attack each of these objects by lookig at subsets of the data (and running LS on the subsets), potentially temporarily eliminating data points which are outliers, and otherwise being clever in your analysis You can begin your analysis in the space below.\n",
"\n",
"You can also download your own data for many variable objects from the American Association of Variable Star Observers' website, aavso.org."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}