variable_lightcurves_vS2D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "749b0ddf",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=170 style=\"padding: 10px\"> \n",
    "<b>Little Demo: Variable Light curves</b> <br>\n",
    "Contact author(s): Ryan Lau <br>\n",
    "Last verified to run: 2023-04-27 <br>\n",
    "LSST Science Pipelines version: Weekly 2023_07 <br>\n",
    "Container Size: medium"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff071fc1-aa10-46b0-a5e4-cc6d3e08eec1",
   "metadata": {},
   "source": [
    "Let's plot a light curve and do a Lomb-Scargle periodogram analysis on an RR Lyrae variable. This is a truncated version of the DP0.2 tutorial notebook 07b, \"Variable Stars in DP0.2\"\n",
    "\n",
    "## 1. Setup and Package Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c128f938-b3ca-43a4-818f-8e1c4143e402",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# General\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('tableau-colorblind10')\n",
    "\n",
    "# Astropy\n",
    "import astropy.units as u\n",
    "from astropy.timeseries import LombScargle\n",
    "\n",
    "# LSST package for TAP queries\n",
    "from lsst.rsp import get_tap_service"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26912cc0-ae49-4b05-a4d9-89457b9a2c85",
   "metadata": {},
   "source": [
    "Set up some plotting defaults so plots will look nice:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dfa2640-9411-4972-9d22-b3082b2b178d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\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",
    "          }\n",
    "\n",
    "plt.rcParams.update(params)\n",
    "\n",
    "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": "3c77197a-75d3-4cdc-8acc-6eb81186edeb",
   "metadata": {},
   "source": [
    "Start the TAP service, which we will use for all data retrieval in this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8b5c0d2-6f26-4579-8025-2f17ca35b25c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "service = get_tap_service()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69599956-3fbe-4851-a9b3-62b80a7d922a",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Create a Lightcurve for a Known RR Lyrae Variable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72433a32-5667-41cb-8e50-b3975403614a",
   "metadata": {},
   "source": [
    "Set the object ID of a known RR Lyrae variable in the `ForcedSource` table and run a search to get the forced source photometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee6e9cac-2b1d-4a48-910a-691d404295bc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sel_objid = 1651589610221899038\n",
    "\n",
    "# Define the query\n",
    "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",
    "\n",
    "results = service.search(query)\n",
    "srcs = results.to_table()\n",
    "del query\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f230c331-11a3-43db-b10e-429c8fb94b0c",
   "metadata": {},
   "source": [
    "Uncomment and execute the following cell to see what this table looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "840b1bbf-d1be-4b42-bf81-187b36b20271",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#srcs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "680bcd70-9cd4-4d4b-b215-0d875b460379",
   "metadata": {
    "tags": []
   },
   "source": [
    "Plot the lightcurve of the variable in all (ugrizy) filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09d47b1a-a0c7-4701-bd9d-5ddf326bd0db",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "pick = {}\n",
    "for filter in plot_filter_labels:\n",
    "    pick[filter] = (srcs['band'] == filter)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9,6))\n",
    "    \n",
    "for filter in plot_filter_labels:\n",
    "    ax.plot(srcs[pick[filter]]['expMidptMJD'], srcs[pick[filter]]['psfMag'],\n",
    "    plot_filter_symbols[filter], ms=10, mew=2, mec=plot_filter_colors[filter],\n",
    "    alpha=0.5, color='none', label=filter)\n",
    "\n",
    "ax.set_xlabel('Modified Julian Date')\n",
    "ax.set_ylabel('Magnitude')\n",
    "ax.set_title('Light Curve of RR Lyrae Variable')\n",
    "ax.legend(loc='center left', ncol=2)\n",
    "ax.invert_yaxis()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21672fd6-4d5a-4006-9c48-1dd257b389f9",
   "metadata": {},
   "source": [
    "## 3 Create a Phased Lightcurve"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d84682c-3242-414f-bd20-c863415cb23b",
   "metadata": {},
   "source": [
    "First, we will extract a periodogram. We start off by defining the parameters of the Lomb-Scargle periodogram analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "097aecbb-243e-48f4-b2c0-1c4170ea8359",
   "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'])\n",
    "    \n",
    "    \n",
    "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": "adc21cfd-3e76-4620-8c50-4819393b956c",
   "metadata": {},
   "source": [
    "Now we run the Lomb-Scargle algorithm the photometry of each filter. We will set only the min/max frequencies, and otherwise use the \"autopower\" method from Lomb-Scargle. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a6a3e94-7e21-40c0-995c-e7efb78d63ee",
   "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",
    "    \n",
    "all_peak_freqs = []\n",
    "\n",
    "for filter in plot_filter_labels:\n",
    "    # find the index with maximum power (= peakbin)\n",
    "    peakbin = np.argmax(power[filter])\n",
    "\n",
    "    # Store the frequency corresponding to the peak power in each filter (band)\n",
    "    all_peak_freqs.append(frequency[filter][peakbin].value)\n",
    "\n",
    "# Convert the frequencies from a list to and array:\n",
    "all_peak_freqs = np.array(all_peak_freqs)\n",
    "\n",
    "# Calculate the mean of the \"best-fit\" frequencies:\n",
    "mean_peak_freq = np.mean(all_peak_freqs)\n",
    "best_period = 1/mean_peak_freq\n",
    "\n",
    "print('Mean frequency:', mean_peak_freq)\n",
    "print('Mean period:', best_period, ' days')\n",
    "print('\\nugrizy frequency results:\\n', all_peak_freqs)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01289bbe-ca34-44d1-ba46-e0e9d2cbd7d2",
   "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34f65e83-dcdd-4e3d-ace3-bc4a4eb0bf5a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "mjd_norm = {}\n",
    "phase = {}\n",
    "\n",
    "# Number of elapsed periods since the first measurement:\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",
    "\n",
    "    # Calculate phase by using the modulus function:\n",
    "    phase[filter] = np.mod(mjd_norm[filter], 1.0)\n",
    "    \n",
    "fig = plt.figure(figsize=(7, 10))     # create a matplotlib.pyplot figure\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 (plus integer offset)')\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d7ef122-4978-40af-86c2-8aae4c4d147d",
   "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.10.9"
  },
  "toc-autonumbering": false
 },
 "nbformat": 4,
 "nbformat_minor": 5
}