{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Light Curve Analysis Module\n",
    "\n",
    "**Lecturer:** Melissa Hayes-Gehrke<br>\n",
    "**Jupyter Notebook Authors:** Melissa Hayes-Gehrke & Cameron Hummels\n",
    "\n",
    "This is a Jupyter notebook lesson taken from the GROWTH Winter School 2018.  For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.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 <module>`.  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": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.stats import LombScargle\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,'Flux')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEmJJREFUeJzt3X2sZHV9x/H3x8UHfARlgwisi5W20lqv9ApYWnOrmIA1rK0PRW2LBrIxlWptmxZrotWkqSamVltiJUrFxogUH1grKSq62sRKd9EbeSp1xSq7oqAI1vqIfPvHDDAu9+6evXdmzpk571dys3POnL3ne+bMPZ/z+/3OnElVIUlSE/druwBJ0uwwNCRJjRkakqTGDA1JUmOGhiSpMUNDktSYoSFJaszQkCQ1ZmhIkho7qO0Cxu2www6rzZs3t12GJM2Uq6666ltVtXF/y81daGzevJmdO3e2XYYkzZQkX22ynN1TkqTGDA1JUmOGhiSpMUNDktSYoSFJaszQkCQ1ZmhIkhozNCRJjRkaUsctLS2xtLTUdhkSYGhI0gHpe4gbGpKmru8H3llmaEiSGjM0JEmNGRqSpMYMDXWCfdzSbDA0dA8P3JL2x9CQJDVmaGiqbM1Is83QkMbEQFQfGBodNs6DkAe07pvUPnLfN+Pr1IyhIWnuGQjjY2hI6gwP7t1naEir8AAm3ZehIUlqzNCQesqWlNbC0JB0wAyc/mo1NJJckOSWJNes8nySvC3JriRfTHL8tGtcC/+gJM2rtlsa7wZO3cfzpwHHDn+2Am+fQk1S73iiM32z+pq3GhpV9Rngtn0ssgV4Tw18DjgkyRHTqW5+LC0tsby83HYZkuZA2y2N/TkSuGlkevdwniSpBV0PjUaSbE2yM8nOW2+9te1yJGludT009gBHj0wfNZz3M6rq/KparKrFjRs3Tq04Seqbg9ouYD+2AeckuQg4Ebijqm5uuSZpqrZv3952CZ1190DyNF+jvu+PVkMjyfuAJeCwJLuB1wH3B6iqfwQuA54F7AK+D7y0nUolSdByaFTVC/fzfAEvn1I5kqT96PqYhiSpQ7o+piFpDvV9XGCWGRq6h3/I6+Prpz6we0qS1JgtDUlzz1bg+NjSkCQ1ZmhIkhqze0qdYPeBNBsMjR7Yvn37VO/b38atHebBpF4v90M3zep+MTSknho9aM3qAUzTZ2hI6owDDS+/XGz6DA1NlWe00mwzNDrMA+z8mMQ4j2NHaoOhIemAdSWoFhYWxva7urJNXWdo9IR/EJLGwQ/3aWqWlpameumvpPEzNCbMA6W6wveixsHuKWkVdulJ92VLQ5LUmC0NTZSXhWqS/HDf9NnSOED2C4+fr6k0O2xpTIBn1ZLmlS2NffAMWJJ+lqEhSWrM7imNnd1zmiQvrmiXLQ1JUmOGhiSpMUND0lh5Acl8MzRmkH+UktpiaEgt8gRA+9LF90eroZHk1CQ3JNmV5NwVnn9JkluTLA9/zm6jzlnUxTebpNnX2iW3STYA5wHPBHYDO5Jsq6rr9lr0/VV1ztQL1MQZatLsafNzGicAu6rqRoAkFwFbgL1DQz3mNfmr8zVRG9oMjSOBm0amdwMnrrDcc5M8Dfhv4FVVddPeCyTZCmwF2LRp0wRKldRF4/yOcDXT9YHwjwCbq+pXgI8DF660UFWdX1WLVbW4cePGqRY4KxzjkDQObYbGHuDokemjhvPuUVXfrqofDSffCfzqlGrTDDMgpclps3tqB3BskmMYhMUZwItGF0hyRFXdPJw8Hbh+uiWuX9+/JMZ+d2m+tBYaVXVnknOAy4ENwAVVdW2SNwA7q2ob8IokpwN3ArcBL2mr3rWyz1XSPGn1LrdVdRlw2V7zXjvy+NXAq6ddlyRpZd4aXZoCu+k0L7p+9ZTm2PLycu/HfKRZY0tD6glbOxoHWxqSpMYMDUlSY3ZPaWrsHpFmny0NSVJjhoYkqTFDQ+oQ75ulrnNM4wDZLz9dvt6zx3023wwNtc6DjDQ7DA1JM2X0JMMTjulzTEOS1JgtjX3wLEbg95SrPV18zxkaao3fNSLNHrunJEmN2dKQWtTF7gdpXwyNGeSBRlJb7J6SJDVmaEiSGjM0JGnMJn0PsTbvUWZoSJIacyB8wtoatHawXNIkGBo9MWsh4qewpW4yNNRJy8vLa/6/Bo00OY3GNJIct8K8pbFXMyf8Ih1J86rpQPjFSf4iAwcn+XvgbyZZmCSpe5p2T50IvAn4LPAw4L3AyZMqSv1gN5I0e5qGxk+AHwAHAw8CvlJVd02sKqln7M7UrGjaPbWDQWg8BfgN4IVJ/mW9K09yapIbkuxKcu4Kzz8wyfuHz1+ZZPN61ylJWrumoXFWVb22qn5SVTdX1RZg23pWnGQDcB5wGnAcgyDae8D9LOA7VfV44C0MusgkSS1pGhq3JNk0+gN8ep3rPgHYVVU3VtWPgYuALXstswW4cPj4EuAZSbLO9UqS1qjpmMZHgQLCYEzjGOAG4JfWse4jgZtGpnczGHBfcZmqujPJHcCjgG+tY72SOsgPdM6GRqFRVU8cnU5yPPCHE6loDZJsBbYCbNq0qeVqNG88iHXbaNgYPJO3phsWVtXnuW+r4EDtAY4emT5qOG/FZZIcBDwC+PYK9ZxfVYtVtbhx48Z1ltUNfkBQUhc1amkk+ZORyfsBxwNfX+e6dwDHJjmGQTicAbxor2W2AWcC/wE8D/hkVdU61ytJWqOmYxoPG3l8J4Mxjg+sZ8XDMYpzgMuBDcAFVXVtkjcAO6tqG/Au4J+T7AJuYxAsGhOb8tPl66150HRM4/WTWHlVXQZctte81448/iHw/EmsW5ImZZ5PDPYZGkk+wuCqqRVV1eljr0iS1Fn7a2m8eSpVSALm+wxV82F/ofGVqvraVCrpMPuipfni3/Ta7e+S2w/f/SDJuga+55GXxUrqm/21NEZv2fG4SRbSBeM6+/DsRdK82l9Lo1Z5LEnqof21NJ6U5LsMWhwHDx8znK6qevhEq9N92BcrzaZ5+dvdZ2hU1YZpFSJJfTcLwdL0E+GSpI5YXl5ubd1rumGhJKmfDA1JUmN2T6mTbr/99rZLkLQCQ0OSZszCwkJr67Z7SpLUmKEhSWrM7impA7p8Xb40ytCQ1AlrDU4Dd7rsnpIkNWZoSJIas3tK0tywq2ryDI2O8s0vqYvsnpIkNWZoSOoEvz55Ntg9tQ52IUnqG0Ojga6Hw1q+uGVpaYnl5eVW72EjafYYGpI0I7rQfWdoSFPS9Rar1IShIal3DPC1MzRGzMIbaZw1LiwszMQ2S+oOQ0OSOmIWTuJa+ZxGkkcm+XiSLw3/PXSV5X6aZHn4s23adUrqBz8j0lxbH+47F7iiqo4FrhhOr+QHVbUw/Dl9euW1a9bewLNWr6S1ays0tgAXDh9fCDynpTokSQegrTGNw6vq5uHjbwCHr7Lcg5LsBO4E3lhVH15poSRbga0AmzZtGnetnbe8vNx2CZLGbC0f2p2GiYVGkk8Aj17hqdeMTlRVJalVfs1jq2pPkscBn0xydVV9ee+Fqup84HyAxcXF1X6XJGmdJhYaVXXKas8l+WaSI6rq5iRHALes8jv2DP+9Mcl24MnAfUJDkjQdbY1pbAPOHD4+E7h07wWSHJrkgcPHhwEnA9dNrUJJGqPt27d3rqtpLdoa03gjcHGSs4CvAi8ASLIIvKyqzgaeALwjyV0Mwu2NVWVoSOq9NsOnldCoqm8Dz1hh/k7g7OHjzwJPnHJpkqR98EuYWuTnGyTNGkNDktSYoSFJaszQ0Eyzi0+aLkNDkg5Qn09WDI116PMbR1I/GRqaG9MM8TZOGMa1Tk92tB6Gxhj5xyhp3hkaknrDE7v18+tee2oe7oEjafpsafSYZ12SDpQtjQ6yFTB+Xf1CG2nWGBpzYGFhoe0SpEYM79lnaMwB/wAlTSuQDQ1JneDJz2wwNKZstbOBWW62z2LN0izqwt+aV09J0pj04YpEWxpST8z7wUzTYWhIK5jl7kLNh66+9wwNSb3X1QN0FxkaknSA+hwyDoRLkhqzpbEOfT7bkNRPtjQkTdXy8rJXcs0wWxqaacvLy62s11am+sqWhiSpMVsaY+TZZ7t8/aXJMzRa5EFuenytpfGwe0qS1FgroZHk+UmuTXJXksV9LHdqkhuS7Epy7jRrlCTdV1vdU9cAvwO8Y7UFkmwAzgOeCewGdiTZVlXXTadEqd+8/9ZsmdZ+aiU0qup6gCT7WuwEYFdV3Thc9iJgC2Bo6B5+1W0zS0tLLC8vs7CwYAhoXbo8EH4kcNPI9G7gxJZqkVrnwV5dMLHQSPIJ4NErPPWaqrp0zOvaCmwF2LRp0zh/tSQ11odgn1hoVNUp6/wVe4CjR6aPGs5baV3nA+cDLC4u1jrX24o23mx9eINLGq8ud0/tAI5NcgyDsDgDeFG7JUmaZZ4orV9bl9z+dpLdwFOBjya5fDj/MUkuA6iqO4FzgMuB64GLq+raNuqVJA20dfXUh4APrTD/68CzRqYvAy6bYmmSpH3ocvfUXLJ5rD7bvn27t0Wfcd5GRJLUmKEh9YQf7NM42D0laaoMrtlmS0OS1JgtDWkFng1LK7OlIUlqzJaGZpotAmm6bGlIkhozNCRJjRkakqTGHNOQtCLHi7QSWxqSpMYMDUlSY4aGJKkxQ0OS1JihIUlqzKunpB7wSiiNiy0NSVJjhoYkqTFDQ5LUmKEhSWrM0JAkNWZoSJIaMzQkSY0ZGpKkxgwNSVJjqaq2axirJLcCX13Dfz0M+NaYy5kFbne/uN39ciDb/diq2ri/heYuNNYqyc6qWmy7jmlzu/vF7e6XSWy33VOSpMYMDUlSY4bGvc5vu4CWuN394nb3y9i32zENSVJjtjQkSY31PjSSnJrkhiS7kpzbdj2TkuToJJ9Kcl2Sa5O8cjj/kUk+nuRLw38PbbvWSUiyIckXkvzrcPqYJFcO9/v7kzyg7RrHLckhSS5J8l9Jrk/y1D7s7ySvGr7Hr0nyviQPmtf9neSCJLckuWZk3or7OANvG74GX0xy/FrW2evQSLIBOA84DTgOeGGS49qtamLuBP60qo4DTgJePtzWc4ErqupY4Irh9Dx6JXD9yPSbgLdU1eOB7wBntVLVZL0V+Leq+kXgSQy2f673d5IjgVcAi1X1y8AG4Azmd3+/Gzh1r3mr7ePTgGOHP1uBt69lhb0ODeAEYFdV3VhVPwYuAra0XNNEVNXNVfX54eP/ZXAAOZLB9l44XOxC4DntVDg5SY4Cfgt453A6wNOBS4aLzN12J3kE8DTgXQBV9eOqup0e7G8GX2N9cJKDgAcDNzOn+7uqPgPcttfs1fbxFuA9NfA54JAkRxzoOvseGkcCN41M7x7Om2tJNgNPBq4EDq+qm4dPfQM4vKWyJunvgD8H7hpOPwq4varuHE7P434/BrgV+Kdht9w7kzyEOd/fVbUHeDPwNQZhcQdwFfO/v0etto/Hcrzre2j0TpKHAh8A/riqvjv6XA0upZury+mSPBu4paquaruWKTsIOB54e1U9Gfg/9uqKmtP9fSiDM+pjgMcAD+G+3Te9MYl93PfQ2AMcPTJ91HDeXEpyfwaB8d6q+uBw9jfvbqIO/72lrfom5GTg9CT/w6D78ekM+voPGXZfwHzu993A7qq6cjh9CYMQmff9fQrwlaq6tap+AnyQwXtg3vf3qNX28ViOd30PjR3AscMrKx7AYMBsW8s1TcSwH/9dwPVV9bcjT20Dzhw+PhO4dNq1TVJVvbqqjqqqzQz27yer6sXAp4DnDRebx+3+BnBTkl8YznoGcB1zvr8ZdEudlOTBw/f83ds91/t7L6vt423AHwyvojoJuGOkG6ux3n+4L8mzGPR5bwAuqKq/brmkiUjy68C/A1dzb9/+XzIY17gY2MTg7sAvqKq9B9bmQpIl4M+q6tlJHseg5fFI4AvA71XVj9qsb9ySLDAY/H8AcCPwUgYninO9v5O8HvhdBlcMfgE4m0Hf/dzt7yTvA5YY3M32m8DrgA+zwj4ehug/MOiu+z7w0qraecDr7HtoSJKa63v3lCTpABgakqTGDA1JUmOGhiSpMUNDktTYQftfRNJKkjyKwQ3hAB4N/JTBrTsAvl9Vv9ZKYdIEecmtNAZJ/gr4XlW9ue1apEmye0qagCTfG/67lOTTSS5NcmOSNyZ5cZL/THJ1kp8bLrcxyQeS7Bj+nNzuFkgrMzSkyXsS8DLgCcDvAz9fVScw+LT2Hw2XeSuD73t4CvDc4XNS5zimIU3ejrvv8ZPky8DHhvOvBn5z+PgU4LjBnR4AeHiSh1bV96ZaqbQfhoY0eaP3OLprZPou7v0bvB9wUlX9cJqFSQfK7impGz7GvV1Vd99sUOocQ0PqhlcAi0m+mOQ6BmMgUud4ya0kqTFbGpKkxgwNSVJjhoYkqTFDQ5LUmKEhSWrM0JAkNWZoSJIaMzQkSY39P1G6gXrJYVxZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGCVJREFUeJzt3X2sZHd93/HPJ2NiDOMyhnXAtb2sXVxRO8BgLiQhbjoUWmxHeEkpkQkVmBi2JHHaCqmKkVUn5Z+4TaRUaa2QLUEYKcIY87QkixyDmZCNa8drNH6uYVlDvBs3Xgw2HUENHr79Y86Y8d2Ze3/zcJ5m3i/pas+cc+bOd8/9nvme3+93HhwRAgAgxU+UHQAAoD4oGgCAZBQNAEAyigYAIBlFAwCQjKIBAEhG0QAAJKNoAACSUTQAAMlOKDuAZduxY0fs2rWr7DAAoFbuvPPOb0XEqdutt3JFY9euXTp48GDZYQBArdj+Zsp6dE8BAJJRNAAAySgaAIBkFA0AQDKKBgAgGUUDAJCMogEASEbRAAAko2jUSKfTUafTKTsMrDByDNuhaKwAdnQUgTyDRNFYKezUKAJ5tt4oGgCAZBQNAEAyigYAIBlFo2D0ByNv5BjyRNGoKHZ85I0cwzwoGgCAZBSNCuHID3kjx7AoikbF9Ho9dmrkihzDIigaFTDt6O/AgQNqtVrPmJe6w3NEiXHkGJaFolEhvV5P/X6/7DCw4vr9vnq9XtlhoKYoGiXrdDrswMhdr9cjz7AUFI0SUTBQhE6nQwsWS0PRAAAkO6HsANYVrQzkjUFq5IGWRsX1+312fuSOM6GQqtSiYfvDth+1fe+U5bb9h7YP2b7b9vlFx1h3fBmgCOTZ+ii7pfERSRdusfwiSedkP3sk/VEBMRUir9NruXAL4/I6a4qzsdZXqUUjIr4s6dtbrLJb0kdj6DZJLdunFRNdvvr9vgaDQdlhYIVx3Q/yUHZLYzunS3p47PWRbB4AoARVLxpJbO+xfdD2wWPHjpUdDgCsrKoXjaOSzhx7fUY27xkiYm9EbETExqmnnlpYcACwbqpeNPZJekd2FtXPSnoiIh4pO6g8DQaDqQOM3W5X7XZ74vxut5tzZFgVg8Fg6ljHtBwbLSPPUOrFfbY/JqkjaYftI5J+W9KzJCkiPihpv6SLJR2S9D1J7yonUgCAVHLRiIi3bbM8JP1GQeEAALZR9e4pAECFUDRWyFb90cCykGfrjaJRMY1GQ+12W91uV81ms+xwsIIajYaazSY5hrlQNAAAySgaFdJutznyQ+6azSbdS5gbRQMAkIyiAQBIxpP7Cja6ovaEE9j0yMcox1qtVrmBYCXxzZWT0TMtpt124YILLtj2eQTtdvu4dbiNA8ZtlWfj4xYHDhyY+P5pt6UBpqF7qiR5nevO/YEwkmcucK3G+qJo5KDT6SQ/1azdbrPzYS6z5BlnTGFZKBoVwtEb8kaOYVGMaeRg1sdsFtmdtN1YC+pjljwb3WWgCOTYaqNoVMz4jsZOhzyQY1gERWNFjb4MRkd9wLJRcNYTYxpL1ul0ZuqaAuZBnqEstDSWaJazWaqC/ud6qWPLkRxbLRSNErEToQjkGZaJ7ikAQDKKxoJarRb3+EGuOp2OWq1WLbumsHooGkvQ7/dXaofudDor9f9ZBf1+v3bjZVshx+qLMY0VR382ikCerQ9aGjPg6Ah5I8dQdRQNAEAyuqfmtN3RYLPZpMmOhaS0OMgzFI2WBgAgGUUDAJCM7ik8jQFY5I0cqz9aGgvq9/saDAZlh4EV1uv1yDFUBkVjDnW8MWGKXq+3kv+vOiLHUFWlFg3bF9p+0PYh21dOWH6Z7WO2e9nPu8uIU1rdnRjVQY6hDkob07DdkHStpH8h6YikO2zvi4j7N6368Yi4ovAA1xDPZ0DeyLH6K7Ol8RpJhyLicET8QNL1knaXGM/MpvU1N5tNtdvtEiLCKpr2RUueoQxlFo3TJT089vpINm+zt9i+2/aNts+c9Its77F90PbBY8eO5RHrTNrtNhdcIXfkGcpQ9YHwz0naFREvl3SzpOsmrRQReyNiIyI2Tj311EIDXDfcGwl5I8eqrcyicVTSeMvhjGze0yLisYh4Mnv5IUmvKig2iJ0XxSDP6qXMi/vukHSO7bM0LBaXSvqV8RVsnxYRj2QvL5H0QLEhbq/RaJQdQm44k6cams2mnnjiibLDyA15Vi+lFY2IeMr2FZJuktSQ9OGIuM/2ByQdjIh9kv6d7UskPSXp25IuKyveVUWfOPJGjq2WUm8jEhH7Je3fNO/qsen3S3p/0XEBACbj3lNrhqM+5I0cW21VP3sKJVu155+jmsiz+qBoLFG32+ViK+Sq2+1yJI9SUTQAAMkoGgCAZAyELxldBygCeYay0NIAACSjaAAAklE0kIyHBCFv5Fj1MaaxgHa7rV6vtzIPlmk2m5J+3F/earVKjAYj7XZbBw4cKDuMpdicY6gfWhpL0Gw22QmQKx64hKqgaMxp9ACcdrvNzoxcjD9kiQcuoSooGgCAZIxpbDK6/82kozqO9rAM5BjqjKKBp9HNhryRY/VH9xQAIBktDTyNbhHkjRyrP4rGgtgJkDdyDFVC9xQAIBlFAwCQjKKBmfT7fe4NNIdOp8PjTBORY/MpKscoGgCAZAyEJ1rHwcjRDRlRjHXMMYk8qxtaGpjJpBvn0fWCZZp2c0byrBooGphqdEPGVOzUmAd5Vi9J3VO2z42I+zfN60REN5eoUAvsuCgCeVYtqS2NG2z/lodOsv3fJf1unoEBAKondSD8ZyT9F0m3SjpZ0p9K+vm8gkJ1rOvgLIo1yjNaFdWX2tL4oaTvSzpJ0rMlPRQRP8otqpLwfGLkjRxD3aUWjTs0LBqvlvRPJb3N9icW/XDbF9p+0PYh21dOWH6i7Y9ny2+3vWvRzwQAzC+1aFweEVdHxA8j4pGI2C1p3yIfbLsh6VpJF0k6V8NCdO7mz5X0nYh4iaQ/0LCLDABQktSi8ajtneM/kv5ywc9+jaRDEXE4In4g6XpJuzets1vSddn0jZJeb9sLfi4AYE6pA+F/LikkWcMxjbMkPSjpvAU++3RJD4+9PqLhgPvEdSLiKdtPSHqBpG8t8LlArrZ6nCuwDGXmWFLRiIiXjb+2fb6kX88lojnY3iNpjyTt3Llzod/FM5pn0+v11O/3p17Fi+ORY7M7cOCAJOmCCy4oORLMde+piPiK7c2tglkdlXTm2OszsnmT1jli+wRJz5P02IR49kraK0kbGxuxYFzAQjg7CnkrM8dSrwh/39jLn5B0vqS/W/Cz75B0ju2zNCwOl0r6lU3r7JP0Tkn/S9K/lnRLRFAUAKAkqQPhJ4/9nKjhGMfmQeuZRMRTkq6QdJOkByTdEBH32f6A7Uuy1f5E0gtsH5L0PknHnZaL8g0GA/X7/bLDqJR+v8+Faks0GAw0GAzKDqNSysqx1DGN/5zHh0fEfkn7N827emz6/0l6ax6fjfnQH4+8kWPVtmXRsP05Dc+amigiLpm2DACwerZrafx+IVGgFjj6Q97Iserbrmg8FBF/W0gkWBmtVkuS9Pjjj5ccCVZVr9dTq9WiK6sE2w2Ef2Y0YfuTOceCGhsMBpxqityRY+XbrmiM37Lj7DwDQf1s9cQ1zh7CsnS7XTUajYnLer0eeVaw7YpGTJkGAKyh7cY0XmH7uxq2OE7KppW9joj4B7lGB9QAR7rI2/i9pjqdTqnXRW1ZNCJicpsQwHHob0cRer2eBoPB1C67vKVeEQ4gUb/f5+pl5K6sHKNoAACSUTQAAMkoGlhYo9FQo9HgQivkbpRjPLulPHM9TwPAbEbFlDOtkJeiDthoaQAAklE0AADJ6J4aQ3888kaOoe5oaWAh3W5XzWaz7DAqpdlslnbh1apimx6Pi/sAAJVH0QAAJKNoYGHtdpsuKuSKHKsOigawZHzBIW+NRqO0HKNoAACSUTSAGXQ6Ha7qRq6qnmMUDSwF951C3sixaqBoYOm269PvdDpqtVqVPppCtaXctLDVaqnVahUU0fqgaCBXVW9qo/56vR45ViBuI4KF0WWAvJFj1UFLA0vH8w5QBPKsHBQN5KLb7a7d0WG73X76S4xB22KsY+EYjRmWlWMUDWBJ1rFQonhlF8pSiobt59u+2fbXsn9PmbLewHYv+9lXdJwAA/nIW91yrKyB8CslfTEirrF9Zfb6tyas9/2IWK+254rrdDrq9Xplh5GsTrFiqE5fwCOjM8Dq0FItq3tqt6TrsunrJL25pDgAADMoq6Xxwoh4JJv+P5JeOGW9Z9s+KOkpSddExGcmrWR7j6Q9krRz585lx4oc9Pv9Wh3F1/Hodd31+/2yQ5hJXa43ya1o2P6CpBdNWHTV+IuICNsx5de8OCKO2j5b0i2274mIr29eKSL2StorSRsbG9N+FwBgQbkVjYh4w7Rltv/e9mkR8Yjt0yQ9OuV3HM3+PWy7K+mVko4rGgCAYpTVPbVP0jslXZP9+9nNK2RnVH0vIp60vUPSz0v6r4VGCSSow+Al6q1KOVZW0bhG0g22L5f0TUm/LEm2NyS9NyLeLemfSPpj2z/ScMD+moi4v6R4gWRV2sGxmsrMsVKKRkQ8Jun1E+YflPTubPpWSS8rODQAwBa4Ihwrq24XTaGe1i3PKBoAgGQUDQBAMooGkFm3bgaUo+55RtFAqeq+A6EeyLPloWigctjBkTdybH487hWlGd1/avRsgGk78Wh+6rnpfBlgZDAYHJcPozstjz+TYtYcG3/PuqGlgbXEkSbytqo5RtEAVL/nfKCeViHPKBoo3OgZx0Cems2mGo1G2WGsHIoGMEWn06ndMxlQP71er1Z5xkA4CjUaaGy1WuUGkqDX62kwGJQdBmZUpxyr44EJLQ2UYp4uqlUdWEQ+yLF8UDRQim63+4xTHlP0er3aDyKiOPPkmFSfx66Whe4pYAY8KwN5q3qOUTSQqyrtACkXcI2OMkfrNBqNuY5WUZy65djootZxzWazNnlG9xRqjT5oFIE8+zGKBjBBHc9qQf3U8Qw9uqewckZN/2U296vUBYJqWHae1SXHaGkAW2g2m7XZmVFPdRs3o6WBUrXbbb6UkTtybHloaQAAktHSQOVwVIi80cKdHy0NrJ1VuD01qm9V72BASwOlST3SW8UdD8WYpTVBnqWhpQEASEZLA5VRZB8zfdrrq6i/e51Oo50FRQMra9EvhzrdDwjlWSTP2u127brF6J4CACQrpWjYfqvt+2z/yPbGFutdaPtB24dsX1lkjACA45XV0rhX0r+S9OVpK9huSLpW0kWSzpX0NtvnFhMeqqjdbid3F/X7fe5Kirmkjnf1+/21vKllKWMaEfGAJNnearXXSDoUEYezda+XtFvS/bkHiFqb1k+83Q3mGL/ALCY9SjZlfKLuY2VVHtM4XdLDY6+PZPOwZmZpYQDzmPfRsOsot5aG7S9IetGERVdFxGeX/Fl7JO2RpJ07dy7zV6PiOG0WRSDPfiy3ohERb1jwVxyVdObY6zOyeZM+a6+kvZK0sbERC34uVtSsR5Jcy4FZzdNaqVueVbl76g5J59g+y/ZPSrpU0r6SY8KK6na7tdpxUU+r0A1W1im3v2T7iKSfk/Tntm/K5v9D2/slKSKeknSFpJskPSDphoi4r4x4AQBDZZ099WlJn54w/+8kXTz2er+k/QWGBgDYArcRQW0U0X00+gyu8VhPRYwv1L0btMpjGgCAiqGlAUxQ96NB1EMd84yigbVTxx0V9bOqeUb3FAAgGS0NrJxut8tANnJX9+st5kVLAwCQjJYGKm9V+4ZRLeRZGloaAIBkFA0AQDKKBgAgGWMaWEn0TyNv65pjtDQAAMkoGgCAZBQNAEAyigYAIBlFAwCQjLOnsDbW9WwXFGcdcoyWBgAgGUUDAJCMogEASEbRAAAko2gAAJJRNAAAySgaAIBkFA0AQDKKBgAgmSOi7BiWyvYxSd+c8+07JH1rieEsC3HNrqqxEdfsqhrbqsX14og4dbuVVq5oLML2wYjYKDuOzYhrdlWNjbhmV9XY1jUuuqcAAMkoGgCAZBSNZ9pbdgBTENfsqhobcc2uqrGtZVyMaQAAktHSAAAkW4uiYftC2w/aPmT7ygnLT7T98Wz57bZ3jS17fzb/QdtvLDiu99m+3/bdtr9o+8Vjywa2e9nPvmXGlRjbZbaPjcXw7rFl77T9teznnQXH9QdjMX3V9uNjy3LbZrY/bPtR2/dOWW7bf5jFfbft88eW5bm9tovr7Vk899i+1fYrxpZ9I5vfs31wmXElxtax/cTY3+zqsWVb5kHOcf3HsZjuzfLq+dmy3LaZ7TNtfyn7TrjP9r+fsE7+eRYRK/0jqSHp65LOlvSTku6SdO6mdX5d0gez6UslfTybPjdb/0RJZ2W/p1FgXK+T9Jxs+tdGcWWv+yVvs8sk/Y8J732+pMPZv6dk06cUFdem9X9T0ocL2ma/IOl8SfdOWX6xpM9LsqSflXR73tsrMa7Xjj5P0kWjuLLX35C0o8Rt1pH0Z4vmwbLj2rTumyTdUsQ2k3SapPOz6ZMlfXXCfpl7nq1DS+M1kg5FxOGI+IGk6yXt3rTObknXZdM3Snq9bWfzr4+IJyPiIUmHst9XSFwR8aWI+F728jZJZyzpsxeObQtvlHRzRHw7Ir4j6WZJF5YU19skfWxJn72liPiypG9vscpuSR+NodsktWyfpny317ZxRcSt2edKxeZYyjabZpH8XHZcRebYIxHxlWz6/0p6QNLpm1bLPc/WoWicLunhsddHdPyGfnqdiHhK0hOSXpD43jzjGne5hkcQI8+2fdD2bbbfvKSYZo3tLVkT+EbbZ8743jzjUtaVd5akW8Zm57nNtjMt9jy316w251hI+gvbd9reU1JMP2f7Ltuft31eNq8S28z2czT84v3k2OxCtpmHXeivlHT7pkW559kJ87wJxbL9byRtSPpnY7NfHBFHbZ8t6Rbb90TE1wsM63OSPhYRT9r+txq21P55gZ+/nUsl3RgRg7F5ZW+zyrL9Og2LxgVjsy/IttdPSbrZ9v/OjsKL8hUN/2Z92xdL+oykcwr8/O28SdJfR8R4qyT3bWa7qWGh+g8R8d1l/u4U69DSOCrpzLHXZ2TzJq5j+wRJz5P0WOJ784xLtt8g6SpJl0TEk6P5EXE0+/ewpK6GRx3Lsm1sEfHYWDwfkvSq1PfmGdeYS7Wp2yDnbbadabHnub2S2H65hn/D3RHx2Gj+2PZ6VNKntbyu2SQR8d2I6GfT+yU9y/YOVWCbZbbKsVy2me1naVgw/jQiPjVhlfzzLI8Bmyr9aNiaOqxhV8Vo0Oy8Tev8hp45EH5DNn2enjkQfljLGwhPieuVGg74nbNp/imSTsymd0j6mpY7EJgS22lj078k6bb48YDbQ1mMp2TTzy8qrmy9l2o4IOmitln2e3dp+qDuL+qZA5R/k/f2Soxrp4Zjda/dNP+5kk4em75V0oXLjCshtheN/oYafvn+bbb9kvIgr7iy5c/TcNzjuUVts+z//lFJ/22LdXLPs6UmQFV/NDyj4KsafgFflc37gIZH75L0bEmfyHaev5F09th7r8re96CkiwqO6wuS/l5SL/vZl81/raR7sp3lHkmXl7DNflfSfVkMX5L00rH3/mq2LQ9JeleRcWWvf0fSNZvel+s20/CI8xFJP9Swv/hySe+V9N5suSVdm8V9j6SNgrbXdnF9SNJ3xnLsYDb/7Gxb3ZX9na/KIce2i+2KsRy7TWOFbVIeFBVXts5lGp4kM/6+XLeZhl2HIenusb/XxUXnGVeEAwCSrcOYBgBgSSgaAIBkFA0AQDKKBgAgGUUDAJCMogHMYeyOuffa/oTt59jeNe3OqMCqoGgA8/l+RLQj4qcl/UDDc+WBlUfRABb3V5Jekk03bP/P7HkHf2H7JEmy/R7bd2Q33/tkdrM72X5r1lq5y/aXs3kN27+XrX93dm8voBIoGsACsnuVXaTh1bfS8IZ610bEeZIel/SWbP6nIuLVEfEKDW9pfXk2/2pJb8zmX5LNu1zSExHxakmvlvQe22fl/78BtkfRAOZzku2epIMa3hPpT7L5D0VEL5u+U8N7GEnST9v+K9v3SHq7hvc1k6S/lvQR2+/R8OFCkvQvJb0j+/23a3ib/ird3RVrjFujA/P5fkS0x2cMn9ulJ8dmDSSdlE1/RNKbI+Iu25dp+FQ6RcR7bf+Mhjeau9P2qzS8f9BvRsRNef4HgHnQ0gCKcbKkR7JbW799NNP2P4qI2yPiaknHNLx99U2Sfi1bV7b/se3nlhE0sBktDaAY/0nDrqZj2b8nZ/N/z/Y5GrYuvqjhHVLv1rBb6yvZY4ePSSr6SYPARNzlFgCQjO4pAEAyigYAIBlFAwCQjKIBAEhG0QAAJKNoAACSUTQAAMkoGgCAZP8f4zx/GRKfE8wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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": [
    "As you can see, in this form we can clearly see this periodic modulation in the lightcurve. This is the exact same data as above, but put in a more useful form. However, for the synthetic data, we know what the period is because we generated the lightcurve with a period=1. In reality though, one usually doesn't know the period beforehand, so now we learn to find it using an algorithm known as lomb scargle.\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 a period=1. Try changing the value of the period, and watch the periodogram change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.9987063533134721"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8W/W5P/DPYzuOs5cdskkCYSRAGBkUOugFWqAt0Mu+XXRxKaWUDu4PCgUKHdANlL33Hk0hG7Ig09mJsxxneO9tybKk5/fHGT6SJVs+tuzk+PN+vfKKJB8dfc98vvuIqoKIiAgAUno7AUREdORgUCAiIhuDAhER2RgUiIjIxqBAREQ2BgUiIrIxKBARkY1BgYiIbAwKRERkS+vtBHRWZmamTp48ubeTQUR0VNm4cWOFqmZ1tNxRFxQmT56M7Ozs3k4GEdFRRUQOJbIcq4+IiMiWtKAgIs+JSJmI7IjzdxGRh0UkV0S2iciZyUoLERElJpklhRcAXNTO3y8GMM38dwOAx5OYFiIiSkDSgoKqrgRQ1c4ilwF4SQ1rAQwXkbHJSg8REXWsN9sUxgPId7wvMD8jIqJeclQ0NIvIDSKSLSLZ5eXlvZ0cIiLP6s2gUAhgouP9BPOzNlT1KVWdpaqzsrI67GZLREQu9WZQmAfgu2YvpLMB1KpqcS+mx1O25tdge0FtbyeDiI4ySRu8JiKvAzgPQKaIFAC4B0A/AFDVJwDMB3AJgFwATQC+n6y09EWXPfoZAODgA1/r5ZQQ0dEkaUFBVa/r4O8K4KfJ+n0iIuq8o6KhmYiIegaDAhER2RgUiIjIxqBAREQ2BgUiIrIxKBARkY1BgYiIbAwKRERkY1AgIiIbgwIREdkYFIiIyMagQERENgYFIiKyMSgQEZGNQcGDjFnJiYg6j0HBgxgTiMgtBgUPYkwgIrcYFDyI1UdE5BaDAhER2RgUPIjlBCJyi0HBg1h7RERuMSh4kLKsQEQuMSh4EEsKROQWgwIREdkYFDyIJQUicotBwYPYpkBEbjEoeBBLCkTkFoOCBzEmEJFbDApERGRjUPAgzn1ERG4xKHgQQwIRuZXUoCAiF4nIHhHJFZHbY/x9kogsE5HNIrJNRC5JZnr6ChYUiMitpAUFEUkF8CiAiwFMB3CdiEyPWuwuAG+p6hkArgXwWLLS06cwKBCRS8ksKcwBkKuqeaoaAPAGgMuillEAQ83XwwAUJTE9fQbHKRCRW2lJXPd4APmO9wUA5kYtcy+AxSLyMwCDAFyQxPT0Gc7qI1WFiPReYojoqNLbDc3XAXhBVScAuATAyyLSJk0icoOIZItIdnl5eY8n8mjjLCewfYGIOiOZQaEQwETH+wnmZ04/BPAWAKjqGgAZADKjV6SqT6nqLFWdlZWVlaTkeoezSypjAhF1RjKDwgYA00Rkioikw2hInhe1zGEA5wOAiJwMIyiwKNBFkSUFhgUiSlzSgoKqBgHcDGARgF0wehntFJH7RORSc7FfAfixiGwF8DqA65V3sS6LaFPovWQQ0VEomQ3NUNX5AOZHfXa343UOgHOTmYa+jiGWiDqjtxuaKQmcXVLZPZWIOoNBwYsiuqT2XjKI6OjDoOBBjANE5BaDggcpSwpE5BKDggexTYGI3GJQ8CCWFIjILQYFD9I4r4mIOsKg4EER01ywqEBEncCg4EEc0UxEbjEoeBwLCkTUGQwKXsegQESdwKDgQZHVR4wKRJQ4BgUPihinwJhARJ3AoOBBbGgmIrcYFDyID9khIrcYFDyIj+MkIrcYFDwosqTQa8kgoqMQg4IHsfcREbnFoOBJbGkmIncYFDyIvY+IyC0GBQ9imwIRucWg4HFsUyCizmBQ8CA+ZIeI3GJQ8KDIx3ESESWOQcGDIksKDAtElDgGBQ9i9RERucWg4EFsXCYitxgUPIglBSJyi0HB41hqIKLOYFDwIJYUiMgtBgUPYpdUInIrqUFBRC4SkT0ikisit8dZ5moRyRGRnSLyWjLT01ewSyoRuZWWrBWLSCqARwFcCKAAwAYRmaeqOY5lpgG4A8C5qlotIqOTlZ6+iiGBiDojmSWFOQByVTVPVQMA3gBwWdQyPwbwqKpWA4CqliUxPX0GJ8QjIreSGRTGA8h3vC8wP3M6AcAJIvKZiKwVkYuSmJ4+I7LKiFGBiBKXtOqjTvz+NADnAZgAYKWInKqqNc6FROQGADcAwKRJk3o6jUcdlhSIyK1klhQKAUx0vJ9gfuZUAGCeqrao6gEAe2EEiQiq+pSqzlLVWVlZWUlLsFfwITtE5FYyg8IGANNEZIqIpAO4FsC8qGU+gFFKgIhkwqhOyktimvoIR5dURgUi6oSkBQVVDQK4GcAiALsAvKWqO0XkPhG51FxsEYBKEckBsAzAbapamaw09RWRJQVGBSJKXFLbFFR1PoD5UZ/d7XitAH5p/qNuwjYFInKLI5o9iNNcEJFbDAoe5OySyuojIuoMBgUPYvUREbnFoEBERLYOg4KIpIjI1T2RGOoebFMgIrc6DAqqGgbwfz2QFuomkVNnMyoQUeISrT5aKiK/FpGJIjLS+pfUlJF7LCkQkUuJjlO4xvz/p47PFMDU7k0OdQdOh0dEbiUUFFR1SrITQt2HD9khIrcSqj4SkYEicpeIPGW+nyYiX09u0sgtPo6TiNxKtE3heQABAOeY7wsB/D4pKaIuY+8jInIr0aBwnKr+GUALAKhqEwBJWqqoS7Sdd0RE7Uk0KAREZADMO4yIHAegOWmpoi6JmOaCMYGIOiHR3kf3AlgIYKKIvArgXADXJylN1EXsfUREbiXa+2ixiGwEcDaMaqOfq2pFUlNG3YIlBSLqjISCgoi8AmAFgFWquju5SaIuY5dUInIp0TaFZwGMBfCIiOSJyLsi8vMkpou6gF1SicitRKuPlonISgCzAXwZwI0AZgB4KIlpI5fYJZWI3Eq0+uhjAIMArAGwCsBsVS1LZsLIPT6jmYjcSrT6aBuMwWunADgNwClmF1U6AmncN0RE7Uu0+ugXACAiQ2B0RX0ewBgA/ZOWMnIt8nGcRESJS7T66GYAXwBwFoCDAJ6DUY1ERyA+jpOI3Ep08FoGgL8D2KiqwSSmh7oB2xSIyK1Eq4/+KiIzAdwoIoAxXmFrUlNGXcBpLojInUSnzr4FwKsARpv/XhGRnyUzYdQ9GBOIqDMSrT76EYC5qtoIACLyIIzuqY8kK2HkHh+yQ0RuJdolVQCEHO9D4NTZRyxOiEdEbiVaUngewDoRed98fzmMqS/oCKSMCkTkUqINzX8XkeUAPm9+9H1V3Zy0VFGXRM59xKhARIlrNyiISAaMeY6OB7AdwGPsknrk49xHRORWR20KLwKYBSMgXAzgr0lPEXUZB68RkVsdVR9NV9VTAUBEngWwPvlJoq7iNBdE5FZHJYUW64WbaiMRuUhE9ohIrojc3s5yV4iIisiszv4GtY9dUomoMzoqKcwUkTrztQAYYL4XAKqqQ+N9UURSATwK4EIABQA2iMg8Vc2JWm4IgJ8DWOdyGyhK5DQXRESJa7ekoKqpqjrU/DdEVdMcr+MGBNMcALmqmqeqAQBvALgsxnL3A3gQgN/VFlAbymkuiMilRAevuTEeQL7jfYH5mU1EzgQwUVU/am9FInKDiGSLSHZ5eXn3p9TTGBWIKHHJDArtEpEUGDOv/qqjZVX1KVWdpaqzsrKykp+4oxy7pBKRW8kMCoUAJjreTzA/swyB8SS35SJyEMDZAOaxsbnr2KZARG4lMyhsADBNRKaISDqAawHMs/6oqrWqmqmqk1V1MoC1AC5V1ewkpqlP4DgFInIraUHB7MJ6M4BFAHYBeEtVd4rIfSJyabJ+l6LHKTAqEFHiEp0QzxVVnQ9gftRnd8dZ9rxkpqUvYUmBiNzqtYZmSiK2KRCRSwwKHhQ5ToFhgYgSx6DgQYwDROQWg4IHsU2BiNxiUPCgyHEKjApElDgGBY9jSYGIOoNBwYM4IR4RucWg4EEMBETkFoOCB2mc10REHWFQ8CLlOAUicodBwYNYUiAitxgUPEgZFYjIJQYFD+IsqUTkFoOCB3FEMxG5xaDgQXzyGhG5xaDgQSwpEJFbDAoexzYFIuoMBgUPimhoZkwgok5gUPA4xgQi6gwGBQ9SNioQkUsMCh4UMUtqL6aDiI4+DAoeFNEllVGBiDqBQcGDImuPGBWIKHEMCh7EwWtE5BaDggfxyWtE5BaDggexpEBEbjEoeBzbFIioMxgUiIjIxqDgQZzmgojcYlDwoMg2BUYFIkocg4IHcZYLInIrqUFBRC4SkT0ikisit8f4+y9FJEdEtonIxyJybDLT01eEldNcEJE7SQsKIpIK4FEAFwOYDuA6EZketdhmALNU9TQA7wD4c7LS05dwmgsiciuZJYU5AHJVNU9VAwDeAHCZcwFVXaaqTebbtQAmJDE9fUZEQzPLCkTUCckMCuMB5DveF5ifxfNDAAuSmJ4+g20KRORWWm8nAABE5NsAZgH4Upy/3wDgBgCYNGlSD6bs6BTRpsCoQESdkMySQiGAiY73E8zPIojIBQDuBHCpqjbHWpGqPqWqs1R1VlZWVlIS6yWqgEjrayKiRCUzKGwAME1EpohIOoBrAcxzLiAiZwB4EkZAKEtiWvqUsAKpZlQIMygQUSckLSioahDAzQAWAdgF4C1V3Ski94nIpeZifwEwGMDbIrJFRObFWR11gkKRkmIFBUYFIkpcUtsUVHU+gPlRn93teH1BMn+/r1KzpCDCNgUi6hyOaPYgVYUIkCLC6iMi6hQGBQ8KqxEQUoTjFIiocxgUPEgVEAAClhSIqHMYFDwobFYfiXi3oTkYCqO6MdDbySDyHAYFjxIRpIh4dpzCH+bvwhn3L0FTINjbSSHyFAYFDwqrIkWAFAHCHq0/WrC9BABQxdICUbdiUPAgo/pIPN37aGB6KgCgpqmll1NC5C0MCi6s2V+JS//1Kfwtod5OSkyqRinBy20KGf2MoMCSAlH3YlBw4R9L92JbQS3WHajq7aTEZJQOBCkp4tnBa6nmiO2mwJEZmImOVgwKLgwf0A8A0NR8pDZyWm0K4vlRCs1BBgXqfXX+FqzZX9nbyegWDApdEAiFezsJMYXDrYPXjrTqo8IaH65/fj0qGmJOiNtpR2oVHvUtN7+2Gdc9vRa1vqO/jYtBwQVrWurmYPygsPlwNUrr/D2UokgKNdOYvIZmf0vI1QVw77ydWL6nHGvzuidX5W/p+cAcDis+2V3q2ao5avXhtiK8uu5Qh8vlFNUBACq7KbPTmxgUXBAYUSFeUAiHFd98bDWufnJNTyar9fed01wk6cb134+txszfLe7093xmG0B3JcvXCyWFN7Pz8YMXsvH+5jaPByGXVu4tR1GNr7eT0cbNr23Gne/v6HC5fqnGPaGsnkGhT2uOc0MqNE/uQ5VNMf/e3YprfZh8+0dYttt4JIV1w00RQThJGemc4jpX3wuZRZc6f9eK2dacTr1RfVRuXvi5ZQ09/tteVNMUwHefW4+bXt3U20lxzbrm6v1Hajtj4hgUuiBem4LVdz49tWd27+pcoyrmPTPnqqpISUGPtCm0dLJdxQoKXb14AmYprTeqjwb1N2acbzhiOxr0jHBYMfn2j/C3xXu6tJ68ikYARubGEgprt1UxdoeOznOrSrk3Sq7djUHBBSuX2hznhtTTUy9YdfsjBhq9ohRGFZf0wOC1zg4eazYvrsYu3lCtYNAbJYX0NOOyaQkdHW0Kt7+7DW9n53f7eq0S8SOf5HZpPf4Y3Yr/smgPrn1qLbYV1HRp3d2lo/M1xYwKsbalsxqag53ObHUnBgUXrJtBvDYFK7cQCIVxuAeqkKzSQIoIAsEwDlc1GV1SU5I/dbavkxdBbZMx2KyrN3Pr+70RFFrM495DBcEue2NDPm57Z1u3r7e72gCs68WZgbGqQjt7fnUnZ3tcR+NhrJKCvxu6SJ9yzyL87LXNXV6PW0fJaX1ksfrGx+sj7zyRv/iXZSip7f5eSOsPVOGRj/cBaD1hfYEQfvHWFmw8VI2mQAgtQcV7mwq7PZfovFg6exFUmyWLRIvZwVAY33l2HTYcjBwo6OsgKORXNeHcBz7BLpdtH+2xftvKHR7Jkjn3VXcNHLT2Z6xOEU29WB3jrB7u6Hy1BlN2NYhZ95SFO0u6tJ6uYFBwwarP7qikYNlRWNvtabj6yTX425K9CATDdtH2zex8fLStGIBR3K1sNBpE/7Z4b8x1uM1lO28G8arQYlFVu4HZF0jse8W1fqzaV9Em59Tc0n6bwtvZ+Sis8SWlh5Czg8GR3i21PontHt0WFAJtSwppZm+ejqpt1uVV4vHl+7slHU6qGnGD7+hmb2UQutqmEGvalkAwjPnbi3vsXGNQcMEOCnHbFCJPjGQ2PhXX+tAYow3DHwzb1VzWBeaUfbAKJ/12IZbklHa6wdR54nZmRHFzMGz30kg0IAXNO4WzRBIKq52Li7dvk3n5WL/50ppDvdbtuD0r9paj1iyR1SVxMFV3tZ357eqj1qOWlpJYULjmqbV4cOFu+5rsLje+shGn37fEft/RNWz9fjKCwpMr9uOmVzdhSU5pl9adKAYFU2dyzc12SSH2d6LX1d31oiFHlqqqMYCm5rbrdy7TL0bl98IdRvH0xy9l48z7l7T5e3ucXTHbG8AXzblfEr14rJtCTVMLfv32VrSEwhH7Pd5xsxrA67vY9TUWZ9o3HKzG5Ns/wvztxd3+O/G0l2NsbA7ie8+tx/dfWA8guW0uvm4qMdltCo5zNsUOCsbfQmG1j3thjQ8F1ZFtdYXdPMZh0c7IG3BHpSIrk9KZhuaqxgD+9cm+iEblWPeKgmpfxP/JxqAA4EBFI0767UJ8sLkwoRxHR9VH0SdQrJx8VxyqbLRf+wKhDtcfjDFYwfmdzuaynCO1O3PTaepEcTzWd97ZWIBDlY0R3/XHSXuV2aDdGCNgdlWsqq+udstM1PztxZhxzyKU1cdup7K6+m46bPTa6UzQXpJTihc+O5Dw8hHViF3IqVv701l9ZJVtrRLiD17YgBPvWggAOPeBT/D5B5dFrCOZJSIjje1fY25KCre+uQV/XbzX7mFVVu+PWRVmBciuju1JFIMCYDdi3vrmFpxy76IOl7cugHg30+ig0N0zeeaVtwaFpkCow/XHqne3GsbccFY3NQfDuOuD7Zj1+6UIdtCNzuempBB1MTY2hyICQbycmfWoTjddX5ftKcMfPsqJ+/dYgbAnBi2pKv6yaA+aAqG4A+camiNvHJ2Zn+vHL2Xj3v/E3m5fIIR5W4vi9sjpyvZb54KzBGgFCOv4rthbDqB14GC07s54RYtbTamKxTtL7GvCl2Ab277SelSY22KVar/19Dp8bPa6crLuMz317JA+GxReXnsIN726EUDkjSMQDHdYFO7t6qNiR069qSXU4Y0v3shrp84U/5257+ZgCK+sPYyKhuYO50Ky9oNI4iWM6KqxxkAw4rvxej9ZdbNuBph9//kNeHrVgbj7JNYNIhijl4+q4qU1B9tUdbhR3RjAlDvm44A50CveDSL65uzMuCTa9z3WPnt8xX7c8vpmLHbUaztzz10ZyGcdz5aQIhAMo6Kh2b75R5cEZ/9hacx1xKpC7Wp6ItYf5xpesbccN7y80a6uTeRaX7G3HBf+Y6U9K0B1Uwu+/cw67IsT6K22m556dkifDQq//WAH5m8vQSisaIi6kOp88U/w1bkV9k14bV4V1sUYdRndANeVXMy7GwtwzZNrsLOotQdTqaOLqy8Q7LCkUOcP4vX1hyM+i64C6czIYOf2OL9X10Fu0brYRg5Mb7ek8P/e2dbaiypq3zU1h+wLb/jAfnGDi3UBxds3W/Nr8NV/rGy3Z1ij47t1/hb84IUNWLWvPOZvxgogOwrrcPe/d+Kef+/ExkPVcX8HMOrMNx+Ov8zm/Mi/lZkZg6IaH7728Cos22PkMKNvzs6gkGjpLFZwt27S0aVUS/Q11Bm+iBJHC2b9fqndRpBo5qE7SwpldW1LI/Fu9vlVkQE/kfRGT5pX0xTAp7kVbZY7aGYArP1c3cSg0K1Kav1YEKNbV2Vjc5tue1Vxdv7+8gb8zzPrIi6ua55a22a56CKk25JCvb8Fv3p7K9YdqMIDC3bbnxfX+tHfHFXblECbAgDc8d72iPddCVwNzcHWmWId+6Kjel1rv40YlN7aDTGsEdVOgWAYb2bn46evGfPgRD+zojEQtEtoIwamxwxmobCisoPqoyU5pdhTWo+lu+L36Ljx5Y14bLkxWnfjoWp8srsM/1y6L+Gb60Gz7efj3WW44vHVWN/OQ5l+/1EOvvnY6rjjKqJLAFau8j9bi7CzqA7PrjrQZrnvPLsO/1rWOto4VlXbs58ewGX/+jTis8bmIO54bxs+3dd6o6oyuzc7p6KIqD5qbkEwZJSy2yt1Pr58P0767QI7Zx0OK4oc64xuMG7v2nHegLuriras3o/SGO018dIRfS4kcm6EokqV8RqQz/vrcny6r8K+VhkUutmr6w7hJ69uwgurD0Z8XlEfaNNDJd7Oj1ef6ew1UV7f3OZG1Bh1QkX3gY7HqioAgN0l9fbr0jo/pmYNBmCchIkWnZ0Xa3SaElmHquLJFfuxr7Qeowb1BxDZwNhRQ5i1zc6SwkMf78Pxdy6wu1BGF5Gj07k6txJXPG50Ax0+sF/Mi7Cysdm+8OIFO+t70Tlc543m09wK/Hmh0YB8wMwhp6ZIwkE+elu2t1MqsR7QcrgqdlVTdAlg1b4KfLyrFJvM0oX1fArn9qzaVxERiHwtITz76QF7mmcAuP/DHGwtqLX3PwBsL6jF6+vz8ROzehUADlYY6XLewJwZi4JqH876/VJMuWM+rnmybUbJ8uDC3fC3hLG31Difn1yZh1WO4FMYdYOM15EAQMSg0K5OmwIYOfY5f/gYv3pra5u/xRtEFx2snedGczAUM0BGf2eP49qOlltWb1fXVjeyTaFb3XrBCZg7ZSSe++xAxE281tfS5iDVxqmvrY5Tp2d9//nPDmD2H5bik6jGouieC4t2luDkuxd2ONrWOukvnH4Myuub7ROupM6PY0cOhEhivY/s9Du2KzoH/rM3Nnc4F/zhqib8acFubDhYjVGD0gFEBYUY1W6HK5vsoNtaUuiHen8Qv3prq52T3WPeJJwP31HVNiWaNx2js4cP6IdAMNwm52UV/yeMGGBfUG+sP4yz//ixvaw1xXGx48ay/oAxdiOaqkYE6ESrNKKDQqx5fHaX1OGmVzfaOd14U0dEn6OHq5rwwxez7a6T1n5rb7DavtIG3P9hjt2W5nSoypkBMc5LXyCEen8Ljv/NfPv4OKtLmgIhDDfn21qaU2pXO60/WIVfv7015hQv1hTTlQ3GvlmcEzlyN7qk4G8JxR2V/c7Ggoi0dMWH24pwyUOrAMQOzNEZgYU7irHpcHWb42K1cZXV+3HiXQvx+vp81Da14P3NBXEng7QCZCyFNT77umFJoZulpggunH4M8qt8EQe91hdoc5Di7fx41UrW8lbf/2hLd5XhG498ih+9uAHPrMrDkhwjaHwWox7RqdS8cc2dMhIAkG82WJbU+jFmWAYG9ktFvT+YcHtARM4q6iTfml+Dh81pM+Jx7rehA9KQliIRF+NPX9uEO95rnWOnvL4ZX/zLMvz2gx0oqvHhyRV5AICRZinj3U2tF0pRjQ/ztxfjz4tau3be/NrmdvtmjxmWAaBtHbjVZXZq1mA0BoJQVdz1wQ6U1Pnxoxc34MXVB+1lnFUX0VNpWOqbg3ZVUFVjIGbpxAqOK/aW2zfO6PNoZ1HbTMADC3Zj/vYSe9/GCwod1dlXNgYQDIXbXW6H2S4Vq8HSOc37x7uM8zMYVmw4WGU3og9MT0VBtc/O/db5gxg/fACA1i6wd33tZADGDfs370dWWQKt84atyi1HMBTGkIx+EX+PPt7+lhAa4mR6nFVjXW1TuPm1zShqZzoaZ1BQVdz4yib892OrI+4dKdJaRbc2zziXnlmVhzc2HMYv3tyKt7Pz8fGuUvs4WCrbaUA+VNlkl4KaAqEemeurzwQFADhutFHlstWRYzNKCtHVR21LCsv3lOGhpbFvmlVNAahqu88Y2F5Yi6W7yvD7j3bZdeIdPZCjtNaP1BTBGZNGADBOkHp/CxqagxgzLAMD0lPbPaHapNOxbKzRqFb3xcIaHx7+eB9CYcXSnFJc8PcV2FFYixV7yu1lB/VPQ/+0FLuu2fL6+tac/L4yIwf0wZYiXPXEGnv/jBwUeSOwfvOmVzdh5d7W3/hoezH+vaUo7vZYN6SqxuaIdolSs6QwNXMQVI0SijUNwbI95Xhl7SG7oba4pvVGEG8uo/L65sigECNX2hQIobjWh+89tx5XPL7aXtYpr7zB/m5pnR8/ejG7zbxYRTWtDcjW092eXpkXcQPMHNw/4jupKQJV4/camlswMD0Vv7nkpDZptHLuViB2Vm04S0J5jtf5Va036WtnT4KvJYSqxgD+smg3tubX2MegoqEZx48ebF9jAFDji39uPrkiD7e/tx2ltX6MGNgPz3x3FoC2JYVV+yrwcNR198h1Z7RZX6zqz9fXH8a7jtJEe6J7aFslIADI6JcSUX3kPK7Oe8ewAa3VmVbGIKNfqn2d3/7edvzwxexOjUw+XNUEXyCEjH7GrbonSgtpSf+FI4h1AjuLazVNLW3qa+//MAdTMwfhyyeNBmBcPNc/vyHueqsbAyio9qHebzTAdtS706pDjK4/jVZS50fW4P6YNHIgAONGsdisLhhrBoWKdgLLxaeMwYIdJZg7ZSTWHahCZWMzVBUiEnNQ16HKJlQ3BnDb21uxen8lZk4cjvc3FyK3rAFffySyMXJQ/zT075dqVwM4hcOKlBRBQVXsBsQRA9PbfMfNaM3xI4zj+c+l+7BwRwn+fOVpOH3icLsUcOwoY781NAcj+uvnljfYx6i03o+WUBjBkOJARewugf9cus++OVY3BZAaJ3g896nR2FtW34yWUDjiAu6XKmgJKfaW1iNzSH/8+q2tWBOj55q1n+58fzuWOYJw9HY7q9mmjx2K7YW1WL63HE+vOoDRQ/pj+IC2+9gqjTQGQnj44324dvZE+285MUrmsGovAAAThUlEQVQxQGtVUvZdF2Dz4Ro899kB7C1twKPLjEFWIwe1/s7Fp4zBBPMaA9p2m43O5S7aUYJgWPE/cydhatYgAIho3LY882nkgLrooAi0bVNoCYXtzhWvrT+MX1xwAj4/LTPmNgJoM8X8sSMHoqbJyNGPGJhuVwH7AiH80xGkanyRQaHEPPes87mhOdipnlk/OHcK5kwZgRtfMTpaHKpsgq8lhGmjB2NfWQMq6gMYO2xAB2vpmqQGBRG5CMBDAFIBPKOqD0T9vT+AlwCcBaASwDWqejBZ6TlmiFHdsL2w9QKo9bWgzhfEZaePw4+/MNW++f1tyR7sK6tHKAw86silxVLd1GJXDXxu6iis3l+JoRlpcbtoWj1HPtpeDH11Ix669oyYU1GU1vlxzLAMjBqUjvS0FPxnaxGyza6NY4cNwKD0NDsHG8sj150BfzCMsCpOu3cxNh+uwc/f2ILHvnVmzJLC6v2VOMMx5cW2/Jq4faNHDUrHkIw0u645It31fowdNiBuo+mwAW1LCgcr4m9HPJNGGjeSD83uq798aysG9EvFjHFDkTk43f6divrIbbACgnWhldT68c3HPkNFjAAHGD18AOPmm1Nch2CcqP/0qtab17Q7F0Rs5xkTR2D9wSrc8d72dkuURTU++AKhiHM02oThA7A1v7W0O2OcERT+z5weu6apBcMGxtjHjiqivy/ZG5HLjZ6V83/mTsJr6w7j9fX5OGnMEGQO7o8JZhD+ZHdrTveU8cOADUbp8MxjR2CcIygUVPvw1oZ8VDUFsLe0HmOGZkT8htX+ceyogfaDixLpvXPM0MigMHJQOmp8LSivb8bA9FQU1fgwID3V/vvGQ9X49rPr8JtLTsINXzyuzfqi26SM38gAYASFrCH9sXRXGT7LNRruX17b+sxm53k7ekgGDlY24Qt//sTORJTV+yNGn594zJCY18yUzEG45xvTcfbUUXbnAaB1f0zNGmScq3V+nIph7e6frkpa9ZGIpAJ4FMDFAKYDuE5Epkct9kMA1ap6PIB/AHgwWekBjHpwABFVFPO3F6OkzriJnTJ+GBbd+kVcNGMMdhTW4Y/zd+PBhbtjDsw5dXzrgfn121vx8tqDSBFgjln/H+uijGX+9hJsyY/9IJHSOj+OGdIfKSmC08YPswMCAJw0dghGDU6PaCh1+tqpY5GWmoLB/dMwNKMfxgzNsHte/WnBroQa5tYeqLSrgKJl9EuFLxCKmcPfUViH7INVdhtItIHpbfMi0bnmR647A7d99cSY3z91/DBcf85kTBzRNsfkawkh+1A1soZk2Deau/8d+xm75xw3CgBwycOr4gYEp7lTR9qvb/vqifjW3En2e2uE+J+vOM3+zNnWYZVa2gsIUzMHoay+GSffvRAVDc34yXnH4cErTsVFM8ZELDchartnjI+8SQRC4ZiBtyKqI8GCOG1gAHC+WUoGgEtOHWukL2sQ0tNS8Fa2USXz1v9+LmIfTBs92N7nlv97dxseWLAb720qxGNxZjP9wrQsDHTcxAFg271fsUvI0ay2JMtxWYOwbE8ZZv9hKWbcswgX/mMl3trQdrr4P87fjZ++uiliMOHOoloc95v5bZYd5SiNWHmAbz2zDg9Ftbs5q4CnZBqZFGeVm78ljP2OsR1fmXFMxPcnm+fF0AH9cN6Jo5HRLxVDM9oeu5PHDgUAlNQmf/6jZLYpzAGQq6p5qhoA8AaAy6KWuQzAi+brdwCcL5K8SepjrdrKPWUONorBJ44Zgt9/85SY379weusBvfWCaVh/5/n2+89yK3Fc1mD7gMbrtjhjnHFwLzi5dV1XPbEGZ92/BGfdvwRz/7gUF/1zJeZtLcLe0gZMNk+03102w17+pR/MwdCMfhHF6FPGD434nYeuPT3i/dWOqgLnSRstc3A6/vDNU3D+SaPxWW5l3Jvl+OED7CeQRfvxS9m48ok1+PeWIgzp3zYARN8AfnHBCW2W+cbMcTjbcRN2euWHc3HvpTMwYlDbKhJLIBjCIDP4ZEcNHJuSadzcrppl7JN4UzQ4T5dpowfj3ONaqx+OHz0Yv/pKa9D68Gefxxs3nI2rZ0/Eytu+jOnmRWzN9pmWwBN5Zk4cHvH+yyeOxjWzJ+GBK06N+NyZGweM3KdT/7SUiDrxeKyAPiLGstYNDgBuOX+aud5UnDFxOGp9LeifloKzjh0RcU1Z1bOPXHcG3r/pHNx/uXEdDeiXio9u+by93D+umYnPTTUC8rWzJ2JK5qA2GYWhGf3sevRoA9PT8NqP5trvxw0f0KbK9uGop8FdceYEAEbp/Py/rcBF/1yJKx9fjZ+/sSXmbwzo13qOWlVC1g08HqsKzDLWDF6Hq5pw+sThuOqsCbj4lLERy1jX9wDHtk4fOxTnHj8Kj3/rTPuzk8YMwds3fg5fP21cu2noDsmsPhoPwBmuCwDMjbeMqgZFpBbAKADtd8vpgpPGDLH7/E8aOdCu4pgwovWAZw7uj/dvOgdvbyzAa+sO4/LTx2Hu1FG4bs4kfLC5ELe+uQVnHTsCwwem47nrZ+GJFXlYf6AKMycOt3M3Gf1S2/z23Ckj8ch1Z+BQVRNSRFDe0GxXA1gNxlMzB2F3ST1ued14fsDsycaNcfrYofjCtEwcMzQDXzwhCwAwedQgezuev34O5m8vxqHKJnz5pKw2N6FLZ45r07to5sTh2Jpfg+evn43V+yvw9KoDWHjrF5E5uD9Onzgcmw5XY8TAdPzyKydgzf5KXDpzHNLTUlDdFMCXThiN40cPxreeWQfAyCmHwor0tJSIUbRfnznOHk19xZkTMCQjDdOOGRyRjuvmTMQ5x4/C6txKpAgwy9zmiSMiL8IPf/Z5TM1qvYH0S03BnCkjYw4KmzRyIEY4GrQnjhyAsUMHIPtQFd6/6Rw0BUIYOywDcyaPRFGtD1+YlhnRSP7Z7f+FwelpmHnfYqQIsOSXX0KdvwWXnz4Op08cjgtPPgYpKYI5U0biSydk2Tk5AJg0aiAuPmUMcorrcMHJx2DhzhJ8/bSxuHrWBDy6bH+bAXOXnz4OH2wpwi3nT0OKCLYV1GB/eYNdGh0+MB1PfPtMrNxXgdfWHcY0szF35KB0XDpzHM6c1BpM/u+iE/HVGWMiAu8xQ/vbje9//Oap+N1/dmLiyIHILWtA5uD+WH7beXhwwW6s3Fdu90IabVb1RD9n/L/PHI91B6ow7ZjBdunolvOnIRAM2wHiGzONG9fkUYOwYHsxrpo1wQ6SAPC1U8chv8qHNXmVOG2CkfbUFMHNXz4e044ZbB//Z783G/d9mIMlOaUYlJ6KG790HFrMqp5zjs/E/Zefgv9sLcJXZ4zBv7cU4fpzJkeMQzrnuFGYNXkkLj5lDEJhxbubCsz9kREx7mf62KFISxXMnjwSOUV1WJNXienjWtP7t6tm4rnPDuD/XXQSLja7rf7XSaPtrudWNeTsKZGZmM8fn4m3zYbuq2ZNwLfmHtumXeW4rMFYvqc8IrCnpAhe/dHZAIyM6JKcUkwcORAzxiW32shmjUDs7n8AroTRjmC9/w6Af0UtswPABMf7/QAyY6zrBgDZALInTZqkXVFS69M1+yv0g80FWljdpB9tK9KHlu7VQDDUZtlAMKSLdhRrMBRud52+QFCfXJGrhdVNGgqF9cXVB7SmMaAvrT6gxTU+XZ1boaW1PvW3BNt8d+XeMn1x9QF9ZlWe7iysVVXV3cV1ujSnRBfvLNFwOP5vl9b59J5/79DtBTUJbfvW/GotqG7Svy3eo5/sKtVQKKy+gJGmcDisLTH2QUfq/S26fE+ZVjU0a0W9X9flVer8bUVa1dCsT67I1fJ6v763KV9X7CmL+N7CHcW6Yk+Zvr7uULvrf+7TPD1c2Rhz36mqNvhb9NN95brxUJU+sypPy+r8Om9Lodb6AhoOh3X5njLNK2/Q8nq/tgRDWtXQHPe31uVVanVjc8R+yCmq1ZJaXyf2iMEXCOrr6w5pc0tI63wB+/NwuHWfl9b5dOGOYg2Hw9rU3Lp9/pZgxHec3914qEpVVfeV1mtpXWu6Nhyo1JV7I/fx4p0lWu9vUVXj2B8ob7DX428J6tMr97c5d/KrGnVfaZ2qqj69cr/uLamL+HswFNYnV+TqlsPVndshqrrlcLUuN8+DqoZmfXZVXkLn3Lq8St3czu+Fw2H7eH+41Tj38sxtdS7zytqDWl7v1wZ/i67cW6af7C7VP36UE3FuBUNh3V5Qo+FwWMvr/Xq4srHNOgqrm1RVNftgpW7Lr9FQKKyNzS0aDof18eW5ujq3Qv/4UY42NrfomxsO6+PLc+1jrqr6ye5SLapp0lV7yzUYCuuzq/LinmMltT5dsL24w32UCADZmsC9WzRJT/MRkc8BuFdVv2q+v8MMQn9yLLPIXGaNiKQBKAGQpe0katasWZqdnZ2UNBMReZWIbFTVWR0tl8w2hQ0AponIFBFJB3AtgHlRy8wD8D3z9ZUAPmkvIBARUXIlrU1BjTaCmwEsgtEl9TlV3Ski98EoxswD8CyAl0UkF0AVjMBBRES9JKnjFFR1PoD5UZ/d7XjtB3BVMtNARESJ61PTXBARUfsYFIiIyMagQERENgYFIiKyMSgQEZEtaYPXkkVEygEc6nDB2DKRxCk0jkDcXm/rS9vbl7YVSM72HquqWR0tdNQFha4QkexERvR5BbfX2/rS9valbQV6d3tZfURERDYGBSIisvW1oPBUbyegh3F7va0vbW9f2lagF7e3T7UpEBFR+/paSYGIiNrhyaAgIheJyB4RyRWR22P8vb+IvGn+fZ2ITO75VHafBLb3ehEpF5Et5r8f9UY6u4OIPCciZSIS88HLYnjY3BfbROTMWMsdLRLY3vNEpNZxbO+OtdzRQEQmisgyEckRkZ0i8vMYy3jm+Ca4vT1/fBN5Es/R9A/GNN37AUwFkA5gK4DpUcvcBOAJ8/W1AN7s7XQneXuvR9RT747WfwC+COBMADvi/P0SAAsACICzAazr7TQneXvPA/Bhb6ezm7Z1LIAzzddDAOyNcS575vgmuL09fny9WFKYAyBXVfNUNQDgDQCXRS1zGYAXzdfvADhfnE8gP7oksr2eoaorYTx7I57LALykhrUAhovI2HaWP6IlsL2eoarFqrrJfF0PYBeM57g7eeb4Jri9Pc6LQWE8gHzH+wK03dH2MqoaBFALYFSPpK77JbK9AHCFWdx+R0Qm9kzSekWi+8NLPiciW0VkgYjM6O3EdAezSvcMAOui/uTJ49vO9gI9fHy9GBSorf8AmKyqpwFYgtZSEh39NsGYvmAmgEcAfNDL6ekyERkM4F0At6pqXW+nJ9k62N4eP75eDAqFAJw54QnmZzGXEZE0AMMAVPZI6rpfh9urqpWq2my+fQbAWT2Utt6QyPH3DFWtU9UG8/V8AP1EJLOXk+WaiPSDcYN8VVXfi7GIp45vR9vbG8fXi0FhA4BpIjJFRNJhNCTPi1pmHoDvma+vBPCJmq06R6EOtzeqzvVSGHWXXjUPwHfNXipnA6hV1eLeTlSyiMgYqz1MRObAuKaPygyOuR3PAtilqn+Ps5hnjm8i29sbxzepz2juDaoaFJGbASyC0TPnOVXdKSL3AchW1XkwDsTLIpILoxHv2t5LcdckuL23iMilAIIwtvf6XktwF4nI6zB6ZGSKSAGAewD0AwBVfQLGM8EvAZALoAnA93snpd0jge29EsBPRCQIwAfg2qM4g3MugO8A2C4iW8zPfgNgEuDJ45vI9vb48eWIZiIisnmx+oiIiFxiUCAiIhuDAhER2RgUiIjIxqBAREQ2z3VJJYpFREIAtjs+ulxVD/ZScoiOWOySSn2CiDSo6uB2/p5mzoNF1Kex+oj6LPM5E/NE5BMAH5uf3SYiG8zJA3/nWPZOEdkrIp+KyOsi8mvz8+UiMst8nSkiB83XqSLyF8e6/tf8/DzzO++IyG4RedUxYnW2iKw2Jz9bLyJDRGSliJzuSMenIjKzp/YR9T2sPqK+YoBj1OgBVf2m+fpMAKepapWIfAXANBjTkQuAeSLyRQCNMEa9nw7jmtkEYGMHv/dDGFMwzBaR/gA+E5HF5t/OADADQBGAzwCcKyLrAbwJ4BpV3SAiQ2GMYH0Wxgj0W0XkBAAZqrq1S3uCqB0MCtRX+FT19BifL1FV63kFXzH/bTbfD4YRJIYAeF9VmwBARKLn0orlKwBOE5ErzffDzHUFAKxX1QJzXVsATIYxfXuxqm4AjInQzL+/DeC3InIbgB8AeCHRDSZyg0GB+rpGx2sB8CdVfdK5gIjc2s73g2iths2IWtfPVHVR1LrOA9Ds+CiEdq5DVW0SkSUwHi5zNbw9wy0dAdimQNRqEYAfmPPbQ0TGi8hoACsBXC4iA0RkCIBvOL5zEK036iuj1vUTc2pkiMgJIjKond/eA2CsiMw2lx9iTusOGNOdPwxgg6pWd2kLiTrAkgKRSVUXi8jJANaYbb8NAL6tqptE5E0Yz78ugzFdueWvAN4SkRsAfOT4/BkY1UKbzIbkcgCXt/PbARG5BsAjIjIARnvCBQAaVHWjiNQBeL6bNpUoLnZJJeokEbkXxs36rz30e+MALAdwkqqGe+I3qe9i9RHREUxEvgvjub13MiBQT2BJgYiIbCwpEBGRjUGBiIhsDApERGRjUCAiIhuDAhER2RgUiIjI9v8BBe9DYOFR1G0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "period=2\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": [
    "Finally, we are ready to combine what we have learned, and use lomb scargle to search for the best period of our lightcurve, and automatically plot the phase folded lightcurve of the object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEKCAYAAADuEgmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF95JREFUeJzt3X+w7HVdx/HXq3NN0CUOeG9K/PBC0RhYrnjQ0n6siQlUXM1sQCswlKzoxzjThMNE5T/Rj5l+MjWEjjrTCIiZmNdBFDdFAjkwy+/Q68WEG8UJBdvRME7v/tjv4pe9u2c/Z3e/P3b3+Zg5c7/7/X53932/573n/f18Pt/vZx0RAgAgxbdVHQAAYH5QNAAAySgaAIBkFA0AQDKKBgAgGUUDAJCMogEASEbRAAAko2gAAJLtqDqAWdu5c2fs3r276jAAYK7cdttt/xURu8btt3BFY/fu3VpfX686DACYK7b/LWU/uqcAAMkoGgCAZBQNAEAyigYAIBlFAwCQjKIBAEhG0QAAJKNoAACSUTRqoNVqqdVqzc3rYv4UmQvk2XKhaFSs1Wqp0+lUHQYWXKfTIc8wExQNAEAyisacoksAZSDPMIiiAQBItnCz3C6D/jhIs9msOhQsKFoXGIWWxpxh4BxFI8ewFYpGTc2qL7nT6XDWiJHIM2wXRaNC/TO6brebfGbX3x9I1Wq11O121e12k/+wbycnsVwqHdOw/W5JPyXpkYh44ZDtlvQXks6U9HVJ50XE7eVGWY7+B7rdbh+0Lf9Bp2BgGv0WwVZ51ul0tLm5WW5gmBtVtzTeI+n0LbafIenE7OcCSX9TQkyl29zcHPkhpX8Zs7BVjknkGdJVWjQi4tOSvrLFLnskvS96bpa0avuocqIDAAyquqUxztGSHsw9fihbBwCoQN2LRhLbF9het72+sbFRdTgAsLDqXjQOSDo29/iYbN3TRMTlEbEWEWu7du0qLTgAWDZ1vyP8WkkX2r5S0sskPR4RD1cc09S4nh1F6+fYsKukgGlUfcnt+yW1JO20/ZCk35P0DEmKiL+VtFe9y233qXfJ7ZuriTQNH1SUgTxDlSotGhFxzpjtIenXSgqnNP1LG4fNHUUrBLMw7vJZ8gyTqvuYxsLijluUgfsvMGsUjSXEdySgDOTZYqr7QPhC4QOEMvTnmmo0GlWHggVE0ZhTzWaTgVAUjhzDILqnKrK5uVnZ5IOdTuepieto/Sy2qmZFzo+lkGeLhaIBAEhG99ScmbSfmmv7sR0rKysTfZ0wLYrFR0tjDjB+gTKQZ0hB0ZiRSa6HbzQaeuyxx7SysrKt5/HhXl6T5tnhhx++rTxrNBrkGIaiaMzAtDdQNRqNbReO7egPfKeuRz3NIs8m6XJKMepmVXJs8TCmUQPNZpMPFgpVVLHA8qGlMUfa7fa2Pvzj9u92u3znOA7SbDa3nWdbXaDBlDmLhZbGnKGfGUUjx7AVWhoAgGQUDUjirl0UjxxbDBSNGZpl32273aabAEORZ6gSRWPBjRsM39zc5OwPUxs3eM5g+OJgIHwCRU7JwVkfpOKnfSHPMCmKRgk4k0fRmFsMZaF7qgLNZrPQO8CBIu/+xnKjpVERPtQoAzmGWaNoTKj/JUaz6A5ot9sju7CK7G5oNBrcEV5zs+p26j+/qjzD4qBo1ETRfdHMjAup+DwbV5ww/xjTqMB255ACtouTBBSForGEtjshHbBdnBgtLorGEms2mwf1NzPVA2Zp1AkKeTa/GNMoEd0FKAN5hiLR0gAAJKOlUZEqzwa5wmV5kGeYNVoaM7SoN+zR/1wv5BmqREtjCdDHjaKRY8uj0paG7dNt3297n+2Lhmw/z/aG7U7285Yq4gQA9FTW0rC9IukySa+W9JCkW21fGxH3Dux6VURcWHqAAICDVNnSeKmkfRGxPyK+KelKSXsqjAcDWq2WVldX+fIcFIo8my9VFo2jJT2Ye/xQtm7Q623fafsa28cOeyHbF9het72+sbFRRKxb4u5XlIE8Qx3U/eqpj0jaHRE/IOl6Se8dtlNEXB4RaxGxtmvXrlIDBIBlUmXROCAp33I4Jlv3lIh4NCKeyB5eIeklJcX2FC4DRBnIM8yLKovGrZJOtH287W+XdLaka/M72D4q9/AsSfeVGB8AYEBlV09FxJO2L5R0naQVSe+OiHtsv1PSekRcK+k3bJ8l6UlJX5F0XlXxAgAqvrkvIvZK2juw7pLc8jskvaPsuJYRUz6gaHzHx2Ko+0D43ODKFpSBPEPVmEZkhjiLQhnIM1SJlsYEFuVsb1H+H4uo3W4vTHEgzxYLRWNCi/ShTtHpdBjvqAB5hrqhaAAAkjGmseSW6SwW1SHPFgctjTE6nQ4TqaFQrVaLHMPcoKWRqN/POskZ07yfZTGIWZ5J82zec0wiz+YFLQ0AQDKKRoJut0v3AQpFjmFe0D2FJPxBQ9H6OUY3Vb3R0kiwubmpbrdbdRi1xtTe0yHH0pBn1aNoAACS0T2FkfJX5KyurlYXCBZafoZlukHrj6KBp0m9dHOaS5AB8mx+0T2VaHNzs+oQsODIMcwDigYAIBndU5gKV7KMRtfK7JBno5WdZ7Q0cricb3ocw/E4RtNh+vRqUTQAAMkoGgCAZBQNAEAyBsIxNeYMQhm48a8eaGmMwR9ClGVlZaXqEICxKBoAgGR0T+Xkm79c0vd0tLhmoz+/UrPZJMeGaDab3NdSc7Q0AADJaGmM0W63tWPHjqWfF4izv+L0c2zZkWPzgZYGACAZpzcJGo0G36qGQjUajapDAJJU2tKwfbrt+23vs33RkO3PtH1Vtv0W27vLjxIA0FdZ0bC9IukySWdIOknSObZPGtjtfElfjYjvkfRnkv6o3CiRqtvtHnT1GVcHYZa63e5BLX7yrHxVdk+9VNK+iNgvSbavlLRH0r25ffZI+v1s+RpJf23bERFlBir1ug8YqEPRyDHUXVJLY0gLQLZbU7730ZIezD1+KFs3dJ+IeFLS45KeM+X7AgAmlNo9dbXt33HPobb/StIfFhnYdti+wPa67fWNjY2qwwGAhZVaNF4m6VhJN0m6VdK/S3rFlO99IHvNvmOydUP3sb1D0uGSHh18oYi4PCLWImJt165dU4YFABgldUzjfyV9Q9Khkg6R9EBE/N+U732rpBNtH69ecThb0hsH9rlW0rmS/kXSz0q6ocjxjFFTZTCFxniNRoPjlKg/VUZ+AJdjl4Y8q15qS+NW9YrGqZJ+RL0rnT4wzRtnYxQXSrpO0n2Sro6Ie2y/0/ZZ2W7vkvQc2/skvV3SQZfloj46nQ7TV6NQg1fpoXypLY3zI2I9W35Y0h7bvzDtm0fEXkl7B9Zdklv+H0lvmPZ9AACzkVo0HrF93MC6f551MJgf/UtDuUYeRSLP6ie1aHxUUkiyemMax0u6X9LJBcVVK1w7j6KRY5gXSUUjIr4//9j2KZJ+tZCIAAC1NdE0IhFxu3qX4QIAlkhSS8P223MPv03SKerdqwEAWCKpYxqH5ZafVG+M44OzDwcAUGepYxp/UHQgWAyDN671/122gd5l/X+XYfB7xDudjlqt1tId66pybMuiYfsj6l01NVREnDVqG4AeCgjKUFaejWtp/Gmh7465129VcJcuitRut7W6uso3aNbAuKLxQER8uZRIMLcoHChDs9kkx2pg3CW3/9hfsM3ANwAsuXFFw7nlE4oMBABQf+O6p2LE8sJj0BLbMckgJDmG7arDRRXjisaLbH9NvRbHodmysscREd9RaHQAgFrZsmhExEpZgQDzjAFaFK1/sUnVX0KVekc4gAF16CrA4qvbtPAUDczE4B/OYd+DwB9ZTGNY3gx+bS45VjyKRg6JhqKRY5h3FA1gCv15j4Ai1WnMbKLv0wAALCeKBgAgGUUDAJCMMQ1gmxjDQBnqmmcUDdTKvFwyWaeBSWwfeTY5igZQsLr/YcJiKCvPGNMAACSjaKBSrVartn23Uv3jQ5q6/x7rHl8e3VPADPA1pChDt9utfJyDlgYAIBlFA9gGvgsdZahznlE0UGt17+utQ3cBplfnPOt2u7Xq/qxkTMP2kZKukrRb0pck/VxEfHXIfpuS7soefjkiziorRswGl5uiaORYuapqaVwk6ZMRcaKkT2aPh/lGRDSzHwoGAFSsqqun9khqZcvvldSW9DsVxYKSzcvduHl17brAaORZMapqaTw3Ih7Olv9D0nNH7HeI7XXbN9t+7agXs31Btt/6xsbGzIMFxmk0Gmo0GlWHgQXXaDQW9zvCbX9C0vOGbLo4/yAiwnaMeJnnR8QB2ydIusH2XRHxxcGdIuJySZdL0tra2qjXQsWGfQUsMEvz1KqYV4UVjYg4bdQ22/9p+6iIeNj2UZIeGfEaB7J/99tuS3qxpIOKBgCgHFWNaVwr6VxJl2b/fnhwB9tHSPp6RDxhe6ekV0j641KjBCbEGS+KVlWOVTWmcamkV9v+gqTTsseyvWb7imyf75O0bvsOSZ+SdGlE3FtJtAAASRW1NCLiUUmvGrJ+XdJbsuWbJH1/yaEByZrNJjf2oXB1u8CCO8IBAMkoGgCAZEyNjtroT9JW9XXoWFxc7j09WhoAgGS0NDA36jwtBK2jxUCOjUdLA5hCs9mszYcZi6vZbNamkNHSQKUGPwh1+nBgcQxOYUOOTY6WBgAgGS0NYJvojkIZ+nlWtxtIaWmg9jqdDpdKonCdTqd2f6DriKIBAEhG9xTm1iwGNUe9xqiWzagBVVpCi2vaPNvq+Sl5VreLQygamDt1+ANdpw8xikGeDUf3FObOdvqeW63Wtj/8rVZLN954I/3bSy41zybJMUlaXV3VjTfeOEFk1aJoAACSUTQAAMkoGlhY/VlzR21L6VJot9u17FdGPWyVYynb++YpzxgIR6212+1aDEiOMy8feAzXz7O6j2PVIc8oGihdHRIfi488KwbdUwCAZLQ0UBucGaJo5Nj0KBpYOnytLIo2D+Nwk6J7CgCQjKIBAEhG9xSQgL5wlGEe8oyWBgAgGUUDtTfN3bLdbrf2N2yhHtrt9sQXRyzTF4XRPYW50+12qw4BS4A8G46WBgAgGUUDAJCskqJh+w2277H9f7bXttjvdNv3295n+6IyYwQAHKyqlsbdkn5G0qdH7WB7RdJlks6QdJKkc2yfVE54AIBhKikaEXFfRNw/ZreXStoXEfsj4puSrpS0p/joMM9WV1e1uro69es0Gg2mGcFIs8yzebg3I6/OYxpHS3ow9/ihbB0AoCKFXXJr+xOSnjdk08UR8eEZv9cFki6QpOOOO26WL40a6Z+R7djBleIoxryd9VehsE9fRJw25UsckHRs7vEx2bph73W5pMslaW1tLaZ8XwDACHU+ZbtV0om2j1evWJwt6Y3VhoRlMC9fMYv5Nq9jZlVdcvs62w9J+iFJH7V9Xbb+u2zvlaSIeFLShZKuk3SfpKsj4p4q4sX86Ha7SXfyNptNuiIwsdQ8m2ZqkrqqpKURER+S9KEh6/9d0pm5x3sl7S0xNADAFurcPQVUhlYIijavOVbnS26BoRqNhhqNxtj92u120n7AMCl5tojdT+PQ0sDcWbYPKapBng1H0cDcSmne88HHtMbl2bJdbUfRwFKb135lzI9FyzHGNAAAyWhpYOks2pkf6meRc4yWBgAgGUUDAJCM7iksFO7LQBmWOc9oaQAAktHSwNxZ5EFG1Ad5NhwtDQBAMooGACAZRQMAkIwxDSw0+qVRhmXKM1oaAIBktDSwUJjVFmVY5jyjpQEASEbRAAAko2gAAJJRNAAAyRgIx0JZpksfUZ1lzjNaGgCAZBQNAEAyigYAIBlFAwCQjKIBAEhG0QAAJKNoAACSUTQAAMkoGgCAZI6IqmOYKdsbkv5tgqfulPRfMw5nVuoaW13jkuobG3FtX11jq2tc0mSxPT8ido3baeGKxqRsr0fEWtVxDFPX2Ooal1Tf2Ihr++oaW13jkoqNje4pAEAyigYAIBlF41surzqALdQ1trrGJdU3NuLavrrGVte4pAJjY0wDAJCMlgYAINlSFA3bp9u+3/Y+2xcN2f5M21dl22+xvTu37R3Z+vttv6bkuN5u+17bd9r+pO3n57Zt2u5kP9fOMq7E2M6zvZGL4S25befa/kL2c27Jcf1ZLqbP234st62wY2b73bYfsX33iO22/ZdZ3HfaPiW3rcjjNS6uN2Xx3GX7Jtsvym37Ura+Y3t9lnElxtay/Xjud3ZJbtuWeVBwXL+di+nuLK+OzLYVdsxsH2v7U9nfhHts/+aQfYrPs4hY6B9JK5K+KOkESd8u6Q5JJw3s86uS/jZbPlvSVdnySdn+z5R0fPY6KyXG9UpJz8qWf6UfV/a4W/ExO0/SXw957pGS9mf/HpEtH1FWXAP7/7qkd5d0zH5U0imS7h6x/UxJH5NkST8o6Zaij1diXC/vv5+kM/pxZY+/JGlnhcesJemfps2DWcc1sO9PS7qhjGMm6ShJp2TLh0n6/JDPZeF5tgwtjZdK2hcR+yPim5KulLRnYJ89kt6bLV8j6VW2na2/MiKeiIgHJO3LXq+UuCLiUxHx9ezhzZKOmdF7Tx3bFl4j6fqI+EpEfFXS9ZJOryiucyS9f0bvvaWI+LSkr2yxyx5J74uemyWt2j5KxR6vsXFFxE3Z+0rl5ljKMRtlmvycdVxl5tjDEXF7tvzfku6TdPTAboXn2TIUjaMlPZh7/JAOPtBP7RMRT0p6XNJzEp9bZFx556t3BtF3iO112zfbfu2MYtpubK/PmsDX2D52m88tMi5lXXnHS7oht7rIYzbOqNiLPF7bNZhjIenjtm+zfUFFMf2Q7Ttsf8z2ydm6Whwz289S7w/vB3OrSzlm7nWhv1jSLQObCs+zHZM8CeWy/fOS1iT9WG718yPigO0TJN1g+66I+GKJYX1E0vsj4gnbv6xeS+3HS3z/cc6WdE1EbObWVX3Masv2K9UrGj+cW/3D2fH6TknX2/7X7Cy8LLer9zvr2j5T0j9KOrHE9x/npyV9NiLyrZLCj5nthnqF6rci4muzfO0Uy9DSOCDp2NzjY7J1Q/exvUPS4ZIeTXxukXHJ9mmSLpZ0VkQ80V8fEQeyf/dLaqt31jErY2OLiEdz8Vwh6SWpzy0yrpyzNdBtUPAxG2dU7EUeryS2f0C93+GeiHi0vz53vB6R9CHNrms2SUR8LSK62fJeSc+wvVM1OGaZrXKskGNm+xnqFYy/j4h/GLJL8XlWxIBNnX7Ua03tV6+roj9odvLAPr+mpw+EX50tn6ynD4Tv1+wGwlPierF6A34nDqw/QtIzs+Wdkr6g2Q4EpsR2VG75dZJujm8NuD2QxXhEtnxkWXFl+71AvQFJl3XMstfdrdGDuj+ppw9Qfq7o45UY13HqjdW9fGD9syUdllu+SdLps4wrIbbn9X+H6v3x/XJ2/JLyoKi4su2Hqzfu8eyyjln2f3+fpD/fYp/C82ymCVDXH/WuKPi8en+AL87WvVO9s3dJOkTSB7IPz+cknZB77sXZ8+6XdEbJcX1C0n9K6mQ/12brXy7pruzDcpek8ys4Zn8o6Z4shk9JekHuub+UHct9kt5cZlzZ49+XdOnA8wo9ZuqdcT4s6X/V6y8+X9LbJL0t225Jl2Vx3yVpraTjNS6uKyR9NZdj69n6E7JjdUf2e764gBwbF9uFuRy7WbnCNiwPyoor2+c89S6SyT+v0GOmXtdhSLoz9/s6s+w8445wAECyZRjTAADMCEUDAJCMogEASEbRAAAko2gAAJJRNIAJ5GbMvdv2B2w/y/buUTOjAouCogFM5hsR0YyIF0r6pnrXygMLj6IBTO8zkr4nW16x/XfZ9x183PahkmT7rbZvzSbf+2A22Z1svyFrrdxh+9PZuhXbf5Ltf2c2txdQCxQNYArZXGVnqHf3rdSbUO+yiDhZ0mOSXp+t/4eIODUiXqTelNbnZ+svkfSabP1Z2brzJT0eEadKOlXSW20fX/z/BhiPogFM5lDbHUnr6s2J9K5s/QMR0cmWb1NvDiNJeqHtz9i+S9Kb1JvXTJI+K+k9tt+q3pcLSdJPSPrF7PVvUW+a/jrN7oolxtTowGS+ERHN/Ire93bpidyqTUmHZsvvkfTaiLjD9nnqfSudIuJttl+m3kRzt9l+iXrzB/16RFxX5H8AmAQtDaAch0l6OJva+k39lba/OyJuiYhLJG2oN331dZJ+JdtXtr/X9rOrCBoYREsDKMfvqtfVtJH9e1i2/k9sn6he6+KT6s2Qeqd63Vq3Z187vCGp7G8aBIZillsAQDK6pwAAySgaAIBkFA0AQDKKBgAgGUUDAJCMogEASEbRAAAko2gAAJL9P/C7rpuaydqNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "period=0.5\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",
    "    return 1/frequency[np.argmax(power)]\n",
    "\n",
    "def plot_best_period(t,y,dy):\n",
    "    phase_fold(t,y,dy,lomb_scargle(t,y,dy))\n",
    "\n",
    "plot_best_period(t,y,dy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's apply our algorithm to real data. Here, we import a lightcurve of a cepheid, and search for its period. Note, that if you are working with a lightcurve that includes errors, you should add the dy column back in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXuQLFd937+/mTuLNLPXlm5fuUIsdhZiKAIUFuhiHqFswdo8rrFwYieFa68iBGHDyHEppeCyqa2ojCpbmCJVzo15KGtM+cIMNs8ksoILI5DLWFgiq7eEbSxg9yIHG917hdHVgiTu/vLHdI96e/pxuud0Tz++n6pTO9t9+syZ09/u33n+jqgqCCGEEABozTsDhBBCygONAiGEkAk0CoQQQibQKBBCCJlAo0AIIWQCjQIhhJAJNAqEEEIm0CgQQgiZQKNACCFkwoF5ZyAthw8f1uXl5XlngxBCKsUdd9xxSlUvSopXOaOwvLyMra2teWeDEEIqhYjsmMRj9xEhhJAJNAqEEEIm0CgQQgiZQKNACCFkAo0CIYSQCTQKhBBCJtAoEEIImUCjQAghZAKNAiGEkAk0CoS4jEYjLC8vo9VqYXl5GaPRaN5ZIjWhStqqnJsLQvJgNBphbW0Nu7u7AICdnR2sra0BAFZXV+eZNVJxqqYtUdV55yEVR44cUfo+IrZZXl7Gzs60a5h+v4/t7e3iM0RqQ1m0JSJ3qOqRpHjsPiIEwMmTJ1MdJ8SUqmmLRoEQAEtLS6mOE2JK1bRFo0AIgI2NDXS73X3Hut0uNjY25pQjUheqpi0aBVJrTGd9rK6uYnNzE/1+HyKCfr+Pzc3NUg4EkvJgoq/KaUtVKxUuvfRSJXYZDofa7/dVRLTf7+twOJx3lqwwHA612+0qgEnodru1+X1VgfoqBwC21OAdO/eXfNpAo2CXqgk7Df1+f9/v8kK/30+8tq4vsqKhvqaZl7ZoFIgRs7w4y46IhP42EYm9rs4vsqKhvvYzT22ZGgWOKTScqk2XS0PWWR/r6+uThUYeu7u7OHbsWOlXo5YN6ms/VdBW7kZBRNoicpeI3BRy7loR+aqI3CsiXxCRft75Ifup2nS5NBw9ehQisu+YyayPuBfWzs4Ojh07hsOHD5fiAS471Nd+krR1xRVXQETmayBMmhOzBADXAvgYgJtCzr0KQNf9PADw8aT02H1kl7p2lYT9LhHRwWAwFS/YvxvV5REMXjlx/CEa6mu/NhzHMdKWV06DwcCatlCGMQUAFwP4AoBXhxmFQNwXAbg1KU0aBfvU8aVm0pcd9cIaDAa6sLBg9OA6jlPLl55NqK+nznc6HWNteYbGlrZMjUKuvo9E5FMA3g3gIIB3qOobYuK+D8Dfq+p/iUuTvo+ICa1WC2HaFhHs7e0BiPdJc+rUKTz22GOZv58+k+rNLPpyHAePPPLIJF5asmpr7r6PROQNAL6jqncYxD0G4AiA90acXxORLRHZevjhhy3nlNQRk77suEHQWQxCXNqkHsyirzNnzmQ2CHHp2iLPgeZ/AeByEdkG8EcAXi0iw2AkEflZAOsALlfVx8MSUtVNVT2iqkcuuuiiHLNM6oKJa4FZB0G73S4cx5kpDVJNitBXcBA77fWZMeljmjUAuAzhA80vAvB1AM82TYtjCsSUpL7suEHQpAHBdrs9GWTmmEIzyVNf/X5fB4OBVW2hDAPNky/xGQUA12PcKgCAmwH8A4C73XBjUlo0CsQmg8FA2+325EXvzR4ZDoexA4L+BUp1HEgldojTV9HaKpVRsBloFIgtkmr5w+Fw8kBH1eZoAEgUSfqKay3koS0aBUISMJlWOBgMYpv57CoiUSTpq2htmRqFRru5qNJm2sQ+Ji4YPvvZz8amsbu7i/X19dBz1FezSdLXLNrKFRPLUaZgq6XAAUJi0lKIcnoWDME+X+qLJOkrq7ayAnYfxVNn743EjMFgkLhiNI1bAv/11FeziZph5NeXqTsVW5UKU6PQqO4jf3M+bKUhwEVHTWE0GuHEiRPjmpGLiODKK68EMF6NKiI4ffp0qnQ9j5fUV3MZjUZYW1ub0o7jONjc3AQQvdo5jqK6kw7k/g0lwbtRQbe1QbjoqBmEuTBWVXziE5/AiRMnYnXSarUyr0ilvupPmLYAYHFxEQCM3kNRpDUkWWhMSyHqRvkp82baxC5RNfbTp08n6mRvb29qNasJ1FcziNLWzs4OrrnmmswGAYhe5WyTxhiFuGZ7JTbTJlaZtcae5sGmvppFnLbSdkcGUdXcZ7E1xihE3ah+v4+9vT1sb2/zgW0QYb5r8oD6ah55ayvvcYXGGAUTB1akOayurk4G/fKC+momeWuryl5SS4V3o/r9PpvzDcebhXbFFVeg3W7n8h3tdpv6aiBFaKsWXlJtBhvrFOjArLmELSqzHby1D9RWs8iqLdNFbJhxrQK4eC0crjRtNmkXDAHjLRRbrVbkecdxJuna3D6RVIss2gKQ6AOJXlJzNgpcadps0tTK/KHX68UaBa5iJlm11e/3dWVlJfTceeedZ61SYWoUGjOm4GHiBI3Ul6z9sbu7uxgMBmi1ph+Z06dPY21tjauYG05WbZ08eRI333wzVlZWps794Ac/wNraWqHOFBtnFKJunKrSk2UDyDpdcGlpCR/4wAdw7tw59Pv9qfO7u7uRA4utVoueUhvALNoCgJtvvjlSW4V6SzVpTpQp5DGm4A/sA642JpMI4rp6TDQR102QNNBIfVUbky0489CWfze2rIBjCtEk3Tj2AVeTNJMI4rZDBMZbJ4a5w07Sjf+lEbVrG/VVTUz1laStoF781+WpGRoFA/K0yqR4otxchz1QSbU5x3H2GYUsLUzqqz4Mh8PIGWhBfSVpq9PppNKXrdYljUIIwaZfmpcIKTdJtbMgcV1ABw4cmHoo4/ZVaLfburKyMtWtwNlI9WA4HOrCwkLk/Q8a+Tht9Xq9qbSS9DUYDKz8jtIYBQBtAHcBuCnk3E8DuBPADwH8skl6WY1ClCWOaq61221rN4PkT1ztrN1uT8VPu3lOlrC4uKidTif0nDeNlZSfpJq/aUvBcZxMuiu6pVDE7KNrAPxVxLmTAN4M4GN5ZyLKdfa5c+fgOE7o8Q9+8IO4+uqrY9PlPrzlIM7P/Llz5/b9PxqN8Oijj07FC5tuOgtnz56FqqLX602dO336NK666qpEvVBf8ydpD4Ogf6ujR4+GxrvkkksyeUnd3d3FNddck/q6zJhYjqwBwMUAvgDg1QhpKfji/QFybinENekcx4ltMUTB1dHlYDgcxt5f05pc0qrlrO4x4tKN60qivuZPkrYcx5m6JkpfUe8YU33VYo9mAJ8CcCmAy+ZtFLIuQR/bzXRpss+4WOLubafTmQzkeX3+WTTgTyOrjsJC3KAz9TV/TLSlqtb0FWc4Zr3vczcKAN4A4APu55mMAoA1AFsAtpaWljIViOk0sTQWmrNLykHcg2gyc8gLUQ9kq9WamiliyyiE1TSTfhf1VRxJ2lI1d4SX9MJP0tas970MRuHdAB4CsA3g7wHsAhhGxI01Cv4wy+yjKP81SYM/UU121uSKIWnBUNJ9MKndd7tdHQwGxovPbA1U+2ubQaiv/JlVW3Fx0uorSVuVbyns+5ISdB+paqg3Qu+hTDvDQJV9vkVgUsZJceJqe1EL1OIWnzmOo4PBwKirwJuTnlZbpr+dZMeGtlSz6Ssqvje9OQ9vu6U1CgCuB3C5+/klGLcmHgNwGsADSWnZnpLqPZRJVjyq6ca9GfLFtLYcdx+y1rjjHvaoqaZRBiSLtpJ+F5kNG9pKk46fNGMPImJlenypjILNkNdAs9e8o2uCcmGjXz2NewL/w2+ri4jaKie2xmxMWxyzaMuGRmgUAphYZq8WwCZ7ebDVr27iyCx4301bA6YPNbVVLmyO2cTpK0pbcaukZzVUYdAoBDAZDPIKnk328lDUizRKH1mnGFJb5Wfe2opbH8WWQgFGIWkRiq2CJ/Yp4kWaxh32wsJC6lYEtVVO5qktETGazkqHeDkZBdXw2Uf+Bz1sBgprc80grishTA/BY1HbKVJbJKmbKqiJwWCQi0ZoFCKIGuDxFhGx37damLywTe6djftObdWLKB2l1VdZ7juNQgRJPnKizrXbbT68JSPsYQvr2jF9AGetxWfVVqvVorZKRtSLPGx6sYm+ytBCpFGIYJYBRdbqyoXJ5IEi+/Rn0VbcymZSPFH3ssrTik2NQhGus0tF2ObaIjK2kAkUvoE2iSXJpbGfkydPJsaZ1U31LNp68sknqa0SEaWtoBt2jyR9VcoFuonlKFOwuUez15SDYW0ToDOyMmE6nQ8GNTlb/b7UVj1Io60kfXFMoQJGIUjZuiFIMmk8lXrdN96DHtana3Mhk0m61FZ5SdJW2BTSxcXFqX2XPfLSVlpoFFIQN1U1zrqXYfCoiZjM7fYMQFR/vqlTs1lr78Ph0GhNQ3BMgdqaDybachwn0uNymL7K4gKdRiEFSbOOwmqXZWkSNhHT2neaxYp51uaS/NwE92umtuaHLW35tcOWQgWNgslGGkHKcqPrQpqasW3XE9735/UijstDGNSWfUz1ZUtbfn2VxcjTKKQgriYX9SCWpUlYB9I+NGkHcE3vbR5dNnH901H7f1NbdkmjL1vaCuqrDN2BNAqGJPX5Rj2IrM3ZI21ZmvT7JoUiamom+QyD2rJLmvK0oa2i9JUWGgVDkmoGUcIJa12UUQhVIEvN2HRyQLBm7t3TIu5TFm15vy2PnbeaSlp9RT3fcRVHx3EiZx+VBRoFQ5K8Ywb9nXgiCMZttVpWdkdqIkk1ubCmdxWmeppqy/8bo/S1uLhYyhdNFchTX1VqvdEoGBJ18/2+jkyblKzNZSOuzzfKv5GpQZjnPTHRVtTvL9tvqTJ56atq94NGwRCTQag0Tckq1RzKRNRAXJaBv1arNXn5evdkHg+vibaGw6HVldkkHJv6KoO2skCjkIK4mQFp+645Q8QuaacILi4uTu5DGWp1SdpK+/uoL7tkmYJaFm2lhUbBAsNh8m5trMnlS5qaXFKzv0z3Jou2yvYb6kBe05vLiKlRyN1Lqoi0ReQuEbkp5NzTROTjIvKgiNwuIst55ycN6+vrY8tpiIjg6NGjOeaoeYR5Hg2j1WrhiSeeiI1j4im1KNJqCwAWFhawsbGRU46aiam+kiiTtmalCNfZ1wD4q4hzbwXwiKr+BIDfAfCeAvJjTBrXzMC41XXixIlyu8WtGKurq9jc3ES/34eIwHEctFrTst3b20tMa2lpKY8sZiKttgCkNiIkmaC+er1epnTKpK1ZydUoiMjFAH4ewIciorwRwAn386cArIiI5JmnNLTb7dTXcM+F2Qn6ngeA7e1t7O3t4dSpU7jwwgtTp9ntdktVy86iLe65YIc4fR0+fDh1emXT1syY9DFlDRi/6C8FcBmAm0LO3w/gYt//XwdwOCTeGoAtAFtLS0s59LZNk8Y1c1iowsBTGQmbsSMi+9aAmPbFe/HKOEMkLt/e7KmoUMbfUxWGw+HU2NPCwkKiR9MqaSsKzHugGcAbAHzA/TyTUfCHIgaaTeaN9/v92EGqqsxIKBtxZZo0jdBxnLn7lzEhaZA5SVvUV3aippc7jqOq1ddWHGUwCu8G8BCAbQB/D2AXwDAQ53MAXu5+PgDgFACJS7cIo2AyI8GrvSY93CQdceXpPbhJtb2yk6QvE21RX9mIK09Vs5ZqVTE1CrmNKajqO1X1YlVdBvAmAF9U1WOBaDcCuNL9/MtuHM0rT6aYzCQ4dOgQTpw4ETv4V6cZCUURN2B3+vRpAMCtt946NdOoBLIxJkkXJtoySYdkIzisqdqsCSRFzD7ah4hcLyKXu//+PgBHRB4EcC2A3yw6P2EkzSQQETz++OPY3d2NjddqtRojJFskDdiNRiPccMMNU8erNAh76NChyHOm2gKoL9uMRiOsra3hsccemzrXqAkkJs2JMoWixhRsbbbBvt/0eKuSg8Hr140q66qs9k3jNoX6skvcmIJJt16Vwby7j6rM6uqqte6IRtUwLHHDDTdgYWFh37GFhQUcP348tsukKjXnM2fOWEuL+krH8ePHM2kLiG/h1QoTy1GmUJSbC1vL31GDGsY8yOrArAo1Z5vaor7Sk1VbVZrMEAbmPfsor1CUUbC1AxPAWSI2MZ0uXGZsaqsKv7cq1EFbcdAoWCCrF8uq1VzLRtJ+tkmboFSh5px2dy/qyx5x+qqDtqKgUbBA2mZ+t9vVwWBQ+UUu8yTNGoSo+1OV2hz1VTymW51WXVth0ChYIE0rwXEcPqAWSFpx6sdkE5syQ30VS9yswuDLvuraCoNGwQJpanJVrkGUibgyDiOpq6nMUF/FkrZbqMraCoNGwQJpBwSj+ifrIqoiSGsUqkxafVFbs9F0lzQ0CpZIGngK6/eN2hC86s3PIkjTfVQH0njjpbZmI+o5FpFGlB2NgmXSGAbHcSI3Y29CjSQKk9rtcDjUTqezr8w6nc4kbh1ryGm7kaIMJ7WVrK04Z3d11JYfGgXL2JpbXuUpbbOQpnYb9XDWtYZMbc0GtWWGVaMA4Hkhxy4zudZ2mJdRUN3flRTVEjCp6TWRuCl+pjW0Ok4T9KC2skNtmWHbKNwP4DcACIDzAfwugL80udZ2mKdRUJ2tVle3mkca4gb5TGtoSXstVL3ZP2uLoaq/e1by1lZd1obYNgo9AO8D8JeugXgngJbJtbbDvI1C2gVHXmi325UVkw2iyi3N2Itp2VfV+GbVlmcUm0qR2qqyvmwbhQUA7wVwN4AHAbzJ5Lo8wryNQhaX2mn7N+tIVJ9tVJlFzRs3rUlXsdmf1V17mL6orfy0VVV92TYK9wC4HkAHwNMB/G8AnzS51naYt1GIm9YW9r9XWwlbw1D3ga0gYS+qOP8/YWVm6i+oioOuptrqdDqTcgjTF7WVXluqYxcYpkahivqybRSOhBy7wuRa22HeRiHqgQv2Ow4Gg9gHswkDWyYkveTj5ubXrSZnqq2ktQrU1hhTbamypeAPpkZhKSyYXGs7zNsoqJo1zZMezKiugirWQLJiunCr3+83os9X1bzbJ05f1FY6balyTMEfTI3CfQDudf/+LYAfAnjA5FrboQxGwYS4/uG4pm0VayAmBF92YS2prKHVammv12tE/7lHkssGaiud+5C48lxZWanF+IxVozB1EfBiAB/Kcu2soSpGIa7m0el09MCBA6HHqyq4OKJWkqZ58SfFr2rtLStx4w8rKyuh57yVu3XChra63W5sV1NdtJWrUVC39ZBw/jwAX8F4kPoBAO8KidMH8AWMWyF/BuDipO+tilHIMue8rtMKZ5lqmeZhr2tNOIy42m2TXKzY0pbjOLHPax3KztQoyDhuPCJyre/fFsYtBUdVXxtzjQDoqepZEekA+AsA16jqbb44nwRwk6qeEJFXA7hKVa+Iy8uRI0d0a2srMc9lYDQa4dixY8bxRQR7e3s55mg+tFotmOjMhH6/j52dndBzdS2/KMaPWLr4dSsfW9oSEXz0ox+NfF7rUHYicoeqHkmK1zJM76AvPA3A/wHwxrgLXON01v2344bg3XsegC+6n29JSrNqrK6uot/vG8dfWlqKPDcajbC8vIxWq4Xl5WWMRiMbWSyEuN+Vhn6/j+3t7cgytfU9VSGqHNrtdujxOurL1j1fWlqKfV4bpS2T5kTWAKCN8YK3swDeE3L+Yxi3HgDgX2FsNJyQeGsAtgBsLS0tWWxQ5U9YN1Kn05nacjKu37Lq885tOHxLmj5YpfKwRdwU1jTlU+XypLbMgY0xBQB/DODGqGDyBW46F2DcEnhB4Pg/BfAZAHcBOA7gIQAXxKVVlTEF1XAnZ/555qYzGuowmyT4e9NsWh9WPk1asRtFcBabf8vOJumL2jLDllH4mbhg8gW+tK4D8I6Y84sAHkpKpypGIaqFkMVxWx3nnZvW8Nrt9ryzWkqor2iorXBsGYXMC9QAXOTV+jH2rPolAG8IxDkM17EegA0A1yelWxWjYFJbMW2WVr0mF8VgMDCaPkimob7iMXWJ0iRsGYU7fZ8/bZKgL/4LMe4Wuhdjz6rXucevB3C5+/mXMV4M9zUAHwLwtKR0y2oU/E3OtM1Xk7Tr2M9pUk5VfjHZhPpKT1I5NU1btozCXWGf5xnKaBRs+MFP6sOsWz+nqRsC7+Gt+u+dhVn1ZaKbJuqrDntwpCGPlsKdJgnmHcpoFGZdQBNcbFSHWloSacusCWUSha0FWk0qR+prGlOjELt4TUTOAXgMT+24tuudAqCq+iORF+dEGRev2Vyc5eHNya8rWcqs7mUShW19NaEcqa9prCxeU9W2qv6Iqh5U1QPuZ+//wg1CWcljYcvJkyetp1kmspRZ3cskCtv6akI5Ul/ZMV3RTGLY2NhAt9u1mmbdV1BubGykvqbuZRJFlrKKownlSH1lh0bBAqurq9jc3Ix0L5CWbrdr/UVQNlZXV+E4TuT5AwcO7Pt/YWEBZ8+erZwbBhsklVUamqAtYFxmvV4v8vzCwsK+/0UEOzs7jdNWKCYDD2UKZRxo9rCx5L7dbtfSxXEYUeW1srIyNQWz0+k0bmDQD7WVnuFwOKUbYOxC3O9tILhWpq7aQt6us+cVymwUVNNNtYwKdRVlGDZ2sWsKcVttUlvhJOmrSdqiUZgjUf7s04Q4UdZtTnkSdXPDMCuz6ovaeoomaYtGYY7MahDiRFn21adpXypsKaSnqdpSTacvams/NApzxKSJH3SdHQxRu7CVWcRpXyqm8avwsiqSJH2FbfXqD+12u5IvyDQ6oLamoVGYI0mDgmEDp1HxgrWcMjd307xUhsNhZDdImEFsWrdGHHFbcfZ6PSNteddXRVuq5vqitsKhUZgzUV4aRUQXFxcTH9pg8GovZa7Nmb5UTGbSpN1zommEeZits7ZUzfRlqi0vbpP0RaNQEkzdQ5sET7hB0YtIKaYaRr1Ugt0VJt1rYRup17VZn5UmaUvVTF8m2or6rXXXF41CSbDpzMyrEYW9DMogaJNus7iuD9OXFxnTJG2pxuvLM14m2vJaBk3TF41CSbBVk/MLtsyCjuvP9V4wafYDiHp5keZpSzVeX6Z7TXhdRk3Tl6lRoJuLnInyp+I4Dvr9PkQEjuNgcXExNh2/e4Iox11lcOi1urqKvb29yPO7u2NHu3G+orrdbqRbB/qneQpTbSW5yKiKtoB4fY3fe8na2tjYiCw76gtsKeRNXN+lf3AvrnZt2idfltqcSbfGYDCY1Ngcx5maadXEPt+0JJWRyQroqmlLNVlffj2FaUu1WVNRPcDuo/IQNsvB1JdNFedWm/y2sJdX2IPbpNkhWYgruzT3wJ9embWlGj8lN6riFaahpumLRqHkpJklEUbZBR01JTfs95X9JVRF0rYQ/JRdW6pmM684g20/pkYhdue1MlLGndeykLQzlIjE9s1XhdFohGPHjoWeExEsLS1hZ2dn6ly/5rtg5U0T9DUajbC+vh6qnziaqi0rO6/NmIHzROQrInKPiDwgIu8KibMkIreIyF0icq+IHM0rP2UjaUAr7YDXaDTC8vJy4fsNJH1v3F4AS0tLkQOYaR90sh+b+pqXtpK+e3V1Fdvb26n3mtjZ2eGeCXGYNCeyBIz3cV50P3cA3A7gZYE4mwAG7ufnAdhOSrfq3Udxftz9Ic2CoXl1wZh8b5RP+4WFhdiBUBFpbDN/Fmzra57de7PqK677sondSCjTmAKALoA7Abw0cPx/APgN9/PLAXw5Ka0qGwXTwWUv9Hq9Sb+uf7ZOsJ93XjNGor7X718mKU7coGGZZrxUgbT68sq4jNqK+25TfSWVR9P0VQqjAKAN4G4AZwG8J+T80wHcB+AhAI8AuDQinTUAWwC2lpaWciu0vLG5AtUvfBsLcbIMLsbVRNM4WYtKo84LifLApr7mrS1VO/oaDqM3vWqavkphFCZfAlwA4BYALwgcvxbAf9KnWgpfBdCKS6vKLQWbK1C9ELdC2LQmlLWLIO4l1G63VURivVV6L4qoOE2ryc2KbX3NU1uqyUau3+9H5s/TX1ycpumrVEZhnB9cB+AdgWMPAHiG7/9vAPixuHSqbBRstxT8L9hZ+n2zdhHE1cKSQtLuYU3s852VPPQ1L22pmumr0+kkaqnT6UztX9JEfc3dKAC4CMAF7ufzAXwJwBsCcf4EwJvdz/8cwP8DxtNko0KVjUKWPl+T4A3KZp1bPsuA5Cx+jMIMRZnnxpedPPSVp7a8Lqo4bOnL3zJtqr7KYBReCOAuAPcCuB/Ade7x6wFc7n5+HoBbAdyD8djDa5LSrbJRUJ3uX7Uh+FmbwUn5iDMMNl9ETevjzQO/vvLeK9yEJG15s9Difo8NfVFbJTAKeYWqG4UgNmpCs/q7T3rw2u124vVlMG5kP7O6KQfMavNJeUh6qefZTUltPQWNQkUYDoeJ+zUX9fDGpZ/ErMatiX28RWBjI55ZN9pJ0lZSLX5W40ZtjaFRqBD+Jn+W7RT9D5dXK8ryEER1N5i0FOI2i+/1elO1RW9xUZP7eIvCxFuqacUj6/2K+/64WnxSpWlhYYHaMoRGoYLYeniz1o4Gg0FoWkm1RJM8R7kwJsVh4qQwL30lrWyPgtqyB41CxbC5365JDSwuH16Lod1uG3UbmOabzfj5UQZ9BY2SSZcntWUPGoUKkKVl0Gq1Uo1BFPGgpKl9csCvOLLqK038vKG27GFqFLgd55wYjUZYW1tL7Q10b28PBw8eNPYMuba2ViqPkGXZ1rHuzKKvXq8HEUmMKyLUVg2hUZgT6+vrk/2K03LmzBmcOnUKw+EQ/X4/Nu7u7i7W19czfU+a/Jhy6NChHHNCPGbR1+7uLj760Y8maktVS6WtVqtVKiNVVWgU5sQstRrPF77nTz7p4c27BpXGN//3vvc9PrgFMKu+yq6tsJbMuXPnStcyriI0CnMi7SY6Hp1OBxsbG5P/r776anzrW9/K5buCRG14srGxgW63uy9ut9tFr9ebSuPJJ5/MvXZJ7OirSG0B4fqK0tbb3/52tNvtqTSKaBnXHpOBhzKFugw0Z1m+H5ytETWF1B9szcpI8nYZNmvJxK1x0/3R5MWs+ipSW1H59dKPmrX7UEF3AAAQJElEQVQUpy9qaxpw9lH5CRPuyspKqNDDpobGzRSx/TDEebuMeqDjXBbPc0evphDU12AwiFwcGdRXkdpStauvWT271hUahQpjulYgrhZn+wGIq5VFPdBhK5m9h3MWl8pkNkz0VaS2VLPpK+rlz/0TwqFRqBFhLYokfzJZakb+l7X30vC+L+tK2MXFxdDVpjZ29CJ2CNNXXtoSkakVyH4jlTaErWamtsKhUagJaZvOwQfGtF81rg96YWEh9aKmpJcIWwrlIEpfSffU05OJvvLaRyRKX9RWODQKFSLu4bLlCymphmfze0weSI4pFEde+grWyKPuX97aCuqL2gqHRqEiJAk4D381YS+JPPaPDr5Awn47Z4jkS9H6Cmud5q2tMH1RW9OYGgUZx60OR44c0a2trXlnwxrLy8uhrgj6/T62t7cjzy8uLuKxxx5DlvvX7Xb3rXbtdrs4//zzcfr06dRpmeL9HlIsWfXVarWwt7c38/cXoS2A+jJBRO5Q1SNJ8bh4bc5ErQj1joct3llYWMDjjz+eySAAmHJ/sLu7i0ceeSRTWiZ0u919C+5IcWTRFwArBgHIX1sA9WUbGoU5E7Ui1O/KYnNzE/1+HyKCfr+PgwcP4sknn7Saj+BLoNWyI41+v4/NzU2srq5aSY+kI62+wlYJz0pQW71ez9ihYxLUVw6Y9DGVKTRtTCGMIvpo4Y4/zLIpS9Nne5SBtPoqg7ZMZ7pRX+nAvAeaAZwH4CsA7gHwAIB3hcT5HQB3u+FrAL6blG7djIJq+kGxWWZzLC4uppoeuLCwELpjVtKDzNke5SGNvmadjZSHtnq93lQ86is9ZTAKAmDR/dwBcDuAl8XE/zUAH05Kt45GIS2zzPt2HEcHg0Gqh99bIJQUh7M9qk+YttK0HrJoKyl+mjURJJq5G4V9XwJ0AdwJ4KUxcb4M4OeS0qJRGBNcIWpS4/LXslZWVlI97ElGKG61KB/oahHmM8m0EpJWW57zurj0qS07lMIoAGhj3DV0FsB7YuL1AXwbQDspTRqFcPwPR5rtOm2FqP5dLiSqB1m29jQJvV6P2iqIUhiFyZcAFwC4BcALIs7/BoDfjbl+DcAWgK2lpaVcCqwu5O1SIKp2mHalNAcJqwe1VW1MjUJhi9dE5DoAu6r6X0PO3QXgV1X1y0np1G3xmm2iFiPlRb/fx8bGRuSUwFarhTCNiYi1ufCkGIrWVrvdxokTJ6gtS8x98ZqIXCQiF7ifzwfwcwD+OiTecwFcCOAv88pLE/B2rSryoe31ejh58iTW19cjt0BMmidPqsFoNCpUW8B4e01qaw6YNCeyBAAvBHAXgHsB3A/gOvf49QAu98X7LQC/bZouxxSmGQ6HM3kxtRGimvns960+JruwUVvlB2UaU7AZaBSmSRpY9m9rmGefcNyAIGeIVBPTvRWorfJDo1BzTGeDBPd1Hg6HmTc0MQmkHpjqi9qqDqZG4QBI5RiNRlhbW5tybBfG97///X3/e4N2b3nLW/DEE09YzZctf0lkvqTRlx9qqx7QdXYFSTugHOZW+PDhw7m4M66ansg0afRFbVWHuc8+IvkR5Q45ip2dHSwvL6PVamF5eRmj0QhnzpxJ/b2O46Df76e+jlSLNPqypS0A1jynktmgUagghw4dShVfRLCzswNVxc7ODq644orUta5ut4vjx49je3s78uHlQ10P0ugrqK2rrroqU43e0xe1NX9oFCrGaDTCo48+GnleRKaOBR/StA+t4zj7fNYfP34cnU5nX5xOp4Pjx4+nSpeUjyR9BQlqKcs+H359UVslwGQ0ukyh6bOPomaEOI6jquMZILPsgRCXth9OBawncfryz0iyue9CcLoptZUPKJubC1s0faDZZGl/Hiubh8Mhd7dqAPPSV9XeQ1WEA801xWRpf9qBaBPW19etp0nKxzz0lccWoCQ7NAoVI2yj9eDG5Xn4fsnD0JDyMQ99nTt3zmp6ZDZoFCpGcKP1sI3Lwx7sWUk744lUk3noizOLygXHFGrKaDTC+vo6Tp48iUOHDs28mMhxHJw6dcpS7kjVsakvaqsYTMcUaBQagI2BQfqoJ1HMqi9qqxg40Ewm2JgpQh/1JIpZ9UVtlQsahRrjbbyTluBskOBAIyGetmZ1VEdtlQ8ahZriebo0qcV5RqDdbmMwGODEiROxA42k2fi1lab72XEcDAYDaqvkcEyhppj284Z5uSQkjlm9qJL5wDGFhmOyrmBhYYFNd5Ia0zUr1Fc1oVGoKUmDd47j4MMf/jCb7iQ1Udryjy9QX9WFRqGmRK1MHQ6HUFWcOnWKDyzJRJS2PvKRj0ycqlFf1SU3oyAi54nIV0TkHhF5QETeFRHv34jIV904H8srP03DZGUqIVmgtupNbgPNMnbs31PVsyLSAfAXAK5R1dt8cZ4N4BMAXq2qj4jIj6nqd+LS5UAzIYSkx3Sg+UBeGXD9d591/+24IWiB3gbg/ar6iHtNrEEghBCSL7mOKYhIW0TuBvAdAJ9X1dsDUZ4D4DkicquI3CYir8szP4QQQuLJ1Sio6jlVvQTAxQB+SkReEIhyAMCzAVwG4FcA/J6IXBBMR0TWRGRLRLYefvjhPLNca/yrUL1N1gmxAbVVHwqZfaSq3wVwC4BgS+AhADeq6pOq+k0AX8PYSASv31TVI6p65KKLLso/wzUkuAp1Z2cHa2trfHjJzFBb9SLPgeaLADypqt8VkfMB/CmA96jqTb44rwPwK6p6pYgcBnAXgEtUNdIPLweasxG1CpUrTsmsUFvVYO4DzQCeDuCEiLQxbpF8QlVvEpHrMd5A+kYAnwPwGhH5KoBzAH49ziCQ7EStQuWOamRWqK16kefso3sBvCjk+HW+zwrgWjeQHFlaWgqtzdFtMZkVaqtecEVzQzDZe5eQLFBb9YJGoSFwFSrJC2qrXtB1NiGENAC6ziaEEJIaGgVCCCETaBQIIYRMoFEghBAygUaBEELIBBoFQgghE2gUCCGETKBRIIQQMoFGgRBCyITKrWgWkYcBTHvfSsdhAKcsZMcmZcwTwHylhflKB/Nlzqx56qtq4oY0lTMKNhCRLZPl3kVSxjwBzFdamK90MF/mFJUndh8RQgiZQKNACCFkQlONwua8MxBCGfMEMF9pYb7SwXyZU0ieGjmmQAghJJymthQIIYSEUCujICKvE5G/EZEHReQ3Q84/TUQ+7p6/XUSWfefe6R7/GxF5bcH5ulZEvioi94rIF0Sk7zt3TkTudsONBefrzSLysO/7/53v3JUi8rduuLLgfP2OL09fE5Hv+s7lUl4i8mER+Y6I3B9xXkTkv7t5vldEXuw7l2dZJeVr1c3PfSLyZRH5Sd+5bff43SJidecqg3xdJiL/6LtX1/nOxd7/HPP067783O9q6ZB7Ls+yeoaI3OK+Ax4QkWtC4hSnL1WtRQDQBvB1AM8CsADgHgDPC8S5GsAN7uc3Afi4+/l5bvynAXimm067wHy9CkDX/Tzw8uX+f3aO5fVmAO8LufYQgG+4fy90P19YVL4C8X8NwIcLKK+fBvBiAPdHnD8K4E8ACICXAbg977IyzNcrvO8D8HovX+7/2wAOz6m8LgNw06z332aeAnF/AcAXCyqrpwN4sfv5IICvhTyLhemrTi2FnwLwoKp+Q1WfAPBHAN4YiPNGACfcz58CsCIi4h7/I1V9XFW/CeBBN71C8qWqt6jqrvvvbQAutvTdM+UrhtcC+LyqnlHVRwB8HsDr5pSvXwHwh5a+OxJV/XMAZ2KivBHAR3TMbQAuEJGnI9+ySsyXqn7Z/V6gOG2ZlFcUs+jSZp4K0RUAqOq3VfVO9/OjAP4KwI8HohWmrzoZhR8H8C3f/w9humAncVT1hwD+EYBjeG2e+fLzVoxrBB7niciWiNwmIr9oKU9p8vVLbnP1UyLyjJTX5pkvuN1szwTwRd/hvMoriah851lWaQlqSwH8qYjcISJrc8jPy0XkHhH5ExF5vnts7uUlIl2MX6yf9h0upKxk3KX9IgC3B04Vpq8Ds1xM7CIixwAcAfAzvsN9Vf07EXkWgC+KyH2q+vWCsvTHAP5QVR8XkX+PcSvr1QV9twlvAvApVT3nOzbP8iotIvIqjI3CK32HX+mW1Y8B+LyI/LVbmy6COzG+V2dF5CiA/wXg2QV9dxK/AOBWVfW3KnIvKxFZxNgQ/UdV/Z7NtNNQp5bC3wF4hu//i91joXFE5ACAHwVw2vDaPPMFEflZAOsALlfVx73jqvp37t9vAPgzjGsRheRLVU/78vIhAJeaXptnvny8CYEmfo7llURUvvMsKyNE5IUY3783qupp77ivrL4D4H/CXpdpIqr6PVU9637+LICOiBxGCcoL8brKpaxEpIOxQRip6mdCohSnrzwGTuYRMG71fAPj7gRvgOr5gTi/iv0DzZ9wPz8f+weavwF7A80m+XoRxoNrzw4cvxDA09zPhwH8LewNupnk6+m+z/8SwG361ODWN938Xeh+PlRUvtx4z8V48E+KKC83zWVED5z+PPYPBH4l77IyzNcSxmNkrwgc7wE46Pv8ZQCvKzBf/8S7dxi/YE+6ZWd0//PIk3v+RzEed+gVVVbu7/4IgP8WE6cwfVkTQBkCxiP0X8P4BbvuHrse49o3AJwH4JPuQ/IVAM/yXbvuXvc3AF5fcL5uBvAPAO52w43u8VcAuM99MO4D8NaC8/VuAA+4338LgOf6rn2LW44PAriqyHy5//8WgN8OXJdbeWFcc/w2gCcx7rd9K4C3A3i7e14AvN/N830AjhRUVkn5+hCAR3za2nKPP8stp3vce7xecL7+g09bt8FntMLufxF5cuO8GeNJJ/7r8i6rV2I8ZnGv7z4dnZe+uKKZEELIhDqNKRBCCJkRGgVCCCETaBQIIYRMoFEghBAygUaBEELIBBoFQkLweVu9X0Q+KSJdEVmO8rBJSF2gUSAknO+r6iWq+gIAT2A8Z5yQ2kOjQEgyXwLwE+7ntoj8nuv3/k9F5HwAEJG3icj/dR28fdp1qgYR+ddua+MeEflz91hbRN7rxr/X9StFSCmgUSAkBtdH1usxXkUKjJ22vV9Vnw/guwB+yT3+GVV9iar+JMauj9/qHr8OwGvd45e7x94K4B9V9SUAXgLgbSLyzPx/DSHJ0CgQEs75InI3gC2M/fL8vnv8m6p6t/v5Dox96QDAC0TkSyJyH4BVjP1pAcCtAP5ARN6G8QYyAPAaAP/WTf92jN23l8VDKGk4dJ1NSDjfV9VL/AfG+zHhcd+hcwDOdz//AYBfVNV7ROTNGO8sBlV9u4i8FGOHZneIyKUY+7H5NVX9XJ4/gJAssKVAiB0OAvi26wJ51TsoIv9MVW9X1esAPIyxm+PPARi4cSEizxGR3jwyTUgQthQIscN/xrgr6GH370H3+HtF5NkYtw6+gLGnzXsx7na6090O9mEARe4SR0gk9JJKCCFkAruPCCGETKBRIIQQMoFGgRBCyAQaBUIIIRNoFAghhEygUSCEEDKBRoEQQsgEGgVCCCET/j+GZJb4vJUKDwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "data=np.loadtxt(os.path.join(data_dir,'cepheid.csv'),delimiter=',',skiprows=1)\n",
    "t=data[:,0]\n",
    "y=data[:,1]\n",
    "def lomb_scargle(t,y):\n",
    "\n",
    "    frequency, power = LombScargle(t, y).autopower()\n",
    "    #plt.plot(frequency, power)\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",
    "\n",
    "def plot_best_period(t,y):\n",
    "    phase_fold(t,y,lomb_scargle(t,y))\n",
    "\n",
    "plot_best_period(t,y)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}