{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Query Pan-STARRS DR2 catalog using MAST API\n",
    "\n",
    "The [new MAST interface](https://catalogs.mast.stsci.edu/) to the Pan-STARRS catalog supports queries to both the DR1 and DR2 PS1 catalogs.  It also has an associated [API](https://catalogs.mast.stsci.edu/docs/panstarrs.html), which is used in this script.\n",
    "\n",
    "This script shows how to query the Pan-STARRS DR2 catalog using the PS1 search API.  The examples show how to do a simple cone search, how to manipulate the table of results, and how to get a light curve from the table of detections.\n",
    "\n",
    "This notebook is available for [download](https://ps1images.stsci.edu/ps1_dr2_api.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import ascii\n",
    "from astropy.table import Table\n",
    "\n",
    "import sys\n",
    "import re\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import json\n",
    "import requests\n",
    "\n",
    "try: # Python 3.x\n",
    "    from urllib.parse import quote as urlencode\n",
    "    from urllib.request import urlretrieve\n",
    "except ImportError:  # Python 2.x\n",
    "    from urllib import pathname2url as urlencode\n",
    "    from urllib import urlretrieve\n",
    "\n",
    "try: # Python 3.x\n",
    "    import http.client as httplib \n",
    "except ImportError:  # Python 2.x\n",
    "    import httplib   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Useful functions\n",
    "\n",
    "Execute PS1 searches and resolve names using [MAST query](https://mast.stsci.edu/api/v0/MastApiTutorial.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ps1cone(ra,dec,radius,table=\"mean\",release=\"dr1\",format=\"csv\",columns=None,\n",
    "           baseurl=\"https://catalogs.mast.stsci.edu/api/v0.1/panstarrs\", verbose=False,\n",
    "           **kw):\n",
    "    \"\"\"Do a cone search of the PS1 catalog\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    ra (float): (degrees) J2000 Right Ascension\n",
    "    dec (float): (degrees) J2000 Declination\n",
    "    radius (float): (degrees) Search radius (<= 0.5 degrees)\n",
    "    table (string): mean, stack, or detection\n",
    "    release (string): dr1 or dr2\n",
    "    format: csv, votable, json\n",
    "    columns: list of column names to include (None means use defaults)\n",
    "    baseurl: base URL for the request\n",
    "    verbose: print info about request\n",
    "    **kw: other parameters (e.g., 'nDetections.min':2)\n",
    "    \"\"\"\n",
    "    \n",
    "    data = kw.copy()\n",
    "    data['ra'] = ra\n",
    "    data['dec'] = dec\n",
    "    data['radius'] = radius\n",
    "    return ps1search(table=table,release=release,format=format,columns=columns,\n",
    "                    baseurl=baseurl, verbose=verbose, **data)\n",
    "\n",
    "\n",
    "def ps1search(table=\"mean\",release=\"dr1\",format=\"csv\",columns=None,\n",
    "           baseurl=\"https://catalogs.mast.stsci.edu/api/v0.1/panstarrs\", verbose=False,\n",
    "           **kw):\n",
    "    \"\"\"Do a general search of the PS1 catalog (possibly without ra/dec/radius)\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    table (string): mean, stack, or detection\n",
    "    release (string): dr1 or dr2\n",
    "    format: csv, votable, json\n",
    "    columns: list of column names to include (None means use defaults)\n",
    "    baseurl: base URL for the request\n",
    "    verbose: print info about request\n",
    "    **kw: other parameters (e.g., 'nDetections.min':2).  Note this is required!\n",
    "    \"\"\"\n",
    "    \n",
    "    data = kw.copy()\n",
    "    if not data:\n",
    "        raise ValueError(\"You must specify some parameters for search\")\n",
    "    checklegal(table,release)\n",
    "    if format not in (\"csv\",\"votable\",\"json\"):\n",
    "        raise ValueError(\"Bad value for format\")\n",
    "    url = f\"{baseurl}/{release}/{table}.{format}\"\n",
    "    if columns:\n",
    "        # check that column values are legal\n",
    "        # create a dictionary to speed this up\n",
    "        dcols = {}\n",
    "        for col in ps1metadata(table,release)['name']:\n",
    "            dcols[col.lower()] = 1\n",
    "        badcols = []\n",
    "        for col in columns:\n",
    "            if col.lower().strip() not in dcols:\n",
    "                badcols.append(col)\n",
    "        if badcols:\n",
    "            raise ValueError('Some columns not found in table: {}'.format(', '.join(badcols)))\n",
    "        # two different ways to specify a list of column values in the API\n",
    "        # data['columns'] = columns\n",
    "        data['columns'] = '[{}]'.format(','.join(columns))\n",
    "\n",
    "# either get or post works\n",
    "#    r = requests.post(url, data=data)\n",
    "    r = requests.get(url, params=data)\n",
    "\n",
    "    if verbose:\n",
    "        print(r.url)\n",
    "    r.raise_for_status()\n",
    "    if format == \"json\":\n",
    "        return r.json()\n",
    "    else:\n",
    "        return r.text\n",
    "\n",
    "\n",
    "def checklegal(table,release):\n",
    "    \"\"\"Checks if this combination of table and release is acceptable\n",
    "    \n",
    "    Raises a VelueError exception if there is problem\n",
    "    \"\"\"\n",
    "    \n",
    "    releaselist = (\"dr1\", \"dr2\")\n",
    "    if release not in (\"dr1\",\"dr2\"):\n",
    "        raise ValueError(\"Bad value for release (must be one of {})\".format(', '.join(releaselist)))\n",
    "    if release==\"dr1\":\n",
    "        tablelist = (\"mean\", \"stack\")\n",
    "    else:\n",
    "        tablelist = (\"mean\", \"stack\", \"detection\")\n",
    "    if table not in tablelist:\n",
    "        raise ValueError(\"Bad value for table (for {} must be one of {})\".format(release, \", \".join(tablelist)))\n",
    "\n",
    "\n",
    "def ps1metadata(table=\"mean\",release=\"dr1\",\n",
    "           baseurl=\"https://catalogs.mast.stsci.edu/api/v0.1/panstarrs\"):\n",
    "    \"\"\"Return metadata for the specified catalog and table\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    table (string): mean, stack, or detection\n",
    "    release (string): dr1 or dr2\n",
    "    baseurl: base URL for the request\n",
    "    \n",
    "    Returns an astropy table with columns name, type, description\n",
    "    \"\"\"\n",
    "    \n",
    "    checklegal(table,release)\n",
    "    url = f\"{baseurl}/{release}/{table}/metadata\"\n",
    "    r = requests.get(url)\n",
    "    r.raise_for_status()\n",
    "    v = r.json()\n",
    "    # convert to astropy table\n",
    "    tab = Table(rows=[(x['name'],x['type'],x['description']) for x in v],\n",
    "               names=('name','type','description'))\n",
    "    return tab\n",
    "\n",
    "\n",
    "def mastQuery(request):\n",
    "    \"\"\"Perform a MAST query.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    request (dictionary): The MAST request json object\n",
    "\n",
    "    Returns head,content where head is the response HTTP headers, and content is the returned data\n",
    "    \"\"\"\n",
    "    \n",
    "    server='mast.stsci.edu'\n",
    "\n",
    "    # Grab Python Version \n",
    "    version = \".\".join(map(str, sys.version_info[:3]))\n",
    "\n",
    "    # Create Http Header Variables\n",
    "    headers = {\"Content-type\": \"application/x-www-form-urlencoded\",\n",
    "               \"Accept\": \"text/plain\",\n",
    "               \"User-agent\":\"python-requests/\"+version}\n",
    "\n",
    "    # Encoding the request as a json string\n",
    "    requestString = json.dumps(request)\n",
    "    requestString = urlencode(requestString)\n",
    "    \n",
    "    # opening the https connection\n",
    "    conn = httplib.HTTPSConnection(server)\n",
    "\n",
    "    # Making the query\n",
    "    conn.request(\"POST\", \"/api/v0/invoke\", \"request=\"+requestString, headers)\n",
    "\n",
    "    # Getting the response\n",
    "    resp = conn.getresponse()\n",
    "    head = resp.getheaders()\n",
    "    content = resp.read().decode('utf-8')\n",
    "\n",
    "    # Close the https connection\n",
    "    conn.close()\n",
    "\n",
    "    return head,content\n",
    "\n",
    "\n",
    "def resolve(name):\n",
    "    \"\"\"Get the RA and Dec for an object using the MAST name resolver\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    name (str): Name of object\n",
    "\n",
    "    Returns RA, Dec tuple with position\"\"\"\n",
    "\n",
    "    resolverRequest = {'service':'Mast.Name.Lookup',\n",
    "                       'params':{'input':name,\n",
    "                                 'format':'json'\n",
    "                                },\n",
    "                      }\n",
    "    headers,resolvedObjectString = mastQuery(resolverRequest)\n",
    "    resolvedObject = json.loads(resolvedObjectString)\n",
    "    # The resolver returns a variety of information about the resolved object, \n",
    "    # however for our purposes all we need are the RA and Dec\n",
    "    try:\n",
    "        objRa = resolvedObject['resolvedCoordinate'][0]['ra']\n",
    "        objDec = resolvedObject['resolvedCoordinate'][0]['decl']\n",
    "    except IndexError as e:\n",
    "        raise ValueError(\"Unknown object '{}'\".format(name))\n",
    "    return (objRa, objDec)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use metadata query to get information on available columns\n",
    "\n",
    "This query works for any of the tables in the API (mean, stack, detection)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "meta = ps1metadata(\"mean\",\"dr2\")\n",
    "meta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple positional query\n",
    "\n",
    "Search mean object table with nDetections > 1.\n",
    "\n",
    "This searches the mean object catalog for objects within 50 arcsec of M87 (RA=187.706, Dec=12.391 in degrees). 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ra = 187.706\n",
    "dec = 12.391\n",
    "radius = 50.0/3600.0\n",
    "constraints = {'nDetections.gt':1}\n",
    "\n",
    "# strip blanks and weed out blank and commented-out values\n",
    "columns = \"\"\"objID,raMean,decMean,nDetections,ng,nr,ni,nz,ny,\n",
    "    gMeanPSFMag,rMeanPSFMag,iMeanPSFMag,zMeanPSFMag,yMeanPSFMag\"\"\".split(',')\n",
    "columns = [x.strip() for x in columns]\n",
    "columns = [x for x in columns if x and not x.startswith('#')]\n",
    "results = ps1cone(ra,dec,radius,release='dr2',columns=columns,verbose=True,**constraints)\n",
    "# print first few lines\n",
    "lines = results.split('\\n')\n",
    "print(len(lines),\"rows in results -- first 5 rows:\")\n",
    "print('\\n'.join(lines[:6]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert the results to an astropy table\n",
    "\n",
    "The CSV results string is easily converted to an [astropy table](http://docs.astropy.org/en/stable/table/).  This table is easily manipulated to extract information on individual columns or rows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab = ascii.read(results)\n",
    "# improve the format\n",
    "for filter in 'grizy':\n",
    "    col = filter+'MeanPSFMag'\n",
    "    try:\n",
    "        tab[col].format = \".4f\"\n",
    "        tab[col][tab[col] == -999.0] = np.nan\n",
    "    except KeyError:\n",
    "        print(\"{} not found\".format(col))\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Query for a single object by objID\n",
    "\n",
    "It is possible to query the catalog directly on the object identifier without any RA/Dec restriction.  This might not be particularly useful when search for objects, but it is very useful when searching the detection table for time-series data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results1 = ps1search(release='dr2',columns=columns,verbose=True,objid=122851876947049923)\n",
    "tab1 = ascii.read(results1)\n",
    "# improve the format\n",
    "for filter in 'grizy':\n",
    "    col = filter+'MeanPSFMag'\n",
    "    try:\n",
    "        tab1[col].format = \".4f\"\n",
    "        tab1[col][tab1[col] == -999.0] = np.nan\n",
    "    except KeyError:\n",
    "        print(\"{} not found\".format(col))\n",
    "tab1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search stack objects at same position\n",
    "\n",
    "There is no need for the `nDetections` limit for stack objects, which can in fact have `nDetections = 0` for objects that are too faint to be detected on single-epoch images.  But we require `primaryDetection=1` in order to eliminate duplicate sources at the edges of the skycell regions used for processing.  (There is another column `bestDetection` that would be better suited for this test but is currently not correct in the database.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sconstraints = {'primaryDetection':1}\n",
    "\n",
    "scolumns = \"\"\"objID,raMean,decMean,nDetections,ng,nr,ni,nz,ny,\n",
    "    nStackDetections,primaryDetection,\n",
    "    gPSFMag,rPSFMag,iPSFMag,zPSFMag,yPSFMag\"\"\".split(',')\n",
    "# strip blanks and weed out blank and commented-out values\n",
    "scolumns = [x.strip() for x in scolumns]\n",
    "scolumns = [x for x in scolumns if x and not x.startswith('#')]\n",
    "\n",
    "sresults = ps1cone(ra,dec,radius,table=\"stack\",release=\"dr2\",columns=scolumns,verbose=True,**sconstraints)\n",
    "stab = ascii.read(sresults)\n",
    "for col in scolumns:\n",
    "    try:\n",
    "        stab[col]\n",
    "    except KeyError:\n",
    "        print(col,\"not found\")\n",
    "# improve the format\n",
    "for filter in 'grizy':\n",
    "    col = filter+'PSFMag'\n",
    "    try:\n",
    "        stab[col].format = \".4f\"\n",
    "        stab[col][stab[col] == -999.0] = np.nan\n",
    "    except KeyError:\n",
    "        print(\"{} not found\".format(col))\n",
    "stab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explore the tables\n",
    "\n",
    "Match the stack and mean tables and look at the subset of sources that are detected in both catalogs.  The commented-out lines show how to restrict the joined table to only stack detections and to only stack non-detections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import join\n",
    "jtab = join(stab,tab,join_type='outer')\n",
    "jtab.sort('objID')\n",
    "jtab_both = jtab[(jtab['primaryDetection']==1) & (jtab['nDetections']>1)]\n",
    "\n",
    "#jtab[jtab['nStackDetections'].mask].show_in_notebook()\n",
    "#jtab[~jtab['nStackDetections'].mask].show_in_notebook()\n",
    "#jtab.show_in_notebook()\n",
    "jtab_both.show_in_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot RA & Dec positions of mean and stack objects\n",
    "\n",
    "Note that `raMean` and `decMean` are defined for all objects, including stack-only objects.   For objects detected only on the stacked images, the `raStack` and `decStack` values are given in the `raMean` and `decMean` columns.  That makes it simple to analyze the positions without testing to see what positions are available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'font.size': 16})\n",
    "plt.figure(1,(10,10))\n",
    "plt.plot(tab['raMean'], tab['decMean'], 'ro', label='Mean (nDet>1)')\n",
    "plt.plot(stab['raMean'], stab['decMean'], 'bo', label='Stack')\n",
    "plt.plot(jtab_both['raMean'], jtab_both['decMean'], 'go', label='Both')\n",
    "\n",
    "plt.xlabel('RA [deg]')\n",
    "plt.ylabel('Dec [deg]')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "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.  A small search radius is used so only a single object is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objname = 'KQ UMa'\n",
    "ra, dec = resolve(objname)\n",
    "radius = 1.0/3600.0 # radius = 1 arcsec\n",
    "\n",
    "results = ps1cone(ra,dec,radius,release='dr2',columns=columns,**constraints)\n",
    "tab = ascii.read(results)\n",
    "# improve the format\n",
    "for filter in 'grizy':\n",
    "    col = filter+'MeanPSFMag'\n",
    "    tab[col].format = \".4f\"\n",
    "    tab[col][tab[col] == -999.0] = np.nan\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the detection information\n",
    "\n",
    "Extract all the objects with the same object ID from the Detection table, which contains all the individual measurements for this source."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def addfilter(dtab):\n",
    "    \"\"\"Add filter name as column in detection table by translating filterID\n",
    "    \n",
    "    This modifies the table in place.  If the 'filter' column already exists,\n",
    "    the table is returned unchanged.\n",
    "    \"\"\"\n",
    "    if 'filter' not in dtab.colnames:\n",
    "        # the filterID value goes from 1 to 5 for grizy\n",
    "        id2filter = np.array(list('grizy'))\n",
    "        dtab['filter'] = id2filter[(dtab['filterID']-1).data]\n",
    "    return dtab\n",
    "\n",
    "objid = tab['objID'][0]\n",
    "dconstraints = {'objID': objid}\n",
    "dcolumns = (\"\"\"objID,detectID,filterID,obsTime,ra,dec,psfFlux,psfFluxErr,psfMajorFWHM,psfMinorFWHM,\n",
    "            psfQfPerfect,apFlux,apFluxErr,infoFlag,infoFlag2,infoFlag3\"\"\").split(',')\n",
    "# strip blanks and weed out blank and commented-out values\n",
    "dcolumns = [x.strip() for x in dcolumns]\n",
    "dcolumns = [x for x in dcolumns if x and not x.startswith('#')]\n",
    "\n",
    "dresults = ps1search(table='detection',release='dr2',columns=dcolumns,**dconstraints)\n",
    "dtab = addfilter(ascii.read(dresults))\n",
    "dtab.sort('obsTime')\n",
    "dtab"
   ]
  },
  {
   "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 = dtab['obsTime']\n",
    "mag = -2.5*np.log10(dtab['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",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.figure(1,(10,10))\n",
    "for i, filter in enumerate(\"grizy\"):\n",
    "    plt.subplot(511+i)\n",
    "    w = np.where(dtab['filter']==filter)\n",
    "    plt.plot(t[w],mag[w],'-o')\n",
    "    plt.ylabel(filter+' [mag]')\n",
    "    plt.xlim(xlim)\n",
    "    plt.gca().invert_yaxis()\n",
    "    if i==0:\n",
    "        plt.title(objname)\n",
    "plt.xlabel('Time [MJD]')\n",
    "plt.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 = dtab['obsTime']\n",
    "mag = -2.5*np.log10(dtab['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",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.figure(1,(10,10))\n",
    "for i, filter in enumerate(\"grizy\"):\n",
    "    plt.subplot(511+i)\n",
    "    w = np.where(dtab['filter']==filter)\n",
    "    magmean = tab[filter+'MeanPSFMag'][0]\n",
    "    plt.plot(t[w],mag[w] - magmean,'-o')\n",
    "    plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))\n",
    "    plt.xlim(xlim)\n",
    "    plt.gca().invert_yaxis()\n",
    "    if i==0:\n",
    "        plt.title(objname)\n",
    "plt.xlabel('Time [MJD]')\n",
    "plt.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": [
    "dtab[ (dtab['filter']=='z') & (np.abs(mag-tab['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(dtab['filter']=='z')\n",
    "zdtab = dtab[w]\n",
    "zdtab['mag'] = mag[w]\n",
    "zdtab['dmag'] = zdtab['mag'] - tab['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 = dtab['obsTime']\n",
    "mag = -2.5*np.log10(dtab['psfFlux']) + 8.90\n",
    "magmean = 0.0*mag\n",
    "for filter in \"grizy\":\n",
    "    magmean[dtab['filter']==filter] = tab[filter+'MeanPSFMag'][0]\n",
    "dmag = mag - magmean\n",
    "dmag1 = dmag[dtab['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",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.figure(1,(10,10))\n",
    "for i, filter in enumerate(\"grizy\"):\n",
    "    plt.subplot(511+i)\n",
    "    w = np.where((dtab['filter']==filter) & (dtab['psfQfPerfect']>0.9))[0]\n",
    "    plt.plot(t[w],dmag[w],'-o')\n",
    "    plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean[w[0]]))\n",
    "    plt.xlim(xlim)\n",
    "    plt.ylim(ylim)\n",
    "    if i==0:\n",
    "        plt.title(objname)\n",
    "plt.xlabel('Time [MJD]')\n",
    "plt.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\n",
    "# convert flux in Jy to magnitudes\n",
    "t = (dtab['obsTime'] % period) / period\n",
    "mag = -2.5*np.log10(dtab['psfFlux']) + 8.90\n",
    "magmean = 0.0*mag\n",
    "for filter in \"grizy\":\n",
    "    magmean[dtab['filter']==filter] = tab[filter+'MeanPSFMag'][0]\n",
    "dmag = mag - magmean\n",
    "dmag1 = dmag[dtab['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",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.figure(1,(10,10))\n",
    "for i, filter in enumerate(\"grizy\"):\n",
    "    plt.subplot(511+i)\n",
    "    w = np.where((dtab['filter']==filter) & (dtab['psfQfPerfect']>0.9))[0]\n",
    "    w = w[np.argsort(t[w])]\n",
    "    plt.plot(t[w],dmag[w],'-o')\n",
    "    plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean[w[0]]))\n",
    "    plt.xlim(xlim)\n",
    "    plt.ylim(ylim)\n",
    "    if i==0:\n",
    "        plt.title(objname)\n",
    "plt.xlabel('Phase')\n",
    "plt.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",
    "ra, dec = resolve(objname)\n",
    "radius = 1.0/3600.0 # radius = 1 arcsec\n",
    "\n",
    "results = ps1cone(ra,dec,radius,release='dr2',columns=columns,**constraints)\n",
    "tab = ascii.read(results)\n",
    "# improve the format\n",
    "for filter in 'grizy':\n",
    "    col = filter+'MeanPSFMag'\n",
    "    tab[col].format = \".4f\"\n",
    "    tab[col][tab[col] == -999.0] = np.nan\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the 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 = tab['objID'][0]\n",
    "dconstraints = {'objID': objid, 'psfQfPerfect.min': 0.9}\n",
    "dcolumns = (\"\"\"objID,detectID,filterID,obsTime,ra,dec,psfFlux,psfFluxErr,psfMajorFWHM,psfMinorFWHM,\n",
    "            psfQfPerfect,apFlux,apFluxErr,infoFlag,infoFlag2,infoFlag3\"\"\").split(',')\n",
    "# strip blanks and weed out blank and commented-out values\n",
    "dcolumns = [x.strip() for x in dcolumns]\n",
    "dcolumns = [x for x in dcolumns if x and not x.startswith('#')]\n",
    "\n",
    "dresults = ps1search(table='detection',release='dr2',columns=dcolumns,**dconstraints)\n",
    "dtab = addfilter(ascii.read(dresults))\n",
    "dtab.sort('obsTime')\n",
    "\n",
    "# add magnitude and difference from mean\n",
    "dtab['magmean'] = 0.0\n",
    "for filter in \"grizy\":\n",
    "    dtab['magmean'][dtab['filter']==filter] = tab[filter+'MeanPSFMag'][0]\n",
    "dtab['mag'] = -2.5*np.log10(dtab['psfFlux']) + 8.90\n",
    "dtab['dmag'] = dtab['mag']-dtab['magmean']\n",
    "dtab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = dtab['obsTime']\n",
    "dmag = dtab['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",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.figure(1,(10,10))\n",
    "for i, filter in enumerate(\"grizy\"):\n",
    "    plt.subplot(511+i)\n",
    "    w = np.where(dtab['filter']==filter)[0]\n",
    "    plt.plot(t[w],dmag[w],'-o')\n",
    "    magmean = dtab['magmean'][w[0]]\n",
    "    plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))\n",
    "    plt.xlim(xlim)\n",
    "    plt.ylim(ylim)\n",
    "    if i==0:\n",
    "        plt.title(objname)\n",
    "plt.xlabel('Time [MJD]')\n",
    "plt.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 = ((dtab['obsTime']-bjd0) % period) / period\n",
    "dmag = dtab['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",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.figure(1,(10,6))\n",
    "plt.plot(t,dmag,'-o')\n",
    "plt.xlim(xlim)\n",
    "plt.ylim(ylim)\n",
    "plt.xlabel('Phase')\n",
    "plt.ylabel('Delta magnitude from mean [mag]')\n",
    "plt.title(objname)\n",
    "plt.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",
    "ra, dec = resolve(objname)\n",
    "radius = 1.0/3600.0 # radius = 1 arcsec\n",
    "\n",
    "results = ps1cone(ra,dec,radius,release='dr2',columns=columns,**constraints)\n",
    "tab = ascii.read(results)\n",
    "# improve the format\n",
    "for filter in 'grizy':\n",
    "    col = filter+'MeanPSFMag'\n",
    "    tab[col].format = \".4f\"\n",
    "    tab[col][tab[col] == -999.0] = np.nan\n",
    "tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objid = tab['objID'][0]\n",
    "dconstraints = {'objID': objid, 'psfQfPerfect.min': 0.9}\n",
    "dcolumns = (\"\"\"objID,detectID,filterID,obsTime,ra,dec,psfFlux,psfFluxErr,psfMajorFWHM,psfMinorFWHM,\n",
    "            psfQfPerfect,apFlux,apFluxErr,infoFlag,infoFlag2,infoFlag3\"\"\").split(',')\n",
    "# strip blanks and weed out blank and commented-out values\n",
    "dcolumns = [x.strip() for x in dcolumns]\n",
    "dcolumns = [x for x in dcolumns if x and not x.startswith('#')]\n",
    "\n",
    "dresults = ps1search(table='detection',release='dr2',columns=dcolumns,**dconstraints)\n",
    "dtab = addfilter(ascii.read(dresults))\n",
    "dtab.sort('obsTime')\n",
    "\n",
    "# add magnitude and difference from mean\n",
    "dtab['magmean'] = 0.0\n",
    "for filter in \"grizy\":\n",
    "    dtab['magmean'][dtab['filter']==filter] = tab[filter+'MeanPSFMag'][0]\n",
    "dtab['mag'] = -2.5*np.log10(dtab['psfFlux']) + 8.90\n",
    "dtab['dmag'] = dtab['mag']-dtab['magmean']\n",
    "dtab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = dtab['obsTime']\n",
    "dmag = dtab['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",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.figure(1,(10,10))\n",
    "for i, filter in enumerate(\"grizy\"):\n",
    "    plt.subplot(511+i)\n",
    "    w = np.where(dtab['filter']==filter)[0]\n",
    "    plt.plot(t[w],dmag[w],'-o')\n",
    "    magmean = dtab['magmean'][w[0]]\n",
    "    plt.ylabel('{} [mag - {:.2f}]'.format(filter,magmean))\n",
    "    plt.xlim(xlim)\n",
    "    plt.ylim(ylim)\n",
    "    if i==0:\n",
    "        plt.title(objname)\n",
    "plt.xlabel('Time [MJD]')\n",
    "plt.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 = ((dtab['obsTime']-bjd0) % period) / period\n",
    "dmag = dtab['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",
    "plt.rcParams.update({'font.size': 14})\n",
    "plt.figure(1,(12,6))\n",
    "plt.plot(t,dmag,'-o')\n",
    "plt.xlim(xlim)\n",
    "plt.ylim(ylim)\n",
    "plt.xlabel('Phase')\n",
    "plt.ylabel('Delta magnitude from mean [mag]')\n",
    "plt.title(objname)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}