{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-14T21:18:49.274149Z",
     "start_time": "2020-03-14T21:18:39.580491Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import py21cmfast as p21c\n",
    "\n",
    "from py21cmfast import global_params\n",
    "from py21cmfast import plotting\n",
    "\n",
    "\n",
    "\n",
    "random_seed = 1605\n",
    "\n",
    "EoR_colour = matplotlib.colors.LinearSegmentedColormap.from_list('mycmap',\\\n",
    "             [(0, 'white'),(0.33, 'yellow'),(0.5, 'orange'),(0.68, 'red'),\\\n",
    "              (0.83333, 'black'),(0.9, 'blue'),(1, 'cyan')])\n",
    "plt.register_cmap(cmap=EoR_colour)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This result was obtained using 21cmFAST at commit 2bb4807c7ef1a41649188a3efc462084f2e3b2e0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook shows how to include the effect of the DM-baryon relative velocities, and the new EOS2021 parameters.\n",
    "\n",
    "Based on Muñoz+21 (https://arxiv.org/abs/2110.13919). See https://drive.google.com/drive/folders/1-50AO-i3arCnfHc22YWXJacs4u-xsPL6?usp=sharing for the large (1.5Gpc) AllGalaxies simulation with the same parameters.\n",
    "\n",
    "It is recommended to do the other tutorials first"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fiducial and lightcones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's fix the initial condition for this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-14T21:18:50.491836Z",
     "start_time": "2020-03-14T21:18:49.277655Z"
    }
   },
   "outputs": [],
   "source": [
    "output_dir = '/Users/julian/Dropbox/Research/EOS_21/'\n",
    "\n",
    "HII_DIM = 64\n",
    "BOX_LEN = 200 #cell size of ~3 Mpc or below for relative velocities\n",
    "\n",
    "# USE_FFTW_WISDOM make FFT faster AND use relative velocities. , 'USE_INTERPOLATION_TABLES': True or code is too slow\n",
    "user_params = {\"HII_DIM\":HII_DIM, \"BOX_LEN\": BOX_LEN, \"USE_FFTW_WISDOM\": True, 'USE_INTERPOLATION_TABLES': True, \n",
    "               \"FAST_FCOLL_TABLES\": True, \n",
    "               \"USE_RELATIVE_VELOCITIES\": True, \"POWER_SPECTRUM\": 5}\n",
    "\n",
    "#set FAST_FCOLL_TABLES to TRUE if using minihaloes, it speeds up the generation of tables by ~x30 (see Appendix of 2110.13919)\n",
    "#USE_RELATIVE_VELOCITIES is important for minihaloes. If True, POWER_SPECTRUM has to be set to 5 (CLASS) to get the transfers.\n",
    "\n",
    "\n",
    "\n",
    "initial_conditions = p21c.initial_conditions(user_params=user_params,\n",
    "                                             random_seed=random_seed, \n",
    "                                             direc=output_dir\n",
    "                                             #, regenerate=True\n",
    "                                            )\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATkAAAEKCAYAAABpDyLyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO29eZhlV3Uf+lt3HurW1NVdXT1oHkACWzIKTozjgEeMnchDbMOXB9jGlnkxH/aL/YLwBHFMPuwXsJ0Xx4lseGADAmyQIX4eUDBEJn4MEgghEEJTS93q6qquuerO99z1/rinag11b3V1zXW1f/3V1+eefc4+e+9z7r5nrf1bv0XMjICAgIB+RWK/GxAQEBCwmwiTXEBAQF8jTHIBAQF9jTDJBQQE9DXCJBcQENDXCJNcQEBAX2PXJjkiOk1EnyKiR4joq0T0C/H+USK6l4gei/8fUee8mYgeJ6JHiej7dqttAQEBzx3QbvHkiGgCwAQzf5GISgAeAPBDAH4SwBwzv52I7gQwwsxvIqKbANwN4MUATgD4HwBuYOZoVxoYEBDwnMCuvckx8yQzfzHeXgbwCICTAG4H8N74sPeiM/Eh3v9BZq4z81MAHkdnwgsICAjYMlJ7cREiugrArQA+B2CcmSeBzkRIRMfiw04C+Kw67Vy8z9d1B4A7AIDSmRdlRzqnJ5vuuIi7bgMAp2htu52UbZA5DO2kbPv6WR3L6jhfB9SlqY3e0M1ouSLu3Zd2KtH1OH9tvR3lbCNZ/dSRe29up1VZQQqTCduZdEKVEfcsS6kLFKhhjmtBBrLF9vc3owalztKo6WrJHMdNdZ4bb93rhLqfiYZtb7Kib4AbU/I3OD6qYftCmayUpW1f9D3T4wu48dfPhLvv+llqp2ybzD3rUR8ARFn1wb3u6PPqz56bYeaj2AbaF25YZzImjn+j+2DuMHZ9kiOiAQAfAfCLzLxEPR4SrJ8egHVPGMDMdwG4CwDy46f5un/1bwAAxUn7RGcX5EFNL9kZqn5E7m6jpB4497DURuXzwHn77Y8yUtYYlG1O2jr0w5kuu770mChzC7YvyZp89n2pHZW+JBv2PP2gUkvKFq7PmONaBWlIZsEOeeWElCW/ZWFteyhfM8edGFhc2y6l6qZsIidlY+nlte1bcs+Y42ajga7bAHAyPbe2faYh37ff//J3muNa03lpb9XeCzXXIj8lZYNP21+V0kNT8qFtx5QzagZRPyrtM2ftta66cm27ecxOxPqelcft7JJdlDqjrLQxN2+fP/1M1I7Yr/HKSakzO6fqcz9uK1eqsoLtZ2pZ6njizl96GocYu7q6SkRpdCa49zPzR+PdU7G/btVvNx3vPwfgtDr9FIDzu9m+gICAvUG7y7+9wq69yVHnle1dAB5h5neqoo8DeC2At8f/f0zt/wARvROdhYfrAXx+w2swkKp2fo2SzuTQb1DlUzlTpt/e6iPy69a0P7hIqBeSyrGkKYtUlZklZU7WbTuyS3IzU1V7Y6NMomtZqmrfLFi9/XrTJ9FSb4qL9i0vPS1vTdFwYW17YNL2ZeW4fG4O2l/7VkHqr52TAeKTvS2NfNq2Y7Epg/XCQXkjmY7sgOdIzss4m/3+8jVr2xMZeaO86cQFc9zD7RNr21HL9aWuzOG8bEc5+zXIXRxe205WnJ9CvRFHg/JGRmM3mcPKY9JntsON2ojcw9I5+4bWKijLQj0unLB9WbpS3ihrR2xZblY9E1XZ9vc2vSyfs7O2kd6M3i6aXdYP98RXtsvXeQmAVwP4ChE9GO/7FXQmtw8T0esAPAPgxwCAmb9KRB8G8DUALQA/H1ZWAwL6A3v55uaxa5McM38G3f1sAPBdPc55G4C37VabAgIC9gfRPkq67dUbY0BAwHMY7fVriHuGQz3JURtIVTvbjnGA5qB0rTFgC1t5ecGM1EJjO21vhC2zL6VRXq1MqZXWzKI5DMUptTK6Yv1MhZmKfGiI76f9lF11TAwNrm3TQNGUZbKqkXMLpoxyynGofHKFcxVzXDMvK5krJdvPZM0QL9a2yhdsOyoFudbY0SVT1lbLyA8sXLG2rX1rADCYqK5tlxJ29Vavyi5GsoK62LD+1sGS9G2lkjVlUVKthqblvqzknI/yGal/+Mt2rJafN7q2rX1r+TnrWdEr262c96PKdm66asqag3I/M2qlPnJ16OexOOmeW9Xt6pHelJqk8jmn/Mr/Di9JRmGSCwgI6GeEN7mAgIC+RjP45LYGYqGOtBzRsawIkZG1WsySvjZRrWlmX9lbA/YmJRTBndUoJh2FJDMvB6bPz5uy6JzQABPKDE04k5TSsp7Pxbwp04RUJK3ZxXnpeKskZlCqbGkRhWn53MpZonDluKLYKH5uouHMJzV0M0+OmrIZVchFsdVSLrziVEFRQwqWIvm8rHxeiGR8nj88ZY77XE3M4VbDPt5cls+ZUTGHOWfdCAvXS/3t1BFT1lDmfFYRp32kSPWI3At2y286IKR2zN1Phfyk2JCNEWuWZ5ak/kbJu1LUtVTXamPuGY4U0dtazaZvO4FgrgYEBPQ1ov2b48IkFxAQsPvYP5ZcmOQCAgL2AFFPyuzuo28mudoR5yNS7qnaceswYeUUyU3JECRtXDkag+q4GXuTRr4hzo7srPi0dCA8ALBSnEDkAu9PSQgSWqqNGRtTw4tCyaCqbWQ0Ko4yylj/UVQQ/5oWDli8rmCOM+4S9yxqtQ4T8B85qklFBnzwKVtHuiwXaJSkb1+Zvc4c99Ap8ZN95dgJU/by419b255Ii29zOGUpHsWMNLics77HxpJcuzkljisecPSPjLR3/vm9KTWZZTlO++AAqxLSsi5WE6aXKtrnNl1RYWMDcv9WTlpfafWYaseStQVJDYmmq2DW9kWf50PP2js8MzS9Y3KTIKIzAJYBRABazHwbEY0C+BCAqwCcAfDjzDzfq44gfx4QELDriEDr/i4DL2PmW5j5tvjznQA+yczXA/hk/LknwiQXEBCw62gzrfvbBnoJ73bFoTZXoxRQPt6Zp9dFPAzJq7jWxgKAZEVpwQ2LecAJe9zo16SOgWfsGnvmqem17faChDm0br3eHKcFOqPxYfRCYknqbw9ZczLZVDZH1UYCJDSV4KStX5s41aO9f8807SWxTrBTtjOLysxylINkTQ4sXnDmX1vRdBT1pFm0NlK9JTSJM0vjpuyeppiar7/mvrXtK7Kz5rhsSnTcshmn5nJUGt1YUeafE/lsKjdFeslp0inq0OI10hcf6TKkNOoaA/ZrptVuilP22jo6Z/mkjIfWlgOA3Ix6vmu2juyijH9tRMY4VXHXUqokPqJnpxkf3d7ctABujLtivUjfkk8QEQP4b3F5L+HdrjjUk1xAQMDhQNTFaNQCuBvgJcx8Pp7I7iWir1/utYO5GhAQsOvYqrnKzOfj/6cB3INO3pdewrtdcbjf5BJAK14kqx530t/qY7JuB7R+VF7nE43eJphnsWvwoFoyW5Lg8fSsjXQuX7eWcRFJF2mgoxASud63IpmWsuiCvZ90XFZUW87807Lm2szSq50eWoQTsIx5LR+eWe6dQKHpVgz1vdCy9JzwQqTKbK7bOqamhta2v3pcUn9knX1dVWZttepk3hs61IW6bwNglcsiMW9XupPKW6Cfj1TVjpsXuTR16HQQPq+IMhu1XLmPpNEisT53SCvf/d3Fy/tbwVhn8s5uy2e2Dg2/fLsJEFERQIKZl+Pt7wXwm+gtvNsVh3uSCwgIOBRob81oHAdwT5wXJgXgA8z8N0T0BXQR3u2FMMkFBATsOrZCBmbmJwF8c5f9s+ghvNsNYZILCAjYdUQ7LVB3GdjNRDbvBvCDAKaZ+QXxvg8BuDE+ZBjAAjPfEudlfQTAo3HZZ5n59Ze6BhPQjt0uXkGkNax8PyPOabEi3U6t9B78paukjJNWBSI7KL4aukL8RZlFm39TUwIyS3a4WwXxU1SuFP+RF2BMX5TzEkMu205ZnETtlE3jp31GSUXjSDv/UVOJiOrEO4Blwtu+WF9YQ4mUerZ8ZkXqzF5YWdseOGIpL40haYdXfYFKQjOgQlOmGoPmMH1WLm/vRZVFlaU9r6JB/EuGopREGe+/lIO1L2xg0t6zVFk+l6LeNJFE05cpP6py4Xp/sR5Tfy+ibPeEPdVx21Hth/PpG3U0x06g3adhXe8B8J8B/MnqDmb+idVtInoHAM0ueoKZb9nF9gQEBOwTGrx/RuNuJrK5L35DW4c4XeGPA/jObuUBAQH9hS0uPOwI9mt6/acAppj5MbXvaiL6EoAlAL/GzH9/qUoSEZCNE6s3raWG/BF5v48il+MhJa/67RUxQ5MLvV+p6y5npc5t6WkXvo29kFkSe0SbhS7lKBBJJewiHlAXk4xuHLPXVuaOFlZ0KWoNHcEz6w3tQJ2Xe8pGGqRHhFJTnXB5KJZVPzPyyA2ct2IDnBRzsuoEF8qihWlM1CeWbJ+jtpw3WrA23kX1HFQK4m5ILtqvQVsxTzzzweT3VZvkTNLMrETJc8qOh6aXJPzNUMO9EdXH5JDIu0aq0yrj0ufaqK2vrSg7OpoFAHJzOyuOFG0vjGtb2K9J7lUA7lafJwFcwcyzRPQiAH9BRDcz85I/UYeCpEsjvjggIOAAolvEw15hz69MRCkAP4KOVAoAgJnr8bIwmPkBAE8AuKHb+cx8FzPfxsy3pfLFbocEBAQcMLQ5se5vr7Afb3LfDeDrzHxudQcRHQUwx8wREV0D4HoAT+5D2wICAnYB+/kmt5sUkrsBvBTAGBGdA/AWZn4XgFfCmqoA8B0AfpOIWuiI472emecudQ0mSVITOeHDek35XJKOFrGok7ooJQbv2tCfXejP0lVSvxYfzC7Zm1kbVklcyGXUUVVqOkKmZvuihTFTbRc+NC8L1Imm7af21eiwo7aNdjINWReepNqoI6jKzz9qjksrGoPPJcpJGZOWEoLMTC2b44rquPqgy5mq7u+XZ0VQc3rOUkgS6l4XcpZC0mx0f9zbR+1xOjcQt1xYlwoRNBFlzuUUFaWftRF73cyS9CXZsPdM+8K0f9SHdZn2uhCyxqAab8V8ap+0/tz2svQtO+9CvoZ3dlJqbiGsa6ewm6urr+qx/ye77PsIgI/sVlsCAgL2F31JBg4ICAhYRb+SgXcfKuKBXB7QRKL3631mRplxSjTEm2qeTqGhFRx0ztf6kD1Hm6Hr6BnqDV5HIfg8EVp0UpuuAJBQJl7WmYmpcTFHNJO+MGPr1/SY8rg1KwbPqjypFaXOUXfcGJVDIipaE08z8DdClFW5cnO9j5uc6i0+qu/gwowbq6y0OTcipluj7r4GKn8FNV0kTQ9lEC8A2hhQOSRc9+tKyDLRsoVZZcrqXBDefNRimOmqvZ/6XjeGVR4H53JJK8HViouG8PlOtovwJhcQENDX6MuFh4CAgIBVbDOnw7ZwqCe5dlLSBnLJhgmk0vLaXy/b5cS8WmSK1CJeO+NMk4JakXQjpdMCtLWplrd1aHO4OG3bmJ2XSIBWTpkfS3a1Lzkrq5Cto3Y1kSrSmQRbcztZl2D+XLN3hEaqprdtWUKt/qWWpV1a8BOwkRGpsgvZyCihA7USWLnKmp2tog5ct1VQS8qSKo1kdMyO1dXjEonx+OPHbRtVDg8dS5nJObEBZb62xl2Qf10emPSyEp0ccKZguXtkBACzEutdJBktxFmTdlXH7Gqzzs8Q5exbUkOv6B9XLow5W4eJmij6yIudnZSa/Ri7GhAQELCKkFw6ICCgr7GXEQ4eYZILCAjYdYQ3uS2Csm0kru44vXKOMtJqio+LW85PpnxQOl9ow9E/jB/O/RDpq6WV381pOBo0C7aSwlnx96QnpVGcsRQMNMRBlahYZ1XjGslPWh13iVsU3UFHOXj/kU7Okp+1/qnMvHAJ9LVzMzZaYeVmSX1ZK1rfT1JTHNQQ+PHQ/jqvglFZ0n49VV3aHqcd3CeutEopM4tCKWkoP226ZCMB0mkZg/KK5bLUx5UYa0pRkVbsmGrfbMK69ZCdVxEyi7b9SUXNaZbkASyfcPlf1WOQcTIWDZ2UZkbuReGCo1mpOuoj9vtTO+mcottEeJMLCAjoa/RlWFdAQEDAKgIZeIvIpFq4YmweAHB21tIR2hfEzMgtuDygOveBYnbnp+0re35KtluOGlI9zqpM9nvqgzZVaiO2HYNp9eumTJ/Ess3dygXpS1SypmBlQj6vTNj666Oqjg0iQHIqx2Y7bR+J2qjIWQ0+KXQEumj1E9LLYsatnLJmszaHdV7RdNkJM2rWhcsROvC0jE/tqKJPpGzkhTZXSxlL2x+fWFbHyViVW7a9Qxnp57mMfa7mkjIeiUnZ9jkYqsdYldm+ZOdlm9q970uU6e5uAICBs8rkXbZj0FA5QRoq9wltQGVpDtt7MXZChB+e6dnCzSPw5AICAvoaIeIhICCgrxHe5AICAvoaz8VENjuCZpTE+YUOZ6O+ZH1VhYsyqMVJl8Cjx0KPX87XPozEjKtDJWSpTPT2qxQuqJymC7b+2pj42jLKP1e70eau0Hk6I6e7qZVNorwvU37DI+IzSxat43B5RCrNXbCDU3patqvHpb3Z7ClznPYtDZyznAmdX7ahQrcGp6wjK9GQNjaHLXWjOK19S1LHymNWAv/MMelLdthSQ7IZqf/moxfWtr1PbrIsPKDlqh3wYkH8fAvHVRnZr1KkksQknJKJpioVLjraUlrltlW+tuKzzq+sQviitK1Di7jqcLv0in1OV07LdvKI9V+eHhTH4RexfTTbYZILCAjoYwSeXEBAQF8jRDxsEe0ogcpibKM5QUCtv59dciKRJflV0UvzKS8+ONjdzAKA8imVu1XlCEhOOXtSWQj5WbvUb3IybPBDpwUkyxOOWa/MolbJ5bJIKyFOJSrarts2ZqcVc7+CntCUhsagy1uwKAOeXrBmYktlVdMugWTFmrWJWaHuR3mbQ6L4jJhT7ZRELlSP2oHLP6OiRZ62kSPL4zL+n6/IGIwOWcpOOinHtZyoJSsHenpY2lR3oq2aNuJpRfqZ8/lac2eFusFJuXa+OGSOqw/1VmzROXY1mm5/47icWCpYc/Xh8xNd69gq9nPhYdfeIYno3UQ0TUQPq31vJaJniejB+O8VquzNRPQ4ET1KRN+3W+0KCAjYe+xnSsLdvNJ7ALy8y/7fZeZb4r+/AgAiugmdLF43x+f8FyLavziQgICAHUUbtO5vr7Brkxwz3wfgkmkFY9wO4INxkumnADwO4MW71baAgIC9RbOdXPe3V9gPn9wbiOg1AO4H8EvMPA/gJIDPqmPOxfvWgYjuAHAHACTHhpCIpVTZqf9qP9bCNXZAtQqJkRNxTBC9hD97s0vOopWIm1ohw1XCOmTKJxKROhINuRYN22tplVmfHzNSSiPZOdvPZgldkbQuM4w8phLUtHrTYVpZ7QfylBpVNmtlMQpN5Yt8VqgblLM0kdYFiaNLlwq2/gE5Vucqzbg8t9oKys3aNpL6YlUy8rwsp62v9NqjM2vbS1XbxqZSt0mmpB1Np4YSRdIQn27X0IDcC432w3FOvp4+/EsrCvuQQ00lWlXOBoCNxHlXFux401y6x5FbQ1/65HrgDwFcC+AWAJMA3hHv7zYCXb9tzHwXM9/GzLclS8VuhwQEBBwwbNVcJaIkEX2JiP4y/jxKRPcS0WPx/yOXqmNPJzlmnmLmiJnbAP4IYpKeA6CoiTgF4Pxeti0gIGD30GZa97dJ/AKAR9TnOwF8kpmvB/DJ+POG2FNzlYgmmHky/vjDAFZXXj8O4ANE9E4AJwBcD+Dzl65PEtY00vbFr6mEAyOn4NAqqNyZisWgKSMAkIi0SeCurSkDynRtb/CW783VZkkObg70TjSqc7fmZ3srdzSLvdUuNKUmP2fNs4GHxUxkp0LCBWVbKZZ9bcLmNK0PKTPuGkv/SJ9XtIimNISGXWKVtLpRTuWkPSq/gaT6Upi241FWSizJhn0mhp7Sn6Wf1ZY11VaGpF0jBRuVcfbcEXQDeQpTo7eCSDQgbZ75JvvA5E7Ii4lWyyk9a3kims5TGbfvKpUrXE7cGJx0xpHKL8uRffbp2M4mXt3KaioRnQLwAwDeBuDfxLtvB/DSePu9AD4N4E0b1bNrkxwR3R03ZoyIzgF4C4CXEtEt6Hw1zwD4OQBg5q8S0YcBfA1AC8DPM3P3OxUQEHDo0OoyyWn/eoy7mPku9fn3APxbANq7PL76osTMk0R0DJfArk1yzPyqLrvftcHxb0Nnxg4ICOgzdDNP4wntrvVHA0T0gwCmmfkBInrpdq59qCMemIFWo/OanXBB51FdmS1VH9ws21o00//YRCqfqhdFbNZVYLwWv8xa86ml8hj4VXMdjK3bVD1iG1I7ItcaOGvrzyjhyVTNPkg6iDtRl+3cBRvWwLPKrnX5JRJVsdNbz06ubecb15jjIiUqUB1zOVmTUpZTgqCc9+aqWk1MuUdTrS6mVV7XtA1WQKoq5+mcEZ0yGauSUoLUuRoA4JnSKHqBVqT+nBKB8O3QYqkNG6yAplp51W4VAEj2iJRoDtg2Vo6qleKb7HL58Ig0Zn5aXoKSi3ZMSeU+SVxlO1DIucQU28QWVldfAuBfxAEDOQCDRPQ+AFOrbi8imgAwfamK9i9qNiAg4DmDy114YOY3M/MpZr4KnUCBv2Pm/w0d//1r48NeC+Bjl7r2oX6TCwgIOBzYQZ7c2wF8mIheh44y+49d6oQwyQUEBOw6thPGxcyfRmcVFcw8C+C7Luf8wz3JtQnt2Cenl8MBgFJKODBj/R6ZRTk2s6ASjtTscZWjvf1pWggxfV6GsVWwdWg/Hyf9jdbCiqodFVtH/XrZ1n48wPobB561/rqB8+JX0f6pdt7e9mRJ6CDcsnlXuaEUVgbVcVnruyueWZGylPMpHhW/XvWbhAqSXrA0heTpE/KhYv1MtXHx32kKSWa5d35QnXcWADihx0r6VT1m6TvVshof52PV6iL5i3KfilOWDJC7KPXP32D5RzmlxOIjGSpHlCKM8qN6FZyaWlNMZey1kzrqRp3Wdt92Lsp5BZcQqN7c2amhFUQzAwIC+hkhx0NAQEBfI0xyWwVjLTg+PZ9cXxYjP20HePAZeTVPKlpBfcTWUT3am0LSVhQS425wb+W5OZ1/0+WQUBEV6YoK1nemQlYFSzeGXWC80vD3ApLUFhNPiw3kli09oD2/IOc46kb0vCvXtllRZZKL1pykc0IvoZPHbRm7cJHV62bseLePCd0hUbfnVNW90TkSck6UQOdybbvcrQllkem8Ez6QvzmgoiFOWvO9NaLyUFyU+xLNO1HVk1pQwOXzvSBjl3r0rCkrXC0me2NU6pi/ztJySNNLlmzZXEu1RUfm+IgHJTBQr1v3w05PSRwmuYCAgH7GXurHeYRJLiAgYNcRzNWAgIC+RhRWV7cGahFyFzpdSLgoFO1nKV6wfpXCORVapPxMKyesPp1W/0i6+ls9mAvpJfuLVZgWR1B6xSWyUeKPqSXx0yRX7AMx+nVp1/lvd363lt62164PK5UJJcaYc4Ew2g9Ho8OmrDohvrHSV9SJdTsgrKgQVLXUEB3apn/QyYUqNZU6TJSxtI6VU1LWKio1FCt4guJZnbfU+qCMIKhJMGSfj2ZBxqMxZMe7nZcTKxNy/zjZ+0tcesb5YluKQlKwCijJi6LYkiiJr8379TSVJTdp/Wnat5yfkXalXOjZyjdJv6MVp4ZyfmdFM4NPLiAgoK8RzNWAgIC+BvOlj9ktHOpJLtFUppf7ocguqrwIc9a2TFTlc3lCTMGEU7DTn71CCStLSyuZaBY8ABTPiuJHVHC5G2YkSoDqqo1ta95kB8RsobY14yjSNBd7bS0OWjkm26UzThxUKY+wUyFpDEjH61eIOkfm2QVzHArSrvopa/KuTMj16kqsOul0GZtKh9PnytA5Zds5NT6eFaFUQsjpi6ZVJElaRUoky/b5aKflmagdsV+RuqJd4IS4GKpk70t6UcYtO2/N4eSCUoFp2mvXbj4lRSU9bk4QVVGTBs7bQZi7SR2rirzwa1JFSiTOWhrKwNM7OyuF1dWAgIC+Rlh4CAgI6GsEc3WLoAjIzXfMB3KDqJnvPgi6nVcRBCX5hamOuRwMarFVp3YDgIxaRc3OS1lmxa3o1cVU8YHrUUlMnPYRuVjmvDUFG4PS3vSybaMR/fRBH+qzzhdQOWFNq8FZZUPOLZqy/KyYqPURZdYmbZKkzIyYYLRBWsOmGsdazgW/q1Xl5qj3HajV24y6t0lbR31UOp1wog3pis7tIX3R6SABIKmEGvxzpU3lgbysMK8MWjM/O6vTVDpTTYmUcmSvraNgtKsg7VaK9b314hF6hX/lOmsqG1SUAKgbK2/abhf7ubq6a++QRPRuIpomoofVvv+LiL5ORA8R0T1ENBzvv4qIqkT0YPz3X3erXQEBAXsPZlr3t1foOcnF+Q0v9Tfc63wA7wHwcrfvXgAvYOZvAvANAG9WZU8w8y3x3+u32qGAgICDh22kJNw2NjJXz8d/G7UmCeCKbgXMfB8RXeX2fUJ9/CyAf7mpVgYEBBxqHFSf3CPMfOtGJxPRl7Zx7Z8G8CH1+eq4viUAv8bMf9/jmmtpzDKF4TUWe/6io4k0FRvd+USoKX4Qo1ThpnMdReEFL7WPJK38cIUpy4vgjAyx94Xp3Jmlp4WO0HY5WLXfMOnUUFjdwbQT24yyil6iysi5u6pXyQt5oWbbnzu/vLbdGpR2NUYs5aByk2RraeXtWEWqO7mLKgojYZ1JmqbTzrpIAyV8ar4weecD1XqRLkpFU0r0M+GjFUxEgqOhaCWPhMmIZNuh1WK0wgwA0ID4X3nFhiEkaypJT0X8fD7ixgtgakQqP5BO8JRM2c40F+XAxqAtS9R31pPVPqCrq/9kE+dv5ph1IKJfRSe/6vvjXZMArmDmWSJ6EYC/IKKbmXnJn6vTmA2Mnt7H34eAgIDNYj+/qD2nV2auAQAR/WMiWhP6IqISEX2rPuZyQESvBfCDAP4Vc+c3mZnrsXY7mPkBAE8AuOFy6w4ICDiY2M+Fh81QSP4QwLeoz+Uu+zYFIno5gDcB+GfMXFH7jwKYY+aIiK4BcD2AJ7VA59kAACAASURBVDdT56rZkZlc99K3hvagXQ+PBuQ1vTCtoh9OunyhWWUiuZ8DLaKpc0NowUUAqA8r2kXvFA+oqVylqbw147TJ5CkNLWWa1AftBfKzcqLOIZGsWdMkqYQCuGBN5cS8RGVgRMaxWbADsnylouIccxEbsyqnhli/SNadqdnUZbb+6jF1bEXlq51z+Srq2hw2RaiNyI6EsiCzzp2hA/nTy6YI+WcUVUSJAxSHrR+hUpZ7GLm8HJxTN61sc+BqColGfci2sTgpPofGoK2/OaTym6g8tKkh64rIjMg7SoPtfdf5MHYEB9QntwpafeMCAGZuE9ElzyOiuwG8FMAYEZ0D8BZ0VlOzAO4lIgD4bLyS+h0AfpOIWgAiAK9n5rnL7UxAQMDBxEFXIXmSiN6IztsbAPxrbOIti5lf1WX3u3oc+xEAH9lEWwICAg4h2u2DTQZ+PYBvA/Bs/PetiFc3AwICAjYFpvV/e4RLvskx8zSAV+5BWy4bnASaxc5geb+bVvLwuUqNYOKTYhUXj42bw+ZuVh9cFToEjEmF8Dg6glYCIRdho1VD9HnrwoBUe9tOyzCpln6qx+x5WolFJ7JpOzqMvl5rxIo4tk4LNaQ2Io9L+YTtZ31UNdI1vz4m96Kp/IZeYFT7yVq2GYgKUkduSvxdzZIPo5NtP1aaWkHqzaJw3h6XmRJHXOrqrCmrj8p5y4vyzBVK1t+l1WEaQ/ZrlntaUVQGS6YssSg+uqTOB+scjDUVvtay7jS004rmUpBBnRixfuvluviB5+YdbWnEc2e2h/3kyV3yTY6IriGi/05EF+MwrY/FiwMBAQEBmwN3+dsjbMZc/QCADwOYAHACwJ8BuHs3GxUQENBfOOgUEmLmP1Wf30dEb9itBl0OmMQEqZy05mp2XjG9VyxdXC/TU0Xsvdy8U6OoK63/MVvWUvkJasrcyyzanyhtIkXWMjFRFG1l7SRc/ohVk7xb/TVlNtfHbRsXI0VjUGKYnsWvcxxEWUujqQ+pfioKRmXciVqOSqOp5igwze4masqxLOsqSkDncQAAVuofjWGVt6Divizqox9HjcpxTWuxptrIgsqLWrPt0Hk0uC79rKwMmOM0laV6xFF7xgfXttPTjqOiFErSOg/tUftVLZ9Q+TCsRY32oHR8ZFDM32zS+kvOl6Ud68RH/bhuFwecQvIpIroTwAfRaepPAPh/iWgUAALVIyAg4FLgfVxd3cwk9xPx/z/n9v80OpNe8M8FBARcAgd4kmPmq/eiIdtF5agzkdSbeWHB2kWJJcUyV2Zcds6vkPUenkhFQ2hzzEMLavrAdc341yuhWvATsOKPc8+3/awdFROVhqxZHqkVs8qEWuV1Afo6SsObstqsqx9TOQFG7FjlMzLg9ZoTkPy6uBKySg/Up9lr5dSK5FEnJlmVfrezKkKj5oRItSiBEzPQq6064F9fFwDK14oZl2zaNppIFxVtEZVcFElNR1644H0lbpBasbYmJY+sbSeUUIA3Sesjqv3D1gwdHltBN0w7k7q5Iu1IL9jnKrP4HDBXiehHNjqRmT+6880JCAjoSxzESQ7AnwN4MP4D7PsmAwiTXEBAwOZwQMO6fhQdf9w3AfgYgLuZ+fE9aVVAQEBfYStkYCLKAbgPnXj3FIA/Z+a3xIueHwJwFYAzAH6cmed71dNzkmPmewDcQ0RFALcDeAcRHQHwq8z8Py+/ybuEePDSZZesRglSRgVLi+C0SnZydnptO33W+oFyMyJuGOV7Uwp1chZO2V+sjMoL41UxSP26pSty7dyM9SFWJsSn1XIikVwSf0yxaP1kVaV2oX06xfNe1FKx+K9zYzAuoo4ZtULWqFq/W105vOiCdSAVJhUDX1XfHHBCpMqV1HQ+Ih3xoPvc9FoRqo2JZm/hTe0zi+zjgWZR50y145FdkDLtY40aXglE2ltz7cguKRWVk0VTlpuWe091rTRixyql3MrRcdvGQlZ8s+mEtGOuYmlW6YtKlLNm6/dJkbaNra2u1gF8JzOvEFEawGeI6K8B/AiATzLz22Pmx53oqBt1xWbIwDUAi+go9hYB5DY+PCAgIMCCeP3fpcAdrP70peM/Ruel673x/vcC+KGN6tkokc3LiOguAA8AeBmA32fmW5n5by/dvICAgACFLmFdRHQHEd2v/tYJfxBRkogeBDAN4F5m/hyAcWaeBID4/2MbXXojn9wnATwE4DPo2MSvIaLXrLWZ+Y2X1cndAIv540UFtYAkufX3wjmrq7+K9ohdYh86I8zxRMsOVWNY6q+p4PQoZ3+iquq91wekaxO7Mib2QfWINWG0GEBryJomCaXbn0pYGkOkqAWJFUV3yNh26OiC1Kij2yj6Q3VOzJ3kcm97Jj9lfzs1VUSbXT5PQVuZjQlHc2krGg23lJiBy91KFWmXFgbo1C+f00tyXGPI3xdllg86yo4K0NdvIz4fbuOIXKtVsM+EjmBJOjM3eVH5N+pidqYqQ+Y4LYjAVTuQKzV53vMZeYZ9noVIuT5KZ0yREX7YEXRZeNCpDHqexhwBuCXODHgPEb3gci+90SS3SvYNCAgI2B62OZMw8wIRfRqdNKdTRDTBzJNENIHOW15PbLTw8J7tNSsgICAgxhaUm+K0CM14gssD+G4Avw3g4wBeC+Dt8f8f26iejXxyb91EIy55TEBAQMAWRTMn0ImdfwjAF9Dxyf0lOpPb9xDRYwC+J/7cExuZqz9DRL2zw3TIwa8E8NauhUTvRicr1zQzvyDe15PfQkRvBvA6dHI8vHEzCxypWhsj3+g4D9opRxdQVI7stI3vSTTEV0UZcQTRkk0qkptUKh4566/TiT4ayl3i/Uwmd6v3QSkVEkqpZDguzEj7X2jQhm7pXJqVmuVCkPJB6V9SH7pllDXOWrXKpmKlFFd0Hldbh3Z7Dj5tHWqlp8QHunyNjGOj6GgL6hZqsVHAKorQgvIvjlipEc6ra7uQL51LVIflpZyIqPYNtlzMZUuxMJLq0pWTzjfY6J1QR/tYo4z1+TX/6cm17bH/Nbm27SlSTeWb1fcPAJZmxadbzslx0Yp9AJP6NjlzsjDtFF63ic2spnow80MA1uV+jjP7fddm69mIQvJHAEob/A3Ex/TCe9CxnzXuRIffcj06Cxt3AgAR3YTOhHlzfM5/IaKdZuoEBATsF/ZRNHMjn9y/207FzHwfEV3ldt+OTgYvoMNv+TQ6JL7bAXyQmesAniKixwG8GMD/t502BAQEBGxGamknYfgtRLTKbzkJ4LPquHPxvnWIuTR3AEAmP4zakY5tkZu1Zpw2E/37amtIeB3tY/Jqnyxb0ye5LOvo+Qt2TZ3aYp9pFQ8d/dBph6ILuGV5rUKiTUivzqHZ5+2avWWsTDDOOjtU0T9y88pUcwx/HWmQXbCmT1OxWfLTysRzOVN1voriGauCkVwUczU7LOOWrNuX9bTKQxHlXW4FRfPQ9BKasp2pqfys0TH7TLSVkkmyrCIXBly0zIxs+3QbWcXw0G8jmTmXW1W5H7y5mlYMJj1ugL03KzcJ/ctHh6RWlCvF9RNKtFSncSUXeZFelM/1UVuFp0xtF1sxV3cKO5xBdsvo5oXsOizMfBcz38bMt6WzA90OCQgIOGho0/q/PcJeT3JTMa8Fjt9yDsBpddwpAC6HUkBAwKHFQfTJrYKIfgfAbwGoAvgbAN8M4BeZ+X1buF4vfsvHAXyAiN6JTrKc6wF8/pK1EdCO38yjrFunaMsolk+45AqujlUUplyRuhG8QVrD0UdlJaox4MyWpDJJHYu/NtrbXNBIKOHG9Iy9ZVqLv+murQPS/aqvaaOy0uuDtkyb2GllWqVXnLBnU6XZi2xHOadEIisyVq2iC8LPSvuTzozLmXSFauV83h6XVO6BStuasm1tQqrxaDtbaulaqWPgGSdm4MQrV1E879wUqr2NkjM1Vd/8SnduSXaUx1XaQZeiUT9LlLKVsIpoIeWy0Ca6v7bPh+FXt7eLg26ufi8zL6FDBzkH4AYA/+elTiKiu9FZOLiRiM4R0evQg9/CzF9FJyPY19CZSH8+DucICAjoBxzkNzl0Iv8B4BXoaMrNEV16lmfmV/Uo6spvYea3AXjbJtoTEBBw2HBAlYFX8d+J6OvomKv/Og612Onw3YCAgD7Gfpqrm0lkcycR/TaAJWaOiKiMDq9t30EtRm6u4/xI1K112xgWf4xPIJNoKRFH5TuJMs6fNih1NIYcW7wmDo3MvCzhc8I6bTKLcoHyCecjUrqTTUVjiAaco0Z9zM5aP1ZK0RGyc7afK2opR1MaCs+6N3HVbe+7G3xKLl6Yln5mpixNhJYlBKJ647gtU+2PcnKxZsH5EBVNJHL3TPukdESQ9x3pBEbFs73FQXU/GyP2GxgpX2btqCkCqVVBUn4s79NKV5Tiicthaig87suvfbpa6abpaC46r2srsvUns/JdaCvaiJ9orGLOLq92HsSUhET0ncz8dzqhjTNTQ46HgICATeGgvsn9MwB/B+CfdykLiWwCAgI2j4M4yTHzW+L/f2rvmnN54BShdqRj8xXP2oBibZr4ZfpUTS2rK+Z+qmJNXk1paDuGSv2oDF1LmV1eXGH2ZomuaDh6hjaTNN1D004AgJQ5Uh+1bdTMd/8gaVM8WZU6cgt2QKpj0v6Ei8vOqCiERE222QkioOQ4DgrlE2KXa8FIL1apx9iLGWgqS25WUVkqvaNDPGVHUznKp2R/q+AjRWSzPm4r0WKhpSflwIyj1GiBCGq7NipxB1+mzWjd53y9t8nripBWOXAjda1W0Q5qQgkY+DHgpZ2l0B5oCgkR/SkRDanPVxLRJ3e3WQEBAX2FfaSQbGa6/gyAzxHRK4joZwHcC+D3drdZAQEB/QRqr//bK2xmdfW/EdFXAXwKwAyAW5n5wq63LCAgIGAHsJmwrlcD+HUAr0En0fRfEdFPMfOXd7txlwIngFbse6uPWeqGDoXKLlq/igmJ0TIN7P074uyoHvEhU9q3JH6aptMM8OKStg6tVKH8hC5JjA7dSi/bOnJzcl59uPcyfVKJX3qVk9I5GZDqiA9LU+ofLZX7NG3b2Fa5bZsl+1iVJ5TPUlWfdvmE6sOy7cctf1GHQimBUafOoek2XimlckJTMuRArRTjwS3bT03d0P7cdtrWUTnW20+rqUMJR//Q46PPy884353yZyZdntua8rGmVDwcjVi1Ep6R83SSHwBobTAmW8JBXHhQ+FEA387M0wDuJqJ70NGCu2VXWxYQENA3OKgUEgAAM/+Q+/x5Inrx7jUpICCg73CQJzkiyqGTe+FmACqLKH56txq1WbSTQG2E4m1nVihz1S/vm2V7bR44WoSOcsisOFqHqkMLOqbK9jVfa/OTM00yS5q6oaks5jBjJjZdXoTqUU0v6f0kpRQdITtneSKpFTFj8lN2HBevlaQGtWGxxQsXXR1V+exNN63coWkRtSO2jY1hZUJesPdC91ObqJ5uo028dtZHQ6hnQglGJpzmZPW4tENTbwAgo0RFc8oN0nLXSi+rXLMbuBHIUXb0kfrtx1OYIpVrwj9z7bSiNzWVWTtgL8ZjyofRcIowUf+Yq5tZXf1TAMcBfB+A/4mO1tvyhmcEBAQEKOzn6upmJrnrmPnXAZSZ+b0AfgDAC3e3WQEBAf0E4vV/e4XNLDyshh4vENELAFxAJ6XgvqOdASonO6PlTYKEfk13qelKZ1UAs2KmV4/Y4YjUKljSB2CX1U+RNivciKZUID/ItiOzJCeWnhYbNVG1F6ucFtFPn8JOr1CSW8ZLKGtk8IyYKtnJ3pkmOWM7kFkWW3NlQq8iu1wTCRksnbavc6z0s3pKL227oPOS9LsW5UxZQqfdU6fZIHMYey9Z82acqk8NcbPkTN60EgAt2/EePCPtL56TVJeL19mIDy0c0LJdMfCpBs39VC4RHaUDAI1haZfPITHwtOxYvlr2c8Gaq4OD0v56097PWrWIHcVB9skBuIuIRgD8GjoKvgPoUEoCAgICNoeDPMkx8x/Hm/cBuGZ3mxMQENCPONCxqxpE9Je71ZCAgIA+xgGXP9fomgv1ckBENwL4kNp1DYDfADAM4GcBXIz3/woz/9Wl6lt1Q2n6gUd2xs7lC9eJPyMlbglUx10yEi2K6Bjg2TmpI7uoBBKrTs2Bum8DlnoykFRJbY5a/45OEpOfs3UkWnJefsaWmbygi9KZ1oitP7Ugg9ActQ41HemhBSRrTp1D/1yS8/0MDYm/8eSQJC4dSNfNcU/Mj61tz25EadACjO62k845mrH304h3qm7qBDcAkLsgX4vsvK1fR8+0iuLk88lqVq6UOlvDdjyoodRtnI81oRLiaP9cK+epLEr41d2K8kR3BZ5Uxh6YTUu7vE8uUdtZCslerqZ6XO4k96XtXpCZH0UcLUFESQDPArgHwE8B+F1m/o/bvUZAQMABw0E2V4noDfHCA5h5pwnA3wXgCWZ+eofrDQgIOEA46BSS4wC+QERfBPBuAH/LzDvVxFcCuFt9fgMRvQbA/QB+iZnn/QlEdAeAOwAgOTICTneasvp/N1RPOlq5mtozKmdCq+hMzbyiEjTs70HtuGyXnpA68hddHQlFZXHMei0UoI/zkRea5uKjCUz9zlTW4gM6aF4H3XuUJ2weisq4iqg4Ig1ODdvODA+KSXq0aPM/3DA4vbZ9RVbs7WVHEymlxHw9Wxg2ZdMrEm3RaElfKss2OF1TPtp5Nx6KhqKjBLKzdrwLF5TJ6PK/Vse6f2W0CCcApK6WMWi5NrKOvHBsnty8im7JK7PWXTar8rPq/L0AEOWV2IMyO/23djQv9+yYu2cPl+1zsG0c5Dc5Zv41dJI9vwvATwJ4jIj+AxFdu50LE1EGwL8A8Gfxrj8EcC06puwkgHf0aM9dzHwbM9+WHNhhLk9AQMDu4ICLZiJ+c7sQ/7UAjAD4cyL6nW1c+/sBfJGZp+JrTDFzxMxtAH8EIIgABAT0CfbTXN2MT+6NRPQAgN8B8L8AvJCZ/3cAL0JHhmmreBWUqUpEE6rshwE8vI26AwICDhAOuk9uDMCP+MUBZm4T0Q9u5aJEVADwPQB+Tu3+HSK6BZ0X2TOurCvSuSaOP6/j75mes1liWkvKp+AHVIXtNJQSRsKFf3FK+eRy1q/XrovvZ+UKJSxZt3XosB3tRwGA3Kz4tXSSkVbe1qEVSjx0SE+Usb42fe3GoLRXC4oCQGNE5Zd1IpTtrDpWsR1adUt9mFsQ10G1kTZlORUTl1WZch5btklNv3ZBHJ0tJ1aZUclZkkkZxxPjC+a4C6kh+VB14qOL0q70ovQz4+QmBi7ItZZP2q9IRdEztKLKukQwk0LTySy7dwk1pOkley/MvVE0F5/YpzyhqTK2LC8uUPN8rKQtdehManRt+7aTz5iya09dXNvekVXBLUxqRHQawJ+gsy7QBnAXM/8+EY2iQ0O7Cp254se7+e9XsZmIh9/YoOyRy2v22nkVAEfcvldvpa6AgIBDgK29ubXQWYD8IhGVADxARPeiszbwSWZ+OxHdCeBOAG/qVcnO5h0LCAgI6IKtmKvMPMnMX4y3lwE8gk5Awu3oqJMj/v+HutfQweWSgQ8UEsQYSHdMvmuvftKUzdXFfJqt2tf0qYti0uix9ioNqMjw5AZrpkjnZK2rqIPGkDd5lfnnVUgW5HP5uFyr7vKR6tyiKUdp2IhJrukm9UFFn3BCkJWj0o7FG10ugaKj069V6HI8qCiEesKaq0/MSSSDNknrjloBRfGgnL1uXUnCHBvrraKSLwgNJVG0fVmCPBNRVsbb53+tTCjXwRX2vpMS6awtKDvRRU1kpqSOpEuMqvNtIGHPaxZ16Ah6oqUeaR9Jo4MX9LWK5+zzV4lE3eZ+XGHKxod2WDKyy6Sm6WAx7mLmu7qdTkRXAbgVwOcAjDPzJNCZCIno2EaXPtSTXEBAwOFAtx/jeELrOqmZc4kGAHwEwC8y8xLR5YWcBXM1ICBg17HV1VUiSqMzwb2fmT8a755aZWPE/0/3Oh845G9yrXYCF8sdE6Togr3b6h2+4VbqMC9mBikzI+FMJJ2tsHa2ZMroqLqeWoVtW0sNCZX+zwhtAkjW5Xp65bU5YNtbOdbb1Eyr3BM+/wMnugfXt/LO5FXb66IE1GqxyafgzLMTJyWS4eTAoim7WJVohfN1cRWkizZqIpsVd0EmZV0HmZQaq6RsP/3MmDkuqVwOQyWbLOOqK2XFsJSR+/fU3Kg5jtWzM5Kzz9XUpERi6GiZ7IId+4FzamXehRo0laCmf16aeSmrHFfRCu44LT4Quftp0j4u65SS7vlQ97M641ZeFzZQ+twKtra6SugEITzCzO9URR8H8FoAb4///9hG9RzqSS4gIOCQYGurqy8B8GoAXyGiB+N9v4LO5PZhInodgGcA/NhGlYRJLiAgYNexFfIvM38GvZdfvmuz9YRJLiAgYNdh0oDuMQ71JBdFSSzETPuHK5aOkM8Jyz6fsYlhUuPiq2kuyXnteUcdV8v7qRX7g9LIy9AliuIHql1p/UycUiz7ihfvzKsy8eFk561frD6i6Cout6rNQWqbX5+QflNW+RtXrIMnodqVqLqIDZUbtnmtiGs+/9QFc9yNpam17WdrVkGkGUnDRkq9lS+emDnS9RzA+uRmyyqaoOSkXRSWV/I9PzcrMgYJJyapfXLlC1YEQvvh0uqZyCza+1K8IL48atj7Of98ab/OJwtYf1r1lDxXpeOO0qF5Ii17z1Jp6U9ajVvFfUci9RyQi9Thy1zBvCQOco6HgICAgO1iP3M8hEkuICBg9xEmuS2iDXCs6d9asCZYeUxev68esYkRtFjgVFHoDXMXhsxxmhrSKtrXeVImQkKZtZlBSzngknyuvNBRVCJF8dAihYu2L6y0+BMD1vSmKRc10AvORDXtKKl2OdJm/YSUfe8NX1/bvsYllPjryZvXti8sWLpNPittLmbFvFyqO9HMvIzVqZINvG+1xUwcyEgd2nQFgJoSB9DnANZEhabGpFxwvYp0KT5tvyI6gkDnAMk5F0PtiNxPn59BR1j4vLGkclloN0Kl7Cgd6tUol7PPxMSQRIQkFQt3JWeflbmMjF2t4lw1LqJluwhvcgEBAf2NMMkFBAT0Mw5Ttq6AgICAy0YwV7eKiJBc7vgOqGn9HtGc+B/mjli/zakB8fckSjL6lZr1WWhxxnJk608oSkZbldU8lWXYqlho5BS1ZUD5o1oj1v9XrYu/5MSIDZl6JikhSa0F61fRSh5J5cvLOh+ODqFqte21nzcmYYEvGXxsbfsrldPmuEpT/F2FnKV1jBSEepJQT3vbyWe8+JjIM6bJ+i8fWxGhCd1Gny9U3zN29UeKKmJoW0vWX5mf1uFa9tupfXJGBNV9iWvDqo0jth21MZWTdcC94ihXGKv7l8vYMLe68i+SC+uqqnuhSxbKllKjfckpV3+zurM+uXVZdPYQh3uSCwgIOBQIb3IBAQH9jefaJEdEZwAsA4gAtJj5tsvVbQc6zszUSscsSDjie0q9bp/NWKWKwnVy8JwS1PRv1KR/fpK2kBfFNFxtAwDkFq1p0iyJ6VAbsqZJQimKRCpHABWt6UCzcq0z5ywD36iGOOGsXEH6ed1RoXy8aNjq+bdVeODFhqV/HFUJECabI2vbC03rArhhRBQ+Es7LPFWR/BtzTsBUo9wSU3+yanN2PHlRRUPU5bEdHi6b48YHJIrCm8PnFiQSo5aQMY3avSkeSRcJoK3o6lhvk66tPAf1I3Y82kfVw+qunc6LK2F8WMa+4SJAltSzSe41aV6ZpZqmpMcNANo6UsJRRpIrO6vCtp8LD/upJ/cyZr6FmW+LP9+Jjm779QA+GX8OCAjoA1B7/d9e4SCJZl6WbntAQMAhAvP6vz3Cfk1yDOATRPRArPMOON12AF1124noDiK6n4juj8rlbocEBAQcMBz0vKu7gZcw8/k4AcW9RPT1S54RQ+vCZ684zc3Yz5Vesr4NUiqo+WcsReDJQfHvaLGFxqKlfzQy8k49MFw1ZZW0HJucVmoiVlgDxfOy3Sw6v4dSDQZJWflk79uSsmK3SDTkd2rlKmsDtJTPRec+vTpr1aJPpsX1eX3aUlTKiq5xXiU+qbQtXeVcVfx1TyxZH+jUopxXW5HzBkdsZx5dkN+1mSXre2wpqojOeTv/zIg5br4kYXojI/ZmDCiV32ZT5S2N7PPRVqq7lRPOj1pTKs0V6noOYNV6I6e2bPLBpl3eVUXrWKhKKFe9btsYKbVrH5KVVX49kyzIKQMnyioXb8N9f7DD2MeFh315k2Pm8/H/0wDuAfBiXKZue0BAwOHBfr7J7fkkR0TFOFEsiKgI4HsBPAzRbQc2odseEBBweEBtXve3V9gPc3UcwD1xWrEUgA8w898Q0RdwGbrtAABitHMdU6DlhANzM71fuJvzStEh03uZZ1XhBADSo5aBPzQoptb8KRnGqGCHVIta5ubsjc0uyrWjrKItuGQ19SO9xTvbOsdmxZa1aoq+okyyb9QmzHFXpEWlxWuVPF8pVTQb4gM9mbHsnk9duGFt++yTR00ZqeQyo2NCi0i6nKNLKuKk5SIZjKCpovN4BRFt7s3NuORDPR4JcjQObbr5sl6vBc56R3tMaCIJ18+2VoSpuWgIkooaaRkDdoKuSSVu6thNqA2psUur8XHJh1JluXbaPVeRDY7YPp5LPDlmfhLAN3fZP4vL0G0PCAg4PAgRDwEBAf2NkONhiyCsvatHg9acbNaVmXjRrRw1tGCiDH7SRRq01WqUN60GcxJ4v5hXLH6X46E6p8UO7eqqzouqV2UTthmGOLlyjS0snJVbOPiUPW+ZxcTRK8qePX//3BVr2y8cPm/Kbi48K+exXOu++RvMcednJJogWXb5ApTAwPyKrIZmZu1xGZW7NDXsclkMahO19xcmcV7Ge2DSm6GyHalFR+0OAOzKfDvrrtXDu+FXXvl+iwAAFRZJREFU97PqvrMzcXVfolNWwEHnZ2iqVfvMnFuZV6ex+xanlOinFuWMXF90ztfGiC/b4UkpvMkFBAT0M4K5GhAQ0NcIKQkDAgL6G+FNbotQPjnvK2mMiG8j0bDdTC+L/6SlkpVGjo6QKQhzXAtLAkC9JXUeG5PEIZETnZxTn6uwNIDaMZW0RPkQOe1oEUoQNL1kfTPaX+fFGRtHpM068c7Ts6PmuEiJfp4u2gQyfzf/vLXtb8xLRMLUGVtHekGx551PMdKMnSlpR8rpiTYlWAGtovMR5eV+ZqZl7D31oTipkgotWz9telkaximVRGjEjmntiNTZLPbOi6q3faRLS7lpG4OuL8pP1nYJhtpLMlhK/xMZG4iC3KzUURvzyjeyTepepBvOV6r8zG3HHaJWbwrWVkBBNDMgIKCvEXI8BAQE9DPCm9xW0QYQv4JrVj0AKwLoxlczvRPqtbzBLt/ppJiXk6fsUOkg6LTKi1rI2PwJzzt5YW07d4Ute3RGzL+VKbHVfG7VpjJpUku9hRqbNqYdVO8umKjFNAHg1KiY21UXrP7IrLRRRxD4SABNtWged6KfKhA8OSXbTtMSTWWi+oB3HUw+9ITsT9XsK0KqKp/zk1ZUoVWS+9ksyNhklm0d2SUxc6tH7H2vHFeNVt30EQ+NIZXH4Zi9n8k5ZW7Pe2qIFgCQ/QlbBbKmzd4MVe6YgnrWnbDsgLCDDKUGAMondjhEP/jkAgIC+hlhdTUgIKC/EczVgICAfkZILr1FUIuQmen4NPjauilrZsVPQW3r98iICwr1YbV/wfo2tA8jccY6LWqnVP2KIlB2jqZqQ/nTkpbSkFO+vLrKz9pcstfSOWV9uE1jSLZ9CJJW00goeozOTQoACeUw+dLkSVNWUZQGNNX4uDy0UUnleM1Zn1zmCXFYmeHxiXfmVP5aV5idVzSassqfmnYhe2qIm4PWUbZySj63lMpGdtGFdW3whdSPkqbKaJoMYH2KhcdsOzTNpXq0t6qMfv4KF22jtPmXW7DPFSlaVFXRSzRFx5clHZ0ntdOi2/v4JneQcjwEBAT0K7jL3yVARO8momkieljtGyWie4nosfj/kY3qAMIkFxAQsAegdnvd3ybwHgAvd/suO6vfoTZXwSIwWZ13a+BZGUTPwE+05Gckq5Qv1pkcWpDSiRtC0TOqF4Xenqw4ZY15dVzBMd+VSRMNqLyrTW+CKRFHa5kYOkxyzp5XO6poDEoVY8V19BslaX9qyppWCWUC07jYNO2m7SepvAXZpywNZfBp1Tf1bHuzcOlKlZPC9WXgrLqfJvepY/sXFEUlYd0Umk6hkxi0nKmm6SC+jYUpHVEhhZx0Y6+iT4aesvyPylF5sPz9HHxW6hx+SKJPOOv6UpLnPbViuSEplRsiysq90AKuAFA+LdfKXbT305uv28YWfHLMfB8RXeV23w7gpfH2ewF8GsCbNqonvMkFBATsOoh5/Z/KvBf/3XHpmjaX1U/jcL/JBQQEHA50WXjQmfd2E3s+yRHRaQB/AuA4Oi+xdzHz7xPRWwH8LICL8aG/wsx/tVFdnAQasbgiudU+rPQ2CVKKCF88L6uy9RE7HPVhJa7pqs8qprpmwSfr/jhl3pyxN3r5tNTfiHoz05P13uzzwrQyed1KY2ZRmblKdNFre7WUidccsHZFe1hMrYSKckg4833oG/L5+CeeNWWcU3kLjkvURPm4NY21eyA/YxuZVpEMlSPSXr+6GqkqvftBB83rQWg7UzPZQ1wTAPIXlfmuhsoHuJfOyUOXXrI3tD0hHT3yNWvK5v7+a1L/iXFp+/CgOS49I+EQdG7Slj3vyrXtgjJDM0vWcKuN6sgI2/76yA6vhu7c6uoUEU0w8+Rms/rtx5tcC8AvMfMX46xdDxDRvXHZ7zLzf9yHNgUEBOwmdo4nt5rV7+3YZFa//UhkMwlg1aZeJqJHAJzc+KyAgIDDjE2uptpziO5GZ5FhjIjOAXgLOpPbZWX121efXLxyciuAzwF4CYA3ENFrANyPztvefO+zAwICDg22YK4y86t6FF1WVr99m+SIaADARwD8IjMvEdEfAvj36NAE/z2AdwD46S7n3QHgDgBIjg6vMe29CgnNKeeMG9/cvPhLsjPKQZewjol0ufevD0WaSiDDWB+0fo+BSX0t65spXBDf0tzzxfmjoxg67VBNdD6/VFXaUZy0zsf5vFIeUSKLqZodkPqQ+GbmX+DY/zrHqRLe5Jw9LiPpVFG7ZsyUJRoyjsunVG5Vl9tT+0rzM27s9eWUCy294iNAFCXI1c/JHv40pyAy+IwULl1h72d9uLuCjaeaJKtyLxrD9gKZFaWU8uScKat/q4iURjm5dnbGcjr4649LO05bQyhRkecsqXLZJuou4mZW2rF4rW1j+eSB9cldNvaFQkJEaXQmuPcz80cBgJmnmDli5jaAPwLw4m7nMvNdzHwbM9+WHCh2OyQgIOCgod3lb4+w55McERGAdwF4hJnfqfbrtO4/DOBhf25AQMDhRDee3F5hP8zVlwB4NYCvENGD8b5fAfAqIroFHcPkDICfu2RNBCDWqeeqEx9UtIuBZ+3Phjaf2hkl6FhxeVdLwgtINF3+h1kxH3IXFL3hdMkcl6yJiZBasiZHoiUmwrH75doLN1qzWdNXitO2jelFoSBUJizfQQst6iiP/LQTcVSCmpUZl+9A5VYojIg9WZu3/Vy4XuVMGLXt0CalDnBPWU1L08baiBNL0E1W45Fx5mqyrj/bOqJM94gHL0iZUaKZeZezV4tJ1oeUOblgn4+kel5aRfdsNqSNK88/YspaykTV40Gjlg+Tv/Hate02uSD/nDy3mmLTGHBRDSpqxdNo2gXHu9ounktSS8z8GZhHbA0bcuICAgIOMaL901oKEQ8BAQG7j+fSm1xAQMBzEGGS2yIiQnKx04XCpLWAKyfk9dj4YgCklJ9MJzfxShLapVMftXE72teWLKvQJ+e70+c1S5YbkqwrKsGzkrizeMEu5zeVLyXpE7csC6ekGNkHKVUV31h6xUmxKOj8odT7MFSXxS/EOduO2kk1pgWX51blRq0dl+M0paNzbRUeV+4tYKr9rT7crjitqBteOUaHpSk/XLLp1GGUH8uHwGUXFBVH+RQbJXutlqLv6OQ6gKUftdO2n8WnRdE1MS+8nNaElU0rXy3PErtvcb2klVh0m+xxOp9v7agty1zsnTBpSwg5HgICAvoaHHxyAQEB/Yyw8LA1UBtIx6KRGafTH6kcD03LdjBIVsU+qx211IdGSb3ru7ft6rgSLaxoYUKXA1NTFVq2knRZrs0ZZd6sWE4DKXmO7LlF1wG5XvrCgi0qK/tEPWO1k1bsv3pERTI4NY3Eklxb54zw6+Oc4q7bANBQAhq5cQnfuGLURu3VI7nW5LxV3ajPi6mcXBFTKu3M2mWVH9el9kBOKcLo++IjHqKsFPqcrNpU1jlfEw3XZ/Xs5OYcNSlFXbcBYOVqeVhTx4Xs3srbfup75sUw9fOeFi8IUhUv2irbPmIjP73TeVeDuRoQENDPCJNcQEBAXyNMclsDtYFUrB2oA7MBoK0sz7Zj1jeL0u38M7KClSxZWy2d6G2eJdSKXKsopoMP0E+rAHodaQEAqUUVbd9S+QKKth3GfE30NiN4Ycl8TrTUSqbS+k+vWPY8aRvVPYuZRelPSnQa1wtSKsu4OWbNs9yIRHocGxT76YZBq3d4Mivm9gP5K0zZg9EpdS2Vcq9uXQw6QiHnhDdTKhpCrzomyr2jJlJly/xP1qRvOi1gctkqJ6xcJ6uf1aP2fraUOexFOZtqlVb3M+1SBGoXQNu5B3Kz3fOWVMfNYSD1WHlz1YtEbBtbkFraKRzqSS4gIOCQILzJBQQE9DXC6mpAQEA/gwNPbmvgJFAf7f4arJPXZKyrCrmL4iPS+Sw9u13LwXiKgF7S1+oOmn4AAGm1bJ+qWP9OfSyvysTXk3nWUkGQVrdp0uXtGBeqevtaK54Y5cUXpCk1mUUr3qnFNpsD9pFoKZ+OTsrTGHS0BSXtR1nbz2ZD6jw3M7y2PZZfMce1WO7FWNY6oV5y9ZNr22Wl3vKFyrXmuIyivJRP2jbmZmS7MKMiVlxEgo6ISZ93lB1tdmVkfDlhfbE5JZC6cJ0LNVBIO3+gjqLQPj+TMxYwHJh1FBLFEGoMK8UdJ3RqNN0c3YZTOzwphYiHgICAvkbwyQUEBPQ1wurq1sApRvNovA7uqQ+TYkr4nAY6EL+tqBVRxpoc2gxNulyoOnohuyg30C/1a3GAxpClEhhRREUhoaalYPCSMuuylnPgjzX116WsOSB25/JVVpQzsYE+oqYgVCZUQPcxe1JmXPgl147ZvAWFlAzedEXo+AnnHziVlQiIc3UbkJ5QN7iUErrG2Clr2s+vjEqb5nvnZyAVDjH6NStmmppSJmrDRp9wQQ2IcqaTOw5DclzBCZ1mluTYKGftxOaAfNaB/Y2SPU6LHpBNu2pcCS1lykYJLx4r40NjdgzS6SCaGRAQELBpcLTDk+ZlIExyAQEBu4+w8BAQENDXCBQSARG9HMDvo7Oo/cfM/Paex7YI6amOnyuzaJfRR74hr8fFs9ZRRk0lrDhhFTkM1I/POuFD9cukRRa9CKIu84KaGs1h8bURu5iatvixdNsBoFUS30+UtX6bViGptlWimRE7VlooM+F8j1Fe0UbySuTzhKV/DBeF+zCctXF0oxnx153IC5/n1oGnzXFJFVt0TdZSZQoq4ez7pv7J2vbMtFUr0UyO5qAT5VzoHhLXTjnlmIq0nxw1RNNGjL9u3lJN0mWpI3nM+hdr6pnzKiSswvZ0chkfdpVdVLSfovM9qgQ7nFAJjCbt1705IGVR09ZRW3LSLNsE7+Ob3L7kXe0FIkoC+AMA3w/gJnQyeN20v60KCAjYNri9/m+PcKAmOXQSSj/OzE8ycwPABwHcvs9tCggI2CY4itb97RWI93Fp14OI/iWAlzPzz8SfXw3gW5n5DeqYOwDcEX98Afo3CfUYgJlLHnX4EPp1+HAjM28gPXuwcdB8ct2cJmYWZua7ANwFAER0PzPfthcN22v0a99Cvw4fiOj+/W7DdnDQzNVzAE6rz6cAnN+ntgQEBPQBDtok9wUA1xPR1USUAfBKAB/f5zYFBAQcYhwoc5WZW0T0BgB/iw6F5N3M/NUNTrlrb1q2L+jXvoV+HT4c6r4dqIWHgICAgJ3GQTNXAwICAnYUYZILCAjoaxzaSY6IXk5EjxLR40R05363ZzsgojNE9BUienB1uZ6IRonoXiJ6LP5/5FL1HAQQ0buJaJqIHlb7evaFiN4c38NHiej79qfVl0aPfr2ViJ6N79uDRPQKVXZY+nWaiD5FRI8Q0VeJ6Bfi/Yf+nq2BmQ/dHzqLEk8AuAZABsCXAdy03+3aRn/OABhz+34HwJ3x9p0Afnu/27nJvnwHgG8B8PCl+oJO6N6XAWQBXB3f0+R+9+Ey+vVWAL/c5djD1K8JAN8Sb5cAfCNu/6G/Z6t/h/VN7rkQ/nU7gPfG2+8F8EP72JZNg5nvAzDndvfqy+0APsjMdWZ+CsDj6NzbA4ce/eqFw9SvSWb+Yry9DOARACfRB/dsFYd1kjsJ4Kz6fC7ed1jBAD5BRA/EYWsAMM7Mk0DnQQRwbN9at3306ks/3Mc3ENFDsTm7atIdyn4R0VUAbgXwOfTRPTusk9wlw78OGV7CzN+CjvrKzxPRd+x3g/YIh/0+/iGAawHcAmASwDvi/YeuX0Q0AOAjAH6RmZc2OrTLvgPdt8M6yfVV+Bczn4//nwZwDzqv/1NENAEA8f/TvWs48OjVl0N9H5l5ipkj7iQV/SOI2Xao+kVEaXQmuPcz80fj3X1zzw7rJNc34V9EVCSi0uo2gO9FR1nl4wBeGx/2WgAf258W7gh69eXjAF5JRFkiuhrA9QA+vw/t2xJWJ4EYPwxRxDk0/SIiAvAuAI8w8ztVUf/cs/1e+djGqtAr0FkJegLAr+53e7bRj2vQWa36MoCvrvYFwBEAnwTwWPz/6H63dZP9uRsd062Jzq/+6zbqC4Bfje/howC+f7/bf5n9+lMAXwHwEDpf/olD2K9vR8fcfAjAg/HfK/rhnq3+hbCugICAvsZhNVcDAgICNoUwyQUEBPQ1wiQXEBDQ1wiTXEBAQF8jTHIBAQF9jTDJBewoiOj1RPSaTR77k0R0kYj+OP78UiJiInqdOubWeN8vb7E9nyKiFSLqyyQzAZdGmOQCdhTM/F+Z+U8u45QPcZyCMsZXAPyE+vxKdDiEW23PywAc6mxTAdtDmOSeoyCifxQHlufiqIuvEtELuhz3z4noc0T0JSL6H0Q0Hu//T0T0G/H29xHRfUSUiDXWfjne/0Yi+lp8nQ9usmnPAMgR0XjMxn85gL9W7fk0Ef0eEf0DET1MRC+O9w8Q0f8T6/I9REQ/ur0RCugXHKhENgF7B2b+AhF9HMBvAcgDeB8zd0vU/RkA/5iZmYh+BsC/BfBL6GiMfYGI/h7AfwLwCmZud+alNdwJ4GpmrhPR8GU0788B/BiALwH4IoC6Ky8y87fFQgbvRifJ+K8DWGTmFwLAYREZDdh9hEnuuY3fRCcOuAbgjT2OOQXgQ3GcZgbAUwDAzBUi+lkA9wH4P5j5iS7nPgTg/UT0FwD+4jLa9WEAHwLwPHTCqb7Nld8dt+E+IhqMJ9DvRse0RVw2fxnXC+hjBHP1uY1RAAPoKMLmAICI3rYq5x0f838D+M/xG9LPrR4X44UAZgGc6FH/DwD4AwAvAvAAEW3qR5WZL6ATI/o96MRNrjuky2fqsj8gIExyz3HchY6Z934Avw0AzPyrzHwLM98SHzME4Nl4e1WVAkR0JTpm660Avp+IvlVXTEQJAKeZ+VPomLjD6Eyom8VvAHgTM0ddyn4ivsa3o2OiLgL4BIA3qOsHczUAQDBXn7OIaR4tZv4AESUB/AMRfScz/5079K0A/oyIngXwWQBXK3meX2bm8zHl4z1E9I/UeUkA7yOiIXTesn6XmRc22z5m/ocNiueJ6B8ADAL46XjfbwH4A+okmokA/DsAH+1xfsBzCEGFJGDfQEQ/CeA2Zn7DpY5V53wancl107SQrZwT0D8I5mrAfqKKjqn7x7t1ASL6FDqafc3dukbAwUZ4kwsICOhrhDe5gICAvkaY5AICAvoaYZILCAjoa4RJLiAgoK8RJrmAgIC+xv8PwWWu95Rm00wAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATsAAAEKCAYAAABt+vLPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9eZydR3Um/Jy7972974vUau27LFmybGxs2WCMCSYQAglkAmYJDiEk4ZtsbDMhk49MQgJMvsyEGRGYkARsNjtmMWDjBeNdm7W2drV635fbffelvj/u1T3nvGpJ7V6kbnU9+t2f6r1Vb731VlXXrefUWcgYAwsLC4vrHa5r3QALCwuLqwG72FlYWCwK2MXOwsJiUcAudhYWFosCdrGzsLBYFLCLnYWFxaLAnC12RLSUiJ4molYiOkpEf5T/vpKIniCiU/n/K8Q9nyKi00R0gojePFdts7CwWHygudKzI6IGAA3GmP1EVAJgH4B3APgAgGFjzN8Q0ScBVBhj/pyINgB4EMBOAI0Afg5gjTEmMycNtLCwWFSYs52dMabHGLM/nx4H0AqgCcDbAXwjX+wbyC2AyH//kDEmYYw5B+A0cgufhYWFxYzhuRoPIaIWANsAvAygzhjTA+QWRCKqzRdrAvCSuK0z/52zrgcAPAAAbpd3e7CoBgCQ8et1O+vjtHEs6cYtLly8s/X7U6pcIsqVUNZRh5+/cEX4AdmQLkhJzjNevYuuCk4U0mP9JYV02q+fBbe4z62ziDjPNeF8UVmQkxlHG5HhTBJpADA+UYnYY1Nal0OA6zSOPFeCr7Ne8X1aVyHznCBR1gSylyznion+djTRFeJK3KLfso6C6bToZEcd7nH+IlMs+sbxzu64qMLRXHd1kuvI6jELejkvleV2xGI+XArk0fPKpLlOT0Q0MeC4UdxHKce4iy5IdnQOGmNqLtmAKSDbu+YiCumqP0mTlZ0rzPliR0TFAL4P4BPGmDDRJd9vsoyLOsgYsxvAbgAoLW4yO7f+HgAgvKJIlZto4gHP6CwkKvmv1gQ5vXZFjyp3el9zIe2J6uZlVkcL6eDLwUI6enNUlaN2fni6IaHy7t/Ka/uP/mFXIT26Vrc3XcZtpKBeITxezgs9H9I3it6TC/7E63QbM2P8h+Qd0atpegm32UR4uvgHdDmzlv+qkiP6r6rkJN8XbeBGBYZ0n8bqxILp+AOWz0utjXE5x0ISPMz9nXGsD4GdQ4V0aYDfK5rSq+xAX1khTR79gMpn+Zdo+Haug4b0wyqO8rt5o/pdSn+ns5AeievJuaO2o5DujnI7Dh5dpspJThaoiqmsRB/Px6p9XPCieVXNP+6BTt3+ZBm/d9sn/uQ8rgPM6WksEXmRW+i+aYx5OP91X16ed0Gu15//vhPAUnH7EgDdc9k+CwuLq4PsJP+uNuZsZ0e5LdzXALQaY74ksn4A4H4Af5P//1Hx/beI6EvIHVCsBvDK5Z6RLnJhaFPuVyzr07uERBX/mhqX/mX1hgXtrOZf51MHl6pyxV1cZ3ijprguQffCW5h6lOwNqnKpYvHc43rH82+DdxTSy85zHZ6Y3mmkA7yrMS4nj+Vk2Tm9cxzayLuQsQ28I1z5Vf0bN7CFp8Hyd55ReYfOsySh+CyXyzioduhnvKscb9F5bn41VB7ltOwbAEiW88t4hnQb4zX8x1H5DPdj9XvbVbk2f2Uhne7QO93EmcLBP0YFPQ/URVS5t205WEg/97UdKq9m7xjXV1leSEc3693VyF2c9p7Ru7exp5gx1BzUO/Wf/Hop39cndls1ev6FKvh50bCeV3JXPLJRiDmadBsrn+a5OrpeLz7ZEoeMYYZITXLOeFVkaFfpebcBeB+Aw0T0av67TyO3yH2HiD4MoB3AuwHAGHOUiL4D4BiANIDftyexFhbXB67FTs6JOVvsjDHPYXI5HAC88RL3fB7A5+eqTRYWFtcGmXngSu5q7yQtLCwWIbIXnzVedSzoxc54gHhVbvOYKtGd6W5htY5UQr9meogFTiEhYxtfp+UiqRDLjSimZUjeHiEnWcfPiu9MqnI4y3KjeKNm5WtW8/lL33mW43gn9Luki3iDXNKpZSmRepbhDa/VgjSprlH3S3Eqt1KXC4zw847sWa7rqGM5YPYWllelzpSqcvGkkLfpw151EhwY5T5IlupxSdVy/zfeOKjy2k/WFdKbP8yCv1+eXanKlTzP4xmv0u0I9gl1Ey83KlKp5aBtEb4xHdTkpH8nn5BGG5maBYu1vPTvNn+/kP56y+tV3pGf8rFo+32ONp7isZGqOO7z+rTUdYyvXfV6vviWjxfSMTfLC7NDWrY3cptoc1jLiYvOX1rVZTrI2MXOwsJiMcDu7CwsLBYFUlZmNzNQFvDmtQbSWuNDUdeyMoeibzmrGgwHmZb4+nV3SM13/7CmOi7BAGgfWz/U7tN0pu8mQe8iuv7OGn52fIWwyEhq6kRZSb90HcE+vm94g8pCw82sJH2+vbqQDrRritLwArd5bJXOI6E2UXvbcCHdVqzVOiKC7tW8oik/Zbj9Y8u5/Z6I/gPwDDKV6j/bqPL8ovt/+fzGQto7rp81tpZpss+hIF1xkt9zdBXTxfSwpndHBlmB19WgTxEzq/l6/VpWDv7mqu+pcm85/P5CemTcoY60WphXRB0iliLuE2+E+zS6Ia7KxYWVRF39qMobi/CYeQe4/nSZfpclD3PewFbdV/H62VWEsDTWwsJiUSBz7dc6u9hZWFjMPa69lp1d7CwsLK4CMpdUub16WNCLXbYkg+QdYQBAqkfLkEr2stwiUqNlMvJI39fC8ru1aztUuSMHWXbT/BP92+RO8nW0litMVOguja4WqihOzxgJvq90OctdRof1u3gGWI6Wdjg16H0T65eQW7exe4hlgptXs3ypYcuYKveMb2shnazWshpXiNVBOgfY3CpQ4zA+75cqH86JPflE92gxFPzCMUDlCa1iE6tgmZJL9OPwRs2PXFXc39L8DADOC3lb2dP8fVGnllf5R4WMcbWuf/vGs4X0siDLMB+eWKHKSW8mTtWnsnKWIY9GtQpPWnhSSTZx329s1k4q2n8oVISeqlZ5JW5+77FV4vvT+j2HhIw3Vazfk5KzuzilnC5orgEW9GJnYWGxMGB3dhYWFosCTp+B1wILerHzuLKoLsnR0KIKffx+toq39hWPa+4XaRAOGIUTx4YiTe8iG5g6tEe1H9F0FdOsskNMWcKrNJVcs6y3kO6f0G4+gn6mXBNxVoVwjWlt9vLj4sKhrxSv5SFs2tCn8oYjTC2PHmjhtENaLF3HecY01ckk+N3KDnO/RZr05PUJB6MurX2D0c3cV95yzkwd0XTdP8JpSVsBoLSd+yrczLTevWJClVPuEo+WqLx0UFDEUuFQ1OHBZeQG4T+wSNPpfUeYru4PstXLjwKbVLm19f2F9B0Np1XeY4/ewvXXaLGBqeT3dAvxxdH0ElXOLVRikmV6LKQ4QHpAGW/RA2+EM9NAh1Y5StTO7pGC3dlZWFgsCmTmQSBDu9hZWFjMOSyNnSFSSQ86zufo6ooVmsJ5jzKFc57YZQNMHQKHmeo8dXarKveJX/1RIf3yLk2Tw0mmxr1NXEf4tLY+T2S4iwM+7Wigu42p9jtv2ltIv1qsKcs5V0Mh7bQKyAoqMvSLBpUnldaruwSdWeqgoGw3jqjDqFzG6Aiv4vtC2mcm4OK8uCNagWec2xw6HLpkudEbuH+8pdqhgusxHs+iER6/8AktGpAn7aEBXX9JF983uElYcjgcF3hHub1uh8vz0HPMeSnD5UbXaC7cWcQn4YdP6/FccoTbMbRBj2dMxNeQp87eTi3aSAlriGSFppyupGh/jOsIrNdzuCrEL97bcVG4l1lF0rivXGiOsaAXOwsLi4WB7Dygsde+BRYWFtc9MqCLPlMBEZUT0feI6DgRtRLR66bbBruzs7CwmHNknPFMp45/APBTY8y7iMgHIHilGy6FuQy483UA9wHoN8Zsyn/3bQAXPBeWAxg1xmzNx5VtBXAin/eSMeajV3yIIVAq14k9T2u5SHy10H9wqHIYoSYhHVymK7SawdND7GRxLKnVV84eYhlH1RoO0eeMDdt+nJ1OfvwNT6i8f+xk7/Q/PsOePALPaZUJrGYZj9thdeAd4iEkh6MKWdYbZbmON+KQ+4nuqTmg5T+jq2SwH/7eqWITGOBMn9bgUfVLOZ0vrMuVtIlQjQ5vJiW/OFFID79lDWes0cFysh0irKVDhLn8N88V0t0vs2lBqFvvMgL9fJ2K6rGIczwfJCvEWDdr2d7oeQ7GgyI9MANb+d1cDl+vRee5s6QsUQaQAoCyE9zGdEj31cQyHpvSk8Jpa5+Wb6Z62CKmbq+WJ0808ljMRhzF7DRUT4ioFMAdAD4AAMaYJIDk5e65HOZyZ/cvAP4ngH+98IUx5jcvpInoiwDkn8UZY4w+IbCwsLgukDQXLzUy4H0eu/NxoS9gBYABAP+XiG4AsA/AHxlj9C/cFDFnMjtjzLMAhifLy4dZ/A0AD87V8y0sLOYPsnBd9DHG7DbG7BCf3Y7bPABuBPAVY8w2ABEAn5xuG66VzO52AH3GmFPiu+VEdABAGMBnjTG/vGItLgNcUL3Yphd7byvTj8Cg3kJHlgjteUGxis9ountwZHUh7V+jORfVMk0ebGNu40o4HG+KGA7/dOgOnRcX1gm/ZIoRrVXF4B9gKinVRAAgJezIZWxVAPCNcP392zkdcKhkSK2A8HJNcSV1rX+Z3yVap/tqaDPTrJjD4aURMVo9Ya4/vkxTp5KjrMWf9enf4eIYc3Ij1FySEa357xHqGimHs8p+QUllnIbxIm3JUXaM2xi5UdNTI+KXBJrYemNrQ5cqd7acVZAqi7Ruy+liVjlqqNKcv62T83xdwgFESNPY8eXCSoJ0XnAZz9VUM5crdpTzHGKqHavSS0HGN7t6cZnp6dl1Aug0xrycv/4eFuBi917oXV0PgGZjzBARbQfwH0S00RgTdt4ot77uqnJntoWFxTzEdCwojDG9RNRBRGuNMSeQC8F6bLptuOqLHRF5ALwTwPYL3xljEgAS+fQ+IjoDYA2Avc7781vd3QDgX75kHvg/tbCwuBKy0z+N/QMA38yfxJ4F8MHpVnQtdnZ3AzhujCk4WCOiGgDDxpgMEa0AsBq5F7OwsLgOMF3bWGPMqwB2zEYb5lL15EEAdwKoJqJOAH9hjPkagPfg4oOJOwD8NyJKA8gA+KgxZtLDDYU0wd2fk2vE4lrWVDLKMgJf2HFs/zyrAgyv5S5IbtceNEyGByjuiJOareET8C2b2grp86MVupyQVdy/6mWV978O7CqkYzUivqcWhyFVKkyD1miXIsEQX1cXa7ll90sctMY3xu2I1er+qNvL9U806H6MLBOmTZtFTFPHzMkEuZxT/ab8EBeWh3LhkH6WHKeiQS1v6/hDPqiPCY8cJZX6nTOnuf9jjVom2DvMY1i0n1VUms5o1ZBYNbejpX5I5Q2KQEPjYR6zowP1qtyOenYE63LIyjr2thTS55dplSZXBc8rygjZckC3cfUWlhHeU9uq8qo9LP05GWf9m+cHHA5Gu7QqikS0fnaXhtT1bC5mjHnvJb7/wCTffR/A9y8ubWFhcT1gBkrFswZrQWFhYTHnmI5S8WxjYS92LoN0aW57X7N0RGUNJlkdpOKEw+uJV6gnCN/7qV5tiWJK2KKidLX2GOH3ct7JATYLiI1oWuItZlpS5tYqCG4P07FVb2ctnENCux8A/EP8q5hO6HgaSfC162ZNw5MNKVFOIKsnXvnHOT7FwHMtKs8T5mePLxdU1a9p5vo1TKtOdNWpvNEtoj7hHNQ7rH/tU8VSnULnRZYJ6xbhkDKb1eXK72BnqdHz2gONv437amIj0/94rVZf8a9iGvi2Gu1486ep9YV0zMv3+TyaZv6yjSljUUDT6cR2Qb0H9Xj6TvP8kRYw6XFNA8+8wo5DH96qPa7cUcdtfujgTVz3eV3O7YgzrNpYOcvOO+3OzsLCYjHAOu+0sLBYFLDOO2cKtymE+hsa1idLgUamCv079ElqzQFBiYQheV2pPtm7q56pZU+iTOU9+xIb7melw8W4/gVbupzp9cHIUt3GIiaXB060FNKlDsN0l2BBKf0q2HAra+i09mr6WF7DtHZsjCn6pmXdqtyK4sFC+uQq7VEz3c+0qqSZ6d14p25I2xCLDcyopoXyVLGkjb+PaN8NCG/k/gid0nWUNrDFg9vF/Z14UVNV3xv4nZct16Yifd18Ou0XMRcSjlPb7Hm2tHgytFbl9Z1hC4ed23h+lHm1pcXjnZv5YlRbaKSLhCMKh/OGxHLBXSf4WH7Jqn5Vbkslj+FgUte/s5jnRHQzv+crP9EaHAPbeFw8MT3nvOHZ3YmlJrGNvdq49i2wsLC47mED7lhYWCwKzMCCYtZgFzsLC4s5h93ZzRgEkxd8urv1sXrCz9cVHVr1pP9G4fDyBMs7epdq04Wf/Oj1XC6j6zAbWW5UfJrrcztcC54rZxlYmU/LdRpLWQaWKWFZ09lxHfyEhBaAqdMWFH1RllXe3Nym8o4PswyvsZZVZ2Jp/Z77B1mWeFvLOZW3ays7zfxi692FdPVybeASS3KdoTatJuF0PHkBQYdscvudZwrpl6JrVJ4vKSxdokLlQ2tuoO0893dxlVb1+ZP3PFxI/9/ztxbS3Se1nDIr1Go6zui8QB3Xub+d+y3brVWOUMlyQM+4/jNLh4Qj1VHdV0k/96MRqkkDL2pPpK+0scXG0FbdvwcaWBi6ZDfX5yrX5Vwp7v9lu7SLzpOHtXx5prA7OwsLi0WB69pczMLCwuICrFLxDOGKEkKHczym8riOH9H5m0wjRjZpiitjrbqFM0lPv1Z3cCeE08mYpgClp/i+lNR6cYgmyl/iOqPNun45ATqG2Ddf5RpNEQc7Oa+kRAeh6G9lmjVgNOWq38jqCu9c8moh/fOBdarcigpWPTk4qCm0VLGRahL33qWdGsQy/G4/3rhJ5VU9J5xQBrmDElprBD1Rh16NhHTGKjRFqg87xr2Kx2ViQKtk/G38zYW0EXpf3jGHJUc508dAv4NmlvO1NPB3Om0t2cdzLukIKQJhmZOt1uNZHuLrO5qY1j/1rZ2qXEDEznVam6QGmFKfu587K3BKv4uMcXHiVKPKK+6c3cXJ6tlZWFgsClgLCgsLi0UBu7OzsLBYFMjand3M4EoBoZ6cfGVsuX6V4KuXfjXj4Y6XclOn6c7wHSw/caq2lJ5leY2Mf5os1b9giQq+PnFOqw94gyxPqa/kSjq6K1U5CrKMp7xIy3iab2SVko5RHZNjYC+rnvxj1xsK6aI2LTs8HVpWSP/9r39D5f3lo/cX0qNr+Z3PRbTA7dQQywv97bqvYiKAkNkuOuuQw+TsHBf0VGsVm6Jl/N7hc/yeHe/Qg1Z6kJ8dXudQF2pnk7ni8zwuo9u1vpDsn2C3rqPkPMvD0u9mx56lN2mvOEM/YdlnxSndxsa7WZY6ntB9FX6e++Cpm3hyTizTdaSLeH77xnQbk81CqDnKqifJdVr1aWdLWyF9LqznXK9vduO7pLJ2sbOwsFgEsHp2FhYWiwLWgmKGyBQBI+tznejTLAIuoZEwsV1v31c0sKpF95OsKd7wki7XA6YsGc02FDXLbGbrh7pyHdi1a4DpgLtPV5ISMU5HAsJp45C2cHALjxQDZzQVju90BIEVyK7g93F3s6mB1xGgcsnPWQfhswMfUHmJpUyRTBXTvYMvrlblPFGhUtKsaaGvm9/n5qb2QvrF0xtVueLTXM5zm+5H6bWl4hg/K1rvoIGb+NnLmgdV3ugPWb0iKVhay3dUMRjhZmaiSf+JROv42e9sPlpIbwtqC4Q/rXufeJbe1WT6WATg8Wp6Gl/KE9c1zmNmAtqZZjrEdaZqtdeW0nIez1i3eFFHXOGzY9yOxI90sGJP8+wG7pvuAQURuZGLMthljLlvJm2Ys70lEX2diPqJ6Ij47nNE1EVEr+Y/vyLyPkVEp4noBBG9efJaLSwsFiKyxnXRZ4r4IwCtVyw1Bcwlkf4XAPdO8v2XjTFb85/HAICINiAXdWxj/p5/yq/oFhYW1wGyoIs+VwIRLQHwVgD/PBttmLPFzhjzLIArh0PM4e0AHjLGJIwx5wCcBrDzCvdYWFgsEKSy7os+RPQAEe0Vnwcct/0PAH8GYFYCYlwLmd3Hiej9yPHwPzbGjABoAvCSKNOZ/+4i5DvkAQBwV5ch05JTSYik9LpN7qy4R9eRSPNrR5tZRjLao11oVB9hWUjXLt1V1ZtZfeD+Zdz0docNVKpGbFAdAU4eaeVYqOk0l8uUajlORphYecN6wzt8nJ/nH3QEn+kT6jETnM749NzxHGc5mm+DNiWLinCooVJW/0h2allZJiBkPEndDt/GsUL61CirqBSt04LW2AmWL6W7HKZjQmYVq+X+iDdpedVbbzhcSLtIv+er7SyXigmzski9lpHG6i6965hYwfPl3w/z7/FD3u2qnAwVW35Sy7/6lvKzf2vjHpX3wqCO7XoBZ7ur1XW6hPvYqUokRbJu8WoyDjIA9LWzuol7pe6rTGh2A+5MJrMzxuwGsHuy8kR0H4B+Y8w+IrpzNtpwtc+DvwJgJYCtAHoAfDH//WSza1IJqTFmtzFmhzFmh7skNFkRCwuLeYZp0NjbAPwqEbUBeAjAG4jo32fShqu62Blj+owxGWNMFsBXwVS1E4B0oLUEQLfzfgsLi4WJrKGLPpeDMeZTxpglxpgW5OT5TxljfnsmbbiqNJaIGowxPfnLXwNw4aT2BwC+RURfAtAIYDWAV65YYdoF05+jUx6H14msYHserVGCbhGYxtPA1Ky0XVPEwc1M1dIVmi55RdCXPeHlhfR/rn9ClXtsgj2AbAh0qbzIGq7/J/s5QEv5YU2rpNpL5XFHDNIybnPvHdoDSDooqHGA+6f8uCqGiTtYjcTB/OALCwr9Mw46VDKiN97DG0VQnXqtNtJcznR1KMYqJM7pLp9d1qqnZqyOnxevZ5rv79PlHnuRRQPumP4tXxZmq4xkMedVHNW6OBPL2U1J/426jrpmFkP3dVUU0hmX7o9QL7+dK6XzyM3XK/19Ku9MEdP8F8/wvAoEtTpPxs9j7V+q50SZj6/7/SwOyCb0/K7cP/n8AIDAvay2046Z47pWKiaiBwHcCaCaiDoB/AWAO4loK3IUtQ3A7wKAMeYoEX0HwDEAaQC/b4zJTFavhYXFwkN6BoudMeYZAM/MtA1zttgZY947yddfu0z5zwP4/Fy1x8LC4trBej2ZIVwpIDCQ+8VIlmqqUHOAryca9a+KeyPTllSKt/IDW4OqnKSPgS5NLTsTfLLXWcynWkMJfWjy5ZbvF9JfG3mdymsdZTpNIu6B516t+R89wCeuyRJNRdKCfix5XE+ogRsEdT3B/THeooohXs3TwKtD56L69T2F9OCzbL1RflrTKv8Qd1b0uDYiD9zGVh4ZYRCeSut3KeMwrBhfprKQKhHjWyIcUp7QJ+juuHCuqVm9ovXK+aVXt0NuQjKOGBcDw0xx3aPcb6Fjen5IcYN3TPdV6BDPkb8cfpfKo1oWq9T9mPt0cKujIcvYSiKe1M9O7OX5SMIRKTn+2kduY1pv0vpvJHZAO4KdKexiZ2FhsShgFzsLC4tFAbvYWVhYLApMxTxsrrGgF7tswCC6KjlpXu/tLIMIdjjkEV0cIYcq+f5oo9a7kEFUAoN6sOJLWeZDYe7G9oe0Bvw9y/+0kG6+UaueDEVYRlheyZ5T7mw8pco9meVnD2e1hYbU1K84rh17lpRz/YM3cXvdpVpVIbCXPa5kHRbJHR38vPpzQq4Y0QKxaBNr8Vcd1PLT1rCIAbudrSnSKf2wiSX8nkX9KgspEbRGBstxeqPxC5UYpxrNwA0s26o4wf2R9TrmRzVfe5ZNqLxEhN8z0MJ5yQkdVSd4gl8g1Vih8srP8rMrT+i+GlnNY2bEC7hW6HZ4REzZyLCOWesR8k1XI+td7WzWSiQeF7fjQO8SlZc5N7vOO9PWeaeFhcVigKWxFhYWiwJ2sZshXHFC8FSOVvhfN6TyRvqZViQqHXymnGlcMMTH756zmg5E72DqsOZ2bb128HlhdSCYSLxaD2r5BlYjGY1p9YHxEaYsf3Lz44X0Or9+ll/oUDzm8CYw2sZ0o/1erTrjiXBb3BNMGYvqNSUaXy6oWZ+DxwqVhPF3sWVE75h+FxrnPp5YouuIrmBRAUWYd5qInn6h7WydMNZRpvKKl7C6kEdYrxTdPaLKdbexwbxnVLfDLcJaTDRxnnHpdwkMi/nyZLHKS27gwW5uZuuH3jFNYwdfz45CU7oKRa99YU1jA8NCRWgZ973rsK4/IahqcETPudr93N/jS3lOHGhar9uxmfs0Nu6QBzRrUcdMYexiZ2FhsRhgDygsLCwWBSyNtbCwWBTI2NPY2cNIp5bxSCeXLodHFDrDMpqI0BTx7dLeOkIBln3IuKgAkPWxzCQbZCHMr9yzX5UbTLDA5uSIrsM9yLKyLz/21kJ62VYts9ta0VlIv3HJSZX36JFbuL6Y4z2FbMg7IdRofqj7KrWM88jhRdAVYtlNiYhZG53QMh53nCezUw4VF2ZVm246V0gf7dLBgyIx4YSyWKu2jA9wP3pLWPg2FtZySt8Qj3uqXMtqAyLPK5yZZj2637wxvq/8lEO1yXAbe881F9LRRsc7C1+bqTKHs9Rq7kf/AW1eSOK1pRqQ2+G5xyu80USW6vrjldzfWWFJlqzQ5bL93HfVe7V8M9I0uzsxK7OzsLBYFLA01sLCYlHAzG5kxmlhQS922YBBdHWOZpQe1n74J8TWfunjmopkvfwr0+1jOpao1N3hFU4iUw5tf0ldPaVc/3BS06qWIKvE7OtaqvKkVYZUTzh7ul6Va9jCKgIDMa3HUMqsEEM3OzyRlDLdy57m+0o6tVpBRFg/VB/S9HFsgtVxercxVfX26f72jfK7FPfoOmK13K8yLsTntv9QlfuLH767kA71OqwaGvi+rKCurozeMbiSgpInHXmiWdLywrh0uaKBjCin2xEc5HaMrOE5kS7R7hfdZTwWTmlVup/7lGjrzF8AACAASURBVLZrNaDseaa1/iExP7TmCVZ+5UwhffajK1Xe0A28slS/yumREj0uoZM8hqEePScGb5rdnZg9jbWwsFgUsAcUFhYWiwKWxs4Qbk8GlTU5iheu1wbynihvm6XTRgDI+EV8AEl1HHEEIu3sv99Tq4/Dbt58upA+H2ZD744JbfT9UhvHEXAf1ydvVceYVgxsEQ40yxKq3KF+PrWc6NAhBl1rheG7I0QiWpn7+ETzx5doZ4+1+5jCOCldQrwOCYuHFTu1UfmJU2wxEBjS04oEwzv8wiq+Z3WtKld6RtA2R+C4C05aAX3KnNqhT9DpIL9ztkj3BwlrkHSReE9tOIN0iOdLtFbPnUSFOAVt4fGrbNJhIdMZUUdUn1x7xKl5ptPhMFaE0Qx2cT/6xvTc7H4nU9dUsSMMYiW3q+du/r70kBY9pMWjwy16zEqPX3+nsXO2tySirxNRPxEdEd/9HREdJ6JDRPQIEZXnv28hohgRvZr//O+5apeFhcXVhzF00edq45KLHRFVTuFzOT8w/wLgXsd3TwDYZIzZAuAkgE+JvDPGmK35z0en+0IWFhbzD681lCIAENFSInqaiFqJ6CgR/dFM2nA5Gtud/1yuVW4AzZNlGGOeJaIWx3ePi8uXAGgH/BYWFtclpimzSwP4Y2PMfiIqAbCPiJ4wxhybTmWXW+xajTHbLnczER2YzkPz+BCAb4vr5fn6wgA+a4z55SWe+QCABwDAU1aB2Ms5VfUih4a5VC3ou0m/Ztlp4dxQaGssWz6gyoXjXMmtDW2XfJGOy2xwPV6WwcSrtXqClNMlhWcWr0fLYKR3FFeFQ73kGAucYg6PK5nVHJTFdZTriNU4Yux6uB0VJ7W80JUSaiMiDutARAvVbt7EqhCjq7UQrP0pjp4T7BHyqlFtyeFK8rgYLZpEfDMPsO841+/Zr3UykjdwxKCbmztU3ivBFq7jHFvRFOlhR/FzLI+lW1epvGSJmEtCxjsyolWC3D08dzJBhyXHOPeBceux8A1zH/uFnC4V0uXG1ohAOo6go/dsOlpIP3GMveSktXgQNYdYVnv+Pl1/zUsO7zczRHYap7H5GNM9+fQ4EbUCaEIu5OprxuUWu9ddJu+1lLkIRPQZ5Fbtb+a/6gHQbIwZIqLtAP6DiDYaY8LOe40xuwHsBoBA49J5cMZjYWFxJcz0DzXPErcBeHm6dVxyuTXGxPMPuSW/hbzw0BIiulmWeS0govsB3AfgPxmT29waYxLGmKF8eh+AMwDWXLoWCwuLhYTJDiiI6AEi2is+D0x2LxEVA/g+gE9MtgGaKqaievIVADeK68gk300JRHQvgD8HsMsYExXf1wAYNsZkiGgFgNUAzl6pPm/EoP6VHK1rf5NWp6Ampj2ZAe2ccXQNb9lrDjIHON+iDfUpwHn+Jq1h/vJASyHd3c5qL2VHHO3YxTEX7rrpqMpzCav7s+Ncx/mDjaochDqC22G5EOriOspP62enA1x2ZB2XKxrQlGX8DUz9UiWajsVrBAUTt3ndmpq9fFjQPYcKjLtM0FNhdB9foim5K8XtdcZ89Qm6ntzAY5tNO5ylhvj390CnjqsgYdax5cJorZ4f6Q+uE21ytEM4OWh5mL/vvFOrl/jX8N9kdbEOxjt6ksc30qT3PHJswsJJhV/7KEXtHtGm9/WpvLEU95UR8Utia7WIoqOG50v5Yb3vCetQKjPHJFs7ydIuBSLyIrfQfdMY8/Dlyl4JU1ns6MIOLN/ALJEz3O4kNxE9COBOANVE1AngL5A7ffUDeIKIAOCl/MnrHQD+GxGlAWQAfNQYMzxpxRYWFgsO01E1odwi8TXkzg++NNM2TGWxO0tEf4jcbg4APoYp7LqMMe+d5OuvXaLs95FbvS0sLK5DZLPT0qu7DcD7ABwmolfz333aGPPYdCqbymL3UQD/H4DP5q9/jvxpqIWFhcWUMI2dnTHmOVxe9e014YqLnTGmH8B7ZuuBs4mshxCryr2C0+mk6RJyi1ItAPKLWKCxDlYbKT2mj9tjddw9D6e0iLKyjmUy9UuZcYertPwnkWC5yEhSq2S8teZwIS1ldk75RnUDy/2i56pVXkyIGb3t+saMEO+lKkSc1GVaEJUd5HZNtGg9hpo9LMtJFgt50hL9nqX1bLZV9F2tUjJ0A6eLermNgQEtfxzdxONUfkRPzYkbWBZXIuRy46MOcyuh4uA+rOWPAWFZFq8WnkG26WBN0Q7u48oTeu5kfELmWMHzxTem/yaXlLP52HBMtzHawM9Ol+n+HhddJ/sgUamKKdO3/l5totiR4rlUVM59FXc4XC3q5zoS2toSgUHMKuaDbewVlV+IaAUR/ZCIBvLmX4/mDxEsLCwspgYzyecqYyqaft8C8B0ADQAaAXwXwINz2SgLC4vrC/PBNnaqp7H/Jq7/nYg+PlcNei1wJ7IoPZM71qes1uiXlGtsjX7N2AhT11QDl5NqHAAQrpYBARxxCoTqRTjKlO79a15R5faMsvWAM47Fk26O4xlNCbURx0QYbOf2+kK6je4bmZsNF2uzg5I2QRl7BSUq0tTJVckqCZkxTS37b+Wy3mGmbSGvpnfhc8KKpEb/hjb+kstONIjYICn9LqXHhUWJw4JCxphd2cIc61hSOzod7+QbA47dQ90eVgEZ3MLUcrhFzx2fGApnG+WOZGyF8GyyRqt1nHuBrSgzK7R5T6bSYfIgIB2OxkRMlHTaIWIZFeMUc3iZSXD/x9xMXYtOaxobGOKXqTym2+g921tIH8YsYB7Q2Kksdk8T0ScBPIRck38TwI+JqBIArIqIhYXFlWCmdxo7q5jKYveb+f9/1/H9h5Bb/Kz8zsLC4gpYAIudMWb5lcpcK6SKXei+I2fJFm10hIkLMXXylOjTx+IXmMKMbuByDc9repER8Snid2knkSvK+AQvVMWWAMcjmladH+NjtNubtHrie6teLKTvb/swP7dcU8Rta9sK6VcP6d+WzHk+caxs01yh7AyfxI2uE4bvJzSduRDHAwD+4Pafq7xdoeOF9KfPvrOQPnlKW3kExcmek4L6Rrn/o9uEQ0rdpSg7x+89uFlPzbplTCAO7mXHlZ6mqCqHUn5WzOuYE36mgsly/uPLJh3OXQX/DS/T7SgSMSikQ01vj6b/RlSZTWla7x3hzGC3wymDoNAbtncV0ge7m3Q5EXvDFdV1hDr5eakQj7VXG3KgaIjfxRXVfyNd7xbz7B8wc8xnGktE77xUHgDM1HTDwsJiEWE+L3YAvgfg1fwH0PtQA8AudhYWFlPDPHDLfrnF7teRk9dtAfAogAeNMacvU97CwsJiUswHpeJLLnbGmEcAPEJEIQBvB/BFIqoC8BljzC+uVgMvC2LZSGBQy0VcIu6oO6a9gchgLsVtQt2hTP/6yBih5eXas0w4xTKwnqgIzOPScqJlZSxr+smJDSrvJ+DrynIWqLxpw3FV7t6yQ4X0waplKu/5UZZfvRJYrfKMi9voDfO7VZzSsslkOffP4/3rVV5dE1tvnOljywLPqMO5o5jMy7+n1e+jy9ksIFMk5GEtWk4UreN2ZBq09zBPkvNk8J1l2/SzWk+zbOu2zadU3guJtVxfDfc3OWICp5tYjWTU4TJxeILb4RV6CNIpKQDEhaGLyyE7TImAOLRBCy4j42zNcnsFt/+VU1p0HujnNre8oU3ldXW18LPFUDf8YkyVS1Xy/Oi41+FI1eF1ZsaYB6exU1EqjgMYQ86DcAhA4PLFLSwsLDTIXPy52rjcAcVdAN4LYCdyxv//YIzZe7UaZmFhcR1hPtNYAE8COATgOeR80L2fiN5/IdMY84dz3LYrwgSzyOTjhqYcVEQe97v6tarFhp3nCumj+1sK6cxOfTYf72O+G5rQWvb9WVb5kM9OD2ljf281a6aXlmgt9ZCfVT66z7B1xdkKbez/qT4+GP9Iy3MqL5oWKg8OFZu4oCmhbp5tUYeFQybENKv9aU2T/36Yr931XEfS4cw0U8RTqf9WbVUuVRxKhdR3qF7/BVAL97/vhDbiTxl+FxLNbz3rcHQqVDKODmg1IOPl58mYCCbtIDhJIQIZ0/OKPFxHaCN71IyQNsb3TIhYG1H9Z7Z5HcfGKPZoywtPLfPOl8eYutbV6bi0w0U8H2+s0LE2Wjdzn7h8wjmtT1NVqR5jPHoslr+O4wIf+wJmjnl+QHFBadjCwsJiZpgHK8nlDij+5Sq2w8LC4npG9spF5hqXC5L9uSvdPJUyFhYWFjB08ecq43I09neI6HKRfAg5p56fmzST6OvIRRHrN8Zsyn9XiVys2BYAbQB+wxgzks/7FIAPIxeD4g+NMT+7UuNdriyCgZzM4zfW6RC2B8JLC+nDxQ0q71gXy3JkzM3iIi0/MbU8IBMDWmZHMRZ4SLWX9CqtMpGKC9UWj1b5GB1kUzJXgp/VMa7j0MZF7NZzCe055fAZDiqzeWWnymvtZplPrI6/r9Rxf1DWyu/ijWi+UdLBcsWMn2WfWZ9W56kQkTydMU6946zHIGPUevq1iZVnnOt3qj4kyrldPtFXgfO6Dul0Mt6rPV4WiTrjQhXCHdW/+aXCqs+4HColVeL6Ca6/NKb7beROIZ+N6L460VtbSGc6tGNP7zKWW35o/QuF9OvLtYrrFw/dXUh/69BNuo4gj1kqzH2aKnOYz4m/fhnABwBOVGrztJniWpy+OnE51ZOvAii5zKc4X+ZS+BcA9zq++ySAJ40xq5E7APkkABDRBuQWzo35e/6JiGY3Sq+FhcW1wzxw3nk5md1fzqRiY8yz+cC2Em9HLuIYAHwDwDPIhVZ8O4CHjDEJAOeI6DRyKi8vwsLCwmIWMBUXT7OJOmNMDwAYY3qI6MJ+vgnAS6JcZ/67i5APpPsAALgrKjB2PKfm8M9H71blpJpB0QrNxitK2VPGoKCZg2c07ak4yhvfzJ0O7xrDTD/KTzE9mIhoneuYUK+IhTWdcVWIuKkiOmV3m1Y9gdDA/9fO23SeYFWHz+ouc4s+CHXwu4yt0j+rARGCYWCjptqD24TlQr9o0kqt+T/sZlWRknOa+iUq+N1G1vCGnUQsEABI9HKf1r2ksjC0ieuUXloqXtF96hHsMZHR7UiV8ntLJ5kuhy/NsdVi7vTrOgLDMs3jQg4BvO8MqyBl3bq/UyIWSUmnJlfjAe6DX9RwnPiRuFZpSgmrII9fc/7iIItjxrq5PqliBACeMI/F8Ebd/tneec13Gns1MZm0ctLuMcbsNsbsMMbscIdCkxWxsLCYb8jSxZ8pgIjuJaITRHQ670R42rjai10fETUAQP7/C3uFTgBLRbklALqvctssLCzmCtOQ2eXl9v8LwFsAbADw3rx8f1q4Io0loi8A+H8BxAD8FMANAD5hjPn3aTzvBwDuB/A3+f8fFd9/i4i+hFxQn9UAXpm0BgFXEgh15X4hUiU6LyPYZPYVfbo5UsXbeX+E13uPNnDA8M1Ml9xZx++C0DhPitNHtz6MVdfZxqTKywpa6xUnjG95/UFV7sgInyZ379Mny+lGcYJ8ueN8keV0GBleJZw4xvV7eiLSgYCg61HtoVPZujuaMbCd6wwM8PexMUdoP2HgH6mDzhMUOjDEJ7BJbRSAUB9zUs8p3RBPjBs5tJGnflG//suTJ7DjLQ4KKuaZpLSjqx39JqQeXodOw9hyngemTdPTDZvYciGc4Ek8MKIneHMjyx7On61VeePtTF1NA88PGtOUX87NjB4KBPpmV8I1TRq7E8BpY8xZACCih5CT7x+77F2XwFR2dvcYY8LIqZF0AlgD4E+vdBMRPYjcAcNaIuokog8jt8i9iYhOAXhT/hrGmKPIRTA7htyC+vvGmEtHJbGwsFhYmGRnR0QPENFe8XnAcVcTAGkLd0lZ/lQwleX7ws/BryDn026Y6Mp82xjz3ktkvfES5T8P4PNTaI+FhcVCwyQ7O2PMbgC7L3PXlGX5U8FUFrsfEtFx5Gjsx4ioBjm3TxYWFhZTwjRp7KzK8qcScOeTRPS3AMLGmAwRRZDjzdcc7iRQej7Hdntfpxl5zT7u3eJ2rTYSXsFykv7XCbY8pPWY6+rZ00T/oENGtZTr9O1lGUnPLVoNwCXU9ouCWma3aimP24mn2Annj/fdoMotW8ECK+8aLQDKJlgOs6Rae8YoXc6/SaeHOYBKRafDkmML91XNiw7vMULMM9HEfVzcpdUYfGP8nqkSXYd3Qo6NCFKzX0+/rLjNE9d/HaOshQFfWKqhaKuXaBPXWXzeEehGeCIpPSfktmHdH0YxF/0uvnFulxHVu7QTGCSFmkukSb+Lp4NlcZEluh8jKZZH9rzKlj4l51QxtG8SHm06dRsTFaKNwoNLYFCXM0IlpvyErr/qADv6PIlZwPScd+4BsJqIlgPoQs7w4Lem24TL+bN7gzHmKRl4x0FfbQwKCwuLKWE6OztjTJqIPg7gZ8j96nw9L9+fFi63s9sF4CkAb5usHbCLnYWFxVQxTUmbMeYxAI/NRhMuZy72F/n/PzgbD5oLkDFwpXK92PS0po/JMt6yx+v0ubonISwL2oUR/LgekcFWtmQwZbr+lOi6vlv4+5aWflWuo5+tMmJtWn3grIhXkfHzswM9eljOG1Yt8DicSWZKmYKVN2ndmTPD7EQztoIpNN2oy62tYPo7+txSlVf1KtOZoa2swjO2QlPEUBfzXd+EpmaDb2CqGTjB9CveoOmjW9BdpyMA6dgzFhcqOw7rgewQ56W0/080PsMigFQZt2O82eGQQFDoVInDGuReriMW4fuyDuexbr94t7SuI+0W1g9h3Y/9T/FhY0UPt8NJ66XoPrFFi2moncU01S+IuB4+pxqNM2AgI1mjVWJmigVhQUFE/0ZEZeJ6GRE9ObfNsrCwuK4wDxwBTEXP7jkALxPRrxDRRwA8AeB/zG2zLCwsridQ9uLP1cZUTmP/DxEdBfA0gEEA24wxvXPeMgsLC4tZxFTMxd4H4L8AeD9yAbMfI6IPGmMOXv7OuQdlDHzjuTP/kdXa24iU+WS8l3aNV3aGZStjyx0qE83slSMz4JBh+PinSVpptTscRlIvt6vsrMpCJMoyMJfw0JEJ6D2+b1jKFbX8xyxlnYcdFedVXvsY1y+djWa6tBrNiUbhhaNWb/bL/DxF0qILkuUOdYoot6vyoI5PuvS7/Lz2XxOywwlHvNZyHgtXSOtyuNrY6UOFUMNIFWt5rG+M2+U0gXKPsNxvYDu3aUyH24VLyNhcWlsIyS4WBO7Yxg4197SuUOVc53nczVKtllpykt87ssQpR+PrcS+3Y8kvdB1BIZeLZh0vWsxzc2Q9j2eoU88daS4Wrdd5Nc9rNaYZYx7I7KaiVPzrAF5vjOkH8CARPYKcL7qtc9oyCwuL6wbz4YBiKjT2HY7rV4ho59w1ycLC4rrDQljsiCiAXGyIjQAkV/zQXDVqqqBECr6TPQCA4mKtMuENC43+Uv2akVq+jgralnaoKiT62DJiyVM6r2sXqx2YANOvspc13U3cyaoKnqNa9aR2v9DiH2HaFq/S3im63sTlUiWOMyVBf795QsciiI/wcLlSXC5Rp9U1PEJdY3yVQ8WmhNtc3M4ztvFZh6pPKVOziZXaFYkUD5RWcqxVOqBjrUqLhHitprjBLmE1ITyiXET5BftKO9wddr6N46lKtRT/sC7nFeGDx9bp9/QP8NzZc5itXioP6PZWnGR1m8HNek5I6pqp1xYg3vNMSeMtzKHbHU5V0z3CU01CU9DgKhYjTHQwXc/6HeoxQgOpuMPh+SWg5+CMMQ8Wu6mcxv4bgHoAbwbwC+Ts08Yve4eFhYWFwHw4jZ3KYrfKGPNfAESMMd8A8FYAm+e2WRYWFtcTyFz8udqYygHFBX41SkSbAPQiFwrxmsP4vUiuzBlLB9u0gXxsKW/fkw7D9KpXuWzHvUy5Wh7WfCbewFwnVqO7quE5Hq3RVbzlv+h0rZ+pcOIGnVezj+nHoDDsTpXqn71AJR+bpc5rbpYVDhnplD6RLhVdUtbGdKzjzZr2yIlX87LTEYCwNunjOvzDmn4F+rnNcYf2feMv+FR7rJf7e2SdIxbGALcrHdJ5iUrOizeLI1JHfIe4CB0YadEUNNTGY+gWVdQc1Eeu402CwhVp+phayeKG5oeEJURU1xFeJsNC6jbKE/Vgjx6zpDgoD55iUYn/dZpMhTMiBsoJPZ7DXu7jol7ezyRLHFRVKA4EB3Sea1CfqM8Y84DGTmWx201EFQA+i5xH4WLkVFEsLCwspoaFsNgZY/45n3wWwIrLlbWwsLCYDPNB9eQ1Bdwhoh/NVUMsLCyuY8wD29jXGlVj2v7fL4CI1gL4tvhqBYD/CqAcwEcAXAjJ8um8e5dLIlFJOP3bOblGUAQZAYDy0yxDKmnTXiEyxSwL8Qgt8oGbtSpE5VG+L9yij+LHm4XKSpBHzjeq5Serbu4spPsmtG5LX4AtHAJ9XF/NxgFdrp9lMKZUy5CqX+IhHL5LezPx7mXZ2eAWMdQ+h1lAC3dC0fNahhRu5vtiVZwO9OtnTSxjWaLzpC26hPOybhGcKObwKCJkSLI/AD1O3qM8fuPLdX9EREzZ4hPam8nEGs4LtHPeud/Q7aUEv0DRaW2d4N/Jct1wC/dvw+NDuhwPO8Zu0HGAYzUsFw1v1JYi77xxXyH9gxNb+J5hPXdKhWPS4k5dR6SB5+qye9oK6bYhbd2DAywg9I1r+ebgG5fxxTcwY1yL01cnXutid2CmDzTGnEDe+iIfKq0LwCMAPgjgy8aYv5/pMywsLOYZFgKNJaKP5w8oYIyZbUXiNwI4Y4w5f8WSFhYWCxYLRfWkHsAeItoP4OsAfmaMma2mvgfAg+L640T0fgB7AfyxMWbEeUM+3NoDAOAtrkDl/hwlSDvs9NMBpkj9N2oKIB0heiLC8LpZ1xFeydR42U2dKq/zl2yxUS6c9I+u1V3T2rqEL7w6r/Q4d7+kcIMHdRxQlDAHWLOuS2W19XE73F2agiaqhKa+zHLEA0j3c+fFqvTvn1Slia9iFZLwcu1MwC8sF1IOywWf0JpwC8epsfXauN0l1EiyQ5qCGuF4MtAr4ky0OeMq8HVaSzaUx4bEChFP1fGXF6pi8UVqSFuDRI+wqMMn6p9YX6XKDW7iNlYe11Q7U8TP85fpPmgN12MyeLt0f0hHs4lyRwyKan7eexo4/PLnzuvQMf6tYmD267njjcwy71wIOztjzGeRC1r9NQAfAHCKiP6aiFZe9sYrgIh8AH4VwHfzX30FwErkKG4PgC9eoj27jTE7jDE7PIHQZEUsLCzmG+bBAcWUTmPzO7ne/CcNoALA94joCzN49lsA7DfG9OWf0WeMyRhjsgC+ilw0cAsLi+sA84HGTkVm94dEtA/AFwA8D2CzMeb3AGxHzv3TdPFeCApLRA0i79cAHJlB3RYWFvMIs73YEdHfEdFxIjpERI8QUfmV7pmKzK4awDudhwjGmCwR3TfNhgYBvAnA74qvv0BEW5Hb4LY58iaFcQHpYE4OU3FSH78H29jcJbxeq5RkfCy76d0lZBOOEdi4luV0x45qgZ5p4Od5x/mov8jhw7nslHDU2OAw0xKPDgjNhYhDwWflgyxfanuz9u4S3MJizVRGy24Cj7FcTXoUGbnNob7yPE+DZJluY1Efp32bhSePJVrGkw5yHZliLe9Jj4pAOmtZ7ndjfY9ur5vVH16iFpXnEsGJXOdZBju+WqtMFHWJdhTp8SQv10GjImhPWL+zifC7VZ/RfTW0SfSxqH5gm/5TKrmJ1Yd6y7XqSTYg3iWtx+zsAMv+0lGus+aUKoakCJYzeLNuoyvO/d3k5fnx+rW6kpfalhfSw+u1alXlcUcg3Jli9ndyTwD4VD7c4t8C+BSAP7/cDVOxoPivl8lrfc1NzN0XBVDl+O5906nLwsJiAWCWFztjzOPi8iUA77rSPa/JgsLCwsJiOpiMxhLRA0S0V3wemGb1HwLwkysVeq1KxfMKxg0k85oBqeJLx5no367X9FSD0LIvZ0uA+CmtZtA2zPogt247ofIOf3tDId34NGvV9+zSWuoyZirGNVWIyVgH1VwuE9XD0nWn1KHQP5HuHzJFH92pKZ2pEZ5C6gSFc6h1yBikThWejLiuL2KVjE1bNAV9+Yfs9cs7rsdCxYB1czuG4vo0vWsPO9f0RTS1DL1usJCOuJjGvmXHIVXuycFt3HanY0/hGHPJM6zyMbBVU/KaA5znSmqK2L1LeKcp5jYax1+SS4hEjEPlSA5h+dO6w2UsiKpeLjhwqyOQrlQfclTf9DT38cfGef2QcU4AwKzh8ZxYpkUPkaViDH+MmWOSnZ0xZjeA3Ze6hYh+jpzqmxOfMcY8mi/zGeQOTb95pSYs6MXOwsJiYWA65mLGmLsvWyfR/QDuA/DGqej+2sXOwsJizjHbqiZEdC9yBxK78mcAV8SCXuxMwCCxMkc5wlFNRbrfyJR01Spt/XD6PAcxqAgyje0xmsZGBpk+vjCuvVsFBfs481tMJTM+/RNmBCX1RDSdzhRx2R0tfNjdMa5P0bsNn+XUPaspYiooaM8rzuHkGUZZfnZkuaZEiXK+z//6QVwKJ08xzeyq1X1Ve4BP72Q8CgBwJ4WjhL08Tr036zroBj6p9b6orV7CB8R51hYu92J3iyrnFfQ3Va7HwiUOGNs+wnnZEU1VvcIB6JBDNFBez6f8SXGSGh3W5hqNxVyu5SbtFPbAM2u5/hv1s/39PBZlZ1jcUrNHO1448WGODVJxRM8rynKdWSGxWPkNPba9u7hPY7Wa4jpPsmeM2T+N/Z8A/ACeICIAeMkY89HL3bCgFzsLC4sFgtk/jV31Wu+xi52FhcWcYz4477SLnYWFxZyDstd+tVvYi12S4O7OyVfSO3VAklWV7Ibj9Gl9el21h197rISt1Iq02ALUI1RFSKuNhIXmvjvCspvqA7oST5yvy/dojyVdv8oeUfZP1zXoogAAIABJREFUsBwnVe1QMxDeQOKVWj4Tq+G8YI9+tk94xshIH5Quh2WBEG2NhbU6SCjEahjlh7nf0nfrOrp+mwVi65s6VF5rF/d/YD/LtjwRVQyxONfvuk2Pp0vE3HV7WSY1flJbx5SOcbuSpbqv4muEk1K/iNObcHqLEWNWr9sR9LMcbXSYn01J/ayDe9lPhqdJy8+l01LfsJZvJiv43fq386AFHOpCJIL4DN+g5X4yOJRHiPooouV+Uk7nlNH5hx1/DDPFtV/rFvhiZ2FhsSBgaayFhcXigF3sZghi7fRETNPMrHDUSM7YnyVctriLOZx3wkEHRMwF34RWY5hYKmhLo9DG92u6sfK7It6DS1OdrGhysXCzkAjrd4kuZ8oV6nE4gvQxDYpre3OM3sB02DsiHF6e1vWLroKrXVM6M8g6Nm7h9HR9rbagkBYDzx/XB2UkaHNasOTgDq0KUfcgW5/03qXpHVXxexuhEvT+e55V5R788R2FdFGfpmLjjeJ6L6u9uB3xVOPC6SnF9HhKGnvb2jOF9L4u7aABh5l2Z5JaNFB2lufSwH3aeafvlFBhEc2V1hoAQLXCyqNfx8lIFXP7pfVKpkar+lQcF3M/quf32PLZXRrszs7CwmJxwC52FhYWiwELMbqYhYWFxWuGpbEzhD+YxMqtOVOwtNHysCKPkHOVaLlIvJJlHI1PsIrK8Dank09Oj67WMiR/c7iQ9j7HTjKTOg4Nzr9FeMmo1IF0WlayKko8zUMxPlKiyklvFfH3a30N8zwL6jxaswDeIa7TO851VB7Tqi0yBuzIBm2m5UrxT7JwKILRpPbWsaa0v5AOntZyrpSQiclhiuzXQsaJ9SIQ0rCempkG9goj5YPfP7NVP6uK3y1V4ZBzRYTqjJBrYblWDUmP8vwo3qf7YyLN1/tDrFLjdVhnyp1MzOFJxi/M00x/wJHH6eIuLjeyzhFYSMjpQh167kuHpu4o5517u55XaTEuxed1f4+v1LLhGWPWYnRNHwt6sbOwsFgYsDs7CwuLxYHFutgRURuAcQAZAGljzA4iqgTwbQAtyMWg+I3J4sZKJNMenBvIqSu4HFYBSaGK0lA3qvIGXXwE33sHqzuMrdNb92JBVSf6HLFnTzMlCE7ws6ONqhjSpVwnBTV9XFPG1K8nxm26q17HCvjOse2F9JYarfLxfDl7rsg6/JcGBoWmvrAsoIzuq4kWVo2oPKItBjrvZl6eFQ5GVTxcAOdquB3xZkdciE4R42K94No9WmVi2dbuQvrmqjaV993WGwtpj7CgyGQ0hVvawuosnSe02IAysj84He3RPNOb4LzqQ0lcCiNrma6Tw+hFqY1U6HnV+QZhKaIlLHCJx42s5QHNbtbj4j/I8y9R6Yh7KyipzHNaSTQ+y9cjq3U73JHZdWI+Hw4orqVb9ruMMVuNMTvy158E8KQxZjWAJ/PXFhYW1wEoe/HnamM+xaB4O4Bv5NPfAPCOa9gWCwuL2YQxF3+uMq7VYmcAPE5E+0SQjTpjTA8A5P+vnexGGaQjE45MVsTCwmKeYT4Eyb5WBxS3GWO6iagWOU+jx6d6owzS4W9ZYlKRnNzENa5fxS9ilfaM1ehKSngPnSliuYgnqtf+yDirBRRVab2OysdZzuWNcH2pkDbF8oq4sWN3afnPoSEW8Ek53U/aN6hy1MXteGZsvcqraON0wqFqkSznGZUSKjGRJt1XpWdxSRR38LslS1lGlXbEhm2uZPHq3etfUHlfeZpDCXjP8ruEtBMYnPOxB5qOGq0GJGOo+oVJVfGuflUuI2zfdu08pvL2dHPs36I1PBYTndozdLaM+21oo1ajiTYJD8dCdmiKHfLeY3yflFkCQHy1MPXy6n5MrhbeWCZYpnnfCv0n8lh6YyHt9E+SOssyyJJz4vuQnt+j7JgFiQq9+rz17j2F9D9iFjAPDiiuyc7OGNOd/78fwCMAdgLoI6IGAMj/33/pGiwsLBYS5sPO7qovdkQUIqKSC2kA9wA4AuAHAO7PF7sfwKNXu20WFhZzA8qaiz5XG9eCxtYBeCQfJMMD4FvGmJ8S0R4A3yGiDwNoB/DuK1VEaYKvN0cb5XYdAGKCubpSjmAi5UJPwFw63mzFM0y5hrdpOhOv4N+J4vMsO6wfdugSyDZVaw32vkamKT9N8VDsbDivyj3dypSu6Un9LoFhVgdpf5NW5XALx6HxBn7nmqVao2ewnNVvovfpSVhfxuogAz2sXlJ0TKtrnDKsijKR0n3lqeY+8bUz/ZeBeIDceF7ApqZulRdJ8budHWQPI6l9WkSRrOX3vKdBU7+xah7PThnUyBFPNSCctrodmifeMI+7EU5Vk6Va9ySynGltqFGrjRS9ws92nkomBZ10i63I21+/X5XrXMJ1HOrQakDJRqbCxsXvsvI7WgUrXsfigLEVWvzy2M9uElcPYcaYBzT2qi92xpizAG6Y5PshAG+82u2xsLCYe8wHC4r5pHpiYWFxvSJrLv7MAojoT4jIEFH1lcoubHMxg8L2mBx2y1WtTCs636Ez3eIEbHw9b/k9Qw7jc79wAFqm+Uyol7f9yXKmWP5BfWo7cBNbRjg12Jsf43ZUf3qokH7mnFZnz3rF6eAmtyOPf6/cSUf8Cxl/QFD5soCm2itu5JPgvS+tUXlb72wtpEeiTF3f9p5XVLlvPnE7l3tOx/xILeU+3nAfO7w8NajnJ7UxzT+4f6XK8whHBhkRDzYd0DywcSn349qAtjZ5doCdit7TxBT3odbbVDkXSwaQ0n43lTigqIPni0npfYOMETzRryuhtfwAf5sWPfjWsNVOVJzGfnz/b6lylSUsOsmOaLGBW1iASBGOnIsAkCzlPLdD+uLVzHvmmIOdHREtBfAm5MReV4Td2VlYWMw55ug09ssA/gxTXErtYmdhYTHnmOw0VhoI5D8PXLmmfH1EvwqgyxhzcKr3LGwaa2FhsTAwyd5LGghMBiL6OYD6SbI+A+DTyKmtTRkLerEzHiBVkZPZjK7Vm9QxqVISceiYC0sJtwiWk3aoIISbWNZEvdrJogxqEq3lZ4WXa/lJOsB11u7X6gkDW1nud/748kI6eFarAaSahOeUzWGV5zrMphFVR7RscmArv2ewi9vY1auDw5ytZrlXYFD3Y1+C6y8WwWZax/UczFTyuwVbHVYk43zdWl5XSBc9rz3JLP3lWCE9tE17QR3ZyP39W7ezhcaxsG6HlPX9dfxeledx8Xs+cmZLIe1rmVDloiFWyfCOahlpaatwsnojz52SQ3p+yEBO4826P6SHlOQtWjgWj3PZ8hdYZje6SdexbinrWvUmtWWlVI8JiJhG482qGDIB7lP/Cj2vqoqFN9K/xoxB07CFNcbcPdn3RLQZwHIAB/MqbEsA7CeincaY3kvVt6AXOwsLiwWCWfRyYow5DGE7n3cZt8MYM3jJm2AXOwsLi6uA6ezsZhsLerHz+NOoW5FbzPtOaTUGadTvHnfEp1jHmuQBL3OKCb9WA0h0Ms1yatLHKrnO8K2s4+Fp03RGqo2E12ma7BZOP0NVTBsiHgclKuf6I13aCiMoVAZSQf2eSRFrFS6mY8bB6iuO8RfJt2rrivZxtt6QsXiTGT11yl5l+u6N6Z9xI4q6D3OfxqtUMbS9Q8RydahCuJewqkVXnK0Huie0OsX6rWx9kspoCnrB0SsAlP6c1UEkRQYAFHH7/SOOMUsIC4c+ni/SOSoAJEtEX23XNDk5zOPrPaupPJp5rEc3cDuMwzntL86welI2qPs7JS+J54Qrrd8l5ef5oWgrgIEXGzCrmMO1zhjTMpVyC3qxs7CwWBi4FrawTtjFzsLCYu5haayFhcViwHyIQXHdLHaeGi3kKXuc1QeiDVpWkUjwMb486k/H9fF+RSvfN7ZKZaFomEdvbJjlVZ6oftbmt7JZks+lVU+ODrLaxNg4txdJh+wtycNEKafXEyETXKnzSk+KOKnCSUmiSs+80TvZfImiDrml6KtkhN/z5hvaVLmDm1ivwdRqGVV4gOVSFBdqEf2OWKii+c42ZkQs12d72HGl29Hf/U2ssrKkXssfPR6uU8rpKo7oOsbWcLuSJU5ZHKcDq4WqTKmWvZWc4r5/yyrtRPTRPRw8KFWp50RlCcvOmppYi+J3Gp9V5f7qxH2F9GBMq+mUrWWZ9NJSbuPhfctVuZYf8LM77tYqPE23steZk5gF2J2dhYXFosC1X+vsYmdhYTH3oOy157ELerFLp93oG8ipHmxt6VB5B27mLbt3WNOl9ABzOl8t0wZvp8PppPAaEhjSVGf0Xaz5XibUV8artdrI4V4+wv/IuudVnlcIMorq2VrjiZ/eqMplS4QzxoqUyhtdKzTuT6gsBATV7nqrsMKY0P1x39oj/Oxz6/SzTzA9CwpLlEfGd6pyFcuZMv71hv9QeR/76Qf42RWsw1O7VtPMrle5r0rX6LxbhEPTZ7/H/VN2l1aYD3h4LPp/pp1aJhq4P6heWM6cC6pyxcKHRlIzRNBNTAuri1kdZteSM6pc6GYWDZyLaB2bUBv/2UVW6PG8vZEDgvznmmcK6bMp3ZChITEuFdrTTriVn3ewglVzgsu0tUbHm7jOdK3WrWo/UYdZxbVf6xb2YmdhYbEwYJWKLSwsFgcW42KXd7j3r8h5M8gC2G2M+Qci+hyAjwAYyBf9tDHmsctWZgCTzp3uHWht0c8pYjqTdoSQk9Q1OSAoTJneaw/sEI/y67xSD9PCmDixzIQdofeI89ri2srDJZx6Pd/NtDtVrp9VEmRKlDitw/75xphaDm/RjgBIODYoOsPtiq/VJ9ePPc0vmq1LqDy3OKAuuZ0Dvu2s7FPldgkOfW9Q1yFjS7hEWMjBw42qXGY933drQ5vK+/HhzdwOwfy623Sf+gaYopcO6T+wrEecBJ9gUQaldTmfsGxJlum5kznItHB4jNM/XqNp35q1fJrZPuwIC7mF59/GRt2PtcJr5rEk09F/7rldlbsw7wEg6nQOWs/9+P4bXi6kU1ktvhhayvetDep27D6mHZrOGItxsQOQBvDHxpj9+Shj+4joiXzel40xf38N2mRhYTGXWIwyO2NMD4CefHqciFoBNF3tdlhYWFw9zIfT2GvqqZiIWgBsA3Bhr/1xIjpERF8noopL3mhhYbGwYMzFn6uMa3ZAQUTFAL4P4BPGmDARfQXAXyGnfvhXAL4I4EOT3PcAgAcAwF1VDsprxXv6tOa/K8nCJo8+mUfjd1h20beT0+PbHK42xHjIID0AsKOeVV2yYJnUs90bVTnfCP+ePHpgq8p721b2KN1Qys4TXS36WcMDrCLgKtWTRAZUKWvVMpmxjUKuuFRo6juclBb183XEp2WO3pUsQxrdw04iX1in1TX+n/onCul3nblP5ZWc4T6I1XL76RYdx7RKqPAsDQyrPM+AEB7ewWopdT6tupFdyu8y2Kw9xJTtYXlhKsTlnB5LPAnu//qXdf3t93D/ZEWT/IO6789X82+1kgsDWLO+s5D+25aHVd7bnv19bof3Fn6Ww1UNxfh5O7aeVnm3VrAazN6xlkL6xTPagiJUwvM926D3PW73LO/E5oHM7prs7IjIi9xC901jzMMAYIzpM8ZkjDFZAF8FsHOye40xu40xO4wxO9ylocmKWFhYzDdkJ/lcZVz1xY5yfpS/BqDVGPMl8b10oPVrAI4477WwsFiYIGMu+lxtXAsaexuA9wE4TESv5r/7NID3EtFW5MhjG4DfvVJFrgihZF+OmjgUzLHk56zdPrxB04ihjax2MH4DH9MXtWrrh2QZD0jWrwfnlSI2fE8kuBu37zilyh14geOw3r5Rm1RL1ZNiL7ejuEKrbhQJetcT1vEGAoMijsCobmPZalYniCQF/XLMs5ElzMcC/fr3zwxxx7oFowsVa9nAH595dyE9Fnf2I6c37mLKFfRoiugXjhL+z1Nv1I0U8RLGh3lHf+cWHVxqMMGWBT63VsXpaxbqIaIPxlfpcnX/ndVG0h2dKi+04dZCemybsDpwOGjwH+d21O4YUHnDMZ6PXx3SKiU0yOPU8BRvf7pvdzhNKOe+ur1C09gyN8/954+xB4tAmZ5XE2H+O3gFOkBF0K/HZsaYBzT2WpzGPgeAJsm6vE6dhYXFwkXm2p/GWgsKCwuLucdi3NlZWFgsQtjFboYg4EJ4WK926ICzv87yiKVPaPlD1y6hMzAhHGM6xmP96zg255mfrlB5rjY220qL+KFbS7WMZ08F3/fL46tVXsPPuB1rP3G0kC5y6/a2jXGgmKpXtQTACBFbtM4hb0vzu1UFWY5z6pCOG1vcwp48TIdWb4w2Mv0oX83qIF9Y931V7oM/+51COlQfUXkpoS7z6lmWDf3v1/+rKvffz/1KIe2O6/dMlwi5mnBu+mKvVqf4tWaW4UXTWo2mx8PtkPFxjccRK/de7h9PTHtOkQF3Gn/K/Tu0yaHOcwOrx4y9rOWsUojzaL3u73KhpjMu1GhCnbr+8SaW2TV6tYeYgTSr3IQqWbZKL+ngRCW7WJYYieu+ipzRZWeMOYhBQUR/AODjyFll/dgY82eXK7+wFzsLC4uFATO7MjsiugvA2wFsMcYkiKj2SvfYxc7CwmLuMfsHFL8H4G+MMQkAMMb0X6H8wl7sTCiLxM25eAfJCb0N/8Zd/1xIf3T0YyovVccqA0Vn9X0Sx59nihTQp/bIiJgOJGJQfPusdrzpDrPKQNav6VK0hq9f/jF79Si+RasqhPexZ4+ikKYzkSamB0XacQWGDtcU0n1CVaHihK5jJMhqEkVaawTNm3sK6Y5+ptMffFobt1Tv4fccu8cxrQSDqXye++pjvv+k6yjn2BXpGk3l657iOgfezINxZ6NW9XlphMfsyOFlKs8l/t5conpyWJSMrhMqRwFNvwK9gma2iAyHfkHsINPTzFodkzUd53fxBrXTzIZ3cVD7M/087tk2rUAfKmbRyaND2jJnbxeLCvxCbSlzq6a7iRS3I9alY2h4l+g2zxiTyOykNVQeu40xu6dY4xoAtxPR5wHEAfyJMWbP5W5Y0IudhYXFAsEki11+Ybvk4kZEP0fOFZwTn0Fu7aoAcAuAmwB8h4hWGHPpkxC72FlYWMw9pnEaa4y5+1J5RPR7AB7OL26vEFEWQDXYH+ZFWNCLnckSktHciWbxcU1HPzT8e4V0eY/u6OafMgXo38EnorF6x4nagLh2jpW4DrawEX8sptvRuIm55dAvGlSeNIpPVTLHivVrcxC/MAJ3RlYv2sDG9BU7tFXDRB9TqbVNLNI4PaFPY10iJkW8RstWpJVHdpCdLfjHnHSar2sr9NH44Cm2GKh7oquQHlmvnXeGj3G5wA26jokm7teWhiFuk8NA/sgRpq6+YS02oPVMk80gn1h6w7oOaY0TbNd1BPu4P2I1fF/W8ZekaLJXW2igk2Ug6Vo9nif3cvuzwmGsy63Ljfcx7Tzlr1F5LheXHRth+muyDmcCwtifKjWdTndpq6MZY/ZdPP0HgDcAeIaI1gDwARi83A0LerGzsLBYIJh9PbuvA/g6ER0BkARw/+UoLGAXOwsLi6uBWT6NNcYkAfz2a7nHLnYWFhZzDjPLenbTwf/f3rnH2FFXcfzz3e0+2u622xYppQh9EUx5BLCIAiFqREqN4jPwj0IEo4kNaiRaQyRoMBETXyjRKKIoCCgiNiQGUEuKVmppKS3l0QeUlrZsael2+9jubu8e/5jZ/uZ3u7e7221779x7PsnNnTvzm5lz5jd7ds6Z3++cfBu7PqH9iQp7Z8VDFSauCKqN6o6fbt98XwjK9GbyOx445WDUbsy2EMvaPzmOd3TPCnG/+kycrrArTiK6c3WImfS0xXJYJg6j7mwSzvgYvZlRB1Mei+vjvnxWGOG/d1I8bqRlRYgNrdsZ4nR942M9R+3IJDrdH+u5+a0Q98vOauiZFo/FsQPhWnUfjG+rbMGdAzNCfKmxM46Hdc3IDAlaGSfeLGQuyYSmMCxiUkM8W6MpU0ypu7Uh2nbeySFu2fKxcB3/+2o8C4Pt4WR9DfH12DEnE+caH+St2xpf+1GZGrv1y2NdxmYy1XTviPs6y7gwgYetc4vuzfXhntu7IS72844rQlx0/KQQxlr75MyoXdfp4W+m7ZTOaFvHztJyHRXHYQbFcMm3sXMcJx/43FjHcWqCCii4k2tjp15o3pq4Tz0TY3ej9Y3w2L9vcqxmXaZOaF1P2K91bdyuK1OStG193FmdvcFtKWQ8mMKkeJhB15TwW23x631l3KWmncGl6y3KNr9vetBl053xsJTmZRlXsH10tC07tGX0ttCuMCd2QQstQa6eg3GSyMunv3poefno4DKPKar9sKcr6DJzQjwCYOm0IPNb+0K7uqJZKXWNmWtXlEO7PnPplq+bdmh5VXucXKHxrOCOqS5+mni5PUyfbGwM1/STs1dG7R7eF+roNm6MBWnO9FPX5HC9x26Oz9WV8SxPfTqejdAxK+zX2xLvV8iMXNo7M9xzDe2xS96cqYnb0xbf+69vC/VmrTv054QdRecaHe73XaNiV/uMx0NfvM4xwJ/sHMepBaxQGLzRccaNneM4xx9/QeE4Tk3gQ08OR9Jc4KdAPXC3mX2/VNumjgJnPJZkcjgwJc7acHBMiK1Mei5+rb72+kx8Qpm4yO44PlMYE/4bdfaVLsTWfVoIKJ09c0u07cVVIQPF1JPjOqkdrZk6pp1hiEf3+Pi/YP24cPy9HXFcLiNiFMcB2H96cB2a3g5d3bsxvlaaGqaZNa6Lj/+fUWFYxqXTwliIZzZPi9od6AyxuKU74yEOjdvDuaNko2fGMcyGzeF6tBTFwPacHuJSo5pDvK2uJ56e1/NyiA8274ljWftnhfMVdoZrsLI1TtDZvC3ExxSP+GDcprBi9yWZjCKN8dCTpkxXHxwb/5l1Zi5P43nxPVEohAt00pgwvKm9qS1q1zEmyNjXHLuIDZl45Ki1Qa5xm+I4a8fZ4VzFMcFC0zEeBFwBT3ZlqRtbCkn1wF3AVcBskopjs8srleM4I8b6Dv+cYCrK2JEUxl5vZq+m00EeJMlG6jhOjrFC4bDPiUaDzJ09oUj6NDDXzG5Mf38WuNjM5mfaZBP+nUP1FtM+iUGyOOQU1yt/nGVmrYM3q2wqLWY3UD3ZyBpnE/5JetbM5gywT+6pVt1cr/wh6dlyy3AsqDQ39g0gm2ztNGBribaO4zhDptKM3TLgTEnTJTUC1wILyyyT4zhVQEW5sWZ2UNJ84HGSoSf3mNmaI+wy1OIceaRadXO98kdV6FZRLygcx3GOF5XmxjqO4xwX3Ng5jlMT5NbYSZor6RVJ6yUtKLc8I0HSRkmrJa3sf80vaaKkJyWtS78nDHacSkDSPZK2p4VQ+teV1EXSt9I+fEXSleWRenBK6HWbpC1pv62UNC+zLS96vVPSIkkvSVoj6Svp+tz32WGYWe4+JC8vNgAzSEqoPQ/MLrdcI9BnI3BS0bofAAvS5QXAHeWWc4i6XA5cCLwwmC4kUwKfB5qA6Wmf1pdbh2HodRtJJfritnnSawpwYbrcCqxN5c99nxV/8vpkVwvTyq4G7k2X7wU+XkZZhoyZLQbeLlpdSpergQfNrNvMXgPWk/RtxVFCr1LkSa9tZrYiXd4DvARMpQr6rJi8GrupQLbyzBvpurxiwBOSlqfT4QAmm9k2SG5I4OSSe1c+pXSphn6cL2lV6ub2u3q51EvSNOACYClV2Gd5NXaDTivLGZea2YUk2V6+LOnycgt0gsh7P/4CmAmcD2wDfpiuz51eklqAvwBfNbPOIzUdYF1F69ZPXo1dVU0rM7Ot6fd24K8kbkG7pCkA6ff20keoeErpkut+NLN2MytYUhT11wR3Lld6SWogMXT3m9kj6eqq67O8GruqmVYmaayk1v5l4MMkmVwWAtelza4D/lYeCY8JpXRZCFwrqUnSdOBM4H9lkO+o6DcGKZ8gZODJjV6SBPwGeMnMfpTZVH19Vu43JCN4izSP5M3RBuCWcsszAj1mkLzdeh5Y068LMAn4J7Au/Z5YblmHqM8DJC5dL8lTwA1H0gW4Je3DV4Cryi3/MPX6A7AaWEViBKbkUK/LSNzQVcDK9DOvGvqs+OPTxRzHqQny6sY6juMMCzd2juPUBG7sHMepCdzYOY5TE7ixcxynJnBj5xxTJH1J0ueG2PZ6SW9Jujv9/X5JJumGTJsL0nU3H6U8iyTtlVSVxXCcoePGzjmmmNkvzez3w9jlIUtLZ6asBq7J/L6WZAzi0crzAaAqqmM5I8ONXY0i6aJ0AntzOotjjaRzBmj3UUlLJT0n6R+SJqfr75R0a7p8paTFkurSHG83p+tvkvRiep4HhyjaJqBZ0uR0dP9c4O8ZeZ6S9BNJSyS9IOk96foWSb9N8wKukvSpkV0hp9qoqII7zonDzJZJWgjcDowG7jOzgQqO/xt4r5mZpBuBbwBfJ8lxtkzS08CdwDwz60vs0yEWANPNrFtS2zDEexj4DPAcsALoLto+1swuSRMm3ENSLP3bwG4zOxcgL8lOnROHG7va5rsk84wPADeVaHMa8FA6D7QReA3AzPZL+gKwGPiamW0YYN9VwP2SHgUeHYZcfwIeAt5FMk3rkqLtD6QyLJY0LjWkHyJxeUm37RrG+ZwawN3Y2mYi0EKSobYZQNL3+tOMp21+Bvw8fWL6Yn+7lHOBncCpJY7/EeAu4N3AcklD+udqZm+SzEG9gmRe5mFNBvitAdY7ziHc2NU2vyJx/+4H7gAws1vM7HwzOz9tMx7Yki73Z8FA0hkk7uwFwFWSLs4eWFId8E4zW0Ti+raRGNahcivwTTMrDLDtmvQcl5G4rruBJ4D5mfO7G+tEuBtbo6TDQw6a2R8l1QNLJH3QzP5V1PQ24M+StgDPANMzaYFuNrOt6VCR30m6KLNfPXCfpPEkT10/NrOOocpnZkuOsHmXpCXAOODz6bqyA5rWAAAAcUlEQVTbgbuUFMQpAN8BHimxv1ODeNYTp2xIuh6YY2bzB2ub2ecpEiM75OEkR7OPU324G+uUky4SF/ju43UCSYtIcgb2Hq9zOPnAn+wcx6kJ/MnOcZyawI2d4zg1gRs7x3FqAjd2juPUBG7sHMepCf4Pewjg4ARrZqUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotting.coeval_sliceplot(initial_conditions, \"lowres_vcb\");\n",
    "plotting.coeval_sliceplot(initial_conditions, \"lowres_density\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run a 'fiducial' model and see its lightcones\n",
    "\n",
    "Note that the reference model has \n",
    "\n",
    "    F_STAR7_MINI ~ F_STAR10\n",
    "    and\n",
    "    F_ESC7_MINI ~ 1%, as low, but conservative fiducial\n",
    "    Also we take L_X_MINI=L_X out of simplicity (and ignorance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-14T21:22:09.380149Z",
     "start_time": "2020-03-14T21:18:50.496755Z"
    }
   },
   "outputs": [],
   "source": [
    "# the lightcones we want to plot later together with their color maps and min/max\n",
    "lightcone_quantities = ('brightness_temp','Ts_box','xH_box',\"dNrec_box\",'z_re_box','Gamma12_box','J_21_LW_box',\"density\")\n",
    "cmaps = [EoR_colour,'Reds','magma','magma','magma','cubehelix','cubehelix','viridis']\n",
    "vmins = [-150, 1e1, 0, 0, 5, 0, 0, -1]\n",
    "vmaxs = [  30, 1e3, 1, 2, 9, 1,10,  1]\n",
    "\n",
    "\n",
    "astro_params_vcb = {\"ALPHA_ESC\": -0.3, \"F_ESC10\": -1.35, \n",
    "                    \"ALPHA_STAR\": 0.5, \"F_STAR10\": -1.25, \"t_STAR\" :0.5, \n",
    "                    \"F_STAR7_MINI\": -2.5, \"ALPHA_STAR_MINI\": 0, \"F_ESC7_MINI\" : -1.35, \n",
    "                    \"L_X\": 40.5, \"L_X_MINI\": 40.5, \"NU_X_THRESH\": 500.0, \n",
    "                    \"A_VCB\": 1.0, \"A_LW\": 2.0}\n",
    "\n",
    "\n",
    "\n",
    "astro_params_novcb=astro_params_vcb\n",
    "astro_params_novcb.update({'A_VCB': 0.0})\n",
    "#setting 'A_VCB': 0 sets to zero the effect of relative velocities (fiducial value is 1.0)\n",
    "#the parameter 'A_LW' (with fid value of 2.0) does the same for LW feedback.\n",
    "\n",
    "\n",
    "flag_options_fid = {\"INHOMO_RECO\":True, 'USE_MASS_DEPENDENT_ZETA':True, 'USE_TS_FLUCT':True, \n",
    "                    'USE_MINI_HALOS':True, 'FIX_VCB_AVG':False}\n",
    "\n",
    "\n",
    "flag_options_fid_vavg = flag_options_fid\n",
    "flag_options_fid_vavg.update({'FIX_VCB_AVG': True})\n",
    "#the flag FIX_VCB_AVG side-steps the relative-velocity ICs, and instead fixes all velocities to some average value. \n",
    "#It gets the background right but it's missing VAOs and 21cm power at large scales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ZMIN=5.\n",
    "\n",
    "lightcone_fid_vcb = p21c.run_lightcone(\n",
    "        redshift = ZMIN,\n",
    "        init_box = initial_conditions,\n",
    "        flag_options = flag_options_fid,\n",
    "        astro_params = astro_params_vcb,\n",
    "        lightcone_quantities=lightcone_quantities,\n",
    "        global_quantities=lightcone_quantities,\n",
    "        random_seed = random_seed,\n",
    "        direc = output_dir,\n",
    "        write=True#, regenerate=True\n",
    ")\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "fig, axs = plt.subplots(len(lightcone_quantities),1,\n",
    "            figsize=(20,10))#(getattr(lightcone_fid_vcb, lightcone_quantities[0]).shape[2]*0.01,\n",
    "                     #getattr(lightcone_fid_vcb, lightcone_quantities[0]).shape[1]*0.01*len(lightcone_quantities)))\n",
    "for ii, lightcone_quantity in enumerate(lightcone_quantities):\n",
    "    axs[ii].imshow(getattr(lightcone_fid_vcb, lightcone_quantity)[1],\n",
    "                   vmin=vmins[ii], vmax=vmaxs[ii],cmap=cmaps[ii])\n",
    "    axs[ii].text(1, 0.05, lightcone_quantity,horizontalalignment='right',verticalalignment='bottom',\n",
    "            transform=axs[ii].transAxes,color = 'red',backgroundcolor='white',fontsize = 15)\n",
    "    axs[ii].xaxis.set_tick_params(labelsize=10)\n",
    "    axs[ii].yaxis.set_tick_params(labelsize=0)\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#also run one without velocities and with fixed vcb=vavg (for comparison) \n",
    "\n",
    "lightcone_fid_novcb = p21c.run_lightcone(\n",
    "        redshift = ZMIN,\n",
    "        init_box = initial_conditions,\n",
    "        flag_options = flag_options_fid,\n",
    "        astro_params = astro_params_novcb,\n",
    "        lightcone_quantities=lightcone_quantities,\n",
    "        global_quantities=lightcone_quantities,\n",
    "        random_seed = random_seed,\n",
    "        direc = output_dir,\n",
    "        write=True#, regenerate=True\n",
    ")\n",
    "\n",
    "lightcone_fid_vcbavg = p21c.run_lightcone(\n",
    "        redshift = ZMIN,\n",
    "        init_box = initial_conditions,\n",
    "        flag_options = flag_options_fid_vavg,\n",
    "        astro_params = astro_params_vcb,\n",
    "        lightcone_quantities=lightcone_quantities,\n",
    "        global_quantities=lightcone_quantities,\n",
    "        random_seed = random_seed,\n",
    "        direc = output_dir,\n",
    "        write=True#, regenerate=True\n",
    ")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plus run one with only atomic-cooling galaxies but same otherwise\n",
    "\n",
    "flag_options_NOMINI=flag_options_fid\n",
    "flag_options_NOMINI.update({'USE_MINI_HALOS': False})\n",
    "\n",
    "lightcone_fid_NOMINI = p21c.run_lightcone(\n",
    "        redshift = ZMIN,\n",
    "        init_box = initial_conditions,\n",
    "        flag_options = flag_options_NOMINI,\n",
    "        astro_params = astro_params_vcb,\n",
    "        lightcone_quantities=lightcone_quantities,\n",
    "        global_quantities=lightcone_quantities,\n",
    "        random_seed = random_seed,\n",
    "        direc = output_dir,\n",
    "        write=True#, regenerate=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compare vcb and novcb\n",
    "\n",
    "fig, axs = plt.subplots(2 ,1,\n",
    "            figsize=(20,6))\n",
    "\n",
    "axs[0].imshow(getattr(lightcone_fid_vcb, 'brightness_temp')[1],\n",
    "               vmin=vmins[0], vmax=vmaxs[0],cmap=cmaps[0])\n",
    "axs[1].imshow(getattr(lightcone_fid_novcb, 'brightness_temp')[1],\n",
    "               vmin=vmins[0], vmax=vmaxs[0],cmap=cmaps[0])\n",
    "axs[0].text(1, 0.05, 'vcb' ,horizontalalignment='right',verticalalignment='bottom',\n",
    "        transform=axs[0].transAxes,color = 'red',backgroundcolor='white',fontsize = 15)\n",
    "axs[1].text(1, 0.05, 'novcb' ,horizontalalignment='right',verticalalignment='bottom',\n",
    "        transform=axs[1].transAxes,color = 'red',backgroundcolor='white',fontsize = 15)                        \n",
    "# axs[0].xaxis.set_tick_params(labelsize=10)\n",
    "# axs[1].yaxis.set_tick_params(labelsize=0)\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot tau\n",
    "\n",
    "tau_vcb=tau_novcb=tau_NOMINI=np.array([])\n",
    "for il,lightcone in enumerate([lightcone_fid_vcb,lightcone_fid_novcb,lightcone_fid_NOMINI]):\n",
    "    z_e=np.array([]);\n",
    "    tau_e=np.array([]);\n",
    "    for i in range(len(lightcone.node_redshifts)-1):\n",
    "        tauz=p21c.compute_tau(redshifts=lightcone.node_redshifts[-1:-2-i:-1], \n",
    "                             global_xHI=lightcone.global_xHI[-1:-2-i:-1]) \n",
    "        tau_e=np.append(tau_e,tauz)\n",
    "        z_e=np.append(z_e,lightcone.node_redshifts[-2-i])\n",
    "\n",
    "    #add lower zs where we manually set xH=1\n",
    "    zlow=np.linspace(lightcone.node_redshifts[-1]-0.1, 0.1, 14)\n",
    "    for zl in zlow:\n",
    "        tauz=p21c.compute_tau(redshifts=np.array([zl]), global_xHI=np.array([lightcone.global_xHI[-1]])) \n",
    "        tau_e=np.append([tauz],tau_e)\n",
    "        z_e=np.append([zl],z_e)\n",
    "        \n",
    "    if(il==0):\n",
    "        tau_vcb=tau_e \n",
    "    elif (il==1):\n",
    "        tau_novcb=tau_e        \n",
    "    else:\n",
    "        tau_NOMINI=tau_e\n",
    "\n",
    "\n",
    "\n",
    "        \n",
    "        \n",
    "linestyles = ['-', '-.',':']\n",
    "colors     = ['black','gray','#377eb8']\n",
    "lws        = [3,1,2]\n",
    "\n",
    "fig, axs = plt.subplots(1, 1, sharex=True, figsize=(8,4))\n",
    " \n",
    "kk=0\n",
    "axs.plot(z_e, tau_vcb, label = 'vcb',\n",
    "                         color=colors[kk],linestyle=linestyles[kk], lw=lws[kk])\n",
    "kk=1\n",
    "axs.plot(z_e, tau_novcb, label = 'no vcb',\n",
    "                         color=colors[kk],linestyle=linestyles[kk], lw=lws[kk])\n",
    "kk=2\n",
    "axs.plot(z_e, tau_NOMINI, label = 'no MINI',\n",
    "                         color=colors[kk],linestyle=linestyles[kk],lw=lws[kk])\n",
    "\n",
    "axs.set_ylim(0., 0.1)\n",
    "axs.set_xlabel('redshift',fontsize=15)\n",
    "axs.xaxis.set_tick_params(labelsize=15)\n",
    "\n",
    "axs.set_xlim(0.,20.)\n",
    "axs.set_ylabel('$\\\\tau$',fontsize=15)\n",
    "axs.yaxis.set_tick_params(labelsize=15)\n",
    "\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.0,wspace=0.0)\n",
    "\n",
    "tauPmin=0.0561-0.0071\n",
    "tauPmax=0.0561+0.0071\n",
    "axs.axhspan(tauPmin, tauPmax, alpha=0.34, color='black')\n",
    "axs.grid()\n",
    "\n",
    "#Planck2020: tau=0.0561±0.0071"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#check that the tau z=15-30 is below 0.02 as Planck requires\n",
    "print(z_e[-1],z_e[55])\n",
    "tau_vcb[-1]-tau_vcb[55]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linestyles = ['-', '-.',':']\n",
    "colors     = ['black','gray','#377eb8']\n",
    "lws        = [3,1,2]\n",
    "labels     = ['vcb', 'no vcb', 'no MINI']\n",
    "\n",
    "fig, axs = plt.subplots(1, 1, sharex=True, figsize=(8,4))\n",
    " \n",
    "for kk,lightcone in enumerate([lightcone_fid_vcb,lightcone_fid_novcb,lightcone_fid_NOMINI]):\n",
    "    \n",
    "    axs.plot(lightcone.node_redshifts, lightcone.global_xHI, label = labels[kk],\n",
    "                         color=colors[kk],linestyle=linestyles[kk], lw=lws[kk])\n",
    "\n",
    "axs.set_ylim(0., 1.)\n",
    "axs.set_xlabel('redshift',fontsize=15)\n",
    "axs.xaxis.set_tick_params(labelsize=15)\n",
    "\n",
    "axs.set_xlim(5.,20.)\n",
    "axs.set_ylabel('$x_{HI}$',fontsize=15)\n",
    "axs.yaxis.set_tick_params(labelsize=15)\n",
    "\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.0,wspace=0.0)\n",
    "\n",
    "axs.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#choose a redshift to print coeval slices and see if there are VAOs. Usually best then T21~T21min/2\n",
    "zz=zlist21[40]\n",
    "print(zz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We plot a coeval box, but we compare the vcb case against the vcb=vavg, since the no velocity (vcb=0) case has a background evolution that is too different.\n",
    "coeval_fid_vcb = p21c.run_coeval(\n",
    "        redshift = zz,\n",
    "        init_box = initial_conditions,\n",
    "        flag_options = flag_options_fid,\n",
    "        astro_params = astro_params_vcb,\n",
    "        random_seed = random_seed,\n",
    "        direc = output_dir,\n",
    "        write=True#, regenerate=True\n",
    ")\n",
    "\n",
    "coeval_fid_vcbavg = p21c.run_coeval(\n",
    "        redshift = zz,\n",
    "        init_box = initial_conditions,\n",
    "        flag_options = flag_options_fid_vavg,\n",
    "        astro_params = astro_params_vcb,\n",
    "        random_seed = random_seed,\n",
    "        direc = output_dir,\n",
    "        write=True#, regenerate=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T21slice_vcb=coeval_fid_vcb.brightness_temp\n",
    "T21avg_vcb = np.mean(T21slice_vcb)\n",
    "dT21slice_vcb = T21slice_vcb - T21avg_vcb\n",
    "\n",
    "T21slice_novcb=coeval_fid_vcbavg.brightness_temp\n",
    "T21avg_novcb = np.mean(T21slice_novcb)\n",
    "dT21slice_novcb = T21slice_novcb - T21avg_novcb\n",
    "\n",
    "\n",
    "\n",
    "sigma21=np.sqrt(np.var(dT21slice_vcb))\n",
    "\n",
    "T21maxplot =  3.0*sigma21\n",
    "T21minplot = -2.0 * sigma21\n",
    "\n",
    "\n",
    "\n",
    "origin = 'lower'\n",
    "extend = 'both'\n",
    "\n",
    "origin = None\n",
    "extend = 'neither'\n",
    "\n",
    "xx = np.linspace(0, BOX_LEN, HII_DIM, endpoint=False)\n",
    "yy = xx\n",
    "\n",
    "\n",
    "\n",
    "indexv=0\n",
    "\n",
    "fig, ax = plt.subplots(2, 2, constrained_layout=True, figsize=(10,8), \n",
    "                       sharex='col', sharey='row',\n",
    "                        gridspec_kw={'hspace': 0, 'wspace': 0})\n",
    "\n",
    "cs0=ax[0,0].contourf(xx, yy, dT21slice_novcb[indexv], extend=extend, origin=origin, \n",
    "                   vmin=T21minplot, vmax=T21maxplot,cmap='bwr') \n",
    "fig.colorbar(cs0, ax=ax[0,0], shrink=0.9, location='left')\n",
    "cs1=ax[0,1].contourf(xx, yy, dT21slice_vcb[indexv], extend=extend, origin=origin, \n",
    "                   vmin=T21minplot, vmax=T21maxplot,cmap='bwr')\n",
    "fig.colorbar(cs1, ax=ax[0,1], shrink=0.9)\n",
    "\n",
    "\n",
    "\n",
    "deltaslice=initial_conditions.lowres_density\n",
    "deltaavg = np.mean(deltaslice)\n",
    "ddeltaslice = deltaslice - deltaavg\n",
    "\n",
    "vcbslice=initial_conditions.lowres_vcb\n",
    "vcbavg = np.mean(vcbslice)\n",
    "dvcbslice = vcbslice\n",
    "\n",
    "print(vcbavg)\n",
    "\n",
    "\n",
    "csd=ax[1,0].contourf(xx, yy, ddeltaslice[indexv]) \n",
    "fig.colorbar(csd, ax=ax[1,0], shrink=0.9, location='left')\n",
    "csv=ax[1,1].contourf(xx, yy, dvcbslice[indexv])\n",
    "fig.colorbar(csv, ax=ax[1,1], shrink=0.9, extend=extend)\n",
    "plt.show()\n",
    "\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "global_quantities = ('brightness_temp','Ts_box','xH_box',\"dNrec_box\",'z_re_box','Gamma12_box','J_21_LW_box',\"density\")\n",
    "#choose some to plot...\n",
    "plot_quantities = ('brightness_temp','Ts_box','xH_box',\"dNrec_box\",'Gamma12_box','J_21_LW_box')\n",
    "ymins = [-120, 1e1, 0, 0, 0, 0]\n",
    "ymaxs = [  30, 1e3, 1, 1, 1,5]\n",
    "linestyles = ['-', '-',':','-.','-.',':']\n",
    "colors     = ['gray','black','#e41a1c','#377eb8','#e41a1c','#377eb8']\n",
    "lws        = [2,2,2,2]\n",
    "\n",
    "textss   = ['vcb','MCGs']\n",
    "factorss = [[0, 1],] * len(textss)\n",
    "labelss  = [['NO', 'reference'],] * len(textss)\n",
    "\n",
    "\n",
    "fig, axss = plt.subplots(len(plot_quantities), len(labelss),\n",
    "                         sharex=True, figsize=(4*len(labelss),2*len(plot_quantities)))\n",
    "\n",
    "for pp, texts in enumerate(textss):\n",
    "    labels  = labelss[pp]\n",
    "    factors = factorss[pp] \n",
    "    axs     = axss[:,pp]\n",
    "    for kk, label in enumerate(labels):\n",
    "        factor = factors[kk]\n",
    "\n",
    "        \n",
    "        if kk==0:\n",
    "            if pp == 0:\n",
    "                lightcone = lightcone_fid_NOMINI\n",
    "            else:          \n",
    "                lightcone = lightcone_fid_novcb\n",
    "        else:\n",
    "            lightcone = lightcone_fid_vcb\n",
    "            \n",
    "        freqs = 1420.4 / (np.array(lightcone.node_redshifts) + 1.)    \n",
    "        for jj, global_quantity in enumerate(plot_quantities):\n",
    "            axs[jj].plot(freqs, getattr(lightcone, 'global_%s'%global_quantity.replace('_box','')),\n",
    "                         color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])\n",
    "            \n",
    "    axs[0].text(0.01, 0.99, texts,horizontalalignment='right',verticalalignment='bottom',\n",
    "                        transform=axs[0].transAxes,fontsize = 15)\n",
    "    for jj, global_quantity in enumerate(plot_quantities):\n",
    "        axs[jj].set_ylim(ymins[jj], ymaxs[jj])\n",
    "    axs[-1].set_xlabel('Frequency/MHz',fontsize=15)\n",
    "    axs[-1].xaxis.set_tick_params(labelsize=15)\n",
    "\n",
    "    axs[0].set_xlim(1420.4 / (35 + 1.), 1420.4 / (5.5 + 1.))\n",
    "    zlabels = np.array([ 6,  7,  8,  10, 13, 18, 25, 35])\n",
    "    ax2 = axs[0].twiny()\n",
    "    ax2.set_xlim(axs[0].get_xlim())\n",
    "    ax2.set_xticks(1420.4 / (zlabels + 1.))\n",
    "    ax2.set_xticklabels(zlabels.astype(np.str)) \n",
    "    ax2.set_xlabel(\"redshift\",fontsize=15)\n",
    "    ax2.xaxis.set_tick_params(labelsize=15)\n",
    "    ax2.grid(False)\n",
    "    \n",
    "    if pp == 0:\n",
    "        axs[0].legend(loc='lower right', ncol=2,fontsize=13,fancybox=True,frameon=True)\n",
    "        for jj, global_quantity in enumerate(plot_quantities):\n",
    "            axs[jj].set_ylabel('global_%s'%global_quantity.replace('_box',''),fontsize=15)\n",
    "            axs[jj].yaxis.set_tick_params(labelsize=15)\n",
    "    else:\n",
    "        for jj, global_quantity in enumerate(plot_quantities):\n",
    "            axs[jj].set_ylabel('global_%s'%global_quantity.replace('_box',''),fontsize=0)\n",
    "            axs[jj].yaxis.set_tick_params(labelsize=0)\n",
    "\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.0,wspace=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " # varying parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "let's vary the parameters that describe mini-halos and see the impact to the global signal.\n",
    "Warning: It may take a while to run all these boxes!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We keep other parameters fixed and vary one of following  by a factor of 1/3 and 3:\n",
    "\n",
    "    F_STAR7_MINI\n",
    "    F_ESC7_MINI\n",
    "    L_X_MINI\n",
    "    A_LW\n",
    "    \n",
    "We also have a NOmini model where mini-halos are not included"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#defining those color, linstyle, blabla\n",
    "linestyles = ['-', '-',':','-.','-.',':']\n",
    "colors     = ['gray','black','#e41a1c','#377eb8','#e41a1c','#377eb8']\n",
    "lws        = [1,3,2,2,2,2]\n",
    "\n",
    "textss   = ['varying '+r'$f_{*,7}^{\\rm mol}$',\\\n",
    "            'varying '+r'$f_{\\rm esc}^{\\rm mol}$',\\\n",
    "            'varying '+r'$L_{\\rm x}^{\\rm mol}$',\\\n",
    "            'varying '+r'$A_{\\rm LW}$']\n",
    "factorss = [[0, 1, 0.33, 3.0],] * len(textss)\n",
    "labelss  = [['No Velocity', 'Fiducial', '/3', 'x3'],] * len(textss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-14T21:28:24.543996Z",
     "start_time": "2020-03-14T21:28:24.540228Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "global_quantities = ('brightness_temp','Ts_box','xH_box',\"dNrec_box\",'z_re_box','Gamma12_box','J_21_LW_box',\"density\")\n",
    "#choose some to plot...\n",
    "plot_quantities = ('brightness_temp','Ts_box','xH_box',\"dNrec_box\",'Gamma12_box','J_21_LW_box')\n",
    "ymins = [-120, 1e1, 0, 0, 0, 0]\n",
    "ymaxs = [  30, 1e3, 1, 1, 1,10]\n",
    "\n",
    "fig, axss = plt.subplots(len(plot_quantities), len(labelss),\n",
    "                         sharex=True, figsize=(4*len(labelss),2*len(global_quantities)))\n",
    "\n",
    "for pp, texts in enumerate(textss):\n",
    "    labels  = labelss[pp]\n",
    "    factors = factorss[pp] \n",
    "    axs     = axss[:,pp]\n",
    "    for kk, label in enumerate(labels):\n",
    "        flag_options = flag_options_fid.copy()\n",
    "        astro_params = astro_params_vcb.copy()\n",
    "        factor = factors[kk]\n",
    "        if label == 'No Velocity':\n",
    "            lightcone = lightcone_fid_novcb\n",
    "        elif label == 'Fiducial':\n",
    "            lightcone = lightcone_fid_vcb\n",
    "        else:\n",
    "            if pp == 0:\n",
    "                astro_params.update({'F_STAR7_MINI': astro_params_vcb['F_STAR7_MINI']+np.log10(factor)})\n",
    "            elif pp == 1:\n",
    "                astro_params.update({'F_ESC7_MINI': astro_params_vcb['F_ESC7_MINI']+np.log10(factor)})\n",
    "            elif pp == 2:\n",
    "                astro_params.update({'L_X_MINI': astro_params_vcb['L_X_MINI']+np.log10(factor)})\n",
    "            elif pp == 3:\n",
    "                astro_params.update({'A_LW': astro_params_vcb['A_LW']*factor})\n",
    "            else:\n",
    "                print('Make a choice!')                \n",
    "\n",
    "            lightcone = p21c.run_lightcone(\n",
    "                redshift = ZMIN,\n",
    "                init_box = initial_conditions,\n",
    "                flag_options = flag_options_fid,\n",
    "                astro_params = astro_params,\n",
    "                global_quantities=global_quantities,\n",
    "                random_seed = random_seed,\n",
    "                direc = output_dir\n",
    "            )\n",
    "            \n",
    "            \n",
    "\n",
    "        freqs = 1420.4 / (np.array(lightcone.node_redshifts) + 1.)\n",
    "        for jj, global_quantity in enumerate(plot_quantities):\n",
    "            axs[jj].plot(freqs, getattr(lightcone, 'global_%s'%global_quantity.replace('_box','')),\n",
    "                         color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])\n",
    "            \n",
    "    axs[0].text(0.01, 0.99, texts,horizontalalignment='left',verticalalignment='top',\n",
    "                        transform=axs[0].transAxes,fontsize = 15)\n",
    "    for jj, global_quantity in enumerate(plot_quantities):\n",
    "        axs[jj].set_ylim(ymins[jj], ymaxs[jj])\n",
    "    axs[-1].set_xlabel('Frequency/MHz',fontsize=15)\n",
    "    axs[-1].xaxis.set_tick_params(labelsize=15)\n",
    "\n",
    "    axs[0].set_xlim(1420.4 / (35 + 1.), 1420.4 / (5.5 + 1.))\n",
    "    zlabels = np.array([ 6,  7,  8,  10, 13, 18, 25, 35])\n",
    "    ax2 = axs[0].twiny()\n",
    "    ax2.set_xlim(axs[0].get_xlim())\n",
    "    ax2.set_xticks(1420.4 / (zlabels + 1.))\n",
    "    ax2.set_xticklabels(zlabels.astype(np.str)) \n",
    "    ax2.set_xlabel(\"redshift\",fontsize=15)\n",
    "    ax2.xaxis.set_tick_params(labelsize=15)\n",
    "    ax2.grid(False)\n",
    "    \n",
    "    if pp == 0:\n",
    "        axs[0].legend(loc='lower right', ncol=2,fontsize=13,fancybox=True,frameon=True)\n",
    "        for jj, global_quantity in enumerate(plot_quantities):\n",
    "            axs[jj].set_ylabel('global_%s'%global_quantity.replace('_box',''),fontsize=15)\n",
    "            axs[jj].yaxis.set_tick_params(labelsize=15)\n",
    "    else:\n",
    "        for jj, global_quantity in enumerate(plot_quantities):\n",
    "            axs[jj].set_ylabel('global_%s'%global_quantity.replace('_box',''),fontsize=0)\n",
    "            axs[jj].yaxis.set_tick_params(labelsize=0)\n",
    "\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.0,wspace=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-14T23:04:29.138653Z",
     "start_time": "2020-03-14T23:04:28.772417Z"
    }
   },
   "outputs": [],
   "source": [
    "# define functions to calculate PS, following py21cmmc\n",
    "from powerbox.tools import get_power\n",
    "\n",
    "def compute_power(\n",
    "   box,\n",
    "   length,\n",
    "   n_psbins,\n",
    "   log_bins=True,\n",
    "   ignore_kperp_zero=True,\n",
    "   ignore_kpar_zero=False,\n",
    "   ignore_k_zero=False,\n",
    "):\n",
    "    # Determine the weighting function required from ignoring k's.\n",
    "    k_weights = np.ones(box.shape, dtype=int)\n",
    "    n0 = k_weights.shape[0]\n",
    "    n1 = k_weights.shape[-1]\n",
    "\n",
    "    if ignore_kperp_zero:\n",
    "        k_weights[n0 // 2, n0 // 2, :] = 0\n",
    "    if ignore_kpar_zero:\n",
    "        k_weights[:, :, n1 // 2] = 0\n",
    "    if ignore_k_zero:\n",
    "        k_weights[n0 // 2, n0 // 2, n1 // 2] = 0\n",
    "\n",
    "    res = get_power(\n",
    "        box,\n",
    "        boxlength=length,\n",
    "        bins=n_psbins,\n",
    "        bin_ave=False,\n",
    "        get_variance=False,\n",
    "        log_bins=log_bins,\n",
    "        k_weights=k_weights,\n",
    "    )\n",
    "\n",
    "    res = list(res)\n",
    "    k = res[1]\n",
    "    if log_bins:\n",
    "        k = np.exp((np.log(k[1:]) + np.log(k[:-1])) / 2)\n",
    "    else:\n",
    "        k = (k[1:] + k[:-1]) / 2\n",
    "\n",
    "    res[1] = k\n",
    "    return res\n",
    "\n",
    "def powerspectra(brightness_temp, n_psbins=50, nchunks=10, min_k=0.1, max_k=1.0, logk=True):\n",
    "    data = []\n",
    "    chunk_indices = list(range(0,brightness_temp.n_slices,round(brightness_temp.n_slices / nchunks),))\n",
    "\n",
    "    if len(chunk_indices) > nchunks:\n",
    "        chunk_indices = chunk_indices[:-1]\n",
    "    chunk_indices.append(brightness_temp.n_slices)\n",
    "\n",
    "    for i in range(nchunks):\n",
    "        start = chunk_indices[i]\n",
    "        end = chunk_indices[i + 1]\n",
    "        chunklen = (end - start) * brightness_temp.cell_size\n",
    "\n",
    "        power, k = compute_power(\n",
    "            brightness_temp.brightness_temp[:, :, start:end],\n",
    "            (BOX_LEN, BOX_LEN, chunklen),\n",
    "            n_psbins,\n",
    "            log_bins=logk,\n",
    "        )\n",
    "        data.append({\"k\": k, \"delta\": power * k ** 3 / (2 * np.pi ** 2)})\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# do 5 chunks but only plot 1 - 4, the 0th has no power for minihalo models where xH=0\n",
    "nchunks = 4\n",
    "k_fundamental = 2*np.pi/BOX_LEN\n",
    "k_max = k_fundamental * HII_DIM\n",
    "Nk=np.floor(HII_DIM/1).astype(int)\n",
    "\n",
    "fig, axss = plt.subplots(nchunks, len(labelss), sharex=True,sharey=True,figsize=(4*len(labelss),3*(nchunks)),subplot_kw={\"xscale\":'log', \"yscale\":'log'})\n",
    "\n",
    "for pp, texts in enumerate(textss):\n",
    "    labels  = labelss[pp]\n",
    "    factors = factorss[pp] \n",
    "    axs     = axss[:,pp]\n",
    "    for kk, label in enumerate(labels):\n",
    "        flag_options = flag_options_fid.copy()\n",
    "        astro_params = astro_params_vcb.copy()\n",
    "        factor = factors[kk]\n",
    "        if label == 'No Velocity':\n",
    "            lightcone = lightcone_fid_novcb\n",
    "        elif label == 'Fiducial':\n",
    "            lightcone = lightcone_fid_vcb\n",
    "        else:\n",
    "            if pp == 0:\n",
    "                astro_params.update({'F_STAR7_MINI': astro_params_vcb['F_STAR7_MINI']+np.log10(factor)})\n",
    "            elif pp == 1:\n",
    "                astro_params.update({'F_ESC7_MINI': astro_params_vcb['F_ESC7_MINI']+np.log10(factor)})\n",
    "            elif pp == 2:\n",
    "                astro_params.update({'L_X_MINI': astro_params_vcb['L_X_MINI']+np.log10(factor)})\n",
    "            elif pp == 3:\n",
    "                astro_params.update({'A_LW': astro_params_vcb['A_LW']+np.log10(factor)})\n",
    "            else:\n",
    "                print('Make a choice!')                \n",
    "\n",
    "            lightcone = p21c.run_lightcone(\n",
    "                redshift = ZMIN,\n",
    "                init_box = initial_conditions,\n",
    "                flag_options = flag_options_fid,\n",
    "                astro_params = astro_params,\n",
    "                global_quantities=global_quantities,\n",
    "                random_seed = random_seed,\n",
    "                direc = output_dir\n",
    "            )\n",
    "            \n",
    "        PS = powerspectra(lightcone, min_k = k_fundamental, max_k = k_max)\n",
    "        \n",
    "        for ii in range(nchunks):\n",
    "            axs[ii].plot(PS[ii+1]['k'], PS[ii+1]['delta'], color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])\n",
    "        \n",
    "            if pp == len(textss)-1 and kk == 0:\n",
    "                axs[ii].text(0.99, 0.01, 'Chunk-%02d'%(ii+1),horizontalalignment='right',verticalalignment='bottom',\n",
    "                        transform=axs[ii].transAxes,fontsize = 15)\n",
    "            \n",
    "    axs[0].text(0.01, 0.99, texts,horizontalalignment='left',verticalalignment='top',\n",
    "                        transform=axs[0].transAxes,fontsize = 15)\n",
    "\n",
    "    axs[-1].set_xlabel(\"$k$ [Mpc$^{-3}$]\",fontsize=15)\n",
    "    axs[-1].xaxis.set_tick_params(labelsize=15)\n",
    "    \n",
    "    if pp == 0:\n",
    "        for ii in range(nchunks):\n",
    "            axs[ii].set_ylim(2e-1, 2e2)\n",
    "            axs[ii].set_ylabel(\"$k^3 P$\", fontsize=15)\n",
    "            axs[ii].yaxis.set_tick_params(labelsize=15)\n",
    "    else:\n",
    "        for ii in range(nchunks-1):\n",
    "            axs[ii].set_ylim(2e-1, 2e2)\n",
    "            axs[ii].set_ylabel(\"$k^3 P$\", fontsize=0)\n",
    "            axs[ii].yaxis.set_tick_params(labelsize=0)\n",
    "\n",
    "axss[0,0].legend(loc='lower left', ncol=2,fontsize=13,fancybox=True,frameon=True)\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.0,wspace=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that I've run these simulations in parallel before this tutorial. With these setup, each took ~6h to finish. Here, running means read the cached outputs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## global properties -- optical depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#defining those color, linstyle, blabla\n",
    "linestyles = ['-', '-',':','-.','-.',':']\n",
    "colors     = ['gray','black','#e41a1c','#377eb8','#e41a1c','#377eb8']\n",
    "lws        = [1,3,2,2,2,2]\n",
    "\n",
    "textss_tau   = ['varying '+r'$f_{*,7}^{\\rm mol}$',\\\n",
    "            'varying '+r'$f_{\\rm esc}^{\\rm mol}$',\\\n",
    "            'varying '+r'$A_{\\rm LW}$']\n",
    "\n",
    "factorss_tau = [[0, 1, 0.33, 3.0],] * len(textss_tau)\n",
    "labelss_tau  = [['No Velocity', 'Fiducial', '/3', 'x3'],] * len(textss_tau)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-14T23:03:40.364465Z",
     "start_time": "2020-03-14T23:03:38.244179Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_quantities = ['tau_e']\n",
    "ymins = [0]\n",
    "ymaxs = [0.2]\n",
    "\n",
    "\n",
    "fig, axss_tau = plt.subplots(len(plot_quantities), len(labelss_tau),\n",
    "                         sharex=True, figsize=(4*len(labelss_tau),3*len(plot_quantities)))\n",
    "\n",
    "\n",
    "\n",
    "for pp, texts in enumerate(textss_tau):\n",
    "    labels  = labelss_tau[pp]\n",
    "    factors = factorss_tau[pp] \n",
    "    axs     = axss_tau[pp]\n",
    "    for kk, label in enumerate(labels):\n",
    "        flag_options = flag_options_fid.copy()\n",
    "        astro_params = astro_params_vcb.copy()\n",
    "        factor = factors[kk]\n",
    "        if label == 'No Velocity':\n",
    "            lightcone = lightcone_fid_novcb\n",
    "        elif label == 'Fiducial':\n",
    "            lightcone = lightcone_fid_vcb\n",
    "        else:\n",
    "            if pp == 0:\n",
    "                astro_params.update({'F_STAR7_MINI': astro_params_vcb['F_STAR7_MINI']+np.log10(factor)})\n",
    "            elif pp == 1:\n",
    "                astro_params.update({'F_ESC7_MINI': astro_params_vcb['F_ESC7_MINI']+np.log10(factor)})\n",
    "            elif pp == 2:\n",
    "                astro_params.update({'A_LW': astro_params_vcb['A_LW']*factor})\n",
    "            else:\n",
    "                print('Make a choice!')                \n",
    "\n",
    "            lightcone = p21c.run_lightcone(\n",
    "                redshift = ZMIN,\n",
    "                init_box = initial_conditions,\n",
    "                flag_options = flag_options_fid,\n",
    "                astro_params = astro_params,\n",
    "                global_quantities=global_quantities,\n",
    "                random_seed = random_seed,\n",
    "                direc = output_dir\n",
    "            )\n",
    "            \n",
    "\n",
    "        z_e=np.array([]);\n",
    "        tau_e=np.array([]);\n",
    "        for i in range(len(lightcone.node_redshifts)-1):\n",
    "            tauz=p21c.compute_tau(redshifts=lightcone.node_redshifts[-1:-2-i:-1], \n",
    "                                 global_xHI=lightcone.global_xHI[-1:-2-i:-1]) \n",
    "            tau_e=np.append(tau_e,tauz)\n",
    "            z_e=np.append(z_e,lightcone.node_redshifts[-2-i])\n",
    "            #print(i,lightcone.node_redshifts[i],tauz)\n",
    "            \n",
    "        #add lower zs where we manually set xH=1\n",
    "        zlow=np.linspace(lightcone.node_redshifts[-1]-0.1, 0.1, 14)\n",
    "        for zl in zlow:\n",
    "            tauz=p21c.compute_tau(redshifts=np.array([zl]), global_xHI=np.array([lightcone.global_xHI[-1]])) \n",
    "            tau_e=np.append([tauz],tau_e)\n",
    "            z_e=np.append([zl],z_e)\n",
    "        \n",
    "            \n",
    "#        freqs = 1420.4 / (np.array(lightcone.node_redshifts) + 1.)\n",
    "        for jj, global_quantity in enumerate(plot_quantities):\n",
    "            axs.plot(z_e, tau_e,\n",
    "                         color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])\n",
    "   \n",
    "        \n",
    "    axs.text(0.01, 0.99, texts,horizontalalignment='left',verticalalignment='top',\n",
    "                        transform=axs.transAxes,fontsize = 15)\n",
    "    axs.set_ylim(ymins[0], ymaxs[0])\n",
    "    axs.set_xlabel('redshift',fontsize=15)\n",
    "    axs.xaxis.set_tick_params(labelsize=15)\n",
    "\n",
    "    axs.set_xlim(0.,20.)\n",
    "\n",
    "    if pp == 0:\n",
    "        for ii in range(nchunks):\n",
    "            axs.set_ylabel('$\\\\tau$',fontsize=15)\n",
    "            axs.yaxis.set_tick_params(labelsize=15)\n",
    "    else:\n",
    "        for ii in range(nchunks-1):\n",
    "            axs.yaxis.set_tick_params(labelsize=0)\n",
    "\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.0,wspace=0.0)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 21-cm power spectra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-15T00:04:49.728969Z",
     "start_time": "2020-03-15T00:04:47.792544Z"
    }
   },
   "outputs": [],
   "source": [
    "# do 5 chunks but only plot 1 - 4, the 0th has no power for minihalo models where xH=0\n",
    "nchunks = 4\n",
    "\n",
    "fig, axss = plt.subplots(nchunks, len(labelss), sharex=True,sharey=True,figsize=(4*len(labelss),3*(nchunks)),subplot_kw={\"xscale\":'log', \"yscale\":'log'})\n",
    "\n",
    "for pp, texts in enumerate(textss):\n",
    "    labels  = labelss[pp]\n",
    "    factors = factorss[pp] \n",
    "    axs     = axss[:,pp]\n",
    "    for kk, label in enumerate(labels):\n",
    "        factor = factors[kk]\n",
    "        \n",
    "        if kk==0:\n",
    "            lightcone = lightcone_fid_vcbavg\n",
    "        else:\n",
    "            lightcone = lightcone_fid_vcb\n",
    "        \n",
    "        PS = powerspectra(lightcone, min_k = k_fundamental, max_k = k_max)\n",
    "        for ii in range(nchunks):\n",
    "            axs[ii].plot(PS[ii+1]['k'], PS[ii+1]['delta'], color=colors[kk],linestyle=linestyles[kk], label = labels[kk],lw=lws[kk])\n",
    "        \n",
    "            if pp == len(textss)-1 and kk == 0:\n",
    "                axs[ii].text(0.99, 0.01, 'Chunk-%02d'%(ii+1),horizontalalignment='right',verticalalignment='bottom',\n",
    "                        transform=axs[ii].transAxes,fontsize = 15)\n",
    "            \n",
    "    axs[0].text(0.01, 0.99, texts,horizontalalignment='left',verticalalignment='top',\n",
    "                        transform=axs[0].transAxes,fontsize = 15)\n",
    "\n",
    "    axs[-1].set_xlabel(\"$k$ [Mpc$^{-3}$]\",fontsize=15)\n",
    "    axs[-1].xaxis.set_tick_params(labelsize=15)\n",
    "    \n",
    "    if pp == 0:\n",
    "        for ii in range(nchunks):\n",
    "            axs[ii].set_ylim(2e-1, 2e2)\n",
    "            axs[ii].set_ylabel(\"$k^3 P$\", fontsize=15)\n",
    "            axs[ii].yaxis.set_tick_params(labelsize=15)\n",
    "    else:\n",
    "        for ii in range(nchunks-1):\n",
    "            axs[ii].set_ylim(2e-1, 2e2)\n",
    "            axs[ii].set_ylabel(\"$k^3 P$\", fontsize=0)\n",
    "            axs[ii].yaxis.set_tick_params(labelsize=0)\n",
    "\n",
    "axss[0,0].legend(loc='lower left', ncol=2,fontsize=13,fancybox=True,frameon=True)\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.0,wspace=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nchunks=5\n",
    "\n",
    "textss   = ['vcb','vcb']\n",
    "factorss = [[0, 1],] * len(textss)\n",
    "labelss  = [['Regular', 'Avg'],] * len(textss)\n",
    "\n",
    "k_fundamental = 2*np.pi/BOX_LEN\n",
    "k_max = k_fundamental * HII_DIM\n",
    "Nk=np.floor(HII_DIM/1).astype(int)\n",
    "\n",
    "PSv= powerspectra(lightcone_fid_vcb, min_k = k_fundamental, max_k = k_max)\n",
    "PSvavg= powerspectra(lightcone_fid_vcbavg, min_k = k_fundamental, max_k = k_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "klist= PSv[0]['k']\n",
    "P21diff = [(PSv[i]['delta']-PSvavg[i]['delta'])/PSvavg[i]['delta'] for i in range(nchunks)]\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "fig, axss = plt.subplots(nchunks, 1, sharex=True,sharey=True,figsize=(2*len(labelss),3*(nchunks)),subplot_kw={\"xscale\":'linear', \"yscale\":'linear'})\n",
    "\n",
    "for ii in range(nchunks):\n",
    "    axss[ii].plot(klist, P21diff[ii])\n",
    "\n",
    "plt.xscale('log')\n",
    "axss[0].legend(loc='lower left', ncol=2,fontsize=13,fancybox=True,frameon=True)\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace = 0.0,wspace=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "   "
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Py3.8",
   "language": "python",
   "name": "py3.8"
  },
  "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.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}