notebooks/MAST/PanSTARRS/PS1_DR2_TAP/PS1_DR2_TAP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MAST Table Access Protocol PanSTARRS 1 DR2 Demo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br> This tutorial demonstrates how to use astroquery to access PanSTARRS 1 Data Release 2 via a Virtual Observatory standard Table Access Protocol (TAP) service at MAST, and work with the resultant data. It relies on Python 3 and astroquery, as well as some other common scientific packages.\n",
    "\n",
    "***\n",
    "### Table of Contents\n",
    "\n",
    "1. [TAP Service Introduction](#TAP-Service-Introduction)\n",
    "2. [Imports](#Imports)\n",
    "3. [Connecting to a TAP Service](#Connecting-to-a-TAP-Service)\n",
    "4. [Use Cases](#Use-Cases)\n",
    "5. [Additional Resources](#Additional-Resources)\n",
    "6. [About This Notebook](#About-this-Notebook)\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TAP Service Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Table Access Protocol (TAP) services allow more direct and flexible access to astronomical data than the simpler types of IVOA standard data services. Queries are built with the SQL-like Astronomical Data Query Language (ADQL), and can include geographic / spatial queries as well as filtering on other characteristics of the data. This also allows the user fine-grained control over the returned columns, unlike the fixed set of coumns returned from cone, image, and spectral services.\n",
    "\n",
    "For this example, we'll be using the astropy affiliated PyVO client, which is interoperable with other valid TAP services, including those at MAST. PyVO documentation is available at ReadTheDocs: https://pyvo.readthedocs.io\n",
    "\n",
    "We'll be using PyVO to call the TAP service at MAST serving PanSTARRS 1 Data Release 2, now with individual detection information. The schema is described within the service, and we'll show how to inspect it. The schema is also the same as the one available via the CasJobs interface, with an additional view added for the most common positional queries. CasJobs has its own copy of the schema documentation, which can be accessed through its own site: http://mastweb.stsci.edu/ps1casjobs/\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the pyvo library as our client to the data service.\n",
    "import pyvo as vo\n",
    "\n",
    "# For resolving objects with tools from MAST\n",
    "from astroquery.mast import Mast\n",
    "\n",
    "# For handling ordinary astropy Tables in responses\n",
    "from astropy.table import Table\n",
    "\n",
    "# For displaying and manipulating some types of results\n",
    "%matplotlib inline\n",
    "import requests\n",
    "import astropy\n",
    "import numpy as np\n",
    "import pylab\n",
    "import time\n",
    "import json\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# For calling the object name resolver. Note these are Python 3 dependencies\n",
    "#import sys\n",
    "#import http.client as httplib \n",
    "#from urllib.parse import quote as urlencode\n",
    "#from urllib.request import urlretrieve\n",
    "\n",
    "# To allow display tweaks for wider response tables\n",
    "from IPython.core.display import display\n",
    "from IPython.core.display import HTML\n",
    "\n",
    "# suppress unimportant unit warnings from many TAP services\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", module=\"astropy.io.votable.*\")\n",
    "warnings.filterwarnings(\"ignore\", module=\"pyvo.utils.xml.elements\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "## Connecting to a TAP Service"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The PyVO library is able to connect to any TAP service, given the \"base\" URL as noted in metadata registry resources describing the service. This is the URL for the PanSTARRS 1 DR2 TAP service."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TAP_service = vo.dal.TAPService(\"https://vao.stsci.edu/PS1DR2/tapservice.aspx\")\n",
    "TAP_service.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### List available tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TAP_tables = TAP_service.tables\n",
    "for tablename in TAP_tables.keys():\n",
    "    if not \"tap_schema\" in tablename:  \n",
    "        TAP_tables[tablename].describe()\n",
    "        print(\"Columns={}\".format(sorted([k.name for k in TAP_tables[tablename].columns ])))\n",
    "        print(\"----\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Use Cases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Positional Query\n",
    "\n",
    "This searches the mean object catalog for objects within .2 degrees of M87 (RA=187.706, Dec=12.391 in degrees). The view used contains information from the [ObjectThin](https://outerspace.stsci.edu/x/W4Oc) table (which has information on object positions and the number of available measurements) and the [MeanObject](https://outerspace.stsci.edu/x/WYOc) table (which has information on photometry averaged over the multiple epochs of observation).\n",
    "\n",
    "Note that the results are restricted to objects with `nDetections>1`, where `nDetections` is the total number of times the object was detected on the single-epoch images in any filter at any time.  Objects with `nDetections=1` tend to be  artifacts, so this is a quick way to eliminate most spurious objects from the catalog.\n",
    "\n",
    "This query runs in TAP's asynchronous mode, which is a queued batch mode with some overhead and longer timeouts, useful for big catalogs like PanSTARRS. It may not be necessary for all queries to PS1 DR2, but the PyVO client can automatically handle the additional processing required over synchronous mode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = TAP_service.run_async(\"\"\"\n",
    "SELECT objID, RAMean, DecMean, nDetections, ng, nr, ni, nz, ny, gMeanPSFMag, rMeanPSFMag, iMeanPSFMag, zMeanPSFMag, yMeanPSFMag\n",
    "FROM dbo.MeanObjectView\n",
    "WHERE\n",
    "CONTAINS(POINT('ICRS', RAMean, DecMean),CIRCLE('ICRS',187.706,12.391,.2))=1\n",
    "AND nDetections > 1\n",
    "  \"\"\")\n",
    "TAP_results = job.to_table()\n",
    "TAP_results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get DR2 light curve for RR Lyrae star KQ UMa\n",
    "\n",
    "This time we start with the object name, use the MAST name resolver (which relies on Simbad and NED) to convert the name to RA and Dec, and then query the PS1 DR2 mean object catalog at that position.  Then we run a spatial query to TAP using those coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objname = 'KQ UMa'\n",
    "coords = Mast.resolve_object(objname)\n",
    "ra,dec = coords.ra.value,coords.dec.value\n",
    "radius = 1.0/3600.0 # radius = 1 arcsec\n",
    "\n",
    "query = \"\"\"\n",
    "SELECT objID, RAMean, DecMean, nDetections, ng, nr, ni, nz, ny, gMeanPSFMag, \n",
    "    rMeanPSFMag, iMeanPSFMag, zMeanPSFMag, yMeanPSFMag\n",
    "FROM dbo.MeanObjectView\n",
    "WHERE\n",
    "CONTAINS(POINT('ICRS', RAMean, DecMean),CIRCLE('ICRS',{},{},{}))=1\n",
    "AND nDetections > 1\n",
    "\"\"\".format(ra,dec,radius)\n",
    "print(query)\n",
    "\n",
    "job = TAP_service.run_async(query)\n",
    "TAP_results = job.to_table()\n",
    "TAP_results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get Repeated Detection Information\n",
    "\n",
    "Extract all the objects with the same object ID from the [Detection](https://outerspace.stsci.edu/x/b4Oc) table, which contains all the individual measurements for this source. The results are joined to the [Filter](https://outerspace.stsci.edu/x/nIOc) table to convert the filter numbers to names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objid = TAP_results['objID'][0]\n",
    "query = \"\"\"\n",
    "SELECT\n",
    "    objID, detectID, Detection.filterID as filterID, Filter.filterType, obsTime, ra, dec,\n",
    "    psfFlux, psfFluxErr, psfMajorFWHM, psfMinorFWHM, psfQfPerfect, \n",
    "    apFlux, apFluxErr, infoFlag, infoFlag2, infoFlag3\n",
    "FROM Detection\n",
    "NATURAL JOIN Filter\n",
    "WHERE objID={}\n",
    "ORDER BY filterID, obsTime\n",
    "\"\"\".format(objid)\n",
    "print(query)\n",
    "\n",
    "job = TAP_service.run_async(query)\n",
    "detection_TAP_results = job.to_table()\n",
    "detection_TAP_results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the light curves\n",
    "\n",
    "The `psfFlux` values from the Detection table are converted from Janskys to AB magnitudes.  Measurements in the 5 different filters are plotted separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert flux in Jy to magnitudes\n",
    "t = detection_TAP_results['obsTime']\n",
    "mag = -2.5*np.log10(detection_TAP_results['psfFlux']) + 8.90\n",
    "xlim = np.array([t.min(),t.max()])\n",
    "xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])\n",
    "\n",
    "pylab.rcParams.update({'font.size': 14})\n",
    "pylab.figure(1,(10,10))\n",
    "\n",
    "#detection_TAP_results['filterType'] is a byte string, compare accordingly:\n",
    "for i, filter in enumerate(['g','r','i','z','y']):\n",
    "    \n",
    "    pylab.subplot(511+i)\n",
    "    w = np.where(detection_TAP_results['filterType'] == filter)  \n",
    "    pylab.plot(t[w],mag[w],'-o')\n",
    "    pylab.ylabel(f'{filter} [mag]')\n",
    "    pylab.xlim(xlim)\n",
    "    pylab.gca().invert_yaxis()\n",
    "    if i==0:\n",
    "        pylab.title(objname)\n",
    "        \n",
    "pylab.xlabel('Time [MJD]')\n",
    "pylab.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot differences from the mean magnitudes in the initial search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert flux in Jy to magnitudes\n",
    "t = detection_TAP_results['obsTime']\n",
    "mag = -2.5*np.log10(detection_TAP_results['psfFlux']) + 8.90\n",
    "xlim = np.array([t.min(),t.max()])\n",
    "xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])\n",
    "\n",
    "pylab.rcParams.update({'font.size': 14})\n",
    "pylab.figure(1,(10,10))\n",
    "\n",
    "#detection_TAP_results['filterType'] is a byte string, compare accordingly:\n",
    "for i, filter in enumerate(['g','r','i','z','y']):\n",
    "    pylab.subplot(511+i)\n",
    "    w = np.where(detection_TAP_results['filterType']==filter)\n",
    "    magmean = TAP_results[f'{filter}MeanPSFMag'][0]\n",
    "    pylab.plot(t[w],mag[w] - magmean,'-o')\n",
    "    pylab.ylabel(f'{filter} [mag - {np.round(float(magmean), 4)}]')\n",
    "    pylab.xlim(xlim)\n",
    "    pylab.gca().invert_yaxis()\n",
    "    if i==0:\n",
    "        pylab.title(objname)\n",
    "        \n",
    "pylab.xlabel('Time [MJD]')\n",
    "pylab.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identify bad data\n",
    "\n",
    "There is one clearly bad $z$ magnitude with a very large difference.  Select the bad point and look at it in more detail.\n",
    "\n",
    "Note that indexing a table (or numpy array) with a logical expression selects just the rows where that expression is true."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detection_TAP_results[ (detection_TAP_results['filterType']=='z') & (np.abs(mag-TAP_results['zMeanPSFMag'][0]) > 2) ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From examining this table, it looks like `psfQfPerfect` is bad.  This flag is the PSF-weighted fraction of unmasked pixels in the image (see the [documentation](https://outerspace.stsci.edu/x/IoOc) for more details). Values near unity indicate good data that is not significantly affected by bad pixels.\n",
    "\n",
    "Check all the `psfQfPerfect` values for the $z$ filter to see if this value really is unusual.  The list of values below are sorted by magnitude.  The bad point is the only value with `psfQfPerfect` < 0.95."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = np.where(detection_TAP_results['filterType']=='z')\n",
    "zdtab = detection_TAP_results[w]\n",
    "zdtab['mag'] = mag[w]\n",
    "zdtab['dmag'] = zdtab['mag'] - TAP_results['zMeanPSFMag'][0]\n",
    "ii = np.argsort(-np.abs(zdtab['dmag']))\n",
    "zdtab = zdtab[ii]\n",
    "zdtab['objID','obsTime','mag','dmag','psfQfPerfect']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Repeat the plot with bad psfQfPerfect values excluded\n",
    "\n",
    "Do the plot again but exclude low psfQfPerfect values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert flux in Jy to magnitudes\n",
    "t = detection_TAP_results['obsTime']\n",
    "mag = -2.5*np.log10(detection_TAP_results['psfFlux']) + 8.90\n",
    "magmean = 0.0*mag\n",
    "for i, filter in enumerate(['g','r','i','z','y']):\n",
    "    magmean[detection_TAP_results['filterType']==filter] = TAP_results[f'{filter}MeanPSFMag'][0]\n",
    "\n",
    "dmag = mag - magmean\n",
    "dmag1 = dmag[detection_TAP_results['psfQfPerfect']>0.9]\n",
    "# fix the x and y axis ranges\n",
    "xlim = np.array([t.min(),t.max()])\n",
    "xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])\n",
    "# flip axis direction for magnitude\n",
    "ylim = np.array([dmag1.max(),dmag1.min()])\n",
    "ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])\n",
    "\n",
    "pylab.rcParams.update({'font.size': 14})\n",
    "pylab.figure(1,(10,10))\n",
    "\n",
    "for i, filter in enumerate(['g','r','i','z','y']):\n",
    "    pylab.subplot(511+i)\n",
    "    w = np.where((detection_TAP_results['filterType']==filter) & (detection_TAP_results['psfQfPerfect']>0.9))[0]\n",
    "    pylab.plot(t[w],dmag[w],'-o')\n",
    "    pylab.ylabel('{} [mag - {:.2f}]'.format(filter,magmean[w[0]]))\n",
    "    pylab.xlim(xlim)\n",
    "    pylab.ylim(ylim)\n",
    "    if i==0:\n",
    "        pylab.title(objname)\n",
    "        \n",
    "pylab.xlabel('Time [MJD]')\n",
    "pylab.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot versus the periodic phase instead of epoch\n",
    "\n",
    "Plot versus phase using known RR Lyr period from Simbad (table [J/AJ/132/1202/table4](http://vizier.u-strasbg.fr/viz-bin/VizieR-3?-source=J/AJ/132/1202/table4&-c=KQ%20UMa&-c.u=arcmin&-c.r=2&-c.eq=J2000&-c.geom=r&-out.max=50&-out.form=HTML%20Table&-oc.form=sexa))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "period = 0.48636 #days, from Simbad\n",
    "# convert flux in Jy to magnitudes\n",
    "t = (detection_TAP_results['obsTime'] % period) / period\n",
    "mag = -2.5*np.log10(detection_TAP_results['psfFlux']) + 8.90\n",
    "magmean = 0.0*mag\n",
    "for i, filter in enumerate(['g','r','i','z','y']):\n",
    "    magmean[detection_TAP_results['filterType']==filter] = TAP_results[f'{filter}MeanPSFMag'][0]\n",
    "    \n",
    "dmag = mag - magmean\n",
    "dmag1 = dmag[detection_TAP_results['psfQfPerfect']>0.9]\n",
    "# fix the x and y axis ranges\n",
    "xlim = np.array([t.min(),t.max()])\n",
    "xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])\n",
    "# flip axis direction for magnitude\n",
    "ylim = np.array([dmag1.max(),dmag1.min()])\n",
    "ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])\n",
    "\n",
    "pylab.rcParams.update({'font.size': 14})\n",
    "pylab.figure(1,(10,10))\n",
    "\n",
    "for i, filter in enumerate(['g','r','i','z','y']):\n",
    "    pylab.subplot(511+i)\n",
    "    w = np.where((detection_TAP_results['filterType']==filter) & (detection_TAP_results['psfQfPerfect']>0.9))[0]\n",
    "    w = w[np.argsort(t[w])]\n",
    "    pylab.plot(t[w],dmag[w],'-o')\n",
    "    pylab.ylabel('{} [mag - {:.2f}]'.format(filter,magmean[w[0]]))\n",
    "    pylab.xlim(xlim)\n",
    "    pylab.ylim(ylim)\n",
    "    if i==0:\n",
    "        pylab.title(objname)\n",
    "        \n",
    "pylab.xlabel('Phase')\n",
    "pylab.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repeat search using eclipsing binary KIC 2161623\n",
    "\n",
    "From [Villanova Kepler Eclipsing Binaries](http://keplerebs.villanova.edu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objname = 'KIC 2161623'\n",
    "coords = Mast.resolve_object(objname)\n",
    "ra,dec = coords.ra.value,coords.dec.value\n",
    "radius = 1.0/3600.0 # radius = 1 arcsec\n",
    "\n",
    "query = \"\"\"\n",
    "SELECT objID, RAMean, DecMean, nDetections, ng, nr, ni, nz, ny, gMeanPSFMag, rMeanPSFMag, iMeanPSFMag, zMeanPSFMag, yMeanPSFMag\n",
    "FROM dbo.MeanObjectView\n",
    "WHERE\n",
    "CONTAINS(POINT('ICRS', RAMean, DecMean),CIRCLE('ICRS',{},{},{}))=1\n",
    "AND nDetections > 1\n",
    "\"\"\".format(ra,dec,radius)\n",
    "print(query)\n",
    "\n",
    "job = TAP_service.run_async(query)\n",
    "TAP_results = job.to_table()\n",
    "TAP_results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get Repeated Detection Information\n",
    "\n",
    "This time include the `psfQfPerfect` limit directly in the database query."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objid = TAP_results['objID'][0]\n",
    "\n",
    "query = \"\"\"\n",
    "SELECT\n",
    "    objID, detectID, Detection.filterID as filterID, Filter.filterType, obsTime, ra, dec,\n",
    "    psfFlux, psfFluxErr, psfMajorFWHM, psfMinorFWHM, psfQfPerfect, \n",
    "    apFlux, apFluxErr, infoFlag, infoFlag2, infoFlag3\n",
    "FROM Detection\n",
    "NATURAL JOIN Filter\n",
    "WHERE objID={}\n",
    "AND psfQfPerfect >= 0.9\n",
    "ORDER BY filterID, obsTime\n",
    "\"\"\".format(objid)\n",
    "print(query)\n",
    "\n",
    "job = TAP_service.run_async(query)\n",
    "detection_TAP_results = job.to_table()\n",
    "\n",
    "# add magnitude and difference from mean\n",
    "detection_TAP_results['magmean'] = 0.0\n",
    "for i, filter in enumerate([b'g',b'r',b'i',b'z',b'y']):\n",
    "    detection_TAP_results['magmean'][detection_TAP_results['filterType']==filter] = TAP_results[filter.decode('ascii')+'MeanPSFMag'][0]\n",
    "detection_TAP_results['mag'] = -2.5*np.log10(detection_TAP_results['psfFlux']) + 8.90\n",
    "detection_TAP_results['dmag'] = detection_TAP_results['mag']-detection_TAP_results['magmean']\n",
    "\n",
    "detection_TAP_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = detection_TAP_results['obsTime']\n",
    "dmag = detection_TAP_results['dmag']\n",
    "xlim = np.array([t.min(),t.max()])\n",
    "xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])\n",
    "ylim = np.array([dmag.max(),dmag.min()])\n",
    "ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])\n",
    "\n",
    "pylab.rcParams.update({'font.size': 14})\n",
    "pylab.figure(1,(10,10))\n",
    "\n",
    "for i, filter in enumerate(['g','r','i','z','y']):\n",
    "    pylab.subplot(511+i)\n",
    "    w = np.where(detection_TAP_results['filterType']==filter)[0]\n",
    "    pylab.plot(t[w],dmag[w],'-o')\n",
    "    magmean = detection_TAP_results['magmean'][w[0]]\n",
    "    pylab.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))\n",
    "    pylab.xlim(xlim)\n",
    "    pylab.ylim(ylim)\n",
    "    if i==0:\n",
    "        pylab.title(objname)\n",
    "        \n",
    "pylab.xlabel('Time [MJD]')\n",
    "pylab.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot versus phase using known period\n",
    "\n",
    "Eclipsing binaries basically vary by same amount in all filters since it is a geometrical effect, so combine the data into a single light curve.  Wrap using known period and plot versus phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "period = 2.2834698\n",
    "bjd0 = 54999.599837\n",
    "t = ((detection_TAP_results['obsTime']-bjd0) % period) / period\n",
    "dmag = detection_TAP_results['dmag']\n",
    "w = np.argsort(t)\n",
    "t = t[w]\n",
    "dmag = dmag[w]\n",
    "xlim = np.array([t.min(),t.max()])\n",
    "xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])\n",
    "ylim = np.array([dmag.max(),dmag.min()])\n",
    "ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])\n",
    "\n",
    "pylab.rcParams.update({'font.size': 14})\n",
    "pylab.figure(1,(10,6))\n",
    "pylab.plot(t,dmag,'-o')\n",
    "pylab.xlim(xlim)\n",
    "pylab.ylim(ylim)\n",
    "pylab.xlabel('Phase')\n",
    "pylab.ylabel('Delta magnitude from mean [mag]')\n",
    "pylab.title(objname)\n",
    "pylab.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repeat search for another eclipsing binary KIC 8153568"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objname = 'KIC 8153568'\n",
    "coords = Mast.resolve_object(objname)\n",
    "ra,dec = coords.ra.value,coords.dec.value\n",
    "radius = 1.0/3600.0 # radius = 1 arcsec\n",
    "\n",
    "query = \"\"\"\n",
    "SELECT objID, RAMean, DecMean, nDetections, ng, nr, ni, nz, ny, gMeanPSFMag, rMeanPSFMag, iMeanPSFMag, zMeanPSFMag, yMeanPSFMag\n",
    "FROM dbo.MeanObjectView\n",
    "WHERE\n",
    "CONTAINS(POINT('ICRS', RAMean, DecMean),CIRCLE('ICRS',{},{},{}))=1\n",
    "AND nDetections > 1\n",
    "\"\"\".format(ra,dec,radius)\n",
    "print(query)\n",
    "\n",
    "job = TAP_service.run_async(query)\n",
    "TAP_results = job.to_table()\n",
    "TAP_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "objid = TAP_results['objID'][0]\n",
    "query = \"\"\"\n",
    "SELECT\n",
    "    objID, detectID, Detection.filterID as filterID, Filter.filterType, obsTime, ra, dec,\n",
    "    psfFlux, psfFluxErr, psfMajorFWHM, psfMinorFWHM, psfQfPerfect, \n",
    "    apFlux, apFluxErr, infoFlag, infoFlag2, infoFlag3\n",
    "FROM Detection\n",
    "NATURAL JOIN Filter\n",
    "WHERE objID={}\n",
    "AND psfQfPerfect >= 0.9\n",
    "ORDER BY filterID, obsTime\n",
    "\"\"\".format(objid)\n",
    "print(query)\n",
    "\n",
    "job = TAP_service.run_async(query)\n",
    "detection_TAP_results = job.to_table()\n",
    "\n",
    "# add magnitude and difference from mean\n",
    "detection_TAP_results['magmean'] = 0.0\n",
    "for i, filter in enumerate([b'g',b'r',b'i',b'z',b'y']):\n",
    "    detection_TAP_results['magmean'][detection_TAP_results['filterType']==filter] = TAP_results[filter.decode('ascii')+'MeanPSFMag'][0]\n",
    "detection_TAP_results['mag'] = -2.5*np.log10(detection_TAP_results['psfFlux']) + 8.90\n",
    "detection_TAP_results['dmag'] = detection_TAP_results['mag']-detection_TAP_results['magmean']\n",
    "\n",
    "detection_TAP_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = detection_TAP_results['obsTime']\n",
    "dmag = detection_TAP_results['dmag']\n",
    "xlim = np.array([t.min(),t.max()])\n",
    "xlim = xlim + np.array([-1,1])*0.02*(xlim[1]-xlim[0])\n",
    "ylim = np.array([dmag.max(),dmag.min()])\n",
    "ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])\n",
    "\n",
    "pylab.rcParams.update({'font.size': 14})\n",
    "pylab.figure(1,(10,10))\n",
    "\n",
    "for i, filter in enumerate(['g','r','i','z','y']):\n",
    "    pylab.subplot(511+i)\n",
    "    w = np.where(detection_TAP_results['filterType']==filter)[0]\n",
    "    pylab.plot(t[w],dmag[w],'-o')\n",
    "    magmean = detection_TAP_results['magmean'][w[0]]\n",
    "    pylab.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))\n",
    "    pylab.xlim(xlim)\n",
    "    pylab.ylim(ylim)\n",
    "    if i==0:\n",
    "        pylab.title(objname)\n",
    "        \n",
    "pylab.xlabel('Time [MJD]')\n",
    "pylab.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eclipsing binaries basically vary by same amount in all filters since it is a geometrical effect, so combine the data into a single light curve.\n",
    "\n",
    "Wrap using known period and plot versus phase.  Plot two periods of the light curve this time.\n",
    "\n",
    "This nice light curve appears to show a secondary eclipse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "period = 3.6071431\n",
    "bjd0 = 54999.289794\n",
    "t = ((detection_TAP_results['obsTime']-bjd0) % period) / period\n",
    "dmag = detection_TAP_results['dmag']\n",
    "w = np.argsort(t)\n",
    "# extend to two periods\n",
    "nw = len(w)\n",
    "w = np.append(w,w)\n",
    "t = t[w]\n",
    "# add one to second period\n",
    "t[-nw:] += 1\n",
    "dmag = dmag[w]\n",
    "xlim = [0,2.0]\n",
    "ylim = np.array([dmag.max(),dmag.min()])\n",
    "ylim = ylim + np.array([-1,1])*0.02*(ylim[1]-ylim[0])\n",
    "\n",
    "pylab.rcParams.update({'font.size': 14})\n",
    "pylab.figure(1,(12,6))\n",
    "pylab.plot(t,dmag,'-o')\n",
    "pylab.xlim(xlim)\n",
    "pylab.ylim(ylim)\n",
    "pylab.xlabel('Phase')\n",
    "pylab.ylabel('Delta magnitude from mean [mag]')\n",
    "pylab.title(objname)\n",
    "pylab.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# Additional Resources\n",
    "\n",
    "\n",
    "## Table Access Protocol \n",
    "* IVOA standard for RESTful web service access to tabular data\n",
    "* http://www.ivoa.net/documents/TAP/\n",
    "\n",
    "## PanSTARRS 1 DR 2\n",
    "* Catalog for PanSTARRS with additional Detection information\n",
    "* https://outerspace.stsci.edu/display/PANSTARRS/\n",
    "\n",
    "## Astronomical Query Data Language (2.0)\n",
    "* IVOA standard for querying astronomical data in tabular format, with geometric search support\n",
    "* http://www.ivoa.net/documents/latest/ADQL.html\n",
    "\n",
    "## PyVO\n",
    "* an affiliated package for astropy\n",
    "* find and retrieve astronomical data available from archives that support standard IVOA virtual observatory service protocols.\n",
    "* https://pyvo.readthedocs.io/en/latest/index.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## About this Notebook\n",
    "**Authors:** Rick White & Theresa Dower, STScI Archive Scientist & Software Engineer\n",
    "**Updated On:** 01/09/2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: right;\" src=\"./stsci_pri_combo_mark_horizonal_white_bkgd.png\" alt=\"stsci_pri_combo_mark_horizonal_white_bkgd\" width=\"200px\"/> "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.12"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}