DP02_07b_Variable_Star_Lightcurves.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "081f3946-bafe-4c43-a766-c08673a2773b",
   "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>Variable Stars in DP0.2</b> <br>\n",
    "Contact author(s): Jeff Carlin and Ryan Lau <br>\n",
    "Last verified to run: 2024-05-01 <br>\n",
    "LSST Science Pipelines version: Weekly 2024_16 <br>\n",
    "Container Size: medium <br>\n",
    "Targeted learning level: intermediate <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ba3304f-ee93-4469-978b-a98ff0f2bbb4",
   "metadata": {},
   "source": [
    "**Description:** Use catalog data to identify variable stars and plot their lightcurves."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58e1fddb-395d-4aad-971f-f120dab415ca",
   "metadata": {},
   "source": [
    "**Skills:** Use various TAP tables, including joining multiple tables. Extract time-series photometry. Measure periods and plot phased lightcurves."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7980fe9c-9d9a-429d-ab4b-4812c6fef3fe",
   "metadata": {},
   "source": [
    "**LSST Data Products:** TAP tables dp02_dc2_catalogs.Object, ForcedSource, CcdVisit, DiaObject, DiaSource, ForcedSourceOnDiaObject<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e317548d-7f06-4b65-928a-bce2bd5a95e6",
   "metadata": {},
   "source": [
    "**Packages:** numpy, matplotlib, astropy.units, astropy.coordinates, astropy.io.fits, astropy.timeseries.LombScargle, lsst.rsp.get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd9d68a6-93b9-42dc-8291-88564cd53563",
   "metadata": {},
   "source": [
    "**Credit:** Originally developed by Jeff Carlin and the Community Science Team for Rubin Data Preview 0 with improvements contributed by Bob Abel and updates contributed by Ryan Lau.<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f2912bf-a172-4f6e-b8ce-38e03a253d68",
   "metadata": {},
   "source": [
    "**Get Support:**\n",
    "Find DP0-related documentation and resources at <a href=\"https://dp0.lsst.io\">dp0.lsst.io</a>. Questions are welcome as new topics in the <a href=\"https://community.lsst.org/c/support/dp0\">Support - Data Preview 0 Category</a> of the Rubin Community Forum. Rubin staff will respond to all questions posted there."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9413598-6c09-45b5-baf9-3b0e9a38da44",
   "metadata": {},
   "source": [
    "## 1. Introduction\n",
    "\n",
    "This notebook demonstrates use of the Table Access Protocol (TAP) service to query the time-domain data products at the location of a known RR Lyrae variable star.\n",
    "\n",
    "In Section 2 a known RR Lyrae star is identified in the `Object` table, and then photometry is retrieved from the `ForcedSource` table.\n",
    "Recall that the `Object` table contains astrometric and photometric measurements for objects detected in coadded images, and that\n",
    "the `ForcedSource` table contains forced photometry on the individual processed visit images (PVIs; direct images)\n",
    "at the locations of all detected `Objects`. \n",
    "A Lomb-Scargle periodogram is used to identify the period and phase-fold the time-series photometry, and then\n",
    "a phase-folded lightcurve is displayed for the single RR Lyrae.\n",
    "\n",
    "In Section 3 the data products of difference image analysis (DIA), in particular the lightcurve summary statistic parameters\n",
    "available in the `DiaObject` table, are used to identify a sample of likely RR Lyrae from the DP0.2 data set. \n",
    "\n",
    "A cautionary tale on the difference between using time-series photometry for `DiaObjects` in the `DiaSource` or `ForcedSourceOnDiaObject` tables is demonstrated in Section 3.4.\n",
    "Both contain forced photometry in the PVI (direct image), but the `DiaSource` table only contains photometry for visits in which the\n",
    "`DiaObject` was *detected in the corresponding difference image*,\n",
    "whereas the `ForcedSourceOnDiaObject` table contains forced photometry for *all* PVIs and difference images.\n",
    "\n",
    "> **Warning:** When variable stars are observed to have nearly the same brightness as they do in the template image,\n",
    "they will not be detected in the difference image,\n",
    "and photometry for that visit's PVI will not be included in the `DiaSource` table.\n",
    "\n",
    "This notebook demonstrates that, for variable stars, it is better to use either the `ForcedSource` or `ForcedSourceOnDiaObject` tables."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "694404da-f100-474f-826a-02b31706d878",
   "metadata": {},
   "source": [
    "### 1.1 Package Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "937c4676-d153-49d0-965d-9e4db0d2c019",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.units as u\n",
    "from astropy.table import unique\n",
    "from astropy.timeseries import LombScargle\n",
    "from lsst.rsp import get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b16a970-4179-4ce7-9e99-a59e807401d2",
   "metadata": {},
   "source": [
    "### 1.2 Define Functions and Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0571eaa5-bcce-4ba2-8dac-ee9e25810f69",
   "metadata": {},
   "source": [
    "Set up some plotting defaults so plots will look nice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12d3c8bb-0365-4d28-823a-ab161433aa51",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.style.use('tableau-colorblind10')\n",
    "params = {'axes.labelsize': 24,\n",
    "          'font.size': 20,\n",
    "          'legend.fontsize': 14,\n",
    "          'xtick.major.width': 3,\n",
    "          'xtick.minor.width': 2,\n",
    "          'xtick.major.size': 12,\n",
    "          'xtick.minor.size': 6,\n",
    "          'xtick.direction': 'in',\n",
    "          'xtick.top': True,\n",
    "          'lines.linewidth': 3,\n",
    "          'axes.linewidth': 3,\n",
    "          'axes.labelweight': 3,\n",
    "          'axes.titleweight': 3,\n",
    "          'ytick.major.width': 3,\n",
    "          'ytick.minor.width': 2,\n",
    "          'ytick.major.size': 12,\n",
    "          'ytick.minor.size': 6,\n",
    "          'ytick.direction': 'in',\n",
    "          'ytick.right': True,\n",
    "          'figure.figsize': [10, 8],\n",
    "          'figure.facecolor': 'White'}\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e94db4b3-8feb-45c3-8f99-3140f85f7d4d",
   "metadata": {},
   "source": [
    "Set up colors and plot symbols corresponding to the _ugrizy_ bands. These colors are the same as those used for _ugrizy_ bands in Dark Energy Survey (DES) publications, and are defined in <a href=\"https://github.com/DarkEnergySurvey/descolors\">this github repository</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8538043b-bdcd-452c-9210-7d168dfecdf7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_filter_labels = ['u', 'g', 'r', 'i', 'z', 'y']\n",
    "plot_filter_colors = {'u': '#56b4e9', 'g': '#008060', 'r': '#ff4000',\n",
    "                      'i': '#850000', 'z': '#6600cc', 'y': '#000000'}\n",
    "plot_filter_symbols = {'u': 'o', 'g': '^', 'r': 'v', 'i': 's', 'z': '*', 'y': 'p'}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dbb728a-3d0b-4930-826c-2c65b1b0780d",
   "metadata": {},
   "source": [
    "Start the TAP service, which we will use for all data retrieval in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8818f0f-45b9-4066-8086-0714f054787e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service(\"tap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "056c4f66-ff8d-4a96-ac54-fc787560f0a2",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Create a Lightcurve for a Known RR Lyrae Variable\n",
    "\n",
    "To start out, we will use a star that is known (from the DC2 truth tables) to be a simulated RR Lyrae star in DC2. This will help us get our feet wet at extracting time-series data and manipulating them.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a80ab89-161f-42ac-bd6d-e14cbdc67e0c",
   "metadata": {},
   "source": [
    "### 2.1 Extract Forced Photometry Measurements for a Pulsating Variable\n",
    "\n",
    "The known RR Lyrae star we will use is at position (RA, Dec) = (62.1479031, -35.799138) degrees. We will initially look for this star in the `Object` table to identify its `objectId`. Then we will use this `objectId` to identify measurements of the star in the `ForcedSource` table.\n",
    "\n",
    "#### 2.1.1 Define star coordinates\n",
    "\n",
    "Initialize two variables with the known position of the star of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3b55053-eb68-49b1-a7d1-2048b8f6a80b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ra_known_rrl = 62.1479031\n",
    "dec_known_rrl = -35.799138"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b39ed8a4-cb0a-4f3f-8bfc-a12717375458",
   "metadata": {},
   "source": [
    "#### 2.1.2 Find the Star in the Object Table\n",
    "\n",
    "Query the `Object` table for this star. We use a spatial query at the RA, Dec position of interest, and specify a very small area to search (0.001 degrees), so that we will retrieve only the star of interest.\n",
    "\n",
    "Notice that in addition to selecting on position, we specify \"detect_isPrimary = 1\" as well. The <a href=\"https://pipelines.lsst.io/getting-started/multiband-analysis.html#filtering-for-unique-deblended-sources-with-the-detect-isprimary-flag\">detect_isPrimary flag</a> is a composite flag that, if TRUE (or set to 1), ensures a unique, deblended source.\n",
    "\n",
    "Notice also the use of the `scisql_nanojanskyToAbMag` function to convert fluxes to magnitudes. (This could also be done after extracting the fluxes if desired, but this step streamlines things.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d63566a-00b3-495d-8448-b138be849adf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(\"SELECT TOP 100 \"\n",
    "                         \"coord_ra, coord_dec, objectId, \"\n",
    "                         \"g_psfFlux, r_psfFlux, i_psfFlux, detect_isPrimary, \"\n",
    "                         \"scisql_nanojanskyToAbMag(g_psfFlux) as gPSFMag, \"\n",
    "                         \"scisql_nanojanskyToAbMag(r_psfFlux) as rPSFMag, \"\n",
    "                         \"scisql_nanojanskyToAbMag(i_psfFlux) as iPSFMag, \"\n",
    "                         \"g_extendedness, r_extendedness, i_extendedness \"\n",
    "                         \"FROM dp02_dc2_catalogs.Object \"\n",
    "                         \"WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), \"\n",
    "                         \"CIRCLE('ICRS',\" + str(ra_known_rrl) + \", \"\n",
    "                         + str(dec_known_rrl) + \", 0.001)) = 1 \"\n",
    "                         \"AND detect_isPrimary = 1\")\n",
    "\n",
    "objs = results.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd8fe632-59cf-4698-877f-683222ac86f4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "objs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29629c26-2a47-4333-a5ea-fe7a01aacee9",
   "metadata": {},
   "source": [
    "Looks good - our limited search retrieved only a single object. Notice that the \"extendedness\" flag is zero for all three bands, meaning this object is likely a point source (i.e., a star). It is around 18th magnitude in all three bands, and has fairly blue colors as expected for an RR Lyrae star.\n",
    "\n",
    "#### 2.1.3 Extract the Measurements from the Forced Source Table\n",
    "\n",
    "Now that we have identified the candidate in the Object table, we can use its `objectId` to identify it in the `ForcedSource` table. First extract the `objectId`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8256ab8-9ada-46d8-9356-61933af7a456",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sel_objid = objs[0]['objectId']\n",
    "print(sel_objid)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8295ce1-e291-4fa4-9313-5d49dab1060a",
   "metadata": {},
   "source": [
    "Now we will use this ID to pick out all measurements of this object from the `ForcedSource` table (by use of a WHERE statement in the query).\n",
    "\n",
    "Note that we are also using a JOIN here to extract the visit info for each entry in `ForcedSource`. This is necessary so we can get the time each exposure was observed (\"expMidptMJD\"), to be used in constructing lightcurves.\n",
    "\n",
    "Define the query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc4fba8f-d9e0-4315-a2c4-3dd9e367e5aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT src.band, src.ccdVisitId, src.coord_ra, src.coord_dec, \"\\\n",
    "        \"src.objectId, src.psfFlux, src.psfFluxErr, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(psfFlux) as psfMag, \"\\\n",
    "        \"visinfo.ccdVisitId, visinfo.band, \"\\\n",
    "        \"visinfo.expMidptMJD, visinfo.zeroPoint \"\\\n",
    "        \"FROM dp02_dc2_catalogs.ForcedSource as src \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.CcdVisit as visinfo \"\\\n",
    "        \"ON visinfo.ccdVisitId = src.ccdVisitId \"\\\n",
    "        \"WHERE src.objectId = \"+str(sel_objid)+\" \"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7017b75d-c8d0-45bb-bf25-6ca5b4f75c17",
   "metadata": {},
   "source": [
    "Execute the query, and extract the results to a table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f9c49a0-6b40-4ca4-b546-a2956bb6e284",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "results = service.search(query)\n",
    "srcs = results.to_table()\n",
    "print(len(srcs))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55fe56c5-f4f4-43c9-8f13-8afd26ff6963",
   "metadata": {},
   "source": [
    "Uncomment and execute the following cell to see what this table looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e804145-7308-4e01-8ef6-925ef22a01ca",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# srcs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e3b0c5e-b0bb-4fcd-ada3-719fc40d36e9",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### 2.1.4 Select Forced Photometry Measurements from Each Band\n",
    "\n",
    "The `ForcedSource` table contains measurements from all bands (i.e., _ugrizy_ filters). It will be useful to extract measurements for each band by selecting on the \"band\" column.\n",
    "\n",
    "In the following cell, and throughout this notebook, we store arrays in a python dict (in this case, one called \"pick\") that is indexed on the \"band\" or \"filter\" name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec6a7816-d415-43c0-906b-7859925f10ce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "pick = {}\n",
    "for filter in plot_filter_labels:\n",
    "    pick[filter] = (srcs['band'] == filter)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83239aaa-574c-40da-81c0-7c0471c2b08d",
   "metadata": {},
   "source": [
    "Plot the observed \"lightcurve\" for the _r_ band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0db9bbb-09cb-4c3f-888b-cc5093e793c5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "plt.plot(srcs[pick['r']]['expMidptMJD'], srcs[pick['r']]['psfMag'],\n",
    "         'k.', ms=10)\n",
    "plt.minorticks_on()\n",
    "plt.xlabel('MJD (days)')\n",
    "plt.ylabel('r')\n",
    "plt.gca().invert_yaxis()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b09713c4-093b-4013-8c96-5019f93c9422",
   "metadata": {
    "tags": []
   },
   "source": [
    "There is a spread of almost 1 magnitude, so clearly we have identified a variable star. Hooray!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cd9de32-b11b-4de3-b0f9-ea566ffcfacb",
   "metadata": {},
   "source": [
    "### 2.2 Create a Phased Lightcurve\n",
    "\n",
    "In the plot above, it is clear that this is a variable star, but we want to know more about it. Assuming it is a periodic variable, let's try to estimate its period so we can create a phased lightcurve.\n",
    "\n",
    "> **Caveat:** The author of this notebook is not an expert in periodic variables or analysis of time-series data. The following is almost certainly not the **best** way of doing these analyses, but is simply one way to do things. It is shown here as a demonstration only.\n",
    "\n",
    "A common way of searching for periodicity in unevenly-sampled time-series data is the Lomb-Scargle Periodogram. In this section we will create such a periodogram using the <a href=\"https://docs.astropy.org/en/stable/timeseries/lombscargle.html\">Lomb-Scargle package from Astropy</a>. We will then estimate the best period, and reference all observations to that period to create a phased lightcurve.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82c43b66-d544-440c-979b-df1faacea52e",
   "metadata": {},
   "source": [
    "#### 2.2.1 Extract and Plot a Periodogram\n",
    "\n",
    "For convenience, we first extract all of the times of observation and measured magnitudes to separate arrays for each bandpass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58cba6ea-6d28-491b-9fe4-349f909866a8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mjd_days = {}\n",
    "mags = {}\n",
    "for filter in plot_filter_labels:\n",
    "    mjd_days[filter] = np.array(srcs[pick[filter]]['expMidptMJD']) * u.day\n",
    "    mags[filter] = np.array(srcs[pick[filter]]['psfMag'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afd1f449-aa9b-4d7c-8b45-a24960220369",
   "metadata": {},
   "source": [
    "The Lomb-Scargle periodogram returns the power at different frequencies. Because we know this is an RR Lyrae star, its period must be between ~0.2-0.9 days. We will use the min/max frequency settings to limit our period search to 0.05-1.05 days. (Recall that frequency is 1/period.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b32afc3b-4c21-4fd1-9cee-8a3e791ee95f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "min_period = 0.05 * u.day\n",
    "max_period = 1.05 * u.day\n",
    "\n",
    "min_freq_search = 1.0 / max_period\n",
    "max_freq_search = 1.0 / min_period"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e4e1999-10ca-44ba-85fe-2965bed31e2a",
   "metadata": {},
   "source": [
    "Now we run the LombScargle algorithm on each of these datasets. We will set only the min/max frequencies, and otherwise use the \"autopower\" method from LombScargle. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36330abc-f5ad-4851-9ab2-e3ab8fa29bcf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "frequency = {}\n",
    "power = {}\n",
    "for filter in plot_filter_labels:\n",
    "    frequency[filter], power[filter] =\\\n",
    "        LombScargle(mjd_days[filter], mags[filter]).autopower(minimum_frequency=min_freq_search,\n",
    "                                                              maximum_frequency=max_freq_search)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1dd6f8b-6036-4758-9000-ed466723b6d5",
   "metadata": {},
   "source": [
    "In a well-behaved and well-sampled scenario, the Lomb-Scargle power will peak at a specific frequency that corresponds to the period of the star's variability.\n",
    "\n",
    "We will assume that the frequency returning the highest power corresponds to the real period of the star (which is not necessarily true in the case of real data with uneven sampling, period aliasing, and other effects). Let's find this frequency for each band.\n",
    "\n",
    "In the following cell, we find the index with the maximum power and assign it as the \"peakbin\" and then store the frequency according to the peak power. Then we calculate the mean frequency and period from the peak frequencies of the results from all the filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b546ca6-3b29-4a36-8ac6-c3bdbbd105a0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "all_peak_freqs = []\n",
    "\n",
    "for filter in plot_filter_labels:\n",
    "    peakbin = np.argmax(power[filter])\n",
    "\n",
    "    all_peak_freqs.append(frequency[filter][peakbin].value)\n",
    "    \n",
    "all_peak_freqs = np.array(all_peak_freqs)\n",
    "\n",
    "mean_peak_freq = np.mean(all_peak_freqs)\n",
    "\n",
    "print('Mean frequency:', mean_peak_freq)\n",
    "print('Mean period:', 1.0/mean_peak_freq, ' days')\n",
    "print('\\nugrizy frequency results:\\n', all_peak_freqs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2b62127-d994-43d0-9a06-b35eb62abd16",
   "metadata": {
    "tags": []
   },
   "source": [
    "Plot a \"periodogram\" -- a figure showing the power vs. frequency from the Lomb-Scargle analysis. We'll also add a panel for the period. (This figure shows only the r-band results. You could try it with others!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc11d6dc-c44a-4d85-8b07-680a7d9410ce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(8, 4))\n",
    "\n",
    "plt.sca(ax[0])\n",
    "plt.plot(frequency['r'], power['r'])\n",
    "plt.vlines(mean_peak_freq, 0, 1, linestyle='--', color='red')\n",
    "plt.minorticks_on()\n",
    "plt.xlabel('frequency (1/d)')\n",
    "plt.ylabel('power')\n",
    "\n",
    "plt.sca(ax[1])\n",
    "plt.plot(1 / frequency['r'], power['r'])\n",
    "plt.vlines(1/mean_peak_freq, 0, 1, linestyle='--', color='red')\n",
    "plt.minorticks_on()\n",
    "plt.xlabel('period (d)')\n",
    "plt.ylabel('power')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d752a584-bcba-4906-8615-d33b17832b76",
   "metadata": {},
   "source": [
    "#### 2.2.2 Create a Phased Lightcurve\n",
    "\n",
    "The red line we overlaid on the previous plot corresponds to a clear peak in the periodogram, so it looks like we have found a good estimate of the period. Now we will use that period to extract and plot phased light-curves of the variable star.\n",
    "\n",
    "Recall that period = 1/frequency:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5cef36d-06c5-4800-a9ed-b9eee942af03",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "best_period = 1/mean_peak_freq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0236bf09-d4eb-4b40-b41e-4968356da2c1",
   "metadata": {},
   "source": [
    "To \"phase\" the lightcurves, we calculate how many periods have passed since some fiducial time \"t0\". Here, we select t0 to be the time of the first g-band observation, but that is completely arbitrary.\n",
    "\n",
    "The second part takes only the non-integer part of the number of elapsed periods. For example, a point that happens 2.75 periods after t0 will have phase 0.75."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e57f91e1-f1a5-4ded-84c4-3e1146dc9ddf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mjd_norm = {}\n",
    "phase = {}\n",
    "\n",
    "t0 = np.min(mjd_days['g'].value)\n",
    "\n",
    "for filter in plot_filter_labels:\n",
    "    mjd_norm[filter] = (mjd_days[filter].value - t0) / best_period\n",
    "    phase[filter] = np.mod(mjd_norm[filter], 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f0cb659-6ae6-4520-875c-aba956005bd4",
   "metadata": {},
   "source": [
    "Now plot them separately for each band (using an integer offset to separate them):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce708016-79f3-4371-acec-5032fc655df7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(4, 6))\n",
    "\n",
    "i = 0\n",
    "\n",
    "for filter in plot_filter_labels:\n",
    "    plt.plot(phase[filter], mags[filter]-np.mean(mags[filter]) + i,\n",
    "             plot_filter_symbols[filter],\n",
    "             color=plot_filter_colors[filter], label=filter)\n",
    "    plt.hlines(i, 0, 1, linestyle=':', color='Gray')\n",
    "    i += 1\n",
    "\n",
    "plt.gca().invert_yaxis()\n",
    "plt.legend()\n",
    "plt.xlabel('phase')\n",
    "plt.ylabel('mag - mean (+offset)')\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32b3a216-8f69-440e-9f98-d1b833b287f3",
   "metadata": {
    "tags": []
   },
   "source": [
    "...and all together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1421c45-60f7-4a35-a87f-2bd23f7bc851",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(6, 4))\n",
    "\n",
    "for filter in plot_filter_labels:\n",
    "    plt.plot(phase[filter], mags[filter]-np.mean(mags[filter]),\n",
    "             plot_filter_symbols[filter],\n",
    "             color=plot_filter_colors[filter], label=filter)\n",
    "plt.hlines(0, 0, 1, linestyle=':', color='Gray')\n",
    "plt.gca().invert_yaxis()\n",
    "plt.legend()\n",
    "plt.xlabel('phase')\n",
    "plt.ylabel('mag - mean')\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1992314a-42fc-4cfc-a1f6-bd4c224e969f",
   "metadata": {
    "tags": []
   },
   "source": [
    "Success! That looks like an RR Lyrae star. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7061ac0e-996e-48e7-8d6b-3a302c2e1a42",
   "metadata": {},
   "source": [
    "## 3. Identify Candidate RR Lyrae From Catalogs\n",
    "\n",
    "Now that we know we can extract a lightcurve, find the star's period, and plot phased lightcurves for a _known_ RR Lyrae star, we move on to exploring how to _find_ candidate variables based on data in the catalogs.\n",
    "\n",
    "Data products from difference image analysis (DIA), where sources are detected and measured on the _difference_ image resulting from subtracting an approriately scaled deep coadd template of the same sky area from a processed visit image (PVI), are a good place to check for transient and variable objects.\n",
    "\n",
    "The three DIA products we will use are:\n",
    "\n",
    "- `DiaObject`: table of all spatially unique objects detected in all difference images (all `DiaSources` associated by coordinate)\n",
    "- `DiaSource`: table of sources detected in the difference images\n",
    "- `ForcedSourceOnDiaObject`: table of forced photometry in all PVI and all difference images at the location of all `DiaObjects`\n",
    "\n",
    "The DP0.2-era `DiaObject` table contains summary statistics for each object, which we can use to identify candidate variables.\n",
    "In the future, this table will contain additional variability measurements as described in the Data Management Tech Note (DMTN) \"Review of Timeseries Features\", <a href=\"https://dmtn-118.lsst.io\">DMTN-118</a>.\n",
    "\n",
    "For simplicity, the RR Lyrae candidate search we will conduct in this section will be a cone search centered at the coordinates of the known RR Lyrae from the previous section (RA, Dec = 62.1479031, -35.799138) with a 5 degree search radius."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23fc9bc8-a06d-4057-8f7d-4c232af5be4e",
   "metadata": {},
   "source": [
    "### 3.1 Query DiaObject and ForcedSourceOnDiaObject to Retrieve Candidates and Photometry\n",
    "\n",
    "The `DiaObject` table has several parameters that characterize the variability of `DiaSource` flux measurements (see the <a href=\"https://dm.lsst.org/sdm_schemas/browser/dp02.html#DiaObject\">DP0.2 Schema Browser</a> for a full list).\n",
    "These include the min/max fluxes, and various statistical measures of the scatter in flux measurements.\n",
    "A particularly useful statistic for variable stars is the \"StetsonJ\" statistic, based on <a href=\"https://ui.adsabs.harvard.edu/abs/1996PASP..108..851S/abstract\">this 1996 paper</a>.\n",
    "Variable stars have flux that varies more than expected based on the measurement uncertainties.\n",
    "High values of the StetsonJ index should identify objects whose fluxes varied much more than expected -- i.e., variable stars.\n",
    "Here we use StetsonJ>20, but users are encouraged to explore this and other statistics.\n",
    "\n",
    "Here are all of the criteria we apply to our selection from `DiaObject`:\n",
    "\n",
    " 1. g-band measurements only\n",
    " 2. coordinates within 5 degrees of our chosen RA, Dec position\n",
    " 3. `gTOTFluxSigma`/`gTOTFluxMean` > 0.25 -- the scatter in measured fluxes is larger than 25% relative to the mean\n",
    " 4. `gTOTFluxSigma`/`gTOTFluxMean` < 1.25 -- the scatter in measured fluxes is no larger than 125% relative to the mean\n",
    " 5. 18 < `diao_gmag` < 23 -- mean g magnitude between 18-23\n",
    " 6. `gPSFluxNdata` > 30 -- at least 30 observations in g band\n",
    " 7. `gPSFluxStetsonJ` > 20 -- StetsonJ index greater than 20 \n",
    "\n",
    "In the ADQL query below, we implement all of these criteria and perform a table JOIN between `DiaObject` and `ForcedSourceOnDiaObject`,\n",
    "the latter of which contains point-source forced-photometry measurements on processed visit images and difference images.\n",
    "Since we are looking for **variable** candidates, which do appear in the template image,\n",
    "we use the forced-photometry measurements from the **non-substracted** processed visit images (PVIs), also called the \"direct\" images.\n",
    "\n",
    "We must perform one more table JOIN using the `ccdVisitId` from `FourcedSourceOnDiaObject` in order to get the MJD of each visit\n",
    "from the `CcdVisit` table.\n",
    "Note that the order of the table JOIN functions is important, and a different order may change the result of the search. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c741ca88-ffd0-4758-ae59-55b1e03457d0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT diao.diaObjectId, diao.ra, diao.decl, diao.nDiaSources, \"\\\n",
    "        \"diao.gPSFluxMin, diao.gPSFluxMax, diao.gPSFluxMean, diao.gPSFluxSigma, \"\\\n",
    "        \"diao.gPSFluxMAD, diao.gPSFluxChi2, diao.gPSFluxNdata, diao.gPSFluxSkew, \"\\\n",
    "        \"diao.gPSFluxStetsonJ, diao.gPSFluxPercentile05, diao.gPSFluxPercentile25, \"\\\n",
    "        \"diao.gPSFluxPercentile50, diao.gPSFluxPercentile75, diao.gPSFluxPercentile95, \"\\\n",
    "        \"diao.gTOTFluxMean, diao.gTOTFluxSigma, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(diao.gTOTFluxMean) as diao_gmag, \"\\\n",
    "        \"fsodo.diaObjectId, \"\\\n",
    "        \"fsodo.ccdVisitId, fsodo.band, fsodo.psfFlux, fsodo.psfFluxErr, \"\\\n",
    "        \"fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, \"\\\n",
    "        \"cv.expMidptMJD, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(fsodo.psfFlux) as fsodo_gmag \"\\\n",
    "        \"FROM dp02_dc2_catalogs.DiaObject as diao \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.ForcedSourceOnDiaObject as fsodo \"\\\n",
    "        \"ON fsodo.diaObjectId = diao.diaObjectId \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.CcdVisit as cv ON cv.ccdVisitId = fsodo.ccdVisitId \"\\\n",
    "        \"WHERE diao.gTOTFluxSigma/diao.gTOTFluxMean > 0.25 \"\\\n",
    "        \"AND diao.gTOTFluxSigma/diao.gTOTFluxMean < 1.25 \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(diao.gTOTFluxMean) > 18 \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(diao.gTOTFluxMean) < 23 \"\\\n",
    "        \"AND diao.gPSFluxNdata > 30 \"\\\n",
    "        \"AND diao.gPSFluxStetsonJ > 20 \"\\\n",
    "        \"AND fsodo.band = 'g' \"\\\n",
    "        \"AND CONTAINS(POINT('ICRS', diao.ra, diao.decl), \"\\\n",
    "        \"CIRCLE('ICRS',\"+str(ra_known_rrl)+\", \"+str(dec_known_rrl)+\", 5)) = 1 \"\n",
    "\n",
    "results = service.search(query)\n",
    "fsodo_sources = results.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32d87585-b680-409f-8aa4-fa12fc9a9357",
   "metadata": {},
   "source": [
    "Uncomment the following cell to view the results table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38012841-131c-44bf-bd01-e9e0ec937dd7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#fsodo_sources"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "debe28ea-a0cd-42cb-a256-6450ab7afb72",
   "metadata": {},
   "source": [
    "We can also check out the number of unique `DiaObjects` and the total number of measurements from the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bf2124b-e87e-4532-a862-d9565994ea51",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(len(fsodo_sources), len(np.unique(fsodo_sources['diaObjectId'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "765b8fe9-6a2a-40fe-b6cd-6374b952cb07",
   "metadata": {},
   "source": [
    "The statement and its output above tells us that the query returned 14 unique `DiaObjects`, with a total of 702 measurements in `ForcedSourceOnDiaObject`.\n",
    "\n",
    "In the following cell, we show the unique `DiaObjects` and their coordinates.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc46254b-99e0-40a4-bb2b-d5b5106031c5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "unique(fsodo_sources,keys = 'diaObjectId')['diaObjectId','ra','decl']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a64c0b9-5a7d-449b-bb6e-5efb0f026e93",
   "metadata": {},
   "source": [
    "**Important Note:** Due to a likely issue with source association, the same variable can be identified with multiple `DiaObject` IDs. As can be seen above, this has occurred in this search, where the `DiaObjects` with IDs 1567762843720221100 and 1568246628836442456 are associated with the same candidate. Scientists on the Rubin Community Science team (CST) are currently investigating this issue."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fecc4941-6812-4251-a15a-0ba41a9e3850",
   "metadata": {},
   "source": [
    "### 3.2 Run the Lomb-Scargle Frequency Algorithm on All Candidates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20949b11-873f-4be9-bd77-ce7d0af03ae3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "objids = np.unique(fsodo_sources['diaObjectId'])\n",
    "LSpower = []\n",
    "LSfreq = []\n",
    "\n",
    "minfreq = 1 / (1.25*u.d)\n",
    "maxfreq = 1 / (0.05*u.d)\n",
    "for objid in objids:\n",
    "    findobj = (fsodo_sources['diaObjectId'] == objid)\n",
    "    obj_mjd_days = np.array(fsodo_sources[findobj]['expMidptMJD']) * u.day\n",
    "    obj_mags = np.array(fsodo_sources[findobj]['fsodo_gmag'])\n",
    "    obj_frequency, obj_power =\\\n",
    "        LombScargle(obj_mjd_days, obj_mags).autopower(minimum_frequency=minfreq,\n",
    "                                                      maximum_frequency=maxfreq)\n",
    "    max_power = np.argmax(obj_power)\n",
    "    LSpower.append(obj_power[max_power])\n",
    "    LSfreq.append(obj_frequency[max_power].value)\n",
    "\n",
    "LSpower = np.array(LSpower)\n",
    "LSfreq = np.array(LSfreq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cfbde2e-96a1-4bd6-ae60-25a2d91db530",
   "metadata": {},
   "source": [
    "### 3.3 Plot Lightcurves of Several Candidates\n",
    "\n",
    "Let's examine several candidates. To see if they look like RR Lyrae (and whether we identified the correct period), plot phased lightcurves based on the frequencies from the Lomb-Scargle analysis.\n",
    "\n",
    "Specifically, we'll be looking at the candidates with `DIAObject` IDs 1567428592185376787, 1567762843720221100, and 1651589610221862935."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe8e3b0b-6090-48cd-9d52-f2d34bd2e963",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "objid_index = [0,1,7]\n",
    "\n",
    "for i in objid_index:\n",
    "    fig = plt.figure(figsize=(6, 4))\n",
    "    findobj = (fsodo_sources['diaObjectId'] == objids[i])\n",
    "    obj_mjd_days = np.array(fsodo_sources[findobj]['expMidptMJD']) * u.day\n",
    "    obj_mags = np.array(fsodo_sources[findobj]['fsodo_gmag'])\n",
    "    obj_period = 1 / LSfreq[i]\n",
    "    print('DiaObjectID: ', objids[i],\n",
    "          '\\n RA: ', fsodo_sources[findobj]['ra'][0],\n",
    "          '\\n Dec: ', fsodo_sources[findobj]['decl'][0],\n",
    "          '\\n period: ', obj_period)\n",
    "\n",
    "    t0 = 0.0\n",
    "    obj_mjd_norm = (obj_mjd_days.value - t0) / obj_period\n",
    "    obj_phase = np.mod(obj_mjd_norm, 1.0)\n",
    "\n",
    "    plt.plot(obj_phase, obj_mags-np.nanmean(obj_mags), plot_filter_symbols['g'],\n",
    "             color=plot_filter_colors['g'], ms=8)\n",
    "    plt.hlines(0, 0, 1, linestyle=':', color='Gray')\n",
    "\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.xlabel('phase')\n",
    "    plt.ylabel('mag - mean')\n",
    "    plt.minorticks_on()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec04b81f-c70d-4887-b934-683d3f69543a",
   "metadata": {
    "tags": []
   },
   "source": [
    "Hooray! A few of these look like actual RR Lyrae variables! Some of the others either aren't RR Lyrae, or have mis-identified periods.\n",
    "\n",
    "Note that the lightcurve with `DiaObjectId` 1651589610221862935 is the same as the one we target in Section 2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83256b00-9b67-477a-b642-cb8f065a7e07",
   "metadata": {},
   "source": [
    "### 3.4 Cautionary tale: comparing lightcurves from ForcedSourceOnDiaObject vs. DiaSource\n",
    "\n",
    "**What if we used `DiaSource` instead of `ForcedSourceOnDiaObject`?** \n",
    "\n",
    "In this section we demonstrate the results when we perform a table JOIN between `DiaObject` and `DiaSource` instead of `ForcedSourceOnDiaObject`.\n",
    "Using `DiaSource` will only provide photometry measurements for visits in which the variable star is detected in the \n",
    "difference image (the detection threshold is a signal-to-noise ratio greater than 5).\n",
    "Variable stars will not be detected in a difference image when they are nearly the same brightness as they are in the template image,\n",
    "leaving a very small flux residual in the difference image.\n",
    "\n",
    "As demonstrated below, using `DiaSource` instead of `ForcedSourceOnDiaObject` can leave some points out of the\n",
    "lightcurve and impact any derived parameters like the period.\n",
    "\n",
    "First, we perform a similar search query as Sec. 3.1 except we use a table JOIN with `DiaSource`.\n",
    "Note that `totFlux` from `DiaSource` provides the forced PSF-flux measurement on the direct image, not the difference image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c648a038-3a21-4839-8464-ac6324022800",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "query = \"SELECT diao.diaObjectId, diao.ra, diao.decl, diao.nDiaSources, \"\\\n",
    "        \"diao.gPSFluxMin, diao.gPSFluxMax, diao.gPSFluxMean, diao.gPSFluxSigma, \"\\\n",
    "        \"diao.gPSFluxMAD, diao.gPSFluxChi2, diao.gPSFluxNdata, diao.gPSFluxSkew, \"\\\n",
    "        \"diao.gPSFluxStetsonJ, diao.gPSFluxPercentile05, diao.gPSFluxPercentile25, \"\\\n",
    "        \"diao.gPSFluxPercentile50, diao.gPSFluxPercentile75, diao.gPSFluxPercentile95, \"\\\n",
    "        \"diao.gTOTFluxMean, diao.gTOTFluxSigma, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(diao.gTOTFluxMean) as diao_gmag, \"\\\n",
    "        \"dias.diaSourceId, dias.diaObjectId, dias.midPointTai, \"\\\n",
    "        \"dias.ccdVisitId, dias.filterName, dias.psFlux, dias.psFluxErr, \"\\\n",
    "        \"dias.totFlux, dias.totFluxErr, \"\\\n",
    "        \"scisql_nanojanskyToAbMag(dias.totFlux) as dias_gmag \"\\\n",
    "        \"FROM dp02_dc2_catalogs.DiaObject as diao \"\\\n",
    "        \"JOIN dp02_dc2_catalogs.DiaSource as dias \"\\\n",
    "        \"ON dias.diaObjectId = diao.diaObjectId \"\\\n",
    "        \"WHERE diao.gTOTFluxSigma/diao.gTOTFluxMean > 0.25 \"\\\n",
    "        \"AND diao.gTOTFluxSigma/diao.gTOTFluxMean < 1.25 \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(diao.gTOTFluxMean) > 18 \"\\\n",
    "        \"AND scisql_nanojanskyToAbMag(diao.gTOTFluxMean) < 23 \"\\\n",
    "        \"AND diao.gPSFluxNdata > 30 \"\\\n",
    "        \"AND diao.gPSFluxStetsonJ > 20 \"\\\n",
    "        \"AND dias.filterName = 'g' \"\\\n",
    "        \"AND CONTAINS(POINT('ICRS', diao.ra, diao.decl), \"\\\n",
    "        \"CIRCLE('ICRS',\"+str(ra_known_rrl)+\", \"+str(dec_known_rrl)+\", 5)) = 1 \"\n",
    "\n",
    "results = service.search(query)\n",
    "dia_sources = results.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "528dcf02-5baf-4443-8508-66fb2ecbfea7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(len(dia_sources), len(np.unique(dia_sources['diaObjectId'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12b6bbf3-e099-479a-a04a-5ebd7f101628",
   "metadata": {},
   "source": [
    "The query returned 14 unique `DiaObjects`, the same as the query from Section 3.1.\n",
    "However, the query with `DiaSource` resulted in a total of 549 measurements, whereas the query with `ForcedSourceOnDiaObject` resulted in a total of 702 measurements.\n",
    "\n",
    "We now repeat the same Lomb-Scargle analysis with the `DiaSource` results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a84134c8-0e50-465c-8c3f-c78dfe1a02af",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "objids_dia = np.unique(dia_sources['diaObjectId'])\n",
    "LSpower_dia = []\n",
    "LSfreq_dia = []\n",
    "\n",
    "minfreq = 1 / (1.25*u.d)\n",
    "maxfreq = 1 / (0.05*u.d)\n",
    "for objid in objids:\n",
    "    findobj = (dia_sources['diaObjectId'] == objid)\n",
    "    obj_mjd_days = np.array(dia_sources[findobj]['midPointTai']) * u.day\n",
    "    obj_mags = np.array(dia_sources[findobj]['dias_gmag'])\n",
    "    obj_frequency, obj_power =\\\n",
    "        LombScargle(obj_mjd_days, obj_mags).autopower(minimum_frequency=minfreq,\n",
    "                                                      maximum_frequency=maxfreq)\n",
    "    max_power = np.argmax(obj_power)\n",
    "    LSpower_dia.append(obj_power[max_power])\n",
    "    LSfreq_dia.append(obj_frequency[max_power].value)\n",
    "\n",
    "LSpower_dia = np.array(LSpower_dia)\n",
    "LSfreq_dia = np.array(LSfreq_dia)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "285d4e4a-336c-4c99-81cc-39e2156d644a",
   "metadata": {},
   "source": [
    "We now plot the resulting phased lightcurves from the `DiaSource` catalog (open markers) and the `ForcedSourceOnDiaObject` catalog (filled markers) of the same 3 selected candidates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dd74942-a250-4292-9adb-4a718c5312a2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "for i in objid_index:\n",
    "    fig = plt.figure(figsize=(6, 4))    \n",
    "    findobj_fsodo = (fsodo_sources['diaObjectId'] == objids[i])\n",
    "    findobj_dia = (dia_sources['diaObjectId'] == objids[i])\n",
    "    obj_mjd_days_dia = np.array(dia_sources[findobj_dia]['midPointTai']) * u.day\n",
    "    obj_mags_dia = np.array(dia_sources[findobj_dia]['dias_gmag'])\n",
    "    obj_mjd_days_fsodo = np.array(fsodo_sources[findobj_fsodo]['expMidptMJD']) * u.day\n",
    "    obj_mags_fsodo = np.array(fsodo_sources[findobj_fsodo]['fsodo_gmag'])\n",
    "    obj_period_fsodo = 1 / LSfreq[i]\n",
    "    obj_period_dia = 1 / LSfreq_dia[i]\n",
    "    print('DiaObjectID: ', objids[i],\n",
    "          '\\n RA: ', fsodo_sources[findobj_fsodo]['ra'][0],\n",
    "          '\\n Dec: ', fsodo_sources[findobj_fsodo]['decl'][0],\n",
    "          '\\n period (DiaSource): ', obj_period_dia,\n",
    "          '\\n period (ForcedSourceOnDiaObject)',obj_period_fsodo)\n",
    "    \n",
    "    t0 = 0.0\n",
    "    obj_mjd_norm_fsodo = (obj_mjd_days_fsodo.value - t0) / obj_period_fsodo\n",
    "    obj_mjd_norm_dia = (obj_mjd_days_dia.value - t0) / obj_period_dia\n",
    "    obj_phase_dia = np.mod(obj_mjd_norm_dia, 1.0)\n",
    "    obj_phase_fsodo = np.mod(obj_mjd_norm_fsodo, 1.0)\n",
    "\n",
    "    plt.plot(obj_phase_dia, obj_mags_dia-np.nanmean(obj_mags_dia), plot_filter_symbols['g'],\n",
    "             markerfacecolor = 'none', mec=plot_filter_colors['g'], ms=10, \n",
    "             label='DIA (N=%i)' % len(obj_phase_dia))\n",
    "    plt.plot(obj_phase_fsodo, obj_mags_fsodo-np.nanmean(obj_mags_fsodo), plot_filter_symbols['g'],\n",
    "             color=plot_filter_colors['g'], ms=6, label='FSODO (N=%i)' % len(obj_phase_fsodo))\n",
    "    plt.hlines(0, 0, 1, linestyle=':', color='Gray')\n",
    "\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.legend(loc='upper right', labelspacing=0.2, handletextpad=0.2, fontsize=11)\n",
    "    plt.xlabel('phase')\n",
    "    plt.ylabel('mag - mean')\n",
    "    plt.minorticks_on()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "780089e3-b3e7-403f-909a-fd3491d016e2",
   "metadata": {},
   "source": [
    "Note that the number of measurements is always greater from the `ForcedSourceOnDiaObject` results.\n",
    "\n",
    "The plots demonstrate the possible deviations in this analysis using the different catalogs.\n",
    "Not only can the periods differ from using the different catalogs (e.g. 1567762843720221100), but the phase may also be offset due to additional measurements provided by the `ForcedSourceOnDiaObject` catalog. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "798f11b3-af08-4c60-91b7-6dc244430cd0",
   "metadata": {},
   "source": [
    "Lastly, we can estimate the fraction of candidates that will have different periods between the `ForcedSourceOnDiaObject` and `DiaSource` catalogs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "331f9016-f9f5-4727-b89f-5943cf0387a2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "1-(LSfreq_dia-LSfreq).tolist().count(0.0)/len(LSfreq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dcc9f22-77b8-454a-9862-3b659538736a",
   "metadata": {},
   "source": [
    "We see that **~70%** of the identified candidates exhibit different periods from both catalogs.\n",
    "\n",
    "Plot a histogram of the absolute difference in the periods measured from `ForcedSourceOnDiaObject` and `DiaSource`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30f20541-aebc-42bc-8e84-e7ab5a000134",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "period_difference = np.abs((1.0 / LSfreq_dia) - (1.0 / LSfreq))\n",
    "\n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "plt.hist(period_difference, bins=20)\n",
    "plt.axvline(1.0 / 24.0, ls='--', color='red', label='1 hour')\n",
    "plt.legend(loc='upper right')\n",
    "plt.xlabel('absolute difference in period [days]', fontsize=12)\n",
    "plt.ylabel('number of candidate RR Lyrae', fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f797101-397f-4332-8cdc-2d355956b729",
   "metadata": {},
   "source": [
    "Above, we can see that only 5 of the 14 candidate RR Lyrae (which have not been confirmed to _be_ RR Lyrae) have\n",
    "periods that change by less than 1 hour when the incomplete `DiaSource` photometry is used instead of `ForcedSourceOnDiaObject`.\n",
    "\n",
    "> Note: because of the source association issue pointed out at the end of Section 3.1, there are two `DiaObject` IDs from\n",
    "the search that correspond to the same candidate.\n",
    "This means that the fraction of unique RR Lyrae stars with different periods should actually be _greater_ than ~70%."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8a6fd31-0176-4250-8831-9aeb7c341c23",
   "metadata": {},
   "source": [
    "## 4. Exercises for the learner\n",
    "\n",
    "There are many further explorations one could try as extensions of this notebook. Some examples:\n",
    "\n",
    "1. Experiment with ways to identify periodic signals in the data.\n",
    "2. Refine the selection criteria for finding candidates.\n",
    "3. Compare results from candidate search to the true star lightcurves from the DESC truth tables.\n",
    "\n",
    "...and many more. Have fun!"
   ]
  }
 ],
 "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.11.7"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}