build_epos_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/epos_transiting_planets_raw_output.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[0]\n",
    "VALUES = obs_data[1:]\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['star_id'] = getkey('# starID')\n",
    "obs['period']  = getkey('Period_days')\n",
    "obs['radius']  = getkey('Radius_earth')\n",
    "obs['tdur']    = getkey('Duration_hours')/24\n",
    "\n",
    "obs['rhostar'] = np.ones_like(obs['period'])\n",
    "obs['rstar']   = np.ones_like(obs['period'])\n",
    "obs['mstar']   = np.ones_like(obs['period'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEvZJREFUeJzt3X2MXuV55/Hvr3ZCA00WClNk/BLTyEELqHXCiLCbBrFLkxg2CqRaZW3tJjQNdaKaVbK7UhW6f5DtCmlf8rKKllI54AW0iSkNoUErN8HNRmErFRKbeMG8NYaYeqYOdk03NCVyYnPtH3MGnvh1Zp6ZORPf34/0aM5zPfc55xoL9PO573Mep6qQJLXp5/puQJLUH0NAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1LDFfTdwMuecc06tXLmy7zYk6WfG9u3b/6aqRqYydsGHwMqVK9m2bVvfbUjSz4wkz011rNNBktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUsJOGQJJNSfYl2TlQ+6MkO7rX7iQ7uvrKJD8a+OwPB/a5JMljSXYl+VySzM2v9Kqly1eSZMG9li5fOde/uiRNyVS+NuIO4L8Dd00WqupfTG4n+TTwg4Hxz1TV6mMc51bgt4GHgS3AGuBPp9/y1P312HNcf8uBuTzFjNy24ey+W5AkYApXAlX1IPDCsT7r/jb/fmDziY6RZAnwhqp6qKqKiUC5dvrtSpJm07BrAu8Anq+q7w7Uzk/ynSTfTPKOrrYUGBsYM9bVJEk9GvZbRNfx01cBe4EVVXUgySXAnyS5aLoHTbIeWA+wYsWKIVuUJB3PjK8EkiwGfgP4o8laVR2sqgPd9nbgGeDNwDiwbGD3ZV3tmKpqY1WNVtXoyMiUvhJbkjQDw0wH/TrwVFW9Ms2TZCTJom77l4FVwLNVtRd4Mcll3TrCB4GvDHFuSdIsmMotopuBvwAuSDKW5MPdR2s5ekH4cuDR7pbRLwEfrarJReXfAW4DdjFxhTCndwZJkk7upGsCVbXuOPXfPEbtXuDe44zfBlw8zf4kSXPIJ4YlqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktSwk4ZAkk1J9iXZOVD7ZJLxJDu619UDn92YZFeSp5O8e6C+pqvtSvKJ2f9VJEnTNZUrgTuANceof7aqVnevLQBJLgTWAhd1+/xBkkVJFgG3AFcBFwLrurGSpB4tPtmAqnowycopHu8a4O6qOgh8L8ku4NLus11V9SxAkru7sU9Mu2NJ0qwZZk3ghiSPdtNFZ3W1pcCegTFjXe14dUlSj2YaArcCbwJWA3uBT89aR0CS9Um2Jdm2f//+2Ty0JGnAjEKgqp6vqsNV9TLweV6d8hkHlg8MXdbVjlc/3vE3VtVoVY2OjIzMpEVJ0hTMKASSLBl4+z5g8s6h+4G1SU5Lcj6wCvgW8G1gVZLzk7yWicXj+2fetiRpNpx0YTjJZuAK4JwkY8BNwBVJVgMF7AY+AlBVjye5h4kF30PAhqo63B3nBuBrwCJgU1U9Puu/jSRpWqZyd9C6Y5RvP8H4m4Gbj1HfAmyZVneSpDnlE8OS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWrYSUMgyaYk+5LsHKj91yRPJXk0yX1JzuzqK5P8KMmO7vWHA/tckuSxJLuSfC5J5uZXkiRN1VSuBO4A1hxR2wpcXFW/AvwlcOPAZ89U1eru9dGB+q3AbwOruteRx5QkzbOThkBVPQi8cETtgao61L19CFh2omMkWQK8oaoeqqoC7gKunVnLkqTZMhtrAr8F/OnA+/OTfCfJN5O8o6stBcYGxox1NUlSjxYPs3OSfw8cAr7QlfYCK6rqQJJLgD9JctEMjrseWA+wYsWKYVqUJJ3AjK8Ekvwm8B7gX3ZTPFTVwao60G1vB54B3gyM89NTRsu62jFV1caqGq2q0ZGRkZm2KEk6iRmFQJI1wO8C762qlwbqI0kWddu/zMQC8LNVtRd4Mcll3V1BHwS+MnT3kqShnHQ6KMlm4ArgnCRjwE1M3A10GrC1u9Pzoe5OoMuB30/yE+Bl4KNVNbmo/DtM3Gn0OibWEAbXESRJPThpCFTVumOUbz/O2HuBe4/z2Tbg4ml1J0maUz4xLEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhk0pBJJsSrIvyc6B2i8m2Zrku93Ps7p6knwuya4kjyZ568A+13Xjv5vkutn/dSRJ0zHVK4E7gDVH1D4BfL2qVgFf794DXAWs6l7rgVthIjSAm4C3AZcCN00GhySpH1MKgap6EHjhiPI1wJ3d9p3AtQP1u2rCQ8CZSZYA7wa2VtULVfW3wFaODhZJ0jwaZk3g3Kra221/Hzi3214K7BkYN9bVjleXJPVkVhaGq6qAmo1jASRZn2Rbkm379++frcNKko4wTAg8303z0P3c19XHgeUD45Z1tePVj1JVG6tqtKpGR0ZGhmhRknQiw4TA/cDkHT7XAV8ZqH+wu0voMuAH3bTR14B3JTmrWxB+V1eTJPVk8VQGJdkMXAGck2SMibt8/hNwT5IPA88B7++GbwGuBnYBLwEfAqiqF5L8R+Db3bjfr6ojF5slSfNoSiFQVeuO89GVxxhbwIbjHGcTsGnK3UmS5pRPDEtSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNMwQkqWEzDoEkFyTZMfB6McnHk3wyyfhA/eqBfW5MsivJ00nePTu/giRpphbPdMeqehpYDZBkETAO3Ad8CPhsVX1qcHySC4G1wEXAecCfJXlzVR2eaQ+SpOHM1nTQlcAzVfXcCcZcA9xdVQer6nvALuDSWTq/JGkGZisE1gKbB97fkOTRJJuSnNXVlgJ7BsaMdTVJUk+GDoEkrwXeC/xxV7oVeBMTU0V7gU/P4Jjrk2xLsm3//v3DtihJOo7ZuBK4Cnikqp4HqKrnq+pwVb0MfJ5Xp3zGgeUD+y3rakepqo1VNVpVoyMjI7PQ4sKyaPFpJFlwr6XLV/b9RyNpns14YXjAOgamgpIsqaq93dv3ATu77fuBLyb5DBMLw6uAb83C+X/mHD50kOtvOdB3G0e5bcPZfbcgaZ4NFQJJzgDeCXxkoPxfkqwGCtg9+VlVPZ7kHuAJ4BCwwTuDJKlfQ4VAVf09cPYRtQ+cYPzNwM3DnFOSNHt8YliSGmYISFLDDAFJapghIEkNMwQkqWGGgCQ1zBCQpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0bOgSS7E7yWJIdSbZ1tV9MsjXJd7ufZ3X1JPlckl1JHk3y1mHPL0maudm6EvgnVbW6qka7958Avl5Vq4Cvd+8BrgJWda/1wK2zdH5J0gzM1XTQNcCd3fadwLUD9btqwkPAmUmWzFEPkqSTmI0QKOCBJNuTrO9q51bV3m77+8C53fZSYM/AvmNdTZLUg8WzcIxfq6rxJL8EbE3y1OCHVVVJajoH7MJkPcCKFStmoUVJ0rEMfSVQVePdz33AfcClwPOT0zzdz33d8HFg+cDuy7rakcfcWFWjVTU6MjIybIuSpOMYKgSSnJHk9ZPbwLuAncD9wHXdsOuAr3Tb9wMf7O4Sugz4wcC0kSRpng07HXQucF+SyWN9saq+muTbwD1JPgw8B7y/G78FuBrYBbwEfGjI80uShjBUCFTVs8CvHqN+ALjyGPUCNgxzTknS7PGJYUlqmCEgSQ0zBCSpYYaAJDXMEJCkhhkCktQwQ0CSGmYISFLDDAFJapghIEkNm42vktYpYtHi0+i+B2pBOW/ZGxnfs7vvNqRTkiGgVxw+dJDrbznQdxtHuW3D2X23IJ2ynA6SpIYZApLUMENAkhpmCEhSwwwBSWqYISBJDTMEJKlhhoAkNWzGIZBkeZJvJHkiyeNJPtbVP5lkPMmO7nX1wD43JtmV5Okk756NX0CSNHPDPDF8CPh3VfVIktcD25Ns7T77bFV9anBwkguBtcBFwHnAnyV5c1UdHqIHSdIQZnwlUFV7q+qRbvvvgCeBpSfY5Rrg7qo6WFXfA3YBl870/JKk4c3KmkCSlcBbgIe70g1JHk2yKclZXW0psGdgtzFOHBqSpDk2dAgk+QXgXuDjVfUicCvwJmA1sBf49AyOuT7JtiTb9u/fP2yLkqTjGCoEkryGiQD4QlV9GaCqnq+qw1X1MvB5Xp3yGQeWD+y+rKsdpao2VtVoVY2OjIwM06Ik6QSGuTsowO3Ak1X1mYH6koFh7wN2dtv3A2uTnJbkfGAV8K2Znl+SNLxh7g56O/AB4LEkO7ra7wHrkqwGCtgNfASgqh5Pcg/wBBN3Fm3wziBJ6teMQ6Cq/hw41j9DteUE+9wM3DzTc0qSZpdPDEtSwwwBSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1LBhHhaT5sWixacx8YD6wnLesjcyvmd3321IQzEEtOAdPnSQ62850HcbR7ltw9l9tyANzekgSWqYISBJDTMEJKlhhoAkNcwQkKSGGQKS1DBDQJIa5nMC0gz5EJtOBYaANEM+xKZTgSEgnWK8QtF0GALSKcYrFE3HvC8MJ1mT5Okku5J8Yr7PL0l61byGQJJFwC3AVcCFwLokF85nD5KkV833lcClwK6qeraqfgzcDVwzzz1IkjrzvSawFNgz8H4MeNs89yCpBwt1wfo1p53OTw6+1HcbR5mvhfRU1Zyf5JWTJf8cWFNV13fvPwC8rapuOGLcemB99/YC4OkZnvIc4G9muO9csq/psa/psa/pORX7emNVjUxl4HxfCYwDywfeL+tqP6WqNgIbhz1Zkm1VNTrscWabfU2PfU2PfU1P633N95rAt4FVSc5P8lpgLXD/PPcgSerM65VAVR1KcgPwNWARsKmqHp/PHiRJr5r3h8WqaguwZZ5ON/SU0hyxr+mxr+mxr+lpuq95XRiWJC0sfpW0JDXslAyBJJuS7Euys+9eJiVZnuQbSZ5I8niSj/XdE0CSn0/yrST/t+vrP/Td06Aki5J8J8n/6ruXQUl2J3ksyY4k2/ruZ1KSM5N8KclTSZ5M8o8WQE8XdH9Ok68Xk3y8774Akvyb7r/7nUk2J/n5vnsCSPKxrqfH5/rP6pScDkpyOfBD4K6qurjvfgCSLAGWVNUjSV4PbAeuraoneu4rwBlV9cMkrwH+HPhYVT3UZ1+TkvxbYBR4Q1W9p+9+JiXZDYxW1YK6vzzJncD/qarbujvwTq+q/9d3X5O6r44ZZ+L5oOd67mUpE/+9X1hVP0pyD7Clqu7oua+Lmfg2hUuBHwNfBT5aVbvm4nyn5JVAVT0IvNB3H4Oqam9VPdJt/x3wJBNPUPeqJvywe/ua7rUg/maQZBnwz4Db+u7lZ0GSfwBcDtwOUFU/XkgB0LkSeKbvABiwGHhdksXA6cBf99wPwD8EHq6ql6rqEPBN4Dfm6mSnZAgsdElWAm8BHu63kwndlMsOYB+wtaoWRF/AfwN+F3i570aOoYAHkmzvnnBfCM4H9gP/o5tCuy3JGX03dYS1wOa+mwCoqnHgU8BfAXuBH1TVA/12BcBO4B1Jzk5yOnA1P/2Q7awyBOZZkl8A7gU+XlUv9t0PQFUdrqrVTDzBfWl3OdqrJO8B9lXV9r57OY5fq6q3MvGNuBu6Kci+LQbeCtxaVW8B/h5YMF/X3k1PvRf44757AUhyFhNfYHk+cB5wRpJ/1W9XUFVPAv8ZeICJqaAdwOG5Op8hMI+6Ofd7gS9U1Zf77udI3dTBN4A1ffcCvB14bzf3fjfwT5P8z35belX3t0iqah9wHxPzt30bA8YGruS+xEQoLBRXAY9U1fN9N9L5deB7VbW/qn4CfBn4xz33BEBV3V5Vl1TV5cDfAn85V+cyBOZJtwB7O/BkVX2m734mJRlJcma3/TrgncBT/XYFVXVjVS2rqpVMTCH876rq/W9pAEnO6Bb36aZb3sXEJXyvqur7wJ4kF3SlK4Febzw4wjoWyFRQ56+Ay5Kc3v3/eSUTa3W9S/JL3c8VTKwHfHGuznVK/vOSSTYDVwDnJBkDbqqq2/vtircDHwAe6+bfAX6ve4K6T0uAO7u7Nn4OuKeqFtTtmAvQucB93dciLwa+WFVf7belV/xr4Avd1MuzwId67gd4JSzfCXyk714mVdXDSb4EPAIcAr7Dwnl6+N4kZwM/ATbM5QL/KXmLqCRpapwOkqSGGQKS1DBDQJIaZghIUsMMAUlqmCEgSQ0zBCSpYYaAJDXs/wN81Gg6Yq/RTQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "npl = np.zeros_like(obs['star_id'], dtype='int')\n",
    "\n",
    "for i, sid in enumerate(obs['star_id']):\n",
    "    npl[i] = np.sum(obs['star_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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3176 planets in OBSERVED catalogue\n"
     ]
    }
   ],
   "source": [
    "print('%d planets in OBSERVED catalogue' %len(obs['star_id']))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/epos_observed_planets.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": "markdown",
   "metadata": {},
   "source": [
    "# Build the underlying physical catalogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# READ IN PHYSICAL PLANETS\n",
    "phys_file = MAINPATH + 'catalogues/epos_transiting_systems_raw_output.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[0]\n",
    "VALUES = phys_data[1:]\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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "phys = {}\n",
    "\n",
    "phys['star_id'] = getkey('# starID')\n",
    "phys['period']  = getkey('Period_days')\n",
    "phys['radius']  = getkey('Radius_earth')\n",
    "phys['tdur']    = getkey('Duration_hours')/24\n",
    "\n",
    "phys['rhostar'] = np.ones_like(phys['period'])\n",
    "phys['rstar']   = np.ones_like(phys['period'])\n",
    "phys['mstar']   = np.ones_like(phys['period'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFC1JREFUeJzt3X+MXeV95/H3pzikCV1hA7OG2maNFDcRWSlAR0CWVbWLG2NoFfNHSonYYCGvvH94m6Sq1IX9By0kKyJVpUHasrLArcmycVmaCCuLQl1DVa20EMyPJQESeQpxPF5jTzGQNmyTQr/7x30MgzvTuQMzc50+75d0dZ/zPc855zlH8nx8ftx7U1VIkvrzM6MegCRpNAwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqeWDdMpyW8C/xYo4NvADcA5wC7gTOAJ4DNV9ZMk7wfuAX4ReBn49ar6flvPTcAW4E3gs1X10D+03bPOOqvWrl07/72SpI498cQTf1lVY3P1mzMAkqwCPgucX1X/L8l9wLXAVcDtVbUryX9l8If9zvb+SlV9KMm1wJeAX09yflvuo8DPA3+a5Beq6s3Ztr127Vr27ds3585Kkt6W5MAw/Ya9BLQM+ECSZcAHgcPA5cD9bf5O4OrW3tSmafPXJ0mr76qqH1fVi8AEcPGQ25ckLbA5A6CqDgG/A/yAwR/+1xhc8nm1qt5o3SaBVa29CjjYln2j9T9zen2GZSRJS2zOAEiygsH/3s9jcOnmNGDjYg0oydYk+5Lsm5qaWqzNSFL3hrkE9MvAi1U1VVV/C3wNuAxY3i4JAawGDrX2IWANQJt/OoObwW/VZ1jmLVW1varGq2p8bGzOexiSpHdpmAD4AXBpkg+2a/nrgeeAR4BPtT6bgQdae3ebps1/uAY/OrAbuDbJ+5OcB6wDvrUwuyFJmq85nwKqqseS3A88CbwBPAVsB/4nsCvJF1rt7rbI3cBXkkwAxxg8+UNVPdueIHqurWfbP/QEkCRpceVk/kWw8fHx8jFQSZqfJE9U1fhc/fwksCR1ygCQpE4ZAJI0zdo1Z5Nk5K+1a85e9H0d6ruAJKkXByaPUPeOehSQ644s+jY8A5CkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnZozAJJ8OMnT014/TPL5JGck2ZNkf3tf0fonyR1JJpI8k+Siaeva3PrvT7J59q1KkhbbnAFQVd+rqguq6gLgF4HXga8DNwJ7q2odsLdNA1wJrGuvrcCdAEnOAG4GLgEuBm4+HhqSpKU330tA64G/qKoDwCZgZ6vvBK5u7U3APTXwKLA8yTnAFcCeqjpWVa8Ae4CN73kPJEnvynwD4Frgq629sqoOt/ZLwMrWXgUcnLbMZKvNVpckjcDQAZDkVOCTwP84cV5VFVALMaAkW5PsS7JvampqIVYpSZrBfM4ArgSerKrjP1R5pF3aob0fbfVDwJppy61utdnq71BV26tqvKrGx8bG5jE8SdJ8zCcAPs3bl38AdgPHn+TZDDwwrX59exroUuC1dqnoIWBDkhXt5u+GVpMkjcCyYTolOQ34BPDvppVvA+5LsgU4AFzT6g8CVwETDJ4YugGgqo4luRV4vPW7paqOvec9kCS9K0MFQFX9CDjzhNrLDJ4KOrFvAdtmWc8OYMf8hylJWmh+EliSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqeGCoAky5Pcn+S7SZ5P8vEkZyTZk2R/e1/R+ibJHUkmkjyT5KJp69nc+u9Psnn2LUqSFtuwZwBfBr5ZVR8BPgY8D9wI7K2qdcDeNg1wJbCuvbYCdwIkOQO4GbgEuBi4+XhoSJKW3pwBkOR04JeAuwGq6idV9SqwCdjZuu0Erm7tTcA9NfAosDzJOcAVwJ6qOlZVrwB7gI0LujeSpKENcwZwHjAF/EGSp5LcleQ0YGVVHW59XgJWtvYq4OC05Sdbbbb6OyTZmmRfkn1TU1Pz2xtJ0tCGCYBlwEXAnVV1IfAj3r7cA0BVFVALMaCq2l5V41U1PjY2thCrlCTNYJgAmAQmq+qxNn0/g0A40i7t0N6PtvmHgDXTll/darPVJUkjMGcAVNVLwMEkH26l9cBzwG7g+JM8m4EHWns3cH17GuhS4LV2qeghYEOSFe3m74ZWkySNwLIh+/0GcG+SU4EXgBsYhMd9SbYAB4BrWt8HgauACeD11peqOpbkVuDx1u+Wqjq2IHshSZq3oQKgqp4GxmeYtX6GvgVsm2U9O4Ad8xmgJGlx+ElgSeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6tRQAZDk+0m+neTpJPta7Ywke5Lsb+8rWj1J7kgykeSZJBdNW8/m1n9/ks2zbU+StPjmcwbwr6vqgqo6/tOQNwJ7q2odsLdNA1wJrGuvrcCdMAgM4GbgEuBi4ObjoSFJWnrv5RLQJmBna+8Erp5Wv6cGHgWWJzkHuALYU1XHquoVYA+w8T1sX5L0HgwbAAX8SZInkmxttZVVdbi1XwJWtvYq4OC0ZSdbbba6JGkElg3Z719W1aEk/xTYk+S702dWVSWphRhQC5itAOeee+5CrFKSNIOhzgCq6lB7Pwp8ncE1/CPt0g7t/WjrfghYM23x1a02W/3EbW2vqvGqGh8bG5vf3kiShjZnACQ5Lck/Od4GNgDfAXYDx5/k2Qw80Nq7gevb00CXAq+1S0UPARuSrGg3fze0miRpBIa5BLQS+HqS4/3/e1V9M8njwH1JtgAHgGta/weBq4AJ4HXgBoCqOpbkVuDx1u+Wqjq2YHsiSZqXOQOgql4APjZD/WVg/Qz1ArbNsq4dwI75D1OStND8JLAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4NHQBJTknyVJJvtOnzkjyWZCLJHyU5tdXf36Yn2vy109ZxU6t/L8kVC70zkqThzecM4HPA89OmvwTcXlUfAl4BtrT6FuCVVr+99SPJ+cC1wEeBjcDvJznlvQ1fkvRuDRUASVYDvwLc1aYDXA7c37rsBK5u7U1tmjZ/feu/CdhVVT+uqheBCeDihdgJSdL8DXsG8HvAbwN/16bPBF6tqjfa9CSwqrVXAQcB2vzXWv+36jMsI0laYnMGQJJfBY5W1RNLMB6SbE2yL8m+qamppdikJHVpmDOAy4BPJvk+sIvBpZ8vA8uTLGt9VgOHWvsQsAagzT8deHl6fYZl3lJV26tqvKrGx8bG5r1DkqThzBkAVXVTVa2uqrUMbuI+XFXXAY8An2rdNgMPtPbuNk2b/3BVVatf254SOg9YB3xrwfZEkjQvy+buMqv/AOxK8gXgKeDuVr8b+EqSCeAYg9Cgqp5Nch/wHPAGsK2q3nwP25ckvQfzCoCq+jPgz1r7BWZ4iqeq/gb4tVmW/yLwxfkOUpK08N7LGYAkLZi1a87mwOSRUQ+jKwaApJPCgckj1L2jHgXkulGPYOn4XUCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkTvlBMKlzfgK3XwaA1Dk/gdsvLwFJUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkTs0ZAEl+Nsm3kvyfJM8m+U+tfl6Sx5JMJPmjJKe2+vvb9ESbv3baum5q9e8luWKxdkqSNLdhzgB+DFxeVR8DLgA2JrkU+BJwe1V9CHgF2NL6bwFeafXbWz+SnM/gB+I/CmwEfj/JKQu5M5Kk4c0ZADXw123yfe1VwOXA/a2+E7i6tTe1adr89UnS6ruq6sdV9SIwwQw/Ki9JWhpD3QNIckqSp4GjwB7gL4BXq+qN1mUSWNXaq4CDAG3+a8CZ0+szLCNJWmJDBUBVvVlVFwCrGfyv/SOLNaAkW5PsS7JvampqsTYjSd2b11NAVfUq8AjwcWB5kuPfJbQaONTah4A1AG3+6cDL0+szLDN9G9uraryqxsfGxuYzPEnSPAzzFNBYkuWt/QHgE8DzDILgU63bZuCB1t7dpmnzH66qavVr21NC5wHrgG8t1I5IkuZnmG8DPQfY2Z7Y+Rngvqr6RpLngF1JvgA8Bdzd+t8NfCXJBHCMwZM/VNWzSe4DngPeALZV1ZsLuzuSpGHNGQBV9Qxw4Qz1F5jhKZ6q+hvg12ZZ1xeBL85/mJKkheYngSWpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdWqYH4Vfk+SRJM8leTbJ51r9jCR7kuxv7ytaPUnuSDKR5JkkF01b1+bWf3+SzbNtU5K0+IY5A3gD+K2qOh+4FNiW5HzgRmBvVa0D9rZpgCuBde21FbgTBoEB3AxcwuC3hG8+HhqSpKU3ZwBU1eGqerK1/wp4HlgFbAJ2tm47gatbexNwTw08CixPcg5wBbCnqo5V1SvAHmDjgu6N9FNm7ZqzSTLSl/q1bD6dk6wFLgQeA1ZW1eE26yVgZWuvAg5OW2yy1Warn7iNrQzOHDj33HPnMzzpp86BySPUvaMdQ64b7fY1OkPfBE7yc8AfA5+vqh9On1dVBdRCDKiqtlfVeFWNj42NLcQqJUkzGCoAkryPwR//e6vqa618pF3aob0fbfVDwJppi69utdnqkqQRGOYpoAB3A89X1e9Om7UbOP4kz2bggWn169vTQJcCr7VLRQ8BG5KsaDd/N7SaJGkEhrkHcBnwGeDbSZ5utf8I3Abcl2QLcAC4ps17ELgKmABeB24AqKpjSW4FHm/9bqmqYwuyF5KkeZszAKrqfwGzPSqwfob+BWybZV07gB3zGaAkaXH4SWBJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnq1DC/CbwjydEk35lWOyPJniT72/uKVk+SO5JMJHkmyUXTltnc+u9PsnmmbUmSls4wZwB/CGw8oXYjsLeq1gF72zTAlcC69toK3AmDwABuBi4BLgZuPh4akqTRmDMAqurPgRN/vH0TsLO1dwJXT6vfUwOPAsuTnANcAeypqmNV9Qqwh78fKpKkJfRu7wGsrKrDrf0SsLK1VwEHp/WbbLXZ6pKkEXnPN4GrqoBagLEAkGRrkn1J9k1NTS3UaiVJJ3i3AXCkXdqhvR9t9UPAmmn9VrfabPW/p6q2V9V4VY2PjY29y+FJkubybgNgN3D8SZ7NwAPT6te3p4EuBV5rl4oeAjYkWdFu/m5oNUnSiCybq0OSrwL/CjgrySSDp3luA+5LsgU4AFzTuj8IXAVMAK8DNwBU1bEktwKPt363VNWJN5YlSUtozgCoqk/PMmv9DH0L2DbLenYAO+Y1OknSovGTwJLUqTnPAKR/jNauOZsDk0dGPQxppAwAdenA5BHq3lGPAnLdqEegnnkJSJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1yq+D1pLye/ilk8eSB0CSjcCXgVOAu6rqtqUeg0bH7+GXTh5LegkoySnAfwGuBM4HPp3k/KUcgyRpYKnvAVwMTFTVC1X1E2AXsGmJxyBJYukDYBVwcNr0ZKtpka1dczZJRv6SdPI46W4CJ9kKbG2Tf53ke6MczwI4C/jLUQ/iJHJWrjs5jsdJch/gpDgeHot3OmmOR/Juj8c/G6bTUgfAIWDNtOnVrfaWqtoObF/KQS2mJPuqanzU4zhZeDzeyePxNo/FOy3F8VjqS0CPA+uSnJfkVOBaYPcSj0GSxBKfAVTVG0n+PfAQg8dAd1TVs0s5BknSwJLfA6iqB4EHl3q7I/SP5nLWAvF4vJPH420ei3da9OORqlrsbUiSTkJ+F5AkdcoAWCRJ1iR5JMlzSZ5N8rlRj2nUkpyS5Kkk3xj1WEYtyfIk9yf5bpLnk3x81GMapSS/2f6dfCfJV5P87KjHtJSS7EhyNMl3ptXOSLInyf72vmKht2sALJ43gN+qqvOBS4Ftfu0FnwOeH/UgThJfBr5ZVR8BPkbHxyXJKuCzwHhV/XMGD4hcO9pRLbk/BDaeULsR2FtV64C9bXpBGQCLpKoOV9WTrf1XDP6Bd/up5ySrgV8B7hr1WEYtyenALwF3A1TVT6rq1dGOauSWAR9Isgz4IPB/RzyeJVVVfw4cO6G8CdjZ2juBqxd6uwbAEkiyFrgQeGy0Ixmp3wN+G/i7UQ/kJHAeMAX8QbskdleS00Y9qFGpqkPA7wA/AA4Dr1XVn4x2VCeFlVV1uLVfAlYu9AYMgEWW5OeAPwY+X1U/HPV4RiHJrwJHq+qJUY/lJLEMuAi4s6ouBH7EIpze/7Ro17Y3MQjGnwdOS/JvRjuqk0sNHtdc8Ec2DYBFlOR9DP7431tVXxv1eEboMuCTSb7P4BtgL0/y30Y7pJGaBCar6vgZ4f0MAqFXvwy8WFVTVfW3wNeAfzHiMZ0MjiQ5B6C9H13oDRgAiySDr768G3i+qn531OMZpaq6qapWV9VaBjf3Hq6qbv+HV1UvAQeTfLiV1gPPjXBIo/YD4NIkH2z/btbT8U3xaXYDm1t7M/DAQm/AAFg8lwGfYfC/3afb66pRD0onjd8A7k3yDHAB8J9HPJ6RaWdC9wNPAt9m8Hepq08FJ/kq8L+BDyeZTLIFuA34RJL9DM6SFvzXE/0ksCR1yjMASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqf+P4G3zRR7FsAfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "npl = np.zeros_like(phys['star_id'], dtype='int')\n",
    "\n",
    "for i, sid in enumerate(phys['star_id']):\n",
    "    npl[i] = np.sum(phys['star_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": [
      "18039 planets in PHYSICAL catalogue\n"
     ]
    }
   ],
   "source": [
    "print('%d planets in PHYSICAL catalogue' %len(phys['star_id']))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/epos_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": [
    "# Isolate the undetected planets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "und = phys.copy()\n",
    "\n",
    "in_observed_catalogue = np.isin(phys['period'], obs['period'])\n",
    "\n",
    "for k in und.keys():\n",
    "    und[k] = und[k][~in_observed_catalogue]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14863 planets in UNDETECTED catalogue\n"
     ]
    }
   ],
   "source": [
    "print('%d planets in UNDETECTED catalogue' %len(und['star_id']))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/epos_undetected_planets.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": "markdown",
   "metadata": {},
   "source": [
    "# Apply limit cuts to the undetected objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# restore all the physical planets\n",
    "und = phys.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1842 long-period planets (P > 400d) discarded\n"
     ]
    }
   ],
   "source": [
    "# first remove long period planets (P > 400)\n",
    "keep = (und['period'] <= 400.)\n",
    "\n",
    "for k in und.keys():\n",
    "    und[k] = und[k][keep]\n",
    "    \n",
    "print('%d long-period planets (P > 400d) discarded' %np.sum(~keep))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "142 systems removed after being reduced to less than 3 planets\n"
     ]
    }
   ],
   "source": [
    "# update number of planets in each system to reflect removals\n",
    "uniquesys = np.unique(und['star_id'])\n",
    "\n",
    "for us in uniquesys:\n",
    "    und['npl'][und['star_id']==us] = int(np.sum(und['star_id']==us))\n",
    "\n",
    "    \n",
    "# remove systems with fewer than 3 planets\n",
    "npl = und['npl'].copy()\n",
    "\n",
    "for k in und.keys():\n",
    "    und[k] = und[k][npl >= 3]\n",
    "    \n",
    "print('%d systems removed after being reduced to less than 3 planets' %np.sum(npl <= 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "656 systems discarded with undetected planets outside valid range\n"
     ]
    }
   ],
   "source": [
    "# elimnate SYSTEMS which do not match period (P > 2) or radius (R < 6) limits\n",
    "bad_planet = (und['radius'] > 6.0) + (und['period'] < 2.0)\n",
    "bad_system = np.zeros_like(bad_planet)\n",
    "\n",
    "\n",
    "uniquesys = np.unique(und['star_id'])\n",
    "\n",
    "for us in uniquesys:\n",
    "    bad_system[und['star_id']==us] = int(np.sum(bad_planet[und['star_id']==us]))\n",
    "\n",
    "for k in und.keys():\n",
    "    und[k] = und[k][~bad_system]\n",
    "    \n",
    "print('%d systems discarded with undetected planets outside valid range' %(len(uniquesys)-len(np.unique(und['star_id']))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3 systems discarded with observed planets outside valid range\n"
     ]
    }
   ],
   "source": [
    "# check for SYSTEMS which do not match limits for detected planets either\n",
    "bad_planet = (obs['radius'] > 6.0) + (obs['period'] < 2.0)\n",
    "bad_system = np.isin(und['star_id'], obs['star_id'][bad_planet])\n",
    "\n",
    "uniquesys = np.unique(und['star_id'])\n",
    "\n",
    "for k in und.keys():\n",
    "    und[k] = und[k][~bad_system]\n",
    "    \n",
    "print('%d systems discarded with observed planets outside valid range' %(len(uniquesys)-len(np.unique(und['star_id']))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# seperate out the undetected planets\n",
    "in_observed_catalogue = np.isin(und['period'], obs['period'])\n",
    "\n",
    "for k in und.keys():\n",
    "    und[k] = und[k][~in_observed_catalogue]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8610 planets in LIMITED UNDETECTED catalogue\n"
     ]
    }
   ],
   "source": [
    "print('%d planets in LIMITED UNDETECTED catalogue' %len(und['star_id']))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/epos_undetected_planets_with_limit_cuts.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": "markdown",
   "metadata": {},
   "source": [
    "# Apply limit cuts to the observed objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 long-period planets (P > 400d) discarded\n"
     ]
    }
   ],
   "source": [
    "# first remove long period planets (P > 400)\n",
    "keep = (obs['period'] <= 400.)\n",
    "\n",
    "for k in obs.keys():\n",
    "    obs[k] = obs[k][keep]\n",
    "    \n",
    "print('%d long-period planets (P > 400d) discarded' %np.sum(~keep))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18 small planets (rp/Re < 0.5) discarded\n"
     ]
    }
   ],
   "source": [
    "# next remove small planets (R < 0.5)\n",
    "keep = (obs['radius'] >= 0.5)\n",
    "\n",
    "for k in obs.keys():\n",
    "    obs[k] = obs[k][keep]\n",
    "    \n",
    "print('%d small planets (rp/Re < 0.5) discarded' %np.sum(~keep))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2499 systems removed after being reduced to less than 3 planets\n"
     ]
    }
   ],
   "source": [
    "# update number of planets in each system to reflect removals\n",
    "uniquesys = np.unique(obs['star_id'])\n",
    "\n",
    "for us in uniquesys:\n",
    "    obs['npl'][obs['star_id']==us] = int(np.sum(obs['star_id']==us))\n",
    "\n",
    "    \n",
    "# remove systems reduced to less than 3 planets\n",
    "npl = obs['npl'].copy()\n",
    "\n",
    "for k in obs.keys():\n",
    "    obs[k] = obs[k][npl >= 3]\n",
    "    \n",
    "print('%d systems removed after being reduced to less than 3 planets' %np.sum(npl < 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "74 systems discarded\n"
     ]
    }
   ],
   "source": [
    "# finally elimnate SYSTEMS which do not match period (P > 2) or radius (R < 6) limits\n",
    "bad_planet = (obs['radius'] > 6.0) + (obs['period'] < 2.0)\n",
    "bad_system = np.zeros_like(bad_planet)\n",
    "\n",
    "\n",
    "uniquesys = np.unique(obs['star_id'])\n",
    "\n",
    "for us in uniquesys:\n",
    "    bad_system[obs['star_id']==us] = int(np.sum(bad_planet[obs['star_id']==us]))\n",
    "\n",
    "for k in obs.keys():\n",
    "    obs[k] = obs[k][~bad_system]\n",
    "    \n",
    "print('%d systems discarded' %(len(uniquesys)-len(np.unique(obs['star_id']))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "380 planets in LIMITED OBSERVED catalogue\n"
     ]
    }
   ],
   "source": [
    "print('%d planets in LIMITED OBSERVED catalogue' %len(obs['star_id']))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/epos_observed_planets_with_limit_cuts.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": "markdown",
   "metadata": {},
   "source": [
    "# Read in masses and combine for OBSERVED MULTIS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# READ IN EPOS DATABASE\n",
    "epos_file = MAINPATH + 'catalogues/epos_observed_planets.csv'\n",
    "\n",
    "epos_data = []\n",
    "with open(epos_file) as infile:\n",
    "    reader = csv.reader(infile)\n",
    "    \n",
    "    for row in reader:\n",
    "        epos_data.append(row)\n",
    "        \n",
    "KEYS = epos_data[0]\n",
    "VALUES = epos_data[1:]\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": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PULL RELEVANT QUANTITIES FROM CATALOGUE\n",
    "epos = {}\n",
    "\n",
    "epos['star_id'] = np.array(getkey('star_id'), dtype='int')\n",
    "\n",
    "epos['period'] = getkey('period')\n",
    "epos['radius'] = getkey('radius')\n",
    "epos['tdur']   = getkey('tdur')\n",
    "\n",
    "epos['rhostar'] = np.ones_like(epos['period'])\n",
    "epos['rstar']   = np.ones_like(epos['period'])\n",
    "epos['mstar']   = np.ones_like(epos['period'])\n",
    "\n",
    "epos['npl'] = np.array(getkey('npl'), dtype='int')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove single-planet systems\n",
    "epos_npl = epos['npl'].copy()\n",
    "\n",
    "for k in epos.keys():\n",
    "    epos[k] = epos[k][epos_npl > 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "r1 = epos['radius']\n",
    "r2 = np.load(MAINPATH + 'catalogues/epos_observed_pred_masses_oct22.npy')[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3WuMZOld3/Hv/5xTp25dfe+5z+zMXrzr9d1MHGOHq01kwMJISQQIkB2Q9k0ChhAhDEp4EREhgQiOiBKtwNgRlklkjLBQAG9MjI3ZtZm115e9787uzM6175e6njrnPHlR1e2enr53dVd19e8j7U53dXXVMz3Vv376f57n/5hzDhER6R9etwcgIiKdpWAXEekzCnYRkT6jYBcR6TMKdhGRPqNgFxHpMwp2EZE+o2AXEekzCnYRkT4THOSTjY+Pu/Pnzx/kU4qIHHpPPvnktHNuYrv3P9BgP3/+PJcuXTrIpxQROfTM7MpO7q9SjIhIn1Gwi4j0GQW7iEifUbCLiPQZBbuISJ9RsIuI9JkDXe4oInIURXFKM0nJ+B5hsP/zaQW7iMg+iuKU63NVHGDA6ZHCvoe7SjEiIvuomaQ4oJgNcO3395uCXURkH2V8DwMqjRhrv7/fVIoREdlHYeBxeqSgGruISD8Jg4MJ9GUKdhHpCwe98qSXKdhF5NDrxsqTXnZ0/+Yi0je6sfKklynYReTQ68bKk16mUoyIHHrdWHnSyxTsItIXDnrlSS/TV0FEpM8o2EVE+syWwW5mHzOzSTP79qrbfsfMnjOzb5rZn5vZ8P4OU0REtms7M/aPA+9bc9tjwBudc28GXgA+0uFxiYjILm0Z7M65LwKza277nHMubr/7BHBmH8YmIiK70Ika+88Bf9WBxxERkQ7YU7Cb2W8AMfDJTe7ziJldMrNLU1NTe3k6ERHZhl0Hu5l9CHg/8NPOObfR/ZxzjzrnLjrnLk5MTOz26UREZJt2tUHJzN4H/Crwfc65ameHJCIie7Gd5Y6fAh4HHjSza2b288AfACXgMTN7ysz+xz6PU0REtmnLGbtz7qfWufmP9mEsIiLSAdp5KiLSZ9QETER6gk5A6hwFu4h0nU5A6ix95USk63QCUmcp2EWk63QCUmepFCMiXacTkDpLwS4iPUEnIHWOvooiIn1GwS4i0mcU7CIifUbBLiLSZxTsIiJ9RsEuItJnFOwiIn1GwS4i0mcU7CIifUbBLiLSZ9RSQEQ6Rj3Ve4OCXUQ6Qj3Ve4e+6iLSEeqp3jsU7CLSEeqp3jtUihGRjlBP9d6hYBeRjlFP9d6w5b+AmX3MzCbN7Nurbhs1s8fM7MX2nyP7O0wREdmu7fxo/TjwvjW3/RrweefcA8Dn2++LyCEUxSmVRkwU62Jnv9gy2J1zXwRm19z8AeAT7bc/Afx4h8clIgdgeYni7cU61+eqCvc+sdti2HHn3M3227eA4x0aj4gcIC1R7E97vsrhnHOA2+jjZvaImV0ys0tTU1N7fToR6SAtUexPu/1XvG1mJwHaf05udEfn3KPOuYvOuYsTExO7fDoR2Q/LSxSPD+a0U7SP7PZf8bPAB9tvfxD4i84MR0QOWhh4FLOBQr2PbGe546eAx4EHzeyamf088NvAD5nZi8B72++LiEgP2HKDknPupzb40Hs6PBaRI02dEaVTtPNUpAeoM6J0kl45Ij1gP5cd7scGJG1q6m2asYv0gP1adrgfvwnot4vep2AX6QFbdUbcbf199W8ClUZMM0n3HML78ZjSWQp2kR6xUWfEvcyQ9+M3AW1q6n0KdpEet5cZ8n70SFff9d6nYBfpcXudIe9Hj3T1Xe9tCnaRHqcZsuyUgl3kENAMWXZCrxQRkT6jYBcR6TMKdhGRPqNgFxHpMwp2EZE+o1UxIj1oJy0E1O5X1lKwi/SYnbQQUEMuWY9eASI9ZictfPez3a8cXgp2kR6zkxYCasgl61EpRqTH7KSFgNoNyHoU7CI9aCctBNRuQNZSsIscAK1ckYOkYBfZB6uDHNDKFTlQe3p1mdkvm9nTZvZtM/uUmeU6NTCRw2p5CeLtxTrX56pUGrFWrsiB2nWwm9lp4BeBi865NwI+8JOdGpjIYbV2CSKglStyoPZaigmAvJk1gQJwY+9DEjm8ojglilPiJF0J8mI2oJgNVGOXA7PrYHfOXTez3wWuAjXgc865z3VsZCKHzOpdoAAjhZBiNlgJcgW6HJS9lGJGgA8AF4BTQNHMfmad+z1iZpfM7NLU1NTuRyrS41aXYIL2zFxhLt2wl1fde4FXnHNTzrkm8BngXWvv5Jx71Dl30Tl3cWJiYg9PJ9LbtAtUesVeauxXgXeaWYFWKeY9wKWOjErkENIuUOkVe6mxf8XMPg18DYiBrwOPdmpgIofJ6nXrxay2h0h37ekV6Jz7TeA3OzQWkUNpu61ztftUDoqmFiJ7tPqiaaUR00zSu4JbfdPlIOmVJX0pilvryKN4/3d5bueiqfqmy0HSjF36zkHPjrdz0VQrZuQgKdil72ynNNJpW61Z14oZOUgKduk7vTo71oYlOSgKduk7250da5WK9CsFu/SlrWbHm9XhtxP4+qEgvUzBLkfSch0+DDwWqk0qjZgwCLd14VVLF6XXKdilq7o18834HnGScm22Cjgyvq201t3qwms3Ls6K7ISCXbqmmzPfMPCYKOVoJo7BfIZm8p0fMFtdeO3Vi7MiyxTs0jXdnvkWswGF0KeZpCsBvfrCq3Pf2Ui0elxauii9TsEuXdPtme9GAb3852a/TWjpovQyBbt0Tbdmvmvr+us9bzNJiZKU0PeI2mUaBbkcFgp26aqDnPku94+ZWqoTtH9b2Kiu7xxMLtRJcXgYp4cLBzJGkU5QsMuRsHyhtholzJTrXJgorVwwXS/YzeDYUI4w8IjiFLMuDFpkl/S7pRxKO+3euHyhdqiQAYzFWnPTun7G9wh9D+cg9D2tfJFDRTN2OXR2c7DF8oXaKE4ZGwgZLmQYyocbloG08kUOMwW79JTtbFhab5nk8u3Ln7de+J8eKVBpxMRJSiNOmVqqEwYbr53Xyhc5rBTs0jO2OxNfu0zSubuXJq4X/q2dpR5B+1xS7RqVfqVgl56x3Q1La8sk632ec1CPktZjrKqRd3vtvMhBULBLz9hJ6C6XSaI4JYpT4iS9YwY/tVTH940kcUwM5+7YfKTaufQ7Bbv0jJ2GbhSnvDpdppk4wDFRyt3RyGu4ELbCfs1Sxe3WztWaVw4rBbv0lM1CdzlonWutM6/UY24t1MiFAfUo5tTwd2ryey23qDWvHGZ7CnYzGwb+EHgj4ICfc8493omBydGz2Qx5OWijJGVyoc6xoVyrhp468kArfls6UW7pdoMykb3Y64z9o8BfO+f+pZmFgPZdy65sNUNeORjD90hxhIFHM3YUMwEZ32NoKEMx+52X816XKuoiqxxmuw52MxsCvhf4EIBzLgKizgxLjpqtZsgrG4ySFA+j0oiZq0SMDIR42B1lmE7QRVY5zPYyY78ATAF/bGZvAZ4EPuycq6y+k5k9AjwCcO7cuT08nfSz9damVxrxHR0Yl4P29HCBatT62PIF0mZ7VUwnQ1gblOSwMufc7j7R7CLwBPBu59xXzOyjwKJz7j9s9DkXL150ly5d2t1Ipe+tvjh6Y75KM2kdWXd+fGDdmvvqFTEZ39uyY6PIYWVmTzrnLm73/nt59V8DrjnnvtJ+/9PA2/fweNKHNmrWtd7tYeCtLFe8NltlptLg2myVSiPe9DmaiVsp4zi+c+qRyFG161KMc+6Wmb1mZg86554H3gM807mhyWG30QXR9W4HVmbrC7WIW4t1CmFAFCdUo5gw8O4oszSTlMD3GCoEzFcjksTpQqdI215XxfwC8Mn2ipjLwL/e+5CkX2zUrGtqqc5Crcl4KUsUp8xXI6aXGqQ45soRudAHoJTLEKceU0sNGnF6xw+H1TX50PeYGM5hhi50irDHYHfOPQVsu+4jR8t6F0RfnS7z6lSV24t1jg2FnBkpsFhrslCLSBxUo5jRgZDjQzmGcgEZPyQX+netltGqFZGNaeepdHzr/OrHW9usq9JIqMVN8jmPRtNRDAMaScpiDa5MlcGMfMbnxGCOU8MFMr7H1FJ93TKLVq2IrE/BfsR1euv82sebKLVKJMthnzhHpZ4QBB650ChkA6JqRLkR43s+p0dyjBWznBouMFIMAQgDzcxFdkLBfsR1euv86sebr0Zcnang+7bSBiDwYCDvE5hPEHgUwoBmnJLL+JwZCfD91u7StbtIlx979fsisj4F+xG3Ud/y3VpdV08Sh+8bYdBuA+B7xL7P2ZEBClmfaiNhoRrxwq1FpsoN0sRxfmKA8YHsHY+phlwiO6NgP8Ki9vFw6/Ut363VFzWX+6JXGjHVRqsTIxgYzFYbZDyPRjMhMXjdiUFuzFUwg3IUU4nilQBXQy6RnVGwH2Fb9S3frdWHYBSzAdPlBqlzPHt7kQBjqBCS8eDCqWGqUcxSo8liLSL0fQaymbsCXA25RHZGwX6E7WdgLpdPqlHC1GINzMh5PoHvMZALSBJHLUoIfY+HTw6BQSETMF+L7hqPljaK7IyC/Qjbz8CsNGKqUUImMOarMfUkpVxrMlzIACGnR/JMlLLMVyMaSdpeQRNQzAXrjkdLG0W2T8F+xO1HYC7X7mfKdepxylgp5PhgnrlynbFSjhNDeYYLIc0kpdyI7yi9FLOBAlxkjxTssi0bbWJae/tyi4AodhwfyrNUa4IzZioNZsoRDiOX8RkuhKqdi+wTBbvcYb0AX30sXZI4zo0VGcgFdy1DHM6HXJ2tMFOu89LtMseGWodLP3C8BEuQ8QzMqLaXVxazgWrnIvtAwS4rNlov3kxSoiRlsR5TbTQBuP946Y5liFNLda7MVLg8Veb6fJVqPWGxHvP6U4PMV5sEPtxaqOOAuJTiXAlQ7VxkPyjYZcXaXaPz1WilZLJUa3J1pkwpl6GRJEwt1cn4HnHSKr3cnK9RayQs1CKcaz1GCoS+Tz70GcgGOGcUsq3OjZ1aWikid1Owy4rlmvd8NWJyoQ606t/D+ZDJxTpXZ6s4Z4zM1zg+mCXwPI4N5hjOh4wUQrJ+wnA+JIodY4UsfgZGihmSxFHIBAwX0pXfBlRPF9k/CnZZsbz8cb7aOpN8eePSQj0i8H1ed3yQmXITvIT5WhMfj1qUEI+4Vt3cpbzpzAjjpdYsP+N5TJcb+L4xX4tWGoKpni6yvxTsctcF00IYML3UYL4atfrHeB7lqMlCuUnsUgIfJheq5IMM2YzPqZEC9x4bYKHa5NRwfqUrY6URU46+s5zRjDuae4nI/tB32RGy2YqX1W12V/ePiUl58pVZYpdQbzYZK4XgjNtpRCmfAXNY+3EK7QMxlmk5o0h3KNiPiM1WvKxusLVQi1ioNcmHPtPlGl9+aYa5akQhDJgoZhgutGrrlWbKieE8vgfnxooUs8G6u0W1nFHk4CnYj4iNOiSunlXXophrcw1evl3BkXJltsL12TrlRkzgGc2xAmdHixTzHqeHcxwrZSmErc1GG4W2ljOKHDwF+xGxXlkkitOVw6CbacpSPebVqSqVZhOXGpnAx/MNPEjSlOOlLKV8hjecHqYQBroQKtKjFOxHxNqySBSnvDS5xPXZGrOViGLOaDQdlWbEzdkaMSm++Zws5Sj4HrU4ppE4rs3WeMPpYQZyeumI9Cp9d/a5tRdMl/u5vDy5xPO3Fnl1qkzqUpoYcdMx32gwu9hguBAyUcpwz3iBgWCIr1+fY7iQYbHWpBrF3f5ricgm9hzsZuYDl4Drzrn3731I0ikb9XiZWqrz9I05Xp6q8sLNBRpximeOU4MFcoHPSC7HscEsI4WQ+yZKjA2E1OIY3/cJ/FirW0R6XCdm7B8GngUGO/BY0kErPV5qrROKalHC2dECX35piq+9usBCPWKpEeN7Hs3UsdRI8ZsJ5hvVRszp0TyjxZBzowNEccrV2SpjQWvTUhTreDqRXrWnYDezM8CPAr8F/LuOjEj2bLn84hwkiWO2XGehkRKnrbNGX7pVBnPtC6lGLvAxz2ekGHDvsQFODRXJZ40Hjg9yari1LPKesQEyvs9gPkMzSXXuqEgP2+uM/feBXwVKHRiLbNNmvdErjZippTpBexXM+ECWl6eWWKpFRFFCGMBMtUGSOgxHPuMxWMjQiBPOjuY5OzJAKR9wYii3EurQWiZZCH2a7dOOVI4R6V27DnYzez8w6Zx70sy+f5P7PQI8AnDu3LndPp20bbTRaPUZozPlBvceG2C+GlFrJgyEPs+XI6I4oZEkkEKapCvNu4aKWabKdYq5kNFShvsmSnetTd/JZqONfvCIyMHYy4z93cCPmdmPADlg0Mz+xDn3M6vv5Jx7FHgU4OLFi24PzyfcvdGotdnIWzljNB/6NNOUy5NlXri1gO95vDC5RK2RUMr7BOZzfiJPlCQM5jNEsSNOEwphQJLC1GKDh04MrRvI29lstNEPHhE5OLsOdufcR4CPALRn7P9+bahL5y33QJ9eagCOuB301+eqNJOUxWqThXoMpFy6Ms9wLsPzkwuM5UNSMhwbCBgrZqg2fU4OZTk2mOPWQp2vLs4xU65TixK+63y80shrpzba4SoiB0fr2A+B1aUNgGbsqDVjMGvV0g2WahFL9ZivXp6lEsWkOKbLEY1SlmbiGCxkyQQ+Z0fzXJgYwPdas+/7j5WYKNWZKtcZzIU04r2tUVfjL5Hu60iwO+e+AHyhE48ldwf56tJGMQyYqbQujk4vRuRDY64a8cLtJW4v1JlcrGOekaYpOd/IeR7jAznOjuY5PpTn/FiRwXzIeCm70kp3fCDH+bEB4gQCP2Qov7vZOqjxl0gv0Iy9x6zdVDReyt5xXF2tUme20qTRTLi2UCWOW+vUPc9IUiNKHHGc4DkYKrRKLeMDGd54ZpSJwSxh+4fF6hl1GHi848I41SimEAZ7bhegxl8i3aVg7zFrD46OkpSMb8yUG9ycr4Fz3Jivs1htUo+bLEUJM4ut5Yu1ZkIh8BgfynN6KEeYDXjLmRESB687UVo5v3T5eVbPqAdyew90EekN+k7uMRnfI0kc1UaTsH2htNZMmVqMuDpbZrEWk6RwY75KNYpJ4pS6g8SlDGYDStmQs6MFzo0VSWJHmjp8jIznrQR5FKfd/muKyD5SsPeYMPA4N1YkSlJmKg2SRkwcp+RCj2tzNV6+XQbAzBHgMVtv4HsejShhcCRgsdGkkA04O1rgvvESzTRlvtKkHMVUonjlhCQtRxTpX/qO7kEDuYB7xoqcGMxzaihPvZnyytQSt+dqVJoxM0t1bi1GzDUiZioRziB1EPhG1EzIBT5z5Yhqs1Uzz2dbR9ZFScp0uU6UpBSzAY5WSUZE+otm7D0q43vUophnbi7QaKbcmK1RbqbEzQQv8PCTFN+B7/tkMNLQqEUpCfD85BJztQiHcX6sSCYw5qsRkwt1RgZC5soRQOugai1HFOk7CvYDdNd69DUXMOfKEdOVOjhYqMfcnK/x3O0lQs944dYiC7WImWqDXOAxkMtSTRIKgU+UOk6Uitx7rEiaOgxjtJBlpJgBax1QvTwzX76AOlII72gboDYAIv1DwX5AVm+1j9shu9yo6/RIgUo95i+/eY1rczVuzVXJZgOaTcelK1PUmzG1CDCIY2j4KcNFR63pOD4U4jnj+x6aIDCf2MUkDo6VssTtFTXFbOufebkFQeh7d4W62gCI9A8F+wFZvdW+1Q4Ahgrf6fdyebrMy5NLTC40ePb2EmFgLJQbzFZb5ZUECIFirv2A5pP1YTyfIwjgzWeGOTtapJmkFDJBK6Hbz7cc0httHFIbAJH+omA/IKu32mf8VurOVyOWak2uTpd5bnKJZ24scGuxQaXexPdgqe6IVj2GAaEHnh9QyBjO+eQLPqXAZ6zY2oy0mY02DqkNgEh/UbB32Ea16uWt9pVGqxdLtdG6MPrizXmevlkhTVPKjRjfDPMczRhWL1gp+HBuNM99E0VODOU5OVrk2myVsyMFxkohQwW1ARCRFgV7B62to0+UcneUQgCmluos1po8dXWea7NlnnptjqlKg0YUEyfg+0YzhkYTfFrlF9+D8WKG150s8f43n2Fyoc5A3qdWzzJeyjKYDfY8y1YbAJH+oWDvoOVadRh4XJut0kwcGd+YKOXI+B5XZyr8w0uT1JoJr85UuLVQ45WZCuVGq4YOkMNhtN7PBpA1uGe8wNvuGeEdFyZ4/akhTg0X8Dy4Z7TEcDFDlKSYde/vLSK9RcHeQcu16oVqE3DkMz7X5ipUo4Qbc1WeuDzNE5dnKNdqVCMoR7C2SW4KZP3WbH18MEs24/EDDx1npJRjptzghduLnB0pcGq4sLKDVOvRRWQ1BXsHra6jZ3yj1kxopo5mkvLabIXLM0tcn65R2eQcqcCDUiFDLgsPjpfaHRoLFHMBI4UMmcBfKe+oLi4i61Gwd9DyhdNiNqCYHWC+EnFleokvXZ3l0iuzfOvmwoahngPCDIyWMhwfzPOW06OcGs0xEAaU6zFDXkCcpMSpo9yIieIqp0cKK2vURUSWKRU6ZO2F02IYcHlmiX+8Mss3Xp3l5kKV9vL1u2SB8aGAQpChVAh4+OQQFybyDOZz3DNe5NZijXtGCxTCgHIjZrgQar25iGxIwd4hqy+cPndjgWduLPD45RkuTy6y1Li7lg5Q9KGYM04OFblnvEjW98hmPL773jGcGWnqmCnXyXoeJ4cK7Za7Va03F5FNKdi3YaO16VGcrqxLXz5k+sp0mUuvzvGNK7M8N7lEZYMjREOgmPfJesZoMWRiIMuxoRwPjA8wWsphwMmhPBh3nGqkurqIbEXBvoUoTnl1uryydPH8+MDKYRUv3lri+nyNarPJQBiQpglffWWev3v+FldmG+vO0qHVK/nCWI7jwzlGclnedH6YKEo5UcoTO0ecOkq5gEoU39W3RevNRWQrCvYtVBoxtxZq5MKAevugijAIma9EPH97kevzFZ56ZZ4YR7ke0UgSppY2DnWA0SyMD2XJhz5DgyGNyFHI+twzXuT2Yuv4O9XRRWS3FOzbYu2eWq3/R3HK1dkKL0+WeeKVKZ65tkSawgbXRlcEwOgA3Ds2zLnRIoZxz2iBgWxAPhOwUIsIPY98JlAdXUR2TcG+hWI24MRQjmqUYCEsVCOmy3WqjYTAg+nFGrVNDiHy2//lszBWyOJ5kMkYmcAoZQNG8hkePDnEqaECzTSlELbWqK+uo6tXuojsxK6D3czOAv8TOA444FHn3Ec7NbBeMpwPWaxVuDFf4UsvTOOb48pchdcWKrw2v1nRBTIGp0ZCRgaynCrlSBLHAydLDOSy3H+syBtOD9/RG32ZeqWLyG7tZcYeA7/inPuamZWAJ83sMefcMx0aW9cth2o1Sri9VGe62mqxe3muwuRcmamKY7MTQ48VW0sZ33FhjDhNGcoHRAmcGikykAt4w+nhLVvtqle6iOzUroPdOXcTuNl+e8nMngVOA30R7FGcMl+NiJKUfMZnsRLx5MvTfPnFaRabbBroAXDfeJZ/ev8Eo4UsJ0dyzFdixgdCYucYCH3ODBcY3karXfVKF5Gd6kiN3czOA28DvrLOxx4BHgE4d+5cJ55u30Vxyou3F7kyU+bVmSqpg6evzvGll6Ypb155YTSE8VKOH3r9Ce49USIbeJweKXJ7oUE+6xElKSO5kPuOlbY181ZPGBHZqT0Hu5kNAH8G/JJzbnHtx51zjwKPAly8eHGT9lfdt7zhaLpc50svTfH4S1O8NlVmttqkHkF9k8/NAhNDHvdNDPHdD0zwXedGGS9leW2minNQCD1GiiEexrmx4sqGo+3Q2nUR2Yk9BbuZZWiF+iedc5/pzJC6o1yPeWlyiZlKg5tzVf7hxSm+eWWO+ebWnzsYwBvPDPGDDx/nXfdPcHKwwHwtwjk4MZRrlVxcnkzg3XXwhohIp+1lVYwBfwQ865z7vc4N6eBFccrVmQrX52q8NlPmK5eneeKluTvOG91IzoP7jw3y9gujvPf1Jzk/MQBAMRfQTFKiZsrNhRq+b4S+p26MIrLv9pIy7wZ+FviWmT3Vvu3XnXP/Z+/DOljNJGW6XOMrL0/zheeuc3Vxs0ujLR5QysDIQMjpkRwPTAwyXvrOCpflWfmV6QpT5TqFMGAwn9GqFhHZd3tZFfP3LG/FPKSWa+pXZ8r8ryeu8OXL08xvtX0UOFH0GcpnePjUEG84O8S54QHees/IXXXzZpLi+0Yhm6HaaFIM9342qYjIVo5UXWD1Dk6Ab7w2y9M3Fnj8+SmefG1my1A/MQD3Hx/mRCnPA8cHmBjMMpQPCX2fQnj3lzLje4S+x2AuoJjxOTdW1GxdRPbdkQn21V0acbBYb/CpJ67yzRszvDq7+RrGELj3WI5/9sAE900MMl4KOT8+QCNunZYUxesfJq2liiLSDUci2KM4ZWqpzmtzVQLPePzyNE+8PMXT1ypstehlKAPvffgkbz83wrnxArkwoFxr0mg65irRyqx8oxKLliqKyEHr+2BfbgswXW7w7I0Frs5WeezbN5jbTi29FPDjbz/Nv3j7OYbyIeVGjBlUo5hi1icT5BgphOv2ehER6Za+D/ZmklKJWhdI//752zx1dWnTjUYA41k4MVrke+6f4AcePMGp4QIZ36PSiIniFA8jSlJC31Ooi0jP6dtgL9djqlFMtRHzxOUpHnv6Jk9eXdr0AAyAB0ZC3nx2lNFSljedGWa4kFmpjy/Xy08PFzBDdXMR6Ul9GezleszXr8yyWIt5/OXbfOXybV6cSTb9nKIPb75niH/1ttM8dHqUoXyGQhis7BRVT3QROSz6MtirUcwr00t8+cVJvvjcLLVN7jtgcGaitcHonfeP8dCpUe5f06BLPdFF5DDpy2B/+tos//X/Psf0Jok+EsKFiQEeODHIQDbgDWeGuXd8YN215uqJLiKHSV8Fe7ke88TLk/ynz35z01A/N5Lhh99whgePF8llA0YKIfdNlBgurn8hVD3RReQw6YtgL9djbi5U+cunrvGnX32V25WNuwPfO+Lzs+++jx9902ky7RDfquOiNhqJyGFy6IO9XI/522dv8Ni3b/PFlyY4FqJ8AAAGnElEQVRZ2GB9+lgO3nZ2jA99zwXeem5sR/3QQRuNROTwONTBXq7HfOHZG3zi71/hmevldS+SnikZP/j6Uzx0cpCHTw/z8KlhBbSI9LVDGezLLQL++ls3+POvX+fFW2XWTtSLwFvvLfEL73k9906UyPg65EJEjoZDF+xz5Yinb8zz5Zem+bsXJlmoNwgz0GifipEF7j2e4/sfPMZP/JMLKwdfiIgcFYcq2CcX6nz2m9f4+qszfOPaErWoQS0C34PjBY+TowW+5/5jvPuBce6dKLWOpANtLhKRI+XQBHu5HvM3z9zg7565za35OuV6g/GBHKVsyj1jBd71wDinh0ucHc2TD33KjZhKI2ailGNqqa7NRSJyZByKYI/ilJsLVW4t1PE8DwsMM48gMM6PDfAT77iHc+MDJIljvJSl3IhXNhNVo1ibi0TkSOn5YF/ezj9fbXVOL+UDUnK88eQwD58Z5F33j3P/scE7TkaqtGfrBhTC4I73tblIRPpdzwf78nb+k8N53nxmmIeODzBfj3nw+BATpexKaWX1LHztZqIw0OYiETk6ej7Yl7fzR3HKhfEBJkq51m2btM1dG/TaXCQiR0nPB7u284uI7MyeUtLM3mdmz5vZS2b2a50a1FphoM1FIiLbteukNDMf+G/ADwMPAz9lZg93amAiIrI7e5kCvwN4yTl32TkXAX8KfKAzwxIRkd3aS7CfBl5b9f619m0iItJF+160NrNHzOySmV2ampra76cTETny9hLs14Gzq94/077tDs65R51zF51zFycmJvbwdCIish17CfZ/BB4wswtmFgI/CXy2M8MSEZHdMuc2PkZuy082+xHg9wEf+Jhz7re2uP8UcGXNzePA9K4H0V0ae3do7N2hsXfHOFB0zm275LGnYO8EM7vknLvY1UHsksbeHRp7d2js3bGbsWvHj4hIn1Gwi4j0mV4I9ke7PYA90Ni7Q2PvDo29O3Y89q7X2EVEpLN6YcYuIiId1LVgP6jOkPvBzM6a2f8zs2fM7Gkz+3C3x7QTZuab2dfN7C+7PZadMrNhM/u0mT1nZs+a2Xd3e0zbYWa/3H6tfNvMPmVmuW6PaTNm9jEzmzSzb6+6bdTMHjOzF9t/jnRzjBvZYOy/037NfNPM/tzMhrs5xo2sN/ZVH/sVM3NmNr7V43Ql2PugM2QM/Ipz7mHgncC/OWTj/zDwbLcHsUsfBf7aOfcQ8BYOwd/DzE4DvwhcdM69kda+j5/s7qi29HHgfWtu+zXg8865B4DPt9/vRR/n7rE/BrzROfdm4AXgIwc9qG36OHePHTM7C/xz4Op2HqRbM/ZD3RnSOXfTOfe19ttLtMLlUDRAM7MzwI8Cf9jtseyUmQ0B3wv8EYBzLnLOzXd3VNsWAHkzC4ACcKPL49mUc+6LwOyamz8AfKL99ieAHz/QQW3TemN3zn3OORe3332CVguUnrPB1x3gvwC/Cmzromi3gr1vOkOa2XngbcBXujuSbft9Wi+QtNsD2YULwBTwx+1S0h+aWbHbg9qKc+468Lu0Zls3gQXn3Oe6O6pdOe6cu9l++xZwvJuD2YOfA/6q24PYLjP7AHDdOfeN7X6OLp7ugZkNAH8G/JJzbrHb49mKmb0fmHTOPdntsexSALwd+O/OubcBFXq3HLCiXYv+AK0fTKeAopn9THdHtTeutZzu0C2pM7PfoFVK/WS3x7IdZlYAfh34jzv5vG4F+7Y6Q/YyM8vQCvVPOuc+0+3xbNO7gR8zs1dplb9+0Mz+pLtD2pFrwDXn3PJvR5+mFfS97r3AK865KedcE/gM8K4uj2k3bpvZSYD2n5NdHs+OmNmHgPcDP+0Ozzrv+2hNCL7R/r49A3zNzE5s9kndCvZD3RnSzIxWnfdZ59zvdXs82+Wc+4hz7oxz7jytr/nfOucOzczROXcLeM3MHmzf9B7gmS4OabuuAu80s0L7tfMeDsFF33V8Fvhg++0PAn/RxbHsiJm9j1YJ8secc9Vuj2e7nHPfcs4dc86db3/fXgPe3v5e2FBXgr19EePfAn9D6wX+v51zT3djLLv0buBnac14n2r/9yPdHtQR8QvAJ83sm8Bbgf/c5fFsqf0bxqeBrwHfovV919M7Ic3sU8DjwINmds3Mfh74beCHzOxFWr+F/HY3x7iRDcb+B0AJeKz9/fo/ujrIDWww9p0/zuH5jURERLZDF09FRPqMgl1EpM8o2EVE+oyCXUSkzyjYRUT6jIJdRKTPKNhFRPqMgl1EpM/8f1UFsxSikSP8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(r1, r2, '.', alpha=0.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.004518297204525764 0.06947544091649573\n"
     ]
    }
   ],
   "source": [
    "err = (r1-r2)/r1\n",
    "\n",
    "print(np.median(err), np.std(err))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in mass\n",
    "epos['mass'] = np.load(MAINPATH + 'catalogues/epos_observed_pred_masses_oct22.npy')[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('%d planets in OBSERVED MULTIS catalogue' %len(epos['star_id']))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/epos_observed_multis_with_masses.csv'\n",
    "\n",
    "    with open(filepath, \"w\") as outfile:\n",
    "        writer = csv.writer(outfile)\n",
    "        writer.writerow(epos.keys())\n",
    "        writer.writerows(zip(*epos.values()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in masses and combine for UNDETECTED MULTIS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# READ IN EPOS DATABASE\n",
    "epos_file = MAINPATH + 'catalogues/epos_undetected_planets_with_limit_cuts.csv'\n",
    "\n",
    "epos_data = []\n",
    "with open(epos_file) as infile:\n",
    "    reader = csv.reader(infile)\n",
    "    \n",
    "    for row in reader:\n",
    "        epos_data.append(row)\n",
    "        \n",
    "KEYS = epos_data[0]\n",
    "VALUES = epos_data[1:]\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PULL RELEVANT QUANTITIES FROM CATALOGUE\n",
    "epos = {}\n",
    "\n",
    "epos['star_id'] = np.array(getkey('star_id'), dtype='int')\n",
    "\n",
    "epos['period'] = getkey('period')\n",
    "epos['radius'] = getkey('radius')\n",
    "epos['tdur']   = getkey('tdur')\n",
    "\n",
    "epos['rhostar'] = np.ones_like(epos['period'])\n",
    "epos['rstar']   = np.ones_like(epos['period'])\n",
    "epos['mstar']   = np.ones_like(epos['period'])\n",
    "\n",
    "epos['npl'] = np.array(getkey('npl'), dtype='int')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove single-planet systems\n",
    "epos_npl = epos['npl'].copy()\n",
    "\n",
    "for k in epos.keys():\n",
    "    epos[k] = epos[k][epos_npl > 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in mass\n",
    "epos['mass'] = np.load(MAINPATH + 'catalogues/epos_undetected_pred_masses_dec11.npy')[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('%d planets in UNDETECTED MULTIS catalogue' %len(epos['star_id']))\n",
    "\n",
    "WRITENEW = True\n",
    "if WRITENEW:\n",
    "    filepath = MAINPATH + 'catalogues/epos_undetected_multis_with_masses.csv'\n",
    "\n",
    "    with open(filepath, \"w\") as outfile:\n",
    "        writer = csv.writer(outfile)\n",
    "        writer.writerow(epos.keys())\n",
    "        writer.writerows(zip(*epos.values()))"
   ]
  }
 ],
 "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
}