crossmatch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Catalog cross-matching\n",
    "\n",
    "Some conventions I use:\n",
    "- The first letter of `astropy.table` object variable names are capitalized\n",
    "- First letters of `pandas.DataFrame` names are lowercase\n",
    "- CXO, CSC = Chandra X-ray Observatory, Chandra Source Catalog\n",
    "- XMM, SSC = XMM-Newton, Serendipitous Source Catalog\n",
    "- XRT, XPS = (Swift) X-ray Telescope, XRT Point Source (catalog)\n",
    "- BeSS = Be Star Spectra database"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy\n",
    "from astropy import units\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=units.UnitsWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np, pandas as pd, matplotlib.pyplot as plt, astropy.units as u, pickle, seaborn as sns, sympy\n",
    "from gdpyc.core import DustMap, GasMap\n",
    "from mpl_toolkits import mplot3d\n",
    "from scipy import stats, optimize\n",
    "from PyAstronomy import pyasl\n",
    "from astropy.io import fits, ascii as ascii_io, votable\n",
    "from astropy.coordinates import SkyCoord, Angle\n",
    "from astroquery.vizier import Vizier\n",
    "from astroquery.simbad import Simbad\n",
    "from importlib.machinery import SourceFileLoader\n",
    "AsymmetricUncertainty = SourceFileLoader(\"AsymmetricUncertainty\", \"../dark-GRBs/code/uncertainty.py\").load_module().AsymmetricUncertainty\n",
    "from tools import hms_to_dec, dms_to_dec, round_scientific\n",
    "\n",
    "keV = 2.418e17 # Hz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['BeSS #', 'Name', 'RA', 'Dec', 'Vmag', 'SpecType',\n",
       "       'Simbad classification', 'Catalog memberships',\n",
       "       'Literature classification', 'Reference', 'epsi', 'sepsi', 'BP-RP',\n",
       "       'Parallax', 'Parallax_err', 'Gflux', 'Gflux_corrected', 'Gflux_err',\n",
       "       'Gmag', 'Gmag_corrected', 'Gmag_err', 'PM', 'PM_dec', 'pmdec_error',\n",
       "       'PM_RA', 'PM_RA_err', 'rgeo', 'rpgeo', 'Fluxb', 'Fluxh', 'Fluxm',\n",
       "       'Fluxs', 'Fluxu', 'Fluxw', 'HRhm', 'HRhs', 'HRms', 'NHGal', 'S_N',\n",
       "       '2CXO', '2SXPS', 'CR0', 'E_CR0', 'e_CR0', 'SwiftHR1', 'SwiftHR2', 'CR1',\n",
       "       'e_CR1', 'CR2', 'E_CR2', 'e_CR2', 'CR3', 'E_CR3', 'e_CR3', 'PcstS0',\n",
       "       'FPO0', 'FPU0', 'Gam1', 'E_Gam1', 'e_Gam1', '4XMM', 'Flux1', 'e_Flux1',\n",
       "       'Flux2', 'e_Flux2', 'Flux3', 'e_Flux3', 'Flux4', 'e_Flux4', 'Flux5',\n",
       "       'e_Flux5', 'Flux8', 'e_Flux8', 'Flux9', 'e_Flux9', 'HR1', 'e_HR1',\n",
       "       'HR2', 'e_HR2', 'HR3', 'e_HR3', 'HR4', 'e_HR4', '2MASS', 'Jmag', 'Hmag',\n",
       "       'Kmag', 'e_Jmag', 'e_Hmag', 'e_Kmag', 'WISE', 'W1mag', 'e_W1mag',\n",
       "       'W2mag', 'e_W2mag', 'W3mag', 'e_W3mag'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Be_xray_cat = pd.read_csv(\"BeStarXrayCatalog.csv\")\n",
    "Be_xray_cat.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize ViZier query instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "v = Vizier(columns=['all', '_RAJ2000', '_DEJ2000', \"_r\"])\n",
    "v.ROW_LIMIT = -1 # no row limit\n",
    "v.TIMEOUT = 99999 # big time limit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the [Be star database](http://basebe.obspm.fr/basebe/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " #    Be star    Category ...       _RAJ2000            _DEJ2000     \n",
      "                          ...         deg                 deg        \n",
      "---- ---------- --------- ... ------------------- -------------------\n",
      "   1 BD+62 2346 Classical ... 0.35291666666666666   63.50436944444444\n",
      "   2  HD 224905 Classical ...  0.4109583333333333   60.44992222222222\n",
      "   3  HD 225095 Classical ...            0.863125  55.550897222222225\n",
      "   4      2 Cet Classical ...  0.9349583333333333 -17.335991666666665\n",
      "   5     10 Cas Classical ...  1.6105833333333333   64.19616944444445\n",
      "   6 BD+59 2829 Classical ...             1.70125  60.600230555555555\n",
      "   7    BD+62 1 Classical ...  1.8881250000000003  63.080327777777775\n",
      "   8   BD+62 11 Classical ...              2.6955   63.17293611111111\n",
      "   9   V742 Cas Classical ...  2.9047916666666667   58.21182777777778\n",
      "  10   EM* AS 2 Classical ...            3.246875   66.32206388888889\n",
      " ...        ...       ... ...                 ...                 ...\n",
      "2255   V817 Cas Classical ...   356.4794166666666  54.735397222222225\n",
      "2256  HD 223044 Classical ...           356.50175   61.99101944444445\n",
      "2257  HD 223387 Classical ...  357.22766666666666   57.20959722222222\n",
      "2258   V818 Cas Classical ...  357.47133333333335   62.21409166666667\n",
      "2259 EM* AS 518 Classical ...          359.383625   72.88435833333334\n",
      "2260 BD+65 1970 Classical ...   359.3920833333333   66.43162222222222\n",
      "2261     LQ And Classical ...  359.69350000000003  46.413177777777776\n",
      "2262  HD 224544 Classical ...   359.7051666666667  32.381705555555556\n",
      "2263  HD 224599 Classical ...  359.81237500000003   60.02248055555555\n",
      "2264    eps Tuc Classical ...  359.97908333333334  -65.57713333333334\n",
      "Length = 2264 rows\n"
     ]
    }
   ],
   "source": [
    "BeSS = ascii_io.read(\"./catalogs/BeSS_catalog.csv\", format=\"csv\", header_start=0, data_start=1)\n",
    "BeSS[\"_RAJ2000\"] = [hms_to_dec(hms) for hms in BeSS[\"RA\"]]*u.degree\n",
    "BeSS[\"_DEJ2000\"] = [dms_to_dec(dms) for dms in BeSS[\"DEC\"]]*u.degree\n",
    "BeSS.pprint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Crossmatch with Gaia data\n",
    "\n",
    "Create `gaia_df`, a `pandas.DataFrame` containing the single closest match (within $1^{\\prime\\prime}$) in the Gaia catalog for every entry in the BeSS database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "Gaia_BeSS_1as = v.query_region(BeSS, catalog=\"I/350/gaiaedr3\", radius=\"0.5s\")\n",
    "Gaia_matched = Gaia_BeSS_1as[0][\"_q\",\"_r\",\"Source\",\"Plx\",\"pmRA\",\"pmDE\",\"RUWE\",\"epsi\",\"sepsi\",\"Gmag\",\"BPmag\",\"RPmag\",\"Tefftemp\"]\n",
    "\n",
    "gaia_df = Gaia_matched.to_pandas()\n",
    "for i in BeSS[\"#\"]:\n",
    "    instances = len(gaia_df[gaia_df[\"_q\"]==i])\n",
    "    if instances == 0:\n",
    "        pass\n",
    "    elif instances >= 2:\n",
    "        dupeframe = gaia_df[gaia_df[\"_q\"]==i][[\"_q\",\"_r\"]]\n",
    "        gaia_df.drop(dupeframe[dupeframe[\"_r\"]>min(dupeframe[\"_r\"])].index, inplace=True)\n",
    "\n",
    "gaia_df.dropna(axis=0, subset=[\"_r\"], inplace=True)\n",
    "F0 = 2.9979246e-5 * 3228.75 / (6217.59**2) # Jy -> erg/cm^2/s/Ang\n",
    "gaia_df[\"Gflux\"] = F0 * 10**(gaia_df[\"Gmag\"]/-2.5)\n",
    "\n",
    "assert(len(gaia_df) == len(np.unique(gaia_df[\"_q\"]))),\"Something went wrong...\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_dists = v.query_region(BeSS, catalog=\"I/352/gedr3dis\",radius=\"0.5s\")[0].to_pandas()\n",
    "gaia = pd.merge(gaia_df, gaia_dists, on=\"Source\", how=\"left\")\n",
    "gaia.drop(['_r_x', 'Source', '_q_y', '_RAJ2000', '_DEJ2000', '_r_y', 'RA_ICRS', 'DE_ICRS', 'Flag'], axis=1, inplace=True)\n",
    "gaia.rename(columns={\"_q_x\":\"_q\"}, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "#### Notes\n",
    "\n",
    "- 2CXO J210335.6+454505 with 92%/100% is a HMXB in our TD\n",
    "- 2CXO J213930.6+565910 with 83%/99% is a HMXB from Simbad, but not in our TD\n",
    "- 2CXO J112057.1-615500 with 96%/100% is a HMXB from Liu+06 in our TD\n",
    "- 2CXO J223920.8+611626 with 81%/100% is a HMXB from Liu+06 in our TD\n",
    "- 2CXO J174445.7-271344 #  not significant with 45%/70% is a HMXB in our TD, but a Be Star in Simbad\n",
    "- 2CXO J024031.6+611345 #  not significant with 58%/95% is a HMXB from Liu+06\n",
    "- 2CXO J113106.9-625648 #  not significant with 36%/35% is not in our TD, is a HMXB from Simbad\n",
    "- 2CXO J130247.6-635008 #  not significant with 53%/56% is not in our TD, is a HMXB from Simbad, \n",
    "- 2CXO J134632.5-625523 #  not significant with 44%/47% is not in our TD, is a Be Star from Simbad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HMXB crossmatch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "HMXBs_BeSS_1as = v.query_region(BeSS, catalog=\"J/A+A/455/1165/table1\", radius=\"1s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WISE and 2MASS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Twomass_xmatch = v.query_region(BeSS, catalog=\"II/246/out\", radius=\"0.5s\")\n",
    "Catwise_xmatch = v.query_region(BeSS, catalog=\"II/365/catwise\", radius=\"0.5s\")\n",
    "Allwise_xmatch = v.query_region(BeSS, catalog=\"II/328/allwise\", radius=\"0.5s\")\n",
    "Unwise_xmatch = v.query_region(BeSS, catalog=\"II/363/unwise\", radius=\"0.5s\")\n",
    "\n",
    "twomass = Twomass_xmatch[0].to_pandas()\n",
    "catwise = Catwise_xmatch[0].to_pandas()\n",
    "allwise = Allwise_xmatch[0].to_pandas()\n",
    "unwise = Unwise_xmatch[0].to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEaCAYAAAAYOoCaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAcL0lEQVR4nO3deZhkVZnn8e+PTbaSothEtpRFFBDRTlGH7hYBW3ZwYGgYxKLFrtGGcXzAbquFXmxxutAZW54BGqrBtkRGKFAWQRhpoBpxBKmCYncBhKEKiiqWYhFkfeePc0JuhZGZkZkRNzLi/D7Pk0/GXeLe90Rk3jfue26cq4jAzMzKs1qvAzAzs95wAjAzK5QTgJlZoZwAzMwK5QRgZlYoJwAzs0I5AZg1kRSSts+Pz5b0N72OaSIk7SlpySjL+7Zt1hlOAH1O0psknSfpYUnPSVosab/K8j3zAe3Spue9O89f0DRfkh6UdG+Lfe0s6UeSnpK0UtIiSftXln9R0q8lPS9piaSLmp6/lqQnJK0vaYGk3+Z1n5F0o6R3jaPdx0n6eW7z45J+KGlaXralpO/lfT0j6W5Jx7a77aqI+HREfLmNeB6StM9E9tErg9w2a48TQP9bA3gE+BCwAXAKMF/SUGWdFcAHJW1UmTcT+GWL7f0xsCmwraT3NS37AXAt8Ja8zmeBZwEkzQSOAfaJiPWBYeC6FtteHBHP5+kT8rozgAXA+e00WNKHgP8OHBUR04B3AtVkcz7pNdkG2CjH9Xg727Z6SVqj1zEULSL8M2A/wJ3AYfnxnsAS4Gzg+DxvdWAp8LfAgqbnfhO4APg+cEZl/sZAANNH2OcZwDfGiOvrwIn58QLgU5VlOwEvV6ZXA2YDDwBPAvOBGXnZ54HLRtnP88Bu43i9/hJ4DHgU+GRu5/Z52beAUyuvwZXASuAp4Mc5zvOB14EX877/Kq9/MbAMeAa4Edi5ss9vAWcCVwHPAbcA21WW70xKtk+RktcXx3pdWrSr8d6fBCzPbfyzphgm2raDgXvy+guAd1a2+17g9tyui0nJ+dSmmL6QX5vzgQ3zvlcAT+fHW1a2twA4Ffi/OYYfkBL7BaQPILcCQ73+v+vHH58BDBhJmwFvJ/1zVn0b+ER+/FHgbtIBr/rcdYHDSf9YFwBHSlorL34SuB/4jqRD836qbgY+IekvJQ1LWr1FePuTDnjNMa8FHJ230fBfgUNJZzZvJR0YzszLbgE+KulLkvaQ9KYWsZwp6UhJW7eIo7rvfUkJ5SPADsBopY6TSAevTYDNgC8CERHHAP8POCgi1o+Ir+b1r87b3BS4jfSaVh0JfIl0ALwf+EqOaRrwb8A1ue3b88bZ1GivSytvIZ0ZbgEcR3pdNpxM2yS9Hfgu8Lm8/g+BH+QS31rApaTkMiOv97EWMc0gnaHNIiWaf83TW5OSzRktXqtjcju2A36anzMDuA/4u1FeAxtJrzOQfzr3A6xJOnCcU5m3J7AkP/4VsCNwIemA+ykqZwDAx0mfwtYA1iZ9cv1YZfmWpH/MB0ifCm8EdqgsPzrv/zekhPGFyrLtgPsr0wuAF0ifIF/K+9q7svy+punNgVeANfL0fqRPgitJnwq/Dqyel20IzCElwdeAxcD7RnjNvgnMqUy/nZHPAP4BuLyxrGk7D5HKXyO9N9PzdjeobPfcyvL9gZ/nx0cBt4+wnVFfl6Z19yQdTNeozFsOfGAybQP+BphfmV6NdEa5J6nMtxRQZflNrHoG8DKw9iiv1W7A001/KydXpv8ncHVl+iBSabHn/4P99uMzgAEhqXG6/jJwwgirnZ+XfZj0Ka3ZTNI/9qsR8Vvge3keABGxJCJOiIjtSJ/WfkM6s2gsvyAi9iEd7D4NfFnSR/Pi/UmfiKs+GxHTgXWAA4FLJO2al20DXJo7m1eSDnyvkT6dEhFXR8RBpE+AhwDHkhIaEfF0RMyOiJ3z+ouByySpRZvfSuovaHi4xToNXyN9Uv9R7iifPdKKklaXNEfSA5KeJR1EIZVaGpZVHr8ArJ8fb0VKsq2M+rq08GREvDrCfqrabhvpNfvd6xQRr5Newy3ysqWRj8zZI6s+nRX57wtIZ56SzskXMjxL+mAxveksstqH82KL6VZtsjE4AQyAfGA7j3QQOCwiXhlh1fOBvwB+GBEvNG1jS2Av4OOSlklaRioH7S9p4+YNRcQjpNLDLi2WvRIRF5P6IhrL9yeVCn5PRLweET8mHYD+JM9+BNgvIqZXftaOiKUtnnsdcP0IsTwB/A/SgWlGi90/RjrgNoxYMoqI5yLipIjYllQDP1HS3o3FTav/Z1Ji2odUghnK81sloWaPANuOsmzM12W8xtm2R0mJCPjd399WpE/+jwFbNCXbrVZ9+u9t7yTSmen7I+LNpLMIaO+1sklwAhgM/0y6EuagiHhxpJUi4tek2vHJLRYfQ7oqaEfSKfhupHLIEuAoSRvmmvv2klbLSeGT5Lq9pGMlHSBpWl6+H6kj85bct7A7cMNIsUn6IKkjuNF3cTbwFUnb5OWbSDokPz4k1/c3zJet7p7b1YjlNEm7SFoj19M/Qyo/Pdli1/OBYyXtlOMcsZYs6cDcfpFKVq+RSmGQPpFWD9rTSKWtJ4F1SVcttetKYHNJn1O6zHeapPeP9bpMxjjbNh84QNLektYkHcBfInXS/jQ/94T8+h9Ceu9HM430KX6lpBm4nl8bJ4A+lw8E/4V0wF6mdF3985KObrV+RNwUEY+2WDQTOCsillV/SAecmaTS0hCpxv8sqRP5JVLphTzvi6QOw5XAV4HPRMRNpDOLn1ZP+7MzGvGSzk5OiYhGmeh04ApSSeI50sG9cRB8GvhzUp/Gs8B3gK9FRKOTdV1SiWsl8CDp0+rBI7weVwPfIJ1B3J9/j2SH3P7nSQe6syKikdT+ETgll2Y+TyqNPUz6VHwvq3ZwjyoiniN1Sh9EKhP9ilS2g9Ffl8lou20R8QtSf9H/Ap7IcR4UES9HxMvAfyR1OK/M611J+lsZyTdIZcAncnuu6UB7rA1atVRn1nmSzgLujoizeh2L1U/SLcDZEfGvvY7FVuUzAKvDYlp3OtsAkvQhSW/JJaCZwK74U/2U5G/hWddFxNxex2C12pHUT7AeqQR3eEQ81tuQrBWXgMzMCuUSkJlZoZwAzMwK1Vd9ABtvvHEMDQ31Ogwzs76yaNGiJyJik+b5fZUAhoaGWLhwYa/DMDPrK5JaDnHiEpCZWaGcAMzMCuUEYGZWKCcAM7NCOQGYmRXKCcDMrFBOAGZmhXICMDMrVF99EWyq22PO9SxdmW7ItcX0dfjJ7L06tp1ObdvMrMEJoIOWrnyRh+YcAMDQ7Ks6up1ObdvMrMElIDOzQjkBmJkVygnAzKxQTgBmZoVyAjAzK5QTgJlZoZwAzMwK5QRgZlYoJwAzs0I5AZiZFcoJwMysUE4AZmaFcgIwMyuUE4CZWaGcAMzMCuUEYGZWKCcAM7NC1XJHMEkPAc8BrwGvRsSwpBnARcAQ8BBwREQ8XUc8ZmZW7xnAhyNit4gYztOzgesiYgfgujxtZmY16WUJ6BBgXn48Dzi0d6GYmZWnrgQQwI8kLZI0K8/bLCIey4+XAZu1eqKkWZIWSlq4YsWKOmI1MytCLX0AwB9GxFJJmwLXSvp5dWFEhKRo9cSImAvMBRgeHm65jpmZjV8tZwARsTT/Xg5cCuwOPC5pc4D8e3kdsZiZWdL1BCBpPUnTGo+BPwHuBq4AZubVZgKXdzsWMzN7Qx0loM2ASyU19ve/I+IaSbcC8yUdBzwMHFFDLGZmlnU9AUTEg8C7W8x/Eti72/s3M7PW/E1gM7NC1XUVkAF7zLmepStfBGCL6evwk9l79eU+zGwwOAHUaOnKF3lozgEADM2+qm/3YWaDwSUgM7NCOQGYmRXKJaAua67Jm5lNFU4AXVatyZuZTSUuAZmZFcoJwMysUE4AZmaFcgIwMyuUE4CZWaGcAMzMCuXLQAvg8YHMrBUngAJ4fCAza8UlIDOzQvkMYEA0yjzV4Sa2mL4OQ7Ov8hAUZtaSE8CAaDXkhGv9ZjYal4DMzArlBGBmViiXgPpY3UNN+3JSs8HiBNDH6h5q2peTmg0Wl4DMzArlBGBmViiXgKaQVtfyTwVTNS4zmxwngClkqt4+cqrGZWaT4xKQmVmhfAYwTtVySCcug6z7Us7G8BCNx76U06xctSUASasDC4GlEXGgpLcBFwIbAYuAYyLi5brimahGOaRTl0HWXV6pHvB9KadZ2eosAf034L7K9GnAP0XE9sDTwHE1xmJmVrxaEoCkLYEDgHPztIC9gEvyKvOAQ+uIxczMkrrOAL4B/BXwep7eCFgZEa/m6SXAFq2eKGmWpIWSFq5YsaLrgfYDD/NsZp3Q9T4ASQcCyyNikaQ9x/v8iJgLzAUYHh6OzkbXn9xxa2adUEcn8B7AwZL2B9YG3gycDkyXtEY+C9gSWFpDLGZmlnW9BBQRfx0RW0bEEHAkcH1EHA3cAByeV5sJXN7tWMzM7A29/CLYF4ATJd1P6hM4r4exmJkVp9YvgkXEAmBBfvwgsHud+zczszd4KAgzs0J5KAjrON85zKw/OAFYx/nOYWb9wSUgM7NCOQGYmRXKJaAO6Nc7ZlWHlHCd3qw8TgAd0K93zGoc9F2nNyuTS0BmZoVyAjAzK5RLQBPUfGtFM7N+4wQwQe40NbN+5xKQmVmhnADMzArlBNAlnbpto2//aGbd4j6ALulUH4H7GsysW3wGYGZWKJ8B2Cqqw1qMdvbhYSTM+p/PAGwVjWEtGuP5j+Qns/dqaz0zm7qcAMzMCuUEYGZWKPcBtKH5Foed4Ms7zazXnADa0I3hnt1xama95hKQmVmhnADMzArlBGBmVignADOzQjkBmJkVygnAzKxQXU8AktaW9DNJd0i6R9KX8vy3SbpF0v2SLpK0VrdjMTOzN9RxBvASsFdEvBvYDdhX0geA04B/iojtgaeB42qIxczMsq4ngEiez5Nr5p8A9gIuyfPnAYd2OxYzM3vDmAlA0jaT3Ymk1SUtBpYD1wIPACsj4tW8yhJgixGeO0vSQkkLV6xYMdlQzMwsa2coiKslbQr8HLgTuKvxOyKebWcnEfEasJuk6cClwDvaDTAi5gJzAYaHh6Pd55mZ2ejGTAARsZOkNwE7Ae8CdgUOAXaV9FJEvK3dnUXESkk3AB8EpktaI58FbAksnVALzMxsQtrqA4iIlyLidtKn91uAZcBvgTvGeq6kTfInfyStA3wEuA+4ATg8rzYTuHy8wZuZ2cSNeQYgaUfgAOBAYBNSDf8CYFZEvNzGPjYH5klanZRw5kfElZLuBS6UdCpwO3DeBNtgk9QYmrrxeKKqt5M0s6mvnT6A+0gH6NOAyyPipfHsICLuBN7TYv6DwO7j2ZZ1R6eGpu7GsNlm1j3tJIDPALsAxwNnSHqS1BF8F6kj+LLuhWdmZt3STifwOdVpSVvyRmfwYcBlXYnMzMy6atx3BIuIJaTr9q/ufDhmZlYXDwZnZlYoJwAzs0L5pvDWUqcuDTWzqcsJwFrq1KWhZjZ1uQRkZlYoJwAzs0K5BDQKD21gZoPMCWAUHtrAzAaZS0BmZoVyAjAzK5QTgJlZoZwAzMwK5QRgZlYoJwDrqsaQEnvMub7XoZhZE18Gal3VGFKiMa6QmU0dPgMwMyuUE4CZWaGcAMzMCuUEYGZWKCcAM7NCOQGYmRXKl4E2aQwBDb4VYjvqvnVk8/vjO5eZTZwTQBMPAT0+dR+Aq++Pv1tgNjkuAZmZFcpnAFaL5lJR85mDSztm9et6ApC0FfBtYDMggLkRcbqkGcBFwBDwEHBERDzd7XisN6oH9FalG5d2zOpXRwnoVeCkiNgJ+ABwvKSdgNnAdRGxA3BdnjYzs5p0PQFExGMRcVt+/BxwH7AFcAgwL682Dzi027GYmdkbau0DkDQEvAe4BdgsIh7Li5aRSkStnjMLmAWw9dZb1xCl1WWsS27dL2DWXbUlAEnrA98DPhcRz0r63bKICEnR6nkRMReYCzA8PNxyHetPY11y634Bs+6q5TJQSWuSDv4XRMT38+zHJW2el28OLK8jFjMzS7qeAJQ+6p8H3BcRX68sugKYmR/PBC7vdixmZvaGOkpAewDHAHdJWpznfRGYA8yXdBzwMHBEDbGYtaXR/+C+BxtkXU8AEXEToBEW793t/ZtNRKP/wX0PNsg8FISZWaGcAGzKaQwb4dFYzbrLYwHZlOOau1k9fAZgZlYoJwAzs0I5AVhfcL+AWee5D8D6gvsFzDrPZwBmZoVyAjAzK5RLQFa7Our5HkrabGxOAFa7Og7GHkrabGwuAZmZFcoJwMysUEWXgFwn7m+NvoTmeX4fzdpTdAJwnbi/tTrQ+300a59LQGZmhSr6DMDK4FKRWWtOADbwXCoya80lIDOzQjkBmJkVygnABoqHjTZrn/sAbKC4Y9esfT4DMDMrlBOAmVmhiiwBNYaAqNaJXTu2TvNQIzbVFZkAqkNANPif0zrNQ43YVOcSkJlZoZwAzMwK1fUEIOmbkpZLursyb4akayX9Kv/esNtxmJnZquo4A/gWsG/TvNnAdRGxA3BdnjYzsxp1PQFExI3AU02zDwHm5cfzgEO7HYeZma2qV1cBbRYRj+XHy4DNRlpR0ixgFsDWW29dQ2hWkuqlmlW+HNhK0PPLQCMiJMUoy+cCcwGGh4dHXM9sIlpdEmxWil5dBfS4pM0B8u/lPYrDzKxYvUoAVwAz8+OZwOU9isPMrFhdLwFJ+i6wJ7CxpCXA3wFzgPmSjgMeBo7odhxmE9Hu7SS7MexDdciSOvZn5el6AoiIo0ZYtHe39202We3eTrIbwz40tlnX/qw8/iawmVmhen4VkJm1r9VItmYT5QRg1kd82ap1kktAZmaFcgIwMyuUS0BWpJLvAOdLSK3BCcCKVPJBz5eQWoNLQGZmhXICMDMrlEtAZuM00vAQ7Wg1/PRYdfjq/jrdZ1HtC5nMcBOjDVvRrpGG5q7GWnLprhucAMzGaTIHoVbX8Y9Vh+/mQa+x7ckONzHasBXtGus7Du6v6DyXgMzMCuUzADPrilZloU6UilrtA6ZuiajTbe4knwGYWVc0SjrVun6reZ3YRye32WmdbnMnOQGYmRXKCcDMrFDF9AE01wrNumUyw0xM5u+0VT287hp5q0tW25k33rg8LHZnFJMAPIyu1aXTl4lO5LmNg2vdwz60avtY8yYSl/+fO8MlIDOzQjkBmJkVqpgSkFkvdXP46Xa3Pdp6Yw1vMdmafSmm8jX/rTgBmNWgjuEcJrPeWNuYbM2+FJ0YEqNOLgGZmRXKCcDMrFAuAZn1WDeHe+6GyQyH3an9jtaPMdbYQ2MNO13dXre/S9HrsYycAMx6rB86C6t6FW+7/RTV5NSqJt/udwjq+C5Fr2/P6RKQmVmhfAZgNoC6edlpp9UR60TKbGPFNdLd3drZX6uhLFqV1qrLunHm1dMEIGlf4HRgdeDciJjTy3jMBkU/lZXqiHUi+xjrOaOVkiby3NGe063yUM9KQJJWB84E9gN2Ao6StFOv4jEzK00v+wB2B+6PiAcj4mXgQuCQHsZjZlYURURvdiwdDuwbEZ/K08cA74+IE5rWmwXMypM7Ar+YwO42Bp6YRLj9yu0ui9tdlvG0e5uI2KR55pTvBI6IucDcyWxD0sKIGO5QSH3D7S6L212WTrS7lyWgpcBWlekt8zwzM6tBLxPArcAOkt4maS3gSOCKHsZjZlaUnpWAIuJVSScA/4d0Geg3I+KeLu1uUiWkPuZ2l8XtLsuk292zTmAzM+stDwVhZlYoJwAzs0INVAKQtK+kX0i6X9LsFsvfJOmivPwWSUM9CLPj2mj3H0u6TdKr+fsXA6GNdp8o6V5Jd0q6TtI2vYiz09po96cl3SVpsaSbBuUb9mO1u7LeYZJC0kBcGtrG+32spBX5/V4s6VNtbzwiBuKH1JH8ALAtsBZwB7BT0zp/AZydHx8JXNTruGtq9xCwK/Bt4PBex1xjuz8MrJsff6ag9/vNlccHA9f0Ou462p3XmwbcCNwMDPc67pre72OBMyay/UE6A2hnaIlDgHn58SXA3pJUY4zdMGa7I+KhiLgTeL0XAXZJO+2+ISJeyJM3k75r0u/aafezlcn1gEG40qPdoWO+DJwG/LbO4Lqoq0PmDFIC2AJ4pDK9JM9ruU5EvAo8A2xUS3Td0067B9F4230ccHVXI6pHW+2WdLykB4CvAp+tKbZuGrPdkt4LbBUR/XFH9va0+3d+WC51XiJpqxbLWxqkBGDWkqSPA8PA13odS10i4syI2A74AnBKr+PpNkmrAV8HTup1LD3wA2AoInYFruWNKseYBikBtDO0xO/WkbQGsAHwZC3RdU+pQ2q01W5J+wAnAwdHxEs1xdZN432/LwQO7WZANRmr3dOAXYAFkh4CPgBcMQAdwWO+3xHxZOVv+1zgD9rd+CAlgHaGlrgCmJkfHw5cH7kXpY+VOqTGmO2W9B7gHNLBf3kPYuyGdtq9Q2XyAOBXNcbXLaO2OyKeiYiNI2IoIoZIfT4HR8TC3oTbMe2835tXJg8G7mt7673u5e5wj/n+wC9JveYn53n/QPpDAFgbuBi4H/gZsG2vY66p3e8j1Q5/QzrjuafXMdfU7n8DHgcW558reh1zTe0+Hbgnt/kGYOdex1xHu5vWXcAAXAXU5vv9j/n9viO/3+9od9seCsLMrFCDVAIyM7NxcAIwMyuUE4CZWaGcAMzMCuUEYGZWKCcAM7NCOQHYQJI0JOlFSYt7HUunSPp7SZ8fYdk6eSjglyVtXHds1p+cAGxgSFq9adYDEbHbOJ6vPKZM34mIF3NbH+11LNY/+vKP3axB0sWSzpF0M/DXY6x7maRFku6RNCvPG8o32/g2cDewlaRP5JEV75B0vqT1JF2Vp++W9KeVbX5c0s/yp+9zGkmoeRt53on5+XdL+lxl//dJ+pcc148krVPZ/smSfinpJmDHPG/EeMzGY41eB2A2Se8C5kfEB9pY95MR8VQ+wN4q6Xt5/g7AzIi4WdLOpNEz/0NEPCFpBrAv8GhEHAAgaYP8+53AnwJ7RMQrks4Cjpa0qHkbkv4A+DPg/YCAWyT9O/B03v9REfHnkuYDhwHfyc85EtiN9L96G7BopHjMxstnANa3JK0NzCCNi9KOz0q6gzRQ2FakAy/AwxFxc368F3BxRDwBEBFPAXcBH5F0mqQ/iohn8rp7k0ZevDX3NexNunNTq238IXBpRPwmIp4Hvg/8Ud7OryNicX68iHQHN/LySyPihUg3eWkMAjZSPGbj4gRg/Wxn4JZIN/cZlaQ9gX2AD0bEu4HbSYMDQhokb0QR8UvgvaQD76mS/raxWWBeROyWf3aMiL+fQDuqw1S/xhhn5qPEYzYuTgDWz94F3NnmuhsAT0fEC5LeQRovvpXrgf8kaSOAXL55K/BCRHyHdFOZ9+Z1rwMOl7RpZd1tWm0D+DFwqKR1Ja0HfCzPG82N+TnrSJoGHJS3N1I8ZuPiPgDrZ+8iDevdjmuAT0u6D/gFqQz0eyLiHklfAf5d0mukM4XvAl+T9DrwCukG80TEvZJOAX6Urx56BTg+9yWsso2IOFbStyrxnhsRt0saGingiLhN0kWkYX6Xk8aGb7T79+IxGy8PB20DKR9Yr4yIXXodS52U7oY13Oh/MBuNS0A2qF4DNhikL4KNpvFFMGBN4PUeh2N9wmcAZmaF8hmAmVmhnADMzArlBGBmVignADOzQjkBmJkVygnAzKxQTgBmZoVyAjAzK9T/BxfaKvvh1opgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEaCAYAAAAYOoCaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAgIUlEQVR4nO3debwcVZ338c+XTSIEIibmiYFLUFRkM+IVUXBEUAcXFl8wCC4ERePuw4COgI6DDy6gAjKP6BiFIbgBoiyiOEZ2VNAEQiAEZTGMxEBACBDIgMBv/jinSaXpvrfv0lW3b33fr1e/bnf1qarfqe5bv65zqk4pIjAzs/pZp+oAzMysGk4AZmY15QRgZlZTTgBmZjXlBGBmVlNOAGZmNeUEYCMiaVtJ8yWp4jhOlPThKmMYLkm7S7qr8HqxpN2ri2j4JB0r6fsDvN+zdRuPnAC6SNI7885xlaTlki6WtFtTmUMlhaR35NevzeVXSXokv7eq8JglaUnTMua1mXZUfr5U0hvy8w3yzvKuvLylkr5emG+ppNVN6/zGANU8DvhaRISkZ0k6TdKdkh6WtFDSm5vi2lPSLZIelXSZpC0L731N0q153lskHdI070xJC/K8CyTNLLz9NeAYSRs0zfNqSb/NzyNv01WS7pP0I0mTBqhbcTmDbbfdJP1W0oOS7pf0G0mv7GTZzSJiu4i4fJB4ZuT6rDecdVRlPNetFzkBdImkI4CvA18CpgJ9wDeBfZuKzgLuBw4BiIirImLjiNgY2C6XmVSY9mtgG0lT8nrWA14GTGia9mrgyhahHQ30AzsDE4HdgeuayuzdWF9+fKxNHacBrwfOz5PWA/4CvA7YFPgscI6kGbn8ZOCnwL8CmwHzgbMLi3wE2DvPOws4RdJr8rwbABcA3weeA8wFLmjs8CNiOXALsE9TmG8FflF4/bK8HV+Ql3Nsq7q10Ha7SdoEuAj4/7le04HPA491uGwrkRNLQUT4McoP0g5sFfBPg5TbEngK2B94Avg/Te/PAAJYr2n67cD++fnOwGWkHWJx2qPA+vn1UuAN+flFwOEDxPR02Q7qeQjw60HKLCrENRv4beG9jYDVwDZt5r0QODI/fxOwDFDh/f8G9iq8/gzwn03LuA7YKT8PYOvCex8BftX0uZ0GLM/r+gKw7mDbjZQYVg7h+zEBOAN4ALgZ+BRwV6vPIH+W84GHgHuAkwp1j/w9W0VK+C8ELgX+BtwH/ID046G43E/mz+RBUvLdsPD+vsDCvK7bG9t2oO3Som7HAucAZwIPA4uB/lGo2zqkHxR3Aivy8jdt+i7emev+r03rORY4l/Tj4SHg/XndvwNW5np9A9igsLzI349bcz2Oy9v3t3kZ5xTL9+rDRwDd8WpgQ+C8QcodAsyPiJ8AS4B3dbj8K4F/yM//AbgKuLpp2jUR8fcW814DHCHpI5J2GGHb/Q7AH9u9KWkq8GLSTgDSEc0Njfcj4hHSjma7FvNOAF7ZNO+iyP+d2aKmeZeQjoYay5hGOvq6vsXynwPsR9oeDWeQEvHWwMtJSef9+b2BttufgCclzZX05rzsgfwbaWfyQuAfSUc77ZwCnBIRm+Ty5+Tpjc96UqSjtN8BAr4MPB94KbAFzzzCORDYC9gK2BE4FEDSzqSd6qeASXn5S/M8Z9B+u7SyD3BWXs6FpJ3rSOt2aH68nnT0tnFjuZK2JR1dvwuYRkpY05vWtS8pCUwiJcYngX8GJpP+X/ck7fCL/hF4BbAL8C/AHODdpO26PXDwANugJzgBdMdzgfsi4olByh0C/DA//2F+3YkrWPNP8lpSAriqadoVbeb9MnAC6Z9lPrBMUvMO6HxJKwuPD7RZ1iTSr6NnkLQ+6R9tbkTckidvTPrlWfQgqUml2X+QksV/DWHeh3NMDW8BftmUNK6TtJL0C7kP+HaOd2ouf3hEPBIRK4CTgYPyfG23W0Q8BOxG+tX4HeBeSRfmZbZyIPDFiLg/Iv4C/HubcgB/B7aWNDkiVkXENe0KRsRtETEvIh6LiHuBk0jNcUX/HhF/jYj7gZ8BM/P0w4DT8/xPRcSyiLilg+3SytUR8YuIeBL4HoWkPNy6kbb7SRFxR0SsIjXJHZSbcw4AfhYRV0fE48DnSJ9F0e8i4vxct9URsSAiromIJyJiKel70LytvhIRD0XEYuAm0tHiHRHxIHAxKRn2NCeA7vgbMHmgtkZJu5J+hZ2VJ/0Q2KGpY7OdK4Ed8y/NXUhf7luAaXnabrRu/ycinoyIUyNiV9LO8ovA6ZJeWii2X0RMKjy+0yaOB2ix85a0Dukf/3Gg2H+wCtikqfgmNCURSV8l/cI6sLDz7mTeiaRD+oa3sHb7P6TmoEmkI7RvAVdJ2pDUHLc+sLyR+Eg7hefB4NstIpZExKERsXmO/fmkPqBWnk/qK2m4s005SDvmFwO3SPqDpLe1KyhpqqSzJC2T9BCpyWNyU7G7C88fJSVWSL9qb2+x2AG3SxvN69iwzf9Cx3UjbbPidrqT1Oc0labtGRGPkv4Hi4rbG0kvlnSRpLvztvoSz9xW9xSer27xemN6nBNAd/yO1AG43wBlZpEO2RdKuhu4tjB9QBFxB/BXUpv6f+dfRI31ziZ9MQf6NdVYzuqIOJW0I992sPItLCL9Az8tN42cRvrH3L+pGWoxazfRbEQ69F9cmPZ54M3Am/Iv6+K8OzY1vexYnJfU7HFDXs76pF9081oFnuP6LikJb0/aQTwGTC4kvk0i4hnNU4Ntt5yMz8jLbWU5aYfb0NemHBFxa0QcTNrhngCcm7dbq2F8v5Sn75CbVd5N+o514i+kz6LV9I62y1ANsW5/JSWjhj5Ss9Q9pO25eeON3Hz43ObVNb3+FumkgRflbXUMnW+rccMJoAvyIeLngFMl7Sfp2ZLWz+3DX8m/OA8k7axnFh4fB97Z4VkKVwFH5L8NV+dp8yNidauZJB2udN75BEnr5WaMibRoJ+/APGCnXJ+Gb5F2xHu3iOE8YHtJ++d5Pkdq178lx3Y08E5S513zL7jLSe22n1A63bRxZHFpoczrSIfmkI6CFjUlkadJWhd4L+mX3B2RziL6FXCipE0krSPphZJel8u33W6StpF0pKTNc9ktSO3D7ZLwOcDRkp6T5/l4m3JIerekKRHxFGuObp4C7s1/X1AoPpF0pPSgpOmk9vxOnQa8V+k03XUkTZe0zWDbZSSGWLcfAf8saStJG5OS3dm5mfVcYG9Jr1E6K+xYBt+ZTyR15q6StA3Qk9eQjJQTQJdExImknfFnSV/ov5CaQ84nHRmsBs6MiLsbD+B00mHtXh2s4grSL6erC9OuytNaNv9kjwInkg7T7wM+SvqlfkehzM+09nUALTuzI+Ie0g54XwClc/o/SEpmdxfmf1cufy/pjKcvkn49v4q125K/RPpld1th3mPyvI+TttshpJ3F+0hNVY/ndU8j/Ro/Py+r+fTPhhskrcrrnwW8PbeHk5e9AenMnAdIO5ZpHWy3h3NdrpX0CGnHfxNwZKvtRjpF9E7gz6Sd6/falIP0XVicYz4FOCgfgTxK2o6/yU0zu+Tl7kTqG/k56ZTbjkTE70kJ8eQ8/xWs+cU90HYZiaHU7XTSdrqStN3+h5w4cxv9x0nNqctJSXAFA5+G+0nSj42HSf02Zw9QdtzS2v1jZkOjdAbGXGDnqPDLJOlE4PaI+GZ+fTNwQETcXFVMVo18hLCS1Lzz54rDGdOcAGzcyc0AR0TE8VXHYuWQtDdwCanp50TSEdlOVf4o6QVOAGbW8yR9l3Q6qEin6X4kItpeo2KJE4CZWU25E9jMrKacAMzMaqqnRsWbPHlyzJgxo+owzMx6yoIFC+6LiCnN03sqAcyYMYP58+dXHYaZWU+R1HK4ETcBmZnVlBOAmVlNOQGYmdWUE4CZWU05AZiZ1ZQTgJlZTTkBmJnVlBOAmVlN9dSFYJbsevylLFu5mumTJvCbo/aoOhwz61E+AuhBy1auZunxb2XZypZ3fTQz60hpCUDSupKul3RRfr2VpGsl3Sbp7HwTDzMzK0mZRwD/F1hSeH0CcHJEbE26z+hhJcZiZlZ7pSQASZuTbtL93fxawB6km0tDuqfsfmXEYmZmSVmdwF8H/gWYmF8/F1gZEU/k13cB01vNKGk2MBugr6+vu1GOYY2OX4DpkyZUHI2ZjQddTwCS3gasiIgFknYf6vwRMQeYA9Df31/b+1c2On7NzEZLGUcAuwL7SHoLsCGwCXAKMEnSevkoYHNgWQmxmJlZ1vU+gIg4OiI2j4gZwEHApRHxLuAy4IBcbBZwQbdjMTOzNaq8DuDTwBGSbiP1CZxWYSxmZrVT6pXAEXE5cHl+fgewc5nrNzOzNXwlsJlZTTkBmJnVlBOAmVlNOQGYmdWUh4MeQzzMs5mVyUcAY4iHeTazMjkBmJnVlBOAmVlNOQGYmdWUO4Er5mGezawqTgAV8zDPZlYVNwGZmdWUE4CZWU05AZiZ1ZQTgJlZTTkBmJnVVNcTgKQNJf1e0g2SFkv6fJ5+hqQ/S1qYHzO7HYuZma1RxmmgjwF7RMQqSesDV0u6OL/3qYg4t4QYzMysSdcTQEQEsCq/XD8/otvrNTOzgZVyIZikdYEFwNbAqRFxraQPA1+U9DngEuCoiHisxbyzgdkAfX19ZYRrbTRftewhq816WymdwBHxZETMBDYHdpa0PXA0sA3wSmAz4NNt5p0TEf0R0T9lypQywrU2Glcte8hqs/Gh1LOAImIlcBmwV0Qsj+Qx4D+BncuMxcys7so4C2iKpEn5+QTgjcAtkqblaQL2A27qdixmZrZGGX0A04C5uR9gHeCciLhI0qWSpgACFgIfKiEWMzPLyjgLaBHw8hbT3YM4itxBa2ZD5eGgx4nisNIzjvp5xdGYWS/wUBBmZjXlBGBmVlNOAGZmNeU+gJpodBK7g9jMGnwEUBONTmJfwWtmDU4AZmY15QRgZlZTTgBmZjXlTuAxaPqkCU9fzDWcTtvG/CPp8PWVxWbjnxPAGFTc2Q7nqt7G/CO5IthXFpuNf24CMjOrKScAM7OacgIwM6spJ4Axrtiha2Y2mtwJPMb57Bsz65Yybgm5oaTfS7pB0mJJn8/Tt5J0raTbJJ0taYNux2JmZmuU0QT0GLBHRLwMmAnsJWkX4ATg5IjYGngAOKyEWMzMLOt6AohkVX65fn4EsAdwbp4+l3RjeDMzK0kpncCS1pW0EFgBzANuB1ZGxBO5yF3A9DbzzpY0X9L8e++9t4xwzcxqoZQEEBFPRsRMYHNgZ2CbIcw7JyL6I6J/ypQp3QrRzKx2Sj0NNCJWApcBrwYmSWqchbQ5sKzMWMzM6q6Ms4CmSJqUn08A3ggsISWCA3KxWcAF3Y7FzMzWKOM6gGnAXEnrkhLOORFxkaSbgbMkfQG4HjithFjMzCzregKIiEXAy1tMv4PUH2BmZhXwlcBDNFo3Vy8up0wjvdeAmY0fTgBD1Bgnf6Rj5BfH2y/TSO81YGbjhweDMzOrKScAM7OacgIwM6sp9wGUYLQ6jps1d+gO9f3RNlg9u7UdzGx4nABKMFodx80G24mWvZMdrJ7d2g5mNjxuAjIzqyknADOzmnICMDOrKScAM7OacgIwM6spJwAzs5pyAjAzqyknADOzmvKFYDVWvFK4OK2TC8ga8/qqXrPe1fUEIGkL4ExgKhDAnIg4RdKxwAeAe3PRYyLiF92Ox9ZotePu9Crdxry+qtesd5VxBPAEcGREXCdpIrBA0rz83skR8bUSYjAzsyZl3BJyObA8P39Y0hJgerfXa2ZmAyu1E1jSDNL9ga/Nkz4maZGk0yU9p8xYzMzqrrROYEkbAz8BDo+IhyR9CziO1C9wHHAi8L4W880GZgP09fWVFa6NksYQ0NB6SGoPEW1WnVISgKT1STv/H0TETwEi4p7C+98BLmo1b0TMAeYA9Pf3R/ejtdE02L2PPUS0WXW63gQkScBpwJKIOKkwfVqh2NuBm7odi5mZrVHGEcCuwHuAGyUtzNOOAQ6WNJPUBLQU+GAJsZiZWVbGWUBXA2rxls/5NzOrkK8EHmXjsVOzeNXvUMo3npvZ2OQEMMrGY6fmUBPZeEl8ZuOdB4MzM6spJwAzs5pyAjAzqyn3AdiIuMPXrHc5AdiIuMPXrHe5CcjMrKacAMzMasoJwMyspgbtA5C0ZUTcWUYwY02nV/UONuSxmdlY1Ekn8MWSngfcAiwCbmz8jYiHuhlc1Tq9qnewIY/NzMaiQRNARGwr6VnAtsAOwI7AvsCOkh6LiK26HKOZmXVBR6eBRsRjwPWSbgNWA5OBF5OOBMzMrAd10gfwEuCtwNuAKcA84AfA7Ih4vLvhmZlZt3RyBLAEuB44AbggHw3YMPiqWTMbSzpJAB8Gtgc+CnxD0t9IHcE3kjqCz+9eeOOLr5o1s7Gkk07gbxdfS9qcNZ3B+wPnDzS/pC2AM4GppNs/zomIUyRtBpwNzCDdEvLAiHhgyDUwM7NhGfJYQBFxF3AXcHGHszwBHBkR10maCCyQNA84FLgkIo6XdBRwFPDpocZjZmbD0/UrgSNieURcl58/TOpTmE46lXRuLjYX2K/bsZiZ2RqljgYqaQbwcuBaYGpELM9v3U1qImo1z2xgNkBfX18JUXZmsA5dd/i2N9rbpvlKbPe1mHWmtAQgaWPgJ8DhEfGQpKffi4iQFK3mi4g5wByA/v7+lmWqMNhOxjuh9kZ72xSvxB5P92I267ZSBoOTtD5p5/+DiPhpnnyPpGn5/WnAijJiMTOzpOsJQOmn/mnAkog4qfDWhcCs/HwWcEG3YzEzszXKaALaFXgPcKOkhXnaMcDxwDmSDgPuBA4sIRYzM8u6ngAi4mpAbd7es9vrt/Gl0w5fdwybDc73BLae0mmHrzuGzQbnO4KZmdWUE4CZWU05AZiZ1ZQTgJlZTTkBmJnVlBOAmVlNOQGYmdWUE4CZWU05AdhaGkM1ewhrs/HPVwLbWjxkgll9+AjAzKymnADMzGrKCcDMrKbcB2BjTmMo57KHcfYQ0lY3PgKwMacxlHNjZ1z2eqtYt1kVyrgl5OmSVki6qTDtWEnLJC3Mj7d0Ow4zM1tbGUcAZwB7tZh+ckTMzI9flBCHmZkVdD0BRMSVwP3dXo+ZmQ1NlZ3AH5N0CDAfODIiHmhVSNJsYDZAX1/fqAbgTr/e5quWzUamqgTwLeA4IPLfE4H3tSoYEXOAOQD9/f0xmkH4vrG9zQnbbGQqOQsoIu6JiCcj4ingO8DOVcRhZlZnlSQASdMKL98O3NSurJmZdUfXm4Ak/QjYHZgs6S7g34DdJc0kNQEtBT7Y7TjMzGxtXU8AEXFwi8mndXu9ZmY2MF8JbGZWU04AZmY15QRgZlZTTgBmZjXl4aBtTGhc1dt43jytWK6T5Yz06uDicrp1wZmvRLeqOQHYmNBq5zecHeJo7UQby+nmFeK+Et2q5iYgM7OacgIwM6spJwAzs5pyAjAzqyknADOzmnICMDOrKScAM7OacgIwM6spXwhmtVK8+rZhPF6F26jneKybjR4nAKuV4tW3DePxKtxGPcdj3Wz0dL0JSNLpklZIuqkwbTNJ8yTdmv8+p9txmJnZ2sroAzgD2Ktp2lHAJRHxIuCS/NrMzErU9QQQEVcC9zdN3heYm5/PBfbrdhxmZra2qvoApkbE8vz8bmBqu4KSZgOzAfr6+koI7ZlaDVVs1s5AHbDNQ0B3W/N31x3CVlR5J3BEhKQY4P05wByA/v7+tuW6yf80NhQDdcC26oTupuJ31x3C1qyq6wDukTQNIP9dUVEcZma1VVUCuBCYlZ/PAi6oKA4zs9oq4zTQHwG/A14i6S5JhwHHA2+UdCvwhvzazMxK1PU+gIg4uM1be3Z73WZm1l7lncBjTdlnadj40Op7U6czcDz0RG9yAmhS9lkaNj60+t7U6QwcDz3RmzwaqJlZTTkBmJnVlBOAmVlNuQ/AbACtOnKbO3xHs9NzOMvutAO2URd31FqDE4DZAFp15BY7fEe703M4y+60A7ZRF3fUWoObgMzMasoJwMysppwAzMxqqpZ9AMVOs4ZiB5mNL6N1P4dufkdaLbsYd6vyo92RO9av5u1m53td1TIBDHbVpo0vo/XZdvM70mrZA62vGx25Y/1q3m52vteVm4DMzGrKCcDMrKacAMzMaqo2fQAe5tna8QkA3eXO27Gr0gQgaSnwMPAk8ERE9HdrXR7m2drxDqm73Hk7do2FI4DXR8R9VQdhZlY37gMwM6upqhNAAL+StEDS7IpjMTOrlaqbgHaLiGWSngfMk3RLRFxZLJATw2yAvr6+KmI0a6vTDuTR7mgeraubOzVWO3JHEtdQ5i2WbRjJ+oa7jNFWaQKIiGX57wpJ5wE7A1c2lZkDzAHo7++P0oM0G0Cn/7yj/U9e9k5jrHbkjiSuoczb6iSSkaxvuMsYbZU1AUnaSNLExnPgTcBNVcVjZlY3VR4BTAXOk9SI44cR8csK4zEzq5XKEkBE3AG8rKr1m5nVXdWdwGY2SobTMdzqCvlW90Fut57itMHeK75fXO5wOkc77VQfyRDXg40eMBr3WB5sW3ebE4DZODGcncdgQ6O36qQcaD2DxdDqvsTD6RzttK4jGeJ6sNEDRuMey4Nt626r+joAMzOriBOAmVlNOQGYmdWU+wDMamYoncV1Gyp7OMPGD3Tv5sb7nfRZjNZyhsIJwKxmhrITGStDPpRlOMPGD7aNOu3cHa3lDIWbgMzMasoJwMysppwAzMxqyn0AZla6bg5nPZyrm+vKCcDMStfNHfBwrm6uKzcBmZnVlBOAmVlNOQGYmdWU+wDMbMwZbFjpgeYZa1c3j+WrqZ0AzGzMGU4ncVX3Zx5r6xuKSpuAJO0l6Y+SbpN0VJWxmJnVTZU3hV8XOBV4M7AtcLCkbauKx8ysbqo8AtgZuC0i7oiIx4GzgH0rjMfMrFYUEdWsWDoA2Csi3p9fvwd4VUR8rKncbGB2fvkS4I9DXNVk4L4RhttrXOf6qGO961hnGFm9t4yIKc0Tx3wncETMAeYMd35J8yOifxRDGvNc5/qoY73rWGfoTr2rbAJaBmxReL15nmZmZiWoMgH8AXiRpK0kbQAcBFxYYTxmZrVSWRNQRDwh6WPAfwHrAqdHxOIurGrYzUc9zHWujzrWu451hi7Uu7JOYDMzq5bHAjIzqyknADOzmhoXCWCwISUkPUvS2fn9ayXNqCDMUddBvf9B0nWSnsjXXfS8Dup8hKSbJS2SdImkLauIc7R1UO8PSbpR0kJJV4+Hq+o7HSpG0v6SQlLPnxrawed8qKR78+e8UNL7R7TCiOjpB6kD+XbgBcAGwA3Atk1lPgL8R35+EHB21XGXVO8ZwI7AmcABVcdcUp1fDzw7P/9wjT7rTQrP9wF+WXXc3a5zLjcRuBK4BuivOu4SPudDgW+M1jrHwxFAJ0NK7AvMzc/PBfaUpBJj7IZB6x0RSyNiEfBUFQF2QSd1viwiHs0vryFdX9LrOqn3Q4WXGwG9fnZHp0PFHAecAPxPmcF1SenD44yHBDAd+Evh9V15WssyEfEE8CDw3FKi655O6j3eDLXOhwEXdzWicnRUb0kflXQ78BXgEyXF1i2D1lnSTsAWETFebvLb6fd7/9zEea6kLVq837HxkADMnkHSu4F+4KtVx1KWiDg1Il4IfBr4bNXxdJOkdYCTgCOrjqVkPwNmRMSOwDzWtGwMy3hIAJ0MKfF0GUnrAZsCfysluu6p41AaHdVZ0huAzwD7RMRjJcXWTUP9rM8C9utmQCUYrM4Tge2ByyUtBXYBLuzxjuBBP+eI+FvhO/1d4BUjWeF4SACdDClxITArPz8AuDRyj0oPq+NQGoPWWdLLgW+Tdv4rKoixGzqp94sKL98K3FpifN0wYJ0j4sGImBwRMyJiBqm/Z5+ImF9NuKOik895WuHlPsCSEa2x6p7vUeo9fwvwJ1IP+mfytP9H+kIAbAj8GLgN+D3wgqpjLqneryS1Iz5COuJZXHXMJdT518A9wML8uLDqmEuq9ynA4lzny4Dtqo6523VuKns5PX4WUIef85fz53xD/py3Gcn6PBSEmVlNjYcmIDMzGwYnADOzmnICMDOrKScAM7OacgIwM6spJwAzs5pyArBxSdIMSaslLaw6ltEi6VhJn2zz3oQ8PPDjkiaXHZv1JicAGzckrds06faImDmE+ZXHmOk5EbE61/WvVcdivaMnv+xmDZJ+LOnbkq4Bjh6k7PmSFkhaLGl2njYj34DjTOAmYAtJh+TRFm+Q9D1JG0n6eX59k6R3FJb5bkm/z7++v91IQs3LyNOOyPPfJOnwwvqXSPpOjutXkiYUlv8ZSX+SdDXwkjytbTxmQ7Fe1QGYjdAOwDkRsUsHZd8XEffnHewfJP0kT38RMCsirpG0HWkkzddExH2SNgP2Av4aEW8FkLRp/vtS4B3ArhHxd0nfBN4laUHzMiS9Angv8CpAwLWSrgAeyOs/OCI+IOkcYH/g+3meg4CZpP/V64AF7eIxGyofAVjPkrQhsBlprJROfELSDaSBw7Yg7XgB7oyIa/LzPYAfR8R9ABFxP3Aj8EZJJ0h6bUQ8mMvuSRqN8Q+5r2FP0t2cWi1jN+C8iHgkIlYBPwVem5fz54hYmJ8vIN3Jjfz+eRHxaKQbvjQGBmsXj9mQOAFYL9sOuDbSTX4GJGl34A3AqyPiZcD1pEECIQ2W11ZE/AnYibTj/YKkzzUWC8yNiJn58ZKIOHYY9SgOWf0kgxyZDxCP2ZA4AVgv2wFY1GHZTYEHIuJRSduQxo9v5VLgnyQ9FyA33zwfeDQivk+6wcxOuewlwAGSnlcou2WrZQBXAftJerakjYC352kDuTLPM0HSRGDvvLx28ZgNifsArJftQBreuxO/BD4kaQnwR1Iz0DNExGJJXwSukPQk6UjhR8BXJT0F/J10s3ki4mZJnwV+lc8e+jvw0dyXsNYyIuJQSWcU4v1uRFwvaUa7gCPiOklnk4b+XUEaL75R72fEYzZUHg7axqW8Y70oIravOpYyKd0dq7/R/2A2EDcB2Xj1JLDpeLoQbCCNC8GA9YGnKg7HeoSPAMzMaspHAGZmNeUEYGZWU04AZmY15QRgZlZTTgBmZjXlBGBmVlNOAGZmNeUEYGZWU/8L0Mq4Sv/b5i8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEaCAYAAAAYOoCaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAcsElEQVR4nO3de7xcZX3v8c+Xm4SIhJAYISRsLAiEO24uFrAURNEgpAeK+KISPGhqi1ULPRirp4da7Al6rHIO2EMEJAoCgQpEPdwOFxFakATCNVwChpJASIBsrimX8Osf69nJYmcme3b2njV75vm+X6/9yqzLrPk9s3fWb57fetYzigjMzCw/G7Q6ADMzaw0nADOzTDkBmJllygnAzCxTTgBmZplyAjAzy5QTQAeTdJGkM9PjQyQtbnVM7UTSrZK+kB6fIOmGVse0viSFpB3qbGvrttn6cwLoAOlEtULSexrc/xuSru2z7vE6645Pj1efQCSNknShpKWSXpH0mKTppeeFpNckvVr6Ob3Psa+X9HFJZ0h6q7TfAknHDKDtB0n6V0kvSXpR0h2S9k3bNpH0fUmL07EXSfpho8cui4hLIuLjDcSzOum2i05um62bE0Cbk9QFHAwEcFSDT7sN+ENJG6ZjbA1sDOzdZ90Oad++fgC8F9gF2CK97sI+++wZEe8t/Xy3FPNIoBv4TVp1ee9+wNeAiyWN668Rkt4H/Ar4P8BoYDzw98AbaZdvpNfZD9gcOAS4p7/jWvUkbdTqGHLkBND+TgTuBC4Cpjb4nLspTvh7peWDgVuAR/useyIinqnx/H2Bn0fEioh4JyIeiYgrBxDzYcAdEfFG3w0RcT3wCvAHveskHSlpvqSe9Gl/j7TpQ+k5l0bEqohYGRE3RMT9pTiviohnorAoIn5aLyhJh0t6JPUmzgFU2naSpNvTY0n6gaRlkl6W9ICk3SRNA04ATk89jl+m/adLeiL1lh6W9Cd9jyvpf6Ve3O8lfbK0fbSkn0h6Jm2/uoH3pZ6PpV5dj6RzJWkI2rZL6oH2SHpI0uoPIZK2kvTLdJy7JZ3Z+zppe0g6RdLjwONp3dmSnk7PmSfp4NL+Z0i6QtLF6b18QNKHVPRol6Xn9duTsTWcANrficAl6ecTjXxyjog3gbuAj6ZVHwV+C9zeZ12tT/9QJJzvSPq8pB3XI+ZPAb/uuzKdfCYDmwAPp3V7AxcCfw5sBZwHzFFR7noMWCVplqRPStqyRpynSvpLSbv3nvBqkTQG+AXwLWAM8ARwYJ3dP07x/nyIogd0HPBCRMyk+D18N/VoPp32f4IioW5B0UO5OPWweu1PkXzHAN8FLijF+jNgM2BX4P0Uva/+3pd6jqRIinukmD8xmLZJ2hj4JXBDiu2vgEsk7ZSOdS7wGvABig8ntT6gTEntn5SW76b4EDIa+DlwhaRNS/t/Or0nWwL3AtdTnMfGA99O74M1KiL806Y/wEHAW8CYtPwI8Nel7RcBZ6bHhwCLS9vOoPh0DHAfsCNwRJ91U0v7B7BDejwC+FtgXnr9hcAn++z7MtBT+vlEafu/AxNKcbyZ9nkNWAWcXtr3n4F/6NPuR4E/So93Se1cDLwNzAHGpW0bAqcAd1CUhZ4pt6nPMU8E7iwtKx3zC2n5JOD29PhQiuRzALBBn+Osfs/X8XubDxxdOu7C0rbN0vv3AWBr4B1gyxrHWOf7UmP/AA4qLc8Gpg+mbRRJbWl5P+DS9DvdMP1t7FTadmbv65RiOrSf92oFRTmx92/lxtK2TwOvAhum5c3TMUe1+v9mu/y4B9DepgI3RMTzafnnNF4Gug04SNJoYGxEPA78K8W1gdHAbtTpAURRavnHiPgwxafP2RSf1EaXdtsnIkaVfq4HkLQ78FJEPF3ad3baZyRF6edESX+etm0HnJZKDD2SeoAJwDYplgURcVJEbJti3gb4Ydq2KiLOjYgDgVHAd4ALJe1So1nbAKtjiuKM8nSN/YiIm4FzKD7hLpM0U8X1iJoknVgq1fSkOMeUdllaOvbr6eF7UztfjIgVNQ67zveljqWlx6+n1xhM27YBno6Id0rrnqL4ND4W2Ih3v4e13s93rZP0NyoGAryU2rQF736vnis9Xgk8HxGrSsvUapfV5gTQpiSNoOie/5GK0ThLgb8G9pS0ZwOH+DeK/1xfpPiETES8TPEp+YvAMxHx+/4Okp7zj8BIYPsGXvdTwP9bx/EWAddSfLqD4gTxnT7JZLOIuLTGcx+h+JS6W41tKyPiXIpPlJP6bgeepTiBAkU5qrxc43j/OyXASRTlkv/Wu6m8n6TtgB8DXwa2iohRwIOUri+sw9PAaEmj6mxr6H0ZqEbbRvG3MkFS+TwyEVgCLKfokW1b2lbr/Vx9zFTvP53i73rL9F69RGPvla0HJ4D2NYWiXDKJoma6F0U55LcU5Yx1ioiVwFzg1PScXrendfXq/0j675L2VTHMclPgqxQlnEcbiLtm/b907G0pSlEPpVU/Br4kaf90jWCkpMmSNpe0s6TT0nOQNAH4LEXtH0lfU3H/wwhJG0maSlEmuLfGS/8a2FXSf1ExIuUrFGWYWjHum+LZmKJs9R8UpRooPqF+sLT7SIqT3PL03M9TI0HVEhHPUiTDH0naUtLGknqv0dR9Xxo5dj0DbNtdFD2J01Nsh1Ak7svSp/JfAGdI2kzSzvT/d7k5RdJYDmwk6e+Auj0rGzwngPY1FfhJRPx7RCzt/aHovp+gxobV/Ybi4t3tpXW/TevqJgCKE9pPgOcpPgUeDkyOiFdL+9ynd98H8MP0SXYSRamp7DO9+1FcBLyD4mIpETGXokdyDsWn94UUNWsoRgvtD9wl6TWKE/+DwGlp++vA9ylKH89TXA84JiKeXKtBRRntT4EZwAsU10TuqNP+91GcgFdQlDxeAL6Xtl0ATEplmasj4uEUw79RnEB3X8dxa/kcRS39EWAZxTDZ/t6XwRhI296kOOF/kuL9/RFwYuqJQdHr2YLi/f8ZxfWBtUZ+lVwPXEdxDeIpiuRTswxnQ0NFqdOs+SQdBxwbEce1OharnqSzgA9ERKPXqazJ3AOwKvWQhjFa50sluj1SiWo/4GTgqlbHZWv47jurTER4vpm8bE5R9tmGovz1feCalkZk7+ISkJlZplwCMjPLlBOAmVmm2uoawJgxY6Krq6vVYZiZtZV58+Y9HxFj+65vqwTQ1dXF3LlzWx2GmVlbkfRUrfUuAZmZZcoJwMwsU04AZmaZcgIwM8uUE4CZWaacAMzMMuUEYGaWKScAM7NMtdWNYJ3gwBk3s6Sn+OrS8aNGcMf0Q1sckZnlygmgYkt6VrJoxmQAuqbX/WZEM7OmcwnIzCxTTgBmZplyAjAzy5QTgJlZppwAzMwy5QRgZpYpJwAzs0xVkgAkjZJ0paRHJC2Q9BFJoyXdKOnx9O+WVcRiZmaFqnoAZwPXRcTOwJ7AAmA6cFNE7AjclJbNzKwiTU8AkrYAPgpcABARb0ZED3A0MCvtNguY0uxYzMxsjSp6ANsDy4GfSLpX0vmSRgLjIuLZtM9SYFwFsZiZWVJFAtgI2Af454jYG3iNPuWeiAggaj1Z0jRJcyXNXb58edODNTPLRRUJYDGwOCLuSstXUiSE5yRtDZD+XVbryRExMyK6I6J77NixFYRrZpaHpieAiFgKPC1pp7TqMOBhYA4wNa2bClzT7FjMzGyNqqaD/ivgEkmbAE8Cn6dIPrMlnQw8BRxXUSxmZkZFCSAi5gPdNTYdVsXrm5nZ2nwnsJlZppwAzMwy5QRgZpYpJwAzs0w5AZiZZcoJwMwsU04AZmaZcgIwM8uUE4CZWaacAMzMMuUEYGaWqaomg8vWgTNuZknPytXL40eNaGE0ZmZrOAE02ZKelSyaMbnVYZiZrcUlIDOzTDkBmJllygnAzCxTTgBmZpnyReBhrjyKaPyoEdwx/dAWR2RmncIJYJgrjyLqmv7rFkdjZp3EJSAzs0w5AZiZZcoloCHiWr2ZtRsngCHiWr2ZtRuXgMzMMlVJD0DSIuAVYBXwdkR0SxoNXA50AYuA4yJiRRXxdIK+k8z1cvnJzBpVZQnojyPi+dLydOCmiJghaXpa/nqF8bS1epPMufxkZo1qZQnoaGBWejwLmNK6UMzM8lNVDyCAGyQFcF5EzATGRcSzaftSYFytJ0qaBkwDmDhxYhWxtsS6SjpmZs1QVQI4KCKWSHo/cKOkR8obIyJSclhLShYzAbq7u2vu0wn8vQFmVrVKSkARsST9uwy4CtgPeE7S1gDp32VVxGJmZoWm9wAkjQQ2iIhX0uOPA98G5gBTgRnp32uaHUtVxo8asfpirEs4ZjZcVVECGgdcJan39X4eEddJuhuYLelk4CnguApiqYSHYZpZO2h6AoiIJ4E9a6x/ATis2a9vZma1+U5gM7NMeS6gFvK1AjNrJSeAFvK1AjNrJZeAzMwy5QRgZpYpJwAzs0w5AZiZZcoJwMwsU04AZmaZcgIwM8uUE4CZWaacAMzMMuUEYGaWKScAM7NMOQGYmWXKCcDMLFNOAGZmmfJ00B3swBk3s6RnJVB834CnnzazMieADrakZyWLZkwGWP3FM2ZmvVwCMjPLlHsAbcRfIWlmQ8kJoI24hm9mQ8klIDOzTFXWA5C0ITAXWBIRR0raHrgM2AqYB3wuIt6sKp715ZE1ZtYpquwBfBVYUFo+C/hBROwArABOrjCW9dY7smbRjMmrE4GZWTuqJAFI2haYDJyflgUcClyZdpkFTKkiFjMzK1TVA/ghcDrwTlreCuiJiLfT8mJgfEWxmJkZFVwDkHQksCwi5kk6ZD2ePw2YBjBx4sShDS4jfYeQru+1C18DMescVVwEPhA4StKngE2B9wFnA6MkbZR6AdsCS2o9OSJmAjMBuru7o4J4O1L5RD2Yu4J9d7FZ52h6CSgivhER20ZEF3A8cHNEnADcAhybdpsKXNPsWMzMbI1W3gfwdeBUSQsprglc0MJYzMyyU+mdwBFxK3BrevwksF+Vr29mZmv4TmAzs0x5LqAO08iEceV9+q73qB6zfDgBdJhGTuD19vGoHrO8uARkZpYpJwAzs0w5AZiZZcoJwMwsU04AZmaZ8iggq8mTvpl1PicAq8mTvpl1PpeAzMwy5QRgZpYpJwAzs0z1mwAkbVdFIGZmVq1GLgJfK+n9wCPA/cADvf9GxMvNDG64a2TiNTOz4arfBBARkyS9B5gE7A7sARwN7CHpjYjYvskxDlseGmlm7ayhYaAR8QZwb/r2rpXAGOBDFD0BMzNrQ/0mAEk7AZOBI4GxwI3AJcC0iHizueGZmVmzNNIDWADcC5wFXJN6A2Zm1uYaSQB/AewGnAKcI+kFigvBD1BcCL66eeGZmVmzNHIR+LzysqRtWXMx+Bjg6qZE1mKeC8fMOt2A5wKKiMXAYuDaoQ9n+PBcOGbW6XwnsJlZppwAzMwy1fQEIGlTSb+TdJ+khyT9fVq/vaS7JC2UdLmkTZodi5mZrVFFD+AN4NCI2BPYCzhC0gEUw0p/EBE7ACuAkyuIxczMkqYngCi8mhY3Tj8BHApcmdbPAqY0OxYzM1ujkm8Ek7QhMA/YATgXeALoiYi30y6LgfF1njsNmAYwceLE5gdbgyd9M7NOVEkCiIhVwF6SRgFXATsP4LkzgZkA3d3d0ZQA++F7AMysE1U6CigieoBbgI8AoyT1JqBtgSVVxmJmlrsqRgGNTZ/8kTQCOJxifqFbgGPTblOBa5odi5mZrVFFCWhrYFa6DrABMDsifiXpYeAySWdSTDZ3QQWxmJlZ0vQEEBH3A3vXWP8ksF+zX9/MzGrzncBmZplyAjAzy5QTgJlZppwAzMwyVcmNYNYefMezWV6cAGw13/FslheXgMzMMuUegK23viUj9yDM2osTgK238gnf35ts1n5cAjIzy5QTgJlZppwAzMwy5QRgZpYpJwAzs0w5AZiZZcoJwMwsU04AZmaZcgIwM8uUE4CZWaacAMzMMuW5gKypDpxxM0t6VgKeMM5suHECsKZa0rOSRTMmA54wzmy4cQnIzCxTTU8AkiZIukXSw5IekvTVtH60pBslPZ7+3bLZsZiZ2RpV9ADeBk6LiEnAAcApkiYB04GbImJH4Ka0bGZmFWl6AoiIZyPinvT4FWABMB44GpiVdpsFTGl2LGZmtkalF4EldQF7A3cB4yLi2bRpKTCuznOmAdMAJk6cWEGU1iyNfoWkRw6ZVaOyBCDpvcC/AF+LiJclrd4WESEpaj0vImYCMwG6u7tr7mPtodGvkPTIIbNqVDIKSNLGFCf/SyLiF2n1c5K2Ttu3BpZVEYuZmRWqGAUk4AJgQUT8U2nTHGBqejwVuKbZsZiZ2RpVlIAOBD4HPCBpflr3t8AMYLakk4GngOMqiGXt4FxvNrNMNT0BRMTtgOpsPqzZr98f15vNLFe+E9jMLFOeC6ik7zBFq07fUlyvRoeOmtnAOQGU+OTSOuVSXFmjQ0fNbOBcAjIzy5R7ANavgZbG6pVzhlIjo7fK+5S5lGRWcAKwfg30ZFmvnDOUGhm9VS8Ol5LMCi4BmZllyj0AGxIDLROV92/0OVXyDYKWAycAGxIDPUEO9xOqbxC0HLgEZGaWKfcAzDqUy1jWHycAsw7lMpb1xyUgM7NMOQGYmWUqyxJQFXeq2tBrxmR967pb2KzTZZkAqrhT1YZeMy5i+m/BcuYSkJlZprLsAVhn8XcGmK0fJwBre/7OALP14xKQmVmm3AOw7PirP9fw3cJ5cwKw7Pgkt4bvFs6bS0BmZplqeg9A0oXAkcCyiNgtrRsNXA50AYuA4yJiRbNjMatas0osLt3YUKiiB3ARcESfddOBmyJiR+CmtGzWcXpLLItmTK55x/FwO67lpekJICJuA17ss/poYFZ6PAuY0uw4zMzs3Vp1DWBcRDybHi8FxrUoDjOzbLX8InBEBBD1tkuaJmmupLnLly+vMDIzs87WqgTwnKStAdK/y+rtGBEzI6I7IrrHjh1bWYBmZp2uVQlgDjA1PZ4KXNOiOMzMslXFMNBLgUOAMZIWA/8DmAHMlnQy8BRwXLPjsDw0+y7fet8f0DeG4TYsc7hNmDfYYayNPN9DZfvX9AQQEZ+ts+mwZr+25afZ/8kb+f6A4XhH7XCbMG+wdyA38nzf5dy/ll8ENjOz1nACMDPLlBOAmVmmnADMzDLl6aDN+jFUI4vKx+m7vlmTxFlhoKOgBjqCqF1HHDkBmPVjqP4z1zvOUI5QaWSUUo4GOgpqoCOI2nXEkUtAZmaZcgIwM8tUNiUg10atHTVSW26X7zj2/8H+VX0tIZsE4NqotaNGasvtcsHR/wf7V/W1BJeAzMwylU0PwKwqgynJ9J1srhmlknrxNTJUshklikaGxzb6vrjMNDBOAGZDbDAnxSrKJPXia2SoZDNKFI0Mj230fXGZaWBcAjIzy5R7AGYt1mjJqFWjfeqVVeqVjIaqTDRcRjc12s525ARg1mKNniBbNdqnXlmlXsloqMpEw2V0UyPtbFcuAZmZZco9ADNrawMtFQ1mUr7hUpYaKk4AZtbWBloqGsykfMOlLDVUXAIyM8uUE4CZWaZcAjKztQymrr6uoaLDWZVf/LM+k/w1o/zkBGBmaxmquno71cyr/OKfgU7y16yJ4VwCMjPLVEt7AJKOAM4GNgTOj4gZrYzHzKxqrSyTtSwBSNoQOBc4HFgM3C1pTkQ83KqYzMyq1soyWStLQPsBCyPiyYh4E7gMOLqF8ZiZZUUR0ZoXlo4FjoiIL6TlzwH7R8SX++w3DZiWFncCHh3gS40Bnh9kuO0mxzZDnu12m/Mw2DZvFxFj+64c9qOAImImMHN9ny9pbkR0D2FIw16ObYY82+0256FZbW5lCWgJMKG0vG1aZ2ZmFWhlArgb2FHS9pI2AY4H5rQwHjOzrLSsBBQRb0v6MnA9xTDQCyPioSa81HqXj9pYjm2GPNvtNuehKW1u2UVgMzNrLd8JbGaWKScAM7NMdUwCkHSEpEclLZQ0vcb290i6PG2/S1JXC8IcUg20+aOS7pH0drrvou010OZTJT0s6X5JN0narhVxDrUG2v0lSQ9Imi/pdkmTWhHnUOqvzaX9jpEUktp+aGgDv+eTJC1Pv+f5kr4wqBeMiLb/obiI/ATwQWAT4D5gUp99/hL4v+nx8cDlrY67gjZ3AXsAPwWObXXMFbX5j4HN0uO/aPff8wDa/b7S46OA61odd7PbnPbbHLgNuBPobnXcFfyeTwLOGarX7JQeQCPTShwNzEqPrwQOk6QKYxxq/bY5IhZFxP3AO60IsAkaafMtEfF6WryT4v6SdtdIu18uLY4E2n10R6NTxfwDcBbwH1UG1ySVT4/TKQlgPPB0aXlxWldzn4h4G3gJ2KqS6JqjkTZ3moG2+WTg2qZGVI2G2i3pFElPAN8FvlJRbM3Sb5sl7QNMiIjmTJZfvUb/vo9JJc4rJU2osb1hnZIAzN5F0p8B3cD3Wh1LVSLi3Ij4A+DrwLdaHU8zSdoA+CfgtFbHUrFfAl0RsQdwI2uqGuulUxJAI9NKrN5H0kbAFsALlUTXHDlOpdFQmyV9DPgmcFREvFFRbM000N/1ZcCUZgZUgf7avDmwG3CrpEXAAcCcNr8Q3O/vOSJeKP1Nnw98eDAv2CkJoJFpJeYAU9PjY4GbI11VaVM5TqXRb5sl7Q2cR3HyX9aCGJuhkXbvWFqcDDxeYXzNsM42R8RLETEmIroioovies9RETG3NeEOiUZ+z1uXFo8CFgzqFVt95XsIr6B/CniM4ir6N9O6b1P8UQBsClwBLAR+B3yw1TFX0OZ9KeqIr1H0dh5qdcwVtPn/A88B89PPnFbHXFG7zwYeSm2+Bdi11TE3u8199r2VNh8F1ODv+X+m3/N96fe882Bez1NBmJllqlNKQGZmNkBOAGZmmXICMDPLlBOAmVmmnADMzDLlBGBmliknAOtIkrokrZQ0v9WxDBVJZ0j6mzrbRqTpgd+UNKbq2Kw9OQFYx5C0YZ9VT0TEXgN4vtIcM20nIlamtj7T6lisfbTlH7tZL0lXSDpP0p3AN/rZ92pJ8yQ9JGlaWteVvoDjp8CDwARJJ6bZFu+T9DNJIyX9Oi0/KOkzpWP+maTfpU/f5/Umob7HSOtOTc9/UNLXSq+/QNKPU1w3SBpROv43JT0m6XZgp7SubjxmA7FRqwMwG6TdgdkRcUAD+/7XiHgxnWDvlvQvaf2OwNSIuFPSrhQzaf5hRDwvaTRwBPBMREwGkLRF+ncX4DPAgRHxlqQfASdImtf3GJI+DHwe2B8QcJek3wAr0ut/NiK+KGk2cAxwcXrO8cBeFP9X7wHm1YvHbKDcA7C2JWlTYDTFXCmN+Iqk+ygmDptAceIFeCoi7kyPDwWuiIjnASLiReAB4HBJZ0k6OCJeSvseRjEb493pWsNhFN/mVOsYBwFXRcRrEfEq8Avg4HSc30fE/PR4HsU3uZG2XxURr0fxhS+9E4PVi8dsQJwArJ3tCtwVxRf8rJOkQ4CPAR+JiD2BeykmCIRisry6IuIxYB+KE++Zkv6u97DArIjYK/3sFBFnrEc7ylNWr6Kfnvk64jEbECcAa2e7A/c3uO8WwIqIeF3SzhTzx9dyM/CnkrYCSOWbbYDXI+Jiii+Y2SftexNwrKT3l/bdrtYxgN8CUyRtJmkk8Cdp3brclp4zQtLmwKfT8erFYzYgvgZg7Wx3iqm9G3Ed8CVJC4BHKcpAa4mIhyR9B/iNpFUUPYVLge9Jegd4i+LL5omIhyV9C7ghjR56CzglXUt41zEi4iRJF5XiPT8i7pXUVS/giLhH0uUUU/8uo5gvvrfda8VjNlCeDto6Ujqx/ioidmt1LFVS8e1Y3b3XH8zWxSUg61SrgC066Uawdem9EQzYGHinxeFYm3APwMwsU+4BmJllygnAzCxTTgBmZplyAjAzy5QTgJlZppwAzMwy5QRgZpYpJwAzs0z9J2KmJ6tQdDDfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEaCAYAAAAYOoCaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAawUlEQVR4nO3de7RcZZnn8e+PhGAEJECQwdwONghyE+kgOGjLcOlR0gKzYBQbneBgR1FHeul0i8q0iG1PaGe0XUtnSbyMQbyAKBelUWkg3EZQYrgFBLlKwjWQcI1I4Jk/9nvIpjiVU5Vz9t61z/v7rFXr7Fvt/by76uyn9vvW+5YiAjMzy88mTQdgZmbNcAIwM8uUE4CZWaacAMzMMuUEYGaWKScAM7NMOQGYbYCkUySdmaZnS3pK0qSm49oYkpZI+kCXda0um20cJ4AWkxSSdupY9uIFa5TnvlnSk+V/eEnf6LLs62n6JRcQSZ+WdHe6cKyQdFZp3RJJf0zrhh8/7YjhU5L+SdKBkl4obbdS0uf6OA8zJf1Y0ipJj0u6WdJxpfXHS/pdKttDkv5V0pa97n9YRPwhIraIiOdHiec4SVf1u/8mTeSyWXdOAPm6juL136e07K3Aio5lfwFc0flkSfOB9wGHRMQWwFzgko7NPpouKsOPd3asnwf8a5q+f3g74C3A8ZKO7LEs3wXuA+YA26a4Hkpxvg34J+A9EbEl8HrgrC77sQap4GtSjXyyJ7D0yXqFpE9IeljSA5LeDxARzwHXUFzgkfRqYApwdsey1zFCAgD2BX4REXem/T0YEYv6iG3rtO9fda6LiLuB/wfsVtp+V0kXS3pM0m2S3tURy3ci4umIWBcRyyLiotK6X0XEsrTvxyJicUQ82SWuHSVdnu4WLgaml9YNpbuuyWn+OEl3pW3vlnSspNcDXwfenO5m1qRt50laJukJSfdJOmWE/c6X9Id0J/OZ0vpJ6W7rznSspZJm9XBeRjJH0tVpP7+UNH0cyraVpDMkPSLpXkknD1/IU+z/O5Xpbkkf7TjOEklfkHQ18AzwWknvl3RrOvZdkj5YOhfD7+m/L72nj5R0mKTb03n49CjnwIZFhB8tfQAB7NSx7BTgzDR9ILAOOBXYFDiM4p9s67T+s8D5afpo4Azg0I5ld5X2vQT4QJp+L/AY8HcUn/4ndcTx4rZdYj8G+EEpzhWldTsDK4GD0vzmFJ/w3w9MBt4IrAJ2S+v/Dbg67XN2x3HeCqwFPgccAGw2yjn9FfAlYDOKRPhk6XwOpXM+OcX0BLBLWrcDsHuaPg64qmO/BwJ7Unzo2oviDuXIjv1+A5gKvAF4Fnh9Wv93wE3ALoDS+m1HOy8jlG0JcCdF4p2a5heOQ9nOAM4Htkz7uR04Pq37EHALMBPYOr1WAUwuxfQHYPd07E0p7gz/LJX1bRTv2X063tP/kLb9G+AR4Pvp+Lun13vHpv8/2/BoPAA/xvDi9ZYA1g7/s6VlDwP7l9Y/mv7RvpL+mbZIF6fhZf+39NwllC7qwLHpH/rptJ9Pdmz7DLCm9Ph8af13gfeV4nghbfNEKtdPgClp/buBKzvKeTrw2TS9NbAQWA48D1wP7Fva9h3AT9P+n6K4wE8a4XzOTheXzUvLvk/3BLAGOAqY2rGf4+i4SI5wrH8Bvtyx35ml9b8GjknTtwFHjLCPDZ6XEbZfApxcmv8w8POxlA2YBPyJUtIBPggsSdOXAh8srTuElyeAU0c5V+cBJ3a8pyel+S3T/vYrbb+UlFz92PDDVUDt9jzFp6CyTYHnSvOPRsS60vwzFBd5KKqAtgD2oPi0e2VEPEXxqXJ42UjVPwBExPci4hBgGsUnvc9L+o+lTT4WEdNKj/8BkKoHDgV+Xtr2/rTNq9L+1gKL07o5wH6S1gw/KJLPv0txrI6IkyJid2B7igRwniSl9RdF0f6wDXAExUVspG/DvAZYHRFPl5bd26XsT1NcgD8EPCDpQkm7djtXkvaTdFmqJnk8PW96x2YPlqbLr9Msik/unTZ4XrrodoyNLdt0ivdc+TzdC8xI06+heD8NK0+PuEzSOyRdk6pz1lDcuZbP1aOxvrF6bfr7UGn92pHKZS/nBNBuf6D45Fa2I10uWp0i4o/Ab4B3AjtExO/SqivTsr3YQAIo7ee5iPgRcCNF4hjNvsC9EfFIl/09TvHJe7jR+D7g8o5kskVEnDDCc1cB/4viwrNNx7oXIuISik+lI8X5ALC1pM1Ly2Z3K0RE/CIiDqWoIvkdRRUOFJ9IO30fuACYFRFbUdSlq9u+O9xHUSUy0vKezku/+ijbKooPHHNKy2ZTVOFBcU5nltbNGulwwxOSNgN+TPEabh8R0yi+KNDrubI+OAG021nAySq+BrmJpEMoLprn9LGPK4ATKRpdh12Vlj0QqZG3U2oknCdpy3Tsd1DUv17bwzEPAy7stlLSFhT1+cvTop8Br5P0Pkmbpse+qVESSadJ2kPSZBVf7zwBuCMiHpV0hKRjJG2twpso6pWv6TxuRNxL8e2oz0maIuktrE9CnTFun/a9OUV9/VMU1VhQfBqdKWlK6SlbAo9FxB9TDH/dw3ka9k2Ku6udUxn2krTtaOdlY/VTtvRJ/GzgC+m9MAf4ODD8VeSzgRMlzZA0DfjkKIefQtH+8giwLr2v/nIs5bHunADa7VSKC/dVwGrgn4FjI+LmPvZxOfDqtI9hV6VlV27geU8An6a4C1mTjn1CRJT381W9tB/A0rS8/PXPYa8Z3o7iDmYbiuoMovjGzl9SJIX7KaoxTqO4UAC8Ejg3xXEXxafRw9O61RRtG79PMZ8JfDEivtelXH8N7EfRwP1ZigbOkWxCcaG7P237NorEA8UdxnLgQUmr0rIPA6dKepKiAfPsLvsdyZfS9r9MZfgWRd38aOdlY/Vbtv9G0Q50F8V75/vAt9O6b6S4bwSWUbzu6yiqL18mleljqbyrKV6PC8ZYHutCEf5BGKuPpO0pLgQzwm++7KRP9F+PiDmjbmyV8x2A1W0r4BO++OdB0tT0Hf3JkmZQ3FWd23RcVvAdgJlVRtIrKaoZd6X4ds6FFF/pfKLRwAxwAjAzy5argMzMMuUEYGaWqclNB9CP6dOnx9DQUNNhmJm1ytKlS1dFxHady1uVAIaGhrjuuuuaDsPMrFUkjTg6gKuAzMwy5QRgZpYpJwAzs0w5AZiZZcoJwMwsU04AZmaZcgIwM8uUE4CZWaZa1RHMmnHAwktZuab46dUZ06Zy9UkHNRyRmY0HJwAb1co1a7ln4TwAhk7q+kuOZtYyrgIyM8uUE4CZWaacAMzMMuUEYGaWKScAM7NMOQGYmWXKCcDMLFNOAGZmmXICMDPLlBOAmVmmnADMzDLlsYBaygO0mdlYOQG0lAdoM7OxchWQmVmmnADMzDLlBGBmliknADOzTDkBmJllygnAzCxTTgBmZplyAjAzy5QTgJlZpmpLAJImSVom6WdpfkdJ10q6Q9JZkqbUFYuZmdV7B3AicGtp/jTgyxGxE7AaOL7GWMzMsldLApA0E5gHfDPNCzgIOCdtshg4so5YzMysUNcdwL8Afw+8kOa3BdZExLo0vwKYMdITJS2QdJ2k6x555JHKAzUzy0XlCUDSXwEPR8TSjXl+RCyKiLkRMXe77bYb5+jMzPJVx3DQBwCHSzoMeAXwKuArwDRJk9NdwExgZQ2xmJlZUvkdQER8KiJmRsQQcAxwaUQcC1wGHJ02mw+cX3UsZma2XpP9AD4JfFzSHRRtAt9qMBYzs+zU+otgEbEEWJKm7wLeVOfxzcxsPfcENjPLlBOAmVmmnADMzDJVaxuAvdQBCy9l5Zq1AMyYNpWrTzqo4YjMLCdOAA1auWYt9yycB8DQSRc2HI2Z5cZVQGZmmXICMDPLlKuABsSMaVNfrAZye4CZ1cEJYECUL/huDzCzOrgKyMwsU04AZmaZcgIwM8uU2wAGUBUNwt06nY2lM5o7spm1mxPAAKqiQbhbp7OxdEZzRzazdnMVkJlZppwAzMwy5SqgATeROoi5zcBssDgBDLiJ1EHMbQZmg8VVQGZmmXICMDPLlBOAmVmm3AbQIp0NwmZmY+EE0CL+1oyZjSdXAZmZZcoJwMwsU64CshF1dtoys4nHCcBGVO60ZWYTk6uAzMwy5QRgZpYpJwAzs0y5DaBmg9C4OpYOZe6MZjZxOAHUbBAaV8fSocyd0cwmDlcBmZllygnAzCxTrgKqQFt/+Wo82wbaUmaznDkBVKCtv3w1Xm0DbSqzWc4qrwKS9ApJv5Z0g6Tlkj6Xlu8o6VpJd0g6S9KUqmMxM7P16mgDeBY4KCLeAOwNvF3S/sBpwJcjYidgNXB8DbGYmVlSeQKIwlNpdtP0COAg4Jy0fDFwZNWxmJnZerW0AUiaBCwFdgK+BtwJrImIdWmTFcCMLs9dACwAmD17dvXBttCgNcAOWjxmNrJaEkBEPA/sLWkacC6wax/PXQQsApg7d25UEmDLDVoD7KDFY2Yjq7UfQESsAS4D3gxMkzScgGYCK+uMxcwsd3V8C2i79MkfSVOBQ4FbKRLB0Wmz+cD5VcdiZmbr1VEFtAOwOLUDbAKcHRE/k3QL8ENJ/wgsA75VQyy1K9eHD8+bmQ2CyhNARNwIvHGE5XcBb6r6+E1zA6iZDSqPBWRmliknADOzTDkBmJllyoPBWaX8C2Jmg8sJwCrlRnCzweUqIDOzTDkBmJllylVANtDa+utqZm3gBGADra2/rmbWBq4CMjPLlBOAmVmmnADMzDLlNoAJpnP00fJyM7MyJ4AJxt+SMbNeuQrIzCxTTgBmZplyFZC1RufAcq7uMhsbJwBrjfIF353CzMZu1CogSXPqCMTMzOrVyx3ARZJeDfwOuBG4afhvRDxRZXBmZladURNAROwmaTNgN2BPYC/gCGAvSc9GxI4Vx2hmZhXoqQ0gIp4Flkm6A1gLTAdeR3EnYGZmLTRqApC0CzAP+CtgO+Bi4HvAgoj4U7XhmZlZVXq5A7gVWAacBpyf7gbMzKzlekkAJwB7AB8BvirpUYqG4JsoGoLPqy48MzOrSi+NwKeX5yXNZH1j8FHAeZVEZmZmleq7I1hErABWABeNfzhmZlYXjwVkZpYpJwAzs0w5AZiZZcqDwdlAOGDhpaxcsxaod6TPfo/bVJxmVXACsIGwcs1a7lk4D6h3pM9+j9tUnGZVcBWQmVmmnADMzDLlKqA+leuAy1wfPBh6qaPv3MYsV04AfSrXAZe5Pngw9FJH3+01NMtN5VVAkmZJukzSLZKWSzoxLd9G0sWSfp/+bl11LGZmtl4dbQDrgE9ExG7A/sBHJO0GnARcEhE7A5ekeTMzq0nlCSAiHoiI36bpJymGl55B8atii9Nmi4Ejq47FzMzWq7UNQNIQ8EbgWmD7iHggrXoQ2L7LcxYACwBmz55dQ5TWtBnTpr5Yf19nI23ncd2obxNdbQlA0hbAj4G/jYgnJL24LiJCUoz0vIhYBCwCmDt37ojb2MTS1IW3fFw36lsOaukHIGlTiov/9yLiJ2nxQ5J2SOt3AB6uIxYzMyvU8S0gAd8Cbo2IL5VWXQDMT9PzgfOrjsXMzNarowroAOB9wE2Srk/LPg0sBM6WdDxwL/CuGmKpTFP11rmqs76+345jdQwY50HpbDxUngAi4ipAXVYfXPXx6+J/wHrVWV/fb8exOgaM86B0Nh48FpCZWaacAMzMMuUEYGaWKQ8GV+KGtYmlLR27/L6zpjgBlLhhbWJpS8cuv++sKa4CMjPLlBOAmVmmXAVkVhPX9dugcQIwq4nr+m3QuArIzCxTTgBmZplyAjAzy5TbAHrQ72iQNrpyJ63hedswvw9tvDkB9KDf0SBtdP4GTP/8PrTx5iogM7NMOQGYmWXKVUCWhSp+sa2zHaO8vIl4OpXbDDqPPVwF10vnNHdgq0/d59oJwLJQxT/SWPZZx0W0W5tBOWn10jnNHdjqU/e5dhWQmVmmnADMzDLlBGBmlim3AZgNsKY6f7nTWR6cAMwGWFOdv9zpLA+uAjIzy5QTgJlZplwFZDaCOjpqjZfOWAepo5Y7kQ02JwCzEbTpQlWOddA6arkT2WBzFZCZWaacAMzMMuUEYGaWKbcBWOu1qcHWDAancdwJwFqvTQ22ZjA4jeOuAjIzy5QTgJlZplwFZFahftsnOn9lzG0aViUnALMK9ds+4fYMq1PlVUCSvi3pYUk3l5ZtI+liSb9Pf7euOg4zM3upOtoAvgO8vWPZScAlEbEzcEmaNzOzGlWeACLiCuCxjsVHAIvT9GLgyKrjMDOzl2qqDWD7iHggTT8IbN9tQ0kLgAUAs2fPriE0s4ltLB3nBnnk0bEalM5ZdWq8ETgiQlJsYP0iYBHA3Llzu25nZr0Zy4VtkEceHatB6ZxVp6b6ATwkaQeA9PfhhuIwM8tWUwngAmB+mp4PnN9QHGZm2aq8CkjSD4ADgemSVgCfBRYCZ0s6HrgXeFfVcfTLA4xZG7X9fdutHr68vHOdbbzKE0BEvKfLqoOrPvZY+M1lbdT29223evjy8s51tvE8FpCZWaacAMzMMuUEYGaWqcb7AZhZ+/XbiaqX7Xtt0O6l4XgsjcadDdDl+Nre5uIEYGZj1m8nql627/Xi2kvD8VgajTsboIdNhIZoVwGZmWXKCcDMLFOuAjKzcTWWX0Gro/Pahur0qzxWv20GdQy85wRgZuNq0H8FrVudftXH6rfNoI6B91wFZGaWKScAM7NMOQGYmWXKbQBm1jrdGo77XT5e2vprYk4AZtY63S6w/S4fL239NTFXAZmZZcoJwMwsU9lXAXXW3ZlZbwa9w1fVxqv8TZ6X7BNAnZ1CzCaSQe/wVbXxKn+T58VVQGZmmXICMDPLlBOAmVmmsmkDqHMEQDNrt7E0zLapsTubBODGXjPr1VgaZtvU2O0qIDOzTDkBmJllKpsqIDPbsDbVXdv4cAIwM6Bdddc2PlwFZGaWKScAM7NMOQGYmWXKCcDMLFNOAGZmmXICMDPLlBOAmVmmnADMzDLlBGBmlqlGE4Ckt0u6TdIdkk5qMhYzs9w0lgAkTQK+BrwD2A14j6TdmorHzCw3Td4BvAm4IyLuiog/AT8EjmgwHjOzrDQ5GNwM4L7S/Apgv86NJC0AFqTZpyTdthHHmg6s0mkb8cx2mw6sajqImrnMeciuzOn6tbHlnjPSwoEfDTQiFgGLxrIPSddFxNxxCqk1ciy3y5yHHMsM41/uJquAVgKzSvMz0zIzM6tBkwngN8DOknaUNAU4BrigwXjMzLLSWBVQRKyT9FHgF8Ak4NsRsbyiw42pCqnFciy3y5yHHMsM41xuRcR47s/MzFrCPYHNzDLlBGBmlqkJlQBGG1pC0maSzkrrr5U01ECY46qHMv+FpN9KWifp6CZirEIP5f64pFsk3SjpEkkjfg+6TXoo84ck3STpeklXTYSe9b0OFyPpKEkhqfVfDe3hdT5O0iPpdb5e0gc2+mARMSEeFA3JdwKvBaYANwC7dWzzYeDrafoY4Kym466hzEPAXsAZwNFNx1xjuf8D8Mo0fUImr/WrStOHAz9vOu6qy5y22xK4ArgGmNt03DW8zscBXx2P402kO4BehpY4Alicps8BDpakGmMcb6OWOSLuiYgbgReaCLAivZT7soh4Js1eQ9HPpM16KfMTpdnNgbZ/w6PX4WI+D5wG/LHO4CpS6xA5EykBjDS0xIxu20TEOuBxYNtaoqtGL2WeiPot9/HARZVGVL2eyizpI5LuBP4Z+FhNsVVl1DJL2geYFREX1hlYhXp9bx+VqjfPkTRrhPU9mUgJwOxlJL0XmAt8selY6hARX4uIPwM+CZzcdDxVkrQJ8CXgE03HUrOfAkMRsRdwMetrNfo2kRJAL0NLvLiNpMnAVsCjtURXjVyH0+ip3JIOAT4DHB4Rz9YUW1X6fa1/CBxZZUA1GK3MWwJ7AEsk3QPsD1zQ8obgUV/niHi09H7+JvDnG3uwiZQAehla4gJgfpo+Grg0UqtKS+U6nMao5Zb0RuB0iov/ww3EON56KfPOpdl5wO9rjK8KGyxzRDweEdMjYigihijaeg6PiOuaCXdc9PI671CaPRy4daOP1nSr9zi3oB8G3E7Riv6ZtOxUijcFwCuAHwF3AL8GXtt0zDWUeV+KesSnKe52ljcdc03l/jfgIeD69Lig6ZhrKPNXgOWpvJcBuzcdc9Vl7th2CS3/FlCPr/P/TK/zDel13nVjj+WhIMzMMjWRqoDMzKwPTgBmZplyAjAzy5QTgJlZppwAzMwy5QRgZpYpJwCbkCQNSVor6fqmYxkvkk6R9N+7rJuahgb+k6Tpdcdm7eQEYBOGpEkdi+6MiL37eL7S+DKtExFrU1nvbzoWa49WvtnNhkn6kaTTJV0DfGqUbc+TtFTSckkL0rKh9OMbZwA3A7Mk/Zc00uINkr4raXNJF6b5myW9u7TP90r6dfr0ffpwEurcR1r28fT8myX9ben4t0r6Rorrl5Kmlvb/GUm3S7oK2CUt6xqPWT8mNx2A2RjtCZwdEfv3sO1/jYjH0gX2N5J+nJbvDMyPiGsk7U4xiua/j4hVkrYB3g7cHxHzACRtlf6+Hng3cEBEPCfp/wDHSlrauQ9Jfw68H9gPEHCtpMuB1en474mIv5F0NnAUcGZ6zjHA3hT/q78FlnaLx6xfvgOw1pL0CmAbinFSevExSTdQDBo2i+LCC3BvRFyTpg8CfhQRqwAi4jHgJuBQSadJemtEPJ62PZhiJMbfpLaGgyl+yWmkfbwFODcino6Ip4CfAG9N+7k7Iq5P00spfsWNtP7ciHgmih97GR4UrFs8Zn1xArA22x24Noof99kgSQcChwBvjog3AMsoBgeEYqC8riLidmAfigvvP0r6h+HdAosjYu/02CUiTtmIcpSHqn6eUe7MNxCPWV+cAKzN9gRu7HHbrYDVEfGMpF0pxo4fyaXAf5a0LUCqvnkN8ExEnEnxwzL7pG0vAY6W9OrStnNG2gdwJXCkpFdK2hz4T2nZhlyRnjNV0pbAO9P+usVj1he3AVib7UkxrHcvfg58SNKtwG0U1UAvExHLJX0BuFzS8xR3Cj8AvijpBeA5ih+ZJyJukXQy8Mv07aHngI+ktoSX7CMijpP0nVK834yIZZKGugUcEb+VdBbFsL8PU4wVP1zul8Vj1i8PB20TUrqw/iwi9mg6ljqp+GWsucPtD2Yb4iogm6ieB7aaSB3BNmS4IxiwKfBCw+FYS/gOwMwsU74DMDPLlBOAmVmmnADMzDLlBGBmliknADOzTDkBmJllygnAzCxTTgBmZpn6/7F+FZm/bXgGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(twomass[\"_r\"],bins=120,histtype=\"step\")\n",
    "#plt.vlines(1,0,100,color=\"k\",linestyle=\"--\")\n",
    "#plt.xlim(-1,2)\n",
    "plt.xlabel(\"$r$ [arcseconds]\")\n",
    "plt.ylabel(\"$N$\")\n",
    "plt.title(\"2MASS/BeSS distance histogram\")\n",
    "plt.show()\n",
    "\n",
    "plt.hist(catwise[\"_r\"],bins=120,histtype=\"step\")\n",
    "#plt.vlines(1,0,1750,color=\"k\",linestyle=\"--\")\n",
    "#plt.xlim(-1,40)\n",
    "plt.xlabel(\"$r$ [arcseconds]\")\n",
    "plt.ylabel(\"$N$\")\n",
    "plt.title(\"CATWISE (2020)/BeSS distance histogram\")\n",
    "plt.show()\n",
    "\n",
    "plt.hist(allwise[\"_r\"],bins=100,histtype=\"step\")\n",
    "#plt.vlines(1,0,1650,color=\"k\",linestyle=\"--\")\n",
    "#plt.xlim(0,1)\n",
    "plt.xlabel(\"$r$ [arcseconds]\")\n",
    "plt.ylabel(\"$N$\")\n",
    "plt.title(\"AllWISE/BeSS distance histogram\")\n",
    "plt.show()\n",
    "\n",
    "plt.hist(unwise[\"_r\"],bins=100,histtype=\"step\")\n",
    "#plt.vlines(1,0,1650,color=\"k\",linestyle=\"--\")\n",
    "#plt.xlim(-1,40)\n",
    "plt.xlabel(\"$r$ [arcseconds]\")\n",
    "plt.ylabel(\"$N$\")\n",
    "plt.title(\"UnWISE/BeSS distance histogram\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filter WISE stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "../dark-GRBs/code/uncertainty.py:249: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = np.log10(self.value)\n",
      "../dark-GRBs/code/uncertainty.py:250: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pos = self.plus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:251: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  neg = self.minus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:62: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  self.sign = self.value/np.abs(self.value) if self.value != 0 else 1\n",
      "../dark-GRBs/code/uncertainty.py:249: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = np.log10(self.value)\n",
      "../dark-GRBs/code/uncertainty.py:250: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pos = self.plus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:251: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  neg = self.minus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:62: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  self.sign = self.value/np.abs(self.value) if self.value != 0 else 1\n",
      "../dark-GRBs/code/uncertainty.py:249: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = np.log10(self.value)\n",
      "../dark-GRBs/code/uncertainty.py:250: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pos = self.plus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:251: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  neg = self.minus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:62: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  self.sign = self.value/np.abs(self.value) if self.value != 0 else 1\n",
      "../dark-GRBs/code/uncertainty.py:249: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = np.log10(self.value)\n",
      "../dark-GRBs/code/uncertainty.py:250: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pos = self.plus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:251: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  neg = self.minus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:62: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  self.sign = self.value/np.abs(self.value) if self.value != 0 else 1\n",
      "../dark-GRBs/code/uncertainty.py:249: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = np.log10(self.value)\n",
      "../dark-GRBs/code/uncertainty.py:250: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pos = self.plus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:251: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  neg = self.minus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:62: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  self.sign = self.value/np.abs(self.value) if self.value != 0 else 1\n",
      "../dark-GRBs/code/uncertainty.py:249: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = np.log10(self.value)\n",
      "../dark-GRBs/code/uncertainty.py:250: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pos = self.plus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:251: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  neg = self.minus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:62: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  self.sign = self.value/np.abs(self.value) if self.value != 0 else 1\n",
      "../dark-GRBs/code/uncertainty.py:249: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = np.log10(self.value)\n",
      "../dark-GRBs/code/uncertainty.py:250: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pos = self.plus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:251: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  neg = self.minus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:62: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  self.sign = self.value/np.abs(self.value) if self.value != 0 else 1\n",
      "../dark-GRBs/code/uncertainty.py:249: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = np.log10(self.value)\n",
      "../dark-GRBs/code/uncertainty.py:250: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pos = self.plus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:251: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  neg = self.minus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:62: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  self.sign = self.value/np.abs(self.value) if self.value != 0 else 1\n",
      "../dark-GRBs/code/uncertainty.py:249: RuntimeWarning: divide by zero encountered in log10\n",
      "  result = np.log10(self.value)\n",
      "../dark-GRBs/code/uncertainty.py:250: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pos = self.plus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:251: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  neg = self.minus/(self.value*np.log(10))\n",
      "../dark-GRBs/code/uncertainty.py:62: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  self.sign = self.value/np.abs(self.value) if self.value != 0 else 1\n"
     ]
    }
   ],
   "source": [
    "wise = pd.DataFrame(columns=[\"_q\",\"W1\",\"e_W1\",\"W2\",\"e_W2\",\"W3\",\"e_W3\",\"W4\",\"e_W4\"])\n",
    "for q in BeSS[\"#\"]:\n",
    "    \n",
    "    W1 = W2 = W3 = W4 = np.nan\n",
    "    e_W1 = e_W2 = e_W3 = e_W4 = np.nan\n",
    "    \n",
    "    if (catwise[\"_q\"] == q).any():\n",
    "        W1 = catwise.loc[catwise[\"_q\"]==q, \"W1mproPM\"].values[0]\n",
    "        e_W1 = catwise.loc[catwise[\"_q\"]==q, \"e_W1mproPM\"].values[0]\n",
    "        W2 = catwise.loc[catwise[\"_q\"]==q, \"W2mproPM\"].values[0]\n",
    "        e_W2 = catwise.loc[catwise[\"_q\"]==q, \"e_W2mproPM\"].values[0]\n",
    "        \n",
    "    elif (unwise[\"_q\"] == q).any():\n",
    "        fw1 = unwise.loc[unwise[\"_q\"]==q, \"FW1\"].values[0]\n",
    "        e_fw1 = unwise.loc[unwise[\"_q\"]==q, \"e_FW1\"].values[0]\n",
    "        FW1 = AsymmetricUncertainty(fw1,e_fw1,e_fw1)\n",
    "        fw2 = unwise.loc[unwise[\"_q\"]==q, \"FW2\"].values[0]\n",
    "        e_fw2 = unwise.loc[unwise[\"_q\"]==q, \"e_FW2\"].values[0]\n",
    "        FW2 = AsymmetricUncertainty(fw2,e_fw2,e_fw2)\n",
    "        w1 = 22.5 - 2.5*np.log10(FW1) - 0.004\n",
    "        w2 = 22.5 - 2.5*np.log10(FW2) - 0.032\n",
    "        W1 = w1.value\n",
    "        e_W1 = w1.plus\n",
    "        W2 = w2.value\n",
    "        e_W2 = w2.plus\n",
    "    \n",
    "    if (allwise[\"_q\"] == q).any():\n",
    "        if pd.isna(W1):\n",
    "            W1 = allwise.loc[allwise[\"_q\"]==q, \"W1mag\"].values[0]\n",
    "            e_W1 = allwise.loc[allwise[\"_q\"]==q, \"e_W1mag\"].values[0]\n",
    "        if pd.isna(W2):\n",
    "            W2 = allwise.loc[allwise[\"_q\"]==q, \"W2mag\"].values[0]\n",
    "            e_W2 = allwise.loc[allwise[\"_q\"]==q, \"e_W2mag\"].values[0]\n",
    "        W3 = allwise.loc[allwise[\"_q\"]==q, \"W3mag\"].values[0]\n",
    "        e_W3 = allwise.loc[allwise[\"_q\"]==q, \"e_W3mag\"].values[0]\n",
    "        W4 = allwise.loc[allwise[\"_q\"]==q, \"W4mag\"].values[0]\n",
    "        e_W4 = allwise.loc[allwise[\"_q\"]==q, \"e_W4mag\"].values[0]\n",
    "        \n",
    "    wise = wise.append(dict(zip(wise.columns, [q,W1,e_W1,W2,e_W2,W3,e_W3,W4,e_W4])), ignore_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XMM, Chandra, and Swift/XRT catalogs\n",
    "\n",
    "| Band      | Energy bin range |\n",
    "|-----------|------------------|\n",
    "| XMM Flux1 | 0.2 - 0.5 keV    |\n",
    "| XMM Flux2 | 0.5 - 1.0 keV    |\n",
    "| XMM Flux3 | 1.0 - 2.0 keV    |\n",
    "| XMM Flux4 | 2.0 - 4.5 keV    |\n",
    "| XMM Flux5 | 4.5 - 12.0 keV   |\n",
    "| XMM Flux8 | 0.2 - 12.0 keV   |\n",
    "| XMM Flux9 | 0.5 - 4.5 keV    |\n",
    "| CXO Fluxu | 0.2 - 0.5 keV    |\n",
    "| CXO Fluxs | 0.5 - 1.2 keV    |\n",
    "| CXO Fluxm | 1.2 - 2.0 keV    |\n",
    "| CXO Fluxh | 2.0 - 7.0 keV    |\n",
    "| CXO Fluxb | 0.5 - 7.0 keV    |\n",
    "| XRT Soft  | 0.3 - 1.0 keV    |\n",
    "| XRT Med   | 1.0 - 2.0 keV    |\n",
    "| XRT Hard  | 2.0 - 10.0 keV   |\n",
    "\n",
    "Load the full catalogs from locally stored copies. Crossmatch from vizier remotely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully matched BeSS with IX/59/xmm4dr9s\n",
      "Successfully matched BeSS with IX/57/csc2master\n",
      "Successfully matched BeSS with IX/58/2sxpscle\n"
     ]
    }
   ],
   "source": [
    "cat_ids = [\"IX/59/xmm4dr9s\",\"IX/57/csc2master\",\"IX/58/2sxpscle\"]\n",
    "Xmatch = [None, None, None]\n",
    "radii = [\"2s\",\"1s\",\"5s\"]\n",
    "for i in range(len(cat_ids)):\n",
    "    Xmatch[i] = v.query_region(BeSS, catalog=cat_ids[i], radius=radii[i])\n",
    "    print(\"Successfully matched BeSS with\",cat_ids[i])\n",
    "[XMM_xmatch, CSC_xmatch, XRT_xmatch] = Xmatch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "Chandra_sources = CSC_xmatch[0][\"_q\",\"_2CXO\",'Fluxb', 'Fluxh', 'Fluxm',\n",
    "       'Fluxs', 'Fluxu', 'Fluxw', 'HRhm', 'HRhs', 'HRms', 'GamPL', 'NHPL']\n",
    "chandra_bands = {\"b\":(0.5,7), \"h\":(2,7), \"m\":(1.2,2), \"s\":(0.5,1.2), \"u\":(0.2,0.5)}\n",
    "\n",
    "XMM_sources = XMM_xmatch[0][\"_q\",\"_4XMM\",'Flux1', 'e_Flux1',\n",
    "       'Flux2', 'e_Flux2', 'Flux3', 'e_Flux3', 'Flux4', 'e_Flux4', 'Flux5',\n",
    "       'e_Flux5', 'Flux8', 'e_Flux8', 'Flux9', 'e_Flux9', 'HR1', 'e_HR1',\n",
    "       'HR2', 'e_HR2', 'HR3', 'e_HR3', 'HR4', 'e_HR4']\n",
    "xmm_bands = dict(zip(list(\"1234589\"),[(0.2,0.5),(0.5,1.0),(1.0,2.0),(2.0,4.5),(4.5,12.0),(0.2,12.0),(0.5,4.5)]))\n",
    "\n",
    "XRT_sources = XRT_xmatch[0][\"_q\",'_2SXPS', 'CR0', 'E_CR0', 'e_CR0', 'HR1', 'HR2', 'CR1',\n",
    "       'e_CR1', 'CR2', 'E_CR2', 'e_CR2', 'CR3', 'E_CR3', 'e_CR3', 'PcstS0',\n",
    "       'FPO0', 'FPU0', 'Gam1', 'E_Gam1', 'e_Gam1']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compile a master catalog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "objs = pd.merge(wise, Twomass_xmatch[0][['_q',\"Jmag\",\"e_Jmag\",\"Hmag\",\"e_Hmag\",\"Kmag\",\"e_Kmag\"]].to_pandas(), on=\"_q\", how=\"outer\")\n",
    "objs = pd.merge(gaia, objs, on=\"_q\", how=\"outer\")\n",
    "xray = pd.merge(XMM_sources.to_pandas(), Chandra_sources.to_pandas(), on=\"_q\", how=\"outer\")\n",
    "xray = pd.merge(xray, XRT_sources.to_pandas(), on=\"_q\", how=\"outer\")\n",
    "\n",
    "master = pd.merge(BeSS.to_pandas(), pd.merge(objs, xray, on=\"_q\", how=\"outer\"), left_on=\"#\", right_on=\"_q\", how=\"left\")\n",
    "master.rename(columns={\"HR1_x\":\"HR1\",\"HR2_x\":\"HR2\",\"HR1_y\":\"SwiftHR1\",\"HR2_y\":\"SwiftHR2\"},inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve for individual-band Swift fluxes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bandshift(fluxA,Gamma,bandA,bandB):\n",
    "    a1,a2 = bandA\n",
    "    b1,b2 = bandB\n",
    "    if Gamma == 2:\n",
    "        integral = lambda lower,upper : np.log(upper) - np.log(lower)\n",
    "    else:\n",
    "        integral = lambda lower,upper : (upper**(2-Gamma) - lower**(2-Gamma))/(2-Gamma)\n",
    "    scalar = fluxA/integral(a1,a2)\n",
    "    fluxB = scalar*integral(b1,b2)\n",
    "    return fluxB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# total flux and power-law spectrum-based\n",
    "for i in master[master[\"_2SXPS\"].notna()].index:\n",
    "    total = master.loc[i,\"FPO0\"]\n",
    "    gamma = master.loc[i,\"Gam1\"]\n",
    "    if pd.isna(gamma):\n",
    "        gamma = 2\n",
    "\n",
    "    master.loc[i,\"XRTsoft\"] = bandshift(total,gamma,(0.3,10),(0.3,1))\n",
    "    master.loc[i,\"XRTmed\"] = bandshift(total,gamma,(0.3,10),(1,2))\n",
    "    master.loc[i,\"XRThard\"] = bandshift(total,gamma,(0.3,10),(2,10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def infer_gamma(HR,freq1a,freq1b,freq2a,freq2b):\n",
    "    #G = sympy.var(\"Gamma\")\n",
    "    #hr_eq = sympy.Eq(((freq2b**(2-G) - freq2a**(2-G)) - (freq1b**(2-G) - freq1a**(2-G)))/((freq2b**(2-G) - freq2a**(2-G)) + (freq1b**(2-G) - freq1a**(2-G))), HR)\n",
    "    hr_func = lambda G : -HR + ((freq2b**(2-G) - freq2a**(2-G)) - (freq1b**(2-G) - freq1a**(2-G)))/((freq2b**(2-G) - freq2a**(2-G)) + (freq1b**(2-G) - freq1a**(2-G)))\n",
    "    solution = optimize.fsolve(hr_func, 1) #sympy.solve(hr_eq, G)\n",
    "    return solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "# symbolic method using hardness ratios (deprecated)\n",
    "# for i in master[master[\"_2SXPS\"].notna()].index:\n",
    "#     # knowns\n",
    "#     SwiftHR1 = master.loc[i,\"SwiftHR1\"]\n",
    "#     SwiftHR2 = master.loc[i,\"SwiftHR2\"]\n",
    "#     SwiftFlux = master.loc[i,\"FPO0\"]\n",
    "\n",
    "#     # unknowns\n",
    "#     F0,F1,F2 = sympy.symbols(\"F_s, F_m, F_h\")\n",
    "    \n",
    "#     # equations\n",
    "#     fluxsum = sympy.Eq(F0+F1+F2, SwiftFlux)\n",
    "#     hr1 = sympy.Eq((F1-F0)/(F1+F0), SwiftHR1)\n",
    "#     hr2 = sympy.Eq((F2-F1)/(F2+F1), SwiftHR2)\n",
    "    \n",
    "#     # solution\n",
    "#     solution = sympy.solve((fluxsum,hr1,hr2), (F0,F1,F2))\n",
    "#     master.loc[i,\"XRTsoft\"] = float(solution[F0])\n",
    "#     master.loc[i,\"XRTmed\"] = float(solution[F1])\n",
    "#     master.loc[i,\"XRThard\"] = float(solution[F2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flux band conversions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "XMM            34\n",
       "CXO            28\n",
       "XRT            23\n",
       "CXO+XMM+XRT    13\n",
       "CXO+XMM        13\n",
       "XMM+XRT        10\n",
       "CXO+XRT         7\n",
       "Name: Flag, dtype: int64"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CSC_fluxs = ['Fluxs', 'Fluxm', 'Fluxh']\n",
    "XMM_fluxs = ['Flux1', 'Flux2', 'Flux3', 'Flux4', 'Flux5']\n",
    "XRT_fluxs = ['XRTsoft', 'XRTmed', 'XRThard']\n",
    "\n",
    "xray = master[['#','_4XMM', 'Flux1', 'Flux2', 'Flux3', 'Flux4', 'Flux5', 'HR1', 'HR2',\n",
    "               'HR3', 'HR4', '_2CXO', 'Fluxb', 'Fluxh', 'Fluxm', 'Fluxs', 'Fluxu',\n",
    "               'Fluxw', 'GamPL', '_2SXPS', 'SwiftHR1', 'XRTsoft', 'XRTmed', 'XRThard',\n",
    "               'SwiftHR2', 'CR1', 'CR2', 'CR3', 'FPO0', 'Gam1']].copy()\n",
    "\n",
    "for i,row in xray.iterrows():\n",
    "    flag = []\n",
    "    if pd.notna(row[\"_2CXO\"]):\n",
    "        flag.append(\"CXO\")\n",
    "    if pd.notna(row[\"_4XMM\"]):\n",
    "        flag.append(\"XMM\")\n",
    "    if pd.notna(row[\"_2SXPS\"]):\n",
    "        flag.append(\"XRT\")\n",
    "    flag = \"+\".join(flag)\n",
    "    if len(flag)<3:\n",
    "        flag = np.nan\n",
    "    xray.loc[i,\"Flag\"] = flag\n",
    "\n",
    "xray.dropna(subset=[\"Flag\"],inplace=True)\n",
    "    \n",
    "xray[\"Flux12\"] = xray.Flux1 + xray.Flux2\n",
    "xray[\"Flux45\"] = xray.Flux4 + xray.Flux5\n",
    "\n",
    "xray = xray.replace(0,1e-17)\n",
    "\n",
    "xray[\"Flag\"].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.multioutput import MultiOutputRegressor\n",
    "from sklearn.ensemble import GradientBoostingRegressor, RandomForestRegressor\n",
    "\n",
    "# convert 2SXPS to CSC bands\n",
    "xrt_to_csc = xray[xray[\"Fluxm\"].notna() & xray[\"Gam1\"].notna()].copy() # sources in common\n",
    "xrt_to_csc.drop([1791,2124],axis=0,inplace=True)\n",
    "\n",
    "X_train,y_train = np.log10(xrt_to_csc[XRT_fluxs].values), np.log10(xrt_to_csc[CSC_fluxs].values)\n",
    "\n",
    "model = MultiOutputRegressor(RandomForestRegressor())\n",
    "model.fit(X_train,y_train)\n",
    "y_pred = model.predict(np.log10(xray.loc[xray[\"Gam1\"].notna(),XRT_fluxs].values))\n",
    "xray.loc[xray[\"Gam1\"].notna(), [\"from_XRTs\",\"from_XRTm\",\"from_XRTh\"]] = 10**y_pred\n",
    "\n",
    "# convert 4XMM to CSC bands\n",
    "xmm_to_csc = xray[xray[\"Fluxm\"].notna() & xray[\"_4XMM\"].notna()].copy() # sources in common\n",
    "xmm_to_csc.drop([1894,1897,1931,2124],axis=0,inplace=True)\n",
    "\n",
    "X_train,y_train = np.log10(xmm_to_csc[[\"Flux12\",\"Flux3\",\"Flux45\"]].values), np.log10(xmm_to_csc[CSC_fluxs].values)\n",
    "\n",
    "model = MultiOutputRegressor(RandomForestRegressor())\n",
    "model.fit(X_train,y_train)\n",
    "y_pred = model.predict(np.log10(xray.loc[xray[\"_4XMM\"].notna(),[\"Flux12\",\"Flux3\",\"Flux45\"]].values))\n",
    "xray.loc[xray[\"_4XMM\"].notna(), [\"from_XMMs\",\"from_XMMm\",\"from_XMMh\"]] = 10**y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in master[master[\"_2CXO\"].notna() | master[\"_4XMM\"].notna() | master[\"_2SXPS\"].notna()].index:\n",
    "    cxoflux, xmmflux, xrtflux = master.loc[i,CSC_fluxs].values, master.loc[i,XMM_fluxs].values, master.loc[i,XRT_fluxs].values\n",
    "    exists = False\n",
    "    if pd.notna(master.loc[i,\"GamPL\"]): # spectral index from Chandra\n",
    "        Gamma = float(master.loc[i,\"GamPL\"])\n",
    "        exists = True\n",
    "    elif pd.notna(master.loc[i,\"Gam1\"]): # spectral index from Swift\n",
    "        Gamma = float(master.loc[i,\"Gam1\"])\n",
    "        exists = True\n",
    "    else: # assume flat spectrum for general case\n",
    "        Gamma = 1.7\n",
    "    \n",
    "    if all(pd.notna(cxoflux)): # good Chandra observations present\n",
    "        converted = cxoflux\n",
    "        source = \"CXO\"\n",
    "    elif all(pd.notna(xmmflux)): # no Chandra, but good XMM\n",
    "        mergedxmm = np.array([xmmflux[0]+xmmflux[1], xmmflux[2], xmmflux[3]+xmmflux[4]])\n",
    "        origbands = [(0.2,1),(1,2),(2,12)]\n",
    "        converted = [bandshift(flux,Gamma,xmmband,cxoband) for flux,xmmband,cxoband in\n",
    "                     zip(mergedxmm,origbands,[(0.5, 1.2), (1.2, 2), (2, 7)])]\n",
    "        source = \"XMM\"\n",
    "    elif all(pd.notna(xrtflux)): # neither Chandra nor XMM, but good Swift-XRT\n",
    "        if Gamma != 2:\n",
    "            coefs = 1/np.array([(1**(2.-Gamma)-0.3**(2.-Gamma))/(1.2**(2.-Gamma)-0.5**(2.-Gamma)),\n",
    "                          (2**(2.-Gamma)-1**(2.-Gamma))/(2.**(2.-Gamma)-1.2**(2.-Gamma)),\n",
    "                          (10**(2.-Gamma)-2**(2.-Gamma))/(7**(2.-Gamma)-2**(2.-Gamma))])\n",
    "        else:\n",
    "            coefs = 1/np.array([(np.log(1)-np.log(0.3))/(np.log(1.2)-np.log(0.5)),\n",
    "                              (np.log(2)-np.log(1))/(np.log(2.0)-np.log(1.2)),\n",
    "                              (np.log(10)-np.log(2))/(np.log(7)-np.log(2))])\n",
    "        converted = xrtflux * coefs\n",
    "        source = \"XRT\"\n",
    "    else:\n",
    "        continue\n",
    "    \n",
    "    master.loc[i,\"Flux (s band)\"] = converted[0]\n",
    "    master.loc[i,\"Flux (m band)\"] = converted[1]\n",
    "    master.loc[i,\"Flux (h band)\"] = converted[2]\n",
    "    master.loc[i,\"X-ray source\"] = source\n",
    "    if exists:\n",
    "        master.loc[i,\"Gamma\"] = Gamma\n",
    "        \n",
    "xray = pd.merge(xray,master[['#',\"Flux (s band)\",\"Flux (m band)\",\"Flux (h band)\",\"X-ray source\",\"Gamma\"]],on=\"#\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No handles with labels found to put in legend.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x238fda662b0>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD8CAYAAACPWyg8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABKvElEQVR4nO3dd3xV9f348dfnjuRmbwhZJGGEvUQZMiMgCCi27taKRdHaXdtv68+21tZWa62rrYMqWvcWBRco083ee2WTPci643x+f5wkJOAgkOTem/t+Ph55kHvOvee+7yU373zW+6O01gghhBDtYfF2AEIIIfyPJA8hhBDtJslDCCFEu0nyEEII0W6SPIQQQrSbJA8hhBDtZvN2AO0RHx+v09PTvR2GEEL4lY0bN5ZqrRM68pp+lTzS09PZsGGDt8MQQgi/opQ62tHX9ItuK6XUXKXUoqqqKm+HIoQQAj9JHlrrpVrrhVFRUd4ORQghBH6SPIQQQvgWvxrzEEKIQOdyucjLy6OhoeGUcw6Hg5SUFOx2e6fHIclDCCF8XG1VI8uf2MGMG4ZQXF5AREQE6enpKKVa7qO1pqysjLy8PDIyMjo9Jr/otpIBcyFEINvwzmEKDlSx4d0jNDQ0EBcX1yZxACiliIuL+8oWSWfwi+QhA+ZCiEBVW9XI7s+KQMPuTwsxDH1K4mj2dcc7g18kDyGECFTLn9yJx2UAoA2Ns97t5YhMkjyEEMJHFR+tpmBfZcttw6NxNXrwuA3vBdVEkocQQviYxjoXALs/KcBiPakrSkNtZeNXPq4rd4aV5CGEED4kf28Fz9/xObs/LaDoUDWGp21CqC1zU15edkqiaJ5t5XA4uiROv5iqq5SaC8zt27evt0MRQohOobVm60e5fPrGQaISQuiZEcWVv0865X7N6zz27NlzyrnmdR5dQXVlM+dsjR49WkthRCFEd+NscLPquT0c2FBM5ogELrhuIEEhHfe3vVJqo9Z6dIddED9peQghRHdWsL+SgxuLGXdpH0bOSOvSKbdnSpKHEEJ4SU15AxGxDtKHxnPNnWOJ7hHq7ZBOmwyYCyFEFzMMzRdvH+K5P35G8dFqAL9KHCAtDyGE6FINtS5WLN5Jzs5yBozvRWxSmLdDOiOSPIQQoouU5Nbw/uPbOV7RyORrshg8Mckvxje+il8kD5mqK4ToDg5vLcXj1lz661EkZvh3rT6ZqiuEEJ3I4zaoLq0nJjEMbWga69w4wjt/v43WOmOqrgyYCyFEJ6mtbGTJ/Zt585+bcNa7URbV5Ymjs/hFt5UQQvibgv2VfPDfHTgbPWRfO6BDF/35gu71aoQQwsu01mxfnccnrx4gIt7BxT8fQVxyuLfD6nCSPIQQooPl7akgbUgc064fRHA3a3E0656vSgghulhVSR3KooiMC2H6gsHYbBaUxT+n4Z4OGTAXQoizdGR7Ka/evYFVz5qVbu1B1m6dOEBaHkIIcca0oVn/7hHWv3OY+JRwpnxvgLdD6jKSPIQQ4gw01rv5cPFOjmwvI2tMIpO/l4U9yOrtsLqMJA8hhDgDFouitsrJpKv6M2Ryst+WGTlTfpE8pDyJEMJXHNpSQurAWOzBVi777TlYrIE5dOwXr1prvVRrvTAqyr9rwQgh/JfHY7DulX2899h2tnyYAxCwiQP8pOUhhBDeVFvVyAf/3UHhgSqGZacwamZvb4fkdZI8hBDiGxQfrebdR7bRWOdm+g8H0f+8RG+H5BMkeQghxDdwhNkJiw5mzk9HEJ/S/cqMnKnA7bATQoiv4XZ62L46D601kfEhXPa70ZI4TiItDyGEaKW6tJ73Ht9Oad5x4lLCSeobHXDTcE+HJA8hhGiSs7OM5Yt3gobZtwwjqW+0t0PyWZI8hBAC2LYql3Wv7CcuKYxZNw8lKiHU2yH5NEkeQggBxCWFk3VeIpOvycIeHDhlRs6UJA8hRMAqL6glf18FQ6ekkJwVQ3JWjLdD8huSPIQQAenAxmI+emY3QQ4r/cckdttNmzqLvFtCiIBieAw+W3KILStySMyMZObCoZI4zoC8Y0KIgKG15p3/bCNnVzlDJydz/uX9sNpkuduZkOQhhAgYSikyRybQ/7yeZI3t5e1w/FqXJg+lVCZwOxCltb6s6dg8YDYQCTyptV7elTEJIbq/nevycYTb6TOyB4MnJns7nG7htNtrSqnFSqlipdSOk47PVErtVUodUEr97puuobU+pLVecNKxJVrrG4GbgSvbE7wQQnwTt8vDymd3s/r5vexfX+ztcLqV9rQ8ngb+DTzTfEApZQX+A0wH8oD1Sqm3AStw90mP/6HW+pv+937fdC0hhDhr1WX1fLBoB8VHaxh9UTrnzsnwdkjdymknD631WqVU+kmHzwMOaK0PASilXgIu0VrfDcw5nesqs2jMPcB7WutNX3F+IbAQIC0t7XTDFUIEsNqqRl792wYMj8FFPxpKxvAEb4fU7ZztNINkILfV7bymY19JKRWnlHoMGKmUuq3p8E+BacBlSqmbT36M1nqR1nq01np0QoL8AAghvl1YVDDDp6Vy+W3nSuLoJF06YK61LsMc22h97GHg4a6MQwjR/Tgb3Kx5YS8jZ6QRnxLB6Fnp3g6pWzvblkc+kNrqdkrTsQ6llJqrlFpUVVXV0ZcWQnQDFUW1vHbPBvZvKKb4aI23wwkIZ5s81gP9lFIZSqkg4Crg7bMPqy2t9VKt9cKoqKiOvrQQws8d3FzMq/dsoKHWxSU/H8Gg85O8HVJAaM9U3ReBz4AspVSeUmqB1toN/AT4ANgNvKK13tk5oQohRFtHtpXy/uM7iO0VxhX/71wpbNiFlNba2zF8K6XUXGBu3759b9y/f7+3wxFCeJnWGqUUHo/BtpV5DJuSgtUuZUa+jlJqo9Z6dEde0y/ebem2EkI0Kz5azZv/3ER9jROr1cLI6WmSOLxAalsJIfzGrk8KWPviPkIi7dRVOwmJCPJ2SAHLL5JHq24rb4cihPACj8tg7Sv72LWugJQBMcy4YTAh4ZI4vMkv2nrSbSVEYPvszYPsWlfAqAt7M/dnIyRx+AC/aHkIIQKTNjTKojhnVm+Ss6JltbgP8YuWhxAisGit2bwih7ce2ozHYxASESSJw8f4RfKQFeZCBA5ng5vlT+zk09cP4Ai1Y7h9fzlBIPKL5CFjHkIEhoqiWl77+0YObipm3KV9uHDhEOzBVm+HJb6CjHkIIXyC1prlT+6kvsbJ3J+PIHVArLdDEt9AkocQwqsMQ6MNjdVmYdr1gwhy2IiIdXg7LPEtJHkIIbym4biLFYt3EhHnYMr3BhCXFO7tkMRp8osxDxkwF6L7Kcmp4ZW715O3r4KEtAhvhyPayS+ShwyYC9G97Pm8kNf/sRFtaL5z6zkMnvi1G5AKHyXdVkKILlVX7WTti/tIzIxkxoIhhEbKanF/JMlDCNElGmpdBIfaCI0M4tJfjyIuKQyL1S86P8RXkP85IUSnK9hfwQt3fsHOdQUAJKRGSOLwc9LyEEJ0Gq0121bm8cnrB4hKCKFXXxm37C78InlISXYh/I+r0cOq5/awf/0xMobHc8H8QQSH+MWvHHEa/KLdKLOthPA/RYeqOLCxmDGXZDLrpqGSOLoZv0geHaGkroT578+ntL7U26EI0a1Vl9UDkDowlu/dOZbRs9JRFuXlqERHC5jk8ciWR9h0bBOPbX3M26EI0S1pQ/Pl0kM8/8fPKTpkLuiNSgjxclSiswREO7KkroQlB5ag0by671UGxw7mwowLCbWHejs0IbqFhloXHz61i6M7yhgwNpH4FCkz0t0FRMvjsW0nWhuGNvjjZ39k0suT+OnKn+IyXF6MTAj/V5p3nFfvXk/u7nImX92f7OsGYguSMurdXbdveZTUlfDWgbdwa3fLMbvFzpzMOVQ1VmG32AF4aNNDRAdHMzV1KmmRad4KVwi/k7OzDI/L4NJbR5GYKZNaAkW3Tx6PbXsMQxunHLdZbDww9QHAbI18WfQl20q2cd+G++gb3ZepqVOZlTGLfjH9ujpkIXyex2NQeayOuKRwRs5IY9D5STjC7d4OS3Qhv0geZ7POY2vx1lO6plyGiy3FW1puW5SF5y96nvzj+azKWcWq3FUs3rEYu8VOv5h+NLgb2HRsE+cmnovdKh8QEdhqqxr54L87KC+s5ft/HocjzC6JIwAprf1nf+DRo0frDRs2dMlzVTZUotHEOGJYlbOKn636GeH2cCamTCQ7NZsJyRMID5JBQRFYCg9W8f6i7Tjr3WRfO5B+5/b0dkjiNCilNmqtR3foNSV5fLsGdwOfF37OypyVrMlbQ3lDOTaLjbcueYu0yDQMbWBRATH3QAQorTXbV+fzyav7iYhzMOvmocQlyx9P/qIzkodfdFt5m8PmYErqFKakTsFjeNhaspXPCj8jNSIVgLs+v4u95XuZmjaV7NRsMqIyUEoWRYnupehgJWmDY5l2/SCCQ6WbKtBJy6MDvLjnRd468BY7y3YCkB6ZzqX9LuWHQ37o5ciEODtVJfVorYnuEYrb5cFqtchqcT8kLQ8fdfWAq7l6wNUU1RaxOnc1K3NWcqz2GGDO5Lpvw32M7TWWMb3GEGwN9m6wQpymozvKWLF4J7FJYVx66yhsdlm7IU6Qlkcn0VqjlCK3JpfLl15OrauWEFsIE5InMDV1KpNTJxMZFOntMIU4hTY0G98/whdLDxOXFM6sm4cQlSDVGPyZtDz8SPOYR2pEKmuvXMuXRV+2TANecXQF/87+N5NTJ1NaX4rbcJMYlujliIUAZ72bFU/t4si2UvqP6cmU7w3ALqvFxVeQ5NEFgqxBTEiewITkCdw+9nZ2lO4gKzYLgJf2vMTj2x5nYOxAstOyyU7Lpl90PxlwF15hsSrqa5xMvLI/Q6cky8+h+FrSbeVludW5fJjzIStzVrK1ZCsaTd/ovrx+8esy/Vd0mUNbSkjJiiEoxIZhaCwyKN6tBGy3VXfeSTA1MpXrh1zP9UOup7S+lDW5ayhrKGtJHDetuIn4kHiy07IZnzSeEJuUuBYdx/AYfPrmQbZ+mMvoi9IZc3GmJA5xWqTl4cPchps/fPIH1uStocZZg8PqYFzSOK4ZeA1je431dnjCz9VVO/ngvzso2F/J0KkpnP/dvlht0trtjgK25RGobBYbd0+8G5fhYuOxjazKWcXK3JXk1+RDL7Ni8LuH3yU7NZvUyFRvhyv8SElODe88so3GWhfTrh9E1hiZsCHaR1oefkZrjUd7sFlsvHPoHX637ncA9I3uaw64p2YzKG6QDHSKb3S8ooH3F+1g8jVZJKRGeDsc0cmktpUkj1Pk1eSxKncVK3NWsql4E4Y2+Ojyj+gR2oOSuhKiHdEte5aIwOZ2eti5roBhU1NQFtWyFkl0f9JtJU6REpHCtYOu5dpB11LRUMGW4i30CO0BwB2f3sGW4i1mJeA0sxJwmD3MyxELb6guq+f9x3dQklNDXHIYKQNiJXGIsyLJoxuJccQwNW1qy+2rBlxFXEgcq3NX8+7hd7Fb7FyZdSW/Pe+33gtSdLncXeUsf3InhqG56JZhpAyI9XZIohuQ5NGNTUqZxKSUSbgNN1uKt7AqdxW9I3sDUOeq48cf/ZgJyRPITjMrAYvuZ8eaPNa+tI+YXmHMumko0T2lzIjoGJI8AoDNYmN04mhGJ57o8iyuK6beXc+Dmx7kwU0Pkh6ZTnZaNldmXUlSeJIXoxUdKT41gv7nJTL5mizswVJmRHQcGTAPcEW1RS0D7huKNvDC7BcYGDeQveV7Ka4rZkyvMQRZg7wdpmiH8oJacnaVMWJamrdDET5CBsxFh0sMS2wpKV/VWNVS6fflvS/z6r5XCbWFtnRtTUyZKJWAfdyBjcWsfGY3tmArA8b1whEmM+1E55CWh/hKTo+TLwq/YGXuSlbnrqa0vpTk8GTe+857KKWoddXKzC0fYngMPl9yiM0rcuiZEcnMhUMJj5G9Y4RJ1nlI8vAKQxtsL91OWX0Z2WnZGNpg2qvTSAhNIDs1m6lpU6USsBdprXnnkW0c3V7GkEnJTLi8H1a7lBkRJ0jykOThExo9jTy36zlW5a5ia8lWAFLCU/jlOb9kRvoML0cXmPZ8VojWmoHjZbKDOJWMeQifEGwNZsHQBSwYuoCSuhJW561mVc4qwoPCAdhTvofndz/P1NSpjEsaJ5WAO8mujwuwBVvof24iA8b18nY4IsB0WfJQSmUCtwNRWuvLmo4NBH4OxAMfaa0f7ap4RMdICE3g8v6Xc3n/y1uOHa0+ykdHP2LJgSUtlYCz07KZmT4Th83hxWi7B4/LYO3L+9j1cQEZw+PpN7qndBmKLnda3VZKqcXAHKBYaz2k1fGZwEOAFXhCa33PaVzrtebk0eqYBXhGa/39b3qsdFv5D5fHxYZjG1iZs5JVuauodlaz9sq1OGwOtpZsJc4RR0pEirfD9Ds15Q28//h2io/WMGpmb9l/Q5wWr415KKUmAccxf8EPaTpmBfYB04E8YD1wNWYiufukS/xQa13c9Lg2yUMpdTHwI+BZrfUL3xSHJA//pLUm73geqRFm2fjL3r6MvRV76R/Tn+y0bKamTmVg7ED56/lb1Nc4eeHOL/C4DaZdN4jMkQneDkn4Ca8OmCul0oFlrZLHOOBPWusLm27fBqC1PjlxnHydU1oeTcff0VrP/orjC4GFAGlpaeccPXr0tOIVviu3OpeVuStZmbOSLSVbMLTBvL7z+Mv5fwHAY3iwWmQ19FfZvDyH9GFxxCTKNGlx+nxtwDwZyG11Ow8Y83V3VkrFAX8FRiqlbtNa362UmgJ8BwgG3v2qx2mtFwGLwGx5nEW8wkekRqZy3eDruG7wdZQ3lLMmdw29ws0B36LaIr779neZlDKJ7LRszk86n1B74NZjcja4WfPCXoZlp9IzPZKRM2TVuPANXTZgrrUuA24+6dhqYHVXxSB8T6wjlkv7Xdpy22W4mJo6lTV5a1h2aBlBliDGJo3l1tG3khmV6cVIu15FUS3vPbadymN1JPWLpme6rO4XvuNskkc+0Hrv05SmYx1OKTUXmNu3b9/OuLzwIakRqdw14S7chpvNxZtZmWOucG8ui7I6dzWHqw4zNXUq6VHp3gy1Ux3aXMKH/9uF1Wbh4p+PkDLqwueczZiHDXPA/ALMpLEeuEZrvbNzQpUB80DVese7u7+4mxf2mPMqMqMymZo6ley0bIYlDPNmiB0qZ2cZS/+1lR69I5h501AiYmV6szg73pxt9SIwBXM9xjHgDq31k0qpi4AHMWdYLdZa/7UjgzuZJA8BUHi80KwEnGtWAh4YO5AX57wIwM6ynfSL7ueXlYCbk6ThMdi2Ko8hk5Ox2WXigDh7AVuepFW31Y379+/3djjCh1Q1VlFSV0LfmL7UueqY+NJE7FY7E5MnMjV1KhNTJhIRFOHtML9VSU4Na1/ax8ybhhAWJQUNRccK2OTRTFoe4pu4DBefFXzWsjCxvKEcm8XGnePv5OI+F3s7vK+1+9NC1rywl5AIO7N/PJz4lHBvhyS6GV+bqiuET7Fb7C1b7/7B+APbS7ezMnclg+MGA+Zg+2NbHyM7LZvs1Gz6RPfx6sJEj9vg41f2s2NtPslZMVx4w2BCIvyvu00EJkkeoluyWqyM6DGCET1GtByzKAtWZeVfm//Fvzb/i9SIVLJTs/nxyB97pXjjl0sPsWNtPiNnpDH2kkwsVimjLvyHX3RbyZiH6EjFdcWszl3NytyV5FbnsuzSZSilePvg20QGRTK219hOLeBoGBqLRdFQ66LwQCUZw6XMiOhcMuYhYx6ig7kNNzaLDa01Fy+5mCPVRwixhTA+aTzZadlMSp5EtCO6Q55La83Wj3I5tLmES34xUjZsEl2mM5KH/PSKgGazmD23SineuPgNHp/2OBf3uZjtpdu5/ePb+dfmfwHmbooFxwvO+HlcjR6WP7mTT147QEhEEB6P0SHxC+Et0vIQ4isY2mBX2S7C7GFkRGWwvWQ717x7DVkxWeaAe1o2WTFZ3zjgXlvVyPIndjD2kj6sfmEvFYW1jLkkk1EX9pYKwqJLBWy3lYx5CG8rrS/lnUPvsDJnJZuLN6PRJIUl8ci0R+gT3ecrH7PmhT3sWFdASLgdw9BcuGAIqYOkzIjoegGbPJpJy0P4grL6MtbmrWVt3lrumXQPwdZg/rfzf+wt30t2Wjbjk8ZTur+edx7ZhuHWWGyK79w6ip4ZUd4OXQQoWechhA+IC4nj0n6XtqkGXOeqY23eWjZv2sfKwgMkV/ZHY6Cw4DE87Pm8SJKH6FYkeQjRAaa5v0PMwVGU5dViOJwYyoNFm3WplGFh1ycFjL4oXUqPiG5DZlsJcYbcTk/L9zk7yjA8kP2DAQwZfeqAuNPj4t9PvMShqkNdHaYQncIvWh6yn4fwJbVVjWxbmcfOdflc/PMR9OgdycSr+mMPsqIsik0fHkEZbf8us2kbjQUW5i2Zx4z0GSwctpD+Mf299AqEOHt+kTy01kuBpaNHj77R27GIwFVeWMuWFTns/bII7dFkjkzAFmR2TQU5TnyUDs34iDf3v4nLcLUcs1vsXJRxEQtCF/DinhdZm7eWjy7/yC8q/grxVfwieQjhbW6Xhzf+sRG3y2DQ+UmMmJZKVMJX762+tXhrm8QBZsXfPeV7uGvCXcwfPJ9tJduICIpAa839G+/ngrQL2tThEsLXBc5U3ZoieO16uOxpiOjZoXGJ7scwNIe3lHB4WykXXDcQpRS5u8qJTw3v0Mq3RbVFXLH0CioaKxiTOIaFwxZybuK5sohQdCgpT3I21twLOZ/Dmr97OxLhw9xODzvW5PHCHZ/z/qIdFB6opLayEYDUQbEdXjI9MSyR97/7Pr8e/WsOVh1kwfIFXPf+deTV5HXo8wjR0fyi5XHWK8xriuDBoeBxgjUIrn8fkkaCJXByp/h2ZfnHeevBzdTXuOjRO4KRM3qTOTIBi6VrWgEN7gbePPAmbx14i6dmPkWILYSi2iJ6hPbAouRnVZw5WWF+pt1Wy34FG58GfWJqJbYQiO8L338DwntA2UFwN0BsH7B3Xjlu4VuqS+upKq0ndUAsHrfB6uf2MGB8L5L6RXut66h5L3O34Wbum3MJtYdy47AbmZ42HatF9jQX7ScrzM9ETRFseb5t4rDYYfhVUF0AIU21hj5/BNY/AcoC0b0hvj8kZMG0O80WitsJNtnlrbsoPlrN5uU5HNxUTEScg+//eRxWm4UL5g/ydmgtSUuhuGXELSzatojfrPkNGVEZ3Dj0RmZlzGqpBiyEt3T/n8A194I+qfy1UmaS+N4rJ46NvQXSxkHpfijda/5bfhBm/MU8/+p8yPvSTCrNX4lDIWNil70UcfaKDlXx+ZKD5O+rJMhhZcT0NIZNTUV1UddUe1gtVub2mctFGRexImcFi7Yt4v99/P8Is4eRnZbt7fBEgOv+ySPvS3OsozWP0zzeWlwf86u11l16A+dAaKyZVHYtgfoKSDkPblhhnn/9RkC3TS5xfcAm5Si8zeM28LgMgkJsOOvdVJXUM/67fRk8IYmgEN//CFgtVmamz2RG7xl8WvAp45PGA/DynpfRaC7tdynBVvk5E10rMMY8OprWUFcGDVUnEs7L10LBFqjKOXG/wZfC5U+b36/6G0SlNnWH9YeQmK6OOuA01rvZuS6fbSvz6HduT87/bl+01hgejdXm/wPQP/noJ6zJW0NCSALzB8/n8qzLvbIXu/B9MmDuK8njmzjroGy/2UIJS4DMydBQDf/o07YFFJYAk38L590I7kY4sg7isyAyWWaBnaXjFQ1sbSof4mrwkJwVwzmzepM6oHvtpaG15suiL3l82+OsL1pPrCOW34/9PdN7T/d2aMLHyIC5PwgKhV7Dza9mjki4vQgqjrQaU9kHUSnm+bID8Nx3ze/toRDfz2yhnHcTpJ5rDtajpQvsNH325kH2byim76gERs7oTUJa9ywBopRiTK8xjOk1hk3HNrFo2yISQhIAc88Rm8VGVLCUgRedwy9aHt1+J0FnHRRsMhNK6X4oaRqwn3M/9JsO+1fAC1dCTPqJbq/4/tB/FoTFeTt6r9Jak7+vks3Lcxh3aR/iU8KpLq0HIDI+cLtw7vzsTt47/B5XD7iaawddS6yje7W6RPtIt5U/dFt1hpJ9sP3VpuSyz2ypeJzwo0+h52DY/pq5jiW+n9n1Fd/PnGYcmWzOLOuGDI/Bwc0lbF6eQ0lODSERdqZ+fwAZwxO8HZpP2Fu+l0XbFrHi6AocNgdX9L+C6wZfR0KovD+BSLqtAlVCf8i+/cRtw2N2gUWnmbeVMsdNdrxuDuI3+7/D5gyxXW+bXWXx/c3kEpvp12tWtKF55W8bKMs/TnTPUKZ8L4ussYnY7LKArllWbBb/nPJPDlYe5IntT/Ds7mepc9fxx3F/9HZoopuQlkd3ojXUlpitk/LDMOpa8/jSX8DGp07cT1nNFsvN68zbOZ+bx+L7QUh0V0d9WuqqnRzYWMzQKckopdi+Oo+w6GAyhsX75BoNX5NTnUOQNYjEsES2lWzjjf1vsGDoAlIjUr0dmugC0vI4C67iYvJ/dSspD9yPLaGbNt2VMkuthPeA9Aknjs99EGbcdWIWWOk+s6XSbPkfTqx7Ce9ptlAyJsHk/zOP1ZWbU4u90AVWeayOLR/msOfzIjxug6R+UcSnRDB0SkqXx+LP0iLTWr7fW7GXtw++zZIDS5idOZsbht5ARlSGF6MT/ihgWh6Ff7qTypdfJvqqq+h1hzTd2yg/BMV7ToyplO4zu7a+s8g8f/9gaKg8MQssvh/0Ph96j++0kGqrGln70j4ObSnBarWQNTaREdNSiUkM67TnDCTFdcU8vfNpXt37Ko2eRq7IuoLfj/29t8MSnUQGzM8webiKizlwwTRwuVB2O31WrMCeKHt6nBatzS6vkqbpxSX7oDoPzr0BZv8TPG54bIKZbJoH6pvrggW3b4qsNjQ15Q1ExofgcRm8/Lf1ZA6PZ+jUFMKiZJpyZyirL+OZXc8QHRzN9UOux9AG+yv2kxWb5e3QRAeS5HGGyaPwT3dS+eqr4DGLIyqHg5irriJyzhwcgwfJxjvt1Xjc7PYKizPLtLz9s6ZZYAeheQe9C/4IE2+F4yWw5p62s8AierXpAnO7POz74hhbPszB7TT4/l/GYrFa0IaW8Ywu9uHRD/nl6l8yMXkiNw2/ieEJw7/9QcLnSfI4g+ThKi7m4PQZ6MZWffwWC1gs2BIS6PvRhyiLBU91NdbIyA6OOMB43E0LIfc1dW/1hfyN8Mw8aKw+cb+gCLjsSRqSs9m5Yg9b1x6jvs5CfGoYI2f0pu85PbtsDw3RVo2zhpf2vMQzu56hsrGSMb3GcNOwmxjdc7T8keXHJHmcQfIo/NOdVL7+Orha7SlttxN18VxirrySkGHD0G43+6dMxZ6cRNScuUTOmoktPr6Dow9gWsPxY03dXk0LIM+9gcOFsbz7yDbSgjYxImwJKY7dqNgMs3Vy4V/NRZH1lWYFZIck9q5U56rjlb2v8PTOp3HYHCy7dJmUgfdjkjzOIHkcmncpjXv2nHI8eMAAMpe8CYBRX0/FCy9QtXSZeV+rlbBx44i/5RZCR43skNgFlObVsHl5DpEJIYyZm4k2NOU5pcTZc08M1Dcnl+uWmnvNr70PVv7F7OpqWQTZH0b9QDbt6gIN7gbyavLoG9OXBncDt665lcv7X87klMnSEvEjAZs8urI8SeP+/VQtXUb1smX0+utdhI0bh/PIERoPHSZ8wvmoIP9dXOcNWmvydlewecVRcndXYA8299A4b85pTg3N3wSHVrfdZ8XdAP+vEKw2c5rxkY+bBuqbkkvz96JDHag4wE9W/oT84/lkxWSxcNhCpvWeJlvk+oGATR7NunKRoNYatEZZLBTf/wBlixZhjYoiYuZMoubOIWTUKJRUv/1Wn7x+gC0rcgiNDGJYdgpDJiUTHGo/8ws2l8MPa+pW/GIR7H3HnAVWU2Aei0qDX243v1/3T3OAvyW59G/3LDBxgstw8d7h9/jvtv9ypPoIfaL6sHjmYqmd5eMkeXhphbl2uTj+ySdUL11GzcqV6Pp6gvv1JeOttySBnMTZ4Gb3J4WkDY4lJjGMktwaSnJqyDovEau9k9+rxhqz66uhGvpMNY89912z5WK4T9wvazZc/YL5/Y7XITTeTCoRid22FlhH8xgeVhxdwcf5H/OX8/9irvov2c6AuAHYLWfxx4HoFJI8fKA8iVFbS83KlbjLyoibP9+s6vqLX+IYPJio2RdhT072anzeUlvVyLZVeexcm09jnZux8zI5Z2a6t8MyeVxmuZbSfWbXV1iCOWZiGPC3JHCbVXgJjjRbJyOuMdexAJQegJjeYJVfiN+koqGCGa/NINYRy4KhC5jXdx5BVuni9RWSPHwgeZzMU1ND7sKbqN+8GYCQc84hau4cIi68EFtMYOwWuPblfexcl4/h0fQZkcCIGWkkZvjBPhJaQ03hicWPzckl6yIY+yOzLMu9GWCxt10EOWAOJI/69uvXFMFr18NlT5uD/92Y1pp1+et4fNvjbCvZRo+QHlw/5Hq+2/+7sruhD5Dk4YPJo5kzL4/qZe9QtWwpzgMHSfrHvUTNnYvneC3KasES0n0+QFprSnOPE58ajlKKT17bj9tpMHxaKtE9Qr0dXsdprIHdS9sml/JDMOcBOOc6OLYTnr/81FL4vUaYU4uX/cpcnX/O9ebeLAFAa80XRV/w+NbH2XBsA29e/CZ9Y/p6O6yAJ8nDh5NHM601jXv3EpSWhiU0lLInF1P6n/8QMX0akXPmEDZuHMrmn/PlDUNzZGspm5Yf5djhai69dRRJ/aK9HVbXcjtBe8AeYk4rXnf/iWnGzuPmfb73OiQOgQeHmvuuWGww615IGgmxGQGzf/3ByoP0ie4DwF8++wsJoQlcM/AaIoNkzU5Xk+ThB8njZPXbt1P5yitUf7Aco7oaa1wcUXNm0+N3v/ObefIel8HuzwrZ8mEOVcX1RMY7GDEtjQHje2EPkj00gBNdYCV7zSTx0Z9h09Pm3isn+/GXZgvl4Co4+qnZJRabYf4bltDtBu3dhptfrf4Vq3JXEW4Pb9ndMMYRGEnUF0jy8MPk0cxwOqldu5aqpcvQTiepjz4CQOXrbxAycgTBmZlejvBUhqGxWBRup4dnbv+UiFgHI2f0JnNkgpQP+SY1RfDQcHM9SjNrsNndVV8O595oLnBccy+svhu0ceJ+QeHw630QFAb7P4Tq/BOJJSLJLK3jp/aU72nZ3TDEFsJ9k+9jUsokb4cVECR5+HHyaE1rjVIKT1UV+yZMBJcLx6BBRM6ZQ+Tsi7D39O7ganVpPVs+zKVgfwVX3H4eFouipryB8Jhgv2ktedWyX8HmZ80uq2bWIBh57aljH24nVOZAxWFzPKW6AKbfaZ579XrY+UarawRDr+Fwwwrz9sFVZhdabKa5tsXqH92hByoO8NTOp/jlOb8kPiSe/RX7iQiKIDEs0duhdVuSPLpJ8mjNVVxMzXvvUbV0GQ07doBSJP3jH0TNmd3lsRQfrWbz8hwObipGWRT9z+vJhMv7nd2ivkD02AQo2n7q8cShcPPHp38dwwNVeScSS/khs3vswr+a55+YBnnrze8tNohKNde3zHnAPJbzhbkzZHRvny7l8sMPfsiW4i3M6zuPBUMXkBwemNPdO5Mkj26YPFprPHyY6mXvEH3F5dh79qT6vfeofvddIufMJXzKZCzBnbenRd7eCt56YDNBDiuDJyYzLDuV8BjZQ8On1RyD8oMnEkv5IYhMPpFc7h9kdnuhzOOxGeY047E3m+eL90BUCgSHe+0lAOQfz2fx9sW8eeBNtNbMzpzNjcNupHdkb6/G1Z1I8ujmyeNkFa++SsnDD+MpKcUSHk7EjBlEzZ1D6NixZ9195HEb7F9/DI/bYPDEZAxDs2tdPv3PSyQoxD+6P8S3yNvQNrGUHza3F77gD2Z32V97muMtYT2aBu0zYfA86H+huYCyoRJCu67syLHaY+buhvte5RejfsH3B32/y567u5PkEWDJA0C73dR+8YVZGmXFCuxJSWQufRsAV0EBtl692pVInPVudq4rYOvKXGorG0nOimbeL09jwZvoXtxOsyZYS2I5Yv479mY4/+dQmQsPDgFH9InEEpsJgy42u98Mw5wV1gljYKX1pYTbw3HYHLy5/03W5q1l4bCFDIwb2OHPFSj8OnkopTKB24EorfVlrY6HAWuAP2mtl33TNQIxebRmNDTgKiggODMTo76e/edPwNazJ5Fz5xA1Zw5BaWnf+Pi9nxey9qV9OBs8JGfFMHJGGmmDYmUQXJygtZkQastg6wtma6U5wVTlwqWLYNjl5njKc9+BmIym2WBNM8L6XQiRvTosnJf2vMTDmx6mxlXDpJRJ3DTsJoYlDOuw6wcKryUPpdRiYA5QrLUe0ur4TOAhwAo8obW+5zSu9dpJyePPwHFglySP02c0NFD19ttUL11G3Xpz0NQxfBg9f/MbQkef+Bkpyz9OUIiNiFgHhQcq2bY6j5HT0+jRWxZqiXZyOwENtmBzxf2GxScSS+VRc3bZ9e9B7/Gw621YeVfbNSyxGZA2zpyG3A7VzuqW3Q2rGqu4ZsA13Dbmts55jd2UN5PHJMxf8M80Jw+llBXYB0wH8oD1wNWYieTuky7xQ611cdPjWpKHUmo6EAc4gFJJHmfGVVhI9bvvUrXU3IPEMWgQh5dvYeu6YgpK7QydnMykq7O8HWbAKa5u4Ccvbubf14ykR4TvznbqEIbHHJwP62HO7Dq0Gr78r7ktcfkhcNWZ9/v5NrPQ5KZnYeebJxJLTNO/8f2/di1L8+6GmdGZTEqZRLWzmh2lOxjXa5y0nr9FZySP0xoZ1VqvVUqln3T4POCA1vpQU3AvAZdore/GbKWcjilAGDAIqFdKvat16xVT4nTYe/UibsEC4hYs4NDmEjbes4HiozXYnQ1kHvuQ1C8t1PScKZtZdbGHP9rP+iPlPPzRAe6aN+TbH+DPLFaIbtVtmjnF/IIT2xCXHzZnd4E5UF9fDts2QGOVeUxZ4PYisASbiadoW5vEEhqbwfwh81ue4vV9r3P/xvsZFj+Mm4bfxMTkiZJEutDZTKtJBnJb3c4DxnzdnZVSccBfgZFKqdu01ndrrW9vOjcfs+VxSuJQSi0EFgKkfUuffqByuzzY7GaZkKM7SmmsdzP5mv6kOY5R+144Ne+9T957bxPUpw+Zy5bKB6yTVTe4qGtw8+rGPLSGl9fnMDI1msQoBw67lRC7lf49w7FZLTS4zPIlwTZL9/1/UcrcKyWi1SLAc64zv7SG+gozsdQUmF1iYHaD7X0PaktOPCa8p7n6HuDzx/heQwVhqbNZXPIFP/7oxwyMHchNw27igt4XdN1rC2CnPWDe1PJY1qrb6jJgptb6hqbb1wJjtNY/6aRYpdvqJHXVTravzmP7mjzm/GQ4iRlROOvd2IKtbcqHNG9m5amsJHrePLRhcPTaHxA6aiSRc+biyOrvxVfhX5xug2PVDSREBOOwW9l4tILXN+VRUFnf9NXA8UY380Yk8c72Qlyer/58bfnjdKJDg/j7+3t4dPVBlIIQu7Uluaz69RSCbBae/Pgwq/cWtxwPsVsJC7bxx7mDAFi1t5icsjrzsUFWHDYLEQ474/rEAXCsugGPoVuuHWyz+FdpmYbqpkWSh81yL8OvMo8/Mw8OrQLABbwTHsYT8T1JTRnLo9MehfVPoK0OVFwfs/US3qPb1QxrD691W32NfCC11e2UpmOik1UW17Hlw1z2fFaIx2WQPiy+pUDhV63RUHY7EVOmtNz2VFVhCQ+jbPFTlP33CYL79ydyzhyiLrkEe88eXfUyfI7WmrJaJwWV9aTEhBIbFsSO/CoeWX2A/MoGCivrKTneiNbw6s3jODc9lryKOj7YUURSdAjpcWGM7xNPpMPG42sPtUkcQVYL/7p6BMF2Kw0ug/Bg8/9pcv8EIhw2Gpwe6l1NX04Du9X8RecxDI43uik97qTB5aHe6cFqUS3J4/WNeSzbVtjmdfSICObL26cB8P/e2M5He4rbnM/qGcEHvzRrSt36ylZ2F1YTEmRtSTB9e4Tzu1kDAHjqk8NU1LmaEpeFkCArydGhTOhnbgO8I78KrSEkyNKS4MKCbTjsHVQw0xFplmTpNbzt8R8sAVc9VBzFXn6IeeWHmBscQdXguQDkrbmHn0RZmV9Vw+zjtdjtYXDOfJj5N/PxW15oWjiZaf7rxzXDvOVsksd6oJ9SKgMzaVwFXNMhUZ1EKTUXmNu3r+wL4HEbvP73jTgb3QwYk8iI6WnEJLZv9ootJoa0xx/HXV5urmJf9g4l999PcP9+2Hv2wF1eDkp1u82s6pxuCiobKKisp3dcKL3jwjhUcpzfL9lhthqqGnC6zZ7Th64awSUjkml0G+wpqiE5OoSsrASSokNIig6hd6y5b8nFw5O4ZETbchq/f3M7xkkteo1m3YGyU8Y+xmbGMTYz7mtjXjipDwsn9fna8/ddPpw7Lx5MvctDg8tDg8vAY5x47h9OyGD6oJ5m4nEZ1Ls8RASf+NgnRTuoqnc1nfdQXuvE0Wq74Dc25bM9v6rNc07sF9+SPG5+biN5FfVtzs8Y1JNFPzD/yM2+bzV1Tg8hQWarJyTIyrSBPfnxVPOzfNsb27FZFA67paX1NCIlmvF94zEMzfJdRSdaXU0JLj48mJiwILTNgSeuP7YeZqKzAs1LGquufR37p3/gD0GHeSwpgx860piXkEUQmHvaL/nRiYCtweYg/thbYPT15qyyw2vNwfzoNNlF8muc7myrFzEHt+OBY8AdWusnlVIXAQ9i/r8t1lr/tfNCDcxuK21ojuwo4+DGYi64biDKosjdVU5schhhUR1XPsSZl4e9Z0+U3U7xAw9S9uSThE+YQOTcOURkZ/v8ZlZuj0FxTSOFVfXkVzaQFhvKiNRoSo838oMnv6Sgqp7KOlfL/W+bNYCbJvehoLKen7ywiV7RISRHh5AU5SApOoQRqdH0iDyzGVIXPbSOXYXVpxwf1CuSd38+8Yxfo7cYhsbpMahvah1ZlCIxynxvvjhURnWD20xeTg8Nbg/J0SFcMNAs7vnXd3ZRVe+i3mU0JTcP5/eN5+bJfdBaM+Hvq1oSX73Lg9awYEIGf5gziDqnm0F//OCUeH6a3ZdbZ2RReryR0Xd9SJDVQnBT8gkJsvKjyX246rw0jlXVc/Mbz3HMsowaDhGqEpgV8wAXD0thVFQt1QX72b1zCzGNuUTW5VGWfhE1fS8hy1pIzFPnA6CVFR2VCrGZWM7/mVk7rPG4WXMsJr3Laoad7cw9v14k2BECKXl4XAZ7vyxiy4ocKorqiIh1MO9XI4mM7/xf4g379lG15C2q33kH97FjqNBQombPJvHPd3plUFdrTXW9m4Kq+paxhcSoEKYP6onH0Ez+xyoKqxra/MU9f3w6f7p4ME63wc3PbSQp2kwKSVFmy6FPQhhx4VK7y5dorWlsavk57FY8hmZvUU2rVpWZYPokhDOwVyQ1DS4Wf3yEBrfZndd8/uLhSVwwsCe55XXc8vwm6lxujrOLRks+rvKJ3DF3MBHxO4hhOFc9vvmUOP5zxUBmxx1j/56tvLfuM9LVMXqrYzzkuZz1tnN4PruOYavmo1GUWeI4ZkuiNCiZlTFXUh2ewU8nJpEZF8Lucs1Hu4/haDWWFRJk5fy+8USF2Ck73kjJ8caWsazgpn/tVnXK5+z3b27n+S9z+N6Y3mc0cy9gk0erbqsb9+/f7+1wOl1FUS1L7t9MXbWT+NRwRs5Io++oHlisXdsvqw2DuvUbqF62FFD0+sufASh74glCR4/GMXx4hySTBpeHoqqGpuRgdivFhNq5dlw6AFP+sYojZXVtHjNzcCKPXXsOAHe8tYMIh72pS8lhtiCiQwgLPpteWdFdHaw8yLy35hEdHM2lfa7mwpTvYiGkJUn16xFBQkQwhVX1fLS7uGWsqb6pW/DaocFkVG8i//BODu/bQbwzn57uAm4Pu4PtRjovjj5Iytpf0xAUy46GeI7qHhw1EnnWM40KIln+8/H0T4zmqU+PcOfSXafEt+7/ppIaG8rTnxzmiY8PY7cojpTXoTU4bBbW/nZqu1sfAZs8mnXnlkdNeQOVx+pIHRiL4TFY+cwessYmkjIgxqemcLpLSzlwwTR0YyP2tDSi5swmcs6cb9zMqux4I7kV9RRW1pPfNCPJblXcdpFZq+iSf3/M1ry2/eoT+sbz3A3mzO8n1h0CaBlvSIp2EB8W7F+zhoRP2Vqylf9u+y9r8tYQYY/gmoHXMH/wfMKDOqDC8LFdsP8DKD+E0VTexVKdT+6CbdTZY8nc9W/sXz5KY2Q6ZUHJVIekUOFIYXf8TOo8Fuafn0F4sI0Vu47x3vZCNhwtp6G8gH8F/YtfeH7OBecOa3frQ5JHN0wepXk1bF6ew/4NxYRFBXHtX8f7/C9FT00NNSs+pHrZUmo//wIMg6rf/pnDg8dQ0JQgahvdPH6t+bN607Mb+GDnsZbHh9itDEmO5NWbxwOwbFsB9U5PS4uheT2EEJ1td9luFm1bxPpj63n/O+8THhSOoQ0sqoNb+a4Gcw2LUnDgQ3MNS3Ol48occz+W24vMWV/v/RYOr4PYDGrDe/P3LxqZwGamWTbxnOcC/sYN7W59BGzy6I7dVsVHq/n8rUPk7irHHmxl0IQkhl+QSkSs75WxOFpWy5bcypYupeaZSUt+PB5LeRmv/uMp/ubJ5HhQKHOOfEZ20Xb2DhnPrXfdTFB0FF8eLqeq3tXSpRQVYvep1pQQVY1VRAVHYWiD77/7fYYnDGf+4Pn0DOuCXT09LqgpPLFCf/0T5hbE5Ydwlx3Cpl0YGiwK6nUQ2e6H2t36CNjk0czfWx4ej4HHaRAUYiNvTzkrFu9iWHYKgycm4wjr2umAWmu0BotFcaS0lrX7S1q6lAoqzS6mlxaOIy0ulEVrD/K3d/cAEB1qp1dUCMnRDv5x2XBiwoI4UHycyjonSdEhOD58l/L/LsJ1NAcVFET45MlEXjyXyOnTu/T1CXEmal213P3F3Sw7tAyLsnBp30tZMHQBSeFJXolnzoOr+WHpvcyxfk6Q8tCobbzsmcJLCb9o18w9SR5+mjycDW52f1LIlo9yyByRwMQr+qO1xnBrrPbOGQRvLnvhsFvJLa/jjU35TS0Gs1upsLKBp64/l7GZcSzdWsBPX9xMkNVCr2hHy4ykX0zrR2psKMU1DVTXu+gVdXqD0FprGnbsoGrpUqrffY+g3r1Jf/45M649ewju1w9llW4p4bvyavJ4cseTLDmwBDQsnrmYkT1Gdn0gNUXw0HBzdX0zm8MsMBlx+q0iSR5+ljxqqxrZtiqPnWvzaaxz06tvFOfMSqf34K9fFHY6DENTerwRu9VCTFgQx6obeGzNwZbyGIVV9ZQed/LAlcO5dGQKG4+W891HPyMhIrhlHUNSdAhXn5dK3x4R1Da6qXN6iAsL6vDxFu124y4rNxcfVlSwf+IkbLGxRF50EZFz5+AYNEi6sITPKqot4rV9r3Hz8JuxWWyszVtLcngyfaK/fuFmh1r2K9j8rFnuvpk1CEZeC3PuP+3LBGzy8Ncxj4+e2c2ezwrpMyKBETPSSMyIOq3H1TS4KKxqwGGzkhYXSk2Dizve2ml2K1XVU1TVgMuj+b+ZWdwypS/5lfVMv38NydEhTYvdzNbD9ME9GZAYidtj4NGaYJt3/9o3nE6Or1xJ1dJlHF+7FlwugjIy6PWXP7fZg0QIX2Rog7lvziW3JpdpvaexcNhCBsQO6NwnfWwCFG0/9XjiULj549O+TMAmj2a+3PLQWlN4sIrNy3M4d3Y6PXpHUl1aj+HRRPcMbbmfy2NQVNVAYVUDwTYLw1Oj0Vpz4zMbyauoI7+ynpoGNwDXju3NX+YNwe0xmPrP1SRGOlpNVw3h3PQYBiRG0vx/6E9/wXsqK6levpzqpj1IgtLSqP3iSxr37SNy1kxs8fHeDlGIU1Q0VPDsrmd5cc+LHHcdZ0rKFH466qf0j/Ht4qKSPHwweRiG5vDWEjYvz+HY4WqCQm30vSgNZy8HVqWYNsjsl7zxmQ1sy6ukuMYsrAcwbWAPnrjuXAB+sPhLgqwWs9XQlBwGJEbQr2eEt15alzt29z2U/+9/YLUSNm4ckXNmEzFtOtbw9tXuEqKzVTureWH3Czy3+znum3wfY3uNxWW4sFt8sw6WJA8vJ48Gl6dlXKGgsp5GtwfH6hKKj9bgdFj41NrIFpsbV1MDoHUtozve2kGt00NSc7dSUxXW1NjQb3jGwNO4fz9Vy96hetkyXPn5BPfvT+bbbwFm686fWlei+6tz1RFiC0Epxb3r72V32W5uGn4TYxLH+NTPqiSPs0ge31ZYzGNoSmoa29RPKqt1ctsscxX0r1/dymsb83AYMNBlZXOQh9jwIJ6ZNJDgMDura49zrKaxpdWQ3LQSWuonnRmtNfWbt+CpriJiyhS008nBuXMJGzuOqLlzCBk1CiVltIUPeWXvKzy+9XGK64sZnjCchcMW+szuhgGbPDpiwPx3r2/j5fW5TM1KIHtgz5YEcc93h+GwW7lr2S6e+Phwm8dEBNvY8IdpBNusvPXJUfI+O4Y6Uot2aybeMoQBg+IJsskvsK7gLivj2N/upmblSnR9PbakXkTNnk3MVVdhT07+9gsI0QUaPY0s2b+EJ3c8SWFtITcOvZGfjfqZt8MK3OTR7ExbHsXVDYy/ZyXuVlVXbRazrPQrN40jKTqEDUfKW/ZtSIoOoVe0g0iHnfoaJ2tf2sfBTcUoi6L/eT0ZMT2NuKQOqIEj2s2oraVm5Uqqli6l9pNP6f3sM4SOGoWroAC0lkQifILL42LZoWUMTxhOZnQmBysPsq9iHzN6z8Bq6fpZj762k6DfePij/dDUcrRZFJeMSOLey4ZjbbWmYXR6LKPTza1ktNbUlDeAw05QiI2yglpGTEtjWHYq4THSDeVNlrAwoubOJWruXNzl5VibNqwqe+IJKl54kZBzziFq7hwiLryw221mJfyH3Wrn0n6Xttx+bd9rPLf7OR6JfIQbht7ARZkX+ezg+unq9i2P4uoGJt67qmWfAPj6ssYet8H+DcfYvDyHxjo31941DqvNgjY0yseLFQY6Z14e1cuWUbV0Gc6DB8FmI2r2bJL+fo+3QxMCQxt8ePRDFm1bxN6KvSSHJ3PLiFu4uM/FXfL80vI4Aw9/tP+ULUE9WvPwRwdaCos5693sXFfA1pW51FY2EpsUxth5mS2tFUkcvi8oJYX4m28m7qabaNyzh6ply7AEm61ErTXF9/ydsPPHEzZ+PMrW7X/shY+xKAsz0mcwvfd01uSt4fGtj3Og4gBg/nw6DSfBVv/q1ej2n6JNOZW4PG2Th8uj2XS0ouV20aEqPn3jAMlZ0Uz9/gDSBsf6xAwJ0X5KKRwDB+IYOLDlmLuwkMo336T8f//DGhtL5KxZRM2d02GbWQlxupRSTEmdwuSUybgMc1vkTwo+4Q+f/IH5g+dzef/LCbX7x/R9v+i26ojZVrVVjSx/YgczbhhCw3EXWz7MISwqmLHzzL2Uy/KPE58SOAvyAo3hdFK7dq1ZGmXVKrTTScoj/yEiOxvt8UihRuE1O0t38sCmB/ii8AtigmP4weAfcFXWVR2zMVUTmW11Fus8Vr+wh51rCwiPDeZ4eSO2IAvDslMZN6+LCpwJn9G8mVXkRbOwOByUPvqoeXvuXCIvmoW9Zxfs4SDESbYUb2HRtkWsy19H3+i+vHHxGy0t45K6En6z9jfcN/k+4kPaX7pHkscZJo/aqkb+d9un6KapuiOmp3LOhek4wv17toPoGFVLl1H+v//RsGMHKEXomDFEzbuE6HnzvB2aCEA7y3ZSVl/GpJRJuDwuFu9YTE5NDksPLuWKrCv4/djft/uanZE8AmKF24Z3DtPctW2xKtxOQxKHaBE1dw4Zr71K5nvvEn/LLbgKC6j5YHnL+drPP8dobPRihCKQDI4bzKSUSQBsKt7Ev7f8m7cPvo1Gs+TAEkrrS70coanbJ4/aqkZ2f1aE0TRobng0uz8tpLZKfhmItoIzMkj46U/o8/77JN37dwBc+fnkzL+e/RMmUnD77dR+/jna4/FypCJQjOk1hlnps7Aqc0zO0AaPbX3My1GZun3y2PDO4Zbuqmba0Gx494h3AhI+TymFNcKcPGHr2ZPUJ58g4oILqHn/A3LmX8+BqdnUbdzo5ShFICipK2Fl7ko82vyDxWW4fKb10e2TR9Gh6pZWRzPDoyk6WOWliIQ/UTYb4eefT9I9d9Pvk49JfuB+HEOHEpSeDkDNhx9S8sgjOHNyvBuo6JYe2/YYhjbaHPOV1odfrPNoNVW33Y+98vfndXxAIiBZHA4iZ80ictaslmN1GzZS/vTTlD78LxzDhxE1Z27LZlau4mLyf3UrKQ/cjy0hwYuRC3+1tXhry3qQZi7DxZbiLd4JqJWAmG0lRGdyFRZS/e67VC1dRuOePTiGDSPjlZcp/NOdVL78MtFXXUmvO+7wdpgigMlUXUkewsc17t+Pp6YGe0oKB6fPQDfP0rLZUEFBKLvdLKNy/XxcxcXkzL8eZbe3+Yr5/veInDEDV2Ehxf/4B8puh1bno+bMIWTYMFzHjlH11tttH2+zETpmDEEpybjLy2nYseOU69vTemMND8Oor8dTXYMKante9knpfqS2lRA+LrhfPwAK/3Qn2mjqq7ZYCO7fn7AxY9AuF8FN3a/KaiU4qz/a5UK7XOByoZ0uaPp7zqivp2HX7pbz2u1Gu1yEDB1qJo/8fEruv/+UGJIffoiglGQadu4id+FNp5xPffIJws8/n+Nr1pL/i1+ccr73Cy8QOmokVUuXcezuu9skJmW3k/LvfxHUuzfVK1ZQ+eKLbRKbstvpedtt2GJiOP7JJ9R+8mmbxyq7nZjvXYPF4aB+x06chw62eSw2G2HjxqEsFlxFRXiqq5vOmYlXBdlbqiUH0s6SvtgFKslDiA7mKi6m6s03wdXUV20YOA8eJO3xx9p88G1xcaQ88MDXXic4M5M+77/3tedDRowga+uWE8nFaf5ri4luOj+c9JdfOnG+6csxYAAAjkEDSfzTn9okJu1yYk/qBYA9OZmIC2ecSGxNX6qp4CRuN0Z9A7q65pTnAGjYtYuKl15qeXyz6CsuB4eD6nffpXzx4lNe14CdOwAofewxKl96uc055XAwYMtmAAr+77dUv/NOm+RjS0ho2ba46G9/o27DhratrqQkkv76V/P6i/6L88iRtud79SL2B9cCUPX223gqK0+5ftj48QDUbdpsvh/N54PsWCMjsfcy3z93RQXKaj2RGK3WM052pY88Sv3GjZQ88ii97vjjGV2jo0m3lRAdrPBPd1L5+uttfmFitxN92WU+88HvalrrlgSkQkNRSuGprMRTVXVK4gk95xwAGnbvxnk0p1XLy2yVxVx5BQDV739Aw949bRObI4Se//cbwEw+9Vu2trm2rUcPUh56EID8X91K3ebNbRNrv36kv/wSAIe+8x0ad+1u8zpCzz2X3s8+A8DBC2fiPHq0zfnwKVNIfexRAPZNnIinpNWUWqWInDuH5HvvBeDAhRe2TT72ICIvvJD4m29Ca03uwptQdjvaMKhduxYMAxUcTN8PV7S79SHdVkL4gfotW9omDgCXi/rNm70Sjy9QSkFQECooqOWYNToaa3T01z7m5OrIJ4uceSGRMy/82vPxN9/8jTEl3//Pbzyf/vzzaKezTXJpXc4/6b77MGpr2yan+LiW8z1+8QuM48fbnG/u1gQIP/98s+XW6rwloqkYosfTklhd+fnQ1AWqDcNnWh/S8hBCCB/lKi5uO/ECzqj1IbWthBAigJQ+8uiJiRdNmlsf3ibJQwghfJQvd4HKmIcQQviozCVvejuEr+UXLQ+l1Fyl1KKqKqlHJYQQvsAvkofWeqnWemFUVJS3QxFCCIGfJA8hhBC+RZKHEEKIdpPkIYQQot38apGgUqoEOPqtd/x6UUBnj7p31nN09HU76nrxgPe3NRMdoSs+H/7C39+Lk+PvrbXu0IqKfpU8zpZSapHWeqE/PkdHX7ejrqeU2tDRK1eFd3TF58Nf+Pt70RXxB1q31VI/fo6Ovm5XvBfCv8jPxAn+/l50evwB1fIQHU9aHkIEpkBreYiOt8jbAQghup60PIQQQrSbtDyEEEK0myQPIYQQ7SZVddtBKZUJ3A5Eaa0v83Y8QnibUioMeARwAqu11s97OSSv8ff3or3x+3TLQykVrZR6TSm1Rym1Wyk17gyvs1gpVayU2vEV52YqpfYqpQ4opX73TdfRWh/SWi84kxgCgVIqTCn1P6XUf5VS3/N2PIFCKZWqlFqllNqllNqplPr5WVyrvZ+V7wCvaa1vBC4+0+ftKEoph1LqS6XU1qb34s6zuJbX3gullFUptVkptewsrtGp8ft08gAeAt7XWg8AhgNtdqNXSvVQSkWcdKzvV1znaWDmyQeVUlbgP8AsYBBwtVJqkFJqqFJq2UlfPTrmJfmXr/sB9IdfJAHEDdyqtR4EjAV+rJQa1PoOnfVZAVKA3Ka7ec7ydXSERiBbaz0cGAHMVEqNbX0HP3kvfs5Jv+9aPb9PxO+zyUMpFQVMAp4E0Fo7tdaVJ91tMrBEKRXc9JgbgX+dfC2t9Vqg/Cue5jzgQFOLwgm8BFyitd6utZ5z0ldxh704//I0J/0A+tEvkoCgtS7UWm9q+r4G85dO8kl365TPCpCH+f8OPvD7RJuON920N32dPKXUp98LpVQKMBt44mvu4hPxe/0/+xtkACXAU03Ntyea+uRaaK1fBT4AXm7qJvkhcHk7niOZE7/swHzzTv7QtVBKxSmlHgNGKqVua8fz+K2v+QH0i18kgUgplQ6MBL5ofbwTPytvAN9VSj2Kj6zKbury2QIUAyu01v72XjwI/B9gfNVJX4nflwfMbcAo4Kda6y+UUg8BvwP+0PpOWut7lVIvAY8CfVr91dHhtNZlwM2ddX0/8lU/fGOAh4F/K6Vm4yO/SAKJUioceB34hda6+uTznfFZ0VrXAtef7XU6ktbaA4xQSkUDbyqlhmitd5x0H598L5RSc4BirfVGpdSUb3gur8fvy38d5gF5rf5qeA0zmbShlJoIDAHeBO5o53PkA6mtbqc0HRNnQGtdq7W+Xmv9I3+baeLvlFJ2zMTxvNb6ja+5T0B9Vpq6uVfx1f3+vvpenA9crJQ6gtmiz1ZKPXfynXwhfp9NHlrrIiBXKZXVdOgCYFfr+yilRmKWx7gEM2PGKaXuasfTrAf6KaUylFJBwFXA22cdfPfnd79IujOllMIcG9yttb7/a+4TEJ8VpVRCU4sDpVQIMB3Yc9J9fPa90FrfprVO0VqnN113pdb6+z4Zv9baZ78wZ0tsALYBS4CYk86fDwxtddsO3PgV13kRKARcmC2aBa3OXQTsAw4Ct3v7NfviF5AO7Gh12wYcwhyXCgK2AoO9HWegfgETMAeFtwFbmr4uOuk+AfFZAYYBm5veix3AH7/iPn7xXgBTgGW+Gr/UthLfSCn1IuYPcTxwDLhDa/2kUuoizIE9K7BYa/1XrwUphOhykjyEEEK0m8+OeQghhPBdkjyEEEK0myQPIYQQ7SbJQwghRLtJ8hBCCNFukjyEEEK0myQPIYQQ7SbJQwghRLtJ8hBCCNFu/x+yIy/wzHiW4AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "log_mean = lambda x1,x2 : 10**np.mean((np.log10(x1),np.log10(x2)))\n",
    "\n",
    "x1,x2,x3 = log_mean(0.3,1.2),log_mean(1.2,2),log_mean(2,7)\n",
    "colors = iter(['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple',\n",
    "               'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive', 'tab:cyan'])\n",
    "\n",
    "for i,row in xray.loc[xray[\"X-ray source\"]==\"XMM\"].sample(5).iterrows():\n",
    "    c = next(colors)\n",
    "    plt.loglog([x1,x2,x3],row[[\"from_XMMs\",\"from_XMMm\",\"from_XMMh\"]],\n",
    "               marker=\"s\",linestyle=\"--\",color=c,label=row[\"Flag\"]+str(i))\n",
    "    plt.loglog([x1,x2,x3],row[[\"Flux (s band)\",\"Flux (m band)\",\"Flux (h band)\"]],\n",
    "               marker=\"^\",linestyle=\"--\",color=c,label=\"\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "master[\"vsini\"] = master[\"vsini\"].apply(pd.to_numeric, errors=\"coerce\")\n",
    "master[\"HR1\"] = (master[\"Flux (m band)\"]-master[\"Flux (s band)\"])/(master[\"Flux (m band)\"]+master[\"Flux (s band)\"])\n",
    "master[\"HR2\"] = (master[\"Flux (h band)\"]-master[\"Flux (m band)\"])/(master[\"Flux (h band)\"]+master[\"Flux (m band)\"])\n",
    "\n",
    "XDBS = master[master[\"X-ray source\"].notna()].drop(['Category', 'RA', 'DEC', 'Nb of spectra in BeSS', '_q',\n",
    "                                             'rgeo', 'b_rgeo', 'B_rgeo', '_4XMM', 'Flux1', 'e_Flux1',\n",
    "                                             'Flux2', 'e_Flux2', 'Flux3', 'e_Flux3', 'Flux4', 'e_Flux4',\n",
    "                                             'Flux5', 'e_Flux5', 'Flux8', 'e_Flux8', 'Flux9', 'e_Flux9', 'e_HR1',\n",
    "                                             'e_HR2', 'HR3', 'e_HR3', 'HR4', 'e_HR4', '_2CXO',\n",
    "                                             'Fluxb', 'Fluxh', 'Fluxm', 'Fluxs', 'Fluxu', 'Fluxw', 'HRhm',\n",
    "                                             'HRhs', 'HRms', 'GamPL', '_2SXPS', 'CR0', 'E_CR0', 'e_CR0', 'SwiftHR1',\n",
    "                                             'SwiftHR2', 'CR1', 'e_CR1', 'CR2', 'E_CR2', 'e_CR2', 'CR3', 'E_CR3',\n",
    "                                             'e_CR3', 'PcstS0', 'FPO0', 'FPU0', 'Gam1', 'E_Gam1', 'e_Gam1',\n",
    "                                             'XRTsoft', 'XRTmed', 'XRThard'],\n",
    "                                              axis=1).rename(columns={'_RAJ2000':\"RA\", '_DEJ2000':\"Dec\",\n",
    "                                                                    \"NHPL\":\"N_H\",\"V\":\"Vmag\"}).copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in XDBS[XDBS[\"Simbad classification\"].isna()].index:\n",
    "    XDBS.loc[i,\"Simbad classification\"] = str(Simbad.query_object(XDBS.loc[i,\"Be star\"])[\"OTYPE\"][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# for i in master.index:\n",
    "    flags = []\n",
    "    if pd.notna(master.loc[i,\"Be star\"]):\n",
    "        flags.append(\"BeSS\")\n",
    "    if pd.notna(master.loc[i,\"EDR3Name\"]):\n",
    "        flags.append(\"Gaia\")\n",
    "#     if pd.notna(master.loc[i,\"Name\"]):\n",
    "#         flags.append(\"HMXB\")\n",
    "    if pd.notna(master.loc[i,\"_4XMM\"]):\n",
    "        flags.append(\"XMM\")\n",
    "    if pd.notna(master.loc[i,\"_2CXO\"]):\n",
    "        flags.append(\"Chandra\")\n",
    "    master.loc[i,\"Memberships\"] = \",\".join(flags)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7,7))\n",
    "plt.scatter(master[\"BP_c\"]-master[\"RP_c\"], master[\"G_c\"]-5*np.log10(master[\"rpgeo\"]/10),marker=\".\")\n",
    "plt.ylim(16,-7)\n",
    "plt.xlim(-1,5)\n",
    "plt.vlines([-0.4769230769230768,-0.03076923076923066],-10,20) # WebPlotDigitized bounds of B type from Gaia HR diagram\n",
    "plt.xlabel(\"$\\mathrm{G}_\\mathrm{BP}-\\mathrm{G}_\\mathrm{RP}$\")\n",
    "plt.ylabel(\"$\\mathrm{M}_\\mathrm{G}$\")\n",
    "plt.title(\"Color-magnitude diagram of Gaia-matched Be stars\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# XDBS = pd.read_csv(\"./XDBS_catalog.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# XDBS.to_csv(\"./XDBS_catalog.csv\",index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots()\n",
    "for classification,data in XDBS.groupby(\"newclass\"):\n",
    "    ax.scatter(data=data,x=\"Flux (0.3-1 keV)\",y=\"Flux (1-2 keV)\",marker=\".\",label=classification)\n",
    "plt.yscale(\"log\")\n",
    "plt.xscale(\"log\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mapping = {'HMXB':\"**\", 'Em*':\"*\", 'Be*':\"*\", 'Star':\"*\", '**':\"**\", 'V*':\"*\", 'YSO':\"YSO\", '*':\"*\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "XDBS[\"newclass\"] = [mapping[cls] for cls in XDBS[\"Simbad classification\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.9.1"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}