09_Single_Star_Lightcurve_with_Butler.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ae30e255",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<br>\n",
    "<b>Single Star Lightcurve with the Butler</b> <br>\n",
    "Last verified to run on 2021-10-13 with LSST Science Pipelines release w_2021_40 <br>\n",
    "Contact authors: Melissa Graham, Jeff Carlin <br>\n",
    "Target audience: All DP0 delegates. <br>\n",
    "Container Size: medium <br>\n",
    "Questions welcome at <a href=\"https://community.lsst.org/c/support/dp0\">community.lsst.org/c/support/dp0</a>. <br>\n",
    "Find DP0 documentation and resources at <a href=\"https://dp0-1.lsst.io\">dp0-1.lsst.io</a>. <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "644abc24",
   "metadata": {},
   "source": [
    "**Credit:** This tutorial notebook was created by Jeff Carlin, Leanne Guy, and Melissa Graham, and draws on material in the fourth notebook \"04 Intro to Butler\". Please consider acknowledging Jeff, Leanne, and Melissa in any publications or software releases that make use of this notebook's contents."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9a2474f",
   "metadata": {},
   "source": [
    "### Caveats\n",
    "\n",
    "The DP0.1 data set is not well suited to large-scale time domain analyses, and one goal of this notebook is to clearly illustrate the limitations.\n",
    "\n",
    "The DP0.1 data set does not contain any data products from difference image analysis.\n",
    "It is recommended that DP0 delegates who would like to do time-domain science wait for DP0.2.\n",
    "\n",
    "A catalog of sources detected in processed visit images (PVIs) is available via the Butler (read more in the [DP0.1 DPDD](https://dp0-1.lsst.io/data-products-dp0-1/index.html#dp0-1-data-products-definition-document-dpdd)).\n",
    "However, this source catalog does not include *associations* of PVI sources by sky coordinate, nor associations between PVI sources and objects detected in the deep coadded images.\n",
    "In this notebook, the user must do the spatial association to build the lightcurve.\n",
    "\n",
    "It is recommended that DP0 delegates work through the Butler notebooks \"01 Intro to DP0 Notebooks\" and \"04 Intro to Butler\" before using this notebook.\n",
    "\n",
    "### Learning Objectives\n",
    "\n",
    "To use the Butler to retrieve all detected sources from all processed visit images (PVIs or `calexps`) at a given sky coordinate, and plot magnitude as a function of time."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35383363",
   "metadata": {},
   "source": [
    "### Set Up\n",
    "\n",
    "Import the packages we will need for this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0270764",
   "metadata": {},
   "outputs": [],
   "source": [
    "### Rubin-specific packages\n",
    "\n",
    "import lsst.daf.butler as dafButler\n",
    "import lsst.geom as geom\n",
    "import lsst.sphgeom as sphgeom\n",
    "import lsst.daf.base as dafBase\n",
    "\n",
    "### If you want to use the TAP Service in Section 2, uncomment.\n",
    "# from lsst.rsp import get_tap_service\n",
    "# service = get_tap_service()\n",
    "\n",
    "### General python / astronomy packages\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from astropy.coordinates import SkyCoord\n",
    "import astropy.units as u\n",
    "from astropy.time import Time\n",
    "from astropy.table import Table\n",
    "from astropy.io import fits\n",
    "from astropy.timeseries import LombScargle\n",
    "\n",
    "from tqdm.notebook import tqdm\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8ede0f0",
   "metadata": {},
   "source": [
    "Create an instance of the Butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf66dada",
   "metadata": {},
   "outputs": [],
   "source": [
    "repo = 's3://butler-us-central1-dp01'\n",
    "collection = \"2.2i/runs/DP0.1\"\n",
    "butler = dafButler.Butler(repo, collections=collection)\n",
    "registry = butler.registry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a4097b9",
   "metadata": {},
   "source": [
    "### 1. Use the TAP Service to choose a variable star.\n",
    "\n",
    "If you're using this notebook for the first time, it is recommended to skip this section and use the example target in Section 2.\n",
    "\n",
    "**Warning:** If you want to use the TAP Service, uncomment the two lines in the cell under Set Up and import the TAP-related packages.\n",
    "\n",
    "Retrieve the RA and Dec of elements in the `truth_match` table that are within a 1 degree radius circle near the center of the DC2 region, and have `is_variable = 1` (is variable), `truth_type = 2` (is a star), and `mag_r > 20` (r-band AB magnitude is brighter than 20).\n",
    "Find more information about the `truth_match` table in the [DP0.1 DPDD](https://dp0-1.lsst.io/data-products-dp0-1/index.html#dp0-1-data-products-definition-document-dpdd)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e892d04d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT ra, dec \"\\\n",
    "#                          \"FROM dp01_dc2_catalogs.truth_match \"\\\n",
    "#                          \"WHERE CONTAINS(POINT('ICRS', ra, dec), \"\\\n",
    "#                          \"CIRCLE('ICRS', 62.0, -37.0, 1.0)) = 1 \"\\\n",
    "#                          \"AND is_variable = 1 AND truth_type = 2 AND mag_r > 20 \",\\\n",
    "#                          maxrec=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4c6f15f",
   "metadata": {},
   "source": [
    "Alternatively, you might want to choose a Type Ia supernova, but note that Section 5 of this notebook assumes you've retrieved data for a variable star.\n",
    "\n",
    "To use the TAP service to find Type Ia supernovae, set: truth_type = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84dfbad5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# results = service.search(\"SELECT ra, dec, host_galaxy \"\\\n",
    "#                          \"FROM dp01_dc2_catalogs.truth_match \"\\\n",
    "#                          \"WHERE CONTAINS(POINT('ICRS', ra, dec), \"\\\n",
    "#                          \"CIRCLE('ICRS', 62.0, -37.0, 1.0)) = 1 \"\\\n",
    "#                          \"AND truth_type = 3 \", maxrec=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4020bc8b",
   "metadata": {},
   "source": [
    "Convert the results to a pandas data frame and show it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5575e0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# data = results.to_table().to_pandas()\n",
    "# data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "750d08db",
   "metadata": {},
   "source": [
    "Clean up, we don't need these anymore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84730721",
   "metadata": {},
   "outputs": [],
   "source": [
    "# del results,data\n",
    "# del service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95b0576a",
   "metadata": {},
   "source": [
    "### 2. Execute spatial search of the `src` catalogs with the Butler.\n",
    "\n",
    "For now, we will select a star that is known to be a simulated RR Lyrae variable. The star of interest is located at a position (RA, Dec) = (60.2837946, -35.4042439) degrees. \n",
    "\n",
    "Use [astropy.SkyCoord](https://docs.astropy.org/en/stable/api/astropy.coordinates.SkyCoord.html) to define a coordinate object for this variable star."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7da198cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "targ_coord = SkyCoord(ra=60.2837946 * u.deg, dec=-35.4042439 * u.deg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b94f0386",
   "metadata": {},
   "source": [
    "Use these coordinates to define the HTM ID spatial search region to pass to the Butler's `queryDatasets` function.\n",
    "\n",
    "Find more details about sky pixelization with HTM, and Butler spatial queries and the `dataId`, in Section 2.5 of notebook \"04 Intro to Butler\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70f53f27",
   "metadata": {},
   "outputs": [],
   "source": [
    "pixelization = sphgeom.HtmPixelization(15)\n",
    "htm_id = pixelization.index(sphgeom.UnitVector3d(sphgeom.LonLat.fromDegrees(\n",
    "                            targ_coord.ra.value, targ_coord.dec.value)))\n",
    "\n",
    "# Obtain and print the scale to provide a sense of the size of the\n",
    "#   sky pixelization being used\n",
    "scale = pixelization.triangle(htm_id).getBoundingCircle().getOpeningAngle().asDegrees() * 3600\n",
    "print(f'HTM ID={htm_id} at level={pixelization.getLevel()} is a ~{scale:0.2}\" triangle.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "502aaf06",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Now we will query the `src` datasets (i.e., measurements from each processed visit image) for this HTM pixel via the Butler. Notice that the ID for the HTM triangle of interest was extracted using the RA, Dec coordinates of the star of interest. When we extract dataset references, we can limit it to return only catalogs overlapping that HTM \"pixel.\"\n",
    "\n",
    "To keep things simple, we will just query for r-band data.\n",
    "\n",
    "Find more details about the `src` catalog in the [DP0.1 DPDD](https://dp0-1.lsst.io/data-products-dp0-1/index.html#dp0-1-data-products-definition-document-dpdd)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "628d2d85",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasetRefs = registry.queryDatasets(\"src\", htm20=htm_id,\n",
    "                                     collections=collection,\n",
    "                                     where=\"band in ('r')\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93c17716",
   "metadata": {},
   "outputs": [],
   "source": [
    "datasetRefs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd2c0468",
   "metadata": {},
   "source": [
    "Extract the datasetRefs into a python list called `refs` to speed up access to the dataIds in Section 3. This may take a little while (~30 seconds)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3a02fb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "refs = list(datasetRefs)\n",
    "totalNrefs = len(refs)\n",
    "\n",
    "print(totalNrefs, ' catalogs matching the requested position.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb94ab48",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "Let's take a look at a `dataId` from the list we have compiled. It should specify the `visit` and `detector` containing the star, and the `band` (r) that we specified in the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55013d81",
   "metadata": {},
   "outputs": [],
   "source": [
    "refs[0].dataId.full"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfe8adb8",
   "metadata": {},
   "source": [
    "### 3. Loop over search results and store `src` data in python lists.\n",
    "\n",
    "The DP0.1 work-around for building a light curve is as follows.\n",
    "\n",
    "**WARNING:** This is *not representative* of how LSST light curves will be built and accessed by users. This is a *temporary work-around* for DP0.1.\n",
    "\n",
    " 1. For each query result (i.e., each element of `refs`), use the `dataId` to retrieve all sources for that query result.\n",
    " 2. Calculate the separation of each source from the target.\n",
    " 3. If the nearest source is within 1.5\" of the target, assume it _is_ the target.\n",
    " 4. Append the following quantities for the nearest source into a python list:\n",
    "   * -- RA and Dec\n",
    "   * -- separation (between the supplied position and the measured position in each visit image)\n",
    "   * -- r-band magnitude\n",
    "   * -- visit and detector IDs\n",
    "   * -- date of observation (as a \"Modified Julian Date,\" or MJD)\n",
    "\n",
    "**The above process can take ~30 minutes: choose option A to do it anyway, and option B to skip to plotting a pre-made lightcurve.**\n",
    "\n",
    "**Option A: retrieve source photometry via the Butler:** <br>\n",
    "Uncomment and execute the following two cells, which perform the steps listed above.  <br>\n",
    "**WARNING**: this takes a _long_ time, sometimes up to 30 minutes, to run. <br>\n",
    "**WARNING:** If you choose option A, scroll down and comment-out the cell for option B.\n",
    "\n",
    "**Option B (RECOMMENDED): read the previously-generated data file (this is the default option for this notebook):** <br>\n",
    "Scroll down and execute the cell under OPTION B. \n",
    "It will read in a FITS file containing the results previously created by OPTION A."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3b11fc0",
   "metadata": {},
   "source": [
    "#### 3.1 OPTION A: Retrive source photometry via the Butler (can take ~30 minutes)\n",
    "\n",
    "Uncomment the following two cells and execute them in order to build the light curve yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f524e85c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# t0 = time.time()\n",
    "# print(t0)\n",
    "\n",
    "# # To retrieve data for all `refs`, set N_refs to be equal to totalNrefs (or to a large value like 1000).\n",
    "# # Alternatively to test the retrieval of a few `refs`, set smaller (e.g., 10).\n",
    "# Nrefs = totalNrefs\n",
    "# Nrefs = 15\n",
    "\n",
    "# # Instantiate empty lists\n",
    "# ra_arr = []\n",
    "# dec_arr = []\n",
    "# sep_arr = []\n",
    "# mag_arr = []\n",
    "# visit_arr = []\n",
    "# detector_arr = []\n",
    "# mjd_arr = []\n",
    "\n",
    "# # Loop over all refs\n",
    "# for i, d in enumerate(tqdm(refs)):\n",
    "#     t1 = time.time()\n",
    "\n",
    "#     if i <= Nrefs:\n",
    "\n",
    "#         # Use the butler to get all sources for this datasetRef's dataId\n",
    "#         did = d.dataId\n",
    "#         src = butler.get('src', dataId=did)\n",
    "\n",
    "#         # Get the separation of all sources from the target\n",
    "#         src_coords = SkyCoord(ra=src['coord_ra'] * u.rad,\n",
    "#                               dec=src['coord_dec'] * u.rad)\n",
    "#         sep = src_coords.separation(targ_coord)\n",
    "\n",
    "#         # If the nearest source is within 1.5\", append source quantities to python lists\n",
    "#         if np.min(sep.arcsecond) < 1.5:\n",
    "#             sx = np.argmin(sep.arcsecond)\n",
    "\n",
    "#             # Append RA, Dec, and separation\n",
    "#             ra_arr.append(src['coord_ra'][sx])\n",
    "#             dec_arr.append(src['coord_dec'][sx])\n",
    "#             sep_arr.append(sep[sx].arcsecond)\n",
    "\n",
    "#             # Append r-band magnitude (AB mag from the calibrated flux in nJy).\n",
    "#             # The calibrated flux in nJy is base_PsfFlux_instFlux times base_localPhotoCalib\n",
    "#             mag_arr.append(-2.5 * np.log10(src['base_PsfFlux_instFlux'][sx] * src['base_localPhotoCalib'][sx]) + 31.4)\n",
    "\n",
    "#             # Append visit and detector information\n",
    "#             visit_arr.append(did['visit'])\n",
    "#             detector_arr.append(did['detector'])\n",
    "\n",
    "#             # Append MJD (from this src's associated image's header)\n",
    "#             visit_info = butler.get('calexp.visitInfo', dataId=did)\n",
    "#             mjd_arr.append(visit_info.getDate().get(dafBase.DateTime.MJD))\n",
    "\n",
    "#             del sx\n",
    "\n",
    "#         t2 = time.time()\n",
    "#         if i == 0:\n",
    "#             print('Each data point takes: ', t2 - t1, ' seconds to retrieve.')\n",
    "#             print('Retrieving ', Nrefs, ' data points will take ',\n",
    "#                   (t2 - t1) * Nrefs / 60.0, ' minutes.')\n",
    "\n",
    "#         del did, src, src_coords, sep\n",
    "#         del t1, t2\n",
    "\n",
    "# print('Total run-time was: ', (time.time() - t0) / 60.0, ' minutes.')\n",
    "# del t0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a67cdd2d",
   "metadata": {},
   "source": [
    "<blockquote>\n",
    "As a side-note, there are options for how the Butler can be used to get the photometric calibration and visit info for a source.\n",
    "One alternative is to retrieve the full `calexp` to get the `photCalib` and `visitInfo`.\n",
    "\n",
    "Retrieve the `calexp`:<br>\n",
    "``calexp = butler.get('calexp', dataId=did, collections=collection)``\n",
    "\n",
    "Apply the `calexp`'s photometric calibration to the flux:<br>\n",
    "``calib = calexp.getPhotoCalib()``<br>\n",
    "``mag_arr.append(calib.instFluxToMagnitude(src['base_PsfFlux_instFlux'][sx]))``\n",
    "\n",
    "Convert the `calexp`'s header date to MJD:<br>\n",
    "``date_obs = calexp.getInfo().getVisitInfo().getDate()``<br>\n",
    "``mjd_arr.append(date_obs.get(dafBase.DateTime.MJD))``\n",
    "\n",
    "In general, the Butler is very flexible software and there are often multiple ways to achieve the same goal.\n",
    "</blockquote>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec877cd9",
   "metadata": {},
   "source": [
    "Make an AstroPy Table of the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83226e3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# tab_timeseries = Table([mjd_arr, ra_arr, dec_arr, sep_arr, mag_arr,\n",
    "#                         visit_arr, detector_arr],\n",
    "#                        names=['mjd', 'ra', 'dec', 'separation', 'mag',\n",
    "#                               'visit', 'detector'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9d5f07b",
   "metadata": {},
   "source": [
    "Save the timeseries data to a FITS file for later use.\n",
    "\n",
    "**WARNING:** This will overwrite the provided file \"data/timeseries_rband.fits\" unless you rename it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "516d24f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# tab_timeseries.write('data/timeseries_rband.fits', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e270d5cf",
   "metadata": {},
   "source": [
    "#### 3.2 OPTION B: Read the previously-generated data file. \n",
    "\n",
    "**WARNING:** Comment this out if you chose to do option A."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca7749b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "tab_timeseries = fits.getdata('data/timeseries_rband.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1bbf927",
   "metadata": {},
   "source": [
    "#### 3.3 After option A OR option B:\n",
    "Confirm that the expected columns are in the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49c1b3a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "tab_timeseries.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7848a3c8",
   "metadata": {},
   "source": [
    "### 4. Plot a time series (magnitude vs. time) of the star's measurements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4c04b2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(tab_timeseries['mjd'], tab_timeseries['mag'], 'k.')\n",
    "plt.xlabel('MJD')\n",
    "plt.ylabel('r magnitude')\n",
    "plt.gca().invert_yaxis()\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1d0d43a",
   "metadata": {},
   "source": [
    "Hooray! The star is clearly a variable, with a peak-to-peak variation of almost a full magnitude. We have successfully extracted a time series of a variable star!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10b7ce59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# What timespan do these measurements cover?\n",
    "time_days = np.max(tab_timeseries['mjd']) - np.min(tab_timeseries['mjd'])\n",
    "time_years = time_days / 365.25\n",
    "print('Time span in days: ', time_days, '; Time in years: ', time_years)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33dedff7",
   "metadata": {},
   "source": [
    "### 5. Bonus: find the star's period and make a phase-folded light curve\n",
    "\n",
    "In the above plot, it's clear that the star is variable, but with ~100 measurements spread over ~1750 days (almost 5 years), its period is unclear. To extract a period from this time series, we will need to construct a _periodogram_, identify the most likely period, then \"fold\" the lightcurve by placing each measurement at the appropriate phase (i.e., between the start of the variability phase at 0.0 and the end of the cycle at 1.0).\n",
    "\n",
    "_Major caveat: the CET members who created this notebook are not variable star experts, so the method shown here is definitely not the \"best\", state-of-the-art way of doing this. Contributions that improve on this are welcome from delegates!_\n",
    "\n",
    "We will use a \"Lomb-Scargle Periodogram,\" which is a statistical method for extracting periodicity from time series data. In particular, we'll use the [Astropy Lomb-Scargle](https://docs.astropy.org/en/stable/timeseries/lombscargle.html) implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60d08aeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make arrays with units attached:\n",
    "\n",
    "mjd_days = np.array(tab_timeseries['mjd']) * u.day\n",
    "mags = np.array(tab_timeseries['mag']) * u.mag\n",
    "\n",
    "# The Lomb-Scargle periodogram returns the power at different frequencies.\n",
    "\n",
    "# Because we know this is an RR Lyrae star, its period must be between\n",
    "# ~0.2-0.9 days. We will use the min/max frequency settings to limit our\n",
    "# period search to 0.05-1.25 days.\n",
    "\n",
    "min_freq_search = 1.0 / (1.25 * u.day)\n",
    "max_freq_search = 1.0 / (0.05 * u.day)\n",
    "\n",
    "# frequency, power = LombScargle(mjd_days, mags).autopower()\n",
    "frequency, power = LombScargle(mjd_days, mags).autopower(minimum_frequency=min_freq_search,\n",
    "                                                         maximum_frequency=max_freq_search)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "598b86c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the frequency spectrum:\n",
    "fig, ax = plt.subplots(1, 2, figsize=(12, 7))\n",
    "\n",
    "plt.sca(ax[0])  # set the first axis as current\n",
    "\n",
    "plt.plot(frequency, power)\n",
    "plt.minorticks_on()\n",
    "plt.xlabel('frequency (1/d)')\n",
    "plt.ylabel('power')\n",
    "\n",
    "plt.sca(ax[1])  # set the second axis as current\n",
    "plt.plot(1 / frequency, power)\n",
    "plt.minorticks_on()\n",
    "plt.xlim(0, 1.25)\n",
    "plt.xlabel('period (d)')\n",
    "plt.ylabel('power')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fb90bbc",
   "metadata": {},
   "source": [
    "In an ideal case, the peak showing the most power in the periodogram represents the most likely period. (This is definitely not always the case in the presence of poor time sampling, large photometric errors, etc.)\n",
    "\n",
    "Let's extract the peak at a period of ~0.5 days, which has by far the largest power (>0.7)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "821676fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select power>0.4 and period (or 1/frequency) within 0.1 days of 0.5.\n",
    "select_peak = (power > 0.4) & (np.abs(0.5 - (1 / frequency.value)) < 0.1)\n",
    "\n",
    "print(1 / frequency[select_peak])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4be087a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Since there are multiple points satisfying the condition, we will take the\n",
    "#   mean of their periods as the \"best\" period.\n",
    "\n",
    "best_period = np.mean(1 / frequency[select_peak])\n",
    "print(best_period)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1df20208",
   "metadata": {},
   "source": [
    "<br>\n",
    "Now we have a period. To place each measurement in phase using this period, we will set the first measurement to be time = 0. Then we calculate how many periods have elapsed since that time for each measurement. The \"remainder\" left over after removing full periods is the \"phase\" of each observation.\n",
    "\n",
    "For example, say the phase is 0.8 days. We have three measurements at [0, 1.0, 2.0] days. This corresponds to [0/0.8, 1.0/0.8, 2.0/0.8] = [0, 1.25, 2.5] periods. Removing the fully completed periods from the last two measurements places them at phases of 0.25 and 0.5.\n",
    "\n",
    "Let's apply this to our data and see what happens:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f934bb9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of elapsed periods since the first measurement:\n",
    "mjd_norm = (mjd_days - np.min(mjd_days)) / best_period\n",
    "\n",
    "# Calculate phase by using the modulus function:\n",
    "phase = np.mod(mjd_norm, 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3b5059c",
   "metadata": {},
   "source": [
    "<br>\n",
    "Plot a phased light curve -- drumroll please..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "639c9d12",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(phase, mags, 'k.')\n",
    "plt.gca().invert_yaxis()\n",
    "plt.xlabel('phase')\n",
    "plt.ylabel('r magnitude')\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd40acf5",
   "metadata": {},
   "source": [
    "Awesome - it looks like a typical RR Lyrae light curve (specifically, an _ab_-type RR Lyrae; see [this Wikipedia page](https://en.wikipedia.org/wiki/RR_Lyrae_variable) for more about these variable stars).\n",
    "\n",
    "Note that by selecting the first measurement as our \"zero-point\" for the phase, we have chosen an arbitrary starting point. Ideally one might choose something more standard, like placing the peak of the light curve at phase = 0.\n",
    "\n",
    "\n",
    "#### Summary \n",
    "You have now seen how to extract and work with time series data for a variable star. As you've seen, this is rather clunky with the DP0.1 dataset and tools. Extracting individual measurements from visits will be much easier in planned data products (e.g., in DP0.2), where the forced photometry on each visit will be tabulated in its own type of table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07cbca6e-e0d3-41a5-8dd1-543c2650d2b5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}