Adding LS periodogram exercise

lessons/07-example.ipynb

@@ -4,10 +4,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Practicing plotting\n",
+    "[Pre-MAP Course Website](http://depts.washington.edu/premap/seminar/seminar-cohort-12/) | [Pre-MAP GitHub](https://github.com/UWPreMAP/PreMAP2016) | [Google](https://www.google.com)\n",
     "\n",
-    "Let's put together some of the you've worked on. There is a file called `data/cepheid.txt`, which contains flux measurements from a variable star over time. In this notebook, your task will be to reproduce some plots. **Bold text** makes clear what your task is in each cell. \n",
+    "# Practicing plotting: Cepheid light curves\n",
     "\n",
+    "Let's put together some of the skills you've learned! You have a file called `data/cepheid.txt`, which contains flux measurements from a variable star over time. In this notebook, your task will be to _measure the approximate distance to this star_. **Bold text** makes clear what your task is in each cell. \n",
+    "\n",
+    "### Read in and plot the data\n",
     "We have used the function `ascii.read` to read in data before (see `06-plotting.ipynb` for example). **Use `ascii.read` to load the file `data/cepheid.txt` into a variable `data_table`**:"
    ]
   },
@@ -46,6 +49,8 @@
    "source": [
     "This type of star is called a [Cepheid variable star](https://en.wikipedia.org/wiki/Classical_Cepheid_variable#Period-luminosity_relation). These stars are useful to astronomers because their flux varies in a predictable sinusoidal pattern. It turns out that the period of the flux oscillation is directly related to how intrinsically bright the star is – so if you know how long the period is, and you measure how dim the star _appears_ to be, you can estimate how far away it is because you know how bright it really is.\n",
     "\n",
+    "### A very rough \"model\"\n",
+    "\n",
     "In the cell below, roughly **measure the period of this Cepheid** by plotting a sinusoidal function over the data. The function that you should play with is: \n",
     "\n",
     "$$\\textrm{flux} = \\textrm{amplitude} \\times \\sin \\left( \\frac{2\\pi}{\\textrm{period}} \\left( \\textrm{time} - \\textrm{offset} \\right) \\right)$$\n",
@@ -69,6 +74,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "### Period-luminosity relation\n",
+    "\n",
     "The brightness of this Cepheid variable star is related to its period by the _period-luminosity relation_: \n",
     "\n",
     "$$ M_{v}=-2.43 \\left(\\log _{{10}}( \\textrm{period})-1\\right) - 4.05 $$\n",
@@ -111,6 +118,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "### Measure the distance\n",
+    "\n",
     "To get the distance to the star, we need to use one more equation, which calculates the distance $d$ in parsecs to a star given its absolute magnitude $M$ and its apparent magnitude $m$: \n",
     "\n",
     "$$ M = m - 5 (\\log_{10}{d} - 1) $$\n",
@@ -130,6 +139,87 @@
    },
    "outputs": [],
    "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### _Phase-folded_ light curve\n",
+    "\n",
+    "When you find the period of some periodic light curve, we often want to line up each cycle of the observations and line them up with one another, so you can see the pattern in detail. This is called _phase folding_. You can **phase fold the light curve replacing the `times` in your plotting commands with the modulus of the times in the light curve and the period**:\n",
+    "\n",
+    "$$\\textrm{folded times} = \\textrm{times} \\,\\%\\, \\textrm{period}$$\n",
+    "\n",
+    "![](files/data/example_plot_5.png)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Making a periodogram\n",
+    "\n",
+    "Astronomers are often looking for periodic signals in various forms of data. One of the diagnostic plots that astronomers often make when looking for periodicities is the _Lomb-Scargle periodogram_, which shows how periodic a signal is over a range of possible periods, allowing you to pick out the approximate period.\n",
+    "\n",
+    "The function `lomb_scargle_periodogram` takes the times and fluxes from your light curve, and returns the strongest period in your light curve and makes a plot of the periodogram. **Run this function on both the red and the blue data, and identify the period of the light curve. How close was your guess from earlier?**"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "from gatspy.periodic import LombScargle\n",
+    "import warnings\n",
+    "\n",
+    "def lomb_scargle_periodogram(times, fluxes):\n",
+    "    \"\"\"\n",
+    "    Calculate the best period for a light curve using\n",
+    "    a Lomb-Scargle periodogram, and plot the periodogram.\n",
+    "    \"\"\"\n",
+    "    # Use gatspy to calculate the LS-periodogram \n",
+    "    model = LombScargle(fit_period=True)\n",
+    "    model.optimizer.period_range = (20, 30)\n",
+    "    model.optimizer.quiet = True\n",
+    "    model.fit(times, fluxes)\n",
+    "    \n",
+    "    periods = np.linspace(10, 40, 10000)\n",
+    "    \n",
+    "    with warnings.catch_warnings():\n",
+    "        warnings.simplefilter(\"ignore\")\n",
+    "        scores = model.score(periods)\n",
+    "        \n",
+    "    # Plot the periodogram\n",
+    "    plt.figure()\n",
+    "    plt.plot(periods, scores)\n",
+    "    plt.xlabel('Period [days]')\n",
+    "    plt.ylabel('Power')\n",
+    "    plt.title(\"Best period: {0:.4f} days\".format(model.best_period))\n",
+    "    \n",
+    "    # Return the best period\n",
+    "    return model.best_period"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {