build_syssim_catalogue.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "import csv\n",
    "\n",
    "MAINPATH = '/Users/research/projects/maiasaurus/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Build the observed catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# READ IN OBSERVED PLANETS\n",
    "obs_file = MAINPATH + 'catalogues/syssim_observed_planets_raw.csv'\n",
    "\n",
    "obs_data = []\n",
    "with open(obs_file) as infile:\n",
    "    reader = csv.reader(infile)\n",
    "    \n",
    "    for row in reader:\n",
    "        obs_data.append(row)\n",
    "        \n",
    "KEYS = obs_data[20]\n",
    "VALUES = obs_data[21:]\n",
    "\n",
    "def getkey(keyname,keys=KEYS,values=VALUES, convert_to_float=True):\n",
    "    '''\n",
    "    keyname = (string) of column definition\n",
    "    '''\n",
    "    kid = keys.index(keyname)\n",
    "    \n",
    "    dummy = []\n",
    "    for row in values:\n",
    "        dummy.append(row[kid])\n",
    "        \n",
    "    if convert_to_float:\n",
    "        return np.array(dummy, dtype='float')\n",
    "    else:   \n",
    "        return np.array(dummy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = {}\n",
    "\n",
    "obs['target_id'] = getkey('target_id')\n",
    "obs['star_id']   = getkey('star_id')\n",
    "\n",
    "obs['period']    = getkey('period')\n",
    "obs['depth']     = getkey('depth')\n",
    "obs['duration']  = getkey('duration')\n",
    "\n",
    "obs['star_mass']   = getkey('star_mass')\n",
    "obs['star_radius']   = getkey('star_radius')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEDtJREFUeJzt3X+snmV9x/H3Z62g4kaxnBBoG9rEzoWZbJATxLAYY6eCEssfajCbdgzSLKkOxxIF9wfZFhPNFlEzQtK0aMkYSFBDY5jaAcb5B8gpMn4VwwmCPR3YI7/8QYQVv/vjXGzH0nLgPOc8Dz3X+5WcPNd93ddz39+7f5xP7+v+cVJVSJL68zujLkCSNBoGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTy0ddwEs5/vjja+3ataMuQ5KOKLt37/5ZVY3NNe5VHQBr165lYmJi1GVI0hElySMvZ5xTQJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1Kk5AyDJVUn2J7l3Vt8/JXkgyd1JvpFkxax1lyaZTPKjJO+Z1X9W65tMcsnCH8qLrVqzliRd/axas3YY/7SSloCX8yqIrwD/Alw9q28XcGlVHUjyOeBS4FNJTgHOA/4QOAn4jyS/375zBfAuYAq4I8nOqrp/YQ7j0P576hEuvOLxxdzFq862LStHXYKkI8ScZwBV9T3giYP6vlNVB9ribcDq1t4IXFdVz1bVj4FJ4PT2M1lVD1XVc8B1bawkaUQW4hrAXwL/3tqrgL2z1k21vsP1v0iSzUkmkkxMT08vQHmSpEMZKACS/B1wALhmYcqBqtpaVeNVNT42NufbTCVJ8zTv10En+QvgHGBDVVXr3gesmTVsdevjJfolSSMwrzOAJGcBnwTeX1XPzFq1EzgvydFJ1gHrgR8AdwDrk6xLchQzF4p3Dla6JGkQc54BJLkWeAdwfJIp4DJm7vo5GtiVBOC2qvqrqrovyfXA/cxMDW2pqufbdj4GfBtYBlxVVfctwvFIkl6mOQOgqj58iO7tLzH+M8BnDtF/E3DTK6pOkrRofBJYkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjo1ZwAkuSrJ/iT3zup7Y5JdSR5sn8e1/iT5UpLJJHcnOW3Wdza18Q8m2bQ4hyNJerlezhnAV4CzDuq7BLi5qtYDN7dlgLOB9e1nM3AlzAQGcBnwVuB04LIXQkOSNBpzBkBVfQ944qDujcCO1t4BnDur/+qacRuwIsmJwHuAXVX1RFU9CezixaEiSRqi+V4DOKGqHm3tx4ATWnsVsHfWuKnWd7j+F0myOclEkonp6el5lidJmsvAF4GrqoBagFpe2N7WqhqvqvGxsbGF2qwk6SDzDYCftqkd2uf+1r8PWDNr3OrWd7h+SdKIzDcAdgIv3MmzCbhxVv9H291AZwBPt6mibwPvTnJcu/j77tYnSRqR5XMNSHIt8A7g+CRTzNzN81ng+iQXAI8AH2rDbwLeC0wCzwDnA1TVE0n+EbijjfuHqjr4wrIkaYjmDICq+vBhVm04xNgCthxmO1cBV72i6iRJi8YngSWpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwMFQJK/SXJfknuTXJvktUnWJbk9yWSSryY5qo09ui1PtvVrF+IAJEnzM+8ASLIK+GtgvKreAiwDzgM+B1xeVW8CngQuaF+5AHiy9V/exkmSRmTQKaDlwOuSLAdeDzwKvBO4oa3fAZzb2hvbMm39hiQZcP+SpHmadwBU1T7gn4GfMPOL/2lgN/BUVR1ow6aAVa29CtjbvnugjV958HaTbE4ykWRienp6vuVJkuYwyBTQccz8r34dcBJwDHDWoAVV1daqGq+q8bGxsUE3J0k6jEGmgP4U+HFVTVfV/wBfB84EVrQpIYDVwL7W3gesAWjrjwUeH2D/kqQBDBIAPwHOSPL6Npe/AbgfuBX4QBuzCbixtXe2Zdr6W6qqBti/JGkAg1wDuJ2Zi7l3Ave0bW0FPgVcnGSSmTn+7e0r24GVrf9i4JIB6pYkDWj53EMOr6ouAy47qPsh4PRDjP018MFB9idJWjg+CSxJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASerUQAGQZEWSG5I8kGRPkrcleWOSXUkebJ/HtbFJ8qUkk0nuTnLawhyCJGk+Bj0D+CLwrar6A+CPgD3AJcDNVbUeuLktA5wNrG8/m4ErB9y3JGkA8w6AJMcCbwe2A1TVc1X1FLAR2NGG7QDObe2NwNU14zZgRZIT5125JGkgg5wBrAOmgS8n+WGSbUmOAU6oqkfbmMeAE1p7FbB31venWp8kaQQGCYDlwGnAlVV1KvAr/n+6B4CqKqBeyUaTbE4ykWRienp6gPIkSS9lkACYAqaq6va2fAMzgfDTF6Z22uf+tn4fsGbW91e3vt9SVVuraryqxsfGxgYoT5L0UuYdAFX1GLA3yZtb1wbgfmAnsKn1bQJubO2dwEfb3UBnAE/PmiqSJA3Z8gG//3HgmiRHAQ8B5zMTKtcnuQB4BPhQG3sT8F5gEnimjZUkjchAAVBVdwHjh1i14RBjC9gyyP4kSQvHJ4ElqVMGgCR1ygCQpE4NehFYrzLLlh9NklGXMTQnrT6ZfXsfHnUZ0hHJAFhinj/wLBde8fioyxiabVtWjroE6YjlFJAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjo1cAAkWZbkh0m+2ZbXJbk9yWSSryY5qvUf3ZYn2/q1g+5bkjR/C3EGcBGwZ9by54DLq+pNwJPABa3/AuDJ1n95GydJGpGBAiDJauB9wLa2HOCdwA1tyA7g3Nbe2JZp6ze08ZKkERj0DOALwCeB37TllcBTVXWgLU8Bq1p7FbAXoK1/uo2XJI3AvAMgyTnA/qravYD1kGRzkokkE9PT0wu5aUnSLIOcAZwJvD/Jw8B1zEz9fBFYkWR5G7Ma2Nfa+4A1AG39scDjB2+0qrZW1XhVjY+NjQ1QniTppcw7AKrq0qpaXVVrgfOAW6rqz4BbgQ+0YZuAG1t7Z1umrb+lqmq++5ckDWYxngP4FHBxkklm5vi3t/7twMrWfzFwySLsW5L0Mi2fe8jcquq7wHdb+yHg9EOM+TXwwYXYnyRpcD4JLEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ1akL8IJo3KsuVHk2TUZQzVSatPZt/eh0ddhpYAA0BHtOcPPMuFVzw+6jKGatuWlaMuQUuEU0CS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpU/MOgCRrktya5P4k9yW5qPW/McmuJA+2z+Naf5J8KclkkruTnLZQByFJeuUGOQM4APxtVZ0CnAFsSXIKcAlwc1WtB25uywBnA+vbz2bgygH2LUka0LwDoKoerao7W/sXwB5gFbAR2NGG7QDObe2NwNU14zZgRZIT5125JGkgC3INIMla4FTgduCEqnq0rXoMOKG1VwF7Z31tqvVJkkZg4ABI8gbga8Anqurns9dVVQH1Cre3OclEkonp6elBy5MkHcZAAZDkNcz88r+mqr7eun/6wtRO+9zf+vcBa2Z9fXXr+y1VtbWqxqtqfGxsbJDyJEkvYZC7gAJsB/ZU1ednrdoJbGrtTcCNs/o/2u4GOgN4etZUkSRpyAZ5G+iZwEeAe5Lc1fo+DXwWuD7JBcAjwIfaupuA9wKTwDPA+QPsW5I0oHkHQFV9Hzjci9g3HGJ8AVvmuz9J0sLySWBJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKRODfI6aEkjsGz50cz8OY5+nLT6ZPbtfXjUZSw5BoB0hHn+wLNceMXjoy5jqLZtWTnqEpYkp4AkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnfJBMEmvej79vDgMAEmvej79vDiGPgWU5KwkP0oymeSSYe9fkjRjqAGQZBlwBXA2cArw4SSnDLMGSdKMYZ8BnA5MVtVDVfUccB2wccg1SJIYfgCsAvbOWp5qfZKkIUtVDW9nyQeAs6rqwrb8EeCtVfWxWWM2A5vb4puBHw2wy+OBnw3w/SNRb8fc2/GCx9yLQY755Koam2vQsO8C2gesmbW8uvX9n6raCmxdiJ0lmaiq8YXY1pGit2Pu7XjBY+7FMI552FNAdwDrk6xLchRwHrBzyDVIkhjyGUBVHUjyMeDbwDLgqqq6b5g1SJJmDP1BsKq6CbhpSLtbkKmkI0xvx9zb8YLH3ItFP+ahXgSWJL16+DI4SerUkgyAJFcl2Z/k3lHXMgxJ1iS5Ncn9Se5LctGoa1psSV6b5AdJ/qsd89+PuqZhSbIsyQ+TfHPUtQxDkoeT3JPkriQTo65nGJKsSHJDkgeS7EnytkXZz1KcAkryduCXwNVV9ZZR17PYkpwInFhVdyb5XWA3cG5V3T/i0hZNZl4NeUxV/TLJa4DvAxdV1W0jLm3RJbkYGAd+r6rOGXU9iy3Jw8B4VXXzHECSHcB/VtW2dsfk66vqqYXez5I8A6iq7wFPjLqOYamqR6vqztb+BbCHJf6Edc34ZVt8TftZev+bOUiS1cD7gG2jrkWLI8mxwNuB7QBV9dxi/PKHJRoAPUuyFjgVuH20lSy+NhVyF7Af2FVVS/6YgS8AnwR+M+pChqiA7yTZ3d4UsNStA6aBL7epvm1JjlmMHRkAS0iSNwBfAz5RVT8fdT2Lraqer6o/ZuaJ8tOTLOnpviTnAPuraveoaxmyP6mq05h5i/CWNsW7lC0HTgOurKpTgV8Bi/LqfANgiWjz4F8Drqmqr4+6nmFqp8e3AmeNupZFdibw/jYnfh3wziT/OtqSFl9V7Wuf+4FvMPNW4aVsCpiadUZ7AzOBsOAMgCWgXRDdDuypqs+Pup5hSDKWZEVrvw54F/DAaKtaXFV1aVWtrqq1zLxG5Zaq+vMRl7WokhzTbmygTYO8G1jSd/dV1WPA3iRvbl0bgEW5oWNJ/knIJNcC7wCOTzIFXFZV20db1aI6E/gIcE+bEwf4dHvqeqk6EdjR/sjQ7wDXV1UXt0V25gTgG+3vAS8H/q2qvjXakobi48A17Q6gh4DzF2MnS/I2UEnS3JwCkqROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXqfwHyCzSH9H24HwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2183 planets in 1622 OBSERVED systems\n"
     ]
    }
   ],
   "source": [
    "npl = np.zeros_like(obs['target_id'], dtype='int')\n",
    "\n",
    "for i, sid in enumerate(obs['target_id']):\n",
    "    npl[i] = np.sum(obs['target_id'] == sid)\n",
    "\n",
    "obs['npl'] = npl\n",
    "    \n",
    "plt.figure()\n",
    "plt.hist(npl, bins=np.arange(1,npl.max()+1), color='cornflowerblue', edgecolor='k')\n",
    "plt.show()\n",
    "\n",
    "print('%d planets in %d OBSERVED systems' %(len(obs['target_id']), len(np.unique(obs['target_id']))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Build the underlying physical catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# READ IN PHYSICAL PLANETS\n",
    "phys_file = MAINPATH +'catalogues/syssim_physical_planets_raw.csv'\n",
    "\n",
    "phys_data = []\n",
    "with open(phys_file) as infile:\n",
    "    reader = csv.reader(infile)\n",
    "    \n",
    "    for row in reader:\n",
    "        phys_data.append(row)\n",
    "        \n",
    "KEYS = phys_data[20]\n",
    "VALUES = phys_data[21:]\n",
    "\n",
    "def getkey(keyname,keys=KEYS,values=VALUES, convert_to_float=True):\n",
    "    '''\n",
    "    keyname = (string) of column definition\n",
    "    '''\n",
    "    kid = keys.index(keyname)\n",
    "    \n",
    "    dummy = []\n",
    "    for row in values:\n",
    "        dummy.append(row[kid])\n",
    "        \n",
    "    if convert_to_float:\n",
    "        return np.array(dummy, dtype='float')\n",
    "    else:   \n",
    "        return np.array(dummy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "phys = {}\n",
    "\n",
    "phys['target_id'] = getkey('target_id')\n",
    "phys['star_id'] = getkey('star_id')\n",
    "\n",
    "phys['planet_mass'] = getkey('planet_mass')\n",
    "phys['planet_radius'] = getkey('planet_radius')\n",
    "\n",
    "phys['clusterid'] = getkey('clusterid')\n",
    "phys['period'] = getkey('period')\n",
    "phys['ecc'] = getkey('ecc')\n",
    "\n",
    "phys['star_mass'] = getkey('star_mass')\n",
    "phys['star_radius'] = getkey('star_radius')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only keep (twin) systems with at least one planet detected\n",
    "keep = np.isin(phys['star_id'], obs['star_id'])\n",
    "\n",
    "for k in phys.keys():\n",
    "    phys[k] = phys[k][keep]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IMPORTANT NOTE\n",
    "\n",
    "* \"target_id\" is a unique identifier for a system in the physical catalogue\n",
    "\n",
    "* \"star_id\" may be shared between systems because stars are reused during simulation (for efficiency)\n",
    "\n",
    "Due to an earlier confusion between these terms, I mistakenly marked all systems for which \"star_id\" had at least one observed planet. This led to high multiplicities (>20) because I was combining multiple systems - all with the same star_id - into a single system. So, the \"physical_systems\" catalogue I sent to Andrew contained ~12000 planets, ~3200 of which were around stars which did not have detected planets but which had stellar twins which did have detected planets.\n",
    "\n",
    "So there are two outputs from this notebook:\n",
    "\n",
    "1. \"syssim_physical_planets_including_twins\" contains the (falsely included) twin systems; ~12000 objects\n",
    "2. \"syssim_physical_planets\" contains only the real systems; ~8800 objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEWNJREFUeJzt3X+s3Xddx/Hny64MBeOGu3ajLbuLKZBqdJCbMcQYdLJfGgsJwc4KDcFU46ZDSczgn+EPDH8oKhFnJqsUrZsLzKwhjbMMEkIisA7nWDfJrrDaW/vj6mSgJGMbb/84n+qhtL0/eu459H6ej+TkfL/v7+d7Pp9v0txXv5/v93tOqgpJUn++a9IDkCRNhgEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnFgyAJM9P8rkk/5zkQJLfbvXLknw2yWySv03yvFY/v63Ptu3TQ5/1zlb/YpJrVuqgJEkLy0JPAicJ8IKq+u8ka4FPAzcDvwncU1V3Jflz4J+r6rYkvwr8SFX9SpKtwBuq6ueTbAbuBK4AXgx8HHhpVT13ur4vuuiimp6eHsFhSlI/Hnzwwf+oqqmF2p23UIMaJMR/t9W17VXATwG/0Oq7gHcDtwFb2jLAR4A/bSGyBbirqp4GvpxklkEY/OPp+p6enmb//v0LDVGSNCTJwcW0W9Q1gCRrkjwEHAf2Af8KfKWqnm1N5oD1bXk9cAigbX8K+P7h+in2kSSN2aICoKqeq6rLgQ0M/tf+8pUaUJIdSfYn2T8/P79S3UhS95Z0F1BVfQX4JPBq4IIkJ6aQNgCH2/JhYCNA2/59wH8O10+xz3Aft1fVTFXNTE0tOIUlSVqmxdwFNJXkgrb83cDrgMcYBMEbW7PtwL1teU9bp23/RLuOsAfY2u4SugzYBHxuVAciSVqaBS8CA5cAu5KsYRAYd1fVx5I8CtyV5PeAfwLuaO3vAP6qXeR9EtgKUFUHktwNPAo8C9x4pjuAJEkra8HbQCdpZmamvAtIkpYmyYNVNbNQO58ElqROGQCS1CkDQJI6ZQBIyzS98WKSjP01vfHiSR+6VonF3AUk6RQOzh2jdo+/32w7Nv5OtSp5BiBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQ/CKNz2vTGizk45w+kSMthAOicNqlf5QLItsn0K42KU0CS1CkDQJI6ZQBIUqcMAEnq1IIBkGRjkk8meTTJgSQ3t/q7kxxO8lB7XT+0zzuTzCb5YpJrhurXttpskltW5pAkSYuxmLuAngXeUVWfT/K9wINJ9rVtf1RVfzDcOMlmYCvwQ8CLgY8neWnb/AHgdcAc8ECSPVX16CgORJK0NAsGQFUdAY605a8leQxYf4ZdtgB3VdXTwJeTzAJXtG2zVfUlgCR3tbYGgCRNwJKuASSZBl4BfLaVbkrycJKdSS5stfXAoaHd5lrtdHVJ0gQsOgCSvBD4KPD2qvoqcBvwg8DlDM4Q/nAUA0qyI8n+JPvn5+dH8ZGSpFNYVAAkWcvgj//uqroHoKqOVdVzVfVN4C/4/2mew8DGod03tNrp6t+iqm6vqpmqmpmamlrq8UiSFmkxdwEFuAN4rKreN1S/ZKjZG4BH2vIeYGuS85NcBmwCPgc8AGxKclmS5zG4ULxnNIchSVqqxdwF9BrgzcAXkjzUau8CbkhyOVDAE8AvA1TVgSR3M7i4+yxwY1U9B5DkJuA+YA2ws6oOjPBYJElLsJi7gD4N5BSb9p5hn/cA7zlFfe+Z9pMkjY9PAktSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnVowAJJsTPLJJI8mOZDk5lZ/UZJ9SR5v7xe2epK8P8lskoeTvHLos7a39o8n2b5yhyVJWshizgCeBd5RVZuBK4Ebk2wGbgHur6pNwP1tHeA6YFN77QBug0FgALcCrwKuAG49ERqSpPFbMACq6khVfb4tfw14DFgPbAF2tWa7gNe35S3Ah2vgM8AFSS4BrgH2VdWTVfVfwD7g2pEejSRp0ZZ0DSDJNPAK4LPAuqo60jYdBda15fXAoaHd5lrtdHVJ0gQsOgCSvBD4KPD2qvrq8LaqKqBGMaAkO5LsT7J/fn5+FB8pSTqFRQVAkrUM/vjvrqp7WvlYm9qhvR9v9cPAxqHdN7Ta6erfoqpur6qZqpqZmppayrFIkpZgMXcBBbgDeKyq3je0aQ9w4k6e7cC9Q/W3tLuBrgSealNF9wFXJ7mwXfy9utUkSRNw3iLavAZ4M/CFJA+12ruA9wJ3J3kbcBB4U9u2F7gemAW+DrwVoKqeTPK7wAOt3e9U1ZMjOQpJ0pItGABV9Wkgp9l81SnaF3DjaT5rJ7BzKQOU9K3OXwuDE/Pxu3TDOp44dHQifWv0FnMGIOk7yNPPQO2eTN/ZdmwyHWtF+FUQktQpA0CSOmUASFKnDABJ6pQXgTUS0xsv5uCcFwilc4kBoJE4OHdsInemZNv4+5RWC6eAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUqQUDIMnOJMeTPDJUe3eSw0keaq/rh7a9M8lski8muWaofm2rzSa5ZfSHIklaisWcAXwIuPYU9T+qqsvbay9Aks3AVuCH2j5/lmRNkjXAB4DrgM3ADa2tJGlCzluoQVV9Ksn0Ij9vC3BXVT0NfDnJLHBF2zZbVV8CSHJXa/vokkcsSRqJs7kGcFOSh9sU0YWtth44NNRmrtVOV/82SXYk2Z9k//z8/FkMT5J0JssNgNuAHwQuB44AfziqAVXV7VU1U1UzU1NTo/pYSdJJFpwCOpWqOnZiOclfAB9rq4eBjUNNN7QaZ6hLkiZgWWcASS4ZWn0DcOIOoT3A1iTnJ7kM2AR8DngA2JTksiTPY3CheM/yhy1JOlsLngEkuRN4LXBRkjngVuC1SS4HCngC+GWAqjqQ5G4GF3efBW6squfa59wE3AesAXZW1YGRH40kadEWcxfQDaco33GG9u8B3nOK+l5g75JGJ0laMT4JLEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjq1rCeBJfXp/LWQZOz9XrphHU8cOjr2flc7A0DSoj39DNTu8febbccWbqQlcwpIkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6tSCAZBkZ5LjSR4Zqr0oyb4kj7f3C1s9Sd6fZDbJw0leObTP9tb+8STbV+ZwJEmLtZgzgA8B155UuwW4v6o2Afe3dYDrgE3ttQO4DQaBAdwKvAq4Arj1RGhIkiZjwQCoqk8BT55U3gLsasu7gNcP1T9cA58BLkhyCXANsK+qnqyq/wL28e2hIkkao+VeA1hXVUfa8lFgXVteDxwaajfXaqerf5skO5LsT7J/fn5+mcOTJC3krC8CV1UBNYKxnPi826tqpqpmpqamRvWxkqSTLDcAjrWpHdr78VY/DGwcareh1U5XlyRNyHIDYA9w4k6e7cC9Q/W3tLuBrgSealNF9wFXJ7mwXfy9utUkSRNy3kINktwJvBa4KMkcg7t53gvcneRtwEHgTa35XuB6YBb4OvBWgKp6MsnvAg+0dr9TVSdfWJYkjdGCAVBVN5xm01WnaFvAjaf5nJ3AziWNTpK0YnwSWJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6tSC3wWkc8v0xos5OHds0sOQdA4wAFaZg3PHqN3j7zfbxt+npLPjFJAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ3yOQBJ3/HOXwtJJtL3pRvW8cShoxPpe6UZAJK+4z39DBN5wBEg21bvk/VOAUlSpwwASeqUASBJnTIAJKlTZxUASZ5I8oUkDyXZ32ovSrIvyePt/cJWT5L3J5lN8nCSV47iACRJyzOKM4CfrKrLq2qmrd8C3F9Vm4D72zrAdcCm9toB3DaCviVJy7QSU0BbgF1teRfw+qH6h2vgM8AFSS5Zgf4lSYtwtgFQwD8keTDJjlZbV1VH2vJRYF1bXg8cGtp3rtUkSRNwtg+C/XhVHU7yA8C+JP8yvLGqKkkt5QNbkOwAeMlLXnKWw5Mknc5ZnQFU1eH2fhz4O+AK4NiJqZ32frw1PwxsHNp9Q6ud/Jm3V9VMVc1MTU2dzfAkSWew7ABI8oIk33tiGbgaeATYA2xvzbYD97blPcBb2t1AVwJPDU0VSZLG7GymgNYBf9e+oOk84G+q6u+TPADcneRtwEHgTa39XuB6YBb4OvDWs+hbknSWlh0AVfUl4EdPUf9P4KpT1Au4cbn9SZJGyyeBJalTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpU2f7baA6hemNF3Nw7tikhyFJZ2QArICDc8eo3ZPpO9sm06+kc49TQJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlT3gYqSWdw/lpoP307VpduWMcTh46uaB8GgCSdwdPPMJHnerJt5R8mdQpIkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROrernAPxhFkk6vbEHQJJrgT8B1gAfrKr3rlRfk/phFn+URdK5YKxTQEnWAB8ArgM2Azck2TzOMUiSBsZ9DeAKYLaqvlRV3wDuAraMeQySJMYfAOuBQ0Prc60mSRqzVNX4OkveCFxbVb/U1t8MvKqqbhpqswPY0VZfBnxxbAMcnYuA/5j0IMbMY+6Dx3xuuLSqphZqNO6LwIeBjUPrG1rt/1TV7cDt4xzUqCXZX1Uzkx7HOHnMffCYV5dxTwE9AGxKclmS5wFbgT1jHoMkiTGfAVTVs0luAu5jcBvozqo6MM4xSJIGxv4cQFXtBfaOu98xO6ensJbJY+6Dx7yKjPUisCTpO4ffBSRJnTIARijJxiSfTPJokgNJbp70mMYhyZok/5TkY5MeyzgkuSDJR5L8S5LHkrx60mNaaUl+o/2bfiTJnUmeP+kxjVqSnUmOJ3lkqPaiJPuSPN7eL5zkGEfNABitZ4F3VNVm4Ergxk6+6uJm4LFJD2KM/gT4+6p6OfCjrPJjT7Ie+HVgpqp+mMENHFsnO6oV8SHg2pNqtwD3V9Um4P62vmoYACNUVUeq6vNt+WsM/jCs6iedk2wAfgb44KTHMg5Jvg/4CeAOgKr6RlV9ZbKjGovzgO9Och7wPcC/T3g8I1dVnwKePKm8BdjVlncBrx/roFaYAbBCkkwDrwA+O9mRrLg/Bn4L+OakBzImlwHzwF+2aa8PJnnBpAe1kqrqMPAHwL8BR4CnquofJjuqsVlXVUfa8lFg3SQHM2oGwApI8kLgo8Dbq+qrkx7PSknys8Dxqnpw0mMZo/OAVwK3VdUrgP9hlU0LnKzNe29hEH4vBl6Q5BcnO6rxq8Etk6vqtkkDYMSSrGXwx393Vd0z6fGssNcAP5fkCQbf7PpTSf56skNacXPAXFWdOLP7CINAWM1+GvhyVc1X1TPAPcCPTXhM43IsySUA7f34hMczUgbACCUJg7nhx6rqfZMez0qrqndW1YaqmmZwUfATVbWq/2dYVUeBQ0le1kpXAY9OcEjj8G/AlUm+p/0bv4pVfuF7yB5ge1veDtw7wbGMnAEwWq8B3szgf8IPtdf1kx6URu7XgN1JHgYuB35/wuNZUe1s5yPA54EvMPi7seqejk1yJ/CPwMuSzCV5G/Be4HVJHmdwJrRiv2A4CT4JLEmd8gxAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1Kn/Be3xlBZkGFF0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "npl = np.zeros_like(phys['target_id'], dtype='int')\n",
    "\n",
    "for i, sid in enumerate(phys['target_id']):\n",
    "    npl[i] = np.sum(phys['target_id'] == sid)\n",
    "\n",
    "phys['npl'] = npl\n",
    "    \n",
    "plt.figure()\n",
    "plt.hist(npl, bins=np.arange(1,npl.max()+1), color='orange', edgecolor='k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "12609 planets in 2475 (TWIN) PHYSICAL systems\n"
     ]
    }
   ],
   "source": [
    "print('%d planets in %d (TWIN) PHYSICAL systems' %(len(phys['target_id']), len(np.unique(phys['target_id']))))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/syssim_physical_planets_including_twins.csv'\n",
    "\n",
    "    with open(filepath, \"w\") as outfile:\n",
    "        writer = csv.writer(outfile)\n",
    "        writer.writerow(phys.keys())\n",
    "        writer.writerows(zip(*phys.values()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in the masses from Andrew"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeMAAAEyCAYAAADJFbiyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3WmMJOl95/dv5BF5Z+VRV9bdXUeze3ruJoccckSBNAmBWmgBUbIsYA1atsA31pq+YO9KMPTKkCEIhveF9gXha4FdGDDkBWwIlrHGSoJkgRySzSFneqZrpqq7rq7K+468IiIz/KIZoazqquo6su7/ByDYnZUVGVnVmF8+z/N//o9iWRZCCCGEuDiui74BIYQQ4qaTMBZCCCEumISxEEIIccEkjIUQQogLJmEshBBCXDAJYyGEEOKCSRgLIYQQF0zCWAghhLhgEsZCCCHEBfOc54sNDw9bc3Nz5/mSQgghxIV5+PBhwbKskZc971zDeG5ujp/85Cfn+ZJCCCHEhVEUZeMoz5NpaiGEEOKCSRgLIYQQF0zCWAghhLhgEsZCCCHEBZMwFkIIIS6YhLEQQghxwSSMhRBCiAsmYSyEEEJcMAljIYQQN46u62iahq7rF30rwDl34BJCCCEumq7rbGxsYFkWiqIwOzuLqqrouo6u66iqiqqq53pPLx0ZK4ryPyuKklMU5VHfYwlFUf5fRVFWfvH/8bO9TSGEEGIwdF3HsizC4TCWZTkhvLGxwc7ODhsbG+c+Yj7KNPX/CvzKnsf+CfBvLctaBP7tL/4uhBBCnLvjTjmrqoqiKGiahqIozqh4b0Cfp5dOU1uW9TeKosztefgfAr/8iz//C+Cvgf96gPclhBDiDF3klOwgHTTlfBhVVZmdnX3h/e8N6PN00jXjMcuy0r/4cwYYO+iJiqJ8F/guwMzMzAlfTgghxKCcJMAuq/4RrT06Psp72fsh5KCAPi+nrqa2LMsCrEO+/n3Lsh5YlvVgZOSlRzoKIYQYgMOmbi96SnaQ9ptyPs21wuHwhXwwOenIOKsoSsqyrLSiKCkgN8ibEkIIcXIvG/kOMsAu2kWPaAflpCPj/wv4zi/+/B3g/xzM7QghhDitl4187QCbmJi40lPUtosc0Q7KS0fGiqL8bzwv1hpWFOUZ8IfAfwf874qi/EfABvDvnuVNCiGEOLqjjHyv0ijyuhSbHeYo1dS/fcCXvj7gexFCCDEAV2nq9mVBe52KzQ4jHbiEEOIauuwhDEcL2pNWS1810ptaCCHEqZy0z/NRqrqvU7HZYWRkLIQQ4sROM4181LXtqzLlfhoSxkIIIU7sNNPIe4MWQNO0fRtyXNcQtkkYCyGE2NdRqphPO41sX/umFGodRMJYCCHEC44ajoOaRr4phVoHkQIuIYQQLzhOy8xBNN24KYVaB5GRsRBCXCMnaZCx3/fsDUd4vp5rG/Q67k0p1DqIhLEQQlwTJ1l3Peh7+sMRIJ1O0+l0SKfTpFIpfD7fwNd1Dwvh696FS6aphRDimjjJaUyHfY89/QxgWRY+n49er4fP5zvX057sDww7OztsbGxc6VOmDiJhLIQQ18RJ1l2PutdXURQ6nQ4ul4tOp3Ou67rX6cjHg8g0tRBCXBMnWXc9yvf0P2dubs557LzC+CYUd0kYCyHENXKSkNzve/au0Z5229Jp1ntvQnGXhLEQQohd7DXaTqeDaZosLCw4a8cnvdZpm3lc1xC2SRgLIcQ1dtKtTp1Oh0qlQqPRwDAM5ufnX7qXeL/XuunNPI5KwlgIIa6Bg4LwJKNSVVUxTZNGo4GqquRyOTweD+Fw+MBrHLZF6rqv9w6ChLEQQuxx1fa0HhSEJx2VqqrKwsKCc23DMIjH47t+Lvvdw36vdRPWewdBwlgIIfpchQML9n5YOCwITzoqDYfD3Lt3D03TSKfT6Lp+6DUOey0J4ZeTMBZCiD6XfY3zoA8L+wXhaUelqqqSSCQIh8NHOr1JRsAnJ2EshBB9Lvsap6ZpaJq2a9rYXsvdLwgHEYzH2bN82X5eV4WEsRBC9LnMIzxd10mn0+TzefL5PJOTk7tGwZfpXsXxSBgLIcQeJwm2QRR9HXQN+3Fd1/F6vSwtLVEul0mlUhLA14SEsRBCnNIgir4Oq4i2HzcMA8MwnC1HextxHPaB4KpViN80EsZCCHFKgyj6OugadgMOn8+HruuYponH8+J/ug/rmnUVKsRvOgljIYQ4pUEUfamqimEYZLPZF0av6XSaXq9Hu91mZmaGsbGxF0J/b9csgHv37p1qv7E4PxLGQghxSmdd9JVKpfD5fE7Y7w19ewq61WrRaDQIhUJ4PJ6B7DcW50PCWAghBuC0IWwXZ8Xj8Rcad/h8Pmdkm0qldr1e//S0YRhEIhGi0Sg+n29g+43F2ZMwFkKIS8AevZZKJUzT3PX4YUG6d3o6kUgwMTHxwqEOEsKXm+uib0AIIa46Xded0exJqapKKpVyCrTsFpT21w46Man/UIdQKEQgEJDgvYJkZCyEECdkh3A6ncbr9Q6kUjkYDB6r0Kr/UAePx7NrelpcHRLGQghxApqmsbq6iq7rVCoVlpaWDj3V6ChOWmhlH+oga8JXl4SxEMIhjSGORtd1VldXnW1IzWaT7e1tYrHYqX5u9lS1pmkHTksf9r3yO7u6JIyFEIA0hjgOXdfxeDyEQiEqlQrdbhev17vv847z4cbuPW1ZFvV6XX4HN4iEsRACuPxHBw7SaWcA7O1GsVgMl8tFIBB4oRHHST7c3KTfgdhNwlgIAVz+owMH5bCQPGpI9283gucdsvZrxHHcYL0pvwPxIgljIQRwcxpDHNYDuv9AhlQqdei6bf/PaL+f20mC9ab8DsSLJIyFEI6bEAAHhaQd0qqqsr6+jq7rhMPhI00v7/dzO2mw3oTfgXiRhLEQ4so6ydrvQSFph3S5XAYgHo8PZKuSBKs4CgljIcSVdJrq78NGspqmOdPWsm4rzouEsRDiVC5qb/JZVB6rqkoikSAcDsu6rThXEsZCiBO7yL3Jg6o83u/DhISwOG+nCmNFUf4z4HcBC/gI+B3LstqDuDEhxOV3kftiB1F5LI1OxGVx4lObFEWZBP4T4IFlWfcBN/DvDerGhBCX30Xviz3sNKN+B52q1P9hwrKsU526JMRpnHaa2gMEFEUxgCCwc/pbEkJcFeexL/Zla9L7fX1vFfTGxgadTgfTNFlYWCAcDjv3b3+YMAzj1NXTQpzUicPYsqxtRVH+BNgEWsC/sSzr3wzszoQQV8JZrq/uN41sP94ftHu/3v9YMpmk0+lQqVRoNBoA3Lt3z7lvu4I6nU5TKBQoFosyXS3O3WmmqePAPwRuARNASFGUf7TP876rKMpPFEX5ST6fP/mdCiFunL3TyJqmsbGxwc7ODhsbG2ia9sI0897vATBNk0ajQSgUwuPx7JqOtkPZ6/WiqiqapqFp2kW9ZXFDnWaa+t8B1izLygMoivKvgXeBf9n/JMuyvg98H+DBgwfWKV5PCHEBLvJYxb1r0sCugjFg3zXr/sfC4TALCwsAeDwefD7fC9PZAIZhsL6+7rzucY8wFOI0ThPGm8AXFUUJ8nya+uvATwZyV0KIS+Giq433rkkDFIvFXUFrj2b7w7P/EAf7e+/du7frQ8Xe95ZMJtF1fSCdt4Q4rtOsGb+vKMqfAT8FTOADfjECFkJcDxd9pN9+o/K94bzf+b8HrSfbhVv7vTd7NCydt8RFOFU1tWVZfwj84YDuRQhxyRy2deksp6/trUh20PZXQfe/Xv+a8d4PCy/7ILH3vYXDYem8JS6MdOASQhzooK1LZzl9rWkaq6ur6LpOPp8nGAw6U852FXT//R30YeFle6APOzBCiPMmYSzECVxkUdN52+89nnT6+ih7hldXV3n27Bkul4t6vY5hGAwPDztV0HsD96B9zkfZA30Tfn/iapAwFuKYLrqoaVBO84HiJJ23jvJzs0O+Xq9Tq9Xw+XykUilisRg+nw/AWd89Sh9pCVtxVUgYC3FMF13UNAin/UBxks5be39u/aHaH6yKouD3++l0OgwNDREMBpmYmEBVVWcN+bQfgm7SzIa4GiSMhTimi+7HPAiD+EBx3CDb23oynU7j9Xp3BauqqiwsLNBsNvF6vcRiMQKBgPM6g/gQdF1mNsT1ImEsxDGdRz/ms7ZfM429079n8Zr2z03XdQqFwr7BGg6Hee2111hdXX2hScegjky86jMb4vqRMBbiBK5qCNv6gxEY2PTvUV7XbrjR37xj75YpwOmaddAe45Pe43WY2RDXj4SxEDeUHWiH7dU9y9fe2yXLtrKy4tzD4uLiCxXSdpifdCR/HWY2xPUjYSzEDTeokeJxi6L265IViUTY3t4mEAjQarVIpVIkEokXXue0a74SwuKyOfGpTUKIq8ceUe49tSiVSuH3+4lEIvs+57Dvtx/vP01pv+896H72nrp0Ft8jxGUnI2MhboiDRpT249vb2xiGgaIozMzM4PP5do06DztbuD8gS6USpVKJRCKx7z7i/tHz3lF5IpFgcnISXddJJBK7eknbXjaSl21L4iqSMBbihjioitgOr0AgAECj0cDn8zmjzoN6PWuaRrFYxLIsDMMAoFQqkU6nAXYd3GB//34fBvau3y4uLr60a9ZBa76ybUlcVRLGQtwQB40o7UBrtVoYhoHP56PT6ezaUrTf98Pufb/Dw8POlHEikTjywQ37tbF8WYAe1k5Tti2Jq0jCWIgb4qARpT0atc/z7d/vuzck9ztbuFQqYZqmc81MJkOpVHppmJ9FSMq2JXFVSRgLcYP0j1L7/w7Pp5U7nQ7pdJpUKvXCNLP9/P6/p1IppznHxsYGAB6PB9M0mZub2zfMNU070/cn25bEVSRhLMQ1cpRTkQ4q4rIsC5/PR6/X23fN+CB2l6xyuQxAPB4/9Pn2OnOxWDyTNV0JYXEVSRgLcU0c51Sk/dZt7endTqeDpmkvTFMf9JpPnjzB7XbT7XYBKJfLuFwu5ubmXvhwIGu6QuxPwliIa+KwoLNDEV7s72x/LZlM0mg0mJ6eRlEUUqnUoUGp6zqbm5u43W5UVSUUCuH1eonH406gZzIZZ+RsfziQNV0hXiRhLMQ1cVDQ7R0xp1Ip5/m6rjtrvqZp4vF4CIfDlMtlZ233sClvj8dDLBaj0WigqirBYNB5na2tLUqlEqFQiFgshq7rhMNhWdMVYh8SxuLSkyYOL2f/jPqD1u47nc1m0TSNsbExJ2DD4bATxNlsllAoRCgUotVqsbOz4xxxuF/zD5uqqvh8PmKxGOFwmIWFhV0j7Z2dHUKhEI1GY9eUt/wehXiRhLG41KSJw8H6p57tU5cMw3CmlzVN4wc/+AGdTod8Ps8bb7yxKxTtkW1/YE5PT+PxeFBVlbW1NUzTdJ67XxgftFVK1/V9g1oIsT8JY3GpScHP3+ufIYC/P2Ch2Ww608vr6+vOdLBdGT0xMQFANBplenoaeN4pS9d1FEV5ITDr9TqappHP5wHw+XxMTEzse0rSSbpkCSFeJGEsLjUp+HlO0zRnbdfn85FMJp0PKbquY5rmrq1Fdgi6XC52dnbodrvO448fP2Z7exuAyclJ5ufnd42YZ2dnKZVKAM6HoM3NTYLB4LFmJySEhTg6CWNxqckIi33Xdn0+H4ZhoGkaPp/P2UZkTxErikI4HObOnTusra0RjUZJp9NUq1WePHlCOp1mbGwMgPn5+ReqrsPhsHMqkqIozsj7ps9OCHFWJIzFpXdTQvigQjV7qt7lcpHP5ykWi3i9XhRFYXh4eNeo1h4pw/N15HK5TLlcJplMsrOzw+rqKpVKhWfPnhGPx/F6vbteZ2/VtX1PxWLxxs9OCHGWJIyFuAReVqhWKBTodDp0Oh1u377tbD8C9l3DtZt3NBoNcrkcm5ubNBoNHj9+TDwexzRNxsbGmJubc44p3G993u6WBbwQ/EKIwZEwFmLATrIVyw5CVVWdPb6JRML5eiqVwufzoWkapmny2WefAbtnDfb+2TRNdF1namqKQqFAJBKhWCwSj8cJBAIsLi7umqJ+2alMN2WGQoiLIGEsxACddCuWqqoYhsH6+rrzd3sUau/ntYMxEong8XiIx+NOYdfe4ipVVVlYWABwRrbNZpOhoSHC4TBzc3NMT08f6VQmmZ4W4uxJGF8AaWJx9Rz1d3acrVh7typFIhFisRijo6MvfC2ZTAI4U8r1et2poj6ouMrerqRpmrP+++abb+L1ekkkEs61+u19fze9eE6I8yJhfM6kicXVc5zf2VG3YvVf0zAM4PkItlqtEgwGnbOA9762PVq2QxKe7zfOZrP7Fn6l0+ldxyL6fD7Gx8eP/G9OQliI8yFhfM6kicXVc5zf2VG2Yum6TqlUotPpkEgkyGazAM5Wo+HhYVRVdRpz2CcoaZrmrCX3n4IEYBiGs0bcv02p0+lgmiadTseZ6u5f/5V/e0JcDhLG5+ymN7E4qyn6s5z6P+7vrP8e+ltW2vpHq/bzAWfPsKqq/PCHP6RarZLJZBgfH6fVajE2NrZrLdm+vn3ucKPRYHV1lXv37jlf73+teDxOKBQinU47W6NkZkaIy0HCeECOGgY3uYnFWU3Rn/XU/0l/Z/Z99U8Tt1otLMtidHSUVCrF8PCwUzVtH+JQLpfZ2tqi0WiwtrZGMBgkGo0yOTnpjJTt77GrphuNBqFQCI/Hs2t0nEqlUBTF6TFtGIY08BDiEpIwHoDjhsFNC2HbWU3Rn8fU/0l+Z/Z9+Xw+Op0OrVaL7e1t50Sk2dnZXUFsj1gLhQLNZpNKpUKz2aRaraLrOuvr6wSDwRcqre2qabtVZv9WJXs7VDAYZHJy0tkadVNnZoS4rCSMB0DWgY/mrKbo97vuZahYt7crVSoV0uk07XabUqnkjHDtdpabm5vouk6lUmFpaYlAIMDo6CjtdptUKsWdO3cIBoOYpsno6KgzLd1fNX3v3r19328ymSQSiTg/E7t1pn1/8u9UiMtBwngAbvo68FGd1RT9fvtjL1PFumEYRCIRZmdnMQyDUqlEIpHAsiyWl5ep1+tOcJfLZRRF4d69e7zyyiusrq4yNDREJpOhVCqxs7NDKpVyAtXWX9BlT3fbxyraP4P+5wkhLhcJ4wG4yevAx3VWP5/+62qadilmKuwZE8uyaDQarK+v43K5qFQqtFotWq0WMzMzqKpKpVJhbGyMVCoF4PSfnp2dpVAo0Ov1iEQiJBIJhoeHD3y9vccqJhIJJ5z321cshLgcJIwHREL48riImYr9psXt4ipd17l9+zaKomAYhnOYQzgcxjRNms0mLpcLr9dLJpOh0WjQarWc0W+326Xb7dJsNolEIoRCoQO3TO09VlFma4S4GiSMxbVzmpmKk/aV3nvakW18fBzDMAgEApimycbGBs+ePSOXy3H//n1CoRAjIyOMjo6Sy+VYXl5me3ubRqPBxsYGn/vc54jFYgDMzMxw584dZ5/xfu/b/hAia8NCXC0SxuJaOmn1c39XrFQqtWuUe1g3LXtEWiqVePz4MZZlOduZVFVlYmICeN5lyzRNDMMgFAoRCARQFMU5g9g0TUzTpNvtOkcdLi0tOa0tD5tqluUSIa4uCWMh2N0VKxwOs76+TqPRoFqtOm0kDyoEs0ekpVKJTCZDPp/H4/Hw6NEjXn31VTweD+Pj4yQSCRRFcVpS9no9FEVxtibZgZ3P5+l0OoyOjnLr1i1SqdSuYw5l25wQ14+Esbjx7BGxvc3ILpAKhUJOdXN/K8q9VFUllUrx+PFjCoUCKysrdLtdCoUCy8vLjIyM8Omnn3Lr1i3Gx8cBnCntmZkZ5xqqqvL6668zNjbG2toa0WiUcDhMOBx+oTJaAleI6+VUYawoSgz4H4H7gAX8h5Zl/WAQNybESR133dfu4dxoNOj1etRqNYaHh51tQisrKy8029j7OoBTvRyPx6nX6ywsLOByuUgmk3z22WeUSiXC4TBvv/32rnaVe7cfJRKJXR22ZB+7ENffaUfG/wz4fyzL+g1FUVQgOIB7EuLETtIa0656zufzlMtl4vG4c5LS+Pg4hULBKYayg9A+R7h/hGuaJqurq7TbbXq9Hnfu3MEwDILBINVqlYmJCSqVijPC3huymqZRLBb3vXfZxy7E9XbiMFYUZQj4JeA/ALAsSwf0w75HiLO23yjSfry/MUb/31VVZWZmhuXlZZrNJt1uF7/fTzKZZHJykkajQT6fR1EUZ8S8urrK2toamUwGv9/P5uYmlUqFarVKs9kkmUxSr9f56le/Cjxv/JHP52m322SzWdrtthP4dsgC+46ApTBLiOvvNCPjW0Ae+F8URXkdeAh8z7KsRv+TFEX5LvBdwFkfE+Ks7N1jDLyw7Wi/9VdVVbl9+zalUon19XXq9Trz8/MEg0HC4TBPnz4lFApRqVS4desW+XyelZUVPvroI3w+H7lcDrfbTbfbJRgMcvv2bUKhEPB8dP3FL36R5eVlLMuiXC4zMjIC/P1xiXbAFovFfUfAEsJCXG+nCWMP8Bbwjy3Lel9RlH8G/BPgv+l/kmVZ3we+D/DgwQPrFK8nLonL0Pf5IHtHkftNBR80cu7fchSLxRgeHiYajbK1tUUmk2F0dJSdnR12dnZ48uQJmqbR6/XodrtUq1WSySSxWAzLev7P3OfzOVXQqqo6wfvZZ59RLped4qz+n6GMgIW4mU4Txs+AZ5Zlvf+Lv/8Zz8NYXGP2Wql9QtBlrOzdG2T2tiPTNBkeHsYwDLLZrPOclZUVZxtTNBrF4/E4pxs1Gg0++ugjNjc3WVlZIZVKEY1GCQaD9Ho9p8grGo3SbrdJJpN8+ctf5itf+Qrj4+O7wtjeTzw5OelsV9r7s5MQFuJmOnEYW5aVURRlS1GUO5ZlfQp8HfhkcLcmLht7rTSbzTqjx8te2WtvO7I/QKTTaarVKo1Gg1gshqZprK+vUywWyWQytFotvF4vpmmysrLCs2fPePToEaOjozQaDbxeL8Vi0fnz4uIipmkyOztLtVrl9ddf59VXX32hOYes+wohDnPaaup/DPyrX1RSPwV+5/S3JC4rXdfxeDyEQiEajca+I7vLyl77XVtb4+HDh4RCIVZWVlAUhVqtxrNnz9jZ2UFRFCKRCIqi8LOf/YxkMkmn0yEUCuH1eolGowQCAZrNJolEgrt372IYBiMjI8RiMRYXFw/8mUgICyEOcqowtizrZ8CDAd2LuORU9flh9bFYzGnPeFbhMsh16f6irlarRbVapdfrsba2RjKZpFqtOqck2VPTPp8Pl8tFr9cjFovh8XiYmZmh0+ngcrnweDzMzc3h8XgYHR1lfn7+hQ8nl3ltXQhxuUgHLnFk5zXVepK9woex79teB04mk06x1vb2Nt1ul0AggN/vxzRNer0eHo/HGSW/8sorDA8PMzk5yaeffkoulyMUCtHpdIhEIty9e/eFaelBvwchxPUmYSyO5binGZ305KTTdJw66HXT6TSNRoN4PE6328Xj8fDs2TN6vR6dTge3241pmk5P6Gq1it/vJxwOY1kWY2NjTqDb0/Ver/elxxlK1ywhxMtIGIszcZqR4WHnEb8s4Pd7XV3XefLkCR9++CGRSIRcLofP52NzcxPDMCiXy0xMTBCJRPB6vfz85z+n0Wigqir37t3j1q1bWJaFqqrOGrHf7ycWixEIBPYN2os4U1kIcXVJGIszcZqR4UHT4UcJeLvPtKIoVCoVDMNgZWWF1dVVPvroI1KpFIVCgZGREVqtFrFYzAnWRqNBIBCg3W4725zsa9rrwaqq8tprr+3a3nXQSU5SPS2EOCoJY3EmTjsy3BvC9v+OEvBPnjxhfX2dWq1GLBYjnU47o95isUiz2aRSqeD1emm1Wvj9ft544w00TcPn8wEwNTWFaZqkUilmZ2d3FWeFw2Hu3bv30qCVEBZCHJWEsTgTgxoZ9o+G9/Zy3hvWAKVSiVqtRrfbpd1u8+GHH7K+vk6z2cTr9TI6OoplWU5QLywsMDk5ydLSEnfu3KHRaPDBBx+g6zq9Xo/79+8feGyiBK0QYlAkjMWxHKco67QhvN9oeG8v51KpRDqdRtd1tra2CAQCrK6uoigKa2tr1Ot1RkdH2dzcdLprdTodZ5vS7Owsi4uLznat8fFxQqEQjx49IhqNUq/XSSQSErxCiDMlYSyO7KhFWafdX3vU0XA6naZcLrO5uQlAJpMhFos5py55vV7cbje5XI56vY6qqvj9ftxuN+FwmEgk4vSTttd+dV0nk8nQ7Xad9WephBZCnDUJY3FkR1mzPe3+Wl3XKZVKdDodEomEMxree2373F/DMFhfX+fJkycEg0HnrGFd12m32/h8PtrtNiMjIwQCAcrlMsPDw6ysrHD//n0sy2J4eJhwOOx8iLiqXcaEEFeXhLE4sqMUZZ2mitoO206nQzqdBnBe79mzZ1QqFba3twmHw+RyOVqtFm63m62tLT777DM8Hg8ul4uhoSEikQimaaIoCqZp4vP5GBoaYmJigvn5eQqFAnfv3sXtdjujbDvcFUU5ly5jQghhkzAWR3aUoqzTVFHb08J2cHa7XfL5PD/60Y94+vQp1WrVOTfY4/EwPj5ONpul2WximibVahXLspxK6WAw6DTtiMViTE1NOfcUCASoVqsMDQ2hqirtdnvfdWkJYiHEeZAwFsfysoA6bRV1Op12RsaBQICVlRUAZyTs8/nIZDLkcjnW19fxeDzE43EqlQrlchm32008HmdsbAy/3+806wgEArzzzjvcvXvXuZ5lWYRCIcLhMPV63fkAIVPTQojzJmEsBm6/EN6vqGu/x1KpFK1Wi1wuh6ZpTgOPZrPprBfbLSvj8TjVahWXy0UgEKDX69FsNimVSoyOjpJMJhkfH2d+fh5N07h3757z516v54yEAWnQIYS4UBLG4sz1rwWbpumsw/ZXTCeTSXRdp1arsbm5yccff8zKygqdTod4PI7P5yMajeJ2uymXy3g8HqfDVq/Xc0a4AC6Xi9u3bzM1NUUgECAajTI+Ps7Y2Biw/1S6hLAQ4iJJGIuB2zvi1TSNcrlMo9FA13UMwyAej1Mul4nFYqyvr/Phhx9SKBTwer0UCgWnqKrb7ZLL5QDwer3E43Hm5ub4pV/PXsk/AAAgAElEQVT6JQqFAoqi0O12WVtbo1Qq4fP58Pv9zM7Ocv/+fUKhEGNjY0xPTzthLa0qhRCXjYSxGChd11lZWXGCbnZ2lnQ6zc7ODsVikenpaba3t9nY2KBWqxEMBqnX61QqFac4yw5t+9zgarVKIBDA5/MRDAaZmJhgYWGBkZERdnZ2KBQKjI+PMzQ0xOzsLH6/n6mpKXq9ntNJa7+DHCSEhRCXhYTxNTGIg+wHcQ1N09je3iYQCDh9n3VdZ3FxEXi+ValerxOJRBgbGyOXy1EoFHj69Cn1et0pnopEIkxPT5PL5XC5XE6TjlgsxsTEBNlsli996UvMzc2xvr7O06dP+fnPf06v12NxcZG5uTna7TbxeHzX+xJCiMtIwvgaGMRB9oO4xl6GYZDNZimXy07RlWEY5PN5PvzwQ1wuF9FolEQiQafToVKpsL6+zvz8PKZp8vrrr7O6uorP52NyctJ5fq/X46OPPiIajTI9PU273abX6zE1NYXX6yUSiRCPx6nX6+i6LkcYCiEuPQnja2AQB9kftbvWy0bO4XCYyclJdF3H7/c7vaAzmQyZTAa/30+n0yGXy7G8vMzQ0BDlcpler4fH4+HJkyf0ej16vR6WZdFutwGIRCLMzs5SKBQolUqUy2Wy2Sx3797l1q1bGIZBIBAA4M6dOyQSCRKJhKwLCyGuBAnjczaIqeC9BnGQ/d5rwPMpZ/s+XzZy7n9fi4uLTqBvbW1RqVR49uyZs0WpXq87/zNNk3w+71y32+3S7Xbp9Xpks1k6nQ69Xo9SqcQ777zj7Amenp4mFotRLpeJx+N861vfcu63/2AHCWEhxFUgYXyOzmIqGAZTHdx/DcBpD2nf52Ej573vK5VKkU6naTQa1Go1NE2jUCgwOjqKx/P8n1yhUKBSqVCpVDBNk263i2VZWJZFuVxmaGjIKbxKJBJEo1FarRaxWAyPx8Pdu3cJhULOCDgcDu971KEQQlwFEsbnaFDTyfuF7iBG2v1bkfa7z4NG3/b7UlXVmT5eWVkhn8+zsbFBt9ulWCzSarXo9XosLy+TyWRotVpYluVcx96WlEwmeffddymXy/j9ftrtttPUQ9d1hoaGME0TgHa7TTqdHtgHGyGEuAgSxufIDrRSqeSEyXGc1cj6oPvc2xTjoNG3qqrO6UkAtVqNx48fs729TTabZXp6GlVVaTabuFwuSqUSbrebbrfrjJR7vR5+v9+pona73fR6Pebm5qhWq9y/f9+Z2p6eniYQCJz6g40QQlwWEsbnSFVVUqkUq6ureDyeY4/oBjGyPup97he8hx0OkUql0HWdeDzO48ePefToEdlslmw2S6FQIJFIOAczZDIZp62ly+Vy9hOHQiEmJyeZmppidHQUn8/H7du3dwU6QDQa3bWuLdXSQoirTsL4AgSDwRMF6iAKtQ6zdwr8sOvbz+2/N/tow5/+9KfkcjmazSa6rlMoFADIZDJUq1U8Hg+RSASv10s4HMbj8eDz+ZidnSUajTIyMkIymcTn82GaJrOzs4RCIe7evesUk/WvV0u1tBDiqpMwPmf2lG42mz12iJxlG8fDpsB1XXcOVAiHw+i6zurqKpZlkU6nnWMI4/E4P/rRj/joo494+vQpuq7TaDTo9XpkMhna7TYjIyNOU49oNEo0GiUWixGLxRgbG3Ne0+VyMTIywtTUFHNzc07w7zdFLoQQV52E8RVzVqPAg6bA7faW29vbTk9pwzAol8sYhsEnn3yCx+PB6/XSbrf56U9/ys7ODoZhYBgGjUYDRVEwTRO3243b7XaONFxcXMQwDF555RXu3LlDIBBga2uLZ8+eUavVeOeddxgaGpLGHUKIa0/C+Jzpuu4ceDDoiurTOGgKXNM08vk8zWaTWq1GLpcjm83icrlYXl5mZWUFv9+PpmlOSNuFVpZl0ev1cLvdu5qA2NPRQ0NDaJrGu+++SyAQoFgsOiPkoaEhPve5zzkjcZmKFkJcZxLG5+y0674vm04+aXDtNwWu6zqbm5t88sknbG1tAc+nqVdWVpwe1OVy2WnkYQdvtVql1+s511YUhXA4zPz8PEtLSywuLtJsNhkeHiYUCjnT2UNDQzSbTVKpFOFw2OlTLSEshLjuJIzP2d7mGv1roUdx2HTyUTtkHfRae7+maRqVSgW/38/4+DiPHz8ml8tRqVRoNBq0220KhQL1eh2Xy0Wv18Pn8zn3YP9/KBRiZmaG+fl551xhwzCYn5/HMAzGxsbodrskEgmCwSDDw8O7umgJIcR1J2F8AeyQscPTMAxnNPiyADpoZH2cDlkHbafae7rR5uYmz549Y319nUAgQK1Wc/pMV6tVCoUCrVZr1weCQCCAoijOyNjtdjM+Ps7ExATf+MY3mJiYcI5FHBkZIRQKMT09TTqdRtM0fD6fBLEQ4saRML4g/V2r1tfX0XWdcDj80n3Hh+0BflmHrP6gth/vn5K2zyFWFIVIJML29jZ+v59isUi326VQKNBsNsnlcui67mxdAvB4PKiqyvz8PJlMBsMwsCyLaDTKW2+9xf3793n99dcplUrOsYrj4+OMj4+faZW4EEJcBRLGF8QOz3K5DHCsc3f3Cyw70OwtSHu/1mw22dnZQVVVhoeHSafTdDodTNNkfHyccrnMysoKgUCAJ0+eMDc3x8bGBjs7O+zs7GBZFp1Oh3w+74yIu90uqqpimiaKouD1elEUhXg87lRYJ5NJvvKVrzA7O8vo6CilUgnA2WP8soYiQghxE0gYX5D+8LRDeRCBVCwWsSyLYrHojLLtaer+k5Oi0Si6rpPP5/nhD39IKBTi0aNHeL1eWq0WgUCAUqnEysoK9Xqd7e1tGo0GpVIJy7Jwu92YpulsWQoGgySTSd577z0AQqEQoVCIe/fucf/+fWcPsX3GsH24gxBCCAnjC2V3rfJ6vei6jmEYaJq2a8R4nArpg9aNNU2j1+sxMjLC+vo6a2trBAIBQqEQjUbDmX5eWVmhUCjg9Xr55JNP6PV6PH36FNM0qVarTgADdLtdgsGg05ZSVVVnNDw9PY1lWczMzLCwsLCrleXi4qJMRwshxB4SxgPystA86Ov2vuNQKMRnn33mhOns7CzAsQ6G2LtuDFAqldB1nV6vx8bGBtvb2wwPDztbkIrFohO89hpws9l0WliWy2V0XafdbjuvE4lE8Hg8xONxp/mH2+1mcnKSt99+mzfffJNKpbJrTbj/HiWEhRBiNwnjATjKtqKDvn7Y2jFw6MEQe9eYdV0nlUo5X7fD1y6mmp6exjRNRkdH6Xa7xGIxp8hqZ2eHQqFANpvF7XZTKBQwDGNX/2lbu90mHo/TaDTodrv4/X5UVWV8fJzp6WnnfewNYiGEEPuTMB6Al52mdNjX+9eO7fVdOwT3q5DuD+p0Ou1sjQKcAqpkMkm5XCaXyzn3YFkWgUCAcDhMJpMhHA6ztbXFo0ePeP/99/F6vVQqFQA6nY5zdvB+3G63c61areZMTd+6dYuxsTFnfViCWAghjkbCeABe1lXrKF+3C5o0TSOdTlMoFFAUxRnp9u9N7nQ61Go1AoEAY2NjZLNZ4PloNJvN8vjxY7rdLh9//DHxeNwJxnQ6jaIoPHr0CMMwyOfzpNNpTNOk0WhQLBadEO6flj7oPcdiMVRVxe/3s7S0hMvlcqqkhRBCHJ2E8QC8bJ/sYduO9j5PVVW8Xi+qqjrrtXZ/Zl3X6XQ6VCoVKpUKvV7PeS4875jVarVwuVz4/X68Xi+3b992KqPts4V3dnbI5/Nks1kymQy9Xg/DMGg2m87o215ztixr1z0Gg0GmpqZ44403SKVSxONxEokE4+PjEsRCCHFCEsbnqH/bUf+Id2+Bk2EYrK+vA+wKRUVRnJOQYrEYoVCIaDS6q/p6eHiYv/zLv+T99993Rs+jo6O4XC6azSblcpm1tTWePn3qdM/yeDwvtOe0LItIJILf78cwDLrdLqZpMjY2xmuvvUYymSQUCjE0NMTbb79NIBBwqsOFEEIcz6nDWFEUN/ATYNuyrH9w+lu6eo7SbrJ/3bhUKrG6uupsDdpb0JVKpdB1nXg8Ti6XI51Ok0gkqNVqTE5O4vf7CYVCzr7ddru96zqhUIhAIMDt27fRNI1MJoPH46FUKpHJZEin09TrdWetudvt4vF4nG1LNlVVefDgAYqiUKlUaLfbPHjwgKWlJTweD6Ojo/j9fqamppwWlrJOLIQQxzeIkfH3gMdAdADXupL2Bm2pVHqhv3L/urFpmng8ngOfb59YZE8720Va29vbTkMOr9dLNBql1WrtKgyzp7K73S4//OEP2dnZodfr0Wg0aDabGIZBqVRygti2N4jheeOOcDjsHGcYCASYnZ3lC1/4Auvr67jdbqeXtIyIhRDi5E4VxoqiTAG/Cvy3wH8+kDu6guygLZVKpNNpAOr1+gsj3v7TmuyOWPZ+Xnvq2p5yTqVSrK6uEo1GaTQauFwuEokEkUgEwzDY2dlxDlyA52Fqt5pstVrs7OywtbVFqVRidHSUcrnsFGu1Wq1D348dstFoFFVVmZiYYHx8HMMwGBkZwePx8N577znr2RLEQghxOqcdGf8PwH8FRAZwL1fK3iYes7OzThgmEol9tzj1H8oAOFuYWq0Wm5ubNBoN4vG40/AjGAw6QWfv5V1eXqbVaqGqKnfv3gWeh6d9wlKj0eDp06c8e/aMbrfLzs4Oq6urzpGHh/F4PPR6PWdEHI1G8fv9xONx3n33XTqdDqOjo877SiQSA/+5CiHETXTiMFYU5R8AOcuyHiqK8suHPO+7wHcBZmZmTvpyl4qmaayuruLxePD5fM4I2N7DWyqV8Pl8B66f2l23RkZG+OSTT3j69CntdptIJEIwGNy1x7hUKlEoFEilUni9Xm7duoXP52NnZ4eVlRWCwSCtVotOp0O73eZHP/oRjx49IpPJ0Gw2abVa1Gq1F6al+9kj4XA47LyfsbExlpaWuHv3Lnfu3GF6eppisehUWsvasBBCDM5pRsZfBn5NUZRvAX4gqijKv7Qs6x/1P8myrO8D3wd48OCB9eJlrhZd11ldXSWbzRIKhYjFYruacNiFUHNzcwcGVv/6cSwWY3R01Dkf2OVyMTEx4bTE7B9tZ7NZisWi8xqpVIpYLMbOzg7b29v84Ac/4OHDh5RKJRRFcUbvhwUxwNDQECMjIySTSe7cuUMikeCVV14hEokwPz/vTEXbW6ykUEsIIQbrxGFsWdY/Bf4pwC9Gxv/l3iC+juytQPYhC/Yar13EZU9Rv0wymSQSiaAoCtvb27hcLrLZLIFAgIcPH/KlL33JCT07uFutFu12m0AgADxfZ37//ffRNI3NzU3W19dpt9sYhoFpmui6jsvlOvQ+fD4fU1NT3L59m6997WssLCywtLTkNAuB3duvJISFEGLwZJ/xMamqis/nIxaLEQ6HWVhYcALKMAyy2eyuUel+W5z6t0HNzMw4+3iz2SyqqlKpVNja2qLT6eDxeFAUBZfLRSaTYXl5GcAp/FpZWXH2D9dqNbLZLL1ez3k9l8uFy+Xa9ZgtHo8TCoW4f/8+ExMTTE5O8s4778hasBBCnLOBhLFlWX8N/PUgrnXZvazblmEYpNNpdF0nFAo5xVj9xVv926DK5TIejwfLsvj000959OgRo6OjNBoNhoeHnZaTjx49Ynt7m7W1NSqVCqurq1SrVWq1GpZlOccw7g3d/ULY5XIRjUYZGhpidnaW6elpp1q7WCzuaiIihBDi7MnI+AT2C+H+oxALhQKBQIB8Pk+z2WRoaGhXUw7DMJxtTalUiqdPn7K6uuocY2iaJltbW/j9fuB5gVWn0wGeH4lYLBYplUq0Wi1nm9LLKqVt9tGHd+7c4fbt2ywuLjI2NkYsFmNubu7AEb0QQoizI2F8Sv2nKNmHL3Q6HQzDcNpNKopCp9PZFXL2sYbxeJxkMsny8jJut5tWq0UulyOfzwMwMjKCaZp8/PHHFItFdnZ2KJfL1Ot1J6CPKhQKkUwmeeutt/jmN7/JvXv3uHXrFqqqSqW0EEJcIAnjU+jf4qQoCs1mE13X6Xa7lEol58+apuFyuZibm3P2H4+MjNBoNCiXy4RCIebn52m1Ws6JTeFwmGw2y89+9jNqtRr1ep1arUaj0aBSqRw7iIPBIK+99hpvv/023/rWt3jzzTd3df1KJBJSKS2EEBdEwvgl9jb36H/c3uKkqiqmaVKr1QiHw/R6PaeXczqd5r333iMQCKBpGltbWzx58gS/38/c3Bzj4+PA85FyvV6nXq+TyWTIZrPk83kqlYozyi4UChiGceQpaXjeLCSVSvG1r32N3/7t32ZxcZHx8fF9T5aSEBZCiIshYXyI/Q6AsB+3tzgpisInn3yCz+ej1Wpx69Yt5zzgcrnsFGm99957tFotPvvsM1wuFzs7O0xNTfHpp58617ErtQuFAul0mkaj4fSUtvtTH5W9L3hubo7XX3+d3/iN3+C9996TwBVCiEtIwngfdtj2Vz5rmoamac4xiPZoNZfL0W63UVUVl8tFpVJheHiYjz/+mA8//JBYLIbL5SIcDtNoNKjX61iWxc7ODn/2Z3+G2+3G7XYTCoWIx+O4XC5M06TT6VCtVo+0Z7lfKBRibGyM119/nWg0ype//GVisRj379+XIBZCiEtKwniP/tGwPRLVNG3XucL2tqRms4nf7ycajdLr9bh16xaKojhT18ViEa/XS6PRYHl5mbGxMUqlEmtra5imSbVaZWhoiI2NDUzTpNvtOvuMTxLEU1NTDA8Pc+/ePe7fv8/09DRLS0uEQiHZOyyEEJeYhDG714X37gOORqO7umzZRxAWi0U6nQ5erxefz4fb7WZ4eJi1tTXW19d58uQJHo+HZrNJLBaj0WhgmiaZTIaPPvoIl8tFLpfD5XLRaDSwLMspyurfyvQyqqri9XqZnp7my1/+MkNDQ4yPj/OFL3yBV1991VkLllGxEEJcXjc+jPeuC6dSqX2PQwTwer00m02q1Sr5fJ7l5WUajQaRSIR3332XSCTiVEQ/efKEQqHgBOHGxgY/+clP+Oijj5zuWvB8Wrndbjvbo47K4/EQCASIRqMkEgm+/vWv841vfIOJiQl0XWdhYUFGw0IIcUVIGO9ZFwZ2HdAQDofZ3t7G6/USDAZ5+PChs+7bbDbZ2toikUjw8OFDstksjx8/JpfL0Ww2MQyDXq9HvV7n2bNnZDIZqtWq87hlWdRqtWPfczAYJBaLMTMzw8jICHfv3uVb3/oWw8PDTuMROWNYCCGujhsfxv0HMdgVzfZZvcVikc8++8zZTrSyssLDhw9pt9tOv+hKpUKpVGJlZcWphC6VSmiaRrfbpdPp8Omnn7K2tka73T5R+Ho8HtxuNyMjI/j9fiYnJxkeHubu3bu88cYbfP7zn3e2SMleYSGEuHokjPf0mgachh1+v9854vDhw4esr6/TarXIZrP4/X5M02RlZQV4PpXt8/mc6eNer0en03FC3u12H3sq2uv1OgViHo+HL3zhC7zyyiuMjo4SDoeZn5/ntdde2zUKlhAWQoir58aHMew+xGFlZYX19XWy2SyRSIRGo8HKygp/93d/x89//nNWVlZot9v4/X48Hg+tVguXy7WrLWatVsPn89Hr9ZyDGrrd7rHuyefzEQwGSSQSjIyMkEql+Pa3v81Xv/pVJ3DlQAchhLgeJIz72IHq9Xpxu93s7Ozws5/9jGq1ys7ODs1mk0ajga7rNBqNA6/T7XZpNpsnvo9AIMDY2Bh+v5+pqSkWFhZ49913uXXrltPMQwghxPUhYdzHHiHb4VsoFFhfX0fTNJ48eYJpmsdqRXkSdpevQCBAPB5nenqaaDRKvV53zjAWQghxvUgY91FVlWQySbVadZp3bG1tUa/XqVQqWJZ1Zq8dDoeJxWIsLS3R6XQYHh5meHiY6elppqameOONN/B45NclhBDX0Y3/r/ve83t//OMf81d/9VeUSiUymQylUunMgtjlchEIBAgGg9y+fZulpSVmZ2ep1Wrcvn0bRVFYXFwklUoRDAbleEMhhLimbnQY72344fP5ePjwIZ9++im5XI5arUatVjuzEbG9d3l0dJQ7d+4QCoWYmZkhHo+TSqVQVZW7d+86xWWyZUkIIa6nGx/G/Q0/Go0GbrfbWTc+yZ7go/J4PMRiMRKJhDMl/corr/DVr36VsbExYPexhhLCQghxfd3oMFZVFcMwnIMdUqkUo6OjZ1Ko5fV6URQFy7IIBoMEAgESiQQzMzPMzMzwuc99jldffZXp6WkJXiGEuGFudBjbDMPANE0AZ9uQvT/4tFwuFz6fD5/Px9DQEIqiMDExwdjYGOPj44yNjTE/P8+dO3ecKWkhhBA3y40OY3uaWtd1MpkMH3zwAX/zN3/D9va2E86n1ev18Pl8hEIhRkdHmZycZGZmhsnJSVRV5c6dO7z55pskEgkJYiGEuKFudBirquoca/iDH/yAtbU1Hj16RKVSGcj13W43Ho8Hr9dLt9ulVCoxMzODruv4fD6mpqZ455135HQlIYS44W5sGNtbmmZmZnj69CnFYpF6ve6cuDQI3W6XbrdLr9cjHA4zOjrK4uIiS0tLvPXWW8zPz0s3LSGEEDcrjPv7R6fTaWdL0/T0NPV6nXQ6TbVaHdjrBYNBQqEQsVjMOS7RNE2i0SjT09MSxEIIIYAbFMb9e4qbzSYej4dEIkGpVEJRFJaWlnj69OlAX9MwDLxeL6Ojo7hcLlRVpV6vEwwGZX1YCCGEw3XRN3Be+vcUezweTNOkVCqRTqfRNI1nz55RKpVO9Rr2+jA8X48eGhpiamrKOV1peHiYcDhMIBAgnU4f+0hFIYQQ19ONGRmrqoqiKGiahs/nY25uDk3TANja2qJYLGIYxqlew+PxYBgGqqrS6/UwDIN8Ps+DBw+wLItXX30Vt9vNyMiIU8UtI2QhhBA3KoxnZ2d3tZVUVZWtrS3+/M//nKdPn1Iul0/1Gp1OB4/H4+wtjsfjRCIRgsEgfr+fyclJPB6Ps1YtQSyEEAJuUBgDL/R27h8tBwKBU3fdsntYK4qCx+NBVVVisRiKohAOhwkGg7z11lu7PgwIIYQQNyqM92OPWgep2+0SCASYn5/nm9/8Jvl8nkAg4EyLSxW1EEKIfjemgOsgMzMzfOlLXxpY4w27tWa73SYWi+Hz+TAMA7/fT6lUkqItIYQQL7jxYRwOh/m1X/s17t27d6rrKIqy6+9+v59arcaTJ08IBAKYpkkkEpGpaSGEEC+48dPUAIlEgtHR0RN/v8vlwu12Y5omiqLQ6/VoNBo8fvwYXdcZHR0lFAoxNjYmYSyEEOIFEsbs7sx1Er1eD5fL5YSyZVlEIhG8Xq9zTGMikeD27dsDvGshhBDXxY2fprZPbDIM44Wp5uMwTROv10symSQajTpT0u12G8uyGBkZcZp/CCGEEP1udBjbLTK3trZYW1tztiadVLvdRlEUJ3gBCoUCLpeLUqlEMpmUMBZCCPGCGx/GnU6HVqvFzs7OQK7Z6/UYGhpyOnEFg0FisRj5fJ5sNivV1EIIIV5w49eM0+k0y8vL1Ov1gVyvWq1Sq9VotVq43W46nQ5DQ0PE43FKpRIej4fZ2VkZIQshhHDc+DBOpVKYpsnw8DAul4ter3eq69mjbL/fTygUIhgM0uv18Pl8hEIhOp2O9KQWQgixy42eplZVFZ/Px/j4OPPz88Tj8YFct16v0263qdVquN1uAoEAlUqF1dVVNjc30TRNpquFEEI4bnwYz87OkkwmmZ+fJxQKDeS6dm9qt9vN+Pg4ExMTpFIpJiYmsCyLdDrNxsaGBLIQQgjgFGGsKMq0oih/pSjKJ4qifKwoyvcGeWPn6dmzZ3S7XYaGhgZyvW63i6ZpKIpCu93m1q1bzM/PO+Ebj8edIxSFEEKI04yMTeC/sCzrHvBF4D9WFOV0PSXPma7rlEolLMsim82Sz+cHdu1Op0O323XWoScmJggGg84pUXKEohBCCNuJC7gsy0oD6V/8ua4oymNgEvhkQPd2puw9xpqm8ezZM+dAh0FqNBp0u128Xi+WZTE1NUUwGGR4eJhEIiFhLIQQAhhQNbWiKHPAm8D7+3ztu8B34fkJSRfNbn1p7zFuNBo0Gg0ePXqE2+0e6GvZQZxIJAiFQmiahs/nkyAWQgixy6nDWFGUMPB/AP+pZVm1vV+3LOv7wPcBHjx4cLoWV3vYoaqq6pHCzR4N240+qtUq6XSazc1NWq2W01/6tNub4PnhEVNTU7z33nt87WtfY3x8/Fj3KoQQ4uY4VRgriuLleRD/K8uy/vVgbulo7GC1LAtFUUilUgCHhp2u62iaRjqdplgsks1mabfbPHr0iGKxiGEYuN3ugYTx9PQ0v/zLv8yv/uqvMj4+LiEshBDiQCcOY+X5qQr/E/DYsqz/fnC3dDS6rmNZFuFwmFKpxOrqqlMgdViHq0ePHvHBBx+g6zrBYNDpjNXr9eh2u6c6LAIgFovx+c9/nm9/+9t8/etfZ2ZmRkJYCCHEoU4zMv4y8O8DHymK8rNfPPb7lmX936e/rZdTVdWpTDZNE4/HQzgcdhpq7BeAmqaxvb3N2toavV4Py7KIxWK43W68Xi+KouBynazAPBwO8/bbb/PgwQN+5Vd+hS9+8YvOYRFCCCHEYU5TTf3/AacbRp6C3bDDXjfe3Nwkm80eOLLVdZ2nT5+SyWRQVZVer0ev10PTNBqNBhsbG9RqNbrd7rHvZXJykj/4gz9gcXGRpaUlZ1paCCGEOIor3ZvaDrx0Ou2E7e3bt0mn07umqnVd5+c//zl/8Rd/wePHj9na2sLlcuHz+ej1etRqNWcb0nF9/etf50/+5E9YWFiQdWEhhBAncqXDGJ5PPZfLZfL5PJqmUSgUAEgmkyQSCec56XQaRVEYGhpifX2dWq2Gruu0Wi263S6dTudYrzs8PMzv//7v853vfMd5HSGEEOIkrnQY2yy8kSkAAAi5SURBVNPTT58+JZvNoqoqKysrZLNZisUir7zyColEgs3NTUqlEisrK3z88cfkcjlqtdqJqqZdLhe/+Zu/yfe+9z3efvttGQkLIYQ4tSsdxpqmkcvliEQiLC8vMzY2RjabpVQq8fDhQx4+fMi9e/coFou0Wi0CgQDlcvlEQez1evn1X/91fvd3f1eKs4QQQgzUlQ5jAMMwMAyDeDxOIpFga2uLer1OPp+nVqvxwQcfYBgGtVqNTCZDOp0+9mu88cYb/Omf/ikPHjyQkbAQQoiBu9JhbG9vajab5PN58vk8pVKJbDZLNpvF7XaztraGaZpUq1Xa7faxrn/nzh1+53d+h+985zuMj4+f0bsQQghx013pMIbn/a4nJiZot9sUCgU8Hg+KolAsFllfX2d7e/vY1xwbG+Of//N/zrvvvit9pIUQQpy5Kx3G9shY13W8Xi8rKytUKhUURaFerzuV1cfxR3/0R/ze7/2erAkLIYQ4N1c6jAGazSbVahVVVfF4PDSbTZaXl8nlcrRarSNf57d+67f44z/+40txspQQQoib5UqHcSaT4cc//jGhUIjV1VVWVlbI5XLkcjk0TTvSNXw+H3/7t3/L5z//+TO+WyGEEGJ/J2vE/P+3d28xVl11HMe//4EhDp3EziiXCqOtkWiIidQQqdEYU3lAbawEoyUSitH4orEajVRfwAdJTIyXJsbQtJI+GC/BJjaSaEgtiU8gyIN20NgUh3KpHRjKRWAO4N+Hs1FqoEyZc85m7fl+XmbW3mfm/GfNmvObfdbee90CWq0WBw8eZHx8nCNHjvx3fvj48eNTDuL169dz+vRpg1iSVKtij4xbrRYDAwP09fUxOjrK0aNHuXjxIidPnpzS12/bto0NGzZ0t0hJkqag2DC+cvLWrFmzuHDhAidOnGBsbIxLly7d8Gt37tzJypUre1ClJEk3VnQYj4yMMDQ0RGZy6tSpKQXx/v37WbZsWQ8qlCRpaoqdMwYYHh5myZIlZOaUric+ceKEQSxJuuUUe2R8xdy5czl27NgNj4onJye9eYck6ZZUdBi3Wi1GR0fZu3fvqz4uM3tUkSRJr13Rb1MD7Nix41WPisfGxnpYjSRJr13RYTwxMcGePXuuu3/Tpk3eUUuSdMsrOoz37dvH0aNHr7t/8+bNvStGkqSbVGwYt1otdu3add39k5OTvStGkqRpKDqMH3nkkWvu6+/v98xpSVIxig3jiYmJ6+5rtVo9rESSpOkpNoy3bt16ze19fcX+SJKkGarY5NqyZcs1t1++fLnHlUiSND3FhrEkSU1RZBhfb064v7+/x5VIkjR9RYbx2bNnr7ndE7ckSSUqMoxHRkbqLkGSpI4pMozPnTtXdwmSJHVMkWEsSVKTNCaMz5w5U3cJkiTdlMaE8eDgYN0lSJJ0UxoTxpIklarIMF6zZs0r2qtXr66pEkmSpq/IMN64cSOLFi1iYGCABQsWsG7dOq8xliQVa3bdBdyMoaEh1q5dy/nz5+nr62NwcJBWq+WyiZKkIhUZxgsXLmTFihUcOnSIgYEB5s+fbxBLkopVZBgPDg6yatUqJiYmmDNnDsPDw4axJKlYRYYxtAPZy5kkSU1Q5AlckiQ1ybTCOCJWRcTfIuK5iHi4U0VJkjST3HQYR8Qs4EfAh4GlwNqIWNqpwiRJmimmc2T8HuC5zHw+M1vAz4H7O1OWJEkzx3TCeBHwwlXtw9U2SZL0GnT9BK6I+HxE7I2IvePj491+OkmSijOdMD4CjFzVXlxte4XMfDQzl2fm8nnz5k3j6SRJaqbphPEfgSURcVdEzAEeAJ7qTFmSJM0cN33Tj8y8FBFfBH4HzAJ+kpnPdqwySZJmiMjM3j1ZxDgw1sFv+UbgeAe/n6bGfq+PfV8f+74+Jff9WzLzhnO0PQ3jTouIvZm5vO46Zhr7vT72fX3s+/rMhL73dpiSJNXMMJYkqWalh/GjdRcwQ9nv9bHv62Pf16fxfV/0nLEkSU1Q+pGxJEnFM4wlSapZkWHsOsq9ExEjEfFMRIxGxLMR8VC1fTgidkbE36uPQ3XX2lQRMSsi9kfEb6r2XRGxuxr/v6jugKcOi4jbI2J7RPw1Ig5ExHsd970REV+pXm/+EhE/i4jXNX3cFxfGrqPcc5eAr2bmUuAe4AtVfz8MPJ2ZS4Cnq7a64yHgwFXt7wDfz8y3ASeBz9ZSVfP9EPhtZr4DeBft34HjvssiYhHwJWB5Zr6T9h0eH6Dh4764MMZ1lHsqM49l5p+qz8/QfkFaRLvPn6ge9gTw8XoqbLaIWAx8FHisagdwL7C9eoh93wUR8XrgA8DjAJnZysyXcdz3ymxgICJmA3OBYzR83JcYxq6jXJOIuBO4G9gNLMjMY9WuF4EFNZXVdD8Avg78u2q/AXg5My9Vbcd/d9wFjAPbqimCxyLiNhz3XZeZR4DvAodoh/ApYB8NH/clhrFqEBGDwK+AL2fm6av3Zfv6OK+R67CIuA94KTP31V3LDDQbeDfw48y8G/gX//eWtOO+O6p5+Ptp/0P0JuA2YFWtRfVAiWE8pXWU1TkR0U87iH+amU9Wm/8ZEXdU++8AXqqrvgZ7H/CxiPgH7emYe2nPY95evX0Hjv9uOQwczszdVXs77XB23HffSuBgZo5n5kXgSdp/C40e9yWGseso91A1R/k4cCAzv3fVrqeAB6vPHwR+3evami4zv5GZizPzTtrj/PeZ+WngGeAT1cPs+y7IzBeBFyLi7dWmDwGjOO574RBwT0TMrV5/rvR9o8d9kXfgioiP0J5Lu7KO8rdrLqmxIuL9wB+AP/O/ectv0p43/iXwZtrLYn4yMydqKXIGiIgPAl/LzPsi4q20j5SHgf3AusycrLO+JoqIZbRPnJsDPA98hvYBjOO+yyLiW8CnaF/NsR/4HO054saO+yLDWJKkJinxbWpJkhrFMJYkqWaGsSRJNTOMJUmqmWEsSVLNDGNJkmpmGEuSVLP/ALwH1kcL4fBSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# read in mass\n",
    "phys['radius'] = np.load(MAINPATH + 'catalogues/syssim_physical_pred_masses_dec17.npy')[0]\n",
    "phys['mass']   = np.load(MAINPATH + 'catalogues/syssim_physical_pred_masses_dec17.npy')[3]\n",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.plot(phys['mass'], phys['radius'], 'k.', alpha=0.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Eliminate the false twin systems\n",
    "keep = np.isin(phys['target_id'], obs['target_id'])\n",
    "\n",
    "for k in phys.keys():\n",
    "    phys[k] = phys[k][keep]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8865 planets in 1622 PHYSICAL systems\n"
     ]
    }
   ],
   "source": [
    "print('%d planets in %d PHYSICAL systems' %(len(phys['target_id']), len(np.unique(phys['target_id']))))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/syssim_physical_planets.csv'\n",
    "\n",
    "    with open(filepath, \"w\") as outfile:\n",
    "        writer = csv.writer(outfile)\n",
    "        writer.writerow(phys.keys())\n",
    "        writer.writerows(zip(*phys.values()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Identify which systems are in the observed catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pull observed/physical periods and IDs\n",
    "obs_per  = np.array(obs['period'], dtype='float')\n",
    "phys_per = np.array(phys['period'], dtype='float')\n",
    "phys_rad = np.array(phys['radius'], dtype='float')\n",
    "phys_mass = np.array(phys['mass'], dtype='float')\n",
    "\n",
    "obs_id  = np.array(obs['target_id'], dtype='int')\n",
    "phys_id = np.array(phys['target_id'], dtype='int')\n",
    "\n",
    "\n",
    "# new lists of TRUE DETECTED periods (slightly different than observed due to simulated measurement uncertainty)\n",
    "true_per = np.zeros_like(obs_per)\n",
    "true_rad = np.zeros_like(obs_per)\n",
    "obs_mass = np.zeros_like(obs_per)\n",
    "\n",
    "\n",
    "# iterate over all OBSERVED periods\n",
    "for i, per in enumerate(obs_per):\n",
    "    \n",
    "    # pull the PHYSICAL periods for this system\n",
    "    phys_per_in_this_system = phys_per[phys_id == obs_id[i]]\n",
    "    phys_rad_in_this_system = phys_rad[phys_id == obs_id[i]]\n",
    "    phys_mass_in_this_system = phys_mass[phys_id == obs_id[i]]\n",
    "        \n",
    "    # select the physical period closest to the observed period and record it as the DETECTED period\n",
    "    loc = np.argmin(np.abs(phys_per_in_this_system - obs_per[i]))\n",
    "    \n",
    "    true_per[i] = phys_per_in_this_system[loc]\n",
    "    true_rad[i] = phys_rad_in_this_system[loc]\n",
    "    obs_mass[i] = phys_mass_in_this_system[loc]\n",
    "    \n",
    "    \n",
    "obs['radius'] = true_rad\n",
    "obs['mass'] = obs_mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make DETECTED catalogue (very similar to OBSERVED)\n",
    "det = obs.copy()\n",
    "det['period'] = true_per\n",
    "\n",
    "\n",
    "# make UNDETECTED catalogue\n",
    "keep = np.isin(phys_per, true_per)\n",
    "und  = phys.copy()\n",
    "\n",
    "for k in und.keys():\n",
    "    und[k] = und[k][~keep]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12aa07c50>]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X10XOV9J/Dv797RCNEERGS1LhDHdHHJsWMaNcLpQNEOa2JeU3RCzjaBxq4xyMZvKAEbTDan2iW1iHbPoiQ2mMHYWEtO0hRaICfQADqeQOJbHPkFDCaAS4zqpFpsEZYTgkeamd/+oZnLzGhG83LvvN35fs7R8Wh0Z+7z2PJXj55XUVUQEVH9M6pdACIicgcDnYjIIxjoREQewUAnIvIIBjoRkUcw0ImIPIKBTkTkEQx0IiKPYKATEXmEr5I3mzVrls6dO7eStyQiqnv79u07oart+a6raKDPnTsXIyMjlbwlEVHdE5G3CrmOXS5ERB7BQCci8ggGOhGRR+QNdBHZISJvi8jLWb52q4ioiMwqT/GIiKhQhbTQHwJweeaTIvJxAEsAjLpcJiIiKkHeQFfV5wC8k+VL9wDYCIAnZBAR1YCS+tBF5BoAv1bVF10uDxHRNJZlob+/H5ZlVbsoNa3oeegiciqAOzHV3VLI9T0AegBgzpw5xd6OiBqcZVlYvHgxJiYm4Pf7MTw8jEAgUO1i1aRSWuj/CcA5AF4UkaMAzgawX0RmZ7tYVUOq2qmqne3teRc6ERGlCYfDmJiYQCwWw8TEBMLhcLWLVLOKbqGr6iEAf5j8PBHqnap6wsVyEREBAILBIPx+v91CDwaD1S5S0SzLQjgcRjAYLOtvF3kDXUS+DyAIYJaIHAPwd6r6YNlKRESUIhAIYHh4uCKBWA6V7DLKG+iq+uU8X5/rWmmIiLIIBAJ1F+RJ2bqMylUXrhQlIiqjZJeRaZpl7zKq6G6LRESNppJdRgx0IqIyq1SXEbtciIg8goFOROQRDHQiIo9goBMReQQDnYjIIxjoREQewUAnIvIIBjoRkUcw0ImIPIKBTkTkEQx0IiKPYKATEXkEA52IyCMY6EREHsFAJyLyCAY6EZFHMNCJiDyCgU5E5BF5A11EdojI2yLycspz/1NEfikiL4nIP4tIa3mLSURE+RTSQn8IwOUZzz0D4FOqej6A1wFscrlcRERUpLyBrqrPAXgn47mnVTWa+PRfAZxdhrIREVER3OhDvwHAUy68DxEROeAo0EXk6wCiAL43wzU9IjIiIiPHjx93cjsiIppByYEuIn8L4GoA16uq5rpOVUOq2qmqne3t7aXejoiI8vCV8iIRuRzARgD/WVV/726RiIioFIVMW/w+AAvAeSJyTERWANgC4KMAnhGRgyKyrczlJCKiPPK20FX1y1mefrAMZSEiIge4UpSIyCMY6EREHsFAJyLyCAY6EZFHMNCJiDyCgU5E5BEMdCIij2CgExF5BAOdiMgjGOhERB7BQCci8ggGOhGRRzDQiYg8goFOROQRDHQiIo9goBMReQQDnYjIIxjoREQewUAnIvIIBjoRkUfkDXQR2SEib4vIyynPfUxEnhGRNxJ/nlHeYhIRUT6FtNAfAnB5xnN3ABhW1XkAhhOfExFRFeUNdFV9DsA7GU9fA2BX4vEuAN0ul4uIiIpUah/6H6nqfyQejwH4I5fKQ0REJXI8KKqqCkBzfV1EekRkRERGjh8/7vR2RESUQ6mB/n9F5I8BIPHn27kuVNWQqnaqamd7e3uJtyMionxKDfQnACxLPF4G4HF3ikNEtcyyLPT398OyrGoXhbLw5btARL4PIAhglogcA/B3AO4G8EMRWQHgLQD/tZyFJKoGy7IQDocRDAYRCASqXZyyKKaOlmVh8eLFmJiYgN/vx/DwsGf/XupV3kBX1S/n+NJil8tCVDMaIbyKrWM4HMbExARisRgmJiYQDoc993dS77hSlCiLbOHlNcXWMRgMwu/3wzRN+P1+BIPBipSTCpe3hU7UiJLhlWy9ejG8UutomiZGR0dhWVbOVncgEMDw8LDnu6HqmUzNOqyMzs5OHRkZqdj9iJxolD70oaEh7Ny5E9Fo1LPdS/VORPapame+69hCJ8ohEAjUdbCl/kACkPWHUyAQQDgcRjQaZd+4BzDQiTwodcDT5/NBVRGLxbK2wBuhe6lRMNCJPCh1wDMejwMAVDVrC5x9497BQCfyoNRWd2YLPVsLPLN7qRHGD7yIgU7kQZmtbiB7H3o2jTAH36sY6EQeldnqLjSUuYCofnFhERGl4QKi+sUWOhGl4SBp/WKgE9E09T4Hv1Gxy4WIyCMY6ER1hnuSUy7sciGqI5xSSDNhC52ojjTCtr5UOgY6UY2ZqUuFUwppJuxyIaoh+bpUOKWQZsJAJ6ohmV0qQ0ND08K7lCmF3JulMTDQqW5lC6l6D65kl0okEoGqIhQKAQCam5vt1nqxdeRAauNgoFNdyhZSAOo+uAKBAAYHB7F69WrEYjH7+UgkYg+AFltH7s3SOBwNiorIV0XkFRF5WUS+LyKnuFUwoplkCymvzAAZHx+39zBPMgwDbW1t6OvrQyQSKaqOHEhtHCW30EXkLADrAcxX1Q9E5IcAvgTgIZfKRpRTrlN2vHDyTjAYRFNTEyYmJgAApmnia1/7Gnp7exGJRBCPx2EYRsF15EBq43Da5eID0CIikwBOBfAb50Uiyi9XSHkhuJLnfA4NDQEAli5dav/2kQzzSy+9FH19fQXXkXuzNAZR1dJfLHILgL8H8AGAp1X1+izX9ADoAYA5c+Z85q233ir5fkRuq+YgajH35sBmYxORfaramfe6UgNdRM4A8CiAvwbwLoB/BPCIqj6c6zWdnZ06MjJS0v2I3OZWSJbyQ6GUe9f7DB4qXaGB7qTL5VIAv1LV44kb/hOACwHkDHSiaksNRTdmf4RCIaxduxaxWCxtamE+pdyb3SaUj5NAHwXwFyJyKqa6XBYDYPObalZmq3hwcNDRIKplWVizZg2i0SiAD6cWFhK6uQZ1iZwoOdBV9QUReQTAfgBRAAcAhNwqGJHbMlvF4+PjRQ2iZnZ5hMPhtOmFpmkWHMyceULl4GhQtFjsQ6dqctJnPtNCpkgkAsMwsHXrVvT09JSzCtSgKtGHTlRXnLSKs/V5b9q0ia1sqikMdGoopQ4s5urz5kAl1RIGOtWNSk3by3Yf9nlTPWCgU12o1MKame7D1jjVOp5YRHWhUhtveWWDL2pMbKFTXajUvO1gMAjTNBGPx+1piJZlpe2rwlY61SoGOtWFYvqwSzkAInk9AKiq/XHo0CGsX78ekUgEALBjx468i4e4RJ+qhYFOdaOQPuxi+9ozr7/gggswOTkJAJicnMTg4KC9jW3yuZkCnZtoUTWxD508ZaY+cMuy0N/fD8uysl4fiUTw/PPPp73fa6+9Bp/vw3ZPU1PTjN097IOnamILnTwlV197rpZz6vUiknbsW9KKFSvsx/n60LlHC1UTA508JVdfe67dDVOvb2trQ29vL06ePAlVhYigubm5qIFQzlenauJeLuSKWh8ITG2hm6aJG264IWtQJ+vR1taG8fHxmq0PNZayH3BRCga6N9XiQGC2qYbJ53bu3InJyUmYpoktW7ZwQy2qedyciyrGjYMiUjlt7VuWhUsuuSTrVMNwOIzJyUnE43HE43GsWbMGCxcurPoPICI3MNCpYLmC1s2BwEJa+5nzxnP1lyel/pBJXTgEALFYDEuXLsWGDRvYUqe6x0CnguTb48StgcDMaYR9fX1pp9unlsPn80FVEYvF7BOIxsfH0dbWBp/PZ88nB4C2tja7rFu2bMGaNWsQi8Wgqjhy5AhWrlwJAAx1qm+pq+LK/fGZz3xGqT5t3rxZTdNUAGqapm7evLks99mzZ4+2tLSoYRgKQA3D0JaWFt2zZ4/u2bNHlyxZYn8NgIqI/dg0TTVNU1taWrS7u9v+mmEY08q7Z88ePffcc+3XAtBFixaVpU5ETgEY0QIylguLqCDJbhXTNMs6vzrZ2r/00kthGAbi8TgmJiYwNDSExYsX49lnn0079k1TBvVjsRhisRhOnjyJ119/HX6/H4ZhwDRNu4Weep8NGzakPbd///60RUdE9YaBTgVJBu1dd91V9lksgUAAfX19aG5utn+AAFN94alhnouq4vDhw4hGo/Ziod7e3mlh3dPTg+7u7rTXcWUn1TNOW6SalTn4mTqPXETs2SqFEBGsXLkS991337R71NqUS6JMFZmHLiKtALYD+BSm+iFvUNWcv7My0MmJUCiERx99FNdeey0WLlyIvr6+aV0wM2lqasJPf/rTGWfNMMypFlVqHvq3AfyLqn5RRPwATnX4flQFtRpomS303t5eTExMYPfu3ejo6MC8efNgGB/2GuYL9mg0mnWOPE8iIq8oOdBF5HQAXQD+FgBUdQLAxEyvodpTq10OoVDInlpomiauvvpqezpjLBbD3r17sXfvXgCAaZq46KKL8POf/zzr5lpJ3CyLvM5JC/0cAMcB7BSRPwOwD8Atqvq+KyWjinB7lacbQqEQbr75ZrvFHY1G8fjjj8Pn82UN7Fgshueeey6ttS4iuOaaazB79myMjY1h9uzZPG2IPM9JoPsA/DmAdar6goh8G8AdAL6RepGI9ADoAYA5c+Y4uB2VQ61t92pZFtasWTOt+0RV0xYKZROPxyEi9i6JGzduZIBTQ3EybfEYgGOq+kLi80cwFfBpVDWkqp2q2tne3u7gdpRPtgMc8il1OmIp9ypEOByesdskH1WFYRgYHBxkmFPDKbmFrqpjIvLvInKeqr4GYDGAw+4VjYrhpC+82EHBcva7t7W1oZiZV8nFR6lUFePj43lfW6uDwUSlcrqwaB2A74nISwA+DWCz8yJRKSp59Fk573XgwAH7sYigq6tr2jUiAgBobm7GbbfdhkWLFqGrq6uolazJH0rf+MY3sHjxYq4QJU9wNG1RVQ8CyDs3ksqvkn3hhd5rphZw8mvvvvsuDh48aM8t3759u31Nrpb6BRdcgO7ubvuEoWQ5vvvd7xZ8KEUtDgYTOcXdFj2ikkefFXKvmbplkl9LHvUGAE8//TS6uroQjUbT3uett96a9t4rVqxAT08P+vv700I5GebJ3xh49ic1Gga6h1RygUy+e83UAh4aGkoL86Rf/vKX095ndHQ07fONGzfaW9xmhnJbW1vBffs8+5O8iIFOZZEatj6fD6Ojo3Y/9YMPPpi1O+XEiRMwTTNtlkvyuuReLN/61rfsr2WGcrHdKFwhSl7DQCfXZPaZDw8PY2BgAD/60Y8QCoWwa9cuLFu2bFq3SpKq2q3v/fv34xe/+IUd6D6fDx0dHejv709rUWeGMrtRqJEx0Mmx5OHLO3bssE8PGh4eBgD8+Mc/tlvckUgEY2NjMAwj61xzwzDQ0dGBnp4eu589EonAMAx89atfTRsAzeyTT/4gYTcKNTIGOjkSCoWwevXqtIBOHkjx5ptvTmuNP/HEE/aKzsxul1gshvXr19uHNhfSnZJt8HXTpk0VqTtRrWGgU1Eyd0DMDHNgqqW9Y8cORKNRqKq9HD/zpCHDMOw55amt+KGhIbsrJV93CqcfEn2Igd5ASlkZmeuQCb/fj2XLlmXtOmltbcWJEyfsFngy1DOdd955+MpXvoK9e/fisccem7EcuWalcPohUYpCDh5164OHRFdP8vDl5CHKe/bsKfo1q1atSjsouqurK+2Q5WI/7r//fvs+fr9fRUT9fn9BZcss5+bNm4t+HVG9QIGHRLOF3iBK6ZrIfA2Q3u1x/Pjxgu4tItMGQru6uuwZLYFAAOFwuOTBTE4/JJrCQG8AlmVhdHQUPt/UP3ehy/UzuzM6OjoAwN5f/PDhw3j11VdnvLfP58ONN96Ijo4OrF+/3n6vu+++O+06hjKRcwx0D8nWR546C8Q0Tdx00005D3rINmMk2W+d3Dcl2wrPmVx99dX2wcwLFy7klEKiMmKge0SuvVNSu02AqUNGcoVp6rWRSAS9vb0488wzMXv2bIyOjiISiRQV5gBw9OhR9Pf3o62treCNs4ioNAx0j8jVR17Mzoijo6MwTROqing8bp/ZCUxNRSw2zAHg4MGDePHFF+1pis3NzTVzbimR1zDQPSJXcBezM+LJkycBZN+2NrkYqBTJ94vH45wrTlRGDPQ6lK2vfKbgLmRnxGL7xkthGAbnihOVEQO9zsy0z3ghM0WS+64AsAdHCw3YQgI/80i4rq4u7NmzB7FYjGd9EpWZ0yPoqMJS+8pPnjxph3MhLMtCMBjEtm3bsG3bNlx88cUIhUIAgI9+9KOOy9bS0oLbbrsNTU1NMAwDLS0tmD9/ftrCh0LO+iSi0rCFXmeCwSB8Ph9isRhUFTt27MDSpUsBIOuUxdTnwuEwJicn7feKxWJYtWpVyV0tyT3KOzo60mawdHd3p20XsGvXLi7NJ6oABnqdCQQCWL58Oe6//36oKmKxGIaGhtJCc3h4GIcOHcLatWsRi8Xg8/lwww034LTTTpu2YtNJv/l1111nzzHPLGNqtwq3tCWqDCn3QFiqzs5OHRkZqdj9vCqzH33ZsmV44IEHEIvF7MVD27dvT9u6NrldbbZta4uR+vqWlhZOQSSqABHZp6qd+a5z3EIXERPACIBfq+rVTt+P8suc0XLo0CF7nrjf78fY2Ni0fciTIez0B/isWbMwPj6eNgURmN7dQ0SV50aXyy0AXgVwmgvvRQVKdmtYloV169ZhcnISpmli3bp1GBwcLNt929vb8bvf/a6kg5mJqLwczXIRkbMBXAVguzvFoWINDQ3ZOyHGYjGEQiH783K45ZZbMDw8jLvuugvDw8MYHx+ftkKViKrDaQt9EMBGADnnvIlID4AeYGofEXJHcgbL2NhY2vPvvvtuWe4nItiwYUPalrdJPGCCqDaUHOgicjWAt1V1n4gEc12nqiEAIWBqULTU+9GHi4LGxsbw1FNPYXJysuTl+MUyDAOtra3Tni9kawEiqgwnLfSLAPyViFwJ4BQAp4nIw6r6N+4UrbFkW8GZfD65fe26devK2p0CAPPmzcORI0egqjBNEz6fD9FodMbWN/cyJ6oNJQe6qm4CsAkAEi302xjmpbEsC5dccgkikQgA4IEHHsC9996LhQsX2gOOhmGkLQoql6amJvh8PkxOTsIwDHznO9/htrdEdYILi2rA0NCQHebA1ODm6tWrcdNNN9kDjm7MIS/E+++/b//gmJycxIEDB7IuHiKi2uPKXi6qGuYc9NJYloUdO3ZMez4Wi2FsbAx+v9/e8KocYZ7aB+/z+exj5oio/rCFXmXhcDhtKX6q2bNnY3BwEKtWrSrLvbu7u7Fo0aK004QA4KmnnrJnrST3iSGi2sdAr7LUgymSpwVFo1EYhmFvelWOlrlpmti4cWPWfvHdu3dz1gpRHWKgV1kgEMDg4CAeffRRXHvttQBgb6rV29uLdevWuX5P0zRx77335gxrzlohqk8M9CpKTlXcvn07YrEYwuEwrrzySkSjUagqPvjgAwwMDLh6z+7u7pwtcyKqbwz0CktdHPTkk09icnLS7lKZmJjAY489VrZ7iwgWLVrEMCfyKAZ6BWXON68kEcEpp5zCpflEHsZAL5NsBzknj4+rlEWLFuHMM8/EO++8g5MnT2LFihVsnRN5GAO9DHId5Jyc0ZLaQjdNM+e0RSf8fj9WrFiB9evX2/c7ePAgFi5cyFAn8igeEl0GqQc5T0xMYGhoCP39/Th06BCuuOIKzJ07177W7TAXEftMz+TWtkmTk5Pc3pbIw9hCL4PMlnjmcXDl5PP50maxpJajqamJfehEHsZAd0nqrojj4+NYt24d7rnnnrJsqCUi9r4uPp8PZ511Fo4ePQoAiMfjCIfD9lzy3bt3Z93FkYi8h4HugmSf+cmTJ+0piOXcSGvlypVYunRp2ra6ST6fz26FJ3/IMMiJGgMD3SHLstDX15cW5oDzw5hzaWpqsgM6EAigv7/f7ocXESxfvtw+a5RnfRI1Fga6A5Zloaurq2L94wBw1VVXpQVz6l4wqZtpZQ7MJrthiMi7GOgO3HHHHRUN82xyHQGXGfQcDCXyPin3gQmpOjs7dWRkpGL3c1vqMXEdHR1YtWpV2Q+cyNTc3Izdu3cX1NrOtriJiOqPiOxT1c5817GFXqBQKITVq1eXZRFQNqtWrbK3z927dy8ef/xxe2vdQrtPuGsiUWNhoBfAsiysWbOmYmFuGAbmzJmDnp4e+/4/+clP7O6TtrY29Pf3s+VNRGkY6AUIh8OIx+Nlv49hTC3cbW5uTuvzTu0nb2trQ29vLyKRCAzDwNatW+3gJ6LGxqX/BQgGg2hubi7rPa6//nr87Gc/wze/+c2sUwwDgQA2bdqE8fFxRCIRxONxRKNRrF27FpZllbVsRFQfSg50Efm4iOwWkcMi8oqI3OJmwWrNZZddhrPPPrss720YBhYsWGCH9kzdKMFg0G7JA7APxiAictJCjwK4VVXnA/gLAGtEZL47xaodt99+Oy688EI89thjOHbsmOvvbxgGfD4fRkdHC2ppBwIBbN26FU1NTTAMY1r3DBE1LtemLYrI4wC2qOozua6pp2mLlmVhYGCgLCcIGYaBnp4edHR04MCBA9i5cyei0Wjais58Uw45JZGocVR02qKIzAXQAeAFN96v2pLL5j/44APX39swDNx33332QGZ/fz+i0Wjaik4AeZftc0oiEWVyHOgi8hEAjwLoVdX3sny9B0APAMyZM8fp7coitbULADfeeKPrYS4iuPjii3H33XfPuHQ/GAxy2T4RlcRRoItIE6bC/Huq+k/ZrlHVEIAQMNXl4uR+5ZC6iZVhGK5ud9vU1ISrrroKs2fPzrnjYa6l+1y2T0TFKjnQRUQAPAjgVVX93+4VqbIGBgbs1rgbC4fOPfdcfOELX0Bra6sdxPlmoWR2n+QKeSKimZQ8KCoifwngeQCHACRX3dypqk/mek2tDIqGQiE8+OCD+O1vf4s33njDtff1+XxpC324hS0RuaHQQdGSpy2q6s9UVVT1fFX9dOIjZ5jXilAohJUrV2Lv3r2uhPnpp5+OjRs3oqmpCfF4HL29vfb0w2x94URE5dJwK0UHBwddfb+BgQG0trYiHo8jHo+nBXdywNM0TfaFE1HZNcReLsltbw8fPoxXX33Vlff8xCc+gTvvvBM9PT2wLCvrICb7womokjy5H3rqNMRDhw65um+5iGDbtm3TNsTiQh8iKpeG3Q89dSDS5/MhEom49t6maeLee+/NurshF/oQUbV5LtBTByLdmIbY3t6OJUuWYMGCBWx9E1FN81ygt7W1ufI+5557LoaGhhjgRFQ3PBXolmVh/fr1jlvm119/PR5++GGXSkVEVBmemrY4NDTkuM/cNE0sWLDApRIREVWOZ1rolmVh27ZtJb1WRNDU1IRYLMb54kRUt+q+hR4KhfDZz34WF154Ycnv8clPfhLhcBh33XUXl+cTUd2q6xZ6chm/U5///Oc57ZCI6l5dttAty8LNN9+M3t7eol87b948LFmyxP5cRNDa2upm8YiIqqLuWui33347BgYGSnqtaZrYtWsXAOD555/nfuNE5Cl1FeilhrmIwOfzYcuWLXa3CvdYISKvqZu9XCzLwkUXXZR3T5bm5macfvrp+NznPocFCxagra0N4+PjDG4iqlue2svFsiz09vbmDXMuCCKiRlbzgZ7cbGumQ5tz7YBIRNRIaj7Qk5ttpTJNE7feeivee+89AMh5ADMRUSOp+UBPnvqT3A53+fLlDHAioixqPtB56g8RUWEcBbqIXA7g2wBMANtV9W5XSpWBqziJiPIreaWoiJgAtgK4AsB8AF8WkfluFYyIiIrjZOn/IgBHVPVNVZ0A8AMA17hTLCIiKpaTQD8LwL+nfH4s8RwREVVB2TfnEpEeERkRkZHjx4+X+3ZERA3LSaD/GsDHUz4/O/FcGlUNqWqnqna2t7c7uB0REc3ESaD/AsA8ETlHRPwAvgTgCXeKRURExXK0OZeIXAlgEFPTFneo6t/nuf44gLcSn84CcKLkm9cGL9QB8EY9WIfa4YV61FodPqGqebs4KrrbYtqNRUYK2T2slnmhDoA36sE61A4v1KNe61CXJxYREdF0DHQiIo+oZqCHqnhvt3ihDoA36sE61A4v1KMu61C1PnQiInIXu1yIiDyiLIEuIpeLyGsickRE7sjy9WYR+YfE118QkbkpX9uUeP41EbmsHOUrRKl1EJHPicg+ETmU+PO/VLrsKWUs+d8h8fU5IvI7EbmtUmXOxuH30/kiYonIK4l/k1MqWfaUcpT6/dQkIrsSZX9VRDZVuuwpZcxXhy4R2S8iURH5YsbXlonIG4mPZZUr9XSl1kNEPp3yvfSSiPx1ZUteAFV19QNTc9L/DcCfAPADeBHA/IxrVgPYlnj8JQD/kHg8P3F9M4BzEu9jul3GMtehA8CZicefAvDrSpffaR1Svv4IgH8EcFs16uDCv4UPwEsA/izxeVsdfj9dB+AHicenAjgKYG6N1mEugPMBDAH4YsrzHwPwZuLPMxKPz6jh76dc9fhTAPMSj88E8B8AWqtRj1wf5WihF7IL4zUAdiUePwJgsYhI4vkfqGpEVX8F4Eji/Sqt5Dqo6gFV/U3i+VcAtIhIc0VKnc7JvwNEpBvArzBVh2pyUo8lAF5S1RcBQFXHVTVWoXKnclIHBfAHIuID0AJgAsB7lSl2mrx1UNWjqvoSgHjGay8D8IyqvqOqvwXwDIDLK1HoLEquh6q+rqpvJB7/BsDbAGpqP5NyBHohuzDa16hqFMD/w1TrqVZ2cHRSh1TXAtivqpEylXMmJddBRD4C4HYA/70C5czHyb/FnwJQEflJ4lfojRUobzZO6vAIgPcx1RocBfC/VPWdchc4Cyf/N2vl/7VrZRGRRZhq4f+bS+VyRc0fQVevRGQBgG9hqpVYb/oA3KOqv0s02OuVD8BfArgAwO8BDIvIPlUdrm6xirIIQAxTv+KfAeB5EXlWVd+sbrEal4j8MYD/A2CZqmb+NlJV5WihF7ILo31N4lfJ0wGMF/jaSnBSB4jI2QD+GcBSVa3WT3AndfgsgAEROQrQ8Q83AAABfUlEQVSgF8CdIrK23AXOwUk9jgF4TlVPqOrvATwJ4M/LXuLpnNThOgD/oqqTqvo2gJ8DqMaSdCf/N2vl/7XjsojIaQB+DODrqvqvLpfNuTIMOvgwNehxDj4cdFiQcc0apA8A/TDxeAHSB0XfRHUGsZzUoTVx/RcqXW636pBxTR+qOyjq5N/iDAD7MTWY6APwLICr6qwOtwPYmXj8BwAOAzi/FuuQcu1DmD4o+qvEv8cZiccfq9Xvpxnq4QcwDKC3GmUvqH5l+ku7EsDrmOpf+nriuf8B4K8Sj0/B1OyJIwD2AviTlNd+PfG61wBcUbW/mBLrAOC/YarP82DKxx/WUx0y3qMPVQx0F76f/gZTA7svAxiotzoA+Eji+VcwFeYbargOF2Dqt6L3MfXbxSspr70hUbcjAJbX+PdT1nokvpcmM/5vf7qadcn84EpRIiKP4EpRIiKPYKATEXkEA52IyCMY6EREHsFAJyLyCAY6EZFHMNCJiDyCgU5E5BH/H4QQ/eZRDrTZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(obs['depth']**0.5, obs['radius']/obs['star_radius'], 'k.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pull observed/physical IDs\n",
    "obs_id = np.array(obs['target_id'], dtype='int')\n",
    "det_id = np.array(det['target_id'], dtype='int')\n",
    "\n",
    "# pull observed periods; make radius & mass arrays\n",
    "obs_per  = np.array(obs['period'], dtype='float')\n",
    "obs_rad  = np.zeros_like(obs_per)\n",
    "obs_mass = np.zeros_like(obs_per)\n",
    "\n",
    "# iterate over all OBSERVED periods\n",
    "for i, per in enumerate(obs_per):\n",
    "    \n",
    "    # pull the DETECTED periods, masses, and radii for this system\n",
    "    det_per_in_this_system  = det['period'][det_id == obs_id[i]]\n",
    "    det_rad_in_this_system  = det['radius'][det_id == obs_id[i]]\n",
    "    det_mass_in_this_system = det['mass'][det_id == obs_id[i]]\n",
    "\n",
    "    \n",
    "    # identify the DETECTED period closest to the OBSERVED period\n",
    "    loc = np.argmin(np.abs(det_per_in_this_system - obs_per[i]))\n",
    "\n",
    "    obs_rad[i] = det_rad_in_this_system[loc]\n",
    "    obs_mass[i] = det_mass_in_this_system[loc]\n",
    "    \n",
    "    \n",
    "# add radius and mass to the OBSERVED planets catalogue\n",
    "obs['radius'] = obs_rad\n",
    "obs['mass'] = obs_mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# count the OBSERVED planets in each system\n",
    "obs_uniquesys = np.unique(obs['target_id'])\n",
    "\n",
    "for us in obs_uniquesys:\n",
    "    obs['npl'][obs['target_id']==us] = int(np.sum(obs['target_id']==us))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6682 planets in UNDETECTED catalogue\n"
     ]
    }
   ],
   "source": [
    "print('%d planets in UNDETECTED catalogue' %len(und['target_id']))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/syssim_undetected_planets_with_masses.csv'\n",
    "\n",
    "    with open(filepath, \"w\") as outfile:\n",
    "        writer = csv.writer(outfile)\n",
    "        writer.writerow(und.keys())\n",
    "        writer.writerows(zip(*und.values()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2183 planets in OBSERVED catalogue\n"
     ]
    }
   ],
   "source": [
    "print('%d planets in OBSERVED catalogue' %len(obs['target_id']))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/syssim_observed_planets_with_masses.csv'\n",
    "\n",
    "    with open(filepath, \"w\") as outfile:\n",
    "        writer = csv.writer(outfile)\n",
    "        writer.writerow(obs.keys())\n",
    "        writer.writerows(zip(*obs.values()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}