galacticHeight.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pylab as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import astropy.units as u\n",
    "from os import path\n",
    "\n",
    "import astropy.coordinates as coord\n",
    "import astropy.units as u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "prefix = '1to13'\n",
    "path = '/astro/store/gradscratch/tmp/scw7/tessData/lightcurves/sec1to13/'\n",
    "log_path = path + 'log/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7495"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(prefix + '_filt_flare_out.csv')\n",
    "df_param = pd.read_csv(prefix + '_filt_param_out.csv')\n",
    "len(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7495 5747 118887\n"
     ]
    }
   ],
   "source": [
    "# Match flare detections to gaia\n",
    "tic_gaia_table = pd.read_csv('ticGAIA1to13.csv')\n",
    "df_g = pd.merge(tic_gaia_table, df, on='TIC', how='inner')\n",
    "df_g = df_g.dropna(subset=['lum'])\n",
    "ed = (df_g['ed'].values*u.day).to(u.s)\n",
    "ed_err = (df_g['ed_err'].values*u.day).to(u.s)\n",
    "energy = (ed*df_g['lum']*u.erg/u.s).value\n",
    "energy_err = (ed_err*df_g['lum']*u.erg/u.s).value\n",
    "df_g['energy'] = energy\n",
    "df_g['energy_err'] = energy_err\n",
    "print(len(df), len(df_g), len(tic_gaia_table))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['Unnamed: 0', 'TIC', 'r_est', 'lum', 'G_BPbrmag', 'G_RPmag', 'G_mag',\n",
       "       'source_id', 'ra', 'dec'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tic_gaia_table.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_fl = coord.SkyCoord(ra=df_g['ra'].values*u.deg, dec=df_g['dec'].values*u.deg, \\\n",
    "                      distance=df_g['r_est'].values*u.pc)\n",
    "c_all = coord.SkyCoord(ra=tic_gaia_table['ra'].values*u.deg, dec=tic_gaia_table['dec'].values*u.deg, \\\n",
    "                      distance=tic_gaia_table['r_est'].values*u.pc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Only care about z, distance from galactic center doesnt matter\n",
    "galcen = c_fl.transform_to(coord.Galactocentric(z_sun=0*u.pc, galcen_distance=8.1*u.kpc))\n",
    "df_g['gal_z'] = np.fabs(galcen.z)\n",
    "galcen = c_all.transform_to(coord.Galactocentric(z_sun=0*u.pc, galcen_distance=8.1*u.kpc))\n",
    "tic_gaia_table['gal_z'] = np.fabs(galcen.z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fab362a4d30>"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7YAAAHnCAYAAACMt8HjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs/XuUnnV9N/q/PyQIKhILpK0mxkSDaGww6Ej0qVqwclJDHpQKxFWBhzb2wO+n7e5qoa2KtLg90NLdLW1JlUNVEpCnaJAgtAVKdSEQMDUQyk6kUQJBIthELKfAd/8xk+whTGBymNxzDa/XWllz39/rcL+vTFauec91qtZaAAAAoKt263UAAAAA2BGKLQAAAJ2m2AIAANBpii0AAACdptgCAADQaYotAAAAnabYAgAA0GmKLQAAAJ2m2AIAANBp43sdYEfst99+berUqb2OAcAYceutt/64tTax1zm6zL4ZgJ1puPvmThfbqVOnZunSpb2OAcAYUVU/6HWGrrNvBmBnGu6+2anIAAAAdJpiCwAAQKcptgAAAHRap6+xBdhVnnjiiaxZsyaPPvpor6OwE+y5556ZPHlydt99915HAWCM8rPDttnRfbNiCzAMa9asyUte8pJMnTo1VdXrOOyA1loefPDBrFmzJtOmTet1HADGKD87DN/O2Dc7FRlgGB599NHsu+++dkxjQFVl33339Rt0AEaUnx2Gb2fsmxVbgGGyYxo7fC8B2BXsb4ZvR/+uFFsAAAA6TbEF6JDLL788VZX/+I//2Dy2evXq/NIv/VKS5Prrr8973/vepy1z9dVXZ9asWZk1a1b22muvHHDAAZk1a1Y+9KEP5frrr8+ECRM2T581a1b++Z//OUly1lln5fWvf30OPPDAzJo1KzfddFOS5Bvf+EYOOuigvOENb8iMGTNy3nnnDZn14Ycfzm//9m/n1a9+dQ466KC86U1vyt///d8/bZ5zzjkne+65Z9avX795bKhtmDt3bt761rdu598aADw/jRs37mn7+NWrVw+5n91RH/7wh/Ptb387SbJx48bst99+Of300582zyGHHJKlS5cmSaZOnZof//jHOzWDm0cBdMjChQvztre9LYsWLcoZZ5wxrGWOOOKIHHHEEUn6dypnn312+vr6kvSXyLe//e35xje+8bRlbrzxxnzjG9/Ibbfdlj322CM//vGP8/jjj+eJJ57I/Pnzc/PNN2fy5Ml57LHHsnr16iE/9zd+4zfyqle9KitXrsxuu+2WdevW5fzzz3/G9rz5zW/O5ZdfnpNOOmnI9fzXf/1Xbrvttuy11175z//8Tzd8AoBheuELX5hly5Y9bWxr++2t2bhxY8aPf/baeNNNN+Vv/uZvkiTXXHNNDjjggFx66aX51Kc+tctOx1ZsAbbRJ6+4Iyvu27BT1znj5XvnE3Ne/6zzPPzww/n2t7+d6667LkcfffSwi+32WLt2bfbbb7/sscceSZL99tsvSfLQQw9l48aN2XfffZMke+yxRw444IBnLP/9738/N998cy6++OLstlv/yUETJ07MH/3RHz1tnocffjif+9zn8qlPfWqrxfZ//+//nTlz5uQXfuEXsmjRomf8BhgARrte/ezwXG6++eZ89KMfzSOPPJIXvvCFueCCC3LAAQfkwgsvzJVXXplHH300P/vZz3Lttdfmc5/7XC699NI89thjOeaYY/LJT34ySXLnnXfmNa95TcaNG5ek/5fWH/nIR/K3f/u3+c53vrPLzrhyKjJAR3zta1/LkUcemde85jXZZ599ctttt+2U9f7bv/3b005T+v73v5/DDz8899xzT17zmtfkd37nd/Kv//qvSZJ99tknRx99dF75ylfmhBNOyFe+8pU89dRTz1jnHXfckTe84Q2bS+1QFi5cmBNOOCFvf/vbc9ddd+WBBx541vlOOOGELFy4cKdsMwA8HzzyyCOb9+/HHHPMM6a/9rWvzQ033JDvfve7OfPMM/PHf/zHm6fdeOONueiii3LttdfmmmuuycqVK3PzzTdn2bJlufXWW3PDDTckSa666qoceeSRmz/vX/7lX/Le9753l++3HbEF2EY7+tvR7bVw4cJ89KMfTZIcf/zxWbhwYd74xjfu8HqHOhU5SW699db827/9W6677rocd9xx+fSnP52TTjopX/jCF7J8+fL88z//c84+++z80z/9Uy688MJn/YyzzjorX/3qV/PAAw/kvvvuS5IsWrQol19+eXbbbbe8733vy1e/+tX87u/+7tOW+9GPfpRVq1blbW97W6oq48ePz+233775mmIA6IJe/eww1KnIg61fvz4nnnhiVq5cmarKE088sXnaYYcdln322SdJ/+nF11xzTQ466KAk/WeRrVy5Mu94xzty9dVX54ILLkjSfx+OQw89NC960Yvy/ve/P3/2Z3+Wc845Z/PR3JGk2AJ0wIMPPphrr702t99+e6oqTz75ZKoqn/3sZ0fsM8eNG5dDDjkkhxxySGbOnJmLLrpo8+nCM2fOzMyZM/Prv/7rmTZtWr74xS/mTW96U5Lk6KOPzoc+9KH8+7//e5566qnstttu+ZM/+ZP8yZ/8Sfbaa68kyfe+972sXLkyhx12WJLk8ccfz6te9apnFNtLLrkkP/nJTzZfV7thw4YsWrQof/7nfz5i2w0Azxcf+9jHcuihh+byyy/P6tWrc8ghh2ye9uIXv3jz69ZaTj/99Hz4wx9+2vL//d//nf/6r//Ky1/+8iT9v4T/9re/nalTpybp//nluuuuy7ve9a4R3xanIgN0wGWXXZYPfehD+cEPfpDVq1fnnnvuybRp0/Ktb31rRD7vrrvuysqVKze/X7ZsWV75ylfm4YcfzvXXX/+M8XHjxmXZsmVZtmxZzjzzzEyfPj19fX350z/90zz55JNJ+h9U31pL0r/jO+OMM7J69eqsXr069913X+6999784Ac/eFqOhQsX5pvf/Obm+W699dYsWrRoRLYZAJ5v1q9fn0mTJiXJs559dcQRR+T888/Pww8/nCS5995788ADD+S6667LoYcemqT/l8/f+ta38sMf/nDzfvvcc8/dZacjd7LYVtWcqlow+PEQAGPZwoULn3FtzPvf//5cfPHFO7zuLa+xveyyy/Lwww/nxBNPzIwZM3LggQdmxYoVOeOMM9Jay2c/+9nNjwz6xCc+sdUd4Re+8IU8+OCDmT59et70pjflXe96Vz7zmc8k6T8NecvtOeaYY55WWlevXp0f/vCHectb3rJ5bNq0adl77703P3oIANh+f/iHf5jTTz89v/zLv7z5F9FDOfzwwzNv3ry89a1vzcyZM3Psscfmpz/96dOur/3Hf/zHvPOd79x848mk/3F9ixcvzmOPPTbi21KbfnveRX19fW3Ts5AARtKdd96Z173udb2OwU401Pe0qm5trfX1KNKYYN8M0O/58LPDG9/4xtx0003Zfffdd8r6dmTf/Ly/xnYkbr29pbmzJmXe7Ckj+hkAwPBcfNMP8/Vl9/Y6xnbxMwUwmuysJzTsDJ08FblLVqzd0NmdJwBjQ1UdWVV3VdWqqjptiOnvqKrbqmpjVR07xPS9q+reqvr8rkk8sr6+7N6sWDuyv9QeCX6mANi65/0R25G+9fZx5904ousHgGdTVeOSnJvksCRrktxSVYtbaysGzfbDJCcl+YOtrObPkvzrSObc1Wa8bO9c8uG39jrGNvEzBXRPay1V1esYnbCjl8g+74stQ3OaFsCYcXCSVa21u5OkqhYlmZtkc7Ftra0emPbUlgtX1ZuS/EKSbyZx/THAMO2555558MEHs++++yq3z6G1lgcffDB77rnndq9Dsd0FVqzd0Lnfst70nw8lSWZP26fHSbbNplPLFFuAzSYluWfQ+zVJZg9nwaraLclfJPn1JL/6LPPNTzI/SaZM8f8vQJJMnjw5a9asybp163odpRP23HPPTJ48ebuXV2xH2NxZk3odYbvMnrZPJ498du0XCAC7wFCHCYZ7vtfvJFnSWrvn2Y42tNYWJFmQ9N8VeZsTAoxBu+++e6ZNm9brGM8biu0Imzd7SufKITA6jRs3LjNnztz8/mtf+1pWr16ds88+O9/4xjd26mfNnj07jz32WB566KE88sgjmx/e/rWvfS2HHHJIXvKSl2TcuHFJkne84x3567/+63znO9/JRz7ykTz22GN57LHHctxxx+WMM87Ij370o5xyyim555578sQTT2Tq1KlZsmTJkJ/75S9/OZ/97Gfz5JNPZvz48Xnzm9+cs88+Oy996UuTJOvWrcvLX/7yfP7zn8+HP/zhzctNnTo1S5cuzX777Zckufzyy/O+970vd955Z1772tfu1L+bDlqT5BWD3k9Oct8wl31rkrdX1e8k2SvJC6rq4dbaM25ABQC9pNgCdMQLX/jCLFu27Gljq1ev3qZ1bNy4MePHP/d//TfddFOS5MILL8zSpUvz+c8//Wa411133eYSucmJJ56YSy+9NG94wxvy5JNP5q677kqSfPzjH89hhx2Wj3zkI0mS733ve0N+5je/+c2cc845ueqqqzJp0qQ8+eSTueiii/KjH/1oc7H96le/mre85S1ZuHDh04rtlhYuXJi3ve1tWbRoUc4444zn3N4x7pYk+1fVtCT3Jjk+ybzhLNha++Cm11V1UpI+pRaA0UixBdhWV52W3L98567zF2cmR316h1Zx880356Mf/WgeeeSRvPCFL8wFF1yQAw44IBdeeGGuvPLKPProo/nZz36Wa6+9Np/73Ody6aWX5rHHHssxxxyTT37ykzu8CQ888EBe9rKXJek/ujxjxowkydq1a3P44Ydvnu/AAw8ccvmzzjorZ5999uajw+PGjcv/+l//62nzLFy4MH/xF3+RefPm5d57790872APP/xwvv3tb+e6667L0Ucf/bwvtq21jVV1apKrk4xLcn5r7Y6qOjPJ0tba4qp6c5LLk/xckjlV9cnW2sg+NgAAdiLPsQXoiEceeSSzZs3KrFmzcswxxzxj+mtf+9rccMMN+e53v5szzzwzf/zHf7x52o033piLLroo1157ba655pqsXLkyN998c5YtW5Zbb701N9xwwzZlOfTQQzdnOeecc5Ikv/d7v5cDDjggxxxzTM4777w8+uijSZLf/d3fzSmnnJJDDz00Z511Vu67b+izYO+444688Y1v3Opn3nPPPbn//vtz8MEH5wMf+EAuueSSIef72te+liOPPDKvec1rss8++4yqh8f3SmttSWvtNa21V7fWzhoY+3hrbfHA61taa5Nbay9ure07VKltrV3YWjt1V2cHgOFwxBZgW+3gkdXtNdSpyIOtX78+J554YlauXJmqyhNPPLF52mGHHZZ99um/y/k111yTa665JgcddFCS/iOcK1euzDve8Y5hZxnqVOSPf/zj+eAHP5hrrrkmF198cRYuXJjrr78+RxxxRO6+++5885vfzFVXXZWDDjoot99+eyZOnLjV9S9fvjy//uu/np/+9Kf51Kc+leOOOy6LFi3KBz7wgSTJ8ccfn1NOOSW///u//4xlFy5cmI9+9KOb51u4cOGzFmYAoPsUW4Ax4mMf+1gOPfTQXH755Vm9enUOOeSQzdNe/OIXb37dWsvpp5/+rNeobq9Xv/rV+e3f/u385m/+ZiZOnLj5+X377LNP5s2bl3nz5uW9731vbrjhhtx222258sorkyTLli3L61//+tx222059NBDM3PmzCxbtiynnnpqHnnkkST9hfVHP/pRvvKVryRJ7rvvvqxcuTL777//5s9/8MEHc+211+b2229PVeXJJ59MVeWzn/2sZwgCwBim2DLmdPG5wUk6+XglRpf169dvvub0wgsv3Op8RxxxRD72sY/lgx/8YPbaa6/ce++92X333fPzP//zO/T5V155Zd797nenqrJy5cqMGzcuL33pS3PttdfmLW95S170ohflpz/9ab7//e9nypQpef/735+zzjpr8/Knn356/uAP/iBf//rXNz/HblOpveuuu/Kzn/0s99577+b5P/GJT2TRokX52Mc+tnnssssuy4c+9KGcd955m8d+5Vd+Jd/61rfy9re/fYe2DwAYvRRbxpSuPjd4xdoNSaLYskP+8A//MCeeeGL+8i//Mu985zu3Ot/hhx+eO++8M29961uTJHvttVe+/OUvb1OxPfTQQzc/7ufAAw/MP/zDP+RLX/pSfu/3fi8vetGLMn78+HzlK1/JuHHjcuutt+bUU0/N+PHj89RTT+U3fuM38uY3v/kZ63z3u9+ddevW5aijjsqTTz6Zl770pfmlX/qlHHHEETnvvPOecV3x+9///hx//PFPK7YLFy7Maaed9oz5Lr74YsUWAMawaq27z1Hv6+trS5cu7XUM2GGbjjBf8uG39jgJW3PnnXfmda97Xa9jsBMN9T2tqltba309ijQmdGHf3NX/c7uaG2BHDHff7K7IAAAAdJpiCwAAQKcptgDD1OVLN3g630sAGFsUW4Bh2HPPPfPggw8qRGNAay0PPvhg9txzz15HAQB2EndFBhiGyZMnZ82aNVm3bl2vo7AT7LnnnpsfKQQAdJ9iCzAMu+++e6ZNm9brGAAADMGpyAAAAHSaYgsAAECnKbYAAAB0mmILAABApym2AAAAdJq7IsMosWLthhx33o29jrHN5s6alHmzp/Q6BgAAz2OKLYwCc2dN6nWE7bJi7YYkUWwBAOgpxRZGgXmzp3SyHHbxCDMAAGOPa2wBAADotFFVbKvqxVV1a1W9t9dZAAAA6IYRLbZVdX5VPVBVt28xfmRV3VVVq6rqtEGT/ijJpSOZCQAAgLFlpI/YXpjkyMEDVTUuyblJjkoyI8kJVTWjqt6VZEWSH41wJgAAAMaQEb15VGvthqqausXwwUlWtdbuTpKqWpRkbpK9krw4/WX3kapa0lp7ast1VtX8JPOTZMqU7t1sBwAAgJ2rF3dFnpTknkHv1ySZ3Vo7NUmq6qQkPx6q1CZJa21BkgVJ0tfX10Y2KgAAAKNdL4ptDTG2uaC21i7cdVEAAADoul7cFXlNklcMej85yX09yAEAAMAY0Itie0uS/atqWlW9IMnxSRb3IAcAAABjwEg/7mdhkhuTHFBVa6rqlNbaxiSnJrk6yZ1JLm2t3TGSOQAAABi7RvquyCdsZXxJkiXbu96qmpNkzvTp07d3FQAAAIwRvbh51A5rrV2R5Iq+vr7f7HUWeL5bsXZDjjvvxl7H2GZzZ03KvNkeGQYAMBZ0stgCo8PcWZN6HWG7rFi7IUkUWwCAMUKxBbbbvNlTOlkOu3iEGQCArevFXZEBAABgp1FsAQAA6LROFtuqmlNVC9avX9/rKAAAAPRYJ4tta+2K1tr8CRMm9DoKAAAAPdbJYgsAAACbKLYAAAB0mmILAABApym2AAAAdJpiCwAAQKd1sth63A8AAACbdLLYetwPAAAAm3Sy2AIAAMAmii0AAACdptgCAADQaeN7HQCgF1as3ZDjzrux1zG22dxZkzJv9pRex6BDqurIJP9XknFJvtBa+/QW09+R5K+SHJjk+NbaZQPjs5L8bZK9kzyZ5KzW2iW7MjsADJdiCzzvzJ01qdcRtsuKtRuSRLFl2KpqXJJzkxyWZE2SW6pqcWttxaDZfpjkpCR/sMXi/53kQ621lVX18iS3VtXVrbX/2gXRAWCbKLbA88682VM6WQ67eISZnjs4yarW2t1JUlWLksxNsrnYttZWD0x7avCCrbX/Z9Dr+6rqgSQTkyi2AIw6nbzG1nNsAWBYJiW5Z9D7NQNj26SqDk7ygiTf30m5AGCn6mSx9RxbABiWGmKsbdMKql6W5EtJTm6tPbWVeeZX1dKqWrpu3brtiAkAO8apyAAwdq1J8opB7ycnuW+4C1fV3kmuTPKnrbXvbG2+1tqCJAuSpK+vb5uK81A+ecUdWXHfhh1dzVatWLshM16294itH4Bdr5NHbAGAYbklyf5VNa2qXpDk+CSLh7PgwPyXJ/mH1tpXRzDjLjfjZXt39iZyAAzNEVsAGKNaaxur6tQkV6f/cT/nt9buqKozkyxtrS2uqjenv8D+XJI5VfXJ1trrk3wgyTuS7FtVJw2s8qTW2rKRzv2JOa8f6Y8AYIxRbAFgDGutLUmyZIuxjw96fUv6T1HecrkvJ/nyiAcEgJ3AqcgAAAB0mmILAABApym2AAAAdFoni21VzamqBevXr+91FAAAAHqsk8W2tXZFa23+hAkTeh0FAACAHutksQUAAIBNFFsAAAA6TbEFAACg0xRbAAAAOk2xBQAAoNMUWwAAADpNsQUAAKDTFFsAAAA6bXyvAwAwfCvWbshx593Y6xjbbO6sSZk3e0qvYwAAY5RiC9ARc2dN6nWE7bJi7YYkUWwBgBHTyWJbVXOSzJk+fXqvowDsMvNmT+lkOeziEWYAoFs6eY1ta+2K1tr8CRMm9DoKAAAAPdbJYgsAAACbKLYAAAB0mmILAABApym2AAAAdJpiCwAAQKcptgAAAHSaYgsAAECnKbYAAAB0mmILAABApym2AAAAdJpiCwAAQKcptgAAAHSaYgsAAECndbLYVtWcqlqwfv36XkcBAACgx8b3OsD2aK1dkeSKvr6+3+x1FgCe24q1G3LceTeOyLpnvHzvfGLO60dk3QBAN3Sy2ALQHXNnTep1BABgjFNsARhR82ZPybzZU3odAwAYwzp5jS0AAABsotgCAADQaYotAAAAnabYAgAA0GmKLQAAAJ3mrsgjbekFyfLLep1i+8w8Nuk7udcpAAAAnpUjtiNt+WXJ/ct7nWLb3b+8u4UcAAB4XnHEdlf4xZnJyVf2OsW2ueA9vU4AAAAwLI7YAgAA0GmKLQAAAJ3mVOSrThvZa2DvX95/KjIAAAAjwhHbkfaLM/vvLgwAAMCIcMT2qE/3OgEAAAA7wBFbAAAAOk2xBYAxrqqOrKq7qmpVVZ02xPR3VNVtVbWxqo7dYtqJVbVy4M+Juy41AAyfYgsAY1hVjUtybpKjksxIckJVzdhith8mOSnJxVssu0+STySZneTgJJ+oqp8b6cwAsK0UWwAY2w5Osqq1dndr7fEki5LMHTxDa211a+17SZ7aYtkjkvxTa+2h1tpPkvxTkiN3RWgA2BadLLZVNaeqFqxfv77XUQBgtJuU5J5B79cMjI30sgCwy3Sy2LbWrmitzZ8wYUKvowDAaFdDjLWduWxVza+qpVW1dN26ddsUDgB2hk4WWwBg2NYkecWg95OT3Lczl22tLWit9bXW+iZOnLjdQQFgeym2ADC23ZJk/6qaVlUvSHJ8ksXDXPbqJIdX1c8N3DTq8IExABhVFFsAGMNaaxuTnJr+Qnpnkktba3dU1ZlVdXSSVNWbq2pNkl9Lcl5V3TGw7ENJ/iz95fiWJGcOjAHAqDK+1wEAgJHVWluSZMkWYx8f9PqW9J9mPNSy5yc5f0QDAsAOcsQWAACATlNsAQAA6DTFFgAAgE5TbAEAAOg0xRYAAIBOU2wBAADoNMUWAACATlNsAQAA6DTFFgAAgE5TbAEAAOg0xRYAAIBOU2wBAADoNMUWAACATlNsAQAA6DTFFgAAgE5TbAEAAOg0xRYAAIBOU2wBAADoNMUWAACATlNsAQAA6DTFFgAAgE5TbAEAAOg0xRYAAIBOU2wBAADoNMUWAACATlNsAQAA6LRRU2yr6nVV9XdVdVlV/Xav8wAAANANI1psq+r8qnqgqm7fYvzIqrqrqlZV1WlJ0lq7s7X2W0k+kKRvJHMBAAAwdoz0EdsLkxw5eKCqxiU5N8lRSWYkOaGqZgxMOzrJt5L8ywjnAgAAYIwY0WLbWrshyUNbDB+cZFVr7e7W2uNJFiWZOzD/4tba/0jywa2ts6rmV9XSqlq6bt26kYoOAABAR4zvwWdOSnLPoPdrksyuqkOSvC/JHkmWbG3h1tqCJAuSpK+vr41cTAAAALqgF8W2hhhrrbXrk1y/a6MAAADQdb24K/KaJK8Y9H5ykvt6kAMAAIAxoBfF9pYk+1fVtKp6QZLjkyzuQQ4AAADGgJF+3M/CJDcmOaCq1lTVKa21jUlOTXJ1kjuTXNpau2MkcwAAADB2jeg1tq21E7YyviTPcoOo51JVc5LMmT59+vauAgAAgDGiF6ci77DW2hWttfkTJkzodRQAAAB6rJPFFgAAADZRbAEAAOg0xRYAAIBOU2wBAADotE4W26qaU1UL1q9f3+soAAAA9Fgni627IgMAALBJJ4stAAAAbKLYAgAA0GmKLQAAAJ2m2AIAANBpii0AAACd1sli63E/ADB8VXVkVd1VVauq6rQhpu9RVZcMTL+pqqYOjO9eVRdV1fKqurOqTt/V2QFgODpZbD3uBwCGp6rGJTk3yVFJZiQ5oapmbDHbKUl+0lqbnuScJJ8ZGP+1JHu01mYmeVOSD28qvQAwmnSy2AIAw3ZwklWttbtba48nWZRk7hbzzE1y0cDry5L8alVVkpbkxVU1PskLkzyeZMOuiQ0Aw6fYAsDYNinJPYPerxkYG3Ke1trGJOuT7Jv+kvuzJGuT/DDJ2a21h0Y6MABsK8UWAMa2GmKsDXOeg5M8meTlSaYl+T+q6lXP+ICq+VW1tKqWrlu3bkfzAsA2U2wBYGxbk+QVg95PTnLf1uYZOO14QpKHksxL8s3W2hOttQeSfDtJ35Yf0Fpb0Frra631TZw4cQQ2AQCe3fheB2AUu395csF7ep1i2808Nuk7udcpAEaLW5LsX1XTktyb5Pj0F9bBFic5McmNSY5Ncm1rrVXVD5O8s6q+nORFSd6S5K92WXIAGCbFlqHNPLbXCbbP/cv7vyq2AEn6r5mtqlOTXJ1kXJLzW2t3VNWZSZa21hYn+WKSL1XVqvQfqT1+YPFzk1yQ5Pb0n658QWvte7t8IwDgOXSy2FbVnCRzpk+f3usoY1ffyd0sh108wgwwwlprS5Is2WLs44NeP5r+R/tsudzDQ40DwGjTyWtsPccWAACATTpZbAEAAGATxRYAAIBOU2wBAADoNMUWAACATlNsAQAA6DTFFgAAgE4bVrGtqrdX1bgtxt44MpEAgC3ZFwPA1g33iO3VSa6tql8YNPaFEcgzLFU1p6oWrF+/vlcRAGA+1cstAAAaxklEQVRXG1X7YgAYTYZbbO9K8rkk11fV/xgYq5GJ9Nxaa1e01uZPmDChVxEAYFcbVftiABhNxg9zvtZa+0ZV3ZXkkqo6P0kbwVwAwNPZFwPAVgz3iG0lSWttZZK3J3lHkgNHKhQA8Az2xQCwFcM6YttaO2jQ658l+UBVTRmxVADA09gXA8DWPWuxrar/O89+mtP/f+fGAQAGsy8GgOf2XKciL01ya5I9k7wxycqBP7OSPDmy0QCA2BcDwHN61iO2rbWLkqSqTkpyaGvtiYH3f5fkmhFPBwDPc/bFAPDchnvzqJcnecmg93sNjAEAu4Z9MQBsxXAf9/PpJN+tqusG3v9KkjNGJBEAMBT7YgDYiue6edT41trG1toFVXVVktkDk05rrd0/8vEA4PnNvhgAnttzHbH9TlWtSfLNJN9srX19F2QCAP4/9sUA8Bye6+ZRfVX1yiRHJfmrqpqU5FtJrkryr621x3ZBxmeoqjlJ5kyfPr0XHw8Au8xo3RcDwGjynDePaq39oLX2d621/5nkfyS5Ism7kvxbVV050gG3kumK1tr8CRMm9OLjAWCXGo37YgAYTYZ786gkycAjBq4d+JOB3xoDALuIfTEAPNOwim1V/XL677z4ysHLtNZeNTKxAIDB7IsBYOuGe8T2i0l+L8mtSZ4cuTgAwFbYFwPAVgy32K5vrV01okkAgGdjXwwAWzHcYntdVX0uyT8m2Xz3xdbabSOSCgDYkn0xAGzFcIvtpofBv2ngayVpSd650xMBAEOxLwaArXjWYltVvz/w8hvp33nWoMltpEIBAP3siwHguT3XEduXDHw9IMmbk3w9/TvUOUluGMFcAEA/+2IAeA7PWmxba59Mkqq6JskbW2s/HXh/RpKvjng6AHiesy8GgOe22zDnm5Lk8UHvH08ydaenAQC2xr4YALZiuDeP+lKSm6vq8vRfz3NMkotGLBUAsCX7YgDYimEV29baWVV1VZK3Dwyd3Fr77sjFAgAGsy8GgK0b7hHbTc/J86w8AOgR+2IAGNpwr7EdVapqTlUtWL9+fa+jAAAA0GOdLLattStaa/MnTJjQ6ygAAAD0WCeLLQAAAGyi2AIAANBpii0AAACdptgCAADQaYotAAAAnabYAgAA0Gnjex0Adrr7lycXvKfXKbbdzGOTvpN7nQIAADpHsWVsmXlsrxNsn/uX939VbAEAYJsptowtfSd3sxx28QgzAACMEq6xBQAAoNMUWwAAADpNsQWAMa6qjqyqu6pqVVWdNsT0ParqkoHpN1XV1EHTDqyqG6vqjqpaXlV77srsADAcii0AjGFVNS7JuUmOSjIjyQlVNWOL2U5J8pPW2vQk5yT5zMCy45N8OclvtdZen+SQJE/sougAMGyKLQCMbQcnWdVau7u19niSRUnmbjHP3CQXDby+LMmvVlUlOTzJ91pr/54krbUHW2tP7qLcADBsii0AjG2Tktwz6P2agbEh52mtbUyyPsm+SV6TpFXV1VV1W1X94S7ICwDbzON+AGBsqyHG2jDnGZ/kbUnenOS/k/xLVd3aWvuXpy1cNT/J/CSZMmXKDgcGgG3liC0AjG1rkrxi0PvJSe7b2jwD19VOSPLQwPi/ttZ+3Fr77yRLkrxxyw9orS1orfW11vomTpw4ApsAAM9OsQWAse2WJPtX1bSqekGS45Ms3mKexUlOHHh9bJJrW2stydVJDqyqFw0U3l9JsmIX5QaAYXMqMgCMYa21jVV1avpL6rgk57fW7qiqM5Msba0tTvLFJF+qqlXpP1J7/MCyP6mqv0x/OW5JlrTWruzJhgDAs+hksa2qOUnmTJ8+vddRAGDUa60tSf9pxIPHPj7o9aNJfm0ry345/Y/8AYBRq5OnIrfWrmitzZ8wYUKvowAAANBjnSy2AAAAsIliCwAAQKcptgAAAHSaYgsAAECnKbYAAAB0mmILAABApym2AAAAdJpiCwAAQKcptgAAAHSaYgsAAECnKbYAAAB0mmILAABApym2AAAAdJpiCwAAQKcptgAAAHSaYgsAAECnKbYAAAB0mmILAABApym2AAAAdJpiCwAAQKcptgAAAHSaYgsAAECnKbYAAAB0mmILAABApym2AAAAdJpiCwAAQKcptgAAAHSaYgsAAECnKbYAAAB0mmILAABAp42qYltV/7Oq/r6qvl5Vh/c6DwAAAKPfiBfbqjq/qh6oqtu3GD+yqu6qqlVVdVqStNa+1lr7zSQnJTlupLMBAADQfbviiO2FSY4cPFBV45Kcm+SoJDOSnFBVMwbN8qcD0wEAAOBZjXixba3dkOShLYYPTrKqtXZ3a+3xJIuSzK1+n0lyVWvttpHOBgAAQPf16hrbSUnuGfR+zcDY/y/Ju5IcW1W/NdSCVTW/qpZW1dJ169aNfFIAAABGtfE9+twaYqy11v46yV8/24KttQVJFiRJX19fG4FsAAAAdEivjtiuSfKKQe8nJ7mvR1kAAADosF4V21uS7F9V06rqBUmOT7K4R1kAAADosF3xuJ+FSW5MckBVramqU1prG5OcmuTqJHcmubS1dsdIZwEAAGDsGfFrbFtrJ2xlfEmSJduzzqqak2TO9OnTdyQaAAAAY0CvTkXeIa21K1pr8ydMmNDrKAAAAPRYJ4stAAAAbKLYAgAA0Gm9eo4tsKX7lycXvKfXKbbdzGOTvpN7nQIAgOcxxRZGg5nH9jrB9rl/ef9XxRYAgB7qZLF1V2TGnL6Tu1kOu3iEGQCAMaeT19i6KzIAAACbdLLYAgDDV1VHVtVdVbWqqk4bYvoeVXXJwPSbqmrqFtOnVNXDVfUHuyozAGwLxRYAxrCqGpfk3CRHJZmR5ISqmrHFbKck+UlrbXqSc5J8Zovp5yS5aqSzAsD2UmwBYGw7OMmq1trdrbXHkyxKMneLeeYmuWjg9WVJfrWqKkmq6n8muTvJHbsoLwBsM8UWAMa2SUnuGfR+zcDYkPO01jYmWZ9k36p6cZI/SvLJZ/uAqppfVUuraum6det2WnAAGC7FFgDGthpirA1znk8mOae19vCzfUBrbUFrra+11jdx4sTtjAkA28/jfgBgbFuT5BWD3k9Oct9W5llTVeOTTEjyUJLZSY6tqs8meWmSp6rq0dba50c+NgAMXyeP2HrcDwAM2y1J9q+qaVX1giTHJ1m8xTyLk5w48PrYJNe2fm9vrU1trU1N8ldJPqXUAjAadfKILQAwPK21jVV1apKrk4xLcn5r7Y6qOjPJ0tba4iRfTPKlqlqV/iO1x/cuMQBsO8UWAMa41tqSJEu2GPv4oNePJvm151jHGSMSDgB2gk6eigwAAACbKLYAAAB0mmILAABApym2AAAAdFoni21VzamqBevXr+91FAAAAHqsk8XWc2wBAADYpJPFFgAAADZRbAEAAOg0xRYAAIBOG9/rAAAADM+KtRty3Hk39jrGNps7a1LmzZ7S6xjAGKbYAgB0wNxZk3odYbusWLshSRRbYEQptgAAHTBv9pROlsMuHmEGusc1tgAAAHRaJ4ttVc2pqgXr16/vdRQAAAB6rJPFtrV2RWtt/oQJE3odBQAAgB7rZLEFAACATRRbAAAAOk2xBQAAoNM87gfYMfcvTy54T69TbLuZxyZ9J/c6BQAAO4FiC2y/mcf2OsH2uX95/1fFFgBgTFBsge3Xd3I3y2EXjzADALBVrrEFAACg0xRbAAAAOk2xBQAAoNMUWwAAADqtk8W2quZU1YL169f3OgoAAAA91sli21q7orU2f8KECb2OAgAAQI91stgCAADAJootAAAAnabYAgAA0GmKLQAAAJ2m2AIAANBpii0AAACdptgCAADQaYotAAAAnabYAgAA0GmKLQAAAJ2m2AIAANBpii0AAACdptgCAADQaZ0stlU1p6oWrF+/vtdRAAAA6LFOFtvW2hWttfkTJkzodRQAAAB6rJPFFgAAADZRbAEAAOg0xRYAAIBOU2wBAADoNMUWAACATlNsAQAA6DTFFgDGuKo6sqruqqpVVXXaENP3qKpLBqbfVFVTB8YPq6pbq2r5wNd37ursADAcii0AjGFVNS7JuUmOSjIjyQlVNWOL2U5J8pPW2vQk5yT5zMD4j5PMaa3NTHJiki/tmtQAsG0UWwAY2w5Osqq1dndr7fEki5LM3WKeuUkuGnh9WZJfrapqrX23tXbfwPgdSfasqj12SWoA2AaKLQCMbZOS3DPo/ZqBsSHnaa1tTLI+yb5bzPP+JN9trT225QdU1fyqWlpVS9etW7fTggPAcCm2ADC21RBjbVvmqarXp//05A8P9QGttQWttb7WWt/EiRO3OygAbC/FFgDGtjVJXjHo/eQk921tnqoan2RCkocG3k9OcnmSD7XWvj/iaQFgOyi2ADC23ZJk/6qaVlUvSHJ8ksVbzLM4/TeHSpJjk1zbWmtV9dIkVyY5vbX27V2WGAC2kWILAGPYwDWzpya5OsmdSS5trd1RVWdW1dEDs30xyb5VtSrJ7yfZ9EigU5NMT/Kxqlo28Ofnd/EmAMBzGt/rAADAyGqtLUmyZIuxjw96/WiSXxtiuT9P8ucjHhAAdpAjtgAAAHSaYgsAAECnORUZAIARtWLthhx33o29jrHN5s6alHmzp/Q6BjAMii0AACNm7qxJvY6wXVas3ZAkii10RCeLbVXNSTJn+vTpvY4CAMCzmDd7SifLYRePMMPzWSeLbWvtiiRX9PX1/WavswAddf/y5IL39DrFtpt5bNJ3cq9TAACMKp0stgA7ZOaxvU6wfe5f3v9VsQUAeBrFFnj+6Tu5m+Wwi0eYAQB2AY/7AQAAoNMUWwAAADpNsQUAAKDTFFsAAAA6TbEFAACg0xRbAAAAOk2xBQAAoNMUWwAAADpNsQUAAKDTFFsAAAA6TbEFAACg0xRbAAAAOk2xBQAAoNMUWwAAADpNsQUAAKDTFFsAAAA6TbEFAACg0xRbAAAAOk2xBQAAoNMUWwAAADpNsQUAAKDTFFsAAAA6bXyvAwAAwGi0Yu2GHHfejb2Osc3mzpqUebOn9DoG7FKKLQAAbGHurEm9jrBdVqzdkCSKLc87ii0AAGxh3uwpnSyHXTzCDDuDa2wBAADoNMUWAACATlNsAQAA6LRRc41tVb0qyZ8kmdBaO7bXeQAAoIvczZnnoxE9YltV51fVA1V1+xbjR1bVXVW1qqpOS5LW2t2ttVNGMg8AAIxlc2dNyoyX7d3rGNtsxdoN+fqye3sdgw4b6SO2Fyb5fJJ/2DRQVeOSnJvksCRrktxSVYtbaytGOAsAAIxp7ubM89WIFtvW2g1VNXWL4YOTrGqt3Z0kVbUoydwkii3AWLT0gmT5ZSO3/l+cmRz16ZFbPwAw6vXiGttJSe4Z9H5NktlVtW+Ss5IcVFWnt9b+z6EWrqr5SeYnyZQp3fttFMDzzvLLkvuX9xdQAKDnLr7phyN66veMl++dT8x5/Yitfyi9KLY1xFhrrT2Y5Leea+HW2oIkC5Kkr6+v7eRsAIyEX5yZnHxlr1MAAEm+vuzerFi7oZPXY29NL4rtmiSvGPR+cpL7epADAAAYJdzNedea8bK9c8mH39rrGDtNL4rtLUn2r6ppSe5NcnySeT3IAQAAjAJzZ03qdYTtsmLthiTpZLEda0a02FbVwiSHJNmvqtYk+URr7YtVdWqSq5OMS3J+a+2OkcwBMGbcvzy54D29TrFtXF8LwHPo8t2cu3ikeaydhpyM/F2RT9jK+JIkS7Z3vVU1J8mc6dOnb+8qALpn5rG9TrB9fnFmd7MDwLPo6pHmGS/bu7PZt6Za6+79l/r6+trSpUt7HQOAMaKqbm2t9fU6R5fZNwOwMw1337zbrggDAAAAI0WxBYAxrqqOrKq7qmpVVZ02xPQ9quqSgek3VdXUQdNOHxi/q6qO2JW5AWC4FFsAGMOqalySc5MclWRGkhOqasYWs52S5CettelJzknymYFlZ6T/6QWvT3Jkkr8ZWB8AjCqKLQCMbQcnWdVau7u19niSRUnmbjHP3CQXDby+LMmvVlUNjC9qrT3WWvvPJKsG1gcAo0oni21VzamqBevXr+91FAAY7SYluWfQ+zUDY0PO01rbmGR9kn2HuWyqan5VLa2qpevWrduJ0QFgeDpZbFtrV7TW5k+YMKHXUQBgtKshxrZ8JMLW5hnOsmmtLWit9bXW+iZOnLgdEQFgx3Sy2AIAw7YmySsGvZ+c5L6tzVNV45NMSPLQMJcFgJ5TbAFgbLslyf5VNa2qXpD+m0Et3mKexUlOHHh9bJJrW/+D7hcnOX7grsnTkuyf5OZdlBsAhm18rwMAACOntbaxqk5NcnWScUnOb63dUVVnJlnaWluc5ItJvlRVq9J/pPb4gWXvqKpLk6xIsjHJ77bWnuzJhgDAs1BsAWCMa60tSbJki7GPD3r9aJJf28qyZyU5a0QDAsAOcioyAAAAndbJYutxPwAAAGzSyWLrcT8AAABs0sliCwAAAJsotgAAAHSaYgsAAECnKbYAAAB0mmILAABAp1VrrdcZtltVrUvyg17n2E77Jflxr0PsRLZndLM9o5vtGT1e2Vqb2OsQXWbfPKrYntHN9oxutmf0GNa+udPFtsuqamlrra/XOXYW2zO62Z7RzfbA6DDW/u3antHN9oxutqd7nIoMAABApym2AAAAdJpi2zsLeh1gJ7M9o5vtGd1sD4wOY+3fru0Z3WzP6GZ7OsY1tgAAAHSaI7YAAAB0mmILAABApym2u0BVvaKqrquqO6vqjqr6yMD4GVV1b1UtG/jz7l5nHa6qWl1VywdyLx0Y26eq/qmqVg58/ble5xyOqjpg0PdgWVVtqKqPdun7U1XnV9UDVXX7oLEhvx/V76+ralVVfa+q3ti75M+0lW35XFX9x0Dey6vqpQPjU6vqkUHfo7/rXfKhbWV7tvpvq6pOH/je3FVVR/Qm9dZtZXsuGbQtq6tq2cD4qP/+8Pxl3zx6jYX9cmLfPJr/77dvHt3fn+3lGttdoKpeluRlrbXbquolSW7N/9ve/YdKVpYBHP8+pIUaWbv9gNxkVVIQ/9CQLMzACmnF3FISRdglF8TMQvzHYP8JQahIoT/EIJJu5o/wR7QEgqSguGCB2+6arKJbhtteVtaFUvyRu/v0x3kvzF5m7s497Z1z3pnvBw4z895zLs/Lc+4888y8cy58E7gKeCszf9ZpgC1ExKvA+Zm5f2Dsp8CBzPxxRPwQ+Fhm3tpVjG1ExAeAfwEXAN+hkvxExJeBt4DfZOY5ZWxoPsoT9feBS2nm+fPMvKCr2BcbMZdLgCcz82BE/ASgzGUt8MeF/fpoxHx+xJBzKyLOBh4APg98GvgTcGZmHppo0EsYNp9FP78D+Hdm3lZDfjS7rM11qLUug7W5z8/91uZ+56ctP7GdgMycz8xt5f6bwC7glG6jWhHrgblyf47mBUJtvgrszsx/dh3IcmTm08CBRcOj8rGe5okvM/NZ4KPlBV4vDJtLZj6emQfLw2eBNRMPrKURuRllPfBgZr6Xmf8AXqEppL2x1HwiImiaggcmGpTUgrW5GlXWZbA295m1eTrZ2E5YeZfkPODPZeimsoTjnhqWBw1I4PGIeC4iri9jn8rMeWheMACf7Cy69q7myD/8WvMDo/NxCvDawH57qOvF3HXAYwOPT4uIv0bEUxFxUVdBtTDs3Ko9NxcB+zLz5YGxWvOjGWJt7rVpqstgbe47a3PFbGwnKCI+DDwC3JyZ/wHuBs4AzgXmgTs6DG+5LszMzwHrgO+VJRBVi4gPApcDD5WhmvOzlBgyVsV3EiJiM3AQuK8MzQOnZuZ5wC3A/RHxka7iW4ZR51a1uSmu4cgXoLXmRzPE2txfM1SXoeLnf2tz781MbbaxnZCIOJ6mcN6XmY8CZOa+zDyUmYeBX9KzZQ1Lycy95fZ14Pc0se9bWDZTbl/vLsJW1gHbMnMf1J2fYlQ+9gCfGdhvDbB3wrEtW0RsBC4Drs1ycYCyLOiNcv85YDdwZndRjmeJc6vK3ABExHHAFcDvFsZqzY9mh7W596atLoO1ubeszfWzsZ2Asrb9V8CuzLxzYHzwuxPfAv62+Ng+ioiTyoU2iIiTgEtoYt8CbCy7bQT+0E2ErR3xjlat+RkwKh9bgA3R+ALNxQTmuwhwXBHxdeBW4PLMfHtg/BPlwiJExOnAZ4G/dxPl+JY4t7YAV0fEhyLiNJr5/GXS8bX0NeDFzNyzMFBrfjQbrM1VmLa6DNbm3rI2T4HMdFvhDfgSzZKFncD2sl0K3As8X8a30FydsfN4x5jP6cCOsr0AbC7jq4EngJfL7aquY13GnE4E3gBOHhirJj80hX8eeJ/mncVNo/JBs6TmLpp36J6nuYJm53M4ylxeofl+y8Lfzy/KvleWc3AHsA34RtfxjzmfkecWsLnk5iVgXdfxjzOfMv5r4IZF+/Y+P26zu1mb+73VXpdLvNbmnj73W5v7nZ+2m//uR5IkSZJUNZciS5IkSZKqZmMrSZIkSaqaja0kSZIkqWo2tpIkSZKkqtnYSpIkSZKqZmMrSZIkSaqaja0kSZIkqWo2ttIUiYi1EfFORGxvefwJEbE9Iv4bER8/1vFJkjRrrM3SZNjYStNnd2ae2+bAzHynHLv3GMckSdIsszZLK8zGVqpMRDxZ3rndHhHvRsS3l9h3bUS8GBFzEbEzIh6OiBPLzzaUsR0Rce/kZiBJ0nSxNkvdO67rACQtT2Z+BSAivgtcDDx6lEPOAjZl5taIuAe4MSIeAzYDF2bm/ohYtaJBS5I0xazNUvdsbKUKRcQGYB1wZWYeOsrur2Xm1nL/t8APgPeAhzNzP0BmHlixYCVJmgHWZqlbNrZSZcrypmuB9Zn5/hiH5JDHMWRckiS1YG2Wuud3bKWKRMRlwI3AFZn57piHnRoRXyz3rwGeAZ4AroqI1eX3utxJkqQWrM1SP9jYSnWZA9YAW8sFKjaNccwuYGNE7ARWAXdn5gvA7cBTEbEDuHPFIpYkabpZm6UecCmyVJHMXN3isMOZecOQ3zVHU4wlSVJL1mapH/zEVpouh4CT/99/Ag8cDxw+ppFJkjSbrM3SBESm31GXJEmSJNXLT2wlSZIkSVWzsZUkSZIkVc3GVpIkSZJUNRtbSZIkSVLVbGwlSZIkSVWzsZUkSZIkVc3GVpIkSZJUtf8BiAhdlzRDFt4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(figsize=(16,8), nrows=1, ncols=2)\n",
    "\n",
    "zmin, zmax, nbins = 0, 200, 10\n",
    "\n",
    "hist_all, bins_all = np.histogram(tic_gaia_table['gal_z'], bins=np.linspace(zmin, zmax, nbins))\n",
    "bins_all = 0.5*(bins_all[1:] + bins_all[:-1])\n",
    "axes[0].plot(bins_all, hist_all, linestyle='steps-mid', label='All TESS-GAIA')\n",
    "\n",
    "hist_fl, bins_fl = np.histogram(df_g['gal_z'], bins=np.linspace(zmin, zmax, nbins))\n",
    "bins_fl = 0.5*(bins_fl[1:] + bins_fl[:-1])\n",
    "axes[0].plot(bins_fl, hist_fl, linestyle='steps-mid', label='Flare TESS-GAIA')\n",
    "\n",
    "axes[0].set_yscale('log')\n",
    "axes[0].set_xlabel('$z$ [{0:latex_inline}]'.format(galcen.z.unit))\n",
    "axes[0].set_ylabel('dn/dz')\n",
    "axes[0].legend()\n",
    "\n",
    "axes[1].plot(bins_all, hist_fl/hist_all, linestyle='steps-mid', label='Flare/All')\n",
    "axes[1].set_xlabel('$z$ [{0:latex_inline}]'.format(galcen.z.unit))\n",
    "axes[1].set_ylabel('dn/dz')\n",
    "axes[1].legend()"
   ]
  },
  {
   "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
}