{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NGC 1275 Limit Plots \n",
    "This notebook allows the reproduction of Figs. 7 and 8 from J. H. Matthews et al. 2022, \"How do magnetic field models affect astrophysical limits on light axion-like particles? An X-ray case study with NGC 1275\".\n",
    "\n",
    "The data is stored in [alpro/docs/source/data/NGC1275_limits.dat](../data/NGC1275_limits.dat), with column headings as follows:\n",
    "\n",
    "- ``logm``: $\\log_{10} [m_a ({\\rm eV})]$, ALP mass\n",
    "- ``logg_abspl``: 99.7% limit on $\\log_{10} [g_{a \\gamma} ({\\rm GeV}^{-1})]$ for absorbed power-law model \n",
    "- ``logg_N``: 99.7% limit on $\\log_{10} [g_{a \\gamma} ({\\rm GeV}^{-1})]$ for model N\n",
    "- ``logg_r20``: 99.7% limit on $\\log_{10} [g_{a \\gamma} ({\\rm GeV}^{-1})]$ for model B from Reynolds et al. 2020\n",
    "\n",
    "The data is given at 0.1 dex resolution as this was the resolution of the grid of ALP survival probability curves. \n",
    "\n",
    "Note that the plots don't match the format in the paper, but the data is identical. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from astropy.io import ascii \n",
    "import matplotlib.pyplot as plt \n",
    "import alpro\n",
    "data = ascii.read('../data/NGC1275_limits.dat') \n",
    "alpro.util.set_default_plot_params(tex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fig 7: Spectral Model Comparison\n",
    "\n",
    "**Figure caption:**\n",
    "     \n",
    "> 99.7 per cent limits obtained in this work using two different approaches for modelling the X-ray spectrum, zoomed in to the $10^{-14}\\geq m_a/{\\rm eV}< 10^{-11}$ region. The partially covering absorber fits use spectral models of the XSPEC form ``tbabs*tbpcf(ALPs*(pow+zgauss))``, whereas the absorbed power-law models use ``phabs*zphabs(ALPs*pow)``. For comparison, we also show the limits from Reynolds et al. (2020) who used the absorbed power-law method. The limits are not very sensitive to the choice of spectral model, and we reproduce very similar results as Reynolds et al. (2020) when using the equivalent spectral model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAERCAYAAACkWKo8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXRb133o++/GzBmE5lkCZVmW7VgmKTu2k0i2yTiJp8QmJadtemMnFpM2aYYmYtV3373vtvcuhXLffatd6SApzU2bOIlE2k5stx4I2ZKnOBFFyYlHWYIGax5ISJyJYb8/zgEEggAJkiBIUL/PWlwmcQ4O9sGR8cPeZ+/fT2mtEUIIIbLJMtENEEIIceWR4COEECLrJPgIIYTIOgk+Qgghsk6CjxBCiKyzTXQDcoFSah8wAzg40W0RQogcshQ4p7W+MXGDBJ/0zMjPc81bvGDOvIluiBBCTBitsdls2B2OtHb/4OAROru6k26T4JOeg0u9i+f93f/864luhxBCTJje3l4KCgpYsGBBWvt/8ctfY+9b7yYdMZJ7PkIIIbJOgo8QQoisk+AjhBAi6yT4CCGEyDoJPkIIIbJOgo8QQoisk+AjhBAi6yT4CCGEyDoJPkIIIbJOgo8QQoisy7ngo5SqUkrtHem2hP28SqkN5v4blFLuzLdUCCFEKjmV200pVQW0AeUj2ZbEFq11tfk8P9AA1GWwqUIIIYaQUz0frbVPa9060m3xlFJewBP3PD+wNnOtFEIIMZyc6vlkSDlGD2kApZTXDETxj+0yf12ZhXYJIcQVI6d6PhniAQIJj7UBct9HCCGy5EoMPmnTWq/RWq8B9k90W4QQYiqZ8GE3pdR6oGyIXZq11r4MvmSyXk6y3pAQQohxMuHBR2u9Ncsv2UrchIO4dviT7CuEEGIcXBHDbua6HjcMDjLm7LcdE9IwIYS4QuVU8DEXhTaYvzeYa3uG3Yaxjid+OnWtubi0BqjTWssaHyGEyKIJH3YbCfPejw+oH+G22oS//cBm88+mzLdUCCHEUHKq5yOEEGJqkOAjhBAi6yT4CCGEyDoJPkIIIbJOgo8QQoisk+AjhBAi6yT4CCGEyDoJPkIIIbJOgo8QQoisk+AjhBAi6yT4CCGEyDoJPkIIIbJOgo8QQoisk+AjhBAi6yT4CCGEyDoJPkIIIbJOgo8QQoisk+AjhBAi6yT4CCGEyDoJPkIIIbJOgo8QQoisk+AjhBAi6yT4CCGEyDrbRDdACCHExOkNQ39EpbVvX9CCPZz+sa3h3pTbJPgIIcQV6rfnHfz2vCPt/cNhF067ldvD/VQttmO3Jg9aKhyk+Fgz9t4LKY8lwUcIIa5Ax7qs/Pa8g3BEY1Hp9XxA0R+GVz8KcjAQpuZqJ3MKB969sXeeoPTgU9i6zw55JAk+QghxhekOKV486SIcAZdNUWRXOKx62Of1BzXd2oJScOJSmH/Z10vVYju3zbdhIULhidcoOv4KlnAf2mJHK2vKY0nwEUKIK4jWsPO0k86gwm6BEgf8+fIeCuzDP7eru5v3e9282l5IZ3+Y3mCY5w71c+rseR6x/gfF3cdQOkLEVkDYWUJf/hygLemxZLabEEJcQX4fsOPvsKEUOGyKBxb1pxV4AJSCm2dG+B9rSlg+3UGh08oNwf3cfupHnDlxlIt9EHKU0DP9Ok7e+jeEHCUpjyU9HyGEuEJc6LPw2lkHEa3Js1u4dUaIpcWRER9ndqGVjassfPjKr+kLvEeECEFt4XiPncOld7J45R+R7xx6IsOIg49SajFQDqwC3IAHo18VAPYAPq31pZEedwSvXwU0aK0rRrItYb8G4BCwA1gL+LXWvvForxBCTAahCLxw0kV/GBxWxbz8CFVzg6M72Nl3sb39JNeEAnSVWjjWZeeknk6T8wuc6pqD55Uj3Ldy7pCHSDv4KKXuBOqAC0Ar4AP8GIHHA3jNn81KqVJgi9b6pdGdWco2VJmvVz6SbSk0mD+bJPAIIaa6N845OddrwWqBfJuiZnEfNgsE+3oJBfvSOkZ/Tyd5594gHDqKVUfA7qIgv4Cl5Z9hb+hT9B5ox9EX4kJXHz978yiB7v6Uxxo2+CilSoD1wF6t9doUu10EDgM7gW3m8+5USn0P2JqpnlA0SKgk0wKH2pbEHq11/XA7KaV2mb+uTLuRQggxyRzutLK/3U5EQ4FN8dn5/cx09HPq1V9gO/MWVjX8TDcArTUfKehyl3D18mugaDZ8/OvY59zAnwDXLQzwf14/zIXOfjr7gnT2hVIeK50JB2u11o+NtBejtd6ptf47oFopVTyS52aLUirdXpIQQuSkrpDCd8pFOKJx2eAad5hb8k/i2v9jnGf3YdEhtI6k9QMalAXsLvCuhs89BnNuiL3WygVu/vbz13HzEg/uPMeQ64eG7florbeN5cS11k+M5fnjxGsO07WY93+2a61bE3fSWq+BWA9odVZbKIQQY6Q1+E656Aoq7FbwOCJ8Of91Sg++yaGPjhKJaKx2FxGrC2VJb/JzR9jOEeetXH3bt43pbwmKXXa+ccdSXv3wPM0pMiDAFTrbTWu9Ofq7UmoL0AyUTVyLhBAi895qt3Ok04pSMFsFqLf/BzPPneDYqTP0BiP0Kye/6V7O1xp+zrTp0zP2ukopPrVsBounF3AyxT4THnyUUusZ+oO/OdMTApRSbq11AEBr7VdKeTN5fCGEmGjnei28fs5JJBLhVss7PKR2MTPYRdulbtq7+jjcncfzgTK+ufFvMhp44tmtqXtTowo+Sqmvaq1/NOoWxdFab83EcdIVnY4NDDkdWwghclXQnFbtCHWyTu/kY9rPbFsfPWELx85e4tmPinmrfyH33HsfH//4xyekjUmDj1LqB0CqpakK44M7I8EnG8yeTZvZ22kB6uO21QBNE9U2IYTItNfOOpne4+fesI9i1c0cV4iQvZD9B8/yT297uGgp5eqrvXz1q1+dsDam6vkcMv/rT7HdMw5tGZbZa6k2f28gbkhuqG0YPZ1mjGnfAaVUm1JqA8bC2DKtdW2WT0UIIdKnNZDedOijlyIsPv8bbgi/jU1pZuYprK4idp8t4u/2dBFRNqZPL+b73/8+TqdzfNs9hFTBZwdQqbXemWxjmmtpMs4MJj7iei5pbqtN+LsVY6GsEEJMWpZIkLltb+Lp+ACLTr1mJqo3DJ4eK2FtlEkocNpwuafzhnUVj+18hmAYSkuLePjhh/F6J/ZWd9Lgo7W+iLFgNKlUQUkIIURm5PeeYdG5nbj621B66F6PBtqDNgL9VjQalMJms+FcWM4Hs6r4x//9Q/r6+igpKWHVqlXcd9992TmJIUz4bDchhBCXKR1mVnsrswN7sUSCRCx2grY8Iilq44TCmvM9YfrDGq3AaoFuq5v+5Z8ltOBGHv/xjzl//jx5eXnMmjWL73znOxM2ehVvJLnd7sh0rjYhhBCXOfsDLDq3k4Le0ygdIWTNo89WzKE5d3Ou5IYBizq11rx3Psiek730WTU2uyLPrljitnPvsgLynBbeeOMN3nrrLZRSFBYW8q1vfYvS0tIJPMPLRtLzKQck+AghRKZpzfSOd5l74Q1s4V4iykrQVkSgsIwDcx+gz+EesHtXf4TXPurl+KWQmTZHkW+3cPviPFbNdaKU4vTp0zz11FOEQiFKS0u57777WLVq1QSd4GAjCT4T308TQogpxhbqYuG5XRR3H8WiQ4QseYRs+RyZWcWJabcaudTiHA4EeeOjXrqDGgUUOizMKbJx37ICZhYYQ3OhUIif/vSndHV1UVBQwNKlS3n44Ycn4OxSG0nwSW+enxBCiLSUdPlZcG439lAXKEW/rYiuvLl8MK+GLtfsAfv2hzVvHu/lw7Yg4YjGaVW47IqPz3fxqYV52CyX+wfPPvssR48exW63U1JSwve+9z0cjqGLu2WbTDgQQogss0T6mX/+dTwd72HRIcIWJyFrHiem3cbRmXcSsQysa32qM8SrR3u51GdUHS1wWJiWZ+XeZfksLBm47/vvv8/LL79MJBKhtLSURx55hMWLF2fr1NImwUcIIbLo9JEPmHNwO5esveRNL8VROI0e5zQOzHuAiwUD196EI5rWU328fa6fYFhjtyrybIqPzXJS7c3DZRs4JHf8+HEef/xxgsEgxcXF3Hzzzdxzzz3ZPL20SfARQogsCIfDvLT7NW5ofxFnQS89IcWRUxc4UzQX/cn14BpY9qytJ8wrR3s53x1GA3l2RYnTymeX5rN8+sAhtEgkwksvvcRzzz1HT08PLpeL2bNn8+1vf3tSTKtOZiTB5+K4tUIIIaaws+fO0/Tr57lOfcjCeX1YLFZCys5PPixhX/tFZr63hXtqvsjsufPRWvPOuX5aTvXRH9LYLMYU6qUeO3dfVUCRY2Bv58KFCzz++ON8+OGHhMNhCgsLcbvdfP/736ekJFWKzomXdvAZa1E5IYS40mit+c3vWnnxpVeZ5+ik+tpurFYLhcUefhe5loOhABZLF6dPfMRP/+UfWHXH3XTPq+RUZzg2hbrAYeHOJXmUz3YO6MVordmzZw9PPPEEnZ2dWCwWPB4PK1as4Lvf/S7z5s2bwDMf3rDBRylVrLW+lI3GCCHEVHHxUgdPPvMCB/1HcBDkT6/vxmGzMGvuAgquuZObKr5HcPfLPPmLf+NSRwcdzpk0n1A4AoeYNWcuxflO5hcbU6in5Q/MbtDZ2UljYyP79u0jFAqRn59PUVERX/ziF1m7di1Wa/JsCJNJOj2fbcC68W6IEEJMFb9/532eeW4nnV1dWCyKP1nez/wSG/MXLMHmWcTJ6x9FWax84vYqFiy9hv/3F88T6LMTCYXo6+ni9DE/y8uK+dL112JLKMj23nvv8Ytf/IILFy6glMLtdrN48WL+8i//kmXLlk3QGY9cOsGnVinlxyhRIBkOhBAihZ6eXp594SXe+sN7hMIhnA4nt8wJUXVVPrOnlxLJK+XctV8m7DKq0vjPdfLr97soXnYTwXOnaT93BkvvRRwHXuOtty/yb0duYN26dRQWFtLf38/TTz/Nq6++SjAYxOVyUVhYyD333MPDDz+My+Wa4LMfmXSCT53WeptSaolS6lGMxaY+rfWR8W2aEELkDv+RYzz5zAu0tbUDUFBQwNVzivnerf2UODVhRzGd82+ne1YlwXCEne+d5beHLxAMR3DaLFy1eCErblzIged/wjFHmEu9EVpbWzly5Ah33XUXr7zyCidPnkRrTXFxMbNnz+Zb3/oWlZWVE3zmozNs8IlONNBaH8YYgkMpdaNSqhq4gBGI5J6QECKr8vrOM//8KwRtRRybsZqIJXMr+C+0tfPE08/THkh/km9HZxfhUAi73Y7L5eKTt1TwzWvaKeg+TsSWT3/RAtqufohTF3t4at8JTl/sRaMpctmZVujgy7cu5saFpfR99gZ+/OMf8+yzz9LZ2cn58+fZsWMHoVAIh8NBUVERt912G9/4xjcm9Wy24YxqnY/Weh+wD0Ap9aBSygMckmE5IUQ2WMJ9LDnzAq7+NrRSRJSNYzNvz8ixg8EQv3ziGY6fPD2i5ykgLy8ft7uEP1n3BVZ7zpJ//B20shK2F3H2+jpeOdzBrgPn6AuGsVstFDrt3LiwlC/ftpiSPCNTgdPp5Otf/zo33XQTf//3f8/p06fp6OigqKgIt9vN1772Ne68885Ju34nXWNeZKq1fgJAKfWoUqoRo1T1xjG3TAghktGaBRdexdnfTsRiR0VCeDre41L+AgKFS8d8+OaXX+XkqTNYrVby8/JIO6eygpXXr+CP1z3ADH2ewnefROkIYaeb4ws/z7++rTly4QzhiKbAaaMkz85DNy3kU1dNTxpIKioq+OEPf8g///M/8/rrr3P99dfzF3/xF8yaNWvM5zgZjCn4KKVWYsyEW49xhbaaP0IIMS5KOw9Q2nEAFASteVgsYeyhThacf4Uu5yyC9qJRH/vAocO88btWIlpT6Mpj7QP3sqpiZVrPtdtt5OfloYLdlP7+V1jC/YQdhRx2XMUP/F66gp1YlcKdb2fZrCIe/aSXmcVDTxIoLi6mvr6e/v7+SZcYdKxGHHyUUsUYweZrwBKgCVgrpbWFEOPNEbzIgvOvYtFhgrZCzpXcQFHPcQojIeyhLhaf28mHc+4bVIYgHZ1dXTz59POxdTMrP3Ytd675xMiGt7TGffhZrL1tBC1OjnTn87/7q+kiTJ7DSoHDxv03zuPu6+dgtaR/3KkWeACGvUJKqe+Z/31AKfUC0I7R22kAPFrrdRJ4hBDjTekwi8/uxBbuIWxx0uWaxYdz7+eDeTUEbfloZaGw5wSzAq0jPrbWmqeefZFLlzpwOBxM83j40y8+OOL7Kvnn9pF3/h06g3D4kuJnfJZLqpCSPDveGYX813tWcN8Nc0cUeKaqdL4e/LVSKgxsxphksFRrvUprvU1rLfnehBBZMbt9LwW9p8wqnwV8MK+GiMVBR/4Cjs24naC1AKUjzGlvoaB3ZJMF3mzZz/sHDqGUwuV08vCX1lJUVDiiY1h7zlPof44zXWE+6rbxG2sFR/NW4M638+lrZ/Pf713BkukFIzrmVJbOsJsfeNSc4SaEEFlX0HOSWYFWlI4QtBdxZGYVnXmXc5d9NP1TuLsOUdpxAHu4h0Vnd/L+vBoiVuewxz599hwv+HYTCYcpKCzk03euZsXVV42sgZEw9nee4FhbL33ayjnrDN5038284jweuW0J18/P3SnR4yWd4LNFAo8QYqJYw30sOrcTS6SfoLWA9sKrjKJrJ8/y1G8PUGC3cN8tK3DOe5AbD/0TVh3CGWxnwYVXOTqzashjB4MhGp/6T/r6+nG6XCxeOJ/7776LtjPH6Xp/J46+9vTaSAhH/0U0ELY62Tnti9y0dDZ/essiilz2YZ9/JUp7kWkipdSDQB2wRGt9lfnYHbLWRwiRMVqz4PxunP0XiVgc9NlLODD3Abr6+vjpawfoCluhF/6p+V1uX1ZKUdl9LD/+SxyhDko7DnApbyHtRanznb2w8xVOnTmL1WalsKCAR760jlPvvErxid0UR0IjayoKi83OHs9n+cLtt3Jr2bScX4sznkY11dpMs1MG1AOx0nta65eUUg9orZ/MUPuEEFcwT+f7uDsPxqZVfzj38/Tbi/n5C2/QFbZiUQplUfSFwrx44CLvn4jwv667niVd+7GHOpl//hW6XLPptxcPOvYHH/p5s2UfaE2eK48H761Cv/c07u6P0FpjsdqwWEaWHTowo5LateuZXpRbedYmwmjX+bRFe0RKqdKEbRLqhRBj5uwPMP/8a1h0iKCtiNOlN9FWfA27W9/jcKcF0DjsNvJsii6Lhf5giKOdij//3Sz+uSyfJY6Qef/Hx8G596PV5UByqaOTJ595gXAohCs/j898vIyV3a9iCfeBUljtDvrcVzGv6us489KbeGB3OFk5fQ5Ibyctow0+8QOhie90YjASQogRUTrMonM7sYV7CVny6HTNxj/7M5w+387OA+1EIuB02PnYnHz+4Zs1fPcfn6D1RBf9oRCBsIP/60g5fzv9RbxuC4W9p5jVvpfTnpsAc1r1My/Q0dFBcVE+6ypKqSg5DmGNstqwWB1Erq9h1af/C1bbmJPAiBRGvhLLUGFmNwAjyzUQy3ggwUcIMSaz2/ZQ0HuaiLISsuXzwfxa+rWVn+1+m6BW2GxW3E4L//DNGmaUFvPT//owf/m56yl22XDYbRwOevhZ4HrePx+ks6ef2YFWCnpOAvDG71o5cOgwV893892P21jp7gYFVruDUNF8Zjz4GBWf+4oEnnE22sSijymldiilbgQCZr0fL8Zw3F0ZbWECpVQV0KC1rkh4vByITm1ZhTE9PJDiGF6gBmgFyjHy0SXdVwiRXYU9x5l1cX9sWvXhWZ+myzWHp15uoa3fhgKcNhv/z0OfYEbp5Xs5j9y7mtvLl/Ptf/o1B9v6eaH3WlbYTxHWZ5jW28ki6052229n1yuvc//KUm6d3oXDYcNmt2O12ugru4uV934dpyt/4k7+CjLq0K61XquUWoLxge/G+AAf10wHZuBpwwgY8Y+7gUqt9Wbz7xpgJ1Ax6CCGLVrranNfP0a2hrrxarcQIj3WcC+Lzr4Um1bdVriMk55bePvQMfafCRHRmjyng89eO5O7brlh0POXzJvFU3/7VR57/Dl+8eZhftr9Sf7a9gwEofPkWcKXnuDrNzmZ4ejEarWSl1dAOM9DyR3fYuUNt03AGV+5RjvsBsRq/DRrrR/LRoodrbVPa50sd0Ylxsy7KB9QbgalAcxejyfumH5gbabbKkSmaK3Z/4d3+befP8Frb7YQDkcyevziriN4T/0Hs9v3oEY4vXg4/iPH+Nn2X/G8bzd9ff1D7qsiQRadexlH8BJhi4M+h5sD8x6gq6eXp1qOEY5EcNhszC+28T/XfyHlcSwWC/Vfupsff/NzFBSV8Mu+2wBFHzbmFcMMRx8WpSgoLKJv/i1c85UtLJXAk3XD9nzMRKIDJBSPa1BKbSBuEkK2i8tprX1Kqdq4h7zm48mG0soxek8DKKW8ZiCKf2yX+Wt6aW2FyLCurm6efs7HO+8dIByOcODQYd59/0MevO+zTPMM+m41IpZIP/MuvM60S+9h0WFKuo9Q0nWYozPuoNc5fUzHDgZD+Ha9xuu/3Us4FEZZFO9+cJCa+z7DwgXzBu2f13eORWd3ktd/AZQiZM3jwNwvELQV8bPnXqc7YsVqUeQ7LPx/dffgsA8/aFO+3Mt//K+v8F+3/YrfHD7JLbYPAI1SCnuhB+fqb1B+y90oy5i+g4tRSuddr8a4N9KOUck0cclwLUYKngCwl7h1P9mU0CNah5GLLhkPRlvjtWEMHQoxaRw46OeH2/6dP7zzvjHclJ+HAg4f/Yh//NG/07Lv92ithz1OMgW9p1h+vJHpl94x7q1Y89EoCnrPcPWJJ5kZ2A96dD2sU6fP8i8/fpzXfrOHcDiMK8+F1Wrl3LnzbPu37fh2vU44HDZ21hFmtbey7MST5PedQ2MhaCvkxLRP0F50NS+1vMvRLmOKtMNuo+6OFazwzk+7LS6ng7/7xlo+/sC3OGWZjcViIZDvZfmj/8p1t90rgWcCpZPh4AmzUumOFIlEm4BHMaZc1zL4gz2rzKG28ug9nbHQWq8xj7kLWD3W4wmRjv7+fp7f+Qq/2/tWrHSyy+nkumuX8977H9LT00NPdw+/evZF3v/Qz+fvrqawIL2ElUqHmd2+h1mB/VgiQSIWO0FbHoHCqyjuPoo91Ikt3Mu8C69T3H2UozPuSLs+TiQS4bU3W3hp9xv09fVhtdkodOVx9bIyjp84xcVLl+jt7eXlV97gw0OH+eN7PkFFeD+FPSdQOkLImke/rRD/7M9xxl3OybMX2HUwEJtWfeO8AuoeuGNU7+mdt60ifPOvOHb0EKu9V0nmgUkgnWG3B4DtQwyl7YkLStuUUl8FfpRuA5RS6zGyJaTSrLX2pXs8jJlwQwWeZL2cZL0hIbLu+IlTND39PGfPnQetyc/PZ5rHw3/5oxquW3E1xz46wY9/up2PTpykt6eHd987wPETp7j/7k9zzbKh/jcCV38bi86+RH7fmVhvp99ezME593K+5Hry+s6z7MQTFHcfwx7uorj7GMuP7+D49E/QXrhsyMWT7YFLPPn0c/iPfkQkHMHlclFQUEDtF+7hU7fdTHvgIj95vJF33z9Ab08Pc3oP4dn3NsEZxWh3CUFbEZcKFvHBvAfpdUwjFArx+CvvEdRW7DYLHpeFv/9mzZjeW6vNypKy1Kl2RHalM9vNM8w9nMTKpSP6SqG1zljlU/PeU735uzvFPZ9W4iYcxLXDn2RfIbIiHI7wyhu/5eVXfkMwGMJms+LKy6N85fV8ad0DFBYaPZuFC+bx19//Jk8+/Rwv7XqN3r4+Ll68xOM7fkXljdfzueo1gwuPac2MS39gTtub2MJ9RkkCexHtBUs5MO8B+u1GxuUe53R+v+SrLDi3m4XndmEL92APd7Ho7E5Kuo/y0fRPEba6Eg5tTIZ49vmX6O7pwWKxUFBYgHfxIh7503XMnjkDAE+pm+/8+Vd59eWdWH//S5bmdaEjYU6fC9Deo+lccRcXFnw6loXgid37aA9aUYDDZuV//NGnmFYy+gqlYvJJJ/gMeS8kyVDchOQON6dXN8UFnLWYgdGc4damtQ5orf3xXW5z245st1eIqAtt7TT9+jmOHT9p9BryXBQVFfLQg/dzy03lg4aIHHY7Dz14Hx+7djk/ebyR8xfa6O3p4Xd73+LwkY948P7PsnD+XADsoU4WnttFcfcxlA4RsuYRtBZwZFY1Jz23DOrNaGXl2Mw7aC+8imUnnqCg7yy2cBelHQco6D3NsRm305G/AICu7h5jMsS7HxAKh3E6neS5XNz9mTv5bPXt2BIWaboCH1JbvJ/gxzwcO9FPd1+Ik91W/n1PPm173+Cu++ey/Lob+P2Hx/j9uTARNHkOB/d+bA5VN10/jldATIR0gs+0ER5zpPunzVznE12f04A5JGcGkEbz8ejufi73yhqA5ri/a81ekh9YpbWWNT4i67TWtOz7A88176Kntxer2Wu4+iovD//JOqZPG9RBH2DF8mX8t7/6Dj/f8RR7Wt+it7eXs+cusO0nv2T1J27mgfIZLGp7DXuoC60s9NuK6MybzwfzHqTbNWvIY3fkL2Bf2Z/hPf08s9t/hy3cgzN4kbJTz3Ku5Hp2B2bx5H+8RCBwEZSioKCAObNn8ciX1uFdvHDAsVS4n+KjL1JwpgUVDmLJL8K7fCavX3CzZecpLoZ66e29yFO/+DeuWXkzf4jMIxyx4LDbWFhi42+++vkxv9di8lHDzZZRSn0fOJROpmrz/lCZ1vqxDLVvUlBK7Vo+v3T1xrWVE92USaUj4uKV3qs4Hs5cRiUrEW5zHeRq+xksjG4m18TSOIMdOMOdw+4ZCoXp7OwiEolgs9mw221cVeZlyeIFg3o753sinO7UhCPJ35OOzi7OX2gjHA4TiWicljBLXQFc4W4UEFY23uycxyuX5hMe4fK+q1zt3F16iALVh0WH6LPmczJYyPmgC6UUFouF4uIipnlKsSS5L1QUDlAQMgYktNVBt83Na9NqOMKEVCAAACAASURBVJF3NW0XzvP2/r10dXYSDAbR9jzCziIsSlHksvHL7z/AskVzR9ReMXmsWbOG3bt3745O3oqXzmy3x5RSLyql2rXWL6faz8zrVjfe6XUmilMFmW8dtDzoymbVLLOf4tnu63i++xoiY1uzzGzrJR4uepNF9jbQmpxNkG4DSzCCpbsNxdDTlWcXGR/eLpeNhfOmk+fqg8DB2PaIhjPd0NGvyB/ii2IRMNuj6e3tIxwOo4F+7aLf5qS9s4+fvmPD39kNHBjx6RwCXrPl8dAyxXWz7WCxUuqM4HF2G2WnXU5s1ovQlWwyrKFLWdDKynvW5fyn7R56OvKhowNwMnPFzZw5eZyLbW1EImHQRrbqP//0dRJ4prB00+usBXxKqQsYw1t+Li/UrMSYYu3FHBKbmowPCRFPYwG+UPg2K12n+UnXrZyPjOamsGa18wBfyNuHXYUBBSo332tt/kQcBWh7Hra+i1hCqVf2KwXTps9g9py5WCyWAaGquz/M6c5++kMRI/GlVQ2Zrd+Kkda/r6+P3t4+tNa80edlR3gVkYUBik7vM0oGjPicFD2eq/in/GV8vNfPQ/l7cFlC2G12Y+1RGtOWgxYXLxXew3vOlTiVIrG4tWfZVbS3t3H06FHCoSCfLPPwlfvWjLitInekFXzMm/iV5rTorzEwt5of40b/lOzxxJTMR3/yLye6FZOG0mHmnXmZos7DEOpleSTAp20vc9TzSc4VXZt2TRN7qJOy8z5Kuo9ARIGtAG3P56NZd9KVn/5iwsnC3x7knTO99IYihCIRrBYLZW4b182wY7UMfk8806ZjKXZzNu6xUCRCy5F29n8UIOiMYM+3kOewUr7AzcqF6Q1xdnR288s/BHi/q4jZ/UZOtMLyT3LbfCdL3OmndLzYF2H3sT5OdxrHaLNfw3bHXTy0NMSCeTPTWy+jFP3Fi7jflsf9w+waCt1ET+dFVi1fnHYbRW4aUWJRc1p0dAbZEjO32xXBPW0mt9+9bqKbMblEHoJ3fwW/3wH9ndDfxbzQb8BxCa79AjiH6QWd/gO88zTYLkFJHjiLYVoZ3PIN5pQuys45ZNjtwFsfBfg/rx/mUmc/nX1BTocUB7tcfOHGecx15w16Tm/c72cv9fLU/hOcDLiJWN0U5FnxFDh46NbFrFo89ASERNdURvjVvhM8/85puvpC9AbD/Oa8hZ48N5+5djYue+oqnVprWo8FaD5ymm6tKCy0UOiysXx2EV/55BJmjmulzpGdp8hNQ044UEoVZyJPW6aOM1GUUrtWr169eteuXRPdlMmpzQ9v/BDaD0Nfh3G/Js8N1z0AM1cM3j/YA+89Ayf2QjgEdhc4CuGae+GGh8Bqz/45ZNil3iD//sYRWo6009kfIhiO4LRaWb1sBrctnT6oF6S15s3Dbex87wy9oTA2i4Uip43r55fwyG1LKC1wpHil4b1/+hI/evUwZy720tEXQgHTC53cv3Iui6YNzozQ2RfimbdO8sHpDkKRCHkOK4VOO5+/cS6fu24OliQ9OCGSGWrCQTqz3R7FmNJ8ZKQvbJZcuFNrnXbGg8lIgk8aQv2w/3H44D+hvwuCvWC1wYKbYfndYDNH+dv88PtG6DoPaKN3VDwPbvkzmHXthJ5Cpmmtef3gBR7/7VEu9gTp6gthtSgWegr4/Mq5TCs03pOLPUF+vf8Eh851EY5EyHfYKM6zUVuxgDuvSXNoaxjd/SF+/ttjvHbwPF19IfpDERxWC7eWTWfN1TOwWY17bB+c7uCZ35/kYk8QpaDIaWPhtHwe/aQ3aaASYihjCj4QC0BejDQ7+9PY/0aM5J6/S2eK9mQnwWcETr0Fb/4zdJw2ekEKKJwF1z0I5z6Aw69AqM8IRo5C8K6ByofBMXU/2M529PKvrx7m/dMddPaGCGkjwNy1YjYOm4Xn/nCKS31BrEpR6LLhnV7I+k95kw7RjdWeI238+xtHaOvqp7M/hFUp5pbk8bmPzWHfsQD7jrUbvTS7lQKHleoVs6mpmI/DlpsTQMTEGnPwAVBKlWDMeqvGmNDThjELM4CRBaEMY4FpCcaCzm0pEpHmHAk+I9TXAXt+BEdeN+4FhfrA6gAdNjIlO4qgYBqsehQW3TLRrc2KSETzn2+f4lf7TtLZF6Sn3xhaQ0EoHMFlt1LgtHH39XO4f+XcWE9kPLR39fPj1w/zh+MX6egLEYpEsFsshCMaDRQ6bcwsdvLIbUu4bt6EJCwRU8SY1vlEmYFkG0by0BKMnpAHI/Acxqgc6p8qAUeMgbMIbvs2zKuAlh9Ddxv0dRr3clwlMPdG+PifQf6Vc2PZYlHc87G5XD+vhG2v+jl2oZuOvhBEoDjPztySPL76ySVcNWv885eVFjj4bvUydr53lsa9H3GpJ0R3fwiHzUKB08bNS6bxpVsWUegcdaFjIYY1qn9dZoDZl+G2iKlEKVjyKZhxDbz5j3D6bbC54MY/hmWfSXsq9lSzaFoB/+2ea2nae5zmd08D8KllM3ho1ULyHKlnn2WaUoqqFbNYMbeYra/4OXKhiwKHjT/++EJu8U6TkgNi3MlXGzG+CmfAnf8dLhyCgunGLLgrnMNm4Y9uXkjVNTMJRfS43NtJ11x3Hv/3PSs4fL6LOSUuCqS3I7JE/qWJ8acUTF860a2YdGYWj+damfRZLYqlMwsnuhniCiNTWIQQQmSdBB8hhBBZJ8FHCCFE1knwEUIIkXUSfIQQQmRdxoKPUupOpVRxpo4nhBBi6srYVGut9U6l1I1KqVKgLZ0ccEIIIa5M47HORwFlSql1wB7z73at9Uvj8FpCCCFyUMaCj1JqE7AnLov1E3HbSpRSK6U3JIQQAjI74aAmVfkErfVFCTxCCCGiMhl8nlBKrczg8YQQQkxRGQs+Wuu/AqYppRZn6phCCCGmpkxOtb4DOBQtt62UusN8TAghhBggk7PdPg14lVJe4AJGNdNpgMxyE0IIMUAmg0+z1nonDCi53ZrB4wshhJgiMjnhoDya4cCc3bYtg8cWQggxhWQyw8FjSqkdZq+nFWOB6Spk2E0IIUSCUfV8zEWjgyYTaK3XAn8FtGHc7/nB2JqX9LWrlFJ7kzxerpTaYP40KqVS1mtWSjUopdYrpdzmf6sy3U4hhBCpjbbn4wF8SimN0cvxYdzzeUlrvQ/Yl6kGxjODRBtQnvC4G6jUWm82/64BdgIVQxyuwfzZpLX2jUd7hRBCJDfaez7lGB/sS4GtQBnQpJQKK6X2KKU2KaVuyFQjo7TWPq11skkMlUB93N8+jHtQqXo/e7TWpebP5ky3UwghxNBG2/PRZg8HYJv5g1JqA0avqAx4WSm1R2t919ibOWxjfEqp2riHvObjgaGep5QqTxHMott3mb9K5gYhhMig0fZ8bkr2oNmL2KO1Xqu19mD0hr436taNQEIQWQcM1aPxmkN4fvP+T/kQ+wohhMiw0QafLebw2heSbCuN/mJOt1ajfI1RMYfayrXW9an20VpvNofwAsAWoDHFfmu01msASYoqhBAZNKphN631YaVUHbBDKdWEcY/FjzHk1gb8KG73Q0MdSym1HmOYLpXmEU4IaNBaVw/zmu7okJzW2m9mZRBCCJElo17nYw5zLTWHr6oAN7BDax1fx+dFYC+QtNSCeZyto21DIvOeU735uzvZPR+zvQ0MPRNOCCHEOBrzIlOzV5KqZ9LKMD2fTDGnVzfFBZy1GDPxMHs2bea2FuJmxkWfl402CiGEMIxHGe0Ys8xCxpi9lmrz9wbMITkzuDSaj0d392MGH4yeTjOwVWsdUEq1mb2kAFCmtY6fKSeEEGKcjWvwybS4XlZ9wuN+hpjYkBhczCFDSXoqhBATJJOJRYUQQoi0SPARQgiRdaNNLLp9mO3/opR6IVsLTIUQQuSW0fZ8orPIVkZr+EQppTZhlNO+C9inlHpgjG0UQggxxYx62E0pdRBjivLhhB5Ojdb6MQCzsum0sTVRCCHEVDOmrNZa66Va62nARaXUYnNb4qwz/yhfQwghxBQ12uDj11pfjPt7B2YmaaA9YV89ytcQQggxRY02+HiUUncopYrNHs9fAS1mCe3ShH0lb5oQQogBRhV8zGzVa4EjGJkF/BiZB34A1CmlvqeUWqyU+ioy7CaEECLBWBKLfg34WsLDTwAopdowekMvaq1fGn3zhBBCTEVjSq9j9mwqMEopNGutfwRgVjlNDExCCCEEMLap1i8ClRjDai1ApVlgrnjoZwohhLjSjarnY/Z4ahNmvEWriK4H/i4DbRNCCDFFjbbn054YeADMejmHx9YkIYQQU91og89Qa3dkXY8QQoghjTb4TIvLaBCjlFoJlI2lQUIIIaa+Ud3z0VpvU0rtUEot4fI6Hi9G5oN1GWudEEKIKWks63zWKqVuxJjx5gZ+YE6xFkIIIYY0pnU+ZrAZEHCUUt/TWstsNyGEECkNG3yUUi+M4HgKY9GpBB8hhBAppdPzUUA9EEhz3x+MqUVCCCGmvHSCT/1I7uUoperH0B4hhBBXgGGnWo90EoHWWhaZCiGEGNKoc7sJIYQQoyXBRwghRNZJ8BFCCJF1EnxyWGtrKxUVFdTX19PU1MTmzZvx+Xzj8lo+n4/q6upBj9fX17N58+YxHz8QCNDU1ERTUxP19ZfnrDQ1NeHz+WKvke5+w2lqaqK0tJRAwJjE6ff7qa6uHtP7V1paGrseQ7VrpG0VYiqS4JPDysvL8Xq9rFu3jpqaGjZs2JA0QGRCVVUVbrd70OPr1qWXTampqWnI7Tt27KCtrY2amhoAtm7dGntO9LV9Pl/a+w3H7XbT0NDA1q1bAfB6vdTX11NVVZXW+STT2NjI3r17aWhoAEjartG0VYipaEwZDsRl995777gd+5lnnklrv6amptgHXyAQYOvWrZSXl+P3+/F4PGzatImdO3fS0tJCY2MjtbW1NDQ0UF9fT2trK1VVVZSXlwPGh3plZSUtLS2sX79+0Gtt3ryZ8vJyWltbY4+1trbS1tYW+3skH+Txr+H3+6mrq2PLli2x4Ob1emltbWXDhg1p7ZfOa69fv56ysjI2bNhAIBDA6/Wm3d5kAoEAfr8/dpw9e/YMateFCxdG1VYhppqcCz5KqSqgQWtdkeTx6FfzVcB2rXVr4vPNfb1ADdAKlANbzVpEOamlpYW2tjaam5tjwWfTpk2sW7eO8vLy2Af09u3bY8/ZsmULQOzbvtfrZcuWLZSXl7N58+ZYIGpra2Pr1q0DgkM0qFVVVeHxeGLf3rdv3051dTVVVVX4/X5GIxoovV5vbEgs6sKFCyPebzhVVVWx9icLAn6/P2XvJDEot7W14fF4Yu93snaNpa1CTCU5FXzMANOGETASNQJLtNYBpRTANoxUP8ls0VpXm8f0Aw1AXeZbnB2VlZWUl5cTCATYtGkTDQ0NtLa2Ul1dTWtrK3V1xqnV1dXFAkdUsm/7zc3NsWGtaFCK/6CND3LxNm7cyKZNm6irq6OxsXHAB/fevXsH9IqS9abA6L1FA6Pb7R7wnNHsl4rH4wGM4FtXVzfgPk08r9ebsq2Jovu53W6ampqStms0bRViKsqp4KO19gGYwSVRRVzvxYMRpAYxez2euGP6lVJrGWPwSXdobDy53e7YMNitt94KEAtKYHyzr6+vH/bDNDpU5/V68fv9rFq1asD2VatW0draitfrHfBB6vP5aGhoYOPGjWzdupUNGzbEXqupqSkW0FJpamqKDav5fD5WrVo1aELASPZLxefzxXo50XNI7JFEpdvz2bp1Kx6Ph5qaGqZNmwaQtF2BQGBEbRViqpoyEw601vHjPLUYvZlkykkSmMyglPjYLqXULmBlJtqYadFA4/P5CAQCsQ/UpqYmvvOd79Da2orP56OlpSX2nIaGhtjEgdbWVlpbW2MfsK2trQQCgVjPKfrYhg0bBuy7YcMG/H5/7LHm5mYCgQB79uyJzeQaLtAk8vl81NfXU1FRQUWF0WGtqakZ8OEfHSJLZ79AIBDbnvieRc8vauPGjQN6g/GiPZ9kP/HWrl07YAJBTU1N0nYle0yIK5HSOveqXiultNZ6UPfHDCB1wCGt9dYUz10PVGuta+MeOwTUJt4jMgMPwMrVq1eX7Nq1i1xVX18f+9BN9UE7ntLp+WRa/M1/IUT2rVmzht27d+/WWq9J3DZlej4Q6/1sAiqUUmP+pNNarzHftP1jPdZEi65hmYjAA2Q98AAph9KEEBNvwu/5mD2RsiF2aY7e60mHOeGgEWhWSpUmmcXWxuVZcVEe0isZkbOuxOGdiQq0QojhTXjwSTU8NhJJpl9Hb3IkCyqtxE04iGvH6OYGCyGEGLGpMuzWBmyJ+7sS8EcDilLKq5Ryw+AgY94n2pGthgohhJgEPZ+RMHs40fU5DZhDclrrVqWUxxzCA2N9T/wc1gagGYj2smqVUhsAP7BKa52za3yEECIX5VTwMe/9+DDKeifblup5tQl/+4FoVsehk44JIYTIuKky7CaEECKHSPCZImprawdMLfb5fNTW1g7xjPSN5liZKrUwVqlKQWRKpkszjKRkhJRmELlMgs8U4ff7Y+UBILNTq0dzrHRLLYy3VKUgMiXTpRnSLRkhpRlErsupez6T2Vd+smfcjv2vX1415Pampia2bdvGo48+OqjkQDT1jtvtpqqqKmnZg8TyCdEcbXV1dQPS6SQeK7Fsw/r165OWWogXTY/T0NAQy0AQ/aCOb0e0BET0vBobG2OpfKIf9vE56BLbm4rP5xvwuk1NTSlfJ5q4dDiZLM2QbmkJKc0gcp0Enykg+k05sZ6Mx+OJfSBVVFSwd+/eQWUPUpVPqKurG5CVINmxEss2pCq1EC+6LXqssrIyDh06NKgd0ZIJ5eXlVFZWApczbNfX1w943WTtTcbv98eCZPQ9qKmpiZWSSHydkchkaYbo/kOVjJDSDCLXybDbFLB37158Ph9er3dAxdD44SaPx4Pf72fjxo00NzdTVlZGIBCgubk5tp/X66W5uRkYnB0g2bGivaho2Ybm5ua0vvXHHyuaOTtZO6I9m+rqapqammI9tsTXTWxvfX099fX1A4Yho8ddv379oA/uVK+TjvjSDMnKTCS+9nAJSqOGKxkhpRlErpOeT4YMNzQ2XuI/pLxeL7W1tbFhp/gPWbfbHQtO8WUPhiufEJXsWNEb+dGyDalKLQx1rOhrJ2tHVVVVLKjV1dXFAkzi6yZKFQSimboTP/DXrVuX9HWGMx6lGSC9khFSmkHkOgk+Oay1tTU2ZBQdnolOPFi/fn3sgzwQCMR6CHv2XL43VVNTg9frjc2WipZPiJZSiM+AnexYGzZsGDDTKvp3fKmF9evXD7rhH+21tLS0xAJnQ0PDoHbA5R5NWVlZbFgs8XWjz0mWsTu+FITf78ftdseCXDTTdnl5edLXCQQC3Hnnnezduzfpe9/Q0BAbGoT0SjMMJ3pPbNOmTbH3paamhs2bNw8a1kv2mBC5IidLKmSbUmrX6tWrV+dySYXJpLa2lsbGxoluRlqkLIMQo3fFlFQQk1+0V+X350YeVynLIMT4kGE3kVVVVVUcOnRoopuRNinLIMT4kJ6PEEKIrJPgI4QQIusk+AghhMg6CT5CCCGyToKPEEKIrJPgk8NaW1upqKhg8+bNNDU1xf6brvh9/X7/kGUT4ksTjHeZgnQM197RmuhSFGPV2toaK8kwlPhyDVHptNfn81FaWjogddHmzZupq6sbNC09/t/JeFwvKWeR2yT45LBoZoNogswNGzbw6KOPprU2JZrXLcrr9Q658DO+NMF4lylIx3DtHa2JLkUxlHS+WGzZsoWamppYRodkopm9E6XT3qqqqkGZGsrLy9myZcugfxPx/07G43pJOYvcJut8MuXn41i/5o+2j/gpiaUDEsskVFVV0dLSEstP1traSn19fSwgJT4/lWhJgp07d9LS0kJjY+OgjNDRjNV+v5+amppBJRySHWO4sglVVVWx9kbPrb6+PlZaILo+Z/PmzbHUQF6vd1Dm62TnOZ6lKCorK0ddDmI4W7dupaKiYsjjpJuxoaKigo0bN8bOM75sRl1dHbW1tbEglM6Xnfh/X6muV7ISHcORcha5S4LPFNDS0hIrQ7Bt27bY7/GlA6JJOuM/fOM/cMvLy2PfUpOVHkilpqaG7dsvB8fEwNPU1BRLHLpp06ZY4Egs4RD95uzxeNIumxDfE4t+442WQygvL4+9dk1NzYCEq1GpznM8S1FUVVWNqhxEOqKLd9euXUtdXR0NDQ2DeiPR938oPp+PxsbGWE8ivr1A7AM+vuzDcOL/faW6XoklOtIl5Sxykwy7TQGVlZWx4ZBostBUpQPSker5qURr+SQTX2ahsbExZQmH2tpafD5f2mUTkrU5UXl5OXv27EmZ9TnVeY5nKYrXXnttROcVTRS7detWmpubY7+ner/Lyspwu91UVFQM2ic+C/dQr7dly5bYe5LsOgCxXoHP5xtVFohk1yvVa6Ui5Sxym/R8MmUUQ2PjJVXpgKH2j/8AGenzo99kk+1fVlYW+x81EAikLOGwdu1aHn30UTZu3AgMXzYhXdHhkWRtS3We41mKYqTnFZ8NO5qFO5VVq1YNeK8Tez3RAn/R7OfJsoBXVlbS0NBAbW0tzc3NKdsbHe5Kdd6jMZL3RspZ5D4JPjksmqBz+/bteL3eAcNmiaUD3G530jIJTU1NsXsaQ5UeiI5tR8sqRH+PfoNNNsQDRvmD6MyhQCCQsnRCdMgt2rbhyibEtze+PdHEpdEPh02bNuHxeGhraxtU8iBViYXxLEVRVVWVdjmIkUosvRD9YK2uro4V5wPj3lCyD2qfzxcbwnW73bEeRXx743tO0eHcVFJdo1TXK/GaV1ZWJi1pIeUspgittfwM8wPsWr16tRaDbdiwQWut9d69eye4JYNF25bq71zT2Ng40U3IukOHDk10E8QYrF69WgO7dJLPVen5iDGJrquYjN/+1q1bF+u1+f3+2BBcrsrEpIRcIyUtpi4JPmJMJmPQiYqvUCqlEXKTXLepS2a7CSGEyDoJPkIIIbJOgo8QQoisk+AjhBAi63JuwoFSqgpo0FpXJHk8utBkFbBda92a4hgNwCFgB7AW8GutJTugEEJkSU71fMwA0wYkmwLTCPi01k3AHmDbMIdrAA4Dbgk8QgiRXTnV84kGCaVUss0VWuvoogAPRpBKZY/WenBBEyGEEFmRU8FnKFrr+AIltRg9myEppcpTDc2Z23eZv96yf/9+1qxZM6Y2CiHElWT//v0AS5NtU9pIH5NTlFJaaz2o+6OU8gJ1wCGtdfK0v8Z+G4BWoAXYSIr7Q3HB5xNAGPjN2Fs/qa00/7t/Qlsxvq6EcwQ5z6kmV89zKXBOa31j4oYpFXzMbW6MXk+zef9nuGN5zX3LhthnF4DWes2oGpwjroTzvBLOEeQ8p5qpeJ4TPuymlFoPpPzgxwgMaU8I0FoHlFKNQLNSqjTuPlD8a7qjj2ut/WYAEkIIkSUTHnyGGh5LV5Lp1y3mfz1AINm+wICp2kIIIbJnSgy7KaXKgcpoIDMDzJboUJrZs2kze0Vuc9/ozLkaYJ3WujbrJyKEEFeoCe/5jIQZVKrN36P3dXxa61allMccwgOjVxNfXrABaAa2mgGozZx0EADKJPAIIUR25WTPRwghRG7LqQwH2aCUqlJK7R1mny3DbPcqpTaYx9pgDvVNKqnO03y8xvxpMIc0Ux2jQSm1XinlNv876Yr7ZOg8c/l6lptt3qCUahyq7ZP9emboHHP2Wg63LWG/SX0tASmjHf8DVGGk7tFD7DPkdnOf5rjfvRj3nyb8/NI5T6AdI+UQQA2wd4jjNJj7twMbJvq8xvE8c/J6YuQ6XB/3d85ezwyeY05ey+G25dK1jLVxohswGX+GCT41QPsQ272J//iH2n+ynSfgjft9ffz/rMnei4k+h/E+z1y+nuaH1aG4v92AjgbdXLyeYznHXL6W6W7LpWspw24joJSq0cMvXC0nSV65XFlLpEeZpmj8WjQ+RnCeOXs9tTGjM34yjdd8fNDat3i5dD1HeI45ey1HazJfSwk+aTL/gfqH3THJ2iKMf/CTbmw5FXNcvAFo1EMv8PWaY8n+4e6bTEZpnmdOX089MG3UOmDzELvn5PUcwTnm9LUcoUl/LSX4pG/IJKRTidkr2ARUmOugUu23WRtT3QPAFoyyFjkj3fOcCswb6+V6iGzuuX490znHK0UuXMucWuczGplI32N+g0g3xU+yb1LJvnFl1JWSpmgCznNKXE+MDCDVQ+2Q7es5Aec4Va5lOq85+VOITfRNp8n4Q/KbmuvjfrT5X2+S5+bsTU3zPPfG/R29eZvsPKuSnKfOVNsm0Xnm7PWMe3wDl2f2pZpskBPXc4znmPPXMp3rkivXUobd0qCN7uvW6I/52FZt3rQ27x24zccH3Bcyv3HsyHqjR6cNo4seVYlRYnzQeWLkz4sNb5jDVsNmEZ8k0j7PHL+eseuiL/fo1sZtmxLXM91zzPVrOZRcvJZTfthtJFSK9D1x290YPZ7o9i3mP+hY+h5z11plpO/xA6u01nXZO4vhpTpPPcXSFGXiPM2/c/J6mh+ujebj0d39XD6vnLmemThH8++cvJbDbSOHrmWUpNcRQgiRdTLsJoQQIusk+AghhMg6CT5CXMEmY2JNcWWQ4CNECma25GalVPNEt2U8KKXW6+FT7axXSu1VSrWbN7BT7aPjM0pHMyqPR7vF1CATDoQYgpmWZNgFmrnGnOnnS5x+nGLfcox1I2qIfRp0QmYBpdQGrfVQ6XzEFUx6PkIMbVxXv08Es0dSlk7ggVjuNH+qFETmdOdkvcPWSVlHRkwKEnyEuPKsBbaP8DlNQKo1MVU6SXoY87FJtY5GTB4SfIS48lTrkSfJ3QJUjeI+TtukzCsmJpxkOBBihMz7JdG6MJ5oyqUk271A9EO+Ltkqc3NfN8bwjhuHiAAAA3RJREFUng8jLxcYK9mj91DKMRJgVgOPDjdJYJi2D1VieoPZXjcJ56WN5JStGBk+Nsc9p5yhk+7uNduf1hCfuHJIz0eIETAzYPu01k3aKCy4I342nFJqC0aeuCaMlC5bzJQ+yQJPtEaUDyM9Snlc/sA9GL0Nr/laWzEC2vrE44yQl+QF1aLn5TPbXpbkHs8WBg+jVQ5z78gPrBpLg8XUJMFHiDSZ3/K98R+2Zi/EH5cnbi1GYsfotqGqSXrN+yLRxKbxyR/LzGPE9yq8jH0ChAc4FP+AGQQT61VtZ3Cg2YFRpCz+fAYFsgR+pmaxNjFGMuwmRPoqST58dAgjOSnm9vgaMZ4Uz4kPLNUMngBQiVHoLvGxsRZKiw7xxasCAgkz09yYJamjzISVTRjVQlvNntFwdWjaMN4DIQaQno8Q6RvqG3z0A3YLUKeUcpu9oU1p3KNJVqxwwL2UaGAYxUSBRAEGn4cbo+fli/tp0lonK4C2hctDf540zs3D8L0jcQWS4CNE+nwk9AZMZVxe59KG0WOpxLiHMuQiS3PIyx0fVMxA40/4YK/FLAsQ30NRSm1QSlUppRrSXFPTxuCqmq0pzmuQuPT+60lvEkEmhgrFFCTBR4g0xS22jN3zMGePVcbNDKs2900rewDJez3VSR6rxCj17cbsuZiTBLaaASE6eWE40WHB+PPyYUyJHnBvaohgthUj60M6pZ+9GJMnhBhA7vkIkYLZK2kAKqOpYrTWtWZvI9pT8AJ3xj2tEdirlIof3mpKTD0Tp4yBVVWjx0x8bAvGUJxXa73VDBTuuN6RN51gZ963GTR8qLWuNs+rEnOYLGECRGJb0lXB2O9TiSlIcrsJkSFmQKhKHGoz18+UZbJqpnnMaVrrejOYNKabf858ri8D94/Sea3GyVhFU0w8GXYTInPqMNLQDGAGo0yv8g9weThrLcaQXLp51LZizFgbV2Z7RtJLElcQCT5CZE4zlzMUxJgfwpkuy7ADWBWdnABMS/eJ5lDdniykvSlP876QuALJsJsQGWQGg/ib/16gbYj7JxPGrOezdfg9R3XsGoyhPZnpJpKS4CPEFUwp5R6PADFexxVThwQfIYQQWSf3fIQQQmSdBB8hhBBZJ8FHCCFE1knwEUIIkXUSfIQQQmSdBB8hhBBZJ8FHCCFE1knwEUIIkXX/PzWic2XvMVFQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# set up some labels and colours to use to match those in the paper\n",
    "labels = [r\"Reynolds+ 2020, $N=500$\", \"Absorbed power-law, $N=200$\", \"Partially covering absorber $+$ $6.4$~keV line, $N=200$\"]\n",
    "mods = [\"r20\", \"abspl\", \"1\"]\n",
    "colors = [\"k\", \"C0\", \"C1\"]\n",
    "for i, label in enumerate(labels):\n",
    "    plt.plot(data[\"logm\"], data[\"logg_{}\".format(mods[i])], c=colors[i], lw=3, ls=\"-\", label=label, alpha=0.7)\n",
    "    plt.fill_between(data[\"logm\"], data[\"logg_{}\".format(mods[i])], y2=-10.8, alpha=0.2, color=colors[i])\n",
    "\n",
    "# aesthetics\n",
    "plt.legend(frameon=False, loc=3)\n",
    "plt.xlim(-14,-11.1)\n",
    "plt.ylim(-13.5,-10.8)\n",
    "plt.ylabel(r\"$\\log g_{a\\gamma}~({\\rm GeV}^{-1})$\", fontsize=18)\n",
    "_ = plt.xlabel(r\"$\\log m_a~({\\rm eV})$\", fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fig 8: B-field Model Comparison\n",
    "\n",
    "**Figure caption:**\n",
    "     \n",
    "> How does the choice of magnetic field model affect the NGC 1275 X-ray limits on light axion-like particles? A comparison of $99.7$ per cent ALP limits obtained using different approaches for modelling the magnetic field. The models used are given in Table 1 and the colours of the lines match those in the table and in Fig. 5. The limits typically vary by 0.1 dex at $m_a<10^{-14}$~eV, with the most pessimistic model (model 3) resulting in weaker constraints on $g_{a \\gamma}$ by 0.3 dex. Overall, the limits obtained are not very sensitive to whether a cell-based or GRF approach is used, and the choices made about the coherence length of the magnetic field only change the limits by $0.1$ dex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAERCAYAAACkWKo8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXSb93kn+u8DgCC4SARBSbYsWQuoxY5lW4aoxHbSSLbBbG42m5RmMkmaaR3SmU6bpjeHtHo725npYcDk9N7OnN4p6fbetvemE4lI0uyJAdmyk7i2ScJOvEoWod2iFoIgJS4gCDz3j3fBix0gQXB7Pj48AvCugCw8/G3PQ8wMIYQQopxMi30DQgghVh8JPkIIIcpOgo8QQoiyk+AjhBCi7CT4CCGEKDvLYt/AckBErwJYD+D0Yt+LEEIsIzsAXGPm+1I3SPApzPoqW+Wm2zdu3LTYNyKEEGXFhHhkbdrLRDFQ5c2ch164fBlT05GM2yT4FOa0c8vtm/5bx58s9n0IIURZRScrcfVNJwDAYotgNlIBsDJic+veUzBXxLIe+9RffAsnz5zN2GMkYz5CCCGyis+a9ccW2wxsdRP688h4zZzPK8FHCCFEVsbgY6qYRVX9Df25BB8hhBALImYIPuaKGKocxuBTjblmaJPgI4QQIqt4UvCZReXaCZBZGeeJzVgRi1TM6bwSfIQQQmSV2u1GJqDKnpjlNteuNwk+Qgghskpt+QCAzTH/cR8JPkIIIbKKzyZW5JjUadXGSQczk7Y5nVeCjxBCiKwytXysNdMAKTMNYhEr4rHiQ4kEHyGEEFmljvkAAJkYFdXT+uuzU9aizyvBRwghRFZJLR/rrP7YWjOlP45OFd/1tuyCDxG5iWiw2G0p+zmJqEPdv4OI7KW/UyGEWN6YAY5pwYdhsiRS6VhrEy2f6GRl0edeVrndiMgNIATAVcy2DHqYuVk9LgjAA6C9hLcqhBDLXlKXmyUGosQ2a62x5VN88FlWLR9m9jNzoNhtRkTkBOAwHBcEcKh0dymEECtDpskGGmPwmZ2qLDrTwbJq+ZSIC0oLKQkROdVAZHzthPpwbxnuSwghlpSklo81OfhYbDMgcwwcMyM+a0F81pwzw3WqZdXyKREHgHDKayEAMu4jhBAG8ZS8bkZE6pRrVbTI9T6rMfgUjJkPMvNBAK8t9r0IIUS5ZZpmbTSfcZ9F73YjojYAjTl28TGzv4SXzNTKydQaEkKIVS3XmA+QPu5TjEUPPszcW+ZLBmCYcGC4j2CGfYUQYtXK2/IxrvUpcrr1quh2U9f12IH0IKPOfju2KDcmhBBLmDGvW76WT7TIGW/LKvioi0I96mOPurYn7zYo63iM06lb1cWlLQDamVnW+AghRIrklk/6TDazNQazNao8YRNmI4Wn2Vn0brdiqGM/fgCdRW5rTXkeBNCtPvWW/k6FEGL5yzfmAyitn6mQUlBudrISFbaZgs69rFo+QgghyqfQ4KMpZsabBB8hhBAZ5ZtwAKQmGJXgI4QQYh6Ycy8y1RgTjBYz3VqCjxBCiDQcMwFQMomSOQYyZZ7KVlEzBUDZNjttBccp436pJPgIIYRIU8h4DwCYzIyKqoj6jAruepPgI4QQIk0h4z2aijlMOpDgI4QQIk0h4z2apHGfAjMdLKt1PkIIIcqjmJaPcbr1zM0qRG5UAQA4bs52iAQfIYQQ6Qod8wFSgs9ENa6/s005x0xt1mOk200IIUSafHndjCqqIjBVRIs6v7R8hBBCpMmX182ICNhw1zmMnd+AeCzRpiFz9oAkwUcIIUSaYrrdAKC6YRzVDeNJr5m/n71MmnS7CSGESFPMhIO5kOAjhBAiTazIlk+xJPgIIYRIU8yYz1zImI8QQog0xY75MDNCfANRJALVDLIHLQk+QgghksRjBLDSMUamOMgcz7k/M+PbM8/jt7FzSa8Pm25kPUa63YQQQiRJnWxAeRJVvxt/Ly3w5CMtHyGEEEmKyevGzPhZ9FX9+VpzLarMNgDABTqf9TgJPkIIIZIktXwsucd7Xo+dw6X4CADAQhb88ZbfQ51lDQDgaxX/AWMYz3icdLsJIYRIktTysWYPPjGO45noa/rzD9a59MCTjwQfIYQQSYx53XItMA3EhnCVxwAANpMVDzkeKPgaEnyEEEIkKWTMZ5Zj8EV/oz//sP0DqDFXFXwNCT5CCCGSFJJa56XZkwjzBACgxlyN36lvSto+PQuMRyqyXkOCjxBCiCT5FphGOIrj0df154/UPwCbKbmC6bPBSoxMWbNeQ4KPEEKIJPmCz69m38IElNLZdsta3F93X9L2iRmC77Qt5zUk+AghhEgSm0l0l5krk2vyTPA0TkTf1J83Oz6IClPyqp1fnK7E1GzulakSfIQQQiSJRRPBxJISfE5E30AEymvrKxzYt/bupO3hacJzweQuuEwk+AghhNDFYwSOqd1uFE+acDAWn8SvZ9/Rn3+04cMwU3IY+ekpG6JxpdVjs3DW60jwEUIIoYtFE11ulspoUl6347O/wayaqXpT5S24u3Z30rHXJkz41bnEJIP1NSUso01E2wC4AOwHYAfgABACEAbQD8DPzJnzKZQAEbkBeJh5XzHbUvbzABgCcAzAIQBBZvYvxP0KIcRyEp9JhAXjeM9IfByvzL6rP/94wwGYUjKO/ugdG+KsvLbnljguVWTPhl1w8CGiRwC0AxgBEADgBxCEEngcAJzqTzcR1QPoYeZnCz1/gffgVq/nKmZbFh71p0sCjxBitZtBFDHEEYlW669R5TQmOQIA+EX0NcShdKM5q27HrurtCN+IYFaNL9cmTHj5zCwYMVCFDW2bLuK/TYxlvV7e4ENEdQDaAAwy86Esu40BOAPgOICn1eMeIaKvA+gtVUtICxKUIb93rm0Z9DNzZ76diOiE+nBvwTcphBDLjN/0Cl6gV8HE2Dv7MO7HVgDAK+Y38OLU99P2f2jth/CH/+mHGD57NuP5PvHHR3DXD48CF25mvWYhYz6HmPmbxbZimPk4M38LQDMRrS3m2HIhokJbSUIIsSINYwQvUABMSqumZsaub5uwhtP2v7OmEa+8ZMoaeADgS3tugs5eynndvMGHmZ/Ot0+e47+7kGNAc+RUu+mCROTJFoSY+SAzHwTwWqbtQgix3PlNr0AdpoEFFqyZcejbopWTqDJV6j+bK2/Fx+zNODGQCCzmCgsqKm2oqLSh0maD85ZKuK78DDPh3OFlVdbzYeZu7TER9QDwAWhcvDsSQojyO49hnDQpFUgJwJNr/y1M8d2YUbf/q+0fRdUtDyUd88N3bLh5+WcAgEoL8Nf/5bO4957E16d5ehz1/+spjM3W57z2ogcfImpD7i9+X6knBBCRnZnDAMDMQSJylvL8Qgix1DEYPvPL+vO7rXfhVsstuDydaLFYqpMzWo9HCM+8GcXs2FUAwJZ6E+7YfXvSPuvf+glmR7Kv79HPPZebJqInmPlv53JsKmbuLcV5CqVNxwaQczq2EEKsZKfpIs7SZQCACSY8UnUAzEAskgg+5qrk4PPzd224OaxMt66qAPbdtRGVlYl1PRUTI6g//SxGwrnzugFZgg8RfQNAXZZjCMoXd0mCTzmoLZuQ2toZANBp2NYCwLtY9yaEEOXGYPhNiVZPU+VeNJgdiE0TtAEgkzUOY8q20CThhbNWzFwJAgA21zHuvHN70nnXv/kDmOKziIxlL6WgydbyGVL/DGbZ7sjy+oJSWy3N6mMPDF1yubZBaen4oEz7DhNRiIg6oCyMbWTm1jK/FSGEAKAEAkJBS0TmJY7Egs+36Azeo+sAlEkGB6s+BACIRRLZrC3VMcQNvWc/PmVDNAbMXD2DWitQX8XYfcc2fbt1/DLsZ34FAIiM5e9Uy7bHMQBNzHw808YC19KUnBpM/DC0XArc1pryPABloawQQiya1+gUfmZ6Ebt5Kz4bP7ggQegGJvCP5p9imEYybr/fth9rTcpqmJhhvOfdScJf/MietG9sIoT4ZBib1zMqK63YunWjvm3D698HMSMeA2ZuaC2f7O8n41w4Zh7LFnjU7Vm3CSGEyI/BeMb0EiZpGq+aTuISri7IdXymV7IGnkpU4sO2B/TnxuCTKTdB9OoZrLUBa22MnTu3wGJRWkqmmUmsvTgIAJgZt0BNhACyZc9uLYlFhRBiEVzDKG7QpP58iHIvypyLqxjFa3RKf06G/6qoCp+s+RiqTYl0Osbgc1NddGoi1n9s40FssyuvG7vcaq6+A2KlW+9GfIv+uqmqKuu9FZPb7eFS52oTQojVKpgSbIJ0CQe4tElXjpv69cwFOyxOfGnt53LuPzOVGPO5aWJ8ef8sPnOXMuONmdH5+hAmKjIEnytv648nZtYDuAAAMFUnAluqYlo+kopGCCFKJLWlc56GEUV6yeq5uoSreMuUmDPWXP1Qjr0VI2OJ4GOqiuMTuxNTrS9evIKJiSkAQO2aGtx223p9W82Vt/THkbDFcI7sLZ9igs/izDIQQogVJoY4ztB7Sa/NUgznabhk1/CZXtEf31VxJzZZNubYG5iIAhFDy+fA+2ZhNfSNnXznnP549+6t+sQzy9QobOPKe4mbKhC7lMgHRyUKPvmXrAohhMjrPVxDhGbSXi/VuE+QLmHIdBGAMs7zSPWBvMccv2RFTTzRxvjQHckLTE+ePKs/vuOOxPqepC63+rtA10eV61qtMNmyLzaVCQdCCFFmxvGeWqo1vH5x3udWFpAmWj33We/BBvO6nMeMRQgvvlcBm5ZhlBi2mkR7Y3Y2htPvntefZxvvGTMlMpVZGxuBHMtyFj23mxBCrDZDhiBzoOpB/HTSBwbjPVzDFCKoQvYpyqkG6R34TC9jBkrVUQZjlpRWixlmPFz14bzn+Pvjb+Hqr1/CwIYmNO14GJW1nBQ3zgQvYWZGOf+6dfVoaFAT4DCjdvhNfb+pqUT+AduuncDLiSwKqaTlI4QQZRTFLC7QFf35XdY7sMmsjMcwIW0sKJdx3MSPTb/EBE0hSrOI0qweeADg/ZX7YDdny5SmuHZ9FK8c/zl4YhQD7x7H5dBZVNYml79+8cVEVZndd2zVH1tvXEHFlNLNFrPWIHot0ZVYuWtXzusWE3yy10MVQghRkPM0rAeI9aYGrDWthbNim759qIiutxOmQFKwMdpovhUPqWlzcvmB70XE43FY1F62l089A2tN4pyXLl1F/yuJ1s0H7r9Hf2yc5Tax6T7gTOLe8wWfgrvd5ltUTgghRHJwaazYrv/5wvSLANLX/2QzgjEM0jv68y/UHsZ2QxCzUv7knleuXMVvX38bAMGiTmgeHj2Hi6FT2AOlVMKPf/gCmJXItGfPDjQ2btaPrzUEn5ub9wHBn+jP593yWaolsIUQYjkyBhenGny2WG6HRW0LXKcwxnAz73meNfUjTkr32HbLVuyq2AErVeg/hTjx3C8xEVWCTgUlplm/MHgczIxg8CJ++9tEhoRPftoway4eR/XVRPCbsO0A3VQyNpjq6mDZsCHntQvpdpMWjxBClMAUIngP1wAoU6C3W5TxkwqyYIsl0aLI1/oZxghep9P68+aqh4pO+Hzx4iWcPDWESXVd64P7Pw+zWQlaV0PDGBx4Cz/8wfP6/vua3ofNm2/Rn9vC52GZmQAARKsbMGMoIGfbuTPv/RQSfFqJqIuIHi7wPQkhhMggSJe0cjnYZN6IKlNiHYzWBQfkX+/jN72in+eOip3YUrE55/6pmBnPPfuC3uqpvH0P7tywA3dvVZKMmiyMo9/5Bd49pSwsNZlM+N1PJs+aq72SGAeauH0fELygP8/X5QYUFnzamfkIgDNE9GUieoKIthVwnBBCCIPkLrdtSdsaDc+DdBGcZV3/eQzjpEkJCgTAXZU/bU6qM8GzOHfuAiZmCTCZUXPXQ1jLhPucH4a1wgaTGZicnNb3f/DBe7FhQ3IZt5phw+LSzU2goUQGhMpdO/PeQ94JB9pEA2Y+A7ULjojuI6JmACMA/Mw8nvdKy9x5ugKP+R8W+zaEEIvEzFGYOT0rQX6EqMkGhglTiOivGls6AHCbeSNsVIlpjuAGTcJj/sf0nGYM3HXpAB4f/Tr6b/8p3md9Df/mXDcicUKFiWHN0JyIzDL+j1/dwNshE0Zqd2DWbEMkoryPm1FC1ba9sNQ4UDFOQEUV9m7/ME5HExMHLBYLPv6J5FlzFIui+vq7+vOJ25uAYCLvdCEtnzktMmXmVwG8CgBE9DgROQAMreSs1zHEcJOmFvs2hBCLhYC5L42MJD2zwIItltuTXjORCdst2/B29CQAYCLD983msd247+LHAQDud38Pn9rwEmp5HLUAEFd/UjxzMorBC8oCUdP0aUzUKlkIZuPADCrguPPDWAPWy2fft+sDuHrhBG7cUMZzDhzcB3v9mqRzVl0/DVNMCWAR++2IhmKgM4Zut535Wz7zXmTKzN9VW0eNRDRCRF3zPacQQqxUJhAOVH0QFZT+u//v2B6AjbLkQ2Pg/ed/V39qjdkwPPbpvNd79XIiIlXGbsI6q8ykm4yZULPnYZir1uJ9NYk2Vq29Av/m859ARYUFt99+Kz76sQfTzmmcYj2xuQn0t0dBas3tmgcfhHnNmrRjUs0rvQ4R7QVwGEAblN8LetWfFWeLZTM67X+y2LchhFgE77v0XWwZ+RcAwOmNj2LotkcLOs55+afY+d6PAQCTFQ78cncHLOYq2Chz+pwtFZvxlP1rmOL0Vk/10DVcntiS9NqrU5/Gn1kewKRJCTB3brCg5a4afXssFsMzx78Jm3UI1TPX8ZcfqYSlwYIXd30F37+wFgOjyr67DWV3rLVx3H3PTvyf/70j6/tKyucW2wp69pj+fP2ffDXrcUZFBx913U8bgCcBbAfgBXBopZfWNpMZa0y1+XcUQqw4u8ZOozamfMG/0fBBWKu35jlCMbz1i9h3+VmltRG7jj2jv8HZ9bkzTFvIjDWU8l0Ti2E8mAg8hDgYJhATdk03wFetdKtNjBIetdbBpE5zvnjuLKZjQLR6M7ZWhLCzwQzgIq7MDiE4mWjRbLGS3jGYmlonlWlmElUhpU4QgzD940QrqNb9CKruuSfbocnnybcDEX1d/fMxIvoFgFEorR0PAAczH17pgUcIsXrZZkKojSi52GImK0bW3lnwsbOWGryz5ZD+fPflH8Mcj+Q4IrPKU+9hfPZWAEAFTWLdI4ltd8+YYY8pwWYyyrh8I5Ea5+xpZVJAzFyJ23bfm7iPS/+MsHobNjOjzjC1obI2d/WcmmsnQWrGg3C0Eeh/Q9lAhA1fLazVAxQ25vNnRBQD0A1lksEOZt7PzE8zs+R7E0KsaOvHE6v4r9XdhbjZWtTxpzd9ElNWZZqyLTqG7VefK+p4ikYxfC7R0tq8/Sys76tFfIMSMMwgfHA60Yn17kiiGuq5ocRC1OqmFkTNSnG3+sh7+LTp1wCA9zkYsalEKKisyd3y0fK5MQPXBxOZFOo+9amCJhpoCgk+QQBNzLyDmZ9Sp1wLIcSqsP5GYnzjav3eoo+PmW14a+vn9Oc7h38Oy+xkwcdb3x7GRKwBAGAzjWP24btARAjelmihvC9qwXq19XN6ROmCm5mJ4NL5xNqb23bfi1ObP6s//5rFiwrMYk8DIzqVaPlY83S71Qwrwefm5UrMnlXbHxYL1v3Rvy/4PQGFjfn0qFOrhRBidWHG+nFD8LEXH3wA4MzGj2L3BS/Ohu7Fu1O/g7sGv4fbat5J2292Ahh9yQSzDbB/II6YqRLfu/Rf9e2bdr2HSLVSfuE3k1GwxQTr5dMInD4B1/od2Lf5XmAYGH9zHNcvvoiqt99S3kNFJWq+9X/hf3//p/D31T+EHTexxXQNh83PYc+6h3Dj7eRuN+v4ZWwMfBuRtZswvPcQYFLyvlmmwrCNvwdm4NpvE6Ua6g8dgnVzcVkW8rZ8smWzVtf3PENE7xpekxQ8QogVY830ZdhmlTX0M5ZajK5pnNN52GTBKw1P4vnxJ/Fe9C48N/IkTJMR1EauJv2M/2oW0xdNmDhtwvSrUxga/QCm48qXfI1lBDMH9wAAIlHG2WuzCJjH8MIb38eV8Hm8+O6zuDA1gahtHaYrG/B61IKIaQ3MMxHsJcae4VP4zMs/wV9HP6Xf11ct38dGzCByw9jtNovNL/5P1A6/iYZTz8AxdELfps1yGz9vQySstF2oqgrrvvJk0Z/JnNb5ENGXAewH0AngKe11Zn6WiB6byzmFEGKpSepys98LGDI/F+v0xUSraZarMDjRkrR94ooVE8OJNT5Xghvw6s3PJO7l3kmwVZmiffpKFHEGmt74Z8RHEyMhL5/y6eUPLo4MYapqPeJkxu5K5bgPXHkbr1zdjcusjEGtpzDG+kfAcaXlU7thFuuuv4yqcGLB6Lo3fwiaVWYn1Fx5CxwHrr2eKHbg+MIXYFm/vujPY66LTEPq+M+rUGa/GRWXWlUIIZaodYbJBlfmMN6jmbkSQ/RK8iyy304/ih/u+Tv87P1P42f7e3EqmHz+s7d9FDNQ1uFU2BnRB3br204Nz8AxNYatp36JmsnLcIy+g5rqKbw2+jb+afRfMDo7jpHxy2Ay4/y6OzC8vUk/9j8On8DlO74AABibvQUXrjj1bdvvn8CG17+fdB8V0+NoOOUHmFFz5S2Eg9WI3lRaPaa1a9HwB78/p89krsHHGHBSg039HM8phBBLBnEM626c1J/PZbIBoGSQvvnb2fQNccK1ofW4Wb0ZNDQC2zmlyBwTIWKtw8VNifVAax6wgUyJr9p3L0fxuZN+DE0pExcsNTbUbbsNNy0VeCPYj5/cSORd2163Dd/f82mwOm5T/e5JxCPbMFWzCa/c/FeIq0P/dZtm0Rh7DpU3ryr3YfhqX/f2T1EVOgPLeAjX30xkL2j4gz+AuS53me5s5hp89qnZDZR7VKmvSfARQix79olzqIgrmZ0nK9fhZtWmOZ1n5r04ZrVaNyagvjkxVXv6bAyz4Rhu+UGitRFp/gyC93xOn9JdjRCqdiS6+25MxRG/eAXusy/h1IzSHWZ1OlFVaUWtrQKz4ct44x0/ZtRv9y0NjXj/mtth+0QiFY/taB/euu0rODWdKJOw675hbHjrB/rzq/c/gUidMonAHJ3E5hf/BqOnazA7pdyLed06OL7w+Tl9JupHUTxm/iaU9T/vAvAQ0VEi6gfgUbctGCJyE9FghtddRNSh/vQRkT3HOZzqfm71z6z7CiFWJ+N4z5X6+4Aii7UBAMeTWz2191hQc2cFKreowYSBmReuwXZJbfVU2jD28GEMr7lLP2bn6/+EykuX9efvDkfx+befwaXINCLxOMx2O27dtQuf+cxn4KhW1t3EJkYxoqbc2dTQiC1jcfAnvwioYz/mM2dx8jdboIWArZUD2PvWf0LFVBiAUhxuZO+/xtX7v6xf1xy+jpG3EpkX1n3lSZiqDXl5ijTnxKLMfAjAR6DkchsA8BQzf3TOd1IAInIDCAFwpbxuh7IWqZuZuwEcBZAr60KPuq8fSnogz0LdsxBieVpXginWkfNxxMaUVg9VAGualNZM3YOJxZk3J+0Yr1UyXE9/rAXvvRsFqxmm7aOn4Bh5Gw3/9D19//DrZ3Hw4qs4GUm0eu699148/vjj2Lgu0QU2SYyoxYJb7FuAOPD2G9OwffowAGBs7XaMjie6z+6v/TaskyH9+bX9XwJX2DC+8xFMrVMWjobeqUVsRgmaFbduQH1r65w+E828Eosy8xki8jHz2XndReHX8wPIVJ61CcrMOy2pqR9AHxHZmTls3JGInAD0qkjMHCSiQwDac107di2OMa+UVBBitfg5/zv9ceysFaCzBR0Xr7Ai6nAgDsZLg36Ex6/jgTs+jtv2r0f9iy+gzv8CaHYWs+s+jVC1Mongt3d/BdbZG4iZtmNqKFEebWvwn/G9sTDCP/kFfvfCJdRVWnH48ihMYLwTiSBud8BcV4e9e/dizZo1+NLnDuOV//JX+oy3yVu2wKyO9Zx7YwRj6z6M2P5NmK5IzFa7dWwQN35zEzewTrl/sxWRwEsgelk5LlKDyrF1iIwnAua6r34NZC0u00OqvMFHTSSaJKV4nIeIOmCYhFDu4nLM7CciYxh2qq+HM+zugtJ6SkJETmYOprx2Qn24d9Zkw42qLamHCSFEujDw5vlXEDh1AgAwMTuOP/lsK2751t/redF2Dh/Fy/v/HCATZirrMFNZB4xG9VM4dtZhYAg4/vYNAMDkb36LP1q3HjYAM/E4hmYiqL77bgDAPWoyz8Mtj+E//o9/wLUR5ev47vtd2LLNgfNvhQAGwtciQE1iMShxDNve+RGmp5IDCV1PTN9mANNIbK/cUI26T31y3h9RId1uzQACUILL0wDcKdtboaTgCQMYhPrFX27MHDA8PQwlF10mDij3ahQCIOM+QoiSiMZmMDiUyOEWwkVc/5sePfAAQM3kFWy5kLn+prnSjI0P1uMZkwlsUVocb09P41REmQBxemYGYft6WNbU4vbbb4fDoeaOs9nwlXalE8dWXYtv/OFhPPj4DlTYMq9P2nLej+qpa4W/MWLc8qf/DmSe+3onTSFltL+rVio9liWRqBfAl6FMuW5F+hd7WanjPy5mbp7vuZj5oHrOE2vXzhz4wP1F/CUJIZa9WMUaxCoL/73U/O3/B7hwDs+OXoEpPgnrpltAJsAyPQHfy2/h/RtuAUwmTPx5F/7rS9cwGY3BhnMYqVmPOJnwxIe2465Ndahy2PDzn3gxEY8Dribw1BRuRmbxx7E12HH/47gwNIDbzZcAAHv3Jo9HPfX7j+Ex9wNoqFuDhjplgsCXuj6I8NVEPjmOzoKuXUJt7RdA9EUADEyMAFV2wJQlLERuoOKWdbBs3lHch5hFId1ujwE4mqMrrd8QlJ4moicA/G2hN0BEbQBy5azwaWM9BfLkCTyZWjmZWkNJTPUOVLo/VsRtCCFWGzI9jshTX8Wvrl1GJTPit9YjbrHCeuky3otGMTA1hX2/+xl4K7fjtzWpX7+Mb787jL/atxETN8bx3DM/VV42m4GaWpwfG8fEzBhmpiax2TaNarUS6r333otUu7ZsTHpurbJgw9aUEZQdDiymQiYcOPKM4aRWLi1qPiIzl6zyqdfKmPsAACAASURBVDr21Kk+TptsoArAMOHAcB/BDPsKIUTBeP/98NU7MPneJRCAjeEx7PzwgxgYeA0A8KMbN7Dl8S/gB89d0o/580fvxF/6TmFyJoazI5N44dQ1XB/8GSJqF9umzVtA9VvQf/aHAICpN/xw3GICLAQiwp49e8r+PkuhkDGfnG3ODF1xc1vuOk9E1ALAawg4hwzbnNpangyTCpwAjkEIIeZpfHwMz1oTJbI/OzmNx988hWqT8lV7zdGAv/zFv2B6VlmDc8eta/D7H9yOP/jQdv2Yf3j2N/jl8Wf055/47GEM1d0DsiiD/hsrpmC1KL/j79q1CzU1ibLZy0khwaehyHMWu3/B1EWhHvWxR133owWQPgBDRMRExFBbQCoPDMEIQKu6uLQFQDsz55xmLYQQhXjmx/+MiM0GrnfgtooK7LdVoe56CM1r1oBNJkRu3YTnfvYDcEyZ1fa/fWQ3TCbCE7/jRF2VMrEg+LIPw+EJAMD2xl24atuM4SkzqnbeD4vJhK0NiYWdmbrclotCut1CRPQYM38v347q+FDaNOZSUcd+/EgOLFprJmt3HzO3pjwPIjEbzlvItd84+S4e+/3iiiUJIRZPhdmEPZvqsNZWkbbtZiSK1y+NYWY2d+G0Yt24cg4ci8FWtQ5fMZn0NYkP1dTiO9Y6nAxFEItOYvS5/xsbG+x4ln6F59R9Np4fxdmh64heP48LBIxMzABN+zD467MAgOqdD2DrxElUGCaarejgw8zfVOv2jDJz1vqval639oXOcrBYZqamcCl4Mv+OQoglY/icGXffVpf8qykDb743jsmZDMk+S2S6YTuu1+4GLiljPVFrNYbv+hgiryuJV2bDw1i/NoI33kjkaDbFGTQ6AjAjxkCkbgvO8Togotzn+vq1+NP2L+I7/98/AgCsVivuuOOOBXsPC63QDAeHAPiJaARK91YQiRZOE5Qp1k4oa4KEEGJJmJ6JYWRiBg21iUWSoYmZBQ08MFtQe89H8P/aqrD3+hAckRv4u7seBW99PywX3sJs+DI2rLWhviZ5YafZRNh5Sy3evDQGMleg5u7kJZWdH78Dn757A/7ll8/j3LlzcLvdsM4zy8BiImbOv5e2szItug3JudWCUAb6n8p81PJHRCfu2HPvgSPf+KvFvhUhRAF+ffo6fvQbJRmno8aK//Do+2AxE2Jxxl/85G1cu6nkRfv4nltxYHfxhdBy2XDrbVhrV5L7041x0MRNxG+9DQAQnZlB6PIFbK2vhMWcecj90ugUTGvXYc3axNytW9fa4FyvrNmJRCIYHh7Gli1bMqUaW1IOHjyI559//nltzaRRUbnd1GnRvQBARNuZ+UyeQ1aMWxrs+OKjB/LvKIRYdJ+djuKV7ucwOhnFTQBn4g14dM9t+MWbwxirHoO1Glhjs+Avn3wYddXpY0Klsy79pTtvy3nEPXnOWFlZia1bt879lpaInLPdMuV10xQTeHKdRwghSm2NrQJ/+FBiJf53Bi5gfCqK7/Sf11978kDjAgcekUvO4MPM40T0ZSLaNpeTE9F2Inqi3IlGhRDi8/dvxcY6GwAgPBnFn33/dVy/OQMAWFdrxb/94LZFvDuRd50PMz8NoJmIugzVS3MiovuI6BsA7mPmglPtCCFEqdgqzPjjR3bqz8+FErnN/v1DO1BtnVdFGTFPBX36zPw0EdUBOEREfwYly3YIwBCUnGh2KPnZGqBkOPAB6MqSiFQIIcqidd9m9L4QxJnrE/prm+xV+NcfkPIoi63g0K8GkqehJA+tgzK12gEl8JyBUjk0KAFHCLFUWMwm/GnzLvzR/3pVf+1rzbtQaZl/SQAxP3Nqd6oB5tW8OwohxCJ79O6N+E7/efz69Aj2b6vHZ+/btNi3JDDPMtpCCLHUmUyEv/u9/Xhn+AbuuHUNzKalvTZmtZDgI4RY8WwVZuy9XYoVLyWFZLUWQgghSkqCjxBCiLKT4COEEKLsJPgIIYQoOwk+Qgghyq5kwYeIHpEEokIIIQpRsqnWzHxczelWDyDEzK+V6txCCCFWloXodiMAjWoi0seI6HEiengBriMA+P1+NDcXV0C2u7sbXq8Xfr8ffr8f3d3dOc+b7Rp+vx+tra1zu/EcOjs7M96TEGLlKGW3WxeA7cx8nJm/y8xHmPl7zPxdAIOFZsQWxXG73fl3MmhubkZLSwtaWlrgdrvhdDoxNDSU8bx2uz3t8XyuXajDhw8vyHmFEEtHKTMctDDzkUwb1Fxw0g1XJs3NzfD5fGmvBwIBhEIhOJ1O/TWn04nW1laEw2H09vbC5XIhGAyira2toGsFg0H4/X6Ew2HY7XY9IPn9fgSDQTidTrjdbv3amqamprTrdXd3w+VyIRAIzPMTEEIsdaXsdvuutG6WhkyBBwAGBgbQ1NSU9rrb7UZXVxfcbjfcbjcGBwcLvpbD4YDb7UZLSws6OzsBKAFJCygejwcAcPToUf1aTqcz7XpaINJeE0KsbKWccPCUOuNtGzOfLdV5l4ttT/1kwc599huPluQ8TqcTPT09aa8Hg0EEAgE0NzcjEAigvb096zm0ANPY2Ii2trak7jiHw6G3dtra2hAOh/VtR44cQVdXF9rb29HX15d2va6uLj1QCSFWvpIFH3VSwZAWeLRJBsz8bKmuIebH7XYnBQhNOBzWJxS4XK6koJEqNUAY97Xb7XA6nQgEAvD7/Uldd36/Hx6PB0eOHEFvb2/a9fbv349AIACn05nUPSeEWJlKOebzEQBOInICGIFSzbQBgASfBeT1ejEwMACv14uWlhYA2cd8AKVLrru7G06nM2lCgcvlSpph5nA4EAgEEAwGEQ6H9cfGoAUorSltH63FFAwGYbfb9f29Xi/6+/v1Y1paWuB0OpOu19HRge7ubgQCAQQCAfh8vrSWlRBi5SBmLs2JiB5h5uPq4zoAh6C0hJZ98CGiEwcOHDhw4sSJxb4VIYRYNg4ePIjnn3/+eWY+mLqtlBMOXFqGA2YeY+anS3huIYQQK0gpJxx8k4iOqa2eAIB+APsh3W5CCCFSzKnlQ0R1mbIWMPMhAE8BCEEZ7/nG/G4v47XdRJQ2F5iIXETUof70EVHWwQIi8hBRGxHZ1T9lbq8QQpTRXFs+DgB+ImIorRw/AB8zP8vMrwJ4tVQ3aKQGiRAAV8rrdgBNzNytPm8BcBzAvhyn86g/XczsX4j7FUIIkdlcx3xcUL7YdwDoBdAIwEtEMSLqV/O63Vuqm9Qws5+ZMy1/bwLQaXjuhzIGla3108/M9eqPJBETQogym2vLh9UWDgA8rf6AiDqgtIoaATxHRP3M/NH532bem/ETkTHDpVN9PfuCFShddVmCmbb9hPpQMjcIIUQJzbXl8/5ML6qtiH5mPsTMDiitoa/P+e6KkBJEDgPI1aJxql14QXX8x5VjXyGEECU21+DTo3avfTbDtnrtgTrdmuZ4jTlRu9pczNyZbR9m7la78MIAegD0ZdnvoDo/fUkmRQ2Hw/B6vfB6vXram0JISQVFsaUoimX83ILB4IJ8VkIsV3MKPsx8BkA7gG+q4zy/IKL/SURHkT7In56v30CdbebJ8VPsTDQPM+f8VjGOBTFzEGo33XJz7NgxhEIhPbNBb29v3mOkpIJCy9gQDAYX7BrGz83pdKKvL+PvOEKsSnNe56N2c+1Qg4MbgB3AMbV+DwCAiJ4BMAjgeznOk/8bs0DqmFOn+tieacxHvV8Pcs+EWxaMudOCwaCe3mY1lFTQcsV1dnYiEAjoKYIAJQg3NTVhYGAg6/vQPi+v14uOjo6810u9f7fbje7ubrjdbgSDQbS0tKS959TjOzs74fP5ct67lvpIO4/2i4UQK828Mxyo3VdPMfOTxsCjCiBPy6dU1OnVXkPAOWTY5jS0dgZgmBmnHVeOe1wowWAQDodDDyqroaSC2+1GKBTSr62dv7u7G01NTXC5XHA6nVlbg3a7Pem4fFLv3+v1wul0wuVy4ejRoxnfs5HL5UpqSWa6d+2cLS0t6O/vl8AjVrRSJhZNw8xPlfJ8aqulWX3sgbK2yK8mM+1TX9d2D0KZBg4oLR0fgF5mDhNRSG0lhQE0MvP8O+P/c928T5H93GM5N3u93oylElKttJIKqUlOASXwal/a2vtNbf0EAgE0NTXBbrcjHA7rrTYg0fLQgle2+/f5fPrnpHWnpb7nYu/d5XKhp6cHLpdrwcejhFhspczttuDUVlYnM5P6p199Pai+ZvxpNBzXauzeY+aAOumgN9fEhOXA2G3k9+deK2ssqWCUWlIh0xejxuPxwOPx6F/o2UoqpE4Y0LqaBgcHkwbitetpJRUAzKukgtaNByhBdf/+/Wn7aFm3ASXD9rFjx/R7NHYd5rr/xsZG/T61rN+lmCShjXcV2u0pxHK1oC0fsbD8fj86OzvR1dUFIFFrZyWVVACARx55JK0rUNtPG3MKBAIIh8PweDz6eQOBQNp4jt/vR09PT9JkAO25z+fDkSNHYLfb0yZYZLp/rRWoBeDU96wF4tTP0PjYeO8A0NXVBYfDgVAohCNHjuhjQUKsNCUrqbCSSUmFxZUp6C0E40QHY3Aul87OzqSux9TnQiw3uUoqzKnlQ0RHmTnrfFgi+hsA26GMyXxrLtcQQlPoOMp8LXYr4/Dhw3qLKRgMlmXKuRCLZa7dbr0AQER7AQSZeVzbQERdUIrIPUlEjxDRY8ycdaq1EPksdlAoF5fLpb/X1fKexeo15wkHRHQayhTlMykpdFqY+ZsAoFY2bZjfLQohhFhp5pXVmpl3MHMDgDEi2qZuS02ns3BLyIUQQixLcw0+QWY2Lj45hkSKmtGUfWVGgxBCiCRzDT4OInqYiNaqLZ6nAAyoJbTrU/ZdlnnThBBCLJy5JhZ9Gkr6mrNQMgsEoWQe+AaAdiL6OhFtI6InIN1uQgghUsx5woGay83BzPuZ+Wlm/i4zf0WdZHAcSmsoxMzPluxuRRqtLIKUVChMvtIG2d5Tts8gk3A4rKcs8noLSxuYaz+tdIYQK8m80usQ0RNaKQW1lQMAYOZX1eAkU6wXkJYJQMsaXUh5gNVeUiFfaYNs7ynbZ5DJwMCAvki1kL8Tv9+fc2q1dt2FLP8gRLnNZ6r1MwCaoHSrDQBoUgvMrS3VzYncXC4XPB4PwuEwnE6nngUg22/o+UoqdHd3w+/3F1QXSKOliNFaUhrtPNprgUBAb2lpJRhSr6c9z5ejDlBaCvv27UM4HIbf79dT+6ReV2uxeL1ePX2P8fNJ3T/Xe9Lk+qwCgYD+mba3txdUrsHn8yX9nXi9XrS2tsLr9epZF1paWgpKHivEssHMRf8AeAJAXYbX7QC+PpdzLuUfACcOHDjAS5XP5+OOjo68+/X09HBbW1vGbR0dHTw4OMjMrO/T0tKibzc+NnK73fpjl8vFzMxDQ0Pc09OTtL2jo4N9Pp++PfV6PT09+vbBwUH2eDx53492bu08ma7LzOx0OpOO095Ltv0zvSfjcZk+K432HpiZPR4Pj46Opt23z+fTr5vpHNoxxnMZry/EcnHgwAEGcIIzfK/ONcPBKCdPtdYCWZiIzszxnMva3f9w94Kd+/Xfez3ndrfbjb6+Pni93pw1YFZaSYXW1taklkmm6wLZswVk2z/be9IU+llp3W7BYBA+nw+tra0IBoM4dOhQzi40u92e8e/S4XBkPUaI5Wau3W651u7Iup4y6ezs1Lt97HZ73lIEK62kwqFDh9DT06N/KRdb1iDb/pnek1Guz8r42Q4ODuopc/bt21dwoT6/34+Wlha91pAQK9FcWz4NRLSNmc8aX1RzvTVmPkSUWnt7e1IZay0grIaSCoASGBwOh96yyXRdu92un1Mr0a29l0z7a/eW+p6Mx2n3qtEmKWgVZbX9Ms1AzBRYjS0t43XtdntSa7DcWbaFWEhzLqlARMegZK7WftVzQsl8sOJS8UpJhcVVrpIK8+X3+zPOlgsGg/B4PPrEDq37U+sC9fv9SRNGsp073z5CLDUlL6kAAMx8iIjugzLjzQ7gG8z86txvU4jMlkvXU64xGa3bTWMce3O73XophUy09y+BR6wk86pkqgabpIBDRF9nqeEjSmi5lBfIdp/BYDDvWE+uiSLZSnsLsZzlDT5E9IsizkcA9gGQ4COEyu12S/AQIkUhLR8C0AmgkL4PgpLfTQghhMiqkODTWcxYDhEVnmRMCCHEqpR3nU+xkwiYeVUuMhVCCFG4eSUWFUIIIeZCgo8QQoiym9dUa7F0dHZ2FpQbLRwOo7e3NynDQSAQQEdHBwKBAL785S/D7XZj//79CAaDcLlceskG47ZQKIS+vr6smRSEECIXCT4rgN/vL7jWS2trK3p6evQFi+FwWK9vo+UqO3z4sL5mhYjAzBm3CSHEXEm32zKXKfVMtno+WgZo4/52uz1rZmav15u1NRUIBGTtihBizqTlUyJv33Hngp37znfezrotGAymBYFsXWFaIs1UqS2ZgYEBhEIh+Hy+tOCjFXuz2+164k8hhCjWsmv5EJGbiNJylaivt6g/HiLK2jdERE4i6lCP6SCiZZkuOFsiy2yampqSuue0zMtaRVDjfm63G83Nzejq6ko6h9vtRltbm9SWEULMy7Jq+RCRG0AIQKbA0gdgu1rQDgCehpLqJ5MeZm5WzxkE4AGQvSrYEuVwOPRyClpBuFzjMdq4jdZV53Q60dzcjHA4nLFFpJUjANJr3OTKRSaEEPksq+DDzH5AGQTPYB8za9+QDihBKg0ROdXt2jmDRHQI8ww+ubrGFooWaHp7e5OCQ656Pn19fXo9H0Bp/TQ2KiWYtECjpe93u93weDzwer1Yt24dgsEgjh49mrHAmhBCFGPO9XwWExExM2eMQOp2HwCPFqxStrUAaNdaPupro1CCVzBl3xPqw70HDhyok3o+QghRuAWp57MUqa2adgB9mQKPyoH0JKkhKDWJhBBClMGym3CQi9py6QKwT23hzPd8B9WI/dp8zyWEECJh0Vs+RNQGoDHHLr4crZg06oSDPgA+Iqo3jANpMrVyMrWGhBBCLJBFDz7M3Dvfc6iz4DzMrM1uG1D/zBRUAjBMODDcR2EpAoQQQszbSul2CwHoMTxvAhDUAoq6rscOpAcZdZzoWLluVAghxBJo+RRDbeFo63M8ULvkmDlARA61Cw9Q1vcYc8x4APgAaK2sViLqABAEsJ+Zl90aHyGEWM6WVfBRx378UMp6Z9qW7bjWlOdBAN3qU28p71EIIUR+K6XbTQghxDKyrFo+Il19fX1SNoJ8ylnPx3hMe3t7WbIidHd3w+Vy6SmHOjo6CjrO7/fD4/FIfSIhykSCzzLX19dXVHLRctbzMR5TjsDj9Xr1YAkoBfYK5Xa70dPTk39HIURJSPApkb9+8tkFO/cf/s3DWbdpv+Ebv9yz5XZbSvV8tAJ4WqtNa3m0t7frLRYtB51xv97eXrhcLgSDwbSSDk6nE52dnXrS1CNHjujburu74Xa7EQwG0dLSknZ9I611mO06Qoj5kzGfZS4UCsHhcCQFkFLU8/H7/fD5fGlfvH6/H93d3RgYGJhzayYYDOpf6lpwMwaGjo4OeL1eOJ1OtLS0oL+/Hy0tLejq6oLb7Ybb7cbgYFpVDbhcLrS2tqK1tRWNjY16+QjtXC6XC0ePHs14faN81xFCzJ8En2Wura0NdrsddrsdXm/uiXtLoZ6PljG7ra0t6ZpAchB0uVzo7+9HMBjUK7MGAgGEQiEEAoGMrbVwOIy2tjYMDg7C5/Pp3W4+n08PlH19fVmvr8l3HSHE/Em3W4nk6hpbKL29vXA4HGhpaUFDQ0Pe/ZdCPZ9AIKDXIcrXnXX48GEA0PfTgpA2oSDVsWPHcOjQIb3kg/aeGhsbEQqF9PcRDAZzXj/fdYQQ8yfBZxk7dOiQ3kUGQA8I5a7nEw6H8cgjj6R1UQUCAQQCAb2ry+fzwW63691/WhD0er16oDMWxOvq6oLD4UAoFMKRI0f0cSBNpjEnrcR3MBjUWy0dHR16K0gLJqnXdzqdCAQCSeNNua4jhJifZVnPp9yI6MSBAwcOSD2f7FInPcxXZ2dn0nhM6nMhxNK3aur5iMVT6u6pw4cP6y2SYDCod8EJIVYGCT6iJApd+1PM+bRzlvrcQojFJ7PdhBBClJ0EHyGEEGUnwUcIIUTZSfARQghRdhJ8hBBClJ0En2UsEAhg37596OzshNfrRXd3t77gtNT8fr++8t+os7MzaUHmXIXDYXi9Xni93qRs1F6vV88nV8x++Xi9XtTX1+tTxLU0PvP5/Orr6/W/j1z3Vey9CrESSfBZxowlC7SEnJkCRCm43e6MKXgKXX+TL+/csWPHEAqF9CwNvb29+jHatf1+f8H75WO32+HxeNDbq1RW1zJizyebQV9fHwYHB/XFsJnuay73KsRKJOt8SuSTn/zkgp37Rz/6UUH7GUsgpJYFcDgc6OrqwvHjxzEwMIC+vj60trbC4/Ggs7NTL5Ggranp7e1FU1MTBgYGMuZA04q2aSl5gERCTk0xX+TGa2ipcXp6evTgpqW/MRaHy7VfIddua2tDY2MjOjo6EA6H552hIbW8RX9/f9p9jYyMzOlehVhpJPisAAMDAwiFQvD5fHrw6erq0gu/aV/QR48e1Y/RCqdpv+07nU709PTA5XLptW9cLhdCoRB6e3uTgoMW1Nxut54kFACOHj2K5uZmvTzCXGiB0ul0pmVNGBkZKXq/fLRaQtrjTPeTrXWSGpSN5S16enoy3td87lWIlUSCzwrQ1NSkZ2Du6uqCx+NBIBBAc3NzUlmA9vZ2PXBoMv227/P59G4tLSgZv2iNQc7oyJEj6OrqQnt7O/r6+pK+uAcHB5NaRdkySnu9Xj0w2u32pGPmsl82WkmIzs5OtLe3Z616qpVfKIS2n1beItN9zeVehViJJPiUSKFdYwvJWALhwQcfBJBcFsDtdqOzszPvl6nWVaflVdu/f3/S9v379yMQCMDpdCZ9kWrVSI8cOYLe3l50dHTo1/J6vXnLMHi9Xr1bze/3Y//+/WkTAorZLxu/36+3crT3kC03XaEtn0zlLTLdVzgcLupehVipZMLBMmYsgRAOh/UvVK/Xi6997WsIBALw+/0YGBjQj/F4PPrEAa2EgfYFGwgEEA6H9ZaT9lpHR0fSvh0dHQgGg/prPp8P4XAY/f39+kyuYuv9+P1+vbDdvn37ACglIoxf/loXWSH7hcNhfXvqZ6a9P82RI0ey5o/TWj6Zfoy0OkLG8haZ7ivTa0KsRlJSoQAroaSCVpLAWC+nnApp+ZRaqcs8CCGKk6ukgrR8VgltDctiZYgud+ABSl/mQQhROjLms0qsxu4dKcUgxNIlLR8hhBBlJ8FHCCFE2UnwEUIIUXYSfIQQQpTdsptwQERuAB5m3pfhdS3z5X4AR5k5kHq8uq8HwBCAYwAOAQgys2R4FEKIMllWLR81wIQAZJrG1AfAz8xeAP0Ans5zOg+AMwDsEniEEKK8llXLRwsSRJRp8z5m1hZ2OKAEqWz6mTlzMi8hhBALblkFn1yY2ZhGuRVKyyYnInJl65pTt59QHz7w2muv4eDBg/O6RyGEWE1ee+01ANiRaduyTK9DRMzMac0fInICaAcwxMy9OY7vABAAMADgCLKMDxmCzz4AYwBOz//uF9Re9c/XFvUuVhb5TEtPPtPSW6qf6Q4A15j5vtQNKyr4qNvsUFo9PnX8J9+5nOq+jSW+zbLTgmWmPEpibuQzLT35TEtvOX6mi97tRkRtAHJ98fuKmRDAzGEi6gPgI6J6wziQ8Zp27XVmDqoBSAghRJksevDJ1T1WqAzTr7UaAg4A4Uz7QulKE0IIsQhWRLcbEbkANGmBTA0wPVpXmtqyCamtIru6rzZzrgXAYWZuLfsbEUKIVWrRWz7FUINKs/pYG9fxM3OAiBxqFx6gtGqMJSI9AHwAetUAFFInHYQBNErgEUKI8lqWLR8hhBDL27LKcCDyIyI3EQ1meN1JRG1E1EJEHWr3oyhAts9U3dZi+FxXX9GkOcr1mRr26SnX/Sx3Of7du9R/7x1E1LeU/t1Ly2cFMaQfGkydik5EHczcbXjukSwP+eX5TFsAOJm5Wx1X7EvNOSjS5fpMDfu4cm0XCdk+TzXQHDKMhbcAOLJU/h+V4LMCZVoHRUQ+Zm42PO9h5vby393ylOUzHVoJ68MWS571ei0Anmbm+jLf1rKVYSJW6sQrO4BRABmXoJSbdLutIkTkIyK7+j9l32Lfz3Km/mae9ljMHxG1FLJAXOSmzug1TqZyqq8veuABJPisGmqrxwElk7dLMnnPmxNAWP0NPaj2qcuYzzyp3ZfBvDuKgqSkDTsMoDvbvuW2rKZai7lTvxg7oXxp9hARjGNAomgOKEHcCwBE1AslsEs30fy4pNVTemqXm8vY9b7YJPgscaVIP6T+NunSgg0R+QEMElHvUmmCl1OJUjoFoSSnBaCndbITkTMlw/qqUKL/T90ApEWO0qcdg5IBZskEHkCCz5JXivRDAJL+Uav57HqRIf3QalCizzSIROVcTRi560itWCX6TAHgkKFel139EvavtoBews9Ty+LfqT62L5VfOCX4rA5+AC0w/KYOpNVAEkVQA3hQ+8esjVUslX/Yy1Hqb/LqjMySfQmvRuqYpNfw/+UhAEviM5Xgs4LkSD+kfVF2QPmN3QFAFvAVINtnqm5uBXCEiIagdJE8sjh3ubzk+Uy18Yk2w/Ye+UUpu2yfp7b2TH1d2z2IJRJ8ZJ2PEEKIspOp1kIIIcpOgo8QQoiyk+AjxCq2lBJNitVFgo8QWagZgX1E5Fvse1kIRNSWb3aemrF7kIhG1Qkr2fZhY9Zk9TUJbCIrmXAgRA5q3rYlt0BvvopZP1NIhulMWdJTM6kLYSQtHyFyW3HrdtQWSWOh05fV/GBBdc1IpvM5oVQKThWQfHciGwk+Qqw+hwAcLfIYL4Bsqykb1gAAA55JREFUJTjcmVK9qK9J2Q6RkQQfIVaf5pRsx4XoAeCewzhOSG0ZCZFEMhwIUSR1vETL4eZITQFj2O5EIqVROzMba6sY97VD6d7zQ8nDBygr1rUxFBeUrBTNAL48nxQ+uYKHOqEgoN5P0vtSs2QEoGQeMFbEdSF3MtBB9f4lQ4FIIi0fIYpARH1QBuq9aur/Y8bZcETUAyXHmxdKGpMeNcVRpsCj1a7xA/BAyTzeq37p90NpbTjVa/VCCWht83wLTmRIfmp4X3713hszjPH0IL0brSnP2FEQwP753LBYmST4CFEg9bf8pJIJaiskqLZgAGU8ZcCwLVelU6c6LtKERMDSNKrnMLYqnJj/BAgHgCHjC4aSG8auuKNIDzTHADhT3k++LN6Zsn8LId1uQhShCZm7j4YA7FMfa4lbtSDhyHKMMbA0I30CQBOArgyvdWJ+tC4+IzeUqqzulP2SxmrU7N1eKBUxA2rLKF9NmRCUz0CIJNLyEaJwuX6D175gewC0q4Xl2gB0FTBGk6mIWtJYihYY5jBRIFUY6e/DDqXl5Tf8eJk5UzGzHiS6/hwFvDcHVmmNI5GbBB8hCudHSmtA1YjEOpcQlBZLE5QxlJyLLNUuL7sxqKiBJrU2UCvUVPjGFgoRdRCRm4g8Ba6pCSG9QmYgy/tKo7XW1MBayCSCUnQVihVIgo8QBTIsttTHPNTZY02GmWHN6r6FVt/M1OppzvBaEwCfej0thU0fgF41IGiTF/LRugWN78sPZUp00thUjmDWCyXrQyFlnJ1QJk8IkUTGfITIQm2VeAA0aalimLlVbW1oLQUnkovI9QEYJCJj95Y3NfWMQSPSC/s5M7zWA6UrzsnMvWqgMJZEdhYS7NRxm7TuQ2ZuVt9XE9RuspQJEKn3Uqh9mP84lViBJLebECWiBgR3alebun6mkZlLttpfPWcDM3eqwaSv0Pxz6rH+EowfFXKtvkzTzIWQbjchSqcdShqaJGowKvUq/zAS3VmHoHTJFZpHrRfKjLUFpd6PlGsXGUnwEaJ0fEhkKNCpX8KlLstwDMB+bXICgIZCD1S76vrLkPbGVeC4kFiFpNtNiBJSg4Fx8N8JIJRj/GTRqPV8evPvOadzt0Dp2pOZbiIjCT5CrGJEZF+IALFQ5xUrhwQfIYQQZSdjPkIIIcpOgo8QQoiyk+AjhBCi7CT4CCGEKDsJPkIIIcpOgo8QQoiyk+AjhBCi7CT4CCGEKLv/HxNiOZj77vNMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# set up some labels and colours to use to match those in the paper\n",
    "labels = [\n",
    "    \"Cell-based\", \n",
    "    r\"Cell-based, no $\\Lambda_c$ scaling\", \n",
    "    r\"Cell-based, variable $\\beta_{\\rm pl}(z)$\", \n",
    "    \"GRF\", \n",
    "    \"GRF, Large Scale\"\n",
    "]\n",
    "colors = [\"C0\", \"C1\", \"C3\", \"C2\", \"C4\"]\n",
    "\n",
    "# cycle over each model and plot the 99.7% excluded region\n",
    "for i, label in enumerate(labels):\n",
    "    plt.plot(data[\"logm\"], data[\"logg_{}\".format(i+1)], label=\"{}: {}\".format(i+1, label))\n",
    "    plt.fill_between(data[\"logm\"], data[\"logg_{}\".format(i+1)], y2=-10.8, alpha=0.2, color=colors[i])\n",
    "\n",
    "# plot the reynolds limits\n",
    "plt.plot(data[\"logm\"], data[\"logg_r20\".format(i+1)], c=\"k\", lw=3, ls=\"-\", label=r\"Reynolds+ 2020, $N=500$\", alpha=0.7)\n",
    "\n",
    "# aesthetics\n",
    "plt.legend(frameon=False, loc=3)\n",
    "plt.xlim(-19,-11.1)\n",
    "plt.ylim(-13.5,-10.8)\n",
    "plt.ylabel(r\"$\\log g_{a\\gamma}~({\\rm GeV}^{-1})$\", fontsize=18)\n",
    "_ = plt.xlabel(r\"$\\log m_a~({\\rm eV})$\", fontsize=18)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
