sample-sim-props.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "import os\n",
    "from scipy.ndimage.filters import gaussian_filter\n",
    "from styles import slides"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "path = os.path.expandvars(\"$phd/inputs/tta_cb13/set_cb13\")\n",
    "sfh_list = [\"SSAG{0:06d}_z0p00\".format(i+1) for i in xrange(140)]\n",
    "\n",
    "mass1gyr, mass1myr = [], []\n",
    "ur_rd, ur_nd, z_nd, t_nd = [], [], [], []\n",
    "for sfh in sfh_list:\n",
    "    t, sfr = loadtxt(\"{}/{}.sfr\".format(path, sfh), unpack=True)\n",
    "    sfh_pr = loadtxt(\"{}/{}nd.tparam\".format(path, sfh), usecols=(1,19))\n",
    "    col_nd = loadtxt(\"{}/{}nd.1ABmag\".format(path, sfh), usecols=(1,2,3,4,5))\n",
    "    col_rd = loadtxt(\"{}/{}.1ABmag\".format(path, sfh), usecols=(1,2,3,4,5))\n",
    "    \n",
    "    mass = (0.5*(sfr[1:]+sfr[:-1])*(t[1:]-t[:-1])).sum()\n",
    "\n",
    "    mass1myr.append((0.5*(sfr[1:]+sfr[:-1])*(t[1:]-t[:-1]))[t[1:]>=t.max()-1e6].sum()/mass)\n",
    "    mass1gyr.append((0.5*(sfr[1:]+sfr[:-1])*(t[1:]-t[:-1]))[t[1:]>=t.max()-1e9].sum()/mass)\n",
    "    z_nd.append(sfh_pr[0])\n",
    "    t_nd.append(sfh_pr[1])\n",
    "    ur_nd.append(col_nd[0]-col_nd[2])\n",
    "    ur_rd.append(col_rd[0]-col_rd[2])\n",
    "\n",
    "mass1myr = array(mass1myr)\n",
    "mass1gyr = array(mass1gyr)\n",
    "t_nd = array(t_nd)\n",
    "z_nd = log10(array(z_nd))\n",
    "ur_nd = array(ur_nd)\n",
    "ur_rd = array(ur_rd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAggAAAIICAYAAAAL/BZjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAOJgAADiYBou8l/AAAIABJREFUeJzt3X2MHfWd5/vPl4Q49G3bQGZ0GdlR6Ka7tRvQAG4yS6S+\nI2WwQ660UkLTThD8Q3BjM6P9Y5CHB4crRatZh4cE5Y5mdieAu5d/QAG3e2ZzpSsTDLma21rYgbZJ\nlCD5dLub3GAlOzsYP/Q248zg7/3jV6e7+jx016lT56nO+yUd1XmoU/Wr6rLre35PX3N3AQAAxF3W\n6gIAAID2Q4AAAADKECAAAIAyBAgAAKAMAQIAAChDgAAAAMoQIAAAgDIECAAAoAwBAgAAKNPWAYKZ\n7TWzeTO7FC0fanWZAADoBm0bIJjZzZJ+IOlDSQ9LOiXpSTMbbWnBAADoAm0bIERc0jFJz0l6Nnrv\nC60rDgAA3cHaOVmTmT2hUHtQ9KGkYXd/rzUlAgCgO3yy1QWoxsx2KgQHb0t6VNJuSbdIOlNh3f9F\n0r+R9GtJ/9zEYgIA0E4ul/R7kv6bu//PejbUtgGCpDGFJoZH3P0nZjarUIPwnKRvlKz7byS91uTy\nAQDQrm6T9Ho9G2jnAKHoqmj5mWi5o8I6v5akY8eO6XOf+1xTCtVu9u/fr6effrrVxWgpzkFzz8Hu\n3dI770g33SQdPtyUXW6Ia4BzIHX3OfjlL3+pnTt3StF9sR7tHCAclrRX0pSZTUnaqVCjUOm/on+W\npM997nMaGBhoXgnbyObNm7v22Is4B809B5/+9OqyXU471wDnQOIcROpubm/bUQzu/ppCM8OspDsl\nfSDpYXf/VksLBgBAF2jnGgS5+7Sk6VaXAwCAbtO2NQiozego80dxDjgH3X78EudA4hxkpa3nQUjK\nzAYkzc3NzdHuBDTJ3XdLP/2pdOON0osvtro0ACRpfn5eg4ODkjTo7vP1bKutmxgAtC+CAiDfaGIA\nAABlCBAAAEAZAgQAAFCGAAEAAJQhQAAAAGUIEAAAQBmGOQJI5dAh6f33pe3bpfHxVpcGQNYIEACk\nMjEhvfmmdOutBAhAHtHEAAAAyhAgAACAMgQIAACgDAECAAAoQ4AAAADKECAAAIAyDHMEkMoLL0gf\nfSRdcUWrSwKgEQgQAKTS39/qEgBoJJoYAABAGQIEAABQhgABAACUIUBogUKhoMnJSRUKhVYXpTaF\ngjQ5GZYAgFyjk2KTFQoFDQ8Pa2lpSb29vZqdndXQ0FCri7WxQkEaHpaWlqTeXml2VuqEcgMAUqEG\noclmZma0tLQkSVpaWtLMzEyLS5TQzEwIDqSw7JRyo2EWFqRf/CIsAeQPAUKTjYyMqLe3V5LU29ur\nkZGRFpcooZGRUHMghWWnlBsNc8890g03hCWA/KGJocmGhoY0OzurmZkZjYyMdEbzghSaE2ZnQ83B\nyAjNCwCQcwQILTA0NNQ5gUHc0BCBAQB0CZoYAABAGQIEAABQhgABAACUIUAAAABl6KQIIJU9e6Tb\nb5e2b291SQA0AgECgFTGx1tdAgCNRBMDAAAoQ4AAAADKECAAAIAyBAgAAKAMAQIAAChDgAAAAMoQ\nIABI5e67peuvD0sA+cM8CABSWVyU3n1X2rKl1SUB0AjUIAAAgDIECAAAoAwBAgAAKEOAAAAAyrRt\ngGBmL5vZpQqPj1tdNgAA8q6dRzH8UNIHsddXS9otabY1xQEQ19cnnT8flgDyp20DBHefljRdfG1m\nz0hyhSABQIu9+GKrSwCgkdo2QIgzs52S7pd0v7v/stXlAQAg79q2D0KJw5Lm3X2i1QUBAKAbtH0N\ngpntlbRF0n0brbt//35t3rx55fXo6KhGR0cbWDoAAFpnenpa09MrrfG6cOFCZts2d89sY41gZq9K\n+iNJV7n7+SrrDEiam5ub08DAQFPLBwBAu5ifn9fg4KAkDbr7fD3b6oQmhtskHa8WHAAAgOy1dYBg\nZn0KIxdOtbosAAB0k7bug+Dui5I+0epyACh36JD0/vvS9u3S+HirSwMga20dIABoXxMT0ptvSrfe\nSoAA5FFbNzEAAIDWIEAAAABlCBAAAEAZAgQAAFCGAAEAAJQhQAAAAGUY5ggglRdekD76SLriilaX\nBEAjECAASKW/v9UlANBINDEAAIAyBAgAAKAMAQIAAChDgAAAAMoQIAAAgDKMYgCQysLC6jBHRjQA\n+UMNAoBU7rlHuuGGsASQPwQIAACgDAECAAAoQ4AAAADKECAAAIAyBAgAAKAMAQIAACjDPAgAUtmz\nR7r9dmn79laXBEAjECAASGV8vNUlANBIdQcIZrZF0i2S+iVdGb19VtKCpLfd/Xy9+wAAAM2VKkAw\ns2sl7ZM0phAYrLfugqTDkp519/fS7A8AADRXTZ0UzewmM3tF0ilJDysEB7bBo1/SI5JOmdlRM7sx\nu+IDAIBGSFSDEDUjPKdQYyCFG78UmhKOSXoren4mev9qheaGL0jaqdWmh12SdpnZYUl7aX4AAKA9\nJW1i+DBamkIg8KykZ9x9McmXzaxP0gOS7lcIFr6uEGzQSRIAgDaUtInBJJ2QtNvdr3b3R5MGB5Lk\n7ovu/oi7X60QHJzQai0EAABoM0kDhF3ufou7H6l3h+4+5e7Dkr5c77YAtM7dd0vXXx+WAPInURW/\nu7+WZD0zu0mhz8EXFIY5vuTu79SzTQDtaXFRevddacuWVpcEQCNk1gfAzB5XGNkgheYDl/SwmT3p\n7t/Kaj8AAKDxMsnFYGb3KwxlNEmLkqaipUl6xMzuyGI/AACgObJK1rQvWr7q7gPu/nV3H1BoZjBJ\nBzLaDwAAaIJEAYKZ/dEGq+xQaFI4XPL+VLQcrrFcAACghZLWIBwzs79fZxbE49FyrOT9MYXA4bgA\nAEDHSBogPCppQNLxaLrkz5V8/qxCU8IuM3vLzP7azOa0mqfhmWyKC6Bd9PVJn/98WALIn6TDHJ+S\n9JSZPSnpIUkL0XTJj7j7L939WTO7LvpsWKHJoTgR0lPufqgBZQfQQi++2OoSAGikmjopuvsjkq6S\ndEhhRsQFM/tPZrY5+uwWSU9JOhIth92dDooAAHSYmkcxuPs5d9+n0OQwrZBj4ayZfUfSfDQN89ej\n5YmMywsAAJog9TBHd19w990KgcLrCv0UFs1sf1aFAwAArVFTgGBmfxZ1QpyLmhY+FwUKuxSaF96T\n9F0z+8DM9jSiwAAAoPESBwjRqIQnFTohXqcwOdJCceijux+PJWF6T9JzZlZgFkUAADpP0omSnlAI\nCkxhSOMjks5Fr9dMjuTux6JA4euSPiFpysz+PstCAwCAxkqarKk44dGz7v7HkmRmr0l6W9J1Znat\nu78X/4K7TykEB3slPZFdkQG0g0OHpPffl7Zvl8bHW10aAFlL2sRQnPBotviGux+v8HkZd3/W3a9O\nUTYAbWxiQvr3/z4sAeRP0gChOFxxb/ENM/uz2OdvZ1aiGDPbamaHzexM9KAmAgCAJkjaxPCIpB9L\nGjazDySdUag1KDY7nG9Q+Y5LularnSMfMrN5ZmYEAKCxkk61fMzMHlXoS3BV9JCkw8U+CVkzs9sk\n9Ul62N2/F723pYHBCAAAiCQe5hjlY7hK0m6FIY7D7v6NRhVMIZ+DJH3ZzC6Z2SWFKZ4BAECDJW1i\nkBSmWVbIs9AMn4mWfQp9H3ZLGjOzx6vld9i/f782b9688np0dFSjo6MNLygAAK0wPT2t6enpldcX\nLlzIbNvm7huvZPaSpO+4+08z2anZTZIOrFcDEXVIfEjSXnefMLOtkj6UNOvuXyhZd0DS3NzcnAYG\nBrIoIoANfPGL0ptvSrfeKr3xRqtLA0CS5ufnNTg4KEmD7j5fz7aS1iAUf72/KukZd/+bNDszs1FJ\nB7TafLCeUwoTMRX7O3wm9j6AFnvhBemjj6Qrrmh1SQA0QtIA4YikOyXtkrTLzM5KeknSlKS3q3Uc\nNLMtCjkadivMrHilwk1f0XercvfnzOxJSU+a2R9I2qkwauKlhGUG0ED9VWc/AZAHSUcx7DazHZKe\nk3Szwq/6fdFDUcCwEPvKlZKujpZFxcBgQdLuhKmgd0h6RiE4OSXpobS1FwAAILnEnRSjmROHzWyn\nwrwIt8U+vkrlzQZW8vqYpCfd/bUa9vmepNuTrg8AALJR0ygGKcyJIOmYmfUpVPvvipZXlqx6ViEo\neFXSMXdfrLOsAACgSWoOEIqiG/5z0QMAAORI4omSAABA90hdgwCguy0srA5zZEQDkD/UIABI5Z57\npBtuCEsA+UOAAAAAyhAgAACAMgQIAACgTF2dFKOplHdK2rCLkrt/r559AQCA5kkdIJjZ45IeruEr\nBAgAAHSIVAGCmd2vMN1yUhvnlAYAAG0jbQ3CvmjpCjkXFrQ2WROAnNuzR7r9dmn79laXBEAjpA0Q\ndmi1VmBXLQmYAOTD+HirSwCgkdKOYjirUHNwluAAAID8SRsgHI6WpRkcAQBADqQNEB5W1OfAzPZn\nVxwAANAO0vZBeFYhQOiX9JSZ3aX1Oym6u9+Vcl8AAKDJ0gYIuxU6KRZHMeyIHpVYtB4BAgAAHaKe\nmRStynMAANDh6mliANDF7r5b+ulPpRtvlF58sdWlAZC1VAGCuz+QdUEAdJbFRendd6UtW1pdEgCN\nQDZHAABQpu4AwczGzewlM5szs4+j5UtmxjxrAAB0qHqyOV6rMGFSfPSCKQx97Jc0ZmZ7Je1291/W\nUUYAANBkqWoQzGyrpONaDQ5MqyMZ4s9vkTRrZrRSAgDQQdI2MezV6jTLJmlKYW6EXdHyiFaDhKsk\n3V9HGQEAQJOlbWLYFXu+292PlHx+xMzGJL0cvf6ypKdT7gtAG+rrk86fD0sA+ZM2QLhFYXbEhQrB\ngSTJ3afMrDgd8y0p9wOgTTH3AZBv9Y5iOFvn5wAAoA2lDRAWFOVgMLMbK61gZn0KnRhd6ydyAgAA\nbSZtgHAsWpqk181sTzTsUWZ2bTQHwtux9d8WAADoGGn7IDwj6WGF2oGrFOVmMFuTsyn+4smU+wEA\nAC2QqgbB3RcVRiaUzn0QnwOh6Ovu/l7aAgIAgOZL3UnR3Y8p1B48p/LOiOcUahWuqjbKAQAAtK/U\nUy1Lkrufk7QvesjM+qLaBQA5d+iQ9P770vbt0jiZV4DcqStAKEVwAHSPiQnpzTelW28lQADyKFGA\nYGavRE/d3b8Se52Uu/tXavwOAABokaQ1CLsURixUe70eq2FdAADQBmppYii90ZeOVgAAADmRNEB4\nZIPXAAAgRxIFCO7+3fVeAwCAfKk3WRMAAMihVMMczWxe0iVJX15vlkQze0jSXkmnGMUA5MsLL0gf\nfSRdcUWrSwKgERIHCGa2JfayX6HDYp+ZnanylaslfUPSddFzADnS39/qEgBopKTzIDwh6aHYW8XR\nDMcqrA4AADpcmmGO8eGOGw119OhBumcAADpIrZ0Ua537wCSdUJSroVZmdtjMLpU8mNQVAIAGSzrM\n8VEzezx6aZLOKNQM7JQ0u873ztVZvpslnZL0qlZrLhbq3CYAANhA4iaG+M3ezF5TdLPOIAhYT7+k\nve5+qIH7AAAAJVLNg+Duu9x93SGO9TKzvujpo7HmhR80an8AarOwIP3iF2EJIH/qmijJzG4ys782\ns5dK3n/JzL5jZtfWsfl+SR8qdHDcq9CUcb+Z/Vkd2wSQkXvukW64ISwB5E+qiZKklUmQnohelv6G\nGJY0JukRMxtz97+pdfvu/pqkz8T2N6sQJOyS9L1K39m/f782b9688np0dFSjo6O17hoAgI4wPT2t\n6enpldcXLlzIbNtpZ1K8WdKTSVaVNGVm19XaHGFm9ysEA4+7+wmtBgsfVvvO008/rYGBgVp2AwBA\nxyr9ITw/P68f/ehHmWw7bRPDgWjpCpMlPVDy+QOS3ou9fkK1O6NQC3E46nvwcrS/Z1JsCwAA1CBt\ngLAjWp6NOiu+Fv/Q3Y+5+3UKTQ+m0ORQE3c/otD3wCXdL+kfJY25+09SlhkAACSUNkAo5mLYaIbE\nYt+EVLO2u/shdx9090+4+1CavgydrlAoaHJyUoVCodVFab5CQZqcDEsAQFOl7aR4XKEW4ZZqK5jZ\n1uhzJjdKqVAoaHh4WEtLS+rt7dXs7KyGhoZaXazmKBSk4WFpaUnq7ZVmZ6VuOXYAaANpaxCKNQdX\nmtnfm9mXitkezWyLmf1RtM6V0XrH6yxnV5qZmdHS0pIkaWlpSTMzMy0uURPNzITgQArLbjr2DrFn\nj/Ttb4clgPxJW4PwjEL/ACn0LzgmSWZVUzXQsTCFkZER9fb2rtQgjIyMtLpIzTMyEmoOijUI3XTs\nHWKcrChArqUKENz9hJk9IKnazIbFrI+S9Ki7v55mP91uaGhIs7OzmpmZ0cjISPc0L0ihOWF2NtQc\njIzQvAAATZZ6oiR3f9bMjinMh3CbVpsTpBAcHJO0z90X6ytidxsaGuquwCBuaIjAAABaJHWAIEnu\nviBpt7TSKfFqSWcanMAJAAA0WF0BQlwUFBAYAACQA3UHCFFCpivXWaVf0pPuPljvvgAAQHPUk6zp\nToXpjwEAQM6kmgchStZ0WKEzYpIHgJy5+27p+uvDEkD+pK1BiCdrspLnZ7Xa5OAK/RI2mpIZQIdZ\nXJTefVfasqXVJQHQCPUma5KkqxSmVDZJp9z9aoXRDM/F3vtyXaUEAABNVW+ypuPufs7djyvUHPSb\n2RZ3P+vu+xRyMOwws/0ZlRcAADRB2gDhbIX3is0IO2PvFdM935VyP0igYzI+kp0RADpG2j4ICwrN\nDDvM7Evu/hOFhEw7FXI0TEcTJxWDhR2VN4N6dUzGR7IzAkBHSVuDUBzeaJKOmdmopJei93aZ2ceS\nzsTWJ91zg3RMxkeyMwJAR0kVILj7U5IWS947EXuvOLzRowfZHBukmPFRUntnfCxmZ5TIzpgTfX3S\n5z8flgDyp56ZFHdIekpRLobIboX5EYr/ZZikV939e3XsB+vomIyPZGfMnRdfbHUJADSSuXv2Gw0T\nKfUrjHJoeDZHMxuQNDc3N6eBgYFG7w4AgLY0Pz+vwcFBSRp09/l6tpWqBsHMxrU6GdKUu78X/zxq\nbjhRT8EAAEDrpG1ieErS1uj5gqT3MikNAABoC2lHMRRnSZRCUwIAAMiRtKMYHpE0pRAkHDCzzZmW\nCgAAtFTaPggvKQxfXFCoQThrZsdUeYZFSXJ3ZzZFAAA6RNo+CLsVAgRpNYvjzirrFudDIEAAcuTQ\nIen996Xt26Xx8VaXBkDW6pkHwTZ4DSDHJiakN9+Ubr2VAAHIo7QBwrOZlgIAALSVVAGCuz+QdUEA\nAED7SDWKwcyuNbNrsy0KAABoF4kCBDO7ZGYfm9l3orcWJM2b2U2NKxoAAGiVWpoYTCGV87HY61vM\n7OokX3b312stHAAAaI2kAcJZhamVd0h6VatDHJOmcfYa9gUAAFos6U37NUl3ajUwKA5pTJIKkuGP\nQA698IL00UfSFVe0uiQAGiFpgFDM3jgcLUsDBQBdpp8sLECuJQoQ3P2cpF3F12Z2SSFI2OHu7zSo\nbAAAoEXS9gt4TSFAqJZ7AQAAdLC0EyXt2ngtAADQqVJNlAQAAPKNAAEAAJRhbgIAqSwsrA5zZEQD\nkD/UIABI5Z57pBtuCEsA+UOAAAAAyhAgAACAMnUFCGa2xcxGzWxL7L3HzewtM3vFzPbXX0QAANBs\nqTspmtmopMPRy2FJ75jZywo5G4p2mtkud/9KHWUEAABNlqoGwcz6JE0ploshem+s+DL22GVmd9RZ\nTgAA0ERpmxgeiZYu6YTClMtjsc+Px9aRpG+l3I8kycyeMbNLZvZBPdsBAADJpG1i2BktF9z9Fkky\ns/j0y+Pu/o6Z7ZPUHz1SMbN+SfcrWWppAE2yZ490++3S9u2tLgmARkgbIPQr3LCPx97bGb13Npbh\n8Xi07pWpSyj9QNIzkvbVsQ0AGRsfb3UJADRS2iaGYhbHfkkys3jHxJdjz68sWb8mZjYm6Tatba4A\nAAANljZAWFDogLjDzMYlHYh9dliSzGyrVmsVFlLu51lJe939fMrvAwCAFNI2MbwsaUf0/JnY+2fd\n/fVoRMO8VvsNHKt1B2b2cPT916J+CMX3t7r7uUrf2b9/vzZv3rzyenR0VKOjo7XuGgCAjjA9Pa3p\n6emV1xcuXMhs2+aeru+fmZ2S1Ffy9pPufiAKEE5pNUC4qtZaADP7sULzQnEopUfPT7n7YMm6A5Lm\n5ubmNDAwUOORAACQD/Pz8xocHJSkQXefr2dbqSdKcvfrzOxJhZv4lZKm3P1AyWrnJN2WsongYa0d\n/TAl6UNJe9OUFwAAJFdXumd3r9h50N0XzWzY3U/Use13JBVHQ8jMorf9J2m3CQAAkmlYsqZ6goMq\n27vM3X8ny212g0KhoMnJSRUKhVYXpT6FgjQ5GZZoC3ffLV1/fVgCyJ+6ahCiJE07JR0rNiOY2ePR\ne2ck/djdn667lB2qUChoZmZGIyMjGhoaynrj0syMNDIiVdl2oVDQ8PCwlpaW1Nvbq9nZ2ezL0QyF\ngjQ8LC0tSb290uxs1WNG8ywuSu++K23ZsvG6ADoPyZoapKE354Q3zJmZGS0tLUmSlpaWNDMz05kB\nwsxMOFYpLGdmCBAAoMFI1tQglW7OGW68/IZZwcjIiHp7eyVJvb29GhkZya4MzTQyEgIhKSw79TgA\noIN0RLKmTtTQm3PCG+bQ0JBmZ2c1MTHRuc0LUqgtmJ2VJiZoXgCAJmn7ZE2dqnhzbkgfhOINc4M+\nCMVydGxgEDc0RGAAAE3UCcmaOlZDb87cMAEADZQ2QDgraasanKwJQPvq65POnw9LAPmTNkBYUMjF\nUEzW9EDssyyTNQFoUy++2OoSAGiktJ0U47UEz0i6OXoeT9Z0RnUkawIAAK2TKkBw96ckLWrtcEYp\npGcustjzx1OVDgAAtEQ7J2sCAAAt0rbJmgAAQOt0TLImAADQPHXVIEiSmV2r9ec56Jf0pLsP1rsv\nAADQHPUka7pTa0czAOgihw5J778vbd8ujY+3ujQAspYqQDCzm7WayTEJ33gVAJ1kYkJ6803p1lsJ\nEIA8StsHoThaIX7jLz4/W/LeWUmvpdwPAABogbQBwo7Y86sk3aIw78Epd79a0tWSnou99+W6SgkA\nAJoqbYCwkqzJ3c+5+3GFmoJ+M9vi7mfdfZ+iKZnNbH9G5QUAAE2QNkColHzp7Wi5M/begkItwl0p\n9wMAAFogbYBQvPHvMLMvRe8dj97bK61J1iStbZJAxgqFgiYnJ1UoFFpdlGwVCtLkZFgCAJoq7TDH\nlxVu+ibpmJntlvSSpIcl7TKzj0vWJ5tjgxQKBQ0PD2tpaUm9vb2anZ3V0NBQq4tVv0JBGh6Wlpak\n3l5pdlbKw3EBQIeoN1lT/L0TsfeKCZw8ejxTRxmxjpmZGS0tLUmSlpaWNDMz0+ISZWRmJgQHUljm\n5bhy5IUXpJ//PCwB5E89Uy3vUBipcC723m5J78Vem6Rj7v69OvaDdYyMjKi3t1eS1Nvbq5GRkRaX\nKCMjI6HmQArLvBxXjvT3S9dfH5YA8qeebI7nJO2LHsX3jku6LppIqV9hlMNilU0gA0NDQ5qdndXM\nzIxGRkby0bwgheaE2dlQczAyQvMCADRZ3bkYKomaG0jW1CRDQ0P5CQzihoYIDACgRerJxbBFYUbF\nHQoTI63H3f0P0u4LAAA0V9pcDH2S5pOuLnIxAADQUdLWIDypcOMHAAA5lDZA2KnVWoHXVFtmRwA5\nsLAgffSRdMUVjGQA8ihtgHBltHQSMQHd6Z57VtM9v/FGq0sDIGtp50EojlColJMBAAB0uLQBwiPR\n8kozuyOrwgAAgPaQqonB3Y+Z2SMKnRWnzGxB6+dbcHf/Spp9AQCA5qtnmOMBhY6KJuk6hZkTK64u\nhjkCANBR0jYxPKlYR8VoaVUeaIBGpXhuaurodkjn3A5lAIA2lMUwx0VJU9kUB0k0KsVzU1NHt0M6\n53YoAwC0qXpyMRSbDna4+/mMyoMEKqV4zuJG3qjtVtlZeTrnZt+c26EMHWzPHun226Xt21tdEgCN\nkDZAeE3SnZLOEhw0XzHFc/GXflYpnhu13So7C7/ai7/eW5HOuR3K0MHGx1tdAgCNZO619x+MOinO\nStoqaae7/yTrgtVYngFJc3NzcxoYGGhlUZqmUCg0JMVzo7ZbZWetT+fcDmUAgIzMz89rcHBQkgbd\nPWnOpIrSBgiPSxrWal+EY9p4mOOfpCphsvJ0XYAAAECpLAOEtE0MjygEBsVhjjsTfKdhAQIAAMhW\nvZ0UKz2vhHkQAADoIGkDhCPipg8AQG6lnWp5d9YFAQAA7SPtTIrIiabOnNgszI7YFHffLV1/fVgC\nyJ96+iCgwzV15sRmYXbEpllclN59V9qypdUlAdAIbV2DYGYPm9m8mV0yszNm9kSry5QnlWZO7HiV\nZkcEANSsbQMEM7tT0hOSzkh6OFo+ZGajLS1YjhRnTpTU+JkTm6U4O6LE7IgAUId2bmJYkLRX0ssK\nwyj/QFKf1p+QCTUYGhrS7OzshjMnNnV2xbU7rn2Ww6Gh0KzA7IgAUJe2DRDc/YSkE9G0zqcUhlU+\n6+7vtLZk+TI0NFR2048HBJJq76eQxfTF9fQlGBoiMACAOqUKEMxsXNKV0cspd38vsxKVcPdFM9up\nMHvjPjN71d2nG7W/blfacfHAgQO1ZXjMqpMgmRYBoKXS1iA8pZCoSQpV/u9lUpoYM9sqqd/dT7j7\n62a2qFCT8A1JFQOE/fv3a/PmzSuvR0dHNTpKl4ValHZclFRbhsesbuxkWmx7fX3S+fNhCaA1pqen\nNT29eku8cOFCZttOm6zpSUkPKVT7P+Lu38usRKv7eEKhc+Jedz8U1SL8WNJhd/9Gyboka8pIpaGP\nkpL3QchymCGZFgGgJi3P5ihJZvaypDFJH0q61t2zC1u0klL6VPRySiEhVMX00gQI2aq7UyI3dgBo\niZZnczSzlxRqDxYk9Us6a2bHJJ2t8hV397tq2Ues78Ezku5UCBYeKg0OUF3aG32ljos12aiTIAEE\nALS9tH3Z4tAWAAAgAElEQVQQdms1WdNGKZ8tWqemAEGS3P11SYNpCtjt2naWRGY6BICOUM9ESRZ7\nlL4u/azrtDrHwYazJCbIV9CQY8h6pkPyLgBAQ6StQXg201LkTDv8ei/Oklhx9EGCX/ENO4YsRydQ\nGwEADZM23fMDWRckTyr9em92gLDuLIkJhiI27BiynOmQuRIAoGHSdlLc4u7nsy5MXqz7672J4p0N\n13RYTPArvqHHkNVMh/Hj6OmRfvObUKtAkNAUhw5J778vbd8ujY+3ujQAspZ2HoSPJc0qjDA43Opg\noR2HObYsf0GVspQ1F0gb/opvp2OoqlCQpqakgwel5WWaGproi1+U3nxTuvVW6Y03Wl0aAFIbDHNU\n6Hw4rNAX4VkzOy7pB2qDYKFd1D1UMEMVmwvuu2/Dm2g7HUNVQ0PSNdeE4ECiqQEAMpLVKIYdCsHC\nh2b2lpntMbMtWRQQ9ctlWuc4UjwDQObSBgjXSdon6Vj0mmChjRU7LE5MTLTPfAhZKnZ8nJigeQEA\nMpJ2FMOipOeih8zsZoUkSjsVAoQimiHaREc0F9SDFM8AkKl6mhhWRBkXH3X3WyRdpZDtUVqdZbG0\nZuE/ZbFfAADQGGk7KZYxs1FJuyR9XdKVKp+KOT6r4j4zO+XuT2e1fwAAkJ3UAYKZXauQzXGX1uZh\nKJ1eeUFhOORxSV9WSBNtCrkZCBCADvXCC9JHH0lXXNHqkgBohLQTJc0pZHFceatklQWFFM3PRP0V\nil6L0jiPlXwfQIfp518wkGtpaxCu02oTQjE4OKvQx+Aldz+xznffUggQrky5bwAA0GD19EEwhaDg\nZYWRCa8l+ZK7f1fSd+vYLwAAaLC0AUKx+SBRUAAAADpL2nkQvp51QQAAQPtIFCCY2St17sfd/St1\nbgNtpiOSOdWjUMgmLTUAdKCkNQi7tNopsVZWx3fRpipmiMzTTbRQkIaHV1NiM4VzmYWF1WGOjGgA\n8qfWmRQtxQM5VClDZK7MzITgQFrNEIk17rlHuuGGsASQP2n6ILjCPAcLIgDoWsUMkcUahNxmiCzW\nIOTt+ABgA7UGCMWmgn5JVytkc3xL0jF3fyfLgqG9FTNE5rYPQjFDJH0QAHSppAHCcwqZGeOZGq9S\nmPBoTJLMTArTKb8t6ceSXiNrY76RIRIA8itRgODu+4rPo9TOtyh0XNyhtVMmF4OIvdG6ZxXVMrj7\n9zIqMwAAaLCa+yBE0yifUKhVkJlt1WrAMCzpttjqV0narVDLQIAAAECHqDvds7ufk/Ra9FCUjGmf\npIe1muoZAAB0kLoCBDPborXNDfG0zwQHQI7t2SPdfru0fXurSwKgEWoKEMzsWoUgYDhalk6PEg8I\nFhQ6Lb6q0HERQI6Mj7e6BAAaqZaplneWvl3y+rjWDns8V3/xAABAK9Qy1XLch2IOBAAAcquWJobS\nfAo7o8eBaA6Edb/r7r9TS8HQWpkmYiLpEQB0nLSdFK+Mlkk7IZKsqYNkmoiJpEcA0JFqSdZEIqYu\nkWkiJpIeAUBHSlqD8FRDS4G2kiQRU+ImiHqTHtE8AQAtkXSq5UcbXRC0j40SMdXUBFFP0iOaJ9ra\n3XdLP/2pdOON0osvtro0ALJW90yKyKf1EjFVaoJYtxYhbdKjSs0TBAhtY3FRevddacuWVpcEQCPU\n0gcBkLTaBCGpahNERjsKNQdhR7U3TwAAUqMGATWr1ASR6bDI1R2lb55Iiz4PACCJAAEpxZsgMh0W\nWb6j5t2o6fMAACtoYkDdMh0W2UoMyQSAFQQIqFvT+iQ0Gn0eAGAFAQJSKxQKmpyclCTNzs5qYmJC\ns7OzkqTJyUkVCoVKX5ImJ8Oy/gJkty1ptc/DxATNCwn09Umf/3xYAsifVH0QzGxe0iVJX3b399ZZ\n7yFJeyWdcvevpCoh2lKlfgf33Xff+v0RYm38l3p6NP3YY/r9sbF0/RXS9hfYqBNiM/s8dDjmPgDy\nLXENgpltKT4k9Uu6TlJf/P2Sx7WSvhGt94VGFB6tU63fwbr9EWJt/JctL+v/fuwxDQ8PV65p2LgA\ntfcXKAYVe/aEZVY1DwCQQ4kCBDN7QiHFc/FRTL50rOT9+OOUpJszLi/aRLV+B+v2R4i18V+QNKM6\nOjWm6S9AJ0QASKyWJgZTCAyKy+J76/Ho8XbtRUM7qzYd87rTNEdt/L+ZmtL/fvCg5paX03dqTDNH\nQr15IQCgi5j7xpmYoxqEh2NvJQ0QJOm4pLH1+iqss98xSU9K6pO0IOkZd/9uhfUGJM3Nzc1pYGCg\n1t10jIZMRpS1hBMNtexYmAgJQI7Nz89rcHBQkgbdfb6ebSUKECTJzLYWn0o6oxAk7JI0W+077n4u\ndcHM+hSaKT6U9LKknQp9H/a5+6GSdXMfIDRyMqLMbtZZTzTEzRwAapJlgJC4iSF+szez1xQChIV6\ngoAN7Iz28bC7T5jZzQrByG5Jh9b9Zg7VnCApoUwDjyyTKzGrIQC0VKphju6+K+uCVPCqQg1Fsf9C\nsTnjwybsu+0UO/8Vb+RZTUaUaeCRZRs/mRzb3qFD0vvvS9u3S+PjrS4NgKzVlYshGvK4U9Ixdz8f\nvfd49N4ZST9296fTbDvqs/BetM0dkg5H23y0njJ3qnU7/9Uh08Ajy+RKdChsexMT0ptvSrfeSoAA\n5FHqAMHM7lToGyBJw5LeMbOXJd0ZW22nme2qZ5IkM9sr6QeS5iX1ufuFauvu379fmzdvXnk9Ojqq\n0dHRtLtuO/EESVluM9PAI6uJhlqRyREAOsz09LSmp6dXXl+4UPUWWbPEnRTXfGm1A6EU+gkMSzoX\nvVccCqnY52Pu/jcp9vOMpPslvezud62zXu47KQLt5otfXK1BeOONVpcGgJRtJ8W0uRgeiZYu6YSk\ns5LGYp8fj60jSd+qdQdmtlMhOHBJZ83sB9Hj/nRFBgAASaVtYtgZLRfc/RZJMrN4x8Vxd3/HzPYp\nDE3sT7mPYvVGPCjok/Rciu1BdQxpZMghAHSVtAFCv8LN+3jsveIN/ay7vxO9dzxa98pad+Duj6pL\nOyQ2Si1DGtcEElJmQw47YrKnJAiY6sP5A9pe2gDhrKStimoGog6LRS/Hnl8ZWx8tlnRIY2kgMXfg\ngK7JYMhhIyd7airmaKgP5w/oCGn7ICwodETcYWbjkg7EPjssrcy8WKxVWKinkMjGuomUYsoCifAF\nqbhMOeRw3UyPnYSkT5KkF16Qfv7zsKwJ5w/oCGlrEF6WtCN6/kzs/bPu/no0ymFea7M+osWSDmks\nnRvh98fGpLGxuquEGzXZU9MxR4MkqT9NzyKJ8wd0iFTDHCXJzE4pdBiMe9LdD8SGQRY3flVxIqVG\nYJhj9hrVV4A+CJDE+QMapCW5GEq5+3Vm9qSk2xT6Gky5+4GS1c5Juq2RwQGqq+dmPDQ0FDonFqt/\nM/pPvBGTPbVEVhNCdSvOH9D26ppq2d0fqfL+opkNu/uJeraP9OruEEhHMgDoamk7KW6I4KC1qnUI\nLBQKmpycVKFQ2GgD6TqSFQrS5GRYAgA6Vl01CJJkZtdq/XkO+hX6JgzWuy8kV6lDYE21Cmk6klHr\nAAC5kVWyJrSZSiMWJicnk6d2TpMsiRTNXWVhQfroI+mKK+oY0QCgbaUKEMzsZkXzHSSUbqgE6lLa\nIbDmYYa1diRj+FpXuecekjUBeZa2BqE4WiGeubH4/KxWmxxcYSTD22kLiOxkntq5fAekaAaAnEgb\nIOyIPb9K0nUKQcApdx80syslPamQZOmUu3+5vmIiKw0fZsjwNQDIhbSjGFaSNbn7OXc/rlBz0G9m\nW9z9rLvvU5hieYeZ7c+ovAAAoAnSBgiVki8VmxF2xt4r5my4K+V+0CUSD79sFoZrAuhyaZsYFhSa\nGXaY2Zfc/ScKqZ13StoraTqWrEla2yQBrNF2WR4ZrgkAqWsQisMbTdIxMxuV9FL03i4z+1jSmdj6\nZHNEVW2X5ZFsgwCQLkBw96ckLZa8dyL2nkUPjx7PCKiiNA31tm3bWtvcUByuGQqU/XDNnDRf7Nkj\nffvbYQkgf+qZSXGHpKck7Y69t1thfoRilkeT9Kq7f6+O/aDDbZQ0Kj78ctu2bRobG0vX3JBVhsBG\nDtfMUfPF+HirSwCgkerJ5nhO0r7oUXzvuKTroomU+hVGOSxW2QS6QNL+BcXhlzXN9rh2R9neeBs1\nXJPZJgF0iIYka3L3E+5+hOAAtfYvKG1u2HC2x9UddUa/gUY3XwBARhLVIJjZK3Xux939K3VuAx2o\n1umdU8/22CnTPDPbJIAOYe4bp0kws0tKn0/BFAKET6T8/sY7MBuQNDc3N6eBgYFG7QYpxfsgSGrc\nVM+lfRCy6pMAAB1ifn5eg4ODkjTo7vP1bKuWPgi28SpAuWL/gobPdxDvN5CjzoAA0ApJA4RHGloK\ndIVK/REaNiESnQEBoC6JAgR3/26jC4L8i/dH2LRpk7Zt29bIna32Sdi0SWrkvrrU3XdLP/2pdOON\n0osvtro0ALLWkFEMQCVDQ0OamprSpk2bdPHiRY2NjTVuMqShIWlqKgQHFy9KY2MdPzFRu1lclN59\nNywB5A8BAhLJKpnS6dOndfHiRUk1TKucdubB06dDcBB21r5DHwGgDdUzkyK6RJadC2sd9lhXZ8NO\nGfoIAG2IGgRsKMtkSsV5DiYmJpIFGvVMgFScc2BiglEMAFAjahCwoZp/9W+gOOwx4c7rqwVo1JTJ\nAJBzBAjYUOrZDbPZOTMPAkAL0MSA9jc0JN13nzQ0pFdeeUXf/OY39corNcz+XUsnx1alYu7AFNB9\nfdLnPx+WAPKHGgRsqOEzICb0yiuv6CtfCSk9nn/+eT344IN64IEH1i9LaSfHqakwuqFSbUSrZl/s\n0FkfmfsAyDdqELChLDsp1uOHP/zhmtff//73NTw8vP7Qy9JOjl/9qrRnT7ghl36vVRkhOyUTJYCu\nQoAASevPc5A6BXP1naWqTr/rrrvK3qsWsBSPZ3HbttX0ysVJk8IXQ21CXKtSMZMCGkAbSpTNsd2R\nzbE+SZoQ4hkZ62peqLM6fWJiQv/xP/5H/eIXv9Bvf/vbiuUtPZ6fTU2p7/TpMN3y6Ki0vBxW/PSn\nw1zB8f23KgMkmScBZKBV2RyRU0mSKNU0NHH9naVOolQoFPSnf/qnWlpaUk9Pj7797W9rbGysrFyl\nx/OT06fVd9994cN9+6Tvfz88/6d/CrUI3/rW6pdbNSyS4ZgA2gxNDMi+CWH9naWuTo/f+JeXl3XN\nNddUDFrWPZ7f+Z30ZQeALkKAgNpnN6xvZ6lnN0wayKx7PGNjUk9PeN7TE14DAMrQBwEdJZO+ELT3\nZ+LQIen996Xt26Xx8VaXBoBEHwR0sUz6QtDen4mJCenNN6VbbyVAAPKIJgYAAFCGAAEAAJQhQEBD\nrTcBU0t0YM4DAGgF+iCgYbLK4dC0SZrovAgAKwgQ0DBJJmDaSKaJotabpKlDEyYBQKN0RBODmT1s\nZpdaXQ7Upt4JmAqFgh5//PHsEkVt2xbyMYQCrZ2kiYRJALBGW9cgmNlWSbskPSGp8yds6EIHDhyQ\npIpTIq+nUCjoGzffrJuXlzUoaU51zvJYKIRJkS5eDEHC1NTaGoLiDI/FGgQSJm3ohRekjz6Srrii\n1SUB0AhtGyBEwcGHIjDoSKVNA2M1zlj4s6kp/d3ysjZLuiDp4B136L4nnsimeeHiRen06bXllfSz\nAwc0IumasTGaFxLo71cIvP6fGelf6LcB5E3bBgjufs7MdkQvX5e0tZXlQbn1Og/W2/9gRNLm6Plm\nSX96yy26JsH3q5ZpZCRMrby8HJaxGoKyfg5jYxoKH9BpcT302wByrW0DBEly93ckycxaXRSU2Kjz\nYLH/QfHzkZIb8kajEq4ZG9Olgwd12fKyLvX0hF/1dZZpjcXFlZt/xWBG4ua3kToycwJof20dINRq\n//792rx588rr0dFRjY6OtrBE+bVRDUExYVJpIJD4Jj40pMtOnJBmZnRZwl/w65ZpZibUHkhh+dWv\nhqaG3l59aWqqPJjh5rcx+m0ALTc9Pa3p6emV1xcuXMhs27kKEJ5++mmSNTXJejUERZXyJtTU9FBj\nzoR1yxS/mW3aFIKDUAj1nT5dMZjh5reBYmZOmmGAlin9ITw/P68f/ehHmWw7VwECmqdaDUFRtWaE\nJIFFPWX62dSU3v/hD7X9rrvUFy9T/Ga2bVsY0RC7+ZcFM9z8kiHxFZBbHZHu2czOSNri7hUDGtI9\nt5eNmhEymxmxdFvS2n4DU1NhtEKlGzwdEOu2sLA6zLG/v9WlASB1Ybpnd7+61WVAckn6JyQNGDb6\nLB6IzB04oGvi/QZi/QzKgoW8/fJtQcBzzz2r6Z7feKMpuwTQRB0RIKCz1NKMsF5tw0Y1EWWBiKSx\nKv0M1gQLeRuRwHBDAA3QEVMto7MU+ydMTExsmDuhUm1Dks+k8qmcf39sLNwcJyak//Jfws1SKg8W\n8jaNMtNEA2gAahDQEJVGMFSyXm3DyMiIenp6tLy8rJ6enrKaiKodJYvLdTolJtIp/RQYbgigAQgQ\n0FIbjYZI8v1EwyRrHZHQSdX2jLgA0AAECGi5ajf5mZkZLUeTGy0vL690dkw8CqK0BqCWG2enTZSU\nt06XAFqOAAFNVcsQx0rND4lmYiwUwqiFgwfDrIlpagCotgfQ5QgQ0DQ15UpQ5eaHycnJ9WdijDcN\nFKWpAaDafkN79ki33y5t397qkgBoBAIENE2aDI+lzQ8bDqGMNw0Upa0BoNp+XePjrS4BgEZimCOa\npnRYYrX5EQqFgiYnJ1UoFMo+23AIZbFpQAppnQ8ebFwHw0JBmpwMSwDIGWoQ0DRJRiwkaYbYcORC\nM5oGOmmUAwCkQICAptpofoQ0zRAVdrLuzTqTXBCdNsqhGTpl3ggAiRAgoK1UG7mQZXKnWjpKxr+3\npgyNGOVQeoPtpBsuNSpA7hAgoK2UNkNISnVDr6a0hmJqakrXXHPNusFH1aAiy6aM0hvs1NTa2R8b\nfcOtNxihRgXIHQIEtJ14M8SGwxprFK+h6Onp0cGDB7W8vLxu8FG12SPLUQ6lN9gf/rB5N9wsfv0z\nbwSQO4xiQFtLOvIhqfgoiMcee2xlpsZKyaAaVYYqO1kdfdHbK91119rXjbzhpkz2dPfd0vXXh+VK\n59CJCZoXgJygBgFtrd5cDdW2WZyy+fHHH98wLXWlMmTZLyLaSXmTxUZNGFn1UUj5639xUXr3XWnL\nltgxEBgAuUGAgLaXNDNkJevdyGsJPuJliPdJ6Onp0WOPPaaxsbHK36/lJl56g13vhlupWUBKFzAw\naySACggQkFt1z6lQRbxPwvLysh577DE9/vjj5dtvZM/+0maBqSnp8cfT74tf/wBK0AcBuVE6A2Ol\nzoWV1ivZyIazI8b7JBSV9WEoFFZv2GGFxG37iZT2WSjuoxH7AtCVqEFALlSqLag5G2TCX/zFpomp\nqak1oyCKfRgWX3lF27/6VV1+8eLql7LuaFjaLCCtrUHohlEEnTRPBNCBqEFALlQbiliat6FarUK0\nkcS/woeGhvStb31LJ06cWLP9QqGg75YGB/fe25ie/UND0n33rTYPdNMogmIwt2dPWJIPA8gcNQjI\nhXhtwaZNm7Rt2zZJsT4GUdPBl7ZtWzMPwm9+8xsVCoXUsyOW9mGYmZnRsYsXdUHSZkn/vGmTLj9w\noDkzIza5H0Ffn3T+fFg2HRMzAY3n7h3/kDQgyefm5hzd6+jRo75p0yaX5L29vX7y5MnwwcmT7r29\n7pJ7b68vHD3qBw8e9J6ensrrTkyEZQonT5703t5eH5T8jzdt8oWjRyuWIe32EeF8AhXNzc25JJc0\n4HXeW6lBQG6cPn1aF6Oq/TUzHpb82uw7fVrXXHNN2SRJlWZHrHW+g9Khk33F7zT6F2+3tcczNBNo\nOAIE5EalTonRB2VNByNS5XVjaprvIKbi0MmREamnR1peDsssOxF2a6IkhmYCDUUnRbSNdYcfJlCp\nU+LKL+upqTUd+CquW6LSfAfDw8Opy9cwKadKBoD1UIOAtpA2DXOpNb/eN/hlXW2SpGKzwrZYh8ai\nagmjNmyKmJkJtQdSWFZpYijbznpNB8XPtm1r/0RJ3dYEAuQAAQLaQtWMifVttOZ2/9JAZWpqSrOz\nsxXnO6j2nYrBTYIREqXb+dnUlPqqpXyulB769On2vAF3axMI0OFoYkBbaEjGxNLZBhNsszRQOX36\ndMX5Dtb7TsWskAnmKSjdzvuVUj6vrrz2s9OnV+dEaDc0gQAdiRoEtIVGZG1M09O9WkfH9XI2VO0c\nWak865ShdDvb77or1AxUqnVImYExS4cOSe+/L23fLo2Pr7NiG5QVQO3MwzwCHc3MBiTNzc3NaWBg\noNXFQYdLk8o5/h1JqQOd0n0vvvKK3v/hD7X9rrvUd/vtpSuHAEKSxsbS1x6k7B/wxS9Kb74p3Xqr\n9MYbjdkHgNrMz89rcHBQkgbdfb6ebREgABnKqrNlom1l0bZfxzZqChAANEWWAQJ9EIAMJeqPkNW2\nsmjbp38AgCoIEIAMZdnZcsNtpeiEWWEn9W8DQC7RSRGp2tyjL9KuXCJtZ8tKf4Nq21qzbj3TDccn\nkWrXIZLNwHUMVESA0OVSt5kztr2q9UY8VAoE1vsblG6r4rr33Vd7ISv9/SRpcrK7bpRcx0BVNDF0\nudRt5rRd16x4c9+zZ8+aKZtr+Rtk1seh9O83NRVulHv2hGW7TSfdKFzHQFUECF0udZs5bdc1q3Zz\nr+VvkFkfh9K/XyjU6jLBjfKFF6Sf/zwsOxbXMVAVwxyR6z4IqY+t+gZTH/N6TQm1lDOzY4ofi9S9\nVe0dcB0DSTEPQgkCBFSS5ZwE0QbrvolmHrBkKbpRLm7bpp+cPt2eZawHgQC6QJYBAp0UkVuZJ4BK\nkfyp1HodGBspUWAyNKSClG1Q1S46KbkV0CYIEJBbiXMkJN9gx+UUKBQKmpqaWpONcr2bfkOyaoaC\ntPbXe2lw99WvShcvdl9zClADAgTkVpYJoIq/wL80NaW+dvjlmeCGG29iKdropp95UFUsa6v7N8SD\nu02bQnAgpa4JAroBAQJyKV6lfl+aeQJKttVW1e4Jb7jx2oCijW76DcmqWcxIKbXuhhzP7LltW0hu\n1UE1QUArECAgd7K+oTes2j19gRLdcOO1AT09PXrsscc0Nja2YdmT9JMoFAr6m7/5qW644Qv61//6\nWvX3V11ROnhw9fWnPqVfuevVycmqAciG/SXqba7o66s5DTjQldy94x+SBiT53NycAxMTEy5p5XHw\n4EGfmJjwkydPptreyZMnvbe31yV5b29vuu2cPOk+MRGW9Tp50r23111y7+31haNHqx7fyZMn13xW\n+jpN2VfPx391yf2mm5arf3diIpQz9jgv+WCFc3ny5Ek/ePCg//6nP+3flPymnp7ycpYce8XzWelc\nJ/kekANzc3PF//sGvN57a70baIcHAUL3qnTDi9/Qe3p6vKenp76be5X91PDl7G9O0U1w4ejRxMFL\nqkCnQtkPHjwY/QcUAoTf/d256tuKfz/2eFTyb0p++OBBd3dfOHrU/3jTJt8VBRAu+QXJ37rjjrXn\nqzTgOHhwTTCwcPSo/3bTptXyHj0aPj94cO33JiYSn2qgk3RNgCBph6RZSZckvS3ptirrdX2AcOTI\nkVYXoelKb3h/+Zd/ueaziYmJ2M0sPCZacWMovallWIbS2pI/+ZM/qbjeyZMn/d577639XJTcWH99\n8OBKwFUMEKT/un7AcfKk/8ODD/qF2I2/+Pzjnh732E19uSSQKAuq4gFHT094ROv83aFDvuuyy9Z+\ntxgsxNfdtCkEDjnVjf8XlOrmc5BlgNC2Uy2b2VaFoOBaSU9KulLSj83splaWq11NT0+3ughNV9o3\n4IUXXlChUNDk5KQk6b777tPY2NiGUxMXCgVNfec7+s13vlNbDoJCISQ32ug7DZzOt3Tq5V/96lcV\nihn6ZDz//PMr7yUaoVDSf+DSFVfo/3z7bS0vL5etWuybUTz/hfg5GRrS/3XDDdoh6T5JByVFZ0OX\nLS/r3DPP6PJoVMEVkv6pfOPS44+H8kSdDf/Hgw9q8V/9K6lYlqUlze7dqzcvXdKF6Gv/8slPro5W\nWF6W9u1bHcEwNia98kqyv59UdlwVjzN+3hJuN9PvRrrx/4JSnIOM1BthNOoh6X5JH0vaH72+TaEm\n4a8rrNv1NQj33HNPw/dRVzV7A5w8edJv6ulZaa/+0pe+VLEKfb1yHz161K+//PKVau2Pe3qSNQHU\n2mxQYx+EhaNH/e/uvdcXEvzSjR9fpeugtJbh3nvv3fhvePKk+733rvk1/u3LL1+zHbM319QgHI01\nd/T09PjBgwfX/A2Kn93U0xPOs+T/vGmTf/+yy1ZqFM5Lvkvy/+MTn/CLl1++Zv8fX365/2znTn/r\nP/yHlb/XpdjnlyTfEX1/7yc/6f/foUMrf6Pfbtrk//Dgg2trFz75yap/v/g5La2pOrpes0616+Lk\nyVAbc/Bg9Wsgo6aoZvxf0O66+Rx0RRODpCeiAGE0et0XBQivVFiXAKHB/yAy6aiXfaFWbjQf9/T4\n/3bjjTVVoZ88edI3bdrk3yyt0k5S9V7SbPDrqC29tuJXDlwWjh5duQGel1aChOL6C0eP+q8PHvTD\nsRtwUaXroOrfrlrQUqHfwG83bfLB2LkdlPx/tTdW+iAcjTpKxs9/tUDt7w4d8rfuuMMX7rjDl2LN\nDo9G2y1u/1HJp83W/m0k/6hSM0T0uCfa1nei4OofHnzQv/2pT/lgSWByqeR7H9x770pAVnq+Spup\n1m2qqdScdPLkavNGsbmj0r+fKk1RR48e9XvvvdePRmUrXgPVAs5q10BmwX2WHW4btM+ODhDqPL9Z\nBgjtPMzxymi5EC3PRMurK6x7uST98pe/bHSZ2taFCxc0P1/XtNvrOnLkyJrq/CNHjmj37t0N21/C\nQoOZ2/QAAAiKSURBVK1WMS8v64qLF9XT06Pl5WX19PTos5/97Lrn5MiRI7p48aJel/SOQrX3pU9/\nWpd99rPSRufys58N6/7TP2lJ0t1//uf6iy98QX19fYmKvri4qK997WsrZf3bv/3ble++9eyz+oKk\n/x6t+9azz2r+ssv0ta99Tb+7vKwdCpHzgKSv/fmf6y9+9KOV71a6Di677DJNT0/r7bff1i233KLL\nLrtM86++Kn3ta+H89fRIf/u3Yfhf8bzG508YHdWv/u2/1el/9++k5WV9TtLLkr7pv9J/1+/o3P84\nodHRe/VXf/VXK+e/KH6tXHbZZeoz0+Y9e9Qr6X9K+nWsnHPRo7j9XknL7pqTZCXn758V/aMvcUHS\nbySdf/55nXv+eX0s6WuSnpc0t7ysvxgd1Vemp8u+688/r9+TdO755/Wjb35zzbX+wQcfrLmu/vAP\n/1Avv/xy5evMSkpqtvY6lcLzI0ek0n8/n/1s+FsU/yaf/az+3//8n1fm8Xj++ef1qU99Sr/3299q\nh6R5qfxvp/JrYL1rrWaLi9Wvm0ZJsc9G/3/YMBmc39h9sNI/kZq0bbImM3tC0kOSdrv7tJn1STol\n6VV3v71k3T+S9FoLigkAQDu6zd1fr2cD7VyDcErhh0NxCpbrouVshXX/m0IfhV8r/LgAAKAbXS7p\n9xTui3Vp5xqErZI+VGhLeU7S1yVtlXSdu7/XwqIBAJB7bTvM0d3PKcyDcFxhRMM/StpJcAAAQOO1\nbQ0CAABonbatQQAAAK3TMQGCme0ws1kzu2Rmb5vZbeusezhaL/4Yb2Z5G8nMHjazSxusk/h8daKE\n5yB314GZjZnZqehY5s3soXXWzeU1UOM5yN01IK1c//PR8ZyJRn1VWzd310GNx5/La6DIzJ6JjumD\nddZJdQ10RBNDrMPih5KelbRbYeKkYXd/p8L68wqdG19VGAnhkqbqHfLRatF52KUwTNzd/RPrrJf4\nfHWSpOcgWjdX10FsqO+HCse/U2GUzz53P1Sybi6vgVrOQbR+rq4BSTKzOyUdVpiK/mVJDyj8bXe7\n+3TJurm7Dmo5/mj93F0DRWbWrzAlhks66+6fqbBO+mug3pmWmvFQDdMuR59fkjTe6nJnfA62Rsf1\ncXGZ1fnqlEct5yCP10Hs77onen2zqs8umtdrIPE5yOM1EDvmcUlbon8TL0fn5KZuuA5qOf68XgOx\nY/uxpL+OjvGDKuukvgY6pYmhOAfCYrQszq7YX7pi9AtDkh6NVSf9oNEFbDRfHdUxLOnsBqsnPl+d\npJZzkNPr4FWF2pPD0evitH0fVlg3l9eAajgHOb0G5O4nPNSWfEbhuO+U9KxX/jWYu+ugluPP6zUg\nhaY2hZv9Ixusmvoa6JQAoZZpl/sVLpq3Je1VmFjpfjP7s4aWsAnc/Z0q/wmUquV8dZQazkHurgN3\nf8/dX3f382a2Q+EmeUbSoxVWz+U1UOM5yN01EOfuiwpNLMck7TOz0Qqr5fI6kBIff56vgWcl7XX3\n8xusl/oa6JQAofhrsRjxFA/sTOmK7v6au3/G3e+Kosy9Cr8ydjW+mG0j8fnKqzxfB2a2V+E/vI8l\n9XnluUFyfQ0kOQd5vQbMbKuZ3SxJHtrRH4g++kaF1XN3HdRy/Dm+Bh5W6HfwWtQPofj+1gqrp74G\nOiVASDztspndb2YvFy8ghWooqXI1bF7VMk11LuX1OjCzZyT9QNLL7j7k7heqrJrbayDpOcjrNSDp\ngKTZWE/869ZZN4/XQeLjz/E1sFOhZmBBq50Ur1IImkulvwZa3ckiYUeMeOe0HyhEPh9LulahN+Yp\nST+O1r0zWne+ZN0vtfo4MjwfZxTroFfhHFQ9X60uexPPQe6ug+g/hfjftfi4v1uugRrPQe6ugei4\n+qLjuqTQQW/luLrhOqjx+PN6DdwkaTT2uCTpg6yvgZYfaI0n5K3owArFP3B0Mj6WVIitO66QOba4\n7h2tLn/G5+KMpH+Jva50Diqer7w8Ep6DXF0Hkp6IjqX0cbRbroEU5yBX10DsuP6o5LiKozq65Tqo\n5fhzeQ2UnI9Lkv4x62ugI+ZBAAAAzdUpfRAAAEATESAAAIAyBAgAAKAMAQIAAChDgAAAAMoQIAAA\ngDIECAAAoAwBAgAAKEOAAAAAyhAgAACAMgQIAACgDAECAAAoQ4AAAADKECAAAIAyBAgAAKAMAQLQ\nBcysz8wuRY+PK3x+Z+zzuVaUcT2xsl0ys2vNbIeZvVrpWABk45OtLgCAptgRLV3SQoXPvxD7/HhT\nSpSQmd0cPfVoOSzpcPR8tvklAroDAQLQHb4Qe14pANgZe/5Wg8tSq3jZzikEB8Vg4e3mFwfoDgQI\nQHeI32RfrfD5Dq3edNuqBkFrg5utkuYlPatQE9JuZQVygwAB6A47Ys/X/OqOVeGbQpDQbr/K48HL\nrLv/QSsLA3QLOikCORcLACRJ7v5OySq3FD+StODu55tSsOT6Y88faVkpgC5DDQKQf/EAoFKV/HDs\neaIq+yjouGXDFSs74+5HEu7ntuJTSe7uP0m5TwA1IkAA8i8eAFRqPkjTQfGApDtTluesmb2asKYi\nPvqC/gZAExEgAPkXDwBOxT8wsx1aW4Wf6Cbs7l/PoFxJxDsoHmvSPgGIPghAN+jXaie/lRuumfVr\n7ZBBufvrzS3ahuKdK9tt+CWQa+buG68FoGOZ2RmF4YEWvXVc0hmFmoXifwAm6ZS7Dza/hJWZ2VZJ\nH2q1jNe5+3utKxHQXahBAPLvWa0GB1L4VX6bQnPDa1obOLSTYidIk3SW4ABoLgIEIOfc/VFJTyoE\nBB4tn4xqC66O3nNJf9+yQlZWnP/ARfMC0HQ0MQAAgDLUIAAAgDIECAAAoAwBAgAAKEOAAAAAyhAg\nAACAMgQIAACgDAECAAAoQ4AAAADKECAAAIAy/z/U6IrJJMhHGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa200bf91d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(1, 1)\n",
    "\n",
    "plot(ur_nd, mass1myr*100., \"ok\", ms=3, mew=0)\n",
    "plot(ur_rd, mass1myr*100., \"or\", ms=3, mew=0)\n",
    "axvline(2.22, ls=\"--\", color=\"b\")\n",
    "\n",
    "xlabel(r\"\\bf $u-r$\")\n",
    "ylabel(r\"\\bf last Myr mass fraction (\\%)\")\n",
    "\n",
    "tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3gAAAGsCAYAAACPXgVeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAOJgAADiYBou8l/AAAIABJREFUeJzt3cuTHdd9J/jvqVsPgDKbD/kRHeEI8QE5HPZiJGoUjnDY\nCw9pd3glDURJ/QeI0niNESX1HyCbw+Z6KGq88c6UByNTshRBSdbCveluQs11CyBpb7yx+RBokURV\n3TOLusUHWAB+CWTWBRKfT0RF4l78KvPkyZPnd383895qvfcAAABw69tYdwMAAAAYhwIPAABgJhR4\nAAAAM6HAAwAAmAkFHgAAwEwo8AAAAGZCgQcAADATCjwAAICZUOABAADMxKQFXmvtodbaudbasrX2\nQmvt4cLvfGsV/69Ttg0AAGBuWu99mhW3dleS11Y/zyT5fJL7k3yq9/7iFX7ngSTnk/Qkr/fePzpJ\n4wAAAGZoyit4X8hBofbN3vs3knwlSVstr+TpJN9axQEAADDAlAXeg6vly6vlS6vlA0cFt9YeTfJw\nkq9N2CYAAIDZ2pxw3XevloeF3aur5b1XiH8myZd7779o7coX8FprH0nye0n+OcnuCO0E4PawleTf\nJ/mvvfd/W1cj5DEAbsA1c9mUBd7rq+UDSV7Me4Xdq5cHttYez8HtnD9ZfQ7v8Pm7eu9vXBb+e0l+\nMn5zAbhNPJzk79e4fXkMgBt1xVw2ZYF3IQefpTss2A5v2Tx3ROwjObjid3i1rye5J8kLST5+Wew/\nJ8mPfvSjfOxjH7vixnd362+K7u3tleLeeeedUtwvfvGLUeOS5N/+rfZm88WLF0txv/zlL0txf/VX\nf5VHH320FHvp0qVSXLW/q+sbss5q3P7+filuuVyW4pLkWl9o9NOf/jR/9Ed/VF5n9QuShnyR0hTr\nXMf6kuRqdwL8wz/8Q/7wD//wmnHXa+z92dio3U0/ZF82N2vT/2KxKMWdPHnyyOd/+MMf5k//9E8/\n8Ny9917pRo4POnXqVCnut3/7t0tx//RP/5TPfOYzySqPrNE/J8mPf/zjq+axpD4PVueN6tyWJG+/\n/XYprpp3qnnsrbfeKsVVt/2Xf/mX+fznP19aX7WN1dcDSf31SDWumseGxFbnrGuNn8M8NmScrTPn\nVa1r21eb0//Lf/kv+YM/+IN3H1fzxFRfrlhRbWM1rpqfqnFJcscdd3zouaPy2D333FNaXzWPJcNy\n2Wc/+9nkKrlsygLv2Rx8YcoTrbVTee9LV55prd2f5MdJLvTe/yTJ4/ngZ/P+JgffvvnlI9a7myQf\n+9jHrtppUxQH1WT3+uuvXzsoyWuvvVaKS+oJ9I03Lr/gebQ333yzFHfixIn8+q//eim2mvCqSWxI\nAq0ew7ET6JgF3vb2du65557yOsdOikNib/UCb3t7O3ffffc1467Xugq8alxSL/CqcUclxSTZ2dnJ\nRz/6wS9E/rVf+7XSOn/zN3+zFPfggw9eO+iD1n1bZCmPJfVcVn1BPeSFd/WNwGreGfsNyOq2T548\nmd/4jd8ora+aG6uvB5L6MazGTfEG9lhvEBzmsSFFqALvyqp5LLk9C7yx81iSfOQjH/nQc0flsV/9\n1V8tra+ax5Jxc9lkBV7v/Y3W2kNJvp3ksRxc0ftc7/2VVYF3X5L9VeyLObiNM8m7A7r33n86VfsA\nAADmZsoreIeF26ePeP7lJFe8Xtp7n/QPsAMAAMyRQoqr+uQnP7nuJtw2fuu3fmvdTbhtXMdtEFyn\n3/md31l3E7jN/f7v//66m3DbkMeOjzx2fG7FPKbA46oUeMdHYjw+Qz70zI353d/93XU3gducAu/4\nyGPHR4F3fG7FPKbAAwAAmAkFHgAAwEwo8AAAAGZCgQcAADATCjwAAICZUOABAADMhAIPAABgJjbX\n3YDrtVwus1wur/j/Gxv12rUae7Xtvd/+/n4pbnd3txSXJG+99VYp7u233x51fdW4JHnnnXfWEpfU\n+3xvb2/UuOp2k6T3XoqrjrOxtztk20PWuS6ttVHj1mmKfVksFqW4zc1amhiy7bHns+pcWj2vj8u1\n8lgyfn6aYj4YO25IbvzlL385atwUufHSpUuluOp+T5Ebx44bksfGzo1T5Kd15bwh8+q6ctkU263O\ne9U8trW1Vd52dX+qc8CQ83XM156u4AEAAMyEAg8AAGAmFHgAAAAzocADAACYCQUeAADATCjwAAAA\nZkKBBwAAMBMKPAAAgJlQ4AEAAMyEAg8AAGAmFHgAAAAzocADAACYic11N2Aqi8Vi9HVubta6q7rt\n6vrWue0hbbx06VI5tqK1NnrsxkbtPY3q+oa0sfc++jrHVu2f5XI5cUuurNrGW8HY42yd58yQ41KN\nrc5nU7TxOOzv72d/f/+qMdU5uHo8h5y7W1tbo8ZN0f9j57zq+oa8xhg77wzZdvV4r3Muqqr2YzXX\nVuOS+v6MneenmNNvhTau6zXdkHWuc66oxN1c2Q4AAIDrpsADAACYCQUeAADATCjwAAAAZkKBBwAA\nMBMKPAAAgJlQ4AEAAMyEAg8AAGAmFHgAAAAzocADAACYCQUeAADATCjwAAAAZkKBBwAAMBOb627A\n9VosFlksFqOsa3d3d5T1HNrYqNXNQ9q/tbVVitve3h417tKlS6W4pN7G/f398jqr9vb2SnHVY9Na\nG3W7SX2/e+/ldY5tuVyW4tbZxuqxWaex21hd35DtVs+FseeeJDlx4kQpbmdnpxS3uVlLZWPljLFs\nbW2V+3cdqmNkXXFJ/ZiOPY6H5LHqvDqFal9WXwdV55gh+1zNJ2PnnSHrG3vbt0IeG/v10hTbHjuP\nJfU5oJqfqnHJuLnMFTwAAICZUOABAADMhAIPAABgJhR4AAAAMzFpgddae6i1dq61tmytvdBae/gq\nsY+21i6sYs+31r46ZdsAAADmZrICr7V2V5IXktyX5Ikkdyd5vrX2iSNi70/y7CrmW0l6kr9orX1p\nqvYBAADMzZRX8L6Qg0Ltm733byT5SpK2Wl7ukVXs4733P1v9bkvy+QnbBwAAMCtT/h28B1fLl1fL\nl1bLB46I/VGSP87BFb/koLhLktemaRoAAMD8TFng3b1aHhZ2r66W914e2Ht/JckrycHn9pJ8ZxX/\n9QnbBwAAMCtTFnivr5YPJHkx7xV2rx4dnrTWvpzk6STnk9zfe794pdgzZ87kzjvvfPfx6dOnc/r0\n6RttMwAz8bd/+7f53ve+9+7jixevmFLWQh4D4FquJ5dNWeBdyMGtloe3ZB7esnnuqODW2reSPJbk\n2d77f7zWyp966qmcOnVqjHYCMEOf+cxn8pnPfObdxxcuXMgPfvCDNbbog+QxAK7lenLZlF+y8uxq\n+URr7enV457kmdba/as/ifB8krTWHslBcdeTvN5ae3r189iE7QMAAJiVya7g9d7fWH2e7ts5KN4u\nJPlc7/2V1Z9FuC/J/ir88Fs0s4o9dP/q9ye1v79/7aAku7u7pbjlclmK29io19ebm7VDtb29Per6\nhrRxsViMuu3W2rWDVqrt3NvbG3V9Q/qnOi6q47Gq937toOuIHdO6trtuQ8Z4xTrP162trfK2T5w4\nUYrb2dkpxVXbWI27FVXnjeocOGSdVdXxOeQ4VcddNa56XlTjkvr+TDEPjp3zqvtdfb2U3J5551bI\nedX8NCSPjb3OsfNYUp8rTp48WYqr5rFk3Fw2abbrvb+Y5NNHPP9yksX7Hn89vlAFAADghkx5iyYA\nAADHSIEHAAAwEwo8AACAmVDgAQAAzIQCDwAAYCYUeAAAADOhwAMAAJgJBR4AAMBMKPAAAABmQoEH\nAAAwEwo8AACAmdhcdwNuBr33UeNaa6W4jY16fb1YLEpxm5u1Q7q1tTVqXJLs7u6WYyv29vZGXd8Q\n1WNTPdZJslwuS3HVYz32uB2iut/rbOPY66wevylUx+OtMKcMiR2yP7e7Kc616pivzgfVMVeNS5Lt\n7e1S3M7OTinu7bffLsUNyXdTzG9V1WMzdtyQc7c6zm6Ffqy2cYp8sq7XskOsa66o5rGknp+q256i\nHytkTwAAgJlQ4AEAAMyEAg8AAGAmFHgAAAAzocADAACYCQUeAADATCjwAAAAZkKBBwAAMBMKPAAA\ngJlQ4AEAAMyEAg8AAGAmFHgAAAAzsbnuBkyl9z76OheLxahxGxv1+rq1tpa46r4kydbWVjl2bGMf\n7/39/VLckGO4XC5H3XZ1n6c4F24F1f6uxg05F8ZWHWdDxmN1f7a3t0txm5v1dFLddnWeYpoxMvY6\np8iN1XVWx2c1jw0Z72PP1dU5a4ixXzsMOYZj73d1fUPml7Hnoiny99j9M+QYVo09fsbOY8n4c8AU\n/Vja7lq2CgAAwOgUeAAAADOhwAMAAJgJBR4AAMBMKPAAAABmQoEHAAAwEwo8AACAmVDgAQAAzIQC\nDwAAYCYUeAAAADOhwAMAAJgJBR4AAMBMKPAAAABmYnPdDZhKa60cu7FRq3MXi0UpbnOz1q1bW1ul\nuCHbrsaNvc9Jfb+rqm1Mkt57KW65XF5vc45NtY37+/uluGrfJMPOm5tdtX+q/T2kb8bux+q5MOSc\nqcZWz+sh5381tjr/VPt7TuP7cmPnpyTZ29sbdZ1j57Fk/HNjin6sWmceG/scGpJ3qrHV8Vhd35D5\nYMiYrKgewyH9OPa4qJ4zU8yrY7+WneL1dnWdQ8bOmOehK3gAAAAzocADAACYCQUeAADATExa4LXW\nHmqtnWutLVtrL7TWHh4jFgAAgA+brMBrrd2V5IUk9yV5IsndSZ5vrX3iRmIBAAA42pRX8L6QpCf5\nZu/9G0m+kqStljcSCwAAwBGmLPAeXC1fXi1fWi0fuMFYAAAAjjBlgXf3anlYrL26Wt57g7EAAAAc\nYco/dP76avlAkhfzXrH26g3GJknOnDmTO++8893Hp0+fzunTp2+kvQDMyHe/+91897vffffxm2++\nucbWfJg8BsC1XE8um7LAu5CDz9Ed3mZ5eBvmuRuMTZI89dRTOXXq1AjNBGCOPvvZz+azn/3su48v\nXLiQ73//+2ts0QfJYwBcy/Xksilv0Xx2tXyitfb06nFP8kxr7f7W2oXW2vPXip2wfQAAALMyWYHX\ne38jyUNJfpbksST/kuSR3vsrq5D7Vj+VWAAAAK5hyls003t/Mcmnj3j+5SSLSuxxWCwW1w5KsrW1\nNWrc5ma9+6vrHHtftre3S3FJ0nsvxbXWSnHVfZnC/v7+Tb/tjY3a+zPL5bK87eqxqR7ratwQY4+f\ndbZxXXFDYqvz1JDzdex1VuOq58ytqLpvQ45Tdf6vzjHV417NT0Niq9uuxu3s7JTiknqfDzl/q8bO\nE9X8NCTvVGOr+zLFnD62dc5F1f6ptnGKvDN23JDX22PPFUPm3DFz2XyzHQAAwG1GgQcAADATCjwA\nAICZUOABAADMhAIPAABgJhR4AAAAM6HAAwAAmAkFHgAAwEwo8AAAAGZCgQcAADATCjwAAICZUOAB\nAADMxOa6G3Az2NysdUPvvRS3s7NTinv77bdLcUmyvb09atzJkydLcdV9TpKNjdr7BYvFohS3u7tb\n3nZVa60Ut7+/X4pbLpflbVf7sto/1TZW44ao7veQ/qmqHsPqeJzCutpY3e4Q1flxa2urvM6x57Nq\nG6txczbkOFVV55hq/1ePe1LPt5cuXSrF7e3tleKGnLtj57Ip8nJ1rq5uu9qPQ2LHbuMQY2/7Vshj\n1W1PkXeqxn7dmdTnyLHzWDJuLnMFDwAAYCYUeAAAADOhwAMAAJgJBR4AAMBMKPAAAABmQoEHAAAw\nEwo8AACAmVDgAQAAzIQCDwAAYCYUeAAAADOhwAMAAJgJBR4AAMBMKPAAAABmYnPdDbgZbGzU6tzN\nzVp3bW1tleK2t7dLcUNiq9ve29srxZ08ebIUlySttXLs2JbL5ajrq46JIdut9nm1H8eOG6L3Pmrc\nkH6sHptqXNWQflznuVBV7Z/FYlGK29nZKW+7GjtkjqRmyNis5pPd3d1R11eNS+pj6dKlS6W46jw9\nhf39/VLc2PluiOq2p5gv19k/Y+eTqur8m4yfd6bY5+prguq+jP36PannnZs9j7mCBwAAMBMKPAAA\ngJlQ4AEAAMyEAg8AAGAmFHgAAAAzocADAACYCQUeAADATCjwAAAAZkKBBwAAMBMKPAAAgJlQ4AEA\nAMyEAg8AAGAmNtfdgFtJ733UuCFaa6W4xWJRitvcrB36vb29UtyQbY8dl9T3p3psqv29XC5LcUPW\nWVU9NkO2O8XYXZex+3vI+qqxY88pGxv19+y2trZKcdvb26W4kydPlrddXeeQOYDxVcdxNa46Pocc\n97FzXjXu0qVLpbhk/P2utnGIseeiseffpN6PQ/LyzW6KvFM1JJ9UVY9Nddtj57EkueOOO0Zd57ry\nmCt4AAAAM6HAAwAAmAkFHgAAwExMVuC11h5qrZ1rrS1bay+01h6+SuyjrbULq9jzrbWvTtUuAACA\nuZqkwGut3ZXkhST3JXkiyd1Jnm+tfeKI2PuTPLuK+VaSnuQvWmtfmqJtAAAAczXVFbwv5KBQ+2bv\n/RtJvpKkrZaXe2QV+3jv/c9Wv9uSfH6itgEAAMzSVH8m4cHV8uXV8qXV8oEjYn+U5I9zcMUvOSju\nkuS1q23gzJkzufPOO999fPr06Zw+ffq6GgvA/Jw9ezZnz5599/HFixfX2JoPk8cAuJbryWVTFXh3\nr5aHhd2rq+W9lwf23l9J8kpy8Lm9JN9ZxX/9aht46qmncurUqRGaCsAcXV4wnT9/Ps8999waW/RB\n8hgA13I9uaxc4K2++OSLObid8kP/vXr+sd77i0leXz3/QJIX815h9+oRv3u4/i8neTrJ+ST3995v\nrrdaAQAAbnLlAq/3/mSSJ4vhF3JQ9B3eknl4y+a5o4Jba99K8liSZ3vv/7HaJgAAAN4z1ZesPLta\nPtFae3r1uCd5prV2/+pPIjyfJK21R3JQ3PUkr7fWnl79PDZR2wAAAGZpks/g9d7fWH2e7ts5KN4u\nJPlc7/2V1Z9FuC/J/ir88Fs0s4o9dP/q9wEAACiY6ktWsvos3qePeP7lJIv3Pf56rvGFKjeL/f39\nawcl2d3dHXV9SdJau3ZQks3N2iFdLBbXDhoQlyRbW1uluOp+L5fL8rZ7P+qjoR9W7cfq+oa0sRpb\njasemyHjbOw2bmzUbhIY0o9jq7ZxiOo4q6qOxyH7Uj1fT548WYo7ceJEedu3wrigPu7WFZfUz7Xq\nmKvm0O3t7VJcMn4+maJ/xm7jkDmwGru3t1eKm2Kc3QrGzjtjry+pn4fV1zdj57Ek2dnZKcVV+2dd\neWyqWzQBAAA4Zgo8AACAmVDgAQAAzIQCDwAAYCYUeAAAADOhwAMAAJgJBR4AAMBMKPAAAABmQoEH\nAAAwEwo8AACAmVDgAQAAzIQCDwAAYCYUeAAAADOxue4G3Ep676W45XJZitvYqNfXi8WiFLe5WTuk\n1fVV45J6/5w4caK8zqpqX1b7Z39/f9S4IbHV8VNd35BxVt322GO8tVaKS8Y/D6dQ3e8hx2Zs1XNh\na2tr1PUl9eO9zv6hbuzjOWQ+qMaOPd6r89AUhvRPNYePnfOmOHeH7HfFFDmi2sZq3K3QxiHHeuxj\nOPZr46Q+B0wxn41J9gQAAJgJBR4AAMBMKPAAAABmQoEHAAAwEwo8AACAmVDgAQAAzIQCDwAAYCYU\neAAAADOhwAMAAJgJBR4AAMBMKPAAAABmQoEHAAAwE5vrbsCtZHt7uxS3t7c3+rZ776W4S5culeJa\na6W4jY36ewCbm7XhVN12NW7Itqv9s7u7W952VXV/9vf3S3FDjk1VdZxVx3h1nxeLRSluiGo/Vvd5\niLH3e8i5UFU9Z7a2tkaNS+pz6ZB1Mr7quKuOpepxf/vtt0txSX0erJ5r1X2pzi9JfRxX96XaxqSe\n8955551SXHVMDHkdVD021f6ptnG5XJbiknqeWGcbq8bOJ0PWN3b/jJ3HhsRW46rz3thcwQMAAJgJ\nBR4AAMBMKPAAAABmQoEHAAAwEwo8AACAmVDgAQAAzIQCDwAAYCYUeAAAADOhwAMAAJgJBR4AAMBM\nKPAAAABmQoEHAAAwE5vrbsCtZGOjVg9vbW2Nvu29vb1SXLWN1bjFYlGKu1VU97u1NmpcMv4xrG57\nyDHsvZfidnd3R11fdZ+nUG3jFOusHsMpztdq7OZmLU1U45Jke3u7FLezs1NeJ+tTPZ7L5bIUN2Qs\njT1XV7ddnc+T8efBIXln7P6pzv1DjJ3Lqv24v79fikvqx3uKvLwu1fN1iLGPYfV8rc5RSb2N1XUO\n2faYXMEDAACYCQUeAADATExW4LXWHmqtnWutLVtrL7TWHi78zrdW8f86VbsAAADmapICr7V2V5IX\nktyX5Ikkdyd5vrX2iav8zgNJHksy/gdiAAAAbgNTXcH7Qg4KtW/23r+R5CtJ2mp5JU8n+dYqDgAA\ngIGm+hbNB1fLl1fLl1bLB44Kbq09muThJI/m6kXgu86cOZM777zz3cenT5/O6dOnr6uxAMzP2bNn\nc/bs2XcfX7x4cY2t+TB5DIBruZ5cNlWBd/dqeVjYvbpa3nuF+GeSfLn3/ovqV8w+9dRTOXXq1PW3\nEIBZu7xgOn/+fJ577rk1tuiD5DEAruV6clm5wGutfTXJF3P0Z+Ta6vnHeu8vJnl99fwDSV7Me4Xd\nqx/6xdYeX/3uT1afwzt8/q7e+xvV9gEAANzuygVe7/3JJE8Wwy/koOg7LNgOb9k8d0TsIzm44nd4\nta8nuScHX9Ly8Wr7AAAAbndT3aL5bA6+MOWJ1tqpvPelK8+01u5P8uMkF3rvf5Lk8Xzws3l/k+S1\nJF+eqG0AAACzNEmB13t/o7X2UJJv5+BPH1xI8rne+yurAu++JPur2BdzcBtnkmT1Gbzee//pFG0D\nAACYq6mu4B0Wbp8+4vmXkyyu8nuT/fF1AACAOVNMAQAAzMRkV/AYV/XPR2xs1Gr2zc3aoe/9qC9N\nvbFt7+/vj7q+JFksrnhReNJtV4/LEHt7e6W46j5Xj3VS35/qtpfLZXnbY6sew7HHbVI/hlXVYzjk\nWFeP4dbWViluyPk6JJb5GDs/JePPWVXV82KI6hwzJO+MnSfeeuutUtwUrx3GPtZDxkQ1trrf1X0Z\nci5U+7Gal6t5bMixru7P2ON2ity4zteJFbIsAADATCjwAAAAZkKBBwAAMBMKPAAAgJlQ4AEAAMyE\nAg8AAGAmFHgAAAAzocADAACYCQUeAADATCjwAAAAZkKBBwAAMBMKPAAAgJnYXHcD5qj3XopbLpej\nr3Nra6u8zorFYlGO3d3dLcXt7e2V4vb398vbrsZW+3Fjo/bex5D+HnK8K1prpbjqvgyxuVmbOqr9\nXY1Lxt/v6hivjttk/P2utnHIsa6usxpXPS5DTLFOxjd2zqvOL0mys7NTipviPK+aYh6sqp5DYx+b\nIftS7fPqOqeYN4aMyYpqG4dstzr/j/16achrm7HzcrV/hryWrR6bsePG5goeAADATCjwAAAAZkKB\nBwAAMBMKPAAAgJlQ4AEAAMyEAg8AAGAmFHgAAAAzocADAACYCQUeAADATCjwAAAAZkKBBwAAMBMK\nPAAAgJnYXHcD5mhjo1Y3t9bK61wsFqW4kydPluK2trZG3e6Q2EuXLo0alyR7e3uluOVyWV5nxeZm\n/RSqHu/9/f1S3Nj7MoXe+6hxSb0fq+Oxur4h52t1XFTXWZ1ThozHseepKeazahtZr+r5W52zhuSd\nj3zkI6W46rz6zjvvlLddNfY8uLu7W952db+r266+dhhy7laPd3X8VOOqfZOMP1ePvb4p1lk9LkPO\n12qOqrZxitfbVWO3cWyyJwAAwEwo8AAAAGZCgQcAADATCjwAAICZUOABAADMhAIPAABgJhR4AAAA\nM6HAAwAAmAkFHgAAwEwo8AAAAGZCgQcAADATCjwAAICZUOABAADMxOZUK26tPZTk20k+meRnSb7W\ne//JFWLvSvL/JHl49dQzvfevT9W2qW1u1rp1f3+/vM7eeymutVaKq7ZxuVyW4pL6/qxz29V+3Nio\nvfdRXV+SbG1tleKq+7K7u1uKG9KP1f2uxu3t7Y0aN0S1jWOfM0PWOXZ/T9HG6rjd3t4ub7u6ziH7\nw/pUx+disSjFVcdmUh9L1Xmwuu0h+bs6v1XH+5C8U40dO27IuVuNHTufVMdtMmxMVoz9mi4Ztj8V\n1XNrSBurc8DYc0p1X5JkZ2enFFfNeevKY5NcwVsVbC8kuS/JE0nuTvJ8a+0TV/iVnyU5neTpJP89\nyVdba1+aom0AAABzNVVZ+YUkPck3e+9PtdZ+nORHSb6S5M/eH9haezjJ/Uke773/59Vz/673/ouJ\n2gYAADBLUxV4D66WL6+WL62WDxwR+9Bq+Settf9r9e+/yUGReEVnzpzJnXfe+e7j06dP5/Tp09fX\nWgBm5+zZszl79uy7jy9evLjG1nyYPAbAtVxPLpuqwLt7tTws7F5dLe89Ivajq+X9Sb6c5PNJHm2t\n/Xnv/RtX2sBTTz2VU6dOjdFWAGbo8oLp/Pnzee6559bYog+SxwC4luvJZeUCr7X21SRfzMGtlx/6\n79Xzj/XeX0zy+ur5B5K8mPcKu1eP+N2sfvcveu9/2Vr7TpLXkjyS5IoFHgAAAB9ULvB6708mebIY\nfiEHRd/hLZmHt2yeu0rsPavHH33f8wAAABRN9Xfwnl0tn2itPb163JM801q7v7V2obX2fJL03r+d\ngyt+T7TWns3Bt2/2JH89UdsAAABmaZICr/f+Rg6+POVnSR5L8i9JHum9v7IKuW/1c+ihJD9O8rlV\n7Jd77//fFG0DAACYq8n++t7qs3ifPuL5l5MsLnvulST/Yaq2AAAA3A6mukUTAACAY6bAAwAAmAkF\nHgAAwExM9hm829nGRq1u3tysd39rrRTX+1F/pvD61zekjdX9Xue2t7e3S3F7e3uluOVyWYobYn9/\nvxRX3echbRz7GE5xLlQtFotrBw2IG9KPY4+LKfqxut/Vc+bkyZPlbVdjt7a2yuvk5rezs1OKq86/\nSf1cGzuVamUeAAARvUlEQVTvVM+fIbHV8T5k29X9qfZ5NT9VX4sk4/f5FPmkauz+meJ1UNUUr2/G\nzmVj57Gknp9OnDhRiltXHnMFDwAAYCYUeAAAADOhwAMAAJgJBR4AAMBMKPAAAABmQoEHAAAwEwo8\nAACAmVDgAQAAzIQCDwAAYCYUeAAAADOhwAMAAJgJBR4AAMBMbK67AbezxWJRjm2tleKWy2Uprvde\nitva2irFJfX92diova9QjRsaW7G/vz9qXFI/NtW4zc3a6TukjdV+HHs8Dmlj1djjcYjq+VWNq/b3\nkDmlem5vb2+X4k6cOFHedjV2yP4wH0POyeo5NPa2h4zNamz1XBuiuj/VOXh3d3fU9Q2J3dvbG3V9\nU4yzsXPZFOOsui/V/D3E2OdX9ZwZ8lp2Z2enFHez5zFX8AAAAGZCgQcAADATCjwAAICZUOABAADM\nhAIPAABgJhR4AAAAM6HAAwAAmAkFHgAAwEwo8AAAAGZCgQcAADATCjwAAICZUOABAADMhAIPAABg\nJjbX3YDb2cZGvb6uxvber7c5R9rf3y/Hbm1tleLGbmOSLBaLUbe9XC5HjUvqfVldZ2utFHfp0qVS\n3JB1Vvu7asiYmGL8VFT7ZmjsmOsbclw2N2vTf/W83t7eLm97SCy3nyG5cexzrTpPV8+fIbHV/R7S\nP9U5obrf1flgb2+vFJfUc1T1WFf7Z8h8Wd32kNcEY253iGobxx63Q2Krx6Y6HqtxSXLixIlS3M7O\nTnmd6+AKHgAAwEwo8AAAAGZCgQcAADATCjwAAICZUOABAADMhAIPAABgJhR4AAAAM6HAAwAAmAkF\nHgAAwEwo8AAAAGZCgQcAADATkxV4rbWHWmvnWmvL1toLrbWHrxL75dba+VXs+dbaV6dqFwAAwFxt\nTrHS1tpdSV5I8lqSJ5J8PsnzrbVP9d5fvCz2k0meTnIuyeNJ/jjJE621C733s1O0b85aa6W4xWIx\n6vqSZGdnpxS3uVkbdsvlsrzt3nspbm9vr7zOiiH9s7Ex7vsp1WNY7e+kvj9D9nvs9VWP9diGjMfq\nsR47rjomkvq42NraKsVtb2+Xt11dJ1zL2HNW9VwbMobvuOOOUtzY+Smpz1vrmvuTYTmqYop9GXuu\nnqKNQ3JURXVfhhy/ao4a+zwccr5Wc9nNnsemuoL3hSQ9yTd7799I8pUkbbU8Sk/y4yTfTvLM6rlP\nT9Q2AACAWZrkCl6SB1fLl1fLl1bLBy4P7L3/j9bak0m+tvpJDq78fetqGzhz5kzuvPPOdx+fPn06\np0+fvpE2AzAjZ8+ezdmz790IcvHixTW25sPkMQCu5XpyWZvidqfW2tNJHkvyqd77i6tbNl9Lcq73\n/unLYh9J8nwObun8eg5u5/xfkzzce//FEes+leTnP//5z3Pq1KnR23672N/fL8W99dZb5XXu7u6O\nuu113qJZbeOQ86e6P9W4an9X4xK3aF6NWzSv7uTJk+Vtv7+oGWPbVefPn8/HP/7xJPl47/38qCsf\nQB47ftW5vzpfXrp0qbzt6jrXeYtmNedV44bknWobq/1TjZviFs1qflrnLZrVY3i73qJZzWXrymNJ\nLZeVj8rqi0++mIPbKT/036vnH1t9xu711fMPJHkxyb2rx68e8buPrn73a733n7bWzuWgGPz2ansA\nAAAUlAu83vuTSZ4shl/IQdF3eEvm4S2b567yO/eslh9dLR+qtg0AAIDpvmTl2dXyidXtms/m4Crd\nM621+1trF1prz69ivpODYvBvWmvP5uBWzb56HgAAgKJJCrze+xs5uAL3sxx8Fu9fkjzSe39lFXLf\n6ie995/k4DbNc0k+l+Rfkzzee/9PU7QNAABgrqb6Fs2sPov3oT910Ht/OcnisufOJvE37wAAAG7A\nVLdoAgAAcMwUeAAAADOhwAMAAJgJBR4AAMBMTPYlK9zcWmujxiXJ5mZtOG1tbZXi9vf3y9ve29sb\ndZ1D9rtqsVhcOyj1fqyurxo3heVyOfo6NzZq70uNfQx776OuLxn/PKz2zZDYKc6FKfoSxlAd79U8\nltTn4Gp+euedd8rbrp5r1fmgur5qHkvqeaK6zurrgSnyU1V1nE2Rv8fe73XmnSHbrqqO8Zs9j7mC\nBwAAMBMKPAAAgJlQ4AEAAMyEAg8AAGAmFHgAAAAzocADAACYCQUeAADATCjwAAAAZkKBBwAAMBMK\nPAAAgJlQ4AEAAMyEAg8AAGAmFHgAAAAzsbnuBnBzWywW5diNjXHfL+i9j77t7e3tUtxyuSxve2zV\nfdnf3y/FbW1tlbdd3e/qsdnd3S1vu6o6JqvHurVWitvb2yvFJfX+Gbu/h4zb6n5Xx9mQ/lnn+cXt\nqTqvbm7WXhYNGcPV2Grc2Lk2qeeJahur80ZSn9+qcdVtT9HGav9U599qHkvGf+0w9nEZEjvFOKua\nYvysgyt4AAAAM6HAAwAAmAkFHgAAwEwo8AAAAGZCgQcAADATCjwAAICZUOABAADMhAIPAABgJhR4\nAAAAM6HAAwAAmAkFHgAAwEwo8AAAAGZic90NYD02Nmq1/fb2dnmdy+Vy1LghWmuluMViUYqrtrHa\nj0m9jb33Utz+/v6o2x2y7Wpctb+r+zJknZubtemtegyHHOvq/lT7sToeh/RjdVyMHQfrMOT8rZhi\nvG9tbY2+zurcMfb+VOe2pD5vVeOqx3pIf1f7cW9vrxRX7e8hbazu99ivg6Y41mOPxynO15s957mC\nBwAAMBMKPAAAgJlQ4AEAAMyEAg8AAGAmFHgAAAAzocADAACYCQUeAADATCjwAAAAZkKBBwAAMBOT\nF3ittcdba8trxDzUWjvXWlu21l5orT08dbsAAADmZrICr7V2V2vt0SR/kaRfLS7JC0nuS/JEkruT\nPN9a+8RUbQMAAJijSQq8VdH2WpK/LoR/IQcF4Dd7799I8pUkbbUEAACgaHOKlfbe32itPbR6+PdJ\n7rpK+IOr5cur5Uur5QNTtI1hNjbq7wEMia3o/YoXfj9kd3d31Ljqtvf390txU1gur3rn87sWi0V5\nndX9HnJsxtZaK8VV97s6boeM7+q2q8dwCtVtV/t7c7OeToaMSThOY+expD7eq/Nq9ZxMkr29vVJc\nNZetM0eMPWdN8fpm7P2eIu+MPR6H5LFqnhi7H4ecM1tbW6Ovcx0mKfCSpPf+YlLqgLtXy8PC7tXV\n8t6r/dKZM2dy5513vvv49OnTOX369PCGAjBLZ8+ezdmzZ999fPHixTW25sPkMQCu5XpyWRvwjtFX\nk3wxR3+erq2ef+ywsHvf772a5K7e+5FvGbTW/iLJV5N8vvd+trV2f5ILSX7Ue/8PR8SfSvLzn//8\n5zl16lSp7dya1nkFr/qu561wBWaIsd+drb4rPORKaPVds52dnVLcFO/MVo/NrTB+qv29vb1d3vav\n/MqvlOKq76RWnT9/Ph//+MeT5OO99/OjrnwAeYyjVOeYd955p7zOW+EKXnUuqubvqiFXx6r7M3Yb\nh8yrY98ZMcUVvHVdAR5yta3a53fccUcprvpaZIhKLitfweu9P5nkyZHa9n4XclAgHt6SeXjL5rkJ\ntgUAADBbx/538Fpr97fWLrTWnl899exq+URr7enV457kmeNuGwAAwK3suAq8y6+13rf6Se/9jSQP\nJflZkseS/EuSR3rvrxxT2wAAAGZhsi9ZOdR7v/eyxy8nWVz23ItJPj11WwAAAObs2G/RBAAAYBoK\nPAAAgJlQ4AEAAMyEAg8AAGAmFHgAAAAzocADAACYCQUeAADATEz+d/DgerXWyrHb29uluP39/VLc\ncrksxb399tuluCHr3Niove+yWCyuHZR63wzZ9u7ubinu0qVLpbhq3yTJ1tZWKa66L9W4IeOx916K\nq+53dX3VfRmyzmrc5mY9nVTXCbeT6hwzZE6v5rxqXHVOH2Ls+bKaI6o5NKnnvKop8k41duxtD8nf\n1WM9ZL/HNmRc3MxcwQMAAJgJBR4AAMBMKPAAAABmQoEHAAAwEwo8AACAmVDgAQAAzIQCj6s6e/bs\nuptw2/jBD36w7ibcNr73ve+tuwm3jeeee27dTeA2J48dn7/7u79bdxNuG9///vfX3YTbxq2YxxR4\nXJXEeHx++MMfrrsJtw2J8fjciomReZHHjo83Ko+PYvr43IpvCivwAAAAZkKBBwAAMBOb627AddhK\nkn/8x39cdztuCxcvXsz58+fX3YzRvP3226W4d955Z9S4JFkul1f9/zfffDMvv/xyNjZq77ssFotS\n3NbWVikuSVprpbi9vb1SXLV/eu+luCTZ3KxNW9vb21f8vzfffDMvvfRSkpT7ewrXGhOHqv0zZF+q\n66zGXem4vPnmm7lw4cIHnrvasbmeuKr35Y36STENeewYzS2PVeeNpJ7zLl26VIrb3d296v8f5rEh\nqnNMdb+rOW/IfFnNedW46ravlu/en8eS+muC6rarrweGjMfqsa5uewpH9c/Fixc/lMd2dnZK6xs7\njyW1XNaGvLC6GbTW/rckP1l3OwC4ZT3ce//7dW1cHgNgBFfMZbdigfeRJL+X5J+TXP2tJAB4z1aS\nf5/kv/be/21djZDHALgB18xlt1yBBwAAwNF8yQoAAMBMKPAAAABmQoEHAAAwEwo8AACAmbgpC7zW\n2pdba+dba8vV8qtXiX2otXZuFftCa+3h42zrrW5gX39nFff+ny8dZ3tvZZf19QuttU9eJda4vgED\n+9q4vk6ttcdba8vLniuP3bmPc7ns+Mhlx0cuOz5y2fTmmsduum/RXA3ec6ufv07yx6ufR3vvZy+L\nvSvJa6ufZ5J8Psn9ST7Ve3/xONt9KxrS16v480l6kh8laat//806/57UraK19kiS55OcT/L/Jnk0\nyQNJHuy9v3JZrHF9A4b09SreuB5oNUb/OMmzSXrvffG+50tjd+7jXC47PnLZ8ZHLjo9cNq3Z57He\n+031k+STSfaT/HmSu5J8LskyyZ8fEfvYKvbM6vHDq9j/e937cSv8DOnrVfwyyZfW3e5b8SfJ06u+\n/sT7+t64XnNfr/7fuB7Wv3et+mz/cPm+/yuP3bmPc7ns5uzrVbxz/vr7Wi67Cft69f/Gdb1vZ5/H\nbrpbNHvv/yPJk0m+loOK+Dur5beOCH9wtXx5tXxptXxgyjbOxZC+bq3dv/rn19936f/pY2vsre9C\nDt5R+z9W7zb/p9XzDx0Ra1zfmHJfG9fD9d7fyEFffirJ65f995CxO+txLpcdH7nsWMllx0cum8jt\nkMduugJvdUn68SQvJHkkBxP0S0lePSL87tXysEMPY+6dso1zMbCvH8hBwnwhyZdzcCvMY621//N4\nWntr670/mVWf5aAPP3GVcOP6Bgzsa+P6OvTeX+xH334yZOzOepzLZcdHLjs+ctnxkcumNfc8trnu\nBhzh0RzcN/y13vtPW2vncjBov53ki5fFHlbdDyR5Me916FGTOh9W7uve+0+SfPTw8Sr2XA7uX/7P\nx9biW1jv/dOttcMJ+qM5uE/+8neOEuP6hlX72rge3ZCxO/dxLpcdH7nsGMllx0cuW4tZ5LGb7gre\n+9yzWh4O2KMu/x9evj68FHp4qfTchO2ao2v2dWvtsdbas+/7BqfD2NembtwctNY+2Vp7LMkbq3eM\nPrX6r/92RLhxfQOG9LVxPbohY/d2Gedy2fGRyyYmlx0fuWxt5pHH1v0hwCM+sHj4AcVlDr7Z5tUc\nfIDxmzn4ZpoLSZ4/4kOST78v9r5178et8DOwrw8/tH7+sr7+o3Xvx63wk/c+HP3qqv+WSfaS/Dvj\neq19bVzfWF+/mg9+OP2KY/d2G+dy2U3b1875G+truezm7Gvj+vr7eZZ5bO0de4XOPp3kv6866ed5\n79tp7l899z/fF/uJ98X+T4N50r7+0irmsK//93W3/1b6uaz//luS/+UqfW1cH19fG9fX38+vJtm7\n7Lkjx+7tOM7lspu2r53zN9bXctnN2dfG9fX18Szz2E33d/AAAAC4PjfzZ/AAAAAYQIEHAAAwEwo8\nAACAmVDgAQAAzIQCDwAAYCYUeAAAADOhwAMAAJgJBR4AAMBMKPAAAABmQoEHAAAwE/8/hu75EtJ1\n4qEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa200c1f810>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ssag_cat = genfromtxt(\"{}/{}\".format(os.path.expandvars(\"$phd/inputs/tta_cb13/ssag/\"),\"true_params.log\"), dtype=None, names=True)\n",
    "\n",
    "mask_nd = array([name.endswith(\"nd\") for name in ssag_cat[\"Galaxy\"]], dtype=bool)\n",
    "mask_rd = ~mask_nd\n",
    "\n",
    "fig, (ax1, ax2) = subplots(1, 2, figsize=(10,5), sharex=True, sharey=True)\n",
    "\n",
    "xl = xlim(8, 10.2)\n",
    "yl = ylim(-1.1, 0.4)\n",
    "\n",
    "H, xe, ye = histogram2d(ssag_cat[mask_nd][\"rFWLAyr\"], log10(ssag_cat[mask_nd][\"ZZsun\"]), bins=(40,30), range=((xl[0],xl[1]+diff(xl)/40.),(yl[0],yl[1]+diff(yl)/30.)))\n",
    "H = gaussian_filter(H, sigma=1)\n",
    "H[H==0] = 1\n",
    "ax1.pcolormesh(xe[:-1], ye[:-1], log10(H.T), cmap=cm.Greys, vmin=0.0, vmax=3.0, edgecolors=\"face\")\n",
    "\n",
    "H, xe, ye = histogram2d(ssag_cat[mask_rd][\"rFWLAyr\"], log10(ssag_cat[mask_rd][\"ZZsun\"]), bins=(40,30), range=((xl[0],xl[1]+diff(xl)/40.),(yl[0],yl[1]+diff(yl)/30.)))\n",
    "H = gaussian_filter(H, sigma=1)\n",
    "H[H==0] = 1\n",
    "ax2.pcolormesh(xe[:-1], ye[:-1], log10(H.T), cmap=cm.Greys, vmin=0.0, vmax=3.0, edgecolors=\"face\")\n",
    "\n",
    "fig.tight_layout()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}