docs/source/tutorials/3-science-examples/periodograms-verifying-the-location-of-a-signal.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "H4cDVgIW1Fgk"
   },
   "source": [
    "# Verifying the location of a periodic signal in the pixel data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PkncygSd1NlF"
   },
   "source": [
    "## Learning Goals\n",
    "\n",
    "By the end of this tutorial, you will:\n",
    "\n",
    "* Understand the causes of signal blending and aperture contamination in _Kepler_ and _K2_ data.\n",
    "* Be able to use Lightkurve's [plot_pixels()](https://docs.lightkurve.org/reference/api/lightkurve.KeplerTargetPixelFile.plot_pixels.html?highlight=plot_pixels) function to visually identify the pixel source of a signal.\n",
    "* Be able to implement difference imaging to find the pixel source of a sinusoidal signal.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OIWLOsJgDVan"
   },
   "source": [
    "## Introduction\n",
    "\n",
    "This tutorial is part of a series on handling _Kepler_ and _K2_ data with Astropy and Lightkurve. To work through this tutorial, you should be familiar with downloading and handling both **light curves** and **target pixel files** with Lightkurve, and you should have experience working with **periodograms**. We'll use light curves and periodograms to detect signal, and follow that up with detailed analysis on the pixel scale to pinpoint the signal's source.\n",
    "\n",
    "Some useful terms to keep in mind when working with signal verification are _contamination_ and _blending_. These terms are often used interchangeably. Here, we'll use blending to refer to any scenario where flux from two or more targets become visible in one target's aperture. We use contamination to refer, more specifically, to the signal that erroneously enters the aperture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "f0yfyjrD2ERM"
   },
   "source": [
    "## Imports\n",
    "\n",
    "We'll use [Lightkurve](https://docs.lightkurve.org/) for downloading and handling _Kepler_ data throughout this tutorial. We'll also use [NumPy](https://numpy.org/) to handle arrays for aperture masks, and [Matplotlib](https://matplotlib.org) for visualizing data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "xOK55DOl3M-2"
   },
   "outputs": [],
   "source": [
    "import lightkurve as lk\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "q7MV3nlB2oxJ"
   },
   "source": [
    "## 1. Background\n",
    "\n",
    "The _Kepler_ space telescope observed with 4x4 arcsecond square pixels. This means that if two stars are within the same _Kepler_ pixel, it is difficult to identify which is the source of an observed signal. Fortunately, these coincidences are unlikely. But the _Kepler_ mission observed areas close to the Galactic plane, and both _Kepler_ and _K2_ observed stellar clusters, all of which are visually crowded regions. Due to crowding, there are still many cases where the true source of a signal is unclear.\n",
    "\n",
    "The process of signal verification begins with a light curve: first, we must detect our signal. In the case of potential exoplanet transits, these will be evident in the light curve; for other phenomena, including binarity and stellar variability, we might need a periodogram to pick up on the signal. To make sure our transits really do belong to a planet, or if there's any uncertainty about the source of a stellar signal (such as, is the star in a crowded region? Is there a bright star nearby? Are we seeing multiple signals?), we need to look at the target pixel files, which will give us a more complete picture, and help us identify nearby contaminating sources."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "YvBqnyzC26v2"
   },
   "source": [
    "## 2. Pixel Level Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BiXxp-vh6NSf"
   },
   "source": [
    "The most basic method of signal verification is to look at the pixels. In this section, we'll use tools included in Lightkurve to examine light curves and periodograms in individual pixels and help us identify the source of a signal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "OQysyDeqdKEO"
   },
   "source": [
    "### 2.1 Initial analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ytsilBSp2liS"
   },
   "source": [
    "KIC 2435971 is a star of [_Kepler_ magnitude](https://keplerscience.arc.nasa.gov/the-kepler-space-telescope.html#flux-calibration) Kp = 14.7, with a massive signal. It's not impossible that this signal comes from the star, but KIC 2435971 is in the field of open cluster NGC 6791, which means it's subject to a  higher level of crowding than much of the _Kepler_ field. This could be a case of contamination; to be absolutely certain, we're going to use this star for our signal verification exercise. Let's begin by downloading a light curve; then we can produce a Lomb-Scargle Periodogram to search for repeating signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Ygpu69Io7aWu"
   },
   "outputs": [],
   "source": [
    "lcs = lk.search_lightcurve('KIC 2435971', author='Kepler', cadence=\"long\").download_all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 401
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1893,
     "status": "ok",
     "timestamp": 1598467661146,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "XZX_vZlD7lxa",
    "outputId": "b8bc444f-c3d2-44ac-a380-44b0753612a3"
   },
   "outputs": [],
   "source": [
    "lc = lcs.stitch().remove_outliers()\n",
    "pg = lc.to_periodogram()\n",
    "pg.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VGIFXMcJ3lpi"
   },
   "source": [
    "These high-amplitude, narrow peaks in the periodogram suggests a compact binary star."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "RzBQ2jG37qd4"
   },
   "source": [
    "### 2.2 Pixel-level analysis\n",
    "\n",
    "Often, the identification of a much brighter target in the frame is enough to confirm contamination. Let's begin by looking at the target pixel file (TPF) we're going to be working with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "uFh5bS-k4Hv4"
   },
   "outputs": [],
   "source": [
    "tpf = lk.search_targetpixelfile('KIC 2435971', author='Kepler', quarter=9).download()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 420
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1141,
     "status": "ok",
     "timestamp": 1598467687980,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "t1yBDMR74aQa",
    "outputId": "4eeb34c0-252f-40f9-e732-ff763371b031"
   },
   "outputs": [],
   "source": [
    "tpf.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DYA35lvHJ79e"
   },
   "source": [
    "As you can see, there's another star very close by, at the top of the frame. This star is brighter than our target in the center, and could conceivably produce a high-amplitude periodic signal. But it's not a huge difference in flux, so we should double-check before we go ahead and claim that this is our contaminant.\n",
    "\n",
    "Luckily, there is a utility in Lightkurve to deal with this scenario. We can look at the light curves — and periodograms — in each individual pixel, using the `plot_pixels()` function.\n",
    "\n",
    "Let's have a look at what the function does at its most basic:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 435
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 2563,
     "status": "ok",
     "timestamp": 1598467725866,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "jTSket5bDpxT",
    "outputId": "bd15f763-c7db-44b7-ffd9-faf83a7b8619"
   },
   "outputs": [],
   "source": [
    "tpf.plot_pixels();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "56bUogh_Dump"
   },
   "source": [
    "These are the light curves in every pixel of this TPF, for Quarter 9. We can't tell much just by looking at this, but we can give [plot_pixels()](https://docs.lightkurve.org/reference/api/lightkurve.KeplerTargetPixelFile.plot_pixels.html?highlight=plot_pixels) a corrector function to flatten these light curves. To create the corrector function, we use a Python construction called a [lambda function](https://docs.python.org/3/howto/functional.html#small-functions-and-the-lambda-expression). This kind of function is known as an \"anonymous\" function, as it has no name. Lambda functions are just one line of code, so they're very useful for when we want a function to take another function as an argument. Here, we take a variable `x` and, assuming it's a `LightCurve` object, we can apply methods like [.remove_outliers()](https://docs.lightkurve.org/reference/api/lightkurve.LightCurve.remove_outliers.html?highlight=remove_outliers) as normal.\n",
    "\n",
    "We can also overlay the TPF colors to remind us of the stars we're looking at:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 435
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 3200,
     "status": "ok",
     "timestamp": 1598467768547,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "MlM_ch_zD9Zu",
    "outputId": "ddb62184-63f9-4c17-b3f3-ab80d6738bca"
   },
   "outputs": [],
   "source": [
    "tpf.plot_pixels(corrector_func=lambda x: x.remove_nans().flatten().remove_outliers(), show_flux=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "guJKCifqEZvB"
   },
   "source": [
    "Looking at it like this, we can start to see some evidence of the signal from our periodogram above, but not enough to make a call. Let's go back to `plot_pixels()`: we can tell the function to compute a periodogram for each pixel, which should make it clear where the signal is strongest. This time, we're also going to overlay the pipeline aperture, which is the selection of pixels that went into producing the light curve for this quarter of data. The pipeline aperture is shown by red boxes around all included pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 435
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 3407,
     "status": "ok",
     "timestamp": 1598467801926,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "n50XGaFPfpsT",
    "outputId": "e4ea36b8-b83a-4d37-a9cd-329af93e693f"
   },
   "outputs": [],
   "source": [
    "tpf.plot_pixels(corrector_func=lambda x: x.remove_nans().flatten().remove_outliers(), periodogram=True, show_flux=True, aperture_mask='pipeline');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "FePm_BlLE0ZT"
   },
   "source": [
    "We can see the two high-amplitude peaks in most of the pixels in this TPF. The key feature is the ratio of that signal's amplitude to the noise level: pixels with higher signal-to-noise ratios contain more flux from the source of the signal. This is often something we can check by eye; here, we can confirm that the signal is coming from our target star, KIC 2435971. The bright star at the top of the frame, in the yellow pixel, shows a significantly higher noise level than the green pixels below it, and the placement of the pipeline aperture reflects this.\n",
    "\n",
    "In fact, KIC 2435971 is a [known eclipsing binary](http://simbad.u-strasbg.fr/simbad/sim-basic?Ident=KIC+2435971). This example challenges the assumption that any high-amplitude signal usually comes from the brightest star in the frame, but nevertheless shows that it's always useful to verify our signal, especially in crowded fields. If we wanted to study the bright star, we could perform **custom aperture photometry** (see the dedicated tutorial on this topic), though we wouldn't be able to exclude all of this signal, as it has such a high amplitude and contaminates so many of these pixels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tFa1QpKc8AZ4"
   },
   "source": [
    "### 2.3 Identifying a contaminant\n",
    "\n",
    "For this exercise, we're going to switch to a different star. KIC 7024511 was flagged as a _Kepler_ Object of Interest (KOI) in Quarter 6 of the mission, and given the designation KOI 311. This means that it's a potential exoplanet host, based on transits detected. Let's have a look at one quarter of data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 387
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 4765,
     "status": "ok",
     "timestamp": 1598467946369,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "pM6T2qGYnCVW",
    "outputId": "4da47c77-594f-40f6-c887-95cdc09edc53"
   },
   "outputs": [],
   "source": [
    "lc_koi = lk.search_lightcurve('KIC 7024511', author='Kepler', quarter=11).download(quality_bitmask='hard').flatten()\n",
    "lc_koi.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "tzI4C8l_oEWy"
   },
   "source": [
    "It's hard to pick out the transits by eye, but you should be able to see one point lower than the others around 1030 BKJD (Barycentric *Kepler* Julian Date). This dip has a depth of less than 0.01%, which is consistent with the signal coming from an exoplanet. Unfortunately, it's also consistent with the diulted signal of a nearby eclipsing binary — dilution here referring to the diffuse flux towards the edges of a star's spread on the detector, where the amplitude of any variable signal can decrease. We can look at the TPF using `plot_pixels()` to confirm whether or not this is the case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 769
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 8222,
     "status": "ok",
     "timestamp": 1598467982200,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "ab2QB-6HlDpi",
    "outputId": "b23f54ca-18d1-455e-ecd0-76e18fba7cec"
   },
   "outputs": [],
   "source": [
    "tpf_koi = lk.search_targetpixelfile('KIC 7024511', author='Kepler', quarter=11).download()\n",
    "tpf_koi.plot_pixels(aperture_mask='pipeline', corrector_func=lambda x: x.remove_nans().flatten());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZA69h-_PmYUJ"
   },
   "source": [
    "In this case, we can clearly see an eclipse to the left of our target. The eclipsing binary is outside of the pipeline aperture, but the signal is strong enough to have contaminated it. And sure enough, if we check on the [NASA Exoplanet Archive](https://exoplanetarchive.ipac.caltech.edu/cgi-bin/DisplayOverview/nph-DisplayOverview?objname=K00311.01&type=KEPLER_CANDIDATE), KOI 311 was designated a false positive in the first data release it was included in.\n",
    "\n",
    "We can double-check this by using Lightkurve's [cone search](https://docs.lightkurve.org/reference/api/lightkurve.search.search_lightcurvefile.html#lightkurve.search.search_lightcurvefile) to look for the eclipsing binary in question. So long as there's a light curve for this contaminant — and we can reasonably expect one, for such a star — we'll be able to find it with Lightkurve's search function.\n",
    "\n",
    "(As an aside, if you run this notebook yourself, you can use the [interact_sky()](https://docs.lightkurve.org/reference/api/lightkurve.KeplerTargetPixelFile.interact_sky.html?highlight=interact_sky#lightkurve.KeplerTargetPixelFile.interact_sky) function, which returns an interactive TPF widget with targets labelled by _Gaia_ ID. This includes many targets which were not collected by _Kepler_ or _K2_.)\n",
    "\n",
    "Making a rough guess from the size of the TPF above, we set a search radius of 20 arcseconds. And we're only going to search Quarter 11, as otherwise the function will return an entry for every available quarter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 100
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 3181,
     "status": "ok",
     "timestamp": 1598468103993,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "1G-JXuYr8Ka-",
    "outputId": "68e10794-1eaf-4ca5-9243-a3d1aae61908"
   },
   "outputs": [],
   "source": [
    "lk.search_lightcurve('KIC 7024511', radius=20, quarter=11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "L0gvkvrrqt9J"
   },
   "source": [
    "Sure enough, there's a nearby target. The **distance** column in the search result tells us that this star is about 1.8 arcseconds away. Let's have a look at it:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 387
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 2558,
     "status": "ok",
     "timestamp": 1598468534973,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "aWody6wy-OOk",
    "outputId": "c4971708-33a0-4476-865f-2b83c9ada2c7"
   },
   "outputs": [],
   "source": [
    "lc_contam = lk.search_lightcurve('KIC 7024530', author='Kepler', quarter=11).download()\n",
    "lc_contam.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZEtDVc18q2FS"
   },
   "source": [
    "That's our eclipse, and it seems like there's a rotational period here too. We can also confirm this by looking in the [_Kepler_ Eclipsing Binary Catalog](http://keplerebs.villanova.edu/overview/?k=7024530), or on [Simbad](http://simbad.u-strasbg.fr/simbad/sim-basic?Ident=KIC+7024530). Simbad shows an image of the contaminant, where we can see KOI 311, fainter, just beside it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xbFgvn-G287I"
   },
   "source": [
    "## 3. Advanced: Difference Imaging"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EcdjoZpPhVDh"
   },
   "source": [
    "The focus of most _Kepler_ and _K2_ signal verification research has been in validating exoplanets. Often, this is done by measuring and tracking image centroids. A centroid is the weighted average along both axes of the image, which measures the point where most flux is concentrated. If a true transit is observed, the centroid of the image will shift in time with it ([Batalha et al., 2010](https://iopscience.iop.org/article/10.1088/2041-8205/713/2/L103/pdf)). Another technique, and the focus of this section, is difference imaging ([Bryson et al., 2013](https://iopscience.iop.org/article/10.1086/671767/pdf)). Difference imaging has been used to validate signals from exocomets ([Rappaport et al., 2018](https://arxiv.org/pdf/1708.06069.pdf)), flaring stars ([Yang et al., 2017](https://iopscience.iop.org/article/10.3847/1538-4357/aa8ea2/pdf)), and binary systems ([Colman et al., 2017](https://arxiv.org/pdf/1705.00621.pdf)).\n",
    "\n",
    "For vetting exoplanet candidates, difference imaging begins with a selection of flux measurements, both in- and out-of-phase with the transit. Using the frames from TPF images, we take an average of each collection of frames and subtract the in-phase fluxes from the out-of-phase fluxes, resulting in a _difference image_. We can then measure the centroid of the difference image, which will move outside of the optimal aperture if the star is contaminated. We can also visually identify the pixel where the variable flux is \"brightest.\" \n",
    "\n",
    "We're going to work through a basic case here, which is applicable to a broad variety of stellar signals. For any periodic and quasi-sinusoidal signal, difference imaging works by selecting flux frames at timestamps around the maxima and minima of the signal, which is equivalent to working in- and out-of-phase with a transit. The result is a pixel image where the brightest pixels are those where the periodic signal is strongest.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "cPLXuC7mH6kz"
   },
   "source": [
    "### 3.1 Determining the signal period\n",
    "\n",
    "Let's go back to KIC 2435971. We're going to use the TPF we downloaded above for difference imaging; as it's a pixel-level method, we can only use one quarter or campaign of _Kepler_ data at a time. But for difference imaging to work well, we need as much data as possible to extract the period for differencing, so we'll use the stitched light curve we prepared earlier. We've already downloaded the data above, so let's remind ourselves of what the periodogram looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 401
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1177,
     "status": "ok",
     "timestamp": 1598468941005,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "pxHdJOXh5JTt",
    "outputId": "8f9ba8bd-7cc5-4928-9983-7f74e9613c22"
   },
   "outputs": [],
   "source": [
    "pg.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nPVQSqedNO47"
   },
   "source": [
    "From this, we can extract the period of the highest amplitude signal; you will have seen the [period_at_max_power](https://docs.lightkurve.org/reference/api/lightkurve.periodogram.Periodogram.period_at_max_power.html?highlight=period_at_max_power) argument in one of the periodogram tutorials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 52
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 403,
     "status": "ok",
     "timestamp": 1598468975624,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "_jZmoWaZ5WOc",
    "outputId": "372bde2f-4eab-4e03-d42d-a814ac997898"
   },
   "outputs": [],
   "source": [
    "peak = pg.frequency_at_max_power\n",
    "period = pg.period_at_max_power\n",
    "print(f'Frequency: {peak:.3f}\\nPeriod: {period:.3f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "iPRPeUtKIkn4"
   },
   "source": [
    "### 3.2 Determining the maximum and minimum phase times\n",
    "\n",
    "To find our maxima and minima, we're going to start by phase-folding the light curve. This will help us identify the maxima and minima of the periodic signal.\n",
    "\n",
    "Note that when we phase-fold the curve below, we use the `epoch_time` argument to shift the phase curve's zero position. This helps us pick out the maxima and minima of the phase curve more clearly; however, it's a matter of trial and error to choose a suitable epoch time for any given target. In many cases, it will suffice not to set the argument at all. An alternative to experimenting with the epoch time is appending the phased light curve to itself, which guarantees at least one clear maximum and one clear minimum. This is a better method for an automated difference imaging pipeline, but is more complicated to implement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 387
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1337,
     "status": "ok",
     "timestamp": 1598469004968,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "S13VV-P49qn3",
    "outputId": "4a401f25-edd4-43de-a0dd-9848060d70f5"
   },
   "outputs": [],
   "source": [
    "folded = lc.fold(period, epoch_time=0.5)\n",
    "folded.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Sm-oj-rZODLF"
   },
   "source": [
    "This looks promising, but there's still a great deal of noise in the light curve. When we smooth the light curve by putting it into 0.001-day bins, we can see the sinusoidal trend more clearly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 387
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 3043,
     "status": "ok",
     "timestamp": 1598469053579,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "pzn2NtWZ40E1",
    "outputId": "e261a9bc-9a08-47f4-ca85-96d0e6751717"
   },
   "outputs": [],
   "source": [
    "folded2 = folded.bin(time_bin_size=0.001)\n",
    "folded2.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ESDkgeNGOe1K"
   },
   "source": [
    "We're going to identify the maxima and minima using the binned phase curve. Once we have identified two ranges in phase space, we'll collect the timestamps within those ranges from the original phase curve. Then we'll compare those timestamps to our TPF data to pick out the corresponding flux frames for each range.\n",
    "\n",
    "It can be tricky to decide how much flux to collect on either side of the maxima and minima. A good value is a tolerance of 5% of the phase curve, which means that we end up using 10% of the light curve at maximum phase and 10% of the light curve at minimum phase for the difference image. This ensures that we're using enough to get a reliable difference image, given some uncertainty in where the maxima and minima actually are, but not so much that the difference image is meaningless.\n",
    "\n",
    "This next part uses some Python tricks to help us quickly determine which timestamps to use in calculating our difference image. If you find the code below unfamiliar, it might be helpful to read up on [NumPy `where()`](https://numpy.org/doc/stable/reference/generated/numpy.where.html) and [list comprehension](https://www.pythonforbeginners.com/basics/list-comprehensions-in-python)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Be-1P4xW87j1"
   },
   "outputs": [],
   "source": [
    "full_phase_range = folded2.phase[-1].value - folded2.phase[0].value\n",
    "tolerance = 0.05 * full_phase_range\n",
    "min_phase = folded2.time[np.argmin(folded2.flux)].value \n",
    "max_phase = folded2.time[np.argmax(folded2.flux)].value\n",
    "\n",
    "min_timestamps = folded.time_original[np.where((folded.time > min_phase - tolerance) \n",
    "                                             & (folded.time < min_phase + tolerance))].value\n",
    "max_timestamps = folded.time_original[np.where((folded.time > max_phase - tolerance)\n",
    "                                             & (folded.time < max_phase + tolerance))].value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-QE794z2c2mh"
   },
   "outputs": [],
   "source": [
    "one_quarter_minima = [f for (f, t) in zip(tpf.flux.value, tpf.time.value) if t in min_timestamps]\n",
    "one_quarter_maxima = [f for (f, t) in zip(tpf.flux.value, tpf.time.value) if t in max_timestamps]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "KR7LD70HQ8cl"
   },
   "source": [
    "### 3.3 Computing the difference image\n",
    "\n",
    "Now that we have identified the maximum and minimum phase timestamps to use, we can calculate our difference image! We're also going to calculate an average of the whole quarter, and we can use Matplotlib to display them side by side for clear comparison. Note that we're also flipping the images, so that they have the same orientation as the Lightkurve TPF plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "IUKOD91Gf0M5"
   },
   "outputs": [],
   "source": [
    "avg_image = np.nanmean(tpf.flux.value, axis=0)\n",
    "diff_image = np.abs(np.nanmean(one_quarter_maxima, axis=0) - np.nanmean(one_quarter_minima, axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 356
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 871,
     "status": "ok",
     "timestamp": 1598469102188,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "YK7EQozigHRi",
    "outputId": "16bb2e37-6835-482f-b4db-790bd8c964b4"
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,2)\n",
    "ax[0].imshow(np.flipud(avg_image))\n",
    "ax[0].set_title('Quarter 9 average')\n",
    "ax[1].imshow(np.flipud(diff_image))\n",
    "ax[1].set_title('Quarter 9 difference image')\n",
    "fig.set_size_inches((15,6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "kGRooweySGkG"
   },
   "source": [
    "Let's take a moment to think about what these images tell us. In the average image, we can see our two stars, the target and the brighter star at the top of the frame. In the difference image, the pixel with the highest difference flux is one of the central pixels, indicating that the signal is coming from the target star, just as we saw when we looked at pixel-level periodograms with `plot_pixels()`.\n",
    "\n",
    "We can still see a little bit of signal in some other pixels. This is down to noise, which we saw in our first folded light curve. Even with a clear signal like the one from KIC 2435971, most difference images will show some degree of noise. But it's good to note that the difference flux in the pixel that hosts the bright nearby star is greatly diminished, as we would expect.\n",
    "\n",
    "So, since we can get this information from `plot_pixels()`, why would we want to go to the effort of difference imaging? Sometimes the signal in a periodogram may not be visually clear from using `plot_pixels()`, but there would still be enough signal for it to show up in a difference image. Difference imaging is also useful for clarifying association in a crowded field, or where bright stars are involved. If a signal is spread out across a lot of pixels, difference imaging can pinpoint where it's strongest; this is particularly evident in the case of KIC 2435971. And of course, it's always good to confirm our conclusions using a different method!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "vbge5j28DVeN"
   },
   "source": [
    "## About this Notebook\n",
    "\n",
    "**Author:** [Isabel Colman](http://orcid.org/0000-0001-8196-516X) (`isabel.colman@sydney.edu.au`)\n",
    "\n",
    "**Updated on:** 2020-09-15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ZANsIso_B_si"
   },
   "source": [
    "## Citing Lightkurve and Astropy\n",
    "\n",
    "If you use `lightkurve` or `astropy` for published research, please cite the authors. Click the buttons below to copy BibTeX entries to your clipboard."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 85
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 94159,
     "status": "ok",
     "timestamp": 1598467254937,
     "user": {
      "displayName": "Geert Barentsen",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Gj8sjdnDeqdejfe7OoouYPIclAQV0KSTpsU469Jyeo=s64",
      "userId": "05704237875861987058"
     },
     "user_tz": 420
    },
    "id": "7vUtrWVjnlY7",
    "outputId": "1bb3c117-fcce-402d-ce33-663c88e581da"
   },
   "outputs": [],
   "source": [
    "lk.show_citation_instructions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "y1HbVfHwDVeN"
   },
   "source": [
    "<img style=\"float: right;\" src=\"https://raw.githubusercontent.com/spacetelescope/notebooks/master/assets/stsci_pri_combo_mark_horizonal_white_bkgd.png\" alt=\"Space Telescope Logo\" width=\"200px\"/>"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "How to verify the location of a signal in Kepler pixel data.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}