{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The warnings filter gets rid of the deprecation and future warnings that sometimes come with imports and running some cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# essential imports\n",
    "import warnings\n",
    "warnings.simplefilter('ignore')\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from astropy import coordinates as coord\n",
    "from astropy import units as u\n",
    "from astroquery.vizier import Vizier\n",
    "from astroquery.gaia import Gaia\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We specify an instance of the VizierClass in order to control the conditions and output columns. In this segment we need to specify that we want all the objects from the catalogs (the dafault value is 50 - the same as in web Vizier) and all the catalog columns. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = Vizier(columns=[\"**\"]) # we want all the catalog columns, akin to checking all fields in the Vizier web interface\n",
    "v.ROW_LIMIT = -1 # unlimited output, akin to selecting unlimited output from the drop-down menu in the web interface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the open cluster data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "astroquery.vizier acts in a similar way to VizieR web interface. It has the same functionality - e.g. querying an objec or querying a region around the target. The choice of catalogs, keywords as well as filters on individual columns can be also specified.\n",
    "\n",
    "We know in advance what catalog we wil be using, so we query using the catalog name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_list = v.find_catalogs('Gaia DR2 open clusters in the Milky Way (Cantat-Gaudin+, 2018)')\n",
    "print({k:v.description for k,v in catalog_list.items()})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The search yielded two catalogs, we only want the latter one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(catalog_list.keys())\n",
    "print(sorted(catalog_list.keys())[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order not to download useless data needlessly, we specify the catalog. This is normally not a concern for most smaller tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaudin = v.get_catalogs(sorted(catalog_list.keys())[-1])\n",
    "print(gaudin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen the catalog contains two tables, with cluster parameters and individual member stars. We will load each of them in a separate table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters = gaudin[0]\n",
    "members = gaudin[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "members"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes, the number of columns is so large, that they will not be displayed in the entirety, just like rows in this example. In order to see all available columns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "members.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basic table info can be obtained as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "members.info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters.info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The astropy tables makes use of numpy structured arrays, so the data can be handled in a classic numpy fashion:\n",
    "    \n",
    "e.g.\n",
    "\n",
    "Get all the data for the first cluster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "get all the distance modes for all clusters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clusters['dmode']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check how the cluster distances are distributed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 10))\n",
    "ax.hist(clusters['dmode'], bins=int(np.sqrt(len(clusters['dmode']))))\n",
    "ax.set_xlabel('distance [pc]')\n",
    "ax.set_ylabel('no. of objects')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also check the distribution in the Galactic X,Y,Z coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 1, figsize=(10, 20))\n",
    "ax[0].scatter(clusters['X'], clusters['Y'], alpha=0.5, marker='+')\n",
    "ax[0].set_xlabel('X [pc]')\n",
    "ax[0].set_ylabel('Y [pc]')\n",
    "\n",
    "ax[1].scatter(clusters['X'], clusters['Y'], alpha=0.5, marker='+')\n",
    "ax[1].set_xlabel('X [pc]')\n",
    "ax[1].set_ylabel('Y [pc]')\n",
    "ax[1].set_xlim(-5000, 5000)\n",
    "ax[1].set_ylim(-5000, 5000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the WD data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We repeat the above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog_list = v.find_catalogs('Gaia DR2 white dwarf candidates (Gentile Fusillo+, 2019)')\n",
    "print({k:v.description for k,v in catalog_list.items()})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fusillo = v.get_catalogs(catalog_list.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fusillo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fusillo = fusillo[0]\n",
    "fusillo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fusillo.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fusillo.info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(3, 1,figsize=(10, 15))\n",
    "ax[0].hist(1000/fusillo['Plx'], bins=int(np.sqrt(len(fusillo))))\n",
    "ax[0].set_xlabel('distance [pc]')\n",
    "ax[0].set_ylabel('no. of objects')\n",
    "ax[0].set_title('WD distances histogram')\n",
    "\n",
    "ax[1].hist(1000/fusillo['Plx'], bins=int(np.sqrt(len(fusillo))), histtype='step', lw=2, \n",
    "           range=(0,4000), label='WDs')\n",
    "ax[1].hist(clusters['dmode'], bins=int(np.sqrt(len(clusters['dmode']))), histtype='step', lw=2, \n",
    "           range=(0,4000), label='OCs')\n",
    "ax[1].set_title('WD vs star cluster distance distribution')\n",
    "ax[1].set_xlabel('distance [pc]')\n",
    "ax[1].set_ylabel('no. of objects')\n",
    "ax[1].legend()\n",
    "\n",
    "ax[2].hist(1000/fusillo['Plx'], bins=int(np.sqrt(len(fusillo))), histtype='step', lw=2,\n",
    "           range=(0,4000), log=True, label='WDs')\n",
    "ax[2].hist(clusters['dmode'], bins=int(np.sqrt(len(clusters['dmode']))), histtype='step', lw=2,\n",
    "           range=(0,4000), log=True, label='OCs')\n",
    "ax[2].set_xlabel('distance [pc]')\n",
    "ax[2].set_ylabel('no. of objects')\n",
    "ax[2].legend()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### To be discussed:\n",
    "\n",
    "We see that we can all needed parameters to compute the matches. We will consider all WDs projected within 1.5$r_{50}$ of the cluster. The potential matches in the parallax space will be WDs that have parallaxes in the interval of $(\\varpi - 3\\sigma; \\varpi + 3\\sigma)$ of the cluster parallax, where $\\varpi$ is the parallax (plx) of the cluster and $\\sigma$ is the standard deviation of parallax (s_plx) of the cluster members. Similarly, corresponding limits will be also applied in the proper motion space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with coordinates\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Astropy coordinates package provides classes for representing a variety of celestial/spatial coordinates and their velocity components, as well as tools for converting between common coordinate systems in a uniform way.\n",
    "\n",
    "The best way to start using coordinates is to use the SkyCoord class. SkyCoord objects are instantiated by passing in positions (and optional velocities) with specified units and a coordinate frame. Sky positions are commonly passed in as Quantity objects and the frame is specified with the string name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = coord.SkyCoord(ra=15.0*u.degree, dec=-10*u.degree, frame='icrs')\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = coord.SkyCoord(15, -10, frame='icrs', unit='deg')\n",
    "c = coord.SkyCoord('1h00m00s', '-10d00m00s', frame='icrs')\n",
    "c = coord.SkyCoord('1h00m00s', '-10d00m00s')\n",
    "c = coord.SkyCoord('1 00 00 -10 00 00', unit=(u.hourangle, u.deg))\n",
    "c = coord.SkyCoord('1:00 -10:00', unit=(u.hourangle, u.deg))\n",
    "c # All of these are equivalent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Coordinates can be manipulated in various ways:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.ra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.ra.deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.ra.hour"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.ra.hms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.ra.radian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.to_string('decimal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.to_string('hmsdms')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.galactic # transformation to galactic coodinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding distance information opens up new possibilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c1 = coord.SkyCoord(ra=15.0*u.degree, dec=-10*u.degree, distance=150*u.pc, frame='icrs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c1.cartesian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c1.cartesian.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distance between the coordinates\n",
    "\n",
    "Various options exist in order to compute the projected distance between the two set of celestial coordinates. One option is to implement the standard equation:\n",
    "$$ \\cos \\theta = \\sin \\delta_{1} \\sin \\delta_{2} + \\cos \\delta_{1} \\cos \\delta_{2} \\cos(\\alpha_{1} - \\alpha_{2}),$$\n",
    "or to use the separation method from Astropy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = coord.SkyCoord(ra=15.0*u.degree, dec=-12.76*u.degree, frame='icrs')\n",
    "separ = c.separation(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "separ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "separ.arcsec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "separ.deg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coordinates can be easily loaded from any table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_cluster = coord.SkyCoord(clusters['RAJ2000'], clusters['DEJ2000'])\n",
    "c_cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(c_cluster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b.separation(c_cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NGC 663 example case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Be/X-ray binary RXJ0146.9+6121 belongs to the class of low-luminosity persistent systems with Be-companions (see Reig, 2011, for details). It may have been formed by an electron-capture supernova and is a runaway from open cluster NGC 663. It would be very interesting to see whether there any any massive white dwarfs in this cluster, as it would provide a direct empirical evidence for multiple supernova channels.\n",
    "\n",
    "From positional constraints, we determined that two WDs, with GDR2 identifiers 511192470770710656 and 511218171855356544 are projected in the cluster vicinity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wd1 = fusillo[fusillo['Source'] == 511192470770710656]\n",
    "wd2 = fusillo[fusillo['Source'] == 511218171855356544]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wd1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wd2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc663_stars = members[members['Cluster'] == 'NGC_663']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc663_stars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 1, figsize=(10,20))\n",
    "ax[0].scatter(ngc663_stars['pmRA'], ngc663_stars['pmDE'], marker='+', alpha=0.5)\n",
    "ax[0].set_xlabel('proper motion in RA [mas/yr]')\n",
    "ax[0].set_ylabel('proper motion in dec [mas/yr]')\n",
    "\n",
    "ax[1].scatter(ngc663_stars['pmRA'], ngc663_stars['pmDE'], marker='+', alpha=0.5)\n",
    "ax[1].scatter(wd1['pmRA'], wd1['pmDE'], marker='s', color='red')\n",
    "ax[1].scatter(wd2['pmRA'], wd2['pmDE'], marker='s', color='green')\n",
    "ax[1].set_xlabel('proper motion in RA [mas/yr]')\n",
    "ax[1].set_ylabel('proper motion in dec [mas/yr]')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ngc663_members = ngc663_stars[ngc663_stars['PMemb'] >= 0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12, 12))\n",
    "ax.scatter(ngc663_members['BP-RP'], ngc663_members['Gmag'])\n",
    "ax.invert_yaxis()\n",
    "ax.set_ylabel('G [mag]')\n",
    "ax.set_xlabel('BP-RP color [mag]')\n",
    "ax.set_title('color-magnitude diagram of NGC 663')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, the observations from Gaia are not deep enough to catch WD members of a cluster ~ 2 kpc away. Any WD studies of NGC 663 are out of reach, at least from Gaia..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some useful documentation:\n",
    "http://docs.astropy.org/en/stable/coordinates/\n",
    "\n",
    "https://astroquery.readthedocs.io/en/latest/vizier/vizier.html\n"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}