{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import uncertainties as un\n",
    "pylab.rcParams['figure.figsize'] = (10, 6)\n",
    "from scipy.interpolate import InterpolatedUnivariateSpline\n",
    "\n",
    "\n",
    "from matplotlib.pyplot import cm \n",
    "from scipy import interpolate\n",
    "from scipy import stats\n",
    "\n",
    "def print_now(msg):\n",
    "    print(msg)\n",
    "    sys.stdout.flush()\n",
    "\n",
    "\n",
    "pylab.rcParams['figure.figsize'] = (12, 8)\n",
    "pylab.rc('text', usetex = True)\n",
    "pylab.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})\n",
    "pylab.rcParams['axes.labelsize'] = 25\n",
    "pylab.rcParams['axes.titlesize'] = 15\n",
    "pylab.rcParams['xtick.labelsize'] = 15\n",
    "pylab.rcParams['ytick.labelsize'] = 15\n",
    "pylab.rcParams['ytick.labelsize'] = 15\n",
    "#pylab.rc('grid', c='0.5', ls='-', lw=0.5)\n",
    "\n",
    "import os\n",
    "import multiprocessing\n",
    "\n",
    "import scipy\n",
    "#from scipy.interpolate import spline\n",
    "from scipy import interpolate\n",
    "import time\n",
    "\n",
    "import matplotlib.colors as colors\n",
    "import matplotlib.cm as cmx\n",
    "from matplotlib import cm\n",
    "from matplotlib.collections import LineCollection\n",
    "\n",
    "def multiline(xs, ys, c, ax=None, **kwargs):\n",
    "    \"\"\"Plot lines with different colorings\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    xs : iterable container of x coordinates\n",
    "    ys : iterable container of y coordinates\n",
    "    c : iterable container of numbers mapped to colormap\n",
    "    ax (optional): Axes to plot on.\n",
    "    kwargs (optional): passed to LineCollection\n",
    "\n",
    "    Notes:\n",
    "        len(xs) == len(ys) == len(c) is the number of line segments\n",
    "        len(xs[i]) == len(ys[i]) is the number of points for each line (indexed by i)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    lc : LineCollection instance.\n",
    "    \"\"\"\n",
    "\n",
    "    # find axes\n",
    "    ax = plt.gca() if ax is None else ax\n",
    "\n",
    "    # create LineCollection\n",
    "    segments = [np.column_stack([x, y]) for x, y in zip(xs, ys)]\n",
    "    lc = LineCollection(segments, **kwargs)\n",
    "\n",
    "    # set coloring of line segments\n",
    "    #    Note: I get an error if I pass c as a list here... not sure why.\n",
    "    lc.set_array(np.asarray(c))\n",
    "\n",
    "    # add lines to axes and rescale \n",
    "    #    Note: adding a collection doesn't autoscalee xlim/ylim\n",
    "    ax.add_collection(lc)\n",
    "    #ax.autoscale()\n",
    "    return lc\n",
    "from astropy.io import fits\n",
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set path to class and change directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_to_bard_class = '/Users/boris/Work/Boris-BIG-MAC/bard_cmb_lab/codes/class_public-2.7.2/'\n",
    "path_to_notebook = '/Users/boris/Work/Boris-BIG-MAC/bard_cmb_lab/notebooks/'\n",
    "\n",
    "os.chdir(path_to_bard_class)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the Planck data using \"FIT\" file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filename: /Users/boris/Work/Boris-BIG-MAC/bard_cmb_lab/notebooks/COM_PowerSpect_CMB_R2.01.fits\n",
      "No.    Name      Ver    Type      Cards   Dimensions   Format\n",
      "  0  PRIMARY       1 PrimaryHDU       4   ()      \n",
      "  1  TTLOLUNB      1 BinTableHDU     46   249R x 4C   [I, E, E, E]   \n",
      "  2  TTHILBIN      1 BinTableHDU     46   83R x 5C   [E, I, I, E, E]   \n",
      "  3  TTHILUNB      1 BinTableHDU     42   2479R x 3C   [I, E, E]   \n",
      "  4  TEHILBIN      1 BinTableHDU     46   66R x 5C   [E, I, I, E, E]   \n",
      "  5  TEHILUNB      1 BinTableHDU     42   1967R x 3C   [I, E, E]   \n",
      "  6  EEHILBIN      1 BinTableHDU     46   66R x 5C   [E, I, I, E, E]   \n",
      "  7  EEHILUNB      1 BinTableHDU     42   1967R x 3C   [I, E, E]   \n"
     ]
    }
   ],
   "source": [
    "hdul = fits.open(path_to_notebook + 'COM_PowerSpect_CMB_R2.01.fits')\n",
    "hdul.info()\n",
    "#To see description of data: hdul[i].data, with i = 1,2,..\n",
    "low_ell = []\n",
    "Dl_low_ell = []\n",
    "err_m_Dl_low_ell = []\n",
    "err_p_Dl_low_ell = []\n",
    "for i in range(len(hdul[1].data)):\n",
    "    low_ell.append(hdul[1].data[i][0])\n",
    "    Dl_low_ell.append(hdul[1].data[i][1])\n",
    "    err_p_Dl_low_ell.append(hdul[1].data[i][2])\n",
    "    err_m_Dl_low_ell.append(hdul[1].data[i][3])\n",
    "    if low_ell[i] == 29:\n",
    "        break\n",
    "        \n",
    "high_ell = []\n",
    "Dl_high_ell = []\n",
    "err_Dl_high_ell = []\n",
    "for i in range(len(hdul[2].data)):\n",
    "    high_ell.append(hdul[2].data[i][0])\n",
    "    Dl_high_ell.append(hdul[2].data[i][3])\n",
    "    err_Dl_high_ell.append(hdul[2].data[i][4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Plc_BF = loadtxt(path_to_notebook + 'COM_PowerSpect_CMB-base-plikHM-TT-lowTEB-minimum-theory_R2.02.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/ipykernel_launcher.py:22: DeprecationWarning: In future, it will be an error for 'np.bool_' scalars to be interpreted as an index\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAFhCAYAAACPlvgAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU9b3/8dc3hLCEJQlBFEUkqBUtSgG76E9FhciI2ygBta5og7sWq2jbW7XeK+LWWlsq2Gst11oVdGwVBwFrXWprRapiq1WMVLEQCCEsYQlJvr8/zgwMwySZ5cxyJu/n45EHOed7lu85GfLJdzfWWkRERMRbCrKdAREREUmcAriIiIgHKYCLiIh4kAK4iIiIBymAi4iIeJACuIiIiAcpgIuIiHhQYbYzEC9jzM1AA1APYK2dH5FWHd4PVFhr74k6N6V0ERGRXOOJAG6MWQxUWWsbQtsbjDFLrLUN4eAbDujGmApjzGxr7dTQdkrpIiIiucjk+kxs4ZK3tXZOxL4Ka21N6Pt3rLWjos751Fo71I10ERGRXOSFNvBbgSWROyKCdwkwMsY5DcaYsammp5hvERGRtMnpKvRQgC0JfT8Rpw18JDAnVJ1eEdoXrT6UVp9iuoiISE7K6QAOjMYJsCURbdRLgXnAOKCM3Z3PIjXgBP5U02MqLi620U0PxcXF7L///h08jmPTpk306dPHlbTofZHb7aVt2LCB0tLSuPLbnvbyG23dunUcdthhrF+/PuX7eoVb79mrOvvzg96Bnt95/vXr1+/1u2/btm1brLW9k764tTZnv4CxgMUJ4JH738EpiY8FPo1x3jzg5lTT28rXqFGjbLQLL7xwr31tae/YRNOi90Vut5d21FFHdZjPeCTy3Pfdd58r9/QSt96zV3X257dW70DP3/bzA0ttCjEy19vAawBsqPd5hHqc4AtOKTpaZOk51XQREZGck9MB3IY6q7WhAVhK7GBbBixzIV1cNGjQoGxnQUQkb+R0AA9ZZoyJ7lBWgVP10ADUhDq7RSqx1i5JNd29RxCASZMmZTsLIiJ5wwsBfHroCwBjzEigxlobLiHPBKqj0iODb6rpcTnllFNcOTbRtOh9kdvtpV100UUd5jMeiTx3Z+TWe/aqzv78oHeg50/f8+f8RC6wawhZuBTez1o7PSq9Gqe9vIS2p0pNOj3a6NGj7dKlS/fYV1dXR3l5eVzP096xiaZF74vcTiQtWYlc59FHH2XKlCkp37Mjt9/ufOUCt96zV3X25we9Az1/288fmkhsdLLXzvVhZMCe8563kT4nnenijg0bNmTkPnfckTsBXEQkXbxQhS4iIiJRFMAlrwQXtXLw8CaCi1qznRURkbRSAJe8EljQzMx56wksaM52VkRE0koBXDKmoiL908v7JxQyvaof/gme6N4hIpI0/ZaTjPH7/em9wWOP4Wto4NPl1+Gr1N+mIpLfFMAlPzz1FNx1F/Tty/VY4LvZzpGISFqpmCIZ8/DDD6fnwtbCnXfCrFnw2GPcWXw37NyZnnuJiOQIBXDJmMbGxvRcePly2LwZTjoJjjiCXsMr4OWX03MvEZEcoQAu3vf738M550BB6ONcVQWBQHbzJCKSZgrg4n2vvw5jxuzePukkeO21rGVHRCQTFMDF25qb4a9/hWOP3b1v+HBYswbWrs1evkRE0kwBXDLmsMMOc/+i770HBx4I/frt3telCxxzDLzxhvv3ExHJEQrgkjETJkxw/6Jvv82yrt/Ye//IkfDuu+7fT0QkRyiAi7e9/z6PvXvU3vuPPNIpnYuI5CkFcMmYn/3sZ+5f9P33eZ8j995/1FHw/vvu309EJEcogEvG7HR7chVrYflyljN877ShQ51ObBs3untPEZEcoQAu3vXvfxPsew5lw3vvvXxoly5wxBHwwQfZyZuISJopgIt3ffghgYor214+9CtfgY8/zny+REQyQAFcvOuTT/D3f6ft5UMPOQQ++YTbb894zkRE0k4BXDJm+PAYbdWp+OQTfMdsY8XyotjLh4YC+B13uHtbEZFcoAAuGVNZWenuBVescIJ0Ww45hOCH+3Pw8Ka928hFRDxOAVy865NP4JBDuO22NtIPOYRA19PabiMXEfEwBXDJmPvvv9+9i+3cCV98AUOGtN3G3bcv/tq5/KiqV+w2chERD1MAF29auRL23x+Kito9zDfkU/otXxa7jVxExMP0W0286dNPnclaOnLQQQzm3+nPj4hIhimAizd9/jkMHtzxcYMHc8mJn6c/PyIiGaYALt70+efOMqIdGTyYk4aqBC4i+UcBXDJm1KhR7l3siy/iC+AHHgj/VgAXkfyjAC4ZM2bMGPculkAJnM9VhS4i+UcBXDKmqanJvYvFG8APPNA51lr37i0ikgMUwCVjHnroIXcu1NoKq1bBAQd0fGyvXtCjB9TVuXNvEZEcoQAu3lNbC6Wl0L17fMcPHqx2cBHJOwrg4j3xVp+HhavRRUTySM7PL2mMGQtMBWYADcBEoMFaOyfimGqgPrRZYa29J+oaKaVLjvn8cxg0KP7jFcBFJA/lfAAHSoAK4B2cAD4nMsCGg6+1dn5ou8IYM9taO9WNdMlBiZbA998f/vOf9OVHRCQLvBDAsda2N4B4amS6tbYmVGp3K11c8q1vfcudC33xRWIl8IEDYflyd+4tIpIjPN0GbowpAUbGSGowxoxNNd3NvAocc8wx7lxo9WonKMdr4ECVwEUk73iiBB4KpiU4VegjI6rQK0L7otWH0upTTBcXbdmyhV69eqV+odWrYb/94j9eAVxE8pAXAvgycKq2AYwx9caYxdbacUAZuzufRWrACfippsdUW1vLiBEj9thXVVXF1KnxNZtv3LjRtbTofZHbiaQlK5HrzJo1i2uuuYbW1taU7lm6ahWbunenJc6x3aZbN0q//JL6LIwFd+s9e1Vnf37QO9DzO88/d+5c5s6dG51cnsq1cz6AhwN3xPYyY8xoY0zWSsgDBgxg6dKle+yrq6ujvDz+n0V7xyaaFr0vcjuRtGTFe53i4mLKyspSvl/TqrWUHn449OkT3wn9+kFrK+XdukHv3infP1FuvWev6uzPD3oHev5ypk2bxrRp0/bYb4xJqVTh1TbwGiDcRh0rIkSWnlNNl1yyeTM7d9rEArExqkYXkbyT0wE8NKRrQzuHLCV2sC3DqXpPNV1yzerVrGY/JygnQgFcRPJMTgfwkBkx9lUAS6y1DUBNqDd5pBJrbcrp7mRf3BRcsIMfDH+c4KIE29E1FlxE8kxOB/BQ+/cevcSNMROBpyPaxmcC1RHpI4HI4JtqurjkhBNOSPkagb+VMHneYAILmhM7USVwEckzXujENic0WxqEqrsjZ0kLp0cMNatwM13cM3r06JSv4S9fxkNVB3PtfQMSO3HgQGcCGBGRPJHzARycIJvNdHHHhg0bKC0tTekavh5/5tXlH+KrPCKxEwcOhLfeSuneIiK5JKer0CW/PProo6lfJNyJLVGqQheRPKMALt6yejVjL0gigO+3nwK4iOQVT1Shi+yyejUXPrBv4ucNGAC1te7nR0QkS1QCF29JdB70sN69oaUFGhvdz5OISBYogIt3bN8OW7Y4U6MmyhiVwkUkryiAS8aMGzcutQusWeME4YIkP7YK4CKSRxTAJWOOPPLI1C6QbPV52IABPPmgAriI5AcFcMmY2lRLv2vWpBzAX3lKAVxE8oMCuGTM448/ntoFamudavAkBXceT83wExOfR11EJAcpgIt3rFsH++yT9OmBLcdy5bzeic+jLiKSgxTAxTvWrk0pgPuP+JK5VavxT9D0ByLifQrg4h0pBnDfSa3cuPx6fJX62IuI9+k3mXhHigGcAQM4vEyd2EQkPyiAS8aceuqpqV3AhQDer1kBXETygwK4ZMywYcNSu0CqAbykxJnNbfv21PIhIpIDFMAlY7744ovkT25uhg0bkptGNcwY5w8AzcYmInlAAVwy5umnn07+5PXrobQUunRJLROaTlVE8oQCuHhDqtXnYQrgIpInFMDFGxTARUT2oAAu3pDiLGy7KICLSJ5QABdvUAlcRGQPCuCSMWeddVbyJyuAi4jsQQFcMmbo0KHJn+xmAF+7NvXriIhkmQK4ZMynn36a/MluBfB99lEAF5G8oAAuGfPcc88lf/LatdC/f+qZUAAXkTyhAC7e4FYJvF8/aGhwZnYTEfEwBXDxBrcCeJcuzoxu69enfi0RkSxSAJfct22bswBJ377uXE/V6CKSBxTAJfeFJ3Exxp3rKYCLSB5QAJeMmTRpUnInulV9HqYVyUQkDyiAS8YMGjQouRPdmkY1TGPBRSQPKIBLxnz44YfJnZiOErgCuIh4nAK4ZMyLL76Y3IlujQEPUwAXkTxQGL3DGPM9oF86b2qtvTXZc40x86y1VVH7qoH60GaFtfYeN9Mly1QCFxHZy14BHJgKTAdc6vK7lzlAUgHcGDMSmBi1rxqot9bOD21XGGNmW2unupEuOWDtWhg2zL3rKYCLSB6IFcCNtfbZdN3QGDMzhdMrYuybaq0dFd6w1tYYY8a6mC7ZphK4iMheYrWBL0nzPZO6vjFmYriUHLGvBBgZ4/AGY8zYVNOTyaekQV2d+23gGkYmIh63VwncWntFOm+YzPWNMRVATYykCqAhxv76UFp9iuniogsuuCC5E9etczeA9+oF1kJjIxQXu3ddEZEMitWJrQ9OR653Y6QdhNNevCn9WdvDyOjSd0gZuzufRWoASlxIj6m2tpYRI0bssa+qqoqpU+NrNt+4caNradH7IrcTSUtWItfp0qUL9fX1tLa2JnSPfuvWUW8Mtq4u0ey1qbS8nI0ffUTr4MGuXTMWt96zV3X25we9Az2/8/xz585l7ty50cnlqVx7jwBujFkEjA19b4GZ1trvRxyyAZhsjLnbWpvWnuoReRpL+qv1EzJgwACWLl26x766ujrKy+P/WbR3bKJp0fsitxNJS1a813n//fc58sgjE7v41q3Q3Ey/gw5ybypVgH33pay5GVx6B+1x6z17VWd/ftA70POXM23aNKZNm7bHfmNMSqWSXW3gxpi7gWVAqbW2ADgEKDPGvBQ+xlq70Vr7CFCayk3jFWqjxlobq5o7rCzGvhIX08UlixcvTvykcPW5m8Eb1JFNRDwvsgReEtk+ba2tAa4wxow0xjwFfCei6txmKH/VsGv42C7GmJtxqrmfJnawLcP5Y2RpiumSbXV1rG7pz35uX1cBXEQ8LjKAfxrrAGvtMkLV5saYuzLZ/h1rQhVjzMzI/caYGmNMSVQpvcRau8SNdMmydetYvkYBXEQkWuQwsvaqqbHW3gJUG2OGpDdLCZtJqKQOu0rrS1xMlywKvt6NXwy/h+CixDq+dUhDyUTE4yID+NPGmO8YY24yxqyPdbC19j6c4VXpmqWtTaFx27ND388Oj9O21s5h97juicDYyFnUUk2X7Ap8eCAXzxtAYEGzuxdWCVxEPG5XFbq1diPwSKiEHWvIVvi4l40xQzORuaj7LsEpGe8VXENBuL1zU0oXd0yZMiXhc/y9/sz9Vd/ixvsGuJsZLSkqIh4XayKXzzo6KZ5jRKKVliY+eMHX7XWeW74NX+Uh7mZGJXAR8bg2lxM1xsxI5cKRw89EgL3Gzsdl3TrW4eIsbGEK4CLice2tB17dTlq7QtXwo5M9X/LTq6++mvhJ69Zx2iVpCODl5bB+PSQ4K5yISK5oL4CXGmNuTPSCxpizaXt8tUhi1q1jyvQ0BPCuXaFPH6iPNZOuiEjuay+AA1QaY06K50LGmD6hCV/mkaGZ2qQTWLcufdOdqhpdRDysvQA+zlp7CjC0oyAeKnV/BlQB00NTsT7jXjalU9q5EzZvhrJYs926QGPBRcTD2gzg1tqXQ/8+QhtBPFTqfgmn1P0ZMDQ0Vhxr7aT0ZFk6jfXrneBd0FFFUZI0lExEPCyu34yxgrgx5nKc1cnGAbdYa0dreJm0J97lVndxex3waKpCFxEPi7toExHELzfGvA3MAV7GKXXfm64MSv7o1atXYidkIIC/Ok8BXES8KaG6yVAQLwBGAdXW2kqVuiVeb775ZmInZCCAf/S6AriIeFPCjYuhaUenAjXuZ0fy2V/+8pfETqirS18PdCC49kgWDq92f6EUEZEMaG8mtl+2lRYqiZe21zu9vfNF4pLmEnig5hAunDfQ/YVSREQyoL0SeLszqVlrnwFK2gnimolNUpPmAO4fs51A1Sf4J+y1JICISM5r7zfXqLaWFY1SYoyJtZa4ZmKT1KxbB8cdl7bL+87qzUlTT6Nb5aa03UNEJF06Knp8BiQz12Q/YEQS54nslu5ObH370o0dsG0b9OiRvvuIiKRBewF8ibW2MtkLG2MWJXuu5Kdrr702sRPSHcCNccaCr1sHBx6YvvuIiKRBe23gi1O8dqrnS54pKipK7IQ090IHNJmLiHhWe1OppjQ5iyZ3kWh/+tOf4j+4tdWZSlUBXEQkpjRNMi2yt3feeSf+gxsaoLgYEi21J0oBXEQ8ytXxM8aYmwCL037+rpvXlk4m3e3fYQrgIuJRCQfw0AQtDcAia+0rkWnhanNjzDnGmNHW2l+5k03pdDIZwNesSf99RERclkwVejkwHVhijGkxxrxtjLkrckKX8CQvbmVSOqG6uswEcC0pKiIelXAJ3FpbZYzpi7OM6DjgZOAWYLoxBpw50mtwSukiyVm3Lv0d2EBV6CLiWUm1gVtrNwLzQ19EBfQqnABe7VIeJU/ceOON8R+sNnARkXa50gvdWrvRWjvfWjsVqAD+jtOZTSQ5CuAiIu1yfRiZtbbBWnsLKoFLlEWLEpicL1MBvH9/515Wf2+KiLekcxy4SeO1xYOWL18e/8GZ6sTWrRv07OmMOxcR8ZCEA7gxZkVEz/MT2zm0LIV8SWeXqU5s4FSj19Zm5l4iIi5JpgR+C07p+hb2HEo2wxhztjHmJGPM93DawkWSk6kqdFA7uIh4UjLDyCJ7n49l91Cy6ezZcW2+MeYoa+17bmRUOhFrMxvANRZcRDwopalUrbVLgCWw11Cyk3GGk00MjQ1fAiwCnrHWrkzlnuJdXbt2je/AxkZnqc/i4vRmKEwlcBHxINc6sUUOJbPWHgyUApOBXwFDgXuBBFazkHxz3XXXxXdgpjqwhSmAi4gHubqYSaQ2JntRu7h0LJPV5+AE8H/8I3P3ExFxwV4lcGPMS+m4UaiE/vd0XV9y34IFC+I7cO1alcBFRDoQqwSe7lJyQtc3xpQAk0KbJTjV8TOttTURx1QD9eHrW2vvibpGSunijo8++ogJEyZ0fODatU7HskzRMDIR8aBYAdyEhoGlgyHx8eEzgenW2gbY1fP9HZw29l3BN9Q7HmNMhTFmdmha15TTJQvWrnWCaqaoBC4iHhQrgE9P8z0TnWJ1dOhrSWi7BigxxpSEgvpUa+2o8MHW2ppQkA9LNV0yrbYW9t8/c/fTMDIR8aC9AnhoLe+cERlcQyqABmttQ6h6fWSM0xpCQXhpKumhYXKSYe8vWcuRN30tczcsKYEtW6CpCYqKMndfEZEUpK0XehpNB74T+r6C2OuO14fS6lNMj6m2tpYRI0bssa+qqoqpU+Ordd+4caNradH7IrcTSUtWItex1lJfX09ra2u7x9Uur2Vjt27srKtLNXtxK+3Xj43/+het++3n6nXdes9e1dmfH/QO9PzO88+dO5e5c+dGJ6c0X7RnArgxZiLOJDEzI0rGZezufBapAafDW6rpMQ0YMIClS5fusa+uro7yBObubu/YRNOi90VuJ5KWrHiv873vddy1Iriolf8bfjvNG/bFl6m50AH23Zey5ua0zL/u1nv2qs7+/KB3oOcvZ9q0aUybNm2P/caYlEop6VyNzFWhTmbTgapQMJc8FFjQzMR5FQTezWAnNlBHNhHxHM8EcNi11vhU4BFjTLjtOlav9sjSc6rp4pJAINDhMf5TuxCsWo7/jG4ZyFEEBXAR8ZicDuDGmJLQMK9oNTjTtC4ldrAtA5a5kC4uqqmp6fAY36gN3L38PHy+DHcm01hwEfGYnA7gwFicceDRSoD1oWFkNaHe6HukW2uXpJruyhNIYmpraS3PcPU5qAQuIp6T6wF8CVHj0o0xFTgl5DmhXTOJGFseqlqPDL6ppksmrV1L/yMyOAtbmMaCi4jH5HQv9NBY7yXGmJtDuxqAUcCo8Mxs1to5xpjq0LjuEpypUKdGXCOldMmw2trMzsIWphK4iHhMTgdwcGZGA9qdm9xaOyed6eKO0tLSjg/K9DzoYQrgIuIxuV6FLnlkypQpHR+U6XnQwxTARcRjFMAlt2SrCr1/fyeAW5v5e4uIJEEBXDLm6aef7vigbFWh9+wJXbvCpk2Zv7eISBIUwCVjvvjii44PylYJHFSNLiKeogAuuSVbJXDQUDIR8RQFcMktKoGLiMRFAVxyR2Oj04msV6/s3F8BXEQ8RAFcMmZAR1Xj4SFkxmQmQ1GCO09g0m/HEFzU/nrlIiK5IGYAN8b0yXRGJP9dcMEF7R+QzepzILDlWM59pA+BBc1Zy4OISLx2BXBjzCfGmKeMMWcDFVnMk3RW2ezABvhHruOJqlX4J6RpgsKPP4YHHoBFizTeXERSFlkC/8xaO9la+6y19t2s5Ujy1uOPP97+AVkugfsmdGPWmkvwVbrfshT8n7e4fPIqgn/uAddfD1dcoSAuIimJ/E0V9/rXxpiDXM+J5L3ajtbbXrMmqyVwBg1inx2r3L9uTQ2BZxs59clhBA64DN5+G/7+dxaO/6n79xKRTiMygNclcN5EtzMiwurVsN9+2bt/SQk0N7s/G9u0afiPrufhm/s51fO9ehG8OsBTq48i+ESWer3X1cGPfgRTpsALL2QnDyKSksjGvnHGmIY4z5sM3JeG/EhntmYNnHRS9u5vDBxwAHz5JfRxpx9nl/feg7ffxrfid/h6FO3aH1jWn9PnFRC4oRbf+RluNqipgTFjYMIEGD2a4A0v8Ot7D+XSHxycluYDEUmPyAA+FLgizvO+loa8SGeX7RI4OAF81SoYNsyVy/WYPRumTYMePfbY759QyIM39uH6Zd+D/9wLAwe6cr8Obd8Op58O06fD1VcDEPigkXOv30Tgx6vxVe6fmXyISMoi/9yeba0dHc8XcG+2MizeNWjQoPYPyKUA7oYNGyh66SW4+OK9knyVBSx8oRjf2X3gkUfcuV887riD4MCLqFz0nV3j3f1n9eDX1xXif+NHsG1b5vIiIimJDOCJdIl92+2MSMfuuadntrOQkkmTJrWdaK1ThZ5PAfyJJ9h50klQXt72MVddBXPmwM6d7tyzPatWwZw5BA6+mivuWb9rvLuvsoDnX+qP7+iN8OCD6c+HiLgiMoAPjfcka+0zaciLdODee70dwNu1YQN0775XVXPGHXAAxLNqWgeCi1q5+DejWDD8u+0fOHw4HHQQvPRSyvfs0IwZcNll+P09d3eoi3T77fDQQ9DUlP68iEjKIgP4ZGNMFnsQSb579NFH207MhepzcK0EHghs5cz/G8zz6+LoLvLtb8MTT6R8z3bV1jr3uOkmfJUFLPp90d4d1r76VfjKV+DZZ9ObFxFxReT/4JOBUmPMOQrkkg4bNmxoOzHPAri/52vMv7QW38ktHR9cVQUvvghbtqR83zb96lcwcSL079/+cddcQ/CB96g8s0lzwovkuF0B3Fr7d2vtM6Hq8XeymCfpjPIsgPuW/5QnbvyUk8bEMa96//4Ej5yG/8wN6QmaLS0wezZceWXHx552GoGWcXu0kYtIbmprLvQhWcyTdEa5EsDLy52hVps3J3+NLVvgzTfhlFPiPiUw8BIunFWYnqC5YAHsvz+MHNnxsUVF+A/8kIW3NKZvTngRcYXmQpfckCsB3BgYMgQ++yz5a7z2GowendC65v7JJbw4+Z/4T3Y/gAd/9k/O7vFk3KV737XDmLOySpO6iOQ4zYUuGVNR0c4id7kSwAEqKpzZypK1ZAmMG5fQKT5/H35VPgNf88Lk7xtLQwOB+qO54JdF8ZfuTzjBaUb497/dzYuIuEpzoUvG+P3+thPzLYCPHZv4eX4/BALJ3zeWZ57BX/a32MPG2tKlC/h8TtW7iOQszYUuuSHXAviKFcmdu2aNM4581KjEzz3rLPiv/3ImdenaNbn7R/vtb/FdfTW+c4o6PjbSaafBY485E82ISE5ybS70UOc3E9rcYK39Y4p5kzzz8MMPc8UVbXzEcimADxkCixYld+7LL8OJJ0JhEh3A9t8fDjkEXn01uRJ8tC+/hHffdRYtSVRlJVx2GTQ2QnFx6nkREddF/paZba2Na45zY8zdMXb3s9Y+Ekrva4z5DrDYWrsy6tyTgU+j90v+a2xsjJ2wZYsz1MmlFcBSlkoVerLV52FnneVUo7sRwH/3O6davnv3xM8tKYGjj4Y//tFZ/EREco6bc6EvDQ1DO8pauzEUzIdGTwpjrX0ZcOG3k+SNcOnbmI6PzYQhQ2DlSmhNcEy2takH8LPPdgJ4oveO5Ykn4Pzzkz+/stJ5HhHJSa7NhW6t/TtQDfzAGPO2MWYGsB6oD1WvR6pPJrOSp3Kp+hycKuO+fZ327DgFF7Vy6qkbCRZUOtXgyfrKV6CsDP761+SvAfCvfznvdcyY5K9x0klOCVxEcpKrc6GHSt6TgEk47eF/xJnVbWaodD7DGPMw0M54Iul0ci2Ag1ON/umncR8eWNDM5T/dRmDY9anXJJxzDjyT4npBTz0FkyY5PcqTNXKkM5ystja1vIhIWsSaC/1sFwL5Z9baW6y1ZcDBwC3AEqAUuNlaqx7sndBhhx0WO+HLL50pTHPJoYfCxx/Hfbh/QiHPXvIf/Ee3M997vMIB3CbSqhXBWqf9+9xzU8tHYSEcfzy88kpq1xGRtIg1F/qzbvYgDwXzZ0Jt4tOBMreuLd4yoa3e0KtW5V4AHzYM/vnPuA/3ndTK4x+ejO/aYanfe/hwgl3Hc5qvPqm50YOPfMaUng8T3Pz11PNy8smqRhfJUXGNdTHG9MHpeHY0UAI04HRkW2Kt3RTvzay1G4GNiWbSGFMd+nZo6P7TrbUNUenhdvUKa+09Mc5POl3SbNUq5n3+DaqynS5dHYMAACAASURBVI9Ihx8Of/pT/Me//TYMHgz77JP6vY0hcNj1TLmvicCsZnyViY3hDjy/g9OeOJTArBZ8lSlUoQPBnmfw9DsjmLSoVVOriuSYDv9HhjqjbQDm4ZSgpwI3h7Y3GGOeNMYMTlcGjTHV1to5oa/pwGIiVksLB19r7Xxr7XxgvjFmtlvp4p6f/exnsRO++IKfPpNjJfDDD4cPP4z/+CSmT22P39eF589fgX9CggHYWvwfPMCjN3RzZTGSwPKBnP7EIQSeSeNSpyKSlHYDuDHmaZyg/QxOO3YVTgB/BPg7Tke1SUCNMeYutzNnjKkgqnd8KMiWGWPC07lODe0Lp9ew5zC1VNPFJTt37oy5P/ifw9k4fHhurT990EEEN32dU07fFl++Fi92Z+x2iO/KQ/h13QX4BibwRwTAm2/i6/4GL7xY6kqJ2T+hkGcuWY2/XCsMi+SaNv+HG2POwQlko6y1k6y194bbsq21V1hrR+N0SrsCeBe4xRjj8koMgDM0LVo9ThAvAWKtkdhgjBmbanryWW7fPff0TNelvae5mUDpufx4XmNurT/dpQuBYdcz9b6GjvO1ZQssWwbHHefe/Y1xxoQn2hv917+GSy5xbUy9r7KA3577Or66J125noi4p70/0W8BvhMa3x1TaNjYHGvtKJxAXmmMce1/urW2xlpbGiOpAlga+jfW/O31obRU09Pi3nsVwHdZswb/+ieYXpXAYhsZ4u/zF/7vqp0d5+vVV51Zy9yecnTSJGcylnh7ozc2OgH/oovczcdxx8Ebb7h7TRFJWXu/mYbGmrClLdbaOcaYJcAiY8xl1tr/TT17ewu1WS+x1i4LlZJjTQrTgNPZrSzF9Jhqa2sZMWLEHvuqqqqYOnVqXM8A5dTVxV78bePG9vr47X1e9PGR24mkJSuR6zQ2NlJfX09rxCxjhR98wNj+7+Ef28zRIzfRxmvJiuNHb2RM3QwaR97Zbr6Kn3+e1m99i20xDkrpPR98MCXWsmXBApq/+c0OD+/25JN0O/poNnXtiqsvcuBAyr74gg0ff4wtS2wQiVufMy/r7O9Az+88/9y5c5k7d250cnkq124vgCc8GbS1tsYYUwm8bYyZl0gP9XiE2sSnhkr8WTNgwACWLl26x766ujrKy+P/WbR3bKJp0fsitxNJS1a81/nmN79JWXQA2LIFhgwJNSvkWM3EmDFw22306Oj53ngDHn2U4jaOS+k9T51Kyfz5zupg7bEW/u//4NZbXfu57uFb36Lfhx/CmWcmfGpa8uMxnf0d6PnLmTZtGtOmTdtjvzEmpb+026tCT2q601AnsEdwOre5bSbOhDORYhUJSlxMF5dUVlbuvTMXx4CHjR7trObV3E4b+MqVzkxlySwfGo8LL4TnnoOGDlb6/ctfYP369C08ctxx8Prr6bm2iCSlvQCe5DRQAMzA6a3uGmPMTKLGf+O0g8cKtmXAMhfSJd2++CJ3A3jfvk7e2pvQ5fnnndJxKlOWtmeffeCMM+AXv2j/uAcegOuvT18+1A4uknPaC+BJz5gWmrDFtaWlQu3es0Ol+/C+saFgXhPqTR6pxFq7JNV0t/Ivjvvvv3/vnblcAgenc9rf/tZ2+u9/7wTYdLr1VoK/+FebQ9qCD6/gspqrCQ6+LH15+PrXYflyp6OciOSE9gL4UGPMiSlc25UVx0Id1ZaGg7cxpiRqiNdMIoaaGWNG4sy77la6pNOqVTBoULZz0bb25gJvaHCCe6ymATcNG0bgiGmxh7RZS+B39Uz43WEEXk5jL/4ePeCoo+Ctt9J3DxFJSHsBvARYYoxpMcYsNMZ8zxgzop3jo6VSBQ/s6rS2GHjHGGONMRZnVrjFONXfWGvnsHtc90RgrLV2V/V9qumSZrleAq+sdCZpibU+97PPOktuuj18LAb/5Qew4NwP8Q/7fM+EX/wC/8anePimsvQPw1M7uEhO6eh//MvAaKAy9GWNM0HEYpxS6hJr7bttnJtyFXqo1N3hdUJBOG3pkiYtLc5SogMHZjsnbRs82Fmf+913neU1I/3mN3DDDRnJhm9yOT5TB9MuImhfIfDRYPw9X8f36//G9/rr+A7plv5MHHccPPhg+u8jInFpL4Avs9ZWAhhjhuDMylaJ0ws8noCecglc8tzq1dCvHxQltlhHxp1+Osyfv2cA//hjZ670tlZYS4dJk5wq83u/YPxvexGY0h3fokVwyCGZuf8xx8D558POndC1a2buKSJtaq8K/anwN6ElQR+x1laF1vgeitPL/Fmc1cUqgXtwqrpbjDEv4ZTcRXYZFT3UauVKGDIkK3lJyCWXOKXtyOFk991H8LSfUllFZudwnzwZ/53/j0du7I3/tm/CkUdm7t5lZXDQQfD3NidnFJEMajOAW2vvbSeto4A+Do2lzopcnmd9zJgxe+5YudIJCLnuiCPgkEMI/uCPVJ7ZRPBXK+G55wh0P5Mr7lmf8Tncfb4igi/0Snmp0KQcfzy89lrm7ysie3Flgd8YAT3N3XJzy7amFppbcqPFIJfnWW9qatpzh1cCOMD99xN4udAJ2L/6HGbOxH9WDxbOyr053NMp2G8yFzxzYm6tHCfSSbkSwKOFxlB3mnq2VmvZsmMn25pasp2VnPbQQw/tucNLAXzUKPznlfHMpWvwn9wCl16Kr7KAOQ8WubJsp1cE1n6Nsx8bSGBB7KVhRSRz0ll06FRjqS2waftOdra20rtbIcal5Rzz2sqVTscsj/DdOALfjdnORXb5/T0JXPx3/JcdAHig/4JIHktb0cFae0u6rp3LtjW1sGHrTlpbc6NKPad5qQQugLM++Jyv/gZf04vZzopIp9d56v4yaGdLK+sbm9jZonbCNrW0OPOgH3hgtnMiiTruOHVkE8kBCuBp0motGxqb2L5T7eIxhceAd++e7ZxIosI90a1qmUSySQE8jSywcdtONm9Xhx+Ab33rW7s3VH3uXQcdBIWFsGJFtnMi0qkpgGfA1qYWGrY2YTt5ieWYY47ZvaEA7l3GOKVwzYsuklUK4Bmyo9lpF2/uxO3iW7Zs2b2hAO5tmtBFJOsUwDOopdVSv7Wp4wPz1OzZs3dvKIB7mwK4SNYpgGdYuBa9s3du++zlGgVwLzvsMNi82RlJICJZoQCeJRu37aRxR2bn0M4VwUWt3N77NoK1w7OdFUmWMc5a6Es61XxNIjlFATyLtuxoZlMn7KEeeL6JM+cdSuDtsmxnRVJRWQmLFmU7FyKdlgJ4lm3rhD3U/V+r5dmqFfhP05rSnlZZCYsXO5PyiEjGKYDngB3NrZ1i+tUTTjgBAF+/dzl3+cxOtQhIXho0CPbZR+uDi2SJfoPmiJ0trdRvbeLumT2ynZW0GT16tPPNihWUf+Pg7GZGXBE88rtU3dhXy4uKZIECeA5pabXcf18xO5r3rpIMLmrl4OFNnv1FefvtMH36BmdjxQq+eYECeD4IdD2N8+b0IrCgc3bIFMkmBfAc1LB17x7qgQXNzJy33vVflLff7url2nTHHXDPPY86GytWwCGHZObGklb+qr4sPO+f+Mdsz3ZWRDodBfActWVHMxu37dzVuc0/oZDpVf3wT3B3Cfc77nD1cvH55BM4WCXwfOA7oydzDngQX+Mfsp0VkU5HATyHbd/prC3e0mrxVRawYnmR9zt+7dgBa9bA4MHZzom4xe+HQCDbuRDpdDweDfLfzpZW6hubaGr2Ztv3Xj77zFkDvNDdmgTJotNPdyZ02bYt2zkR6VQUwD2g1Voa8mUOdVWf55/ychg1Cl56Kds5EelUFMA9IjxCvL6xiRbPjhcf53RgUwDPP+edB3PnZjsXIp2KArjH7GxpZf2WHSkvhnLPPT1dylEijoR//QsOPTQL95a0mjwZ/vhHWLcu9Ws1NcFvf0tw/H1UH/dXgmPvhl//2uk/ISK7KIB7kMVZDGXj1p0kWxi/997MBfCtTeGhb7Xw4YcwbFjG7i0Z0qcPnHEG/Pa3SV8iuKiV6otrCR52DTz2GIG+VYz/1WAC5efDU08R/OoNVF9a59m5EETcpgDuYdubW9i0rYnGHc05OZd6S6tlQ2MTm7eHA/jjCuD5bMoUmDMHWpMIsNYSmLWC8d9vJXDsj2HxYvyXDeLhm/vhn3IALFxI4Os/ZPwtOwk8ttr9vIt4kAK4x1mcMePrG5tyZo1xay2NO5pZv2UHTS1Rv8x37ID99stOxiS9TjiBYJdKzjylNrFScmsrXHst/pqHeeTG3vgv3AcAX2UBi36/e+ik/+L9+N8buuF/9fvw4ovpeAIRT1EAzxMtrZaN25ylSd0O5InM1rZ9ZwvrG5vYsqN5V8e7V16BIYc3ceiIc+Gww5y1pCX/GEPgq9/j4p8T/4yBzc1w6aXw3nv4Xr+N4Au92pzrwFdZwIJgGb5nroSLL4bly13MvIj3KIDnkVdegYOHN/Hcizup27KDHc2trlStR87WFiuYW2vZ1tRC3ZYdbNy2c69e8gsXw1HHNjHjyQJVn+c5/0UDWHjpp/j7vd3xwU1NBMffx4UfXUPwppegb9/4bvLNb8IDDxA8exbV12xVm7h0WgrgeWThYpg5bz0LFzsl8q1Nzazb7ATVWAukJCMymDc1t7K1qYV1W3awafvegRvg3hmFjB8H7/25iFvPbVUAz3M+X1fmzGjB96vzYNOmtg/ctg3OOotA6zj8cw8g8HKCE/tceCGBoVcx/tqNWkhFOq2cD+DGmBJjTLUxZl4b6dXGmImhr5vdTo8rjytXQg50Ihs/DqZX9WP8uN37LE61dsPWnTRsbaJha1NEr/D2RZe2w+uVb96+k3Wbd7BhaxM7mlvaffT77y5k6V8L+eyfRXz87iYF8M7ghBNgwgS48MLYHdpqa+GUU6C0FP9Nw1k4K7k5/v1XDeWF8z/Bv++7LmRaxHtyOoAbY0YCY4F6oCJGejVQb62db62dD8w3xsx2Kz0uLS10P+M0how6ir4XnkeP2bMoXP5+cj1xU3TiibBieREnnhg73QI7mlt39Qpfu8lZQWpDYxMbt+3k1h+2sHn7TrbsaGbLjmbuuINd7ep1W3awboszDndrUwut1nLvjPh+6d5/dyGGVmCYAnhn8eCD0NjIX6pmU33VFqeau6UFfvc7GD3aCfJz5+LzFTHnweTm+Ped0ZNH72/F94uzob4+DQ8hkttyOoBba5eFAmtNG4dMDaWHj6/BCfhupXesSxe2vf8BXzwfZMdpZ1D4wQf0veQC+g85gJJzJ9LzoZ9SuOwdp7NOjgkXnJtaWtm+s4W7/6cLW5taaNzRvGs503CHuFjV4/ff3X4ADwf4w/kH7zKCQbwBQ4a49wCSu4qKYMECft/Tz/jrNxO4/S+w774E7/07px/+DsHj7oAuXVK/z5gxzmIq06alfi0Rj8npAN4eY0wJMDJGUoMxZmyq6Ynmp3n//dk++Tw2PzSL9e+8z/q33mFb1WS6rFxJ36uq6T9kf0rOOZOeP7mPr/BRopf3lHDgvv/uLlzOI7zKCTzANL7gLS1i0pl060blZQew8KES/BcNgHffJXDcf3Ppz1rcbbeeMQP+9CdYtMi9a4p4gJd/m1YADTH2h6vb61NMT0nrvvux45wqdpxTBYBZX0fRm3+m6LU/8TInU/b/ytnun8iOsyfSkqVSabjX+iuv0Ga1ezLuv7uQ6Ret5HluYH++5Hhe40MOB+537ybiCSeNaWbSxBLAmf/eP6GVQJJt3m3q1QsefhimTnWGlvXq5d61RXKYlwN4GU6wjdYAlLiQ3qba2lpGjBixa7vVWs70n8NFUy5v+6QuBXDccXDccRw25+f8+4fP0ev5P1Ay9nia9z+AzaedzpbTTgeOZlNDW+15A/dKa9yyuc1j7p9ZxI3T69s4diALXtzJzHlbCDzQi1Ff2xzj+uHtyP17fr/X/ZubuZbnKDv2dpZyHWfzLDsp2pVcX19Paxb6B2TLxo0bs52FrIp+/qNHOl8AdXUu3mj0aN489Er+ULWSyssO4KQxudNkpc+Anh9g7ty5zN17wZ/yVK7t5QCeNQMGDGDp0qW7tht3NPOf2rX0KSmL6/xWutB1/AR2jJ/AjuZmil57lV7PzqffqeN5g8MY8NQ57DjTT+u+e89YFuse0fvC27NndeeOGdvbPHbCqV2ZXtWPO+9s3bW/rWtF7o/8/n8fGcgP7+gK1lK0ZBHbr/wBZzKADS/9kTu+PoJoZWXxvaN8Ul6e0v9Rz8vU878w9GpOvX4TC+9uZtLE3Hrn+gzo+adNm8a0qL4axpiU/oz1bBt4SKxoUOJievoVFtJ00sls+vkvWffxZ9zF9+m67B36fX0kpaedQo9Hf4VZ705R5Wc/6b3Hdke91uPx85/2ouiPL1PqG0vv70/ninX/w1iW0PKVw/Y6dty4GBcQcYn/rB7MvaoZ/1t3OCuaieQ5LwfwpcQOtmXAMhfSM6+oiBeZwKbZ/8u6jz9j6xVXU/T6q5R/7auU+E/nSmZR8MXnSV/+5z/t3fFBEHN42F77tm/nYh5jGSPpfetNbLv0ctb/9R3+wJmAM1Vq9+LW0BSqOzj/4mZ++cuzks67SEd8lQU8t+QAfENXwt13Zzs7Imnn2QBurW0AakK9ySOVWGuXpJqernxHdhxrV/fu7DjtDDb++v9Y9+GnbLvoEr7JX+l3/DGUHft1iu+8na5vv8XP7+/heh4jh4eF8/uLhwp45Y+Wrw1fz18uf5z+RxzKuTzJD/gf1v9lKdsnn7fXsKADKppDU6jW0613M0OHDnU9ryJ7MAZ++Ut46CH45z+znRuRtPJKAG+r4XQmUB3eCE38ssTFdNdFTncat+JidvjP4WLmsm7Fv9l8/08xLS30vv4a7vzZYErOnEDxzLvo+tqr9KHjDiNx/xEBLHyxmZnz1jOm4iNenvE2P5zXxPONxzLznNfwsZAgp0JB7I/RqprC0BSqZZxWWcCnn36awEOLJOmAA+DOOwle8BjV1+3QXOmSt3I6gBtjKkLTm04HRhpjZoZmTwPAWjuH3eO6JwJjrbVT3UpPh1jTncbrmhs2Q5cu7PzmMWy5/U7q33ybg1jJ1iuuwmzZQq8f/4j/MJDyrx5KyaSzmcEt9PjVHIoWvkjh8vc5kH9j6utZ+FKr80fEIktPGilY9QWF77/HBF6gxyOzuZfv0fe8KmoYwuRnLuN3VV9QXrOcyvP3YXpVP8ZVD+HHjx/KwcOb6F7s/HKMrmJ/5RWnBL7634V8/G43KisLeO6559x4hSIdq64m0Pscxl9dr7nSJW/ldC/00Mxo94S+2jpmTgfXSCndbbs7jm3v+OAo1313M9GVERsoo8k3gSbfBAAG9u1K7e8/pPAfH7DlpRUUvv8e3Ra+SMF/vuR1NlI+YjOT7Yk8+8YVTF75S2bxEt3GltBaWsbVDKLwg/1Zy6EsGHEGP/7iSG69o4j5Z/cEvs59l25nxQ1O3g+ocErm06v6AU61+023NtO9uJUDKpp58ml2pa9YXkSvbjn9UZN8U1CA/6qhBC74AP/ZxcDR7ly3pQVWrYItW5zV0/bfX8vjStbot2qeaaULLUMPpmXowfwP3bn2Z7v/UBjctztrPt/OSODUvt35743HU9y3O2s+co45tW931jy4nXsf687FG7dz27wN/OGh7jHvs6qmkOlV/VhVUwjs7vEbDuy/+mFpRLpI5vkml+M7oCv4J8CEJXDkkcldaPNmgre/zhNv9Of8T+7F1+PPBHufxfyeE5m45lZ8Yw1cdZWzzKlIBuV0FbqXRA/R8prodvGOqvq3NxawYnkR2xsL9rhGSwvcem4Z37l473SRjDv2WPj5z53Vz5Yvj/u04KJWqq/cTPCiuTBkCIF393WWPT3/MfjySwKn/IQJvzuMwBmzYMQIguf9L9UnLyP4xNr0PYtIFP12dUm8Q7RyVXTnunjGiHcvbt0j6C9cDPc9u57WFoP/tKK2TxTJpEmT4Cc/IXjKvVSf/2XHndqWLydw/weMv2ELgbqR8NZb+KePcJY9PcOpkfJPKHS2zy6GadMInPYQ42ftR+CBf8DChRl4KBEF8JySSO/wtlxzQ/TUqvFJpnNduLo8HPRbtnRtt9p80qRJSeVNJGXnnktg/AOMv62AwN1/J/jTD/booR6c30D1uZ8TPPZHcMop+Ad/yJwb++C/4XAYOhRfZcEey55Gb/tPL3IC+gX94fLLCV71LNXXN6kHvKSVAngOSWqIWRSno1vi2itxR5e0w8Lt4GeOdz5Gj/+my17V5rfd5nwBDBo0KKm8ibjBf24ZC39eiv+rXxJ4eqPTQ/2/3oCBAwnc9Q7j7+hKYNB34LPP8M2ZzMIXiuNep3xXQL/hq/DWWwT+vg/jr1qvHvCSVuphlEPCpeA778zcX+033tLxL5hwSfsPD5XusT/cDu4/re1zb7999/cffvghw4YNSzKnIqnxVRbgq+wOnAGLWgk8tB3/VRVw8t/w/3M/ArNa8E8phG4plmv23x//jcU8ecE/OffQhWDb+Q8ikgIF8BySyhCzZN10654BPFY1frikfeedrfTv3cybrxVw8PAmVtUUsr2xgOCiVg4e3rxrOyxc8g578cUXFcAlJzjBvCfQ09k+AHyVXdo/KZHrTyzBd9JhMOYKGh/42FmzXMRlqkLPU9GBON729VjV+Fdf28qK5UWc6evKPXd1IRg65riTneAfWOCU0A+ocLbDgTuy9C3S6ZSVwaJFvPrUeqrPqVF7uLhOATxPRQfieNvXozuzFXcr5O7/dj4mPYq6YIzBP8EpkX/3aqcCJ7wd7rymwC0Ssu++zDtxBuPv6kHg0VXZzo3kGQXwPBUdiKO32+qtPn6c065ddXo3brsNenUrpLDLnh8TX6VzTGSP3FU1hRxQ0czFU1TKEInk83Vh4X098P/5R/CHP2Q7O5JHFMDzVHSv8ujt6767GQMUdSmguFshpT2dcdsloX+7FJi4StI33bQV2N3RraiXet2KRDppTDNzHinBF7gGLr8c/vSnbGdJ8oQCeCdREJqvuWdRF/r26Eqf7l3Zp093SoudecqLCgv26nQWj5tvdgJ4uKObf0Lb/SIvuOCCpPIukhdGj4annnImllm6NNu5kTygAO5hkR3TDE6pubCggB5FXXYtHlJWXMQ+vbvRv3c3AHp370r3rl3oUrD3AgyptF2Hh5S1N252wIAByd9AJB+ceCLB6/7AZdWbCT76ebZzIx6nAO4BxkBhgaFrFyc4F3crpE/3rix52TBz3nr++McC9unTnfJe3ejd3UkrDgXwrl0KMC6slhRdOu9oO5b3338/5XyIeF1g3UhnHvWHP4OVK7OdHfEwBfAcYXCCLTgdx/r26EpZcRH9e3Vjn97d6derG71CgbtXt0J6FHXhnNO6dlhtnYxYwTi6dN7RdiyLF6cwxZxIntg1j/qxW+C44xJaZEUkkiZyyZKuXQoo7GLoWhD6N6Knd3Gca2fv7g3ubt40DEwkfZxJZIqACfCNewmOv4/A2Hvxf7s87qlbRUAl8IwpLDAUdyukrNjp5V1WXESf7l3pUdRlj+DtZcl0ghPp1M49l8DYexj//RYCsz511uMViVN+RI4cVVhg6NWtkH7FRbuqwDMZrMNDvNradptK7iKJ83+7Pwt/0gv/f34NxxyjHuoSN1Whu6xrlwK6dy2gW2Hsnt6Z5Azx6tnmtohkn1Ol3hta/xvmHkrQ/3N+u/81fPss8E07EoqKsp1FyVEK4C4IB+ryXt1SDtpOKTk/g+yUKVOynQWR3FVQAJdcQuCdczn7mg0ELl+Jb+YAgodfS6DkXPyDPsR31FqCHw8isPIw/CVv4StZSrBmCIHG/4d/+++hcRtPF5zBpMan8JX/nWC3CQR6noV/+Gp8lw2GigpnWIvkBVWhu6B7V2cVIzdK3OGJUfJRaWlpxweJdHL+04ucXur/9Q346CMCQ65k/AOlBFYdDu+9R+Cjgxh/dzGB+qPhwAMJdDuT8b/Yj8CQKwmMuInTnxhK4P/9GGbMIDDgIsbfV0Lg7VI47jiCFdVU+z4g+JPl0Kppj71OJXDJmKVLlzJ69OhsZ0Mkp+3upQ4wAP9FrQRmNeO/bgBUPox/UWj76gFQ+V38R4S2L3R+nQdmNeOfVAjHHYd/WyjtpgEw7ksC1Q2M/952ApfW4HtgMMGT7iLQ/Uz85/RSD3gPUgD3OC9Vub/66qsK4CIJ2jOgx7Pddpq/qi+BWcX4bx8A+71I4JYCxj/QSODKT/C1roPKSqcqXzxBAdzj1DFNROK1Z0Afjv/6VgI/34F/9Ptw660Er3mOwPAb8Z9bhm9yv6zmVTqmAC4i0kk5Ab0HcArYSgJT1jP+lp0Ezn8P3xMPEvz6rQTWjMB/evvrHEh2KICLiAgYg/+8MqfN/LZvwLozCPwGxj+ygcDVa/DVfgDHHgtDhqgne45QAHdBcFErBw9vJrioUH+liohn7a5iLwIuwz+olZ/duJXrRi2H3/+e4A9eYX6/C5lY/Crsty+BlnH4h32O7xubCX46hMD7++EfvR7f1zcRfLuEwNIyKo/4komnrCG4tC+BN/vgP7EJ31m9oEePbD+u5ynauCCwoJmZ89YTWNCc7azktKlTp2Y7CyKSAF9lAcEXeuGbWQnz5xPwP8yEJw8jcNBUAraS8TN6EvhbKTzyCIHAVqf6/Xfr4corCfx2HeOnN/HSM5th8mQCD69k/I1bCTzwDygtJVh+LpeOeo3gqOkwZQrBK5+h+oL/EPzdumw/tmeoBO4C/4RCrqvqx8/u0+tsT69evbKdBRFJgX9CIYFZ/fBfFDFk7aYBUPn87uFtPxoAlUt3bY+b1gMm/gP/olZ+cnMz3/2f42HcNgLXbuP0azcSuOtGfN98jsDvhzH+AUPg2x/g+94FBIddQ6B3Ff7jt+GbOhR6qrNuNEUcF6RrVbB88+abb3LMMcdkOxsikqREhqyFt+vqNsVM95/RncCsQvzfkuFTSAAACaVJREFULoTKavwHhf4AuOt4OOQNAt8vZPzthQQuWwf/fSlPVtzAuUP/gW9iidMWv99+GXji3KYALhnzl7/8RQFcRICOxrMPwX/p7hJ94PnfcOY1GwhML8f32HcJXh0gcFA1/n7v4Du9B8EupxBYPrDT9ZZXABcRkZyzZ0AvcqrurxoAlS8QuH4H46+qJ/BfB+H72+0EPvoq4x/rQeCSFfh++QDBXmcR2HY8/qPWOAF9332hTx/o3RsK8yfs5c+TiIhIXtqr+n1CVyegX14Ilb/Gv6iVn97UxA2XDYKSbxN4/CjGz+xO4MoifM9fSbB2OM+Uns85NT/H1/wSwbIqnh5wGZOansVX8g5BxvJ48+lcUP469OjBY6uP55Jh7+E7fDXB/wxj9juHM/Xk1fiO2Qq9ekFZGQwYAMXFBBe1EljQjH9C5kchKYCHGGOqgfrQZoW19p5EzvfSlKYiIl4Wu/q9O3AQcBD+PqEOc7cMgMq/Ebi+iVOvWk9g1lP47tpJ4HuW02/YTOCe/8J3yT8JzDmYc34IgR+dBS0tTP5VDwI3HYWvx5sE/nUSF80qJHBlAb6XbiNYP5J5RWdRVXMzvp0LCRzxIuPnHkLgyv/ge2aOE9j32cf5d8AAunTtCl27Qt++rr8HBXB2B29r7fzQdoUxZra1Nu5xT5rSVEQkN+xdYg/1np9QCMVF+P2tBGYV4Z8ctejLZRG9668ZAJUP7u5df8sAqHyVwPVNnHbVegKznsB353b8gU3Muq6Qq8bsgP4jCC4t4Td/HM7F3V+Ajct4utXHpP1+je+FG1x/TmOtdf2iXmOMecdaOypq36fW2qGxjh89erRdunTpHvvq6uooLy+P837Q1mtv7zqx0qL3RW4nkpasRK7T1NREUVFRxwfmEbfes1d19ucHvYN8e/6Oqsyrr29i/FXrWTjLmUs+/P2cB/f+3ReKPUmv8NTpS+DGmBJgZIykBmPMWGvtkkznKV91tuAtIvknunQfbY/SPnDfjaV877r0hNpOH8CBCqAhxv76UJq45E9/+hNjxozJdjZERNImOsAfPXJL2mogFMChjN2d1yI1ACWxTqitrWXEiBF77Bs9ejR33313XDc8++zV1NXFnoRg8eLFjBs3Lu606H2R2+2lPf3000yaNCmu/LanvfxGe+211zjyyCNpbW1N+b5e4dZ79qrO/vygd6Dnd55/7ty5zJ07Nzo5tchure3UX8BY4NMY++cBN8c6Z9SoUTbahRdeuNe+trR3bKJp0fsit9tLO+qoozrMZzwSee777rvPlXt6iVvv2as6+/Nbq3eg52/7+YGlNoX41XmmrGlfWYx9MUvfIiIiuUABHJYSO1iXAcsynBcREZG4dPoAbq1tAGpCvdEjlVj1QBcRkRylceDsmsilxIZmXzPGjASm2jYmcjHGrAP+HbW7L7Axzlu2d2yiadH7IrfbSysH6uLMb3sSeW6AQ4BPXLivV7j1nr2qsz8/6B3o+dt+/sHW2v7JXlgBPCQUxGtwqtMTnkpVREQkkxTARUREPKjTt4GLiIh4kQK4iIiIBymAi4iIeJACuIiIiAcpgAvgrMpmjKk2xszMdl5ERKRjCuASFl6TVlPIioh4gAK4ABCada4m2/kQEZH4KICLiIh4kNYDzyOh+dwnAeOstVUx0qvZvfa5ZptLkjFmLDAVmIGzbvxEoMFaOyfimHbftZd+Fql+rrz+Ltp7/s7yWQjlEWAoTjPb9NA6EpHpefsZgPbfQdY+B6msRaqv3PkCRoY+NBOBd2KkVwMTI7YrgNlRx4yN3qevmO96IvAOYIENwMxE3nU8P4tc+Ur1c+X1dxHH8+f9ZwGojvHMn3aWz0Cc7yArn4Osvxh9ufsV+oUT6xdNrH2fRm0rgMf3jid2kN7uu47nZ5FrX8l+rvLlXbTz/Hn9WQgFkpkx9m8IP3u+fwbifAdZ+RyoDbwTCFUBjoyR1BCq+hGXdPSu8+lnkeqz5tO7iCWPnr86xr56oKwTfQbafAcdnZjOd6A28M6hAqddJlp9KC2yDafCGFNtI9puZG+h91WC815H2t3tVR296/oO0r0k1WfNi3eRz58Fa20NUBojqQJYSif4DMTxDoDsfA4UwDuHMnZ3jojUQGjct3WGkS3JZKY8bBns+o+NMabeGLPYWjuOjt91hz8LD0n1WfPhXXS6z0Kos9USa+2yUNDqdJ+ByHcQ2pWVz4ECuEiCwv9JI7aXGWNGG2NyosQgmdPZPguh55pqrR2V7bxkS6x3kK3PgdrAO49YbTU5+xeuB9XgdAKEjt91Pv0sUn3WfHoXYfn8WZgJnBy1r7N9BmK9g1jS/jlQAO8clhL7w1BGqOpH4mOMqTDGbGjnkI7edT79LFJ9Vk+/i872WQitk7DH+G862Wcg1jvI5udAAbwTCH3YakK9HSOVhNq+JTEzYuyrwGkTa/dd59PPItVnzZN30Sk+C6E239mRVcXGmLGd6TPQ1jsIfZuVz4ECeP5pa1jDTCKGQhhjRqJOawkL/efdo8eoMWYi8HTEf+yO3rUXfxbJfq7y5V3s9fyd5bMQClJLIzpolUQNb8r7z0B77yCbnwMTGjAuHhfqLDERGIfT7nIPzkQA0VP51eBU1+TkdIVeETGtYrgXf6xpEdt81175WbjxufLyu0jg+SEPPwuh5/+0jeRSu3sq0Xz/DMT7DiCDnwMFcBEREQ9SFbqIiIgHKYCLiIh4kAK4iIiIBymAi4iIeJACuIiIiAdpLnQRybjQpBUzccbPvm2tnZ/lLIl4jkrgIpINLwMzrbXTgaOznRkRL1IAF5GMMsbMw5liMjxLVS4vXCGSs1SFLiIZEzGzWWlou4SoaShFJD4qgYtIJs0E5kSs5jQJeCqL+RHxLAVwEcmIiNL3vIjdQ621ObdspIgXKICLSKZMBRrCSySGl2fMbpZEvEsBXEQypZrQEomhtu+SyLWVRSQx6sQmImkXqj4vARaHds201k7NYpZEPE8lcBHJhLGhf5cYY2YD07OZGZF8oAAuIpkwCme42FRgdkQvdBFJkgK4iGRCGU4V+lPqdS7iDgVwEcmECmC+greIexTARSStjDE3h76tz2pGRPKMAriIpI0xZiwwH2e2tbI20kUkCcZam+08iEieC437fsdaOzRi32yc4WQaCy6SBAVwEckIY8xEYBy7Fy+Zod7oIslTABcREfEgtYGLiIh4kAK4iIiIBymAi4iIeJACuIiIiAcpgIuIiHiQAriIiIgHKYCLiIh4kAK4iIiIBymAi4iIeJACuIiIiAcpgIuIiHjQ/wderxkTbqJNwwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "label_size = 18\n",
    "title_size = 25\n",
    "legend_size = 18\n",
    "handle_length = 1.5\n",
    "xscale = 'log'\n",
    "yscale = 'linear'\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "fig, (ax1) = plt.subplots(1,1,figsize=(7,5))\n",
    "ax = ax1\n",
    "\n",
    "\n",
    "\n",
    "verr = err_p_Dl_low_ell\n",
    "v = Dl_low_ell\n",
    "verr2 =err_m_Dl_low_ell\n",
    "verr2[verr>=v] = v[verr>=v]*.999999\n",
    "\n",
    "ax.errorbar(low_ell, Dl_low_ell, \n",
    "            yerr=[verr2,verr],elinewidth=1,capthick=1,capsize=0,\n",
    "             label=r'$\\mathrm{Planck\\,\\,data}$',color = 'b',ls='None',\n",
    "            fmt = 'o',\n",
    "            markerfacecolor='lightblue',markersize=2,linewidth='0.2',markeredgewidth=0.5)\n",
    "\n",
    "\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlim((1.7, 30.))\n",
    "ax.set_ylim(-100,6e3)\n",
    "\n",
    "divider = make_axes_locatable(ax)\n",
    "axLin = divider.append_axes(\"right\",size=4.5, pad=0, sharey=ax)\n",
    "\n",
    "\n",
    "axLin.errorbar(high_ell, Dl_high_ell, \n",
    "            yerr=err_Dl_high_ell,elinewidth=1,capthick=1,capsize=0,\n",
    "            label=r'$\\mathrm{Planck\\,\\,data}$',color = 'b',ls='None',\n",
    "            fmt = 'o',\n",
    "            markerfacecolor='lightblue',markersize=2,linewidth='0.2',markeredgewidth=0.5)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "axLin.set_xlim((30., 2550.))\n",
    "\n",
    "\n",
    "ax.set_yscale(yscale)\n",
    "ax.set_ylabel('$D_{\\ell}^\\mathrm{TT}\\,\\,\\,\\,[\\mu\\mathrm{K}^2]$',size=title_size)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "if yscale == 'log':\n",
    "    axLin.yaxis.set_label_coords(-.5,0.5) \n",
    "if yscale == 'linear':\n",
    "    axLin.yaxis.set_label_coords(-.5,0.5) \n",
    "\n",
    "\n",
    "\n",
    "ax.tick_params(axis = 'x',which='both',length=5,direction='in', pad=10)\n",
    "axLin.tick_params(axis = 'x',which='both',length=5,direction='in', pad=10)\n",
    "\n",
    "ax.plot(Plc_BF[:,0][Plc_BF[:,0]<30],Plc_BF[:,1][Plc_BF[:,0]<30],color = 'r',ls='-',lw=1)\n",
    "axLin.plot(Plc_BF[:,0][Plc_BF[:,0]>=30],Plc_BF[:,1][Plc_BF[:,0]>=30],color = 'r',ls='-',lw=1)\n",
    "\n",
    "\n",
    "\n",
    "cosmiv_variance_low_ell_low = np.asarray(Plc_BF[:,1][Plc_BF[:,0]<30])-np.sqrt(2./(2.*np.asarray(Plc_BF[:,0][Plc_BF[:,0]<30])+1.))*np.asarray(Plc_BF[:,1][Plc_BF[:,0]<30])\n",
    "cosmiv_variance_low_ell_high = np.asarray(Plc_BF[:,1][Plc_BF[:,0]<30])+np.sqrt(2./(2.*np.asarray(Plc_BF[:,0][Plc_BF[:,0]<30])+1.))*np.asarray(Plc_BF[:,1][Plc_BF[:,0]<30])\n",
    "\n",
    "ax.fill_between(Plc_BF[:,0][Plc_BF[:,0]<30],cosmiv_variance_low_ell_low,cosmiv_variance_low_ell_high,alpha=0.1)\n",
    "\n",
    "#unbined cosmic variance\n",
    "#cosmiv_variance_high_ell_low = np.asarray(Plc_BF[:,1][Plc_BF[:,0]>=30])-np.sqrt(2./(2.*np.asarray(Plc_BF[:,0][Plc_BF[:,0]>=30])+1.))*np.asarray(Plc_BF[:,1][Plc_BF[:,0]>=30])\n",
    "#cosmiv_variance_high_ell_high = np.asarray(Plc_BF[:,1][Plc_BF[:,0]>=30])+np.sqrt(2./(2.*np.asarray(Plc_BF[:,0][Plc_BF[:,0]>=30])+1.))*np.asarray(Plc_BF[:,1][Plc_BF[:,0]>=30])\n",
    "#axLin.fill_between(Plc_BF[:,0][Plc_BF[:,0]>=30],cosmiv_variance_high_ell_low,cosmiv_variance_high_ell_high,alpha=0.1)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#multiline(xs, ys, np.log(1e10*param_values),ax=axLin, cmap='jet', lw=2)\n",
    "\n",
    "#ax.legend(loc=2,ncol=2,prop={'size':14},handlelength=handle_length,frameon=False)\n",
    "\n",
    "axLin.set_xlabel('$\\ell$',size=title_size)\n",
    "axLin.xaxis.set_label_coords(.4,-0.1) \n",
    "axLin.spines['left'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "ax.xaxis.set_ticks_position('both')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "ax.tick_params(axis = 'y',length=5,direction='in', pad=5)\n",
    "\n",
    "axLin.yaxis.set_ticks_position('right')\n",
    "axLin.tick_params(axis = 'y',length=5,direction='in', pad=5)\n",
    "\n",
    "\n",
    "\n",
    "axLin.xaxis.set_ticks_position('both')\n",
    "\n",
    "\n",
    "ax.axvline(30,ls='--',alpha=0.6,c='k')\n",
    "\n",
    "plt.setp(ax.get_yticklabels(), fontsize=label_size)\n",
    "plt.setp(ax.get_xticklabels(), fontsize=label_size)\n",
    "plt.setp(axLin.get_xticklabels(), fontsize=label_size)\n",
    "plt.setp(axLin.get_yticklabels(), visible=False)\n",
    "ax.grid( b=True, which='both', alpha=0.3, linestyle='-')\n",
    "axLin.grid( b=True, which='both', alpha=0.3, linestyle='-')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "\n",
    "#figure_name = 'bard_planck_data.pdf'\n",
    "#plt.savefig(path_to_notebook + figure_name)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define parameter X and range (\"one by one\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "param = 'h'\n",
    "col_label = r'$h$'\n",
    "param_value_min = 0.5\n",
    "param_value_max = 0.9\n",
    "param_spacing = 'lin'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open bard_parameters_varying_X.ini and run class for many different values of parameter X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing for h =5.0000e-01\n",
      "Computing for h =5.4000e-01\n",
      "Computing for h =5.8000e-01\n",
      "Computing for h =6.2000e-01\n",
      "Computing for h =6.6000e-01\n",
      "Computing for h =7.0000e-01\n",
      "Computing for h =7.4000e-01\n",
      "Computing for h =7.8000e-01\n",
      "Computing for h =8.2000e-01\n",
      "Computing for h =8.6000e-01\n",
      "Computing for h =9.0000e-01\n",
      "--- 16.365 seconds ---\n"
     ]
    }
   ],
   "source": [
    "N_param_values = 11\n",
    "param_name = param + ' ='\n",
    "str_param = 'varying_' + param\n",
    "\n",
    "if param_spacing == 'log':\n",
    "    param_values = np.logspace(np.log10(param_value_min),np.log10(param_value_max),N_param_values)\n",
    "if param_spacing == 'lin':\n",
    "    param_values = np.linspace(param_value_min,param_value_max,N_param_values)\n",
    "\n",
    "\n",
    "param_file_content = []\n",
    "with open(path_to_bard_class + 'bard_parameters.ini', 'r') as file:\n",
    "    param_file_content = file.readlines()\n",
    "    file.close()\n",
    "\n",
    "index_root = 0\n",
    "i = 0\n",
    "root_name = 'root = output/bard_'\n",
    "for line in param_file_content:\n",
    "    if root_name == line[:len(root_name)]:\n",
    "        index_root = i\n",
    "    i+=1\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "i = 0\n",
    "\n",
    "for line in param_file_content:\n",
    "    if param_name == line[:len(param_name)]:\n",
    "        index_param = i\n",
    "    i+=1\n",
    "\n",
    "\n",
    "\n",
    "new_param_file_content = []\n",
    "\n",
    "start_time = time.time()\n",
    "for p in param_values:\n",
    "    new_param_file_content =  param_file_content\n",
    "    new_param = param_name + str(p)+'\\n'\n",
    "    p3 = \"%.4e\" % p\n",
    "    print('Computing for ' + param_name +str(p3))\n",
    "    new_root = root_name + str_param + '_'+str(p3[0])+'d'+str(p3[2:])+'_\\n'\n",
    "    new_param_file_content[index_param] = new_param\n",
    "    new_param_file_content[index_root] = new_root\n",
    "    with open(path_to_bard_class + 'bard_parameters_' + str_param + '.ini', 'w') as file:\n",
    "        file.writelines(new_param_file_content)\n",
    "        file.close()\n",
    "    str_cmd = './class bard_parameters_' + str_param + '.ini'\n",
    "    os.system(str_cmd)\n",
    "print(\"--- %.3f seconds ---\" % (time.time() - start_time))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot several quantities for all these values of parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/ipykernel_launcher.py:188: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAFiCAYAAAAJAoquAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd1wUZ/7A8c9soQosHQEBQayAXWOssXeNUVNMjClnYkzP/TTJ3SW55O6SaHKXchej6eqlqEksSRQl9qhRBMXeUATpLEtbYNv8/lj0LIDALizI83695vViZ2af+T4T4pdn5imSLMsIgiAIgtCyKBwdgCAIgiAI9ScSuCAIgiC0QCKBC4IgCEILJBK4IAiCILRAIoELgiAIQgskErggCIIgtEAqRwfQnHh4eMhOTk7X7FOr1QQFBdW5jOLiYjw9PRt8Tk3Hqtt//b6rP19/rLCwEG9v7zrXoy7qUter5eXl0blzZwoKCuwaR3PWGPe9JWnt9QdxD0T9rfUvLi6mqKjommNarbZYlmWvBhcuy7LYqrbevXvL13vggQdu2Febupxf2zk1Hatu//X7rv58/bHu3bvfNK76qu+9WbFihd1jaO4a4763JK29/rIs7oGof831BxJlG3KWeIQuNJn777/f0SEIgiDcMkQCFwRBEIQWSCRwocmsXLnS0SEIgiDcMkQCv4kxY8bY/fzazqnpWHX7r9939efrj82ePfumcdVXfe9NTk6O3WNo7hrjvrckrb3+IO6BqH/j1V+SxWImV/Tp00dOTEy8Zl9+fj5+fn51LqMu59d2Tk3Hqtt//b6rP9d2zF7qW+a7777LCy+8YNcYmrvGuO8tSWuvP4h7IOpfc/0lSTooy3KfhpYtWuCCIAiC0AKJBC4IgiAILZBI4EKTadeunaNDEARBuGWImdiuotVqb+hwMGjQIKZNm1bnMq6faae+59R0rLr91++7+nNtx+ylvmUOHz4crVaLxWKxeyzNVWPc95aktdcfxD0Q9bfWPz4+nvj4+OsPN3wWNkQCv4aPjw/Lly+/Zl9DOmDU5fzazqnpWHX7r9939efajtlLa+6cUlet/R619vqDuAei/n7MmjWLWbNmXbN/xYoVNv11Ix6hC03m888/b/RrbNxsYfQUAxs3t55WviAIrZNI4EKTKSwsbJRyS1NTKRp/B2ZfT358aTuPLyrgxw2VjXItQRCE5kIk8Juo0GodHYJQA7PBwE+D+3ApKoqcjdtZqi2hW+4HrJt3gu4bXuLxR/ZjMomWuCAItyaRwGtx6JNPWB4Tw/mEBEeHIlxHlmW+7xuHx+6DOAMJCsgDdBnrmJU0grvjvmTo+XeJ6buPigqTo8MVBEGwO5HAa5H22ssYjUZWTRhP4fnzjg5HuMr+F/+PvJRT6Nyn8MqgTZxY9QGHh7+MBOwtAu25Iiad/YmIHlomTt6FmHFQEIRbjUjgtRh4dzkRSqgwGPl6QH8MpaWODqlFi4yMtEs5pbm5bH3nXUYDP3Sax9RP44j/9FF6D5/B4AetS5ZuOgGqIRXcnbietPL+vPjyRbtcWxAEobkQCbwW61+ezNQp4C1Bfk4eP4wehdyKxjDb25133mmXcnY/PBsfCwQDI+/+mr/8zR2XvFRefrk7Q5Z9hqaNK4UyHN5qYXebmbz9aRHrNrtRWmqwy/UFQRCaA5HAaxGt+JSkDwdwV29wBk7t3cf2p59ydFitWkVREUkb47kdqIyEPfO7c/xQMRvXtUOSJJROTgz925sA7M+GqZr/8P3cc6S5efDmm4m1Fy4IgtCCiARei8E+bhwsXknJ5+2Z2A4kYOd/PuL4V185OrQW6eOPP7a5jORX/4zaAh2B4n8H8c3aB7i/ZxohIZ5Xzol94gk8XJ3RytDx3Fo+PzGczncf5Z//saDXG22OQRAEoTkQCfwm5vh5stR1PQGfeDO4Kkf8+PBDZB886NjAWqCysjKbvi/LMklffEYMUOYJnw59mKJ3TLyzuPs15ynVano9+BAASedBHgILT75FRXRPPvssxaYYBEEQmguRwOvgXx1iWBC9kj5vqumiBpNF5puhg9Hn5zs6tFYla98+8ovL6Qmo7lWzMvlBhvmeJzCwzQ3n9nzpJQDOmcCkNzLhx02opyhYtDhP9EgXBOGWIOZCv0pti5ks9ujHn6Yu5h9nnqPwA5nssnL+27cXk/f8jlKtvnK+WMykZmVlZTYtZpL06ssUuE/hz1HziBy2k/QlwXz0vIH86v6QcnMjJDyUS2kZHE+GWC8T40J/Zn1FHxISjtOzZ2CDYqgvsZBD664/iHuwePFili5dSkJCAhEREY4Op8k15mImyLIstqqtd+/e8vXy8vKu/FxhssiL8+bLeRORF4P8GsjrJ42v8fya1HZOTceq23/9vqs/13bMXupb5jvvvNPga1ksFvk9V7U8PXaT/P3JTPn2J9Jkz8hDssViqfE7v//97/JrIC9XIlfMUMrr5o2TmWiQH3xwS4PjqK/GuO8tSWuvvyyLe5CXlydrNBpHh+Ewtf33BxJlG3KWeIReD85KiZlu77Fr2RjujAUlkLThFw688bqjQ2sROnfu3ODv5h46hK7cyJTUJfz0/FHOBrrz6FQDkiTV+J0uc+YAkGYGs97M8HU7cBphZPUPFjHFqiA0kR07djBy5EhHh3FLEgm8nsLcVLgrVlP0TVfGBFn3bXrlVS5s3OjYwFqACRMmNPi7p95/FzdgZtk6Zv7n7+QdVPL009G1fscjOJi2wUGYgXPHQOViYnDb7ejdurJ9+4UGxyIIQt3t2LGDvn37OjqMW5JI4A0wOtCDbc6bafulP33cwAJ8N3kihefOOTq0W9aZjT/REdB7wDeae/A9cYHwcM1Nvxc91Tp5zLlM2BwxAZevQ3DpEMDy5WcbOWJBEABSUlLQaDQkJCSwaNEili1b5uiQbhkigTfQCx1CeLPdzwx824UIFVSYLHzdtydGG4dK3co++OCDBn2vsqSEzNwiogHlICVbD47m3rF1G8/d4QFrp8QLRtggPcbDiwJo52th3Xqj6I0uCE1gx44dREZGMnLkSBYsWMDbb7/t6JBuGSKB2+A/Xfvy/NjlTJ4rWadbLSzh14ljxHSrNTAaGzaJyvnvv8cCRAIZ88PI+s6fxx6LqNN3Q/r2xVmtokiGEYVL2HjvEfR9dRRLHUhKympQPIIg1E1SUhJxcXFX3oEnJSXZbU0EQSRwm30WNYO3X/4r00dYp1s9d/QUB//ykqPDuqWkLv+EQMCsgJ8GT0C9O5uuXf3q9F2FUkl4XCwAERfX8b7bZO5quxLC2rF+vXiMLgiNKSEhgalTp17zedSoUQ6M6NYiEriNJAleCfwzm768j5Htrft2LF6MqbLSsYHdQtIOJhEBGP1hQ8FkhrTX1tr7/Hrtp1jfg6fpQNnexJjkzdDLie+/b93jcwWhsW3ZsoUhQ4Zc+fzdd98xffp0UlNTHRjVrUMkcDtwU0kM8vgC5bJO+Cmg1Cizb/4fHB3WLaGiqIjc4goiAcVwNSlbe/DIg971KqP9lCkAXDKB0mSm12+HUPc3cey0BzpdRSNELQgCQGJiIt27/2+q49TUVCIjI0lISHBgVLcOkcDtpKOnExsD1zKkh/Xzni9XUFlS4tigmpnY2Nh6f+fi+vVIQDhwbm40xWsVjBsXXq8yAmJicFarKAG0F8HvrJaQqIvg01EMJxOERqLT6W4Y/z1z5kzWrFkjxoXbiUjgdvRybGeyPhhIkBLKzbBr9kxHh9SsjB49ut7fSfv2K4IAgwJ+6TUGv/PZeHg416sMSaEgtGsXAC6eA0Osiom5P0G4Dxs3Xqx3TIIg3JxGo2H16tXX7Fu6dCnTp08XHdnsRCRwO0v1+IIht1vfzyau34S+oMDBEbVs6Qf2EwoYvCE+dxyje+kbVE7YSOsfDxmFoIgyMzx5K8Sp2bhRPEIXBKFlEgnczqYHe5Py3hTC1FBpgV/vGu/okJqNd999t17nW0wmsvNLCAVUvZQc2RXH/ffVrff59ULHW/87ZBlBpTbRd3cS6p5G0rO8yc0VY/cFQWh5RAJvBMEhXzF4pBKAIzv3U5KR4eCIWqacxESMMoQBBfcFUbRBxZAhoQ0qK7RfPyQgT4bKXAhKycE7Kh98Itm1K82ucQuCIDQFsZzoVWpbTrSuioqK6OLlxXdvPkTUr59yzgC/TLiDEb/uveac2r5f1/23+nKiZ1Z8gRvgBvwwbjDe7+ag18voG/YUHd8AP/Jz88k8DSEd1AzTb2dV0HR++ukgQ4cGNKzQm2jtS0m29vqDuAei/o23nKhI4Ffx8fFh+fLl1+zLz8/Hz69+j239/PwY6/xvmPol51aZOJNyljFaLZqOHa85p7bv13X/9fuu/lzbMXupT5nu7u74+PjU+fyC37ahc5/C49HzuLhdxdAuepvqENq3P/k//0xmLoSOkBlyZgerut3L7t1yo9ybyxqz7JagtdcfxD0Q9fdj1qxZzJo165r9K1assOmvG/EIvZFEeTiz6k8v0tkVzMDPE+5wdEgtTta5C5yNnMfEb+PI+bk3d9/laVN5oWPGWMstBaWfif7JB5A6y5xNdaGkREy8IwhCyyISeCOa3+U1+t7vggScO5tJ3sGDjg7JoXr37l3nc43l5eSXGpmQuoQNzxzlgknByJENe/99WfCgQQDkWEBRYSH691RcY8vBI5L9+y/ZVLYgCEJTEwm8EXmqlXz15LvEtAEZ+Hlq654DeNiwYXU+N3vPHmTg7rJ1zF/4HMbzhfj6utl0/YCYGJQKBToZ9FngptUT5n0evIPYvTvdprIFQRCamkjgjewvMfOIm+eBEkjLKOTS1tY7haDBYKjzuZk/rsYbMEuwLW4YXVQ6m6+vVKsJCA0GIOssVMY6MTIzAaKc+fVX28sXBKHuOkgSwTZsY8eOdXQVHE50YmtkKoXE1w8up+cnd5Kkg033TGHK8dY5bOnDDz/khRdeqNO5WbsSaAvoXSBBN4pJd9R98ZLaBPXsTdbFDLLy4VTX8aR/OwwXLyVJ8SDLcr0WSREEoeEqgGdt+P6q/Hx7hdJiiRZ4E3il21SiF/ijBjLy9KT/+KOjQ2r2si+k0xagLaQc7M6UKf52KTd4+HAAcvSwQZ7L7MVtaYeRMkMAZ89q7XINQRBuTgl42LAJIoE3mfjpP9KrKgftX/CkY4Np5iwmE/nFBoIABriQv6UNcXH2Gafd9rbbAMizwFSnj/h51jHy+prAK4zffxcd2QShqSgAVxs2QSTwJvNM9ECCXgvHGcjWGUj56N+ODqnZyktJwQwEA5dmRuB+tBi1WmmXsgNiY5EArQx3GNfxceo4QgecAm8f9uzJtMs1BEG4OSXgacMmiATepE5M3ESvttafD/71j44NphnL3rAWD6y/nNv7D6a7pthuZatdXfEN8EMGss9CZYyaoVk7oYOKXbsaOMWbIAj1JiFa4LYSCbwJ3RvWGd9nw1EA6bmV5BxIdHRITWrAgAF1Oi97+ybaAqVOsLV8BJNGqO0aR2DXGABy80ARYaHvmQMQqeTkSSVGo9mu1xIEoXqiBW47kcCb2LEZP9C5alz4lgcmOTqcJnX77bfX6byckycJAsw+cOhQT8aPt+885YFVceSUgNLXRPfDR1BGWTAp2nLsWJ5dryUIQvXEO3DbiQTexGa370XnO90ByDidTUUrmui/tLS0TuflFZQSBEhdlGT96k2nTr52jSOoaka2fBMoLWYiki/iEq0Hr1ASE8V7cEFoCqIXuu1EAneAzU+/TbAKKmXYcNcER4fTZJYuXXrTc8oLCyk1ygQCZZN8cT1chlJp31/TwLg4wLq0qCUP2uSWEeCdCRpv9u3Ltuu1BEGongKurDbYkE0QCdwh7g+fQc9+1l7VWbt+Q67j8pqtQc6OHaixvuNKntCLGDf7dWC7zCM4GBdnZyqBolSo6OLEkNydEKFiz55yu19PEIQbKSTwcGn4JogE7hBKScGpv96HhwSFBti+YKGjQ2o2cjatww8oVcCvHsMZN9g+w8euJkkSAe3bW693CeSO0K+qI9vp0yrRkU0QmoBCAa7ODd8EkcAdZtTgZfSIsP6c+vl7Do2lOck9sJsAoNIF9p4fyPhx9n3/fVlA9+7W6xWCOshE96MpSGFglgJFRzZBaAKSAtTuDd8EkcAdJsDZBdWCOFRARqGJsxs2OjqkZiHvQjoBgBwocXZ3e2Ji7DOF6vUCqzqy5VWA0s1Eh+RUnDuWg1coSUlZjXJNQRCuogDcbdgEsZjJ1bRaLbNnz75m36BBg5g2bVqdyyiqQ6/yy+fox3xNV00MKTrY88y9aAacrfH71e2/ft/Vn2s7Zi/1LbN79+5otVostbzz328cy7bYeYz3/wTjbxaKHy60NcxquURFAVBgAalIxvdcIZrAArI9fdm9+wCTJ7ezy3Ua4763JK29/iDuQY31VwA2vMs2mUzkt4AFTS7XPz4+nvj4+OsPe9lStkjgV/Hx8WH58uXX7MvPz8fPz69e5dTlfD8/P8b6+XH4fl9S/l1A+vki1BUGvLy8avx+dfuv33f159qO2Ut9yhw5cmStx8vy80mNmMfjq+P4/pWnicgtxs+vi60hVstj4EDAOqWq6QIYOzjTr2g/69tOJinJYtd71Rj3vSVp7fUHcQ+qrf/lFngDqYyqFnNf/fz8mDVrFrNmzbpm/4oVK2z66048QnewQ0+tINwJTMCGabUnuJausLD21nTeju1EpC5h44wU2k7cz9BupkaLxdnTE08vT8xAQSqYuygYkLoXopScOCFhNouRAYLQqCTA2YZNEAnc0R7sOI64odapQnOTT2A2Gh0cUeP5/PPPaz2eu3k9A8rW8daxseQM1jBhTONO1xAQHQ1ATi6owkx0P5ECYRIGoy9nzoilRQWhUYl34DYTCbwZOPTXP6JRQIkJfn3qaUeH4zC5B/dYe6C7wv7jfbn9dvtOoXq9gF69AcgrBpWPic6HT6FubwSPEJKTRUc2QWhUl9+BN3QTRAJvDh7o/zd6dpQAKN38g4OjcZz8tEtXeqBnbA/E17dx51sKqJoTPd8AKkyEHMnGLawUPAM4eFAkcEFoVKIFbjORwJsBtUJB2St34ARklVjY8+/aHzXfqgp0FfgDdHfG5Yih0a8XUDWlaoEMXASDv4pupiPg68zevfafAU4QhKsoEO/AbSR6oTcTMVN/wOinITkfUt95mtuffNjRITWpypISSk0QAJybHEzMjyWNfk3/Ll2QJNDJYDgHhhg1t6f/xp72g0hJsCDLMpIkNXocgtAq2dgL/WZ0Oh3Lli0jMjKS1NRURo4cSa9evWo9V6PRABAZGXnTUTPNgUjgzUQ7Vy9yHggk+V85XEwro+hSNl4hQY4Oy65GjRpV47G8fftog7Vj6oHB/Rid1fgPh1QuLnj7+aPNyyM3DRJ7j+fM+nG4tJEoLdWQnl5MWJhNwzQFQaiJjePAb2bGjBksXbqUyMhIwPrvz+rVq68k6astW7aMBQsWXPm8cOFC+vTpU+25zYl4hN6MXHxmFe2cwAisv3uco8Oxu7iqR9bVydu8Hn+gRAHbuYOxo32aJKaAzp2s1y+AjeUPM/vdtrQzmMAl2C4d2Uzl5ViMRpBl2LMH3n8f1qwBvd7msgWhRWvEd+A6nY7U1NQryRusreqEhIRqz//uu++u+ezr60tqamr969TERAJvRqaFDyGmv/WhSFHiIWRZdnBE9pWTk1Pjsfzfd+ALVDrD78f6ERPTOHOgXy+gbz8AcktgSuRH/DzrGLk9zeDVluTkhi8tem7TJj5uG8TSsDDec3LiW687+cOTJfzwyx5Sl7zA152foO+YHLy75TNk6C42bjzbaP+9zQYDxWlpmDLSkUu1YK5A5tb63RJaoEYcB56YmHhD61mj0bBly5Zqz4+MjKR3796kpqaSmppKQUFBjY/bmxORwJuZ0y8/hJsEBZWw+U9vODocu1q5cmWNx/JSz+MHyL6QuSsAZ+emebsTcNttABSYYIJmLf++OJHgPmfB0529exs2TePBJUtYOW4cOdk5OAF3A2sj5jHumzg2lD5M2+ezWD/4D7z0gYVO0UZ2uQ5g/JQy7rrrR/R6+80DoC8oYP3s2bzl5sa/IiJ4p10Ym/180d/hiuV+BeW3OaPtoWHxpJfo+2A2ni9WorzHhFMHLRrvvYwb9yObNp21WzyCcA0bW+B5eXn06dPnyrZs2bIrRet0Onx8rn2K5+vri1Zb/fwOq1evJjIykqioKB577DHefvtt+9a1kYgE3swM6fU6ccHWn3O+/Idjg2lC2oIy/AApWk2b443fA/2ywNhYwDonuuUclMW5MiB7L4QpSUqq/7Ki57du5ef58wG4zRkWRIH77So0+Z/w8YwUvFOWIL9t5N5e/+Lbv1/gjwdeoPP00zA8hh/Xd2HcuG+prLR9BrrC1FQ+6dGD5BUrMJnNuAKVwNFKKN0F8g/gEm3AaWElR5RjeOllmX4X0nB5qhzjvT4Ude7LpoSBjBu3kzlz1tolJkG4hhKbEri/vz+JiYlXtrlz515TfE3JujpJSUmMGjWKLVu2kJqaSu/evW2uXlMQCbyZcVOo8Hq8AxKQllVJ1pGTjg6p0VnMZnQVMn5A+Qhfeno03fthnw4dUCoVFAPlqSB3leiXuh/aK8nPdycvr6zOZVWWlPDjzJnIskw/NYwcAbJBwZLNT7B701P0PTcLr9J17NgLE49+z3/+NJ7AF1JZ+cUslNMt0KMDO3dGMn/+LzbVqVyrZcWIEegyMgiU4OFAeGEe3Pc3BbMVsNx9CtOjN7EhaQptMiqYMnkJqxan8sRvr/Kg9lNcRpXBcBX09wHFVL76Ko977vleTC8r2JdEo03kotFo0Ol01+wrKCi4oVUOkJqaynfffcfcuXMZOXIk586dIzIykoULF9pUvaYgEngzs2iRG/q562nvChYgYc54R4fU6HSnT6MA3ICj43owdmDTDY5QqFT4BVsfeeSmgyrcSOyxoxAigap+78F3vf46JQUFBEpg7DGFeXmbeOexBayRZnJi3yBWBq4C4IAJ0n8Gv806vOfk42su4Lncd2GMGkJ78dlnJaxefaxB9ZFlmfWPPELhhQsESHB3MIQ+CJUPOuG/1Z2T98cR32E+s1fF8Y16HsZtcJfTKha/cRceqy7x1iOv0N9/F4rbjdBVBV3bIEmTWLs2nVdf3d6gmAShWja2wGvTp0+fG1rgOp2u2pEwSUlJ9O3b95p9n3zyyQ1/ADRHIoE3M4sXuzE8oAtd77D+iZmfch6z6dZ+fJn/azy+QAmwzesORo1smg5slwV0iwEgRwtO3kaiD5/FKbISvILrvDZ4aXY2v7/3HgAjfOAnzTzGrohjQ+48TiXGYvlZ4mJaMAVhdyAD2wvAuA6ijqVzalkEf3rvLQJuy4T+zqAewbx5CRQWlte7LsdXr+bk2rU4AZPagNd4MI9VoNviRbHRkz9/tghtXAkfz0hBc3YJJftA3gfBp3MJ7JVN8vOx3PWHeLqsL8fF3wwdnJB9PIAR/OMfu9i3L6PeMQlCtRqxE5tGo6FPnz7X9CRPTEy8Mrb7cmc1sK6SeH3ntsTERGbMmGFjBRufSODNVMpfXsJTgmITrH1snqPDsYvx46t/mpC3cyN+gF4Nu1MH0rGjd5PGdXlK1Tw9qOVKvHKK8ffPAS8N+/fn1qmMvYsXYzKZiFJAxBCYevpjlv8jk0uDJMoSPOH3ctTK7cxf+wnOrq6kWeDsfnDZZaBj0Dm+7v8AvZdqcQmSoZMPBQWxvPLKtnrVw6jXE//sswAMVkHgQFAMhyw3f3zeK+aNlQs5m92J5Dum0v7iSwSVrSO5BAzJoDwqE5yeR8ULKhItk3n9+TLCMowQA4ruLiBFIcvtefzxn8SjdME+Gnkq1dWrV7NmzRrWrFnDokWL+OSTT670TF+6dOmVjmoajYbHHnuMRYsWsWzZMpYtW4ZOp2sRE7mIBN5Mze7/Z3pEWmcBK177pWODsZMuXapf2zv/2GH8AHMbuLArDKWyaX8tg6oen+WbQXkKimPa0Fe7H4KU7NtXedPvG8rKOLhkCWB9baz0h4EhCchfnSFH4wu7TJBzlGeeiaRTzygGvPgiAHvLwLgZgo/msjV4Fo++60t0USn0UoNbLz766DSnTtW9J/y+996jJCuLAAl6hILqdiiKdWPvC8O5d8huth0dT86RUOTdCn51fx2AgyYoOwqWI+CbriPUnEnUhJ2sfeg05ePykHqAxUcJoWrU6qEcPpzLihUp9b3FgnCjRl7MRKPRsGDBAqZPn86CBQuuGRb29ttvs3Tp0iufe/XqxYIFC5g7dy5z585l+vTpdqliYxMJvJlSSgqk53ugBC5qTaTt3O/okGyWnp5e7X5tlhY/gBAFrkebfnxyYNUEM/kymM+BKUbFgIvWtcGzstzQamt/lH3k66+pLC+nrQThfcB0RMme+f04mhdHxS43OF6Ji8tBFiwYCED/Z57B2d2NDAucTwanAyYmjVzGqufOM7riC+gDdHTDYunPX/5St1Z4hU7HnkWLAGvr26Uf0Ad0JzzYaH6UB94MQf1DEOYjKrhg5lzmIC4RTjlwTA+mkyCdBv9MLX1f3sYblgeYoViJk68euoKqqzNGkzfQkb/+dQdGY/176AvCNcRiJjYTCbwZ877vJzq2sf584s8POTYYO1i1alW1+7UlRmsLvKc7PVyafoayNm3b4urqSiWgOwuKjiZ6nDgE7SRQtuXgwcxav59U9e67uwuoe4PiJPz3zrvJOR0C+82Qe5LHHovE39/6r46Llxd9nnwKgMRSMG6FmcH/ZfbSP/HK1r8Q2jEN4tTg1o3Vqy9x+PDNO9L9/uGHVBQVEaqAyChQ9YWCKC80b5ThMeEEr33owoUwFziIdSwZx0lxt746OGyGyuMgnwLfLB3B5iyS/i+WR977Cp+oXBSxYHJVQFsVzs63ceGCjm++Odrg+y0IQKO+A28tRAJvxrprgome5AFA5qlcKvX179TU3BnKyigzgQ+QOaU9YwcomzwGSZII6NABgOwccPGvpMvhU6jCTeAZSmJizQk899gxMo8fxwno2hmMl5Scubc9Ry3dKdvjCQfkcFsAACAASURBVCcNSFISzz572zXf6//00yiUSs5ZIP8wuCQb6OBxntOTo3hm+3vQQ4IOLkBf3nhjZ63xG8rKrnSgu00JTj3AEgtlia6UObmT+KfOXJgDFbISMkyQVwHsYcYrj+MeGEiBDJeKwHwelOkyPnlFmKdAYGYugzJ/Q2pjhChQd3KmstIHaMc//7n3lpspUGhiSqCNDZsgEnhzd2jhP/FTQrkFfrznLkeHY3f5+/fjBZQDe2MHMWpk08yBfr22/axTquYUgbOTkaDTOXiG60DjzZ49NU8Be/hz69KvnZTg1g8UCbDukfGkX4iAI0BGJpMnexIRce20jh7BwXSdPh0ZSM4H824IvJhH9gI/Hv7iKzRd86CLGpy78f336Zw4kVdjDMmff065VkuQBBHhoOoF2ggvvP9WwnevTyO7JJiSCxprPCcrQT5EaKiS+U/fTvcHHwTgmBmMZ4Cz4JujI0yVwYEne/LQv7/Eq50WZTcwtlGClwJn594cPpzD7t0XbbjjQmsnS2BybvgmiATe7D3U/VF6xFj/MxVt2+TgaOyvIGE9fkCpAn7VjWjyHuiXBQ0ZAkCuAVRpJsqi3OlRllTVka36meFki4UjX30FgNZ/CvMyNrHeczq/dh9B0WEfOGQEcwpPPNGn2u/3r+oxftQMZYfAK1lPaHgGGf2DmXdwCcRJ0N4V6MmiRXuqLcNiNrPv3XcB6KsC5ziQu0LFHmeKNZ6sGTuN/NNBcFSCSyakwgrgAH/5yxBcXFT0qErgZ8xQfgHki+BWUIlHWQnZc/wZsuE3wrzTMPvK4ANShBMGQwTgzscfH2zIrRYEACxKCb27U4M3QSTwFiH/T2NxBrJLZQ5++qWjw7GrvN+344d1EZMjv8c0eQ/0y4J69gQgVwbLCajo6cQd57ZBtJL8fE8yM29cnzxt1y5KCgrwBPZFz2Psu3F82+VZTuV0xpyihot62rW7xMiRkTd8FyCkf3+CYmKoAE6kgiUJQjJzOL+wHfM/XYJzrB46OoGyOytXniQjo/iGMk6uXYsuLQ0vCaKDQNUTdOFt0LxZwqo37iS7LJiSNA0cA05XIsuHCA93Zs6cHgD4d+1KUI8eVAJpBjBdBNLAJ6+YAL9czt8RwYRNP+Gi0UM0KMPUyJIEdOH77483aKy6IADISBiUTg3eBJHAW4TeE78jxt/68/E3nnRsMDaYOnXqDfu0qan4ARZvMCW6Nn1QVfy7dEGlVFIsQ9kZkHqauT1pH4QpQB1a7QQmR5cvB6CjEu4KXsLGe49QeV8OhWmBcMgMRSeZNasDCoVU7TUlSaLvM88AkGIC4wHwOabDu18Bm3yn0P+LHFy8ZQhzx2SK4b339t1Qxr5//QuAnkpwiQG6wMbPpzK/7S98Y36IgrOBcFKCTBNSYSVwkJdfHoyT0//6GnSdOROA02YwXgDSQJNfTKCcw7e9HiD1yzF4p+lRdQKTJEGgCnf37lRWmvn2W9GZTWgYCwr0uDV4E0QCbxFCXdsQPCcIgLSLZRRn1N4rurmKioq6YZ82r5Rj7lP4a+RmwnMc9+uoUKkIaB8BQGYGOIcZ6L7nCOoOBtCEsnfvtUPgLCYTx1evBqBrAIxxW8cr0XPIm+hBaaInnDECx7jnnk61Xjfm3ntxdncnS4aMY6A4CpG56fzU6VGeedOF0EojdHACqSdLlyaj01Vc+W5mYiLpv/2GMxDjCao4KA13Zutvs5myrBOVG6Iovuhtffd91oAsp9CunfpK6/uyrlVjXs9boDId5Exw1hlpU6rncO44pn8SSehOEyYnGdqCU6SasjJPwJ+VK4/YcNeF1syawF0bvAkigbcYybM/I9wJzMBPM8c6OpwGOXfu3A37dHoLOyLnMfXTGLxdHdu1NPi2AQBkFYCruZzdxhF03VyGi8aDrVuvnVDl/LZtlJeU4C1BcGew/K5g16O3kXqpA5wEMrSMGOFESIhHrdd0cnenxyOPAHCoDEy/gf/ZAno9sI61D58m0v836KaAthpKSyP5+OPEK9/9/f33AYhRQpsYkGKgfLsrQ9Xf8PZHMtk9ZOSTCsgyQ14lkMjChQOvaX0D+EZH49+tG5VAhhFMmUAGaApKGDBpNauevUAv9QZcfcogCkw+KlCBStWFPXvSuXCh+c8ZLTQ/MhIGnBu8CSKB292iRY3zaGdKUD9iB6kByE08gmxpedNZrl279prP5YWFmCxwZ+oSVv4jnSljHFun0KqpE7MNoDpp4YfOT/DKK5WEhsDhw0rKy/+3VvexFSsA6+NzdRRYihWsHTGZ4lRfSDKB4SgPPdS9Ttft88QTAJw0Q9ERcD5i4r6+XzL+/97nnzvnQ6wFopyBPvzrX/uoqDBRfOkSR7/5Bgnr+HOnOKgIV+L5Xhn5/3Km6J1cCoK8IAU4WwnyMYKCJB55pFe1MXSaMgWAVEvVe/B08NKWMGb0OsY/9z5PJf8Lz2Ad6kiwSBIEqXF2ts6st2pVwxZeEVo38Qjddk237FMLoNVqmT179jX7Bg0axLRp0+pcxuLFfixYUPv0l0VFRbUc9SM//8bvFxUVceilR/HctoQiI8Q/+zTRzzxXY7nXX6P2azZMfcssKytDq9ViqfrjIzf+F3yAAWXr2LqwH1M14dXWvam4dbI+7s60gPkIjOrwFSueeYqMwADMlmA2bz7OwIEhmI1GTnz/PQCdA4EMOP1wJKmVHSg/4AYXy3FzO8PgwbF1u0e+voQNHcrFHTs4lA137IGA2wtQTDHi+loFMwu/ZVW3++BEALm5gXzwwS66ZPyIxWwmWgH+nUGKhfJtbhT18WJr3+FoTwdgOa2ETDNkGoBE5s/vQWmpjtLSG0MIHDQIsCZwYzo4Z4FLiQHX8nLKhzkRfiqdTpZT7FSHQJCEU7iKsgw3IIBvvknh4Yerf1XQGL93LU1rvwc11d/aAm94ZzSTqdKh/17U1eX6x8fHEx8ff/1hL1vKFgn8Kj4+Piyv6ph0WX5+Pn5+fvUqpy7n13ZOTcdmDvuAA+FL2HUBMr9eRp9XXrvh3Ks/13bMXupTpru7+zXr8Wbs344voFfCjrND+OfkQLvHVx++vr64uLqiLy9HexLG/H0ddyyIJ+yJDNgUyaFDZ5gypTtnN22iQq/HR4LgTsA+iS1LhnPxYnvr4/PsdO6ZE0pYWNs6//4Mfukl/rtjB4dM0O8geB4qp1v0CXa92p9XXn2DNfNmYjnsDPn9+fj9Nfyh+AsAequtrW9jiESbP5fzxa/3k2qKRJvqD4eBMwbgBEFBZl54YQiurupqr+8zejS/+PpSXFCAthzc84Bs8NSWEhiSy8mJHRkZv4Xfuw2gMswNS7YKlKBWdOTQod2UlakID9dUW3Zj/N61NK39HlRX/8st8IZSqcwt5r76+fkxa9YsZs2adc3+FStW2PTXnXiE3oJ4qFSon++GEkgvMJKflOTokGySn/wbPoDRBTL2hDo6HCRJIrS79bH3xUxwV+vxLC4hxjcF/HzZuNHaefDol18C0FEFai8o7eXGzojBFJ30hWQjWI7x0EM963XtqNGj8e/UiTLgSCpY9kLwxRyUU43IXhJPnv8QYlTgE0JoZjKG0lJCJQjrCIreYPjWmTP3tGdn7BDy0wIxn1LDRTNcqgR+5+WXB9WYvAEUSiVRo0cDkGYB0yUgE7wKS2lLDiendWDsD5tpE1iMcwcwyRL4qXB17QzAunWn6lVfQbCgwIBTgzdBJPAWp/iutXSq6uuV+Mzs2k9u5grSMvABZB+JoPTmMS1nxMSJAFwsBpcTBnLG+jLpzDrooOLAARWFeUWcXLcOgK5tQT4JSU/Gcqa4I4ZkF8ioIDIyl4ED29XrupIkMfiVVwDYbwL9HvBOLCVOf4TN/xnKoJUH6ZZXgmeYiUFYH8P1U4NzD/j5xFSePbmORXf8ibOmDhScCYQk4FQlyEeIiJCYO7f3TWOIGjMG+F8ClzPBvUiPi7mcgtEauh04SZhbOpWugBeo26koLvYAvEQCF+pNvAO3nUjgLcyE4A50vtOawdNP5VDRgt6vzawab3yZrqAcX0DuoGZYO5NjgrpO+7HWHv7pZjDtB+6QmbJrA8QosUjtWffOV1RWVBAgQVAEyIUKfpk0hsy0MEixgPYEc+fGIEnVj/2uTbe778Y3MpJiGZJOg3kbRJy9RGhsBv/t+zyvv1VKZ7dLOFNCiARRHUHRCdbum8e4lbGc2j6CnHMhmI+qIdUMWXpgH2+9NQJn55u/LYsaNQqADAsYcgEtKIqhTVEZvm5azg1uz9Dft6F0MkIYKNtay5SkSHbsuHDNEDdBuBkZiUqcGrwJ1bwDlyTpj4BvY15UluWXGrP8W13C04sJ+WYel0zwy/RxTNtS/TSbzU27dte2SnXlMj7A2YH+TBjcPP6ibtuzJ65urhTry8k9BB4LSwg6VED4pPOk+UVzbvmjqIDOTqAqh7RHgkhR9KD4oA8cMaBUHmPOnJk3vU51FEolI955h1XTprHHCJ1/Bf9eenqHHCJhXiLfzY3hD4cXc8B9Cmei5xHiu4QxS9cxcNJ3fPhRe/JGKSk4HgYHZDhaAfJehg0LYObMbnW6vkdwMP5du5J3/DhZFmiTDeos8NCVEeSTw/lxYQzduJPPn3yUivY+VBxRgKcCjbILhYXJxMef5e67YxpUd6H1saCgXLSkbVJdC/wx4ACQ2Ejb3EarTSvxVK+59K5atSttx14spubRer2ZEydOXPlZX1CAxQJOwLHRfbjttubRGUVSKIgaNgyA0+mgOVHM+RlhPHruE9q0LUKVfQQJiAuBTRem8FThz5xY3RtzshLSc5k+3YvAwIaPZ+88dSpRw4dTCfxyCSq+hYhdl3gw7n1iLk0itmwtlsh5TPs2jh+188j7ow8F/3bG9P4Fcvx9kfep4JgBdBm4uh5l2bKJ9Xoa0H7ECMD6+tycBWSDh06PD1pSx0cwZONvePjpMPlj/fM/UEVJifXDzz+faXC9hdZHtMBtV91zNUmW5R8a64KSJL3dWGW3FiqFguTXZqMZ+QU6I+x78Rluf+c/jg7rpn755Re6dLGOHdZu34ovUAxscRvNLM/mMzFDl4ce5ugvGzldAUO2WTA+qmDelGXsCzUDMlEK8AZW9XuKh54O5OV/aOCECSqTePbZ/jZdW5IkJn35JUtjupFWXMJ3+6DLczLJRcex6OCABGOLlrDpuXkMf3I5CQ8PYw8DOXO+M6W/e8E+M5woAjaydOl4oqPr9zCt/YgR7P/wQ9ItYMoCcsC1rBJ1pQFVpJFyjQsDMg+wRmoPbSU89GpKzhiAdmzadBaLRa5x6lhBuJpogduuuhZ4QiNfs7HLt5vmvN7x6Ns+one09eeUpR87NpgGyN+2Fh+sQ8j2HbYt6dlb9IQJODs7kSPDpe0Q7pzB2eFhjDlzlt9iNyF5T0Eqh5C5e3njQyfSFE5wupgBA4q57Tbbe9N7tWvHfZu34OLuRpoFNqVBgQ58JLirLUx7YB1vvTMD9UOF/MI4fr80gJyEUPjZDPt1YFnLW2/14YEH6jaRzNXChwwBSSJbhspCsBQBheBRpCeQHL7pOYeyt7vheaYCQqHMSQEqaNOmI3l5epKSsmyuv9A6WFBQiXODN6GaBC7L8uONecHGLt+eTBaZonIjBlPzm/Us1M2Fir/0wxnIKbWQ+u1KR4dUL/mH9+ELGF2hMql5PD6/TO3qSvd77wPgtwug+aGYE975JEY8xuOr49jebh6Z3/iRM9mTs3erqDhngaJ9/O1vg+0WQ2j//jxx+gwDH36YLmGB3NHOmQenqPF8x5Uj8zvzeacH+I/+SX7ZN5lLS8Lho0rYdpwOob+wbt1QFi4c1KDrunp707ZnTyxUTWiTjfUxelEZgeSxQ303j77tS8gJM67twYIEviqUyggANm06a7d7INzaZCTKcW3wJlTfic0TiJRl+VA1xyIArSzLN65reIuyyDKFegOuTko8nFUN6l3cWO776geWB4Tyey7sfG4uE49cdHRIdaZNy6Ij1iFkscbm9w7/9tdfJ2nFcs6YLfz4T7ioz0DlvoS1981jxsSP2TL4DnYXDKH0Jw9IzmHSJD3Dh7e3awwewcGM/OyzG/bHVW1/dAduA1MfmcqXVJhM0Xh5xdp83fBhw8hKSuKSBTpkgzoHPIrK8aCYwLsT+eXecOSnIyh3iwY3cApRUZTkBngSH3+OP/95iM0xCLc+WydyEa5rgUuStBnQAUmSJJklSfrHdecXAndLklTQVAE2F+UGM/mlBiqMZkeHckVeQghtnwpGAVzMLqf4eMuZk1qnrcAXMHdwZlL/5jchoFe7dox5exEAx/QgA/eY1/H5xLH4LrzI16ZZnF3fFdaU4++awLJl4xwWq0ol4e6uwsvLxS7lRQwdCliHk5mzgTxwKjfgVGGgz4TdPB3wHKPbbABJhmBrAgeQpHD27cuguLjSLnEItzbrI3TRic0WVxK4JElvYZ3+wVuWZQUQDfhIknRl8lZZlotkWf4Eax+eVsciWx+pF+mNWCzN4/342vFf0bmNNcH8/vi9jg6nVvfff/+Vn4uqhpCVDAxg9Ijm+evU54UXuH/VKvrHhDOttxN9PvDix+cn8mLZW2z/cDSmN7T00mxn794JBAU5diU1ewobNOjKe3CDDizFQCG0KdLjTx5po8IZ8OseXDV6VGFQigJcJHx8OmEyWdi+/YKjqyC0AJc7sTV0E65tgWtkWX5RluUiAFmWU6veV78kSdJ3VY/WL2se2ctBKkxm8ssqKTc4vjX+aq+RxNxlfR907kQWJWlpDo6oZoGB1rnOK3Q6zBbr+5uU4f0JC6t9yU1Hipoxg7FHLhCbWIn3H3TM8NrAzsBRGF5wRU4N4uDBCURF+dy8oBbE1ceHwNhYzEC2DOYcIAfaFOsJJJ/zo0IZuGUfbn4lqC732fNTode3BWDLlhuXjRWE64lhZLa7OoFX+3+dLMtJsizfDbx8XRJv1WQZiiuMFJYZMJkd28ntv/PfoIOzda3whLtGOTSW2qSkpABQsGsHvkARsFkxpln1KxCswgZbO+RdslyVwIvKcaeEolgP3EvLiNEfp8IJ8AS3cBXl5SrAny1bUh0ZutBCiBa47a5O4LraTpRl+UVgriRJ9u2p08IZzBa0ZQbKKh3XEWtB92fpe6d1WMXxg2coychwWCy12bJlCwAFW7+/MoQs+dDN5+gWmt7VCdyUA+RJOFcYUFca8ZfyOTGyE8MObENSWCAY8LO+B3d2juTUqQIyMlpNP1ehgawtcDGMzBZXJ/BVkiT9QZKk/6upk5osy+8AkYBoMl1FBkorTRSUOqbzjqeTkq+ff4koZzAB22aOcUgcdZV/aI91CJkLuJ5qnu+/W7uwqvXBsyxgKgC5VAbd/96DZ4xqx4Bf9+Ki0eMcBnqLBO4KPD07AJCQIFrhQu2svdBdG7wJVyXwqzqorQH61PQFWZZ/BaKaILYWx1TVsU1vaPrW+GOd/0S/SdblIo/sO44+O7vJY6irggtZV1YhG+rTqrtTNFueISFo2rfHAOTLYM4Fcq3vwQMo4MLIEPpvPYCHd5G1BQ7gp6SoyA9QiAQu3JR1OVHnBm9C9RO5nJdl+XxtX7rZ8daupMJEkd7YpDO5tfVQsXrBs7R3ApMMhx67p8muXV9FBeX4AOYoJ+4caZ+hT4L9hV/9GP1KAi+nDcVUhjhRFOTFwMx9VKoADXi0V2EwKIAgtm4936xnMhQcz4IkWuA2qnE5UUmS3rSl4KuHn7VGFSYzBU3cwW2Y5/PcNt76LvLo3mOUNtNW+EEm81rsJhLazaJf30Zd+E6wQbuqx+iZlzuy5Um46CtRGUz4kc+JUR0Zsn+79T14WzB5WX/33N2jyMoq5eTJfAdGLzR3smiB26y29cAbvGpYVUe3Gh/DtxZmi4y2rOkmf+nkq2LDq48T4QRGGbZPH90k162rhx9+GINez/mIeYxfHcdmeT5ubmpHhyXUIGzgQKCqBZ4HcokMpdbH6P7kcWl0OwYk7MNFo8clHMotEngocHe3vmETj9GF2lyeia2hm1B7AveWJOmF+hYoSdI0rMuGahoc1S1EBorKjRRXNM0j9Z4ebzF4orUllPzbEfKSkxv9mnXl7e2Nds9uuqQuYcOMFCp6tLoJ/VoUv86dcfXxoQwoNoG5AMgF96r34OlDA4lJPIafay5y0OUvKdFqNYCKX38Vb9qEmslIGHBq8CbUnsABRkuSNLwuBUmS5ClJ0nfAalrpTG21KTeYKdQbMTfyDG7Do9z57o0niXUHCxA/rfm0whMTE9Fu+57bytbx12NjyZSCbv4lwWEkhYJ2l1vhlyd0yb08L7oO3GFF10eIeEeJdMIC3uAVpcZkkoBgtm+/gKkZLgQkNA8WiwK93q3Bm1B7Ah8ly/IYIOpmSbyq1X0emAEsrJqK9Xv7hXlrMJotFJRV8uZb1XfA2LjZQodYAxs32/aPXi/F3+n7hCtq4NyFfM59941N5dnLjh07yDu4E2/A6AR+ac13BjbBKuz69+C54KKvQGUy40c+m9rO4dnXnAg9YYS2UOmuBECjiaaoqJLDh3MdGL3QnFnMCvQlrg3ehFoSeNVwMaqGllWbxKta3fFYW93ngaiqseLIsjyzcUJu2WQZ/vmue7UTv/z4s4m3Vxfw48+2DUMb09mN1+99h/4B1s/xf5iDxfy/9/CvvWZT8fV2+Xrx8aBNzUADWLwkxrYX0wk0d5cT+JX34MUg6a2P0f3JI2TMQdbOOU1eXzOuEVBhkcBTgVodAcDOnc1zUiGhGZAlLJXODd6Emz9CB6pP4pIkPYp1dbJRwIuyLPcRw8vqrrTShE5vuGZRlDsnqFg4w5c7J9i+Otds9V2E/yMADwnySgycfPP1K8f++lebi6+Xy9fbsgWK8vTWBN5OzZTRXk0biFBvwX36oHJxQStDeQVYCoE864QuARTg/Vg6/zw/jQ5dUzBbp7pHEagiP98DcBYJXKiZRYJSVcM3oW4JHK5J4o9KknQAWAb8irXVvbixAryVVZosaPUGjFVDzcaNVnD2iBPjRtf5P0uNugQpeS7yc4Z2tH7e8+FHlFy6ZHO5tioqteANlMd506mTmFq/uVM6ORHSvz9w1WP0HOt4cC8KUSpMJI3tzqTkjRiUgC94RamQZQkIZf/+LPR6oyOrIDRXFqDChk2oewKHK0lcAfQG5sqyPFq0um1jtsgUNtJQs897jcX8SQwRKqi0wM+jBtv9GvVVbLIOTzg1vCdKpe1/qAiN74Z50XPArbQCpck6HvzU+GiGbd96ZV70MmclSODn1wWDwcLu3RcdWwGhebIApTZsQv0SOIAsy8uAxwAxyNNOrh5qZk/+XkoW5KxkzCwlauDUifMc+/cHdr1G/TyGG2AEfgma7MA4hPoIHzIEgIyqFrisA6lCpk1JOQHkkjcqiD67kvB2z8c9EgwWCbyVWCzWtUbF8qJCtUQL3Ga1zcS2pKZjVS1x79p6p9f2faF6l9cXt+fsbVumxPHp355haNV81b889yyuOGr8dRs0QJkEh47e7qAYhPpqN2AACpWKXBkq9GApomp50TICKcDsLXO6RwdGnd2MwRdQgFOICq3WFfBg82bxt75QjUZuget0OhYtWsSaNWtYtGgRSUlJtZ6flJR05fxly5Y1tFZNqrYWeK0zqcmy/D2gqSWJt/qZ2BpKW2a4ksxtpVZLlO99mcClIYQoQW+SWage4KB5qvfgDVSoISxTzPPTUji1aUNw377IVL0HzwZyrOPBPdHiRCWHx8cy8WA8RgXgD+4R1k5GTk6RpKTkkJVV4sgqCM1RI7fAZ8yYwfTp05k+fToLFixg4cKF6HTVr5qdlJTEm2++yYIFC5g+fTpLly69acJvDmrryte7pmVFr6ORJKm6uyL+hW4gGSiuMGIwWfB0VSFJtg23+vsMX26LX84Pd47gqzVgMJ7ht2eeZNAH/7FPwHW2Fw1gbgOTu4pepC1JxLBhZOzdS7oForPAKRtcS63jwf1VeZybGMnjkz5FMdmI3E6NLkcBatBoOpGbe5j4+HPMmdPjSnlFFy+Ssnw52u3bcC8uoKu/goBuChjgSnmc7/+zd97xUVXpH37O9GRSJo2EQCgJHRJ6L9KrFJG2WFZZV3Rde8H96dp3FcvuWlYFG6KiUkQUMAiCiMAqofeWhFTS6yRT7/39cScQIAklndzHz3zCzD1z7jnXO/POec/7fl8yw1qSneZHdqo/aUFdON18ENp0C+50LV4ZEtp8N35yKVrtWUJCnPTq1Zx27QLr8QqpXDVuam0vOz8/n/j4eCIjI8+9FhkZyaZNm5gxY8Yl7f/85z+zYsWKc89/+uknLJaGb8Iu9y2aAOReQ79BQI/LtlKpEpvLjdMq4Weqvl74b1MGszbxdibtXcrq07D57XdpMWY8bSdProGRXp6yIL0AQArTcaOaQtaoaDtyJL++/DLJnhW4XACiRManwEqzoEySurXFbjJwQ8bPxLUbQ1GcgGAdhXmhgIYffjjFHXf0QJZlfnvrLTY98QRuh+Nc/zuAfhtgSAswR0PUUThum8aqqPu4MelDxvnfyzP/+jvfhU0nLyMYDunhkID94WA9BSxh7NhmvPvuRKKiVEPeKJCptb3suLi4SwywxWJh48aNlxjw/Px89uzZQ2RkJHv27MFisVxg+BsyVbnQN3lyu8dew6M3SoqZSjVxSzJ5JQ5KHO7qub0dRpZ98TLiq3b09VE+Oytvvom8hNpNInB4pDTLUuUsQElXC2Fh5lo9r0rNEjF4MFqjkUwZSmwg5QDp4JtvJZxMtMLB1llDmb1jOcVmwAR+7fXYbDqgORs2nMLpdLP5qafY8NBDuB0O2mtgtC9EewSHvjZNZb4llg1bp+JuDd8OvYcbP+jKiin3ceTZDky542f6vphNc7JgggT9dDDUGwLbA7fw44/F9O69mG3bztTfRTW5NAAAIABJREFUhVK5cspW4LWwB56fn09g4IU/5IKCgsjNvXQ9GhcXR2RkJCtXriQyMpL4+Hjmz59franVFVUZ8I3V7Lu671cph91TntTuuva98Z8fC2d26gr6v+5FhBZKnG4+6xFNcUZGDY70PMV2F3kljgteCwBODouu9raASt2i9/I6V50sSQJXGpAOfnmleMlF+FDMvlnRTIldC7IbWoIrUHHwBQZ2paDAzvIFr/Dryy+jASbq4abOMHAeTP8IJizxoTTqXm5eHsOnkfdii4ebnO8R+8IBbu71HsHTMvh44pPMf8lC53X5BHXMRAyXobUG3SBv8DGj18+goMCLiROXsW9fwyylq1KOshX4NT6ysrLo06fPucfFgWcVGeuKKHO3jx49GovFwujRo4mPj2flypU1Ms3apCop1WqJs6jiLjWPW5LJL3FSUOq8QMHtaljdrRP3T3qb6X+CEAF5hVY+794NW0FBjY3T5ZbItToqlIv1Ab5tM7PGzqVSd0SNGwdAogSuVCBdYCy1Yyx10Jx0pK56Npgn0/vzbCyyRIlbgEUDRNKcNE6++RwAo/XQNRJMo0AzHtJjggnJtjAtfwmLZx7AHP8eOekwMmMN74dPYEbbLwnJymPEtKV8+Ho29/3yLEMcW/EKLsYwFlzeAvNwb5ySAX//GRQXu7nppq/Jz1dzjRo01VyBh4SEEBcXd+5x993nK2BbLJZLAtZycnIuWZWDsjdusVgucLlHRkaycWPDX4OqShqNEJvTTbbVDnDVbvWOUSa+umEKH786nz/cDBYBGRnZLO3a+YKV+LXopcuyTLHdRa71vLrchdyPU8Cho8OuvnOVeqfd+PGAx4BneuqD54JfXjHhnMVP5LG88wP87R8SzTJdoAFjax2FuYLprEFILrprIToQvIaAGAZ5bXyx7vUi6F95xP3ajebd36JV8Rq2uKB0n4CzMhyCsFPZjBn9PX9c9BTuV2z8a+4COkSfxmmQ8B0NVqHBMsBEQYEvoaFjSEzM54EHfqjnK6ZSJRJQVI1HFfTp0+eSFXh+fj5jxoy5pG1kZGSF0emNIYitRg24EOJxIcRjQgg1gK2WKbPb2cVXn3Imx4fwzovPsO6D6cwZDX4C0lPT+bBDFNnHjgGX6qVfzqDbnG6yi5VVd/mfFJs3y/SKTuQff8+lXTQYOgi6l/hc1XhVGgbNoqPxa9mSEhky3J5VeAr45xYTTDZ6HDjnnmXtH46Q09YG4WBorWcorxNCBgEChunBayiIHuBoZeCbbbN55YWPeOr/XmdTxCi2PvASNm8/EiVIccjY94D7sMBodRJ2Jpv2Ip7Tt7XC1sbMbe9+RPP2mRSFgqYVFAbpERYNWVmdMBhC+eyzA2zapOagN1gkwF6NRxVYLBb69OlDfPz5//9xcXGMHj0agPj4+HPHLBYLM2bMuCBtLC4ujtmzZ9fAJGuXqzbgQoj3hBAvCyFGXHxMluXXPNXIojzFTuoVIcRoIUQvIcQMIcQT9T2emmTLFmgX7eCnzTKFNifZxXbsLumKV+SnXg3jobffIu6LkdwxEUIFFBRa+TC6GwcWLwYu7KeyAiiK4bYrbv2Lzl2UdIYvX9/LUyuM/LxTsHBFDof7d2LuULUUYGNECEEHT9ZCvBucSUCyoouudyhu9LYzjzO98xvcmfU+hg5gLDzJEP4DwBg9eHcDbaSAGEhuFsKOj2czcVk3dp6YisPqR2JGd3YOXQDAry6B/TSQKyMdgJCMAox2G93FCda/PZI2b2bT8l9gPuogYCRICFqM8UKSBBERyjgffDBWrUneUKllIZcVK1awcuXKc0IuH3zwwblV9aJFi1i4cOG5th988AFff/01ixcvZsGCBSxcuJBevXrV2FRri2tZgQcDC4BNQgi3EGKXEOKf5QVdykReamqQ14IQwgIskmV5jyzLK4H5QojGkRtwBcRuhIUrcoj1bNO4JZkSh4usYjvFdheX2yLXaDScuC+EO779hF8/G82t86CdTgmWWz1/Pn8VQyuNUJc85yooVfbj3RedzJabw28P/ZUvY7oQtO9VFs08QJ8zH/HkTCvfOacydLCa5tNY6XTTTQCclMCZDHIOiAIZ/5wiIkghmEy+eHgWj3z7b0xhJUw4/ABanHTTQgs/MPYA+snkNvMl6MEieoV/z7/ftZM/yUri3laYSjTs1D9AiS6IDFkmVQLbARDxAmGVCE3JJ4wzSJGCL/s8zIJHZCLj7eTowLsLpNi0mNtoOX3al7Cwzhw5ksWnn+6r34umUjG1LORisVjOCbM88cQTFxjkhQsXsmjRogvaLly4kLvvvpuFCxeeW6k3dK5FC30mSjDxbOBDz7+fBDZ6DPpJT43wvjU60qtEluV8WZajLnrtuvGnjR8DC2YGMf6iLR1ZBqvdRUGpg1yrgxKH6xIDW0ZgoIHYPhZuW/EZy9+czc3PCcYFgB4IlrfzTlQU38+YQfrevYCMwyVRUKqs9otsrgtW3LIskxH3O7/c/gc+b9eavZ98hCRLjLOv4Qu/8bz0wcucPryK+A3BGI2qiEtjpc3w4XgFBZErQ7bD40ZPhMDsQsLIwIiNgj5+HOrfjacX30tU1kYMwBAdePUH0Q5c4TqcX+lxZuhIWxZA0ZtnKe1ooDTXh9KfQKMz85usBCTFucGV6NlvPwpBmfnoHC66EY/fvMOsnXsEy8DDaA0uvPujFFEZpHh4fH2HA/DCC7/gqCFlQ5UaRC1mUm2u6ZtUluUCYKXngRDCH6Uu+BhgJkqhk7sr7aCOEULcDSy8bMOL3/fpEnzzCxHTbkYODqmFkV07I0bAqYMGRoyo/Keo0y3hdEsU4QJMFJQ6MWg1gBa3JPPiC4LnnvNjs4Ch/3mH5Mcj+OugRUT+XxE/x8FRl8yeVavYs2oVfyeYNTOGEDBgMJt+i2DaLA1FWZm409PJ3bmdtLhdWK3Wc+duo4EhkRDxZ8Gx29qzOPDPSLLMw/9oW/sXR6XW0Or1dJk5k93vv89RNzQ/Bfp4gU90CcZSG628ksimOW++dBdf9ZjLJ0B3HfhGgLaNgF4yhdu98P/Uypvb57Pb2AOXw0z64XAMCeBIB/t+G3HuPzGU10iUXBS4wXgMTAE6RHcXIen5OFsn0vamw/Q5vAr7Zkh79gNOH26OuSOcOabBt7WOkydNtGrVlaSkw3z++QHmzetZ35dPpTwSUFrfg2jc1EgQmyzLBbIsr5RleT4QCShLthpACGERQtwthFhRyfG7PXvcFe5zCyHKfCFXryhnseCzcSPBPbsRMH4U3v95A+2J4+cjyBoZNqf7XMWz7GI7zz8PmYU2IqKMxI6Ff9z3FLf1WkLi6mimvaFhXlfooQcToCWbI2u+ZfvfHsfr2zlsmDuLHQ/+ld9e+Qcnt/6M1WrFDPQywLyucMszYPg+lMWP3sk8w4d88M5fQBZETxlQr9dApfr0uOMOAI64wZYMUqYMaRB8No8oEvAjh8CVa/lZlgg3T2VJ11i2tp+KJlrG8Y0On7dK+WrDNM6ER5FEK5IOtERn0+H4FXxKXchJTjS+mRykGwAH3OCM1yGXuiAeQjKK0UguOpDIrr92YcTaX2ivPYTQSAR4bq8Ww4wAWCxDAHj11e3XnHqpUktIgLUaD5VrW4FXhSzL+cCTQojHgNer05cQohfKD4Jcz9+Lj98N5Hr2uBFCRAohFnl+SJSNZ5Pn2G4hBGVtr2guU6eR3mcAfiZvDNu2YvxhPQFTJyJ7eWOfMAn7xBtx9h8AusbrEpZR3N8dOvqw///sdBnUm0NLvuHWvyzh9qnLGL8liXFr3KTuh6QMyLeDVVbepwPMGgj2hZahENYH5LFaUoY2Z0nEMFbaZ7A9djjBL2ewbVd/hjELwsLqecYq1aVFv36E9ezJ2b17OeqC3kfA1BaCWxeQEZFHyLHD6F/5jt1AQtS9TPk6htW338v4Z9ewPnIan018HFPCKU52aElphj9FaRb4BfROmeLtpfj6llBUtJmijuPh+D4OCy2DrC5cZ0DfzIS2QymWnCIiQxI5aOnIrj8OZOZ/P2PfLb1JSQpE1wyOZ2kwBgoOHPAmLCyC48eTWb/+JDfe2KG+L59KGWV74NeKd00NpPFSm3ng1ZbaKheAVtne9fzyBtmzxz0azq3MF5VrG8+17subTDjGjKPoX2+SfeQUBR8vRfbxwffJxwlp1xq/u+dhXL0KUVh4Td03FAKDjHCkFX3fsfPqA09yk/sbnrz9BbZ+Pgzdhub0W25i0mINc16Euc/BrFdh3Cdaold641rfgp8W3cCLty3gjsCPeXTT22yYdyMtZ29j6/bhdP7Hn7Dzd9Co0gONHSEEgx57DIDfXVByDKQE0Ca7CTiZjeuWDQinm2gtTM19j9j7D3BTl/fIWeLHl30e5paHI4hfP4xip4WU/RHokgTEg3OvDWGTKCr6DpMJ3lz5EKHdu2NzuYmXwBGvh1wb5EDIWQdaSmhBKnsf6sr0j9fQLVQJumwzHGRZ0HGsCYC2bZWEmTff/K1erpdKJdRiHnhT4aqXjkKIU0AeilTqRlmWt1TStFZDjT1R5hXF+ed73ObxQFy51yOBl2vgxLh69MTVoyfWvz2NJiUZY+wPeH2+FL+/3sNahmHYNB/HqDHQKOVCBe++HcU3/nasvxl4s9fDfHbLnUR2O0XX8QdpoU3FQj5a3DjRk4+FRHcbEoqjSDwVSc4vzXB+J+B/x+lXtIM1PMsTvMZ///oHeLy+56ZSU3SdPZvtCxeSceAAW20wZhsY/WHHW7lYDygCQcMMEDBwDTeNXkNJPxPHo1vTTLOXN98NI3uSm4y9rSBPj2sr6DKduJKcaLQ7cbvTefvtyXTrFop13jxiH3yQw7KGDqlOpCLQnNZhDs7GUOpPF6+zbGjdigPjYhj95TK29exCsvACA5yRdKCBgweD8PLyZtOmeI4dy6ZTp+D6vnwqcD4PXOWauRbf75PlHgs8mtZ7gE3ALiAfxbDWdspWpOdcF5MLRMqyvNizLx7pabtIluUaL/AqtYyg9K67Kb3rbkRREStarmPsM08hnv4b1vsewDZrDhiNNX3aOsCfXVuMHNifzVNPZ7Aroy27evdH2wH0LZ1gkJEdGlwZWtyntXDKCcfSae37P86c9mMyCXzIM9zJJ6xnEv/FCfwIjK3vianUABqtlikffcRHgwZxwOkkLwlK3occFxiByXrwGwzanuDsriGhUwsOamPQjimiYGwm+cdbU3TQH2IFItuNK64UjeY4bvdvPPBAP+66S/ltHj13Lj8++ihn3G6sWjCeBmOwgN4QlFmMo3US3nTl6P19mHP7ClZueZC4Pa2I6A/J2wQdhxs4vtnB4MEj2L59He+/H8d//jO+fi+eikJZFLrKNXPVBtzjsi7bcx6NEnk+CiU3vHyUyEohRHdZlvfXxEArIJCKA9Py8eSgX81+N0BGRgY9epwXkZOBKdOmc/u8K9ek+ZQ7eGrdGLx/3Yblg0UEPf8M+bf/kYLbbkcKKHNKhFOYX1lMXcXHrMUV+YwubPvGQgOPLsitpL3SdtuvBtpFa1i3thgwVXCusj6Vv21aa/jis+Z0aB3Gp8/uZecOK+8/p6NXHxsul4vAIA3t2gl63eZN776hBFjas7jNMh7iP0xiHXGeXYu83FzgILm5fZCkpiOsUVCDGvMNDUObNoz78EM2zZ9Pss0GkqLqN9kILQeDYSzYbtAT3y2C3aZe/E5fdopBpB5vTfa+MNgo4KQLeVcpuI8gsZH587vz1FN9yM7OPneeVqNHkxgby3E39D5jwNDdgUiFIO8S0lsFECWncHBAB2w+RvoeXMduzXyk9jJs0yJFaAHIzGwDwKef7uPRR3vg5VV3cSvX8z1wJVQ5/2rEFbpcrgvuk4ZK2fyXLl3K0qVLLz5cLXdQte5iT4BYWZBY+VSyUSjpZDM8K/RNKMuvVbIsJ1bnnLVJaGgocXHnve5Ot0Ri6ln8LFe3G+AXEASTp1E8eRqlRw5j/u9bBN4wFNuMWZTcdz8QXmWflR2r6PXyry1618TzL9sqbe9nCeSXHYoAzHdvB1y2zwuPaRg3qSvjJsH775hY/5ONl5518vTz52uV/+cFB88kz2cOJxjA/0gh4tyxAE8RgYqKCVzvBAdfvy7b4Ntvp+v48Zz69lt0Rw7S1qsQU/sSpI4SOe1NHPMN4oRXZ466o0gq6IjhZDBhB70J+E3CHe9ASktH1+oY3boV8Ze/zGXMmKhLztHnzjtJjI3lmEZPrzwHUjZoE73Qty7Cp7CESP80DtKBJaMfIOu1XoTNKiY13A9TGJw8q8evuY2TJw106dKFI0eOsHVrJrfeGlO31+k6vgeuhIrnLwPOa+5Tp9M1musaHBzMI488wiOPPHLB60KIav0CqbGfoZfJDR8FvAb8HxBUU+ek4n32BqVA7+7SlcL/LkLzzPN4LX6fwDEj+IYh6OIextWnX52Pp0wA5sUXJT79b/X6euc/vjz9vPKDQXvyBHPeuA3mdGMo2yhVQ0SbDOZmzeh+94WyDxqUD3rH7GwGewWDFuXT2t/zuEsLeKHsblW929bhxhsx+PiQUVxMvgEMCQJtqB3sEJTpoNg/hyByiLXdxj3v+ZL2mo30KX60HQJHV0LHUV7s+ryEkJD+wBE+/HBPnRtwlYpwo0ajVY9aCwkunxsuy3I7FMW2mtSni6NiYx2IsiffoJBCw7D+/TmyDh5jMyOx3HErlskTMGzZXCt55W/927fC188LwNTQiWQZ0+dLCRw3inf5C4Xvf6gab5UaRe/tTadp0wA4LoEzWY/skuAMWHLyEZJEBzIJnbyLpQ+nMcDxFVqDk7OBgIBkFDf6/v0WvLyMbN16hlOnrl4WQqWmkYCSajxULjHgHhnUGsdj0PfWVP+efPN4TzR6eSxlud8NErOZd7if7L2Hsf3hFnyfeITAkUMxfr8GanBv+J3/VGzAr5bXXq7cSRPKWfxvnYP3u2+T9/0PLGZ+hZH3jz6p1AUfP15/yTEVlSuhq6cy1AmtAbnIgTsDSDajcdvwy7PSnDO0HnuUgf/+kb9teJ4uLc+SJ0GzrnC2QBDVX09+PvTrNxSAJUtUffT6R80jqy4VrcDrInr8aqls43Qh5SRbPcIvDdd4l0evxzb3VnJ+24P14ccwv/EqQQN6cwefgK066gbXRmWG+o1XLn1dIGH67FP20x13+w7kbvkVV9dulfb9+N8UA/7DDw/UzGBVmhxRY8disljItjvIkcB5RgPpJWCF4Cw3GqyEkkFxRycnurdn5q6vAAjprbw/MMYAgN3eEYBPP92Pu8Ka9Sp1R5mW6rU+VCr61hYeFbXaQHAV+eGeFLAZKPvovYQQC4HTsiwvBvCkit3tiYa3oKSPza+8xwaIRoN9yjTsk6di2LKZ2Te9TXC3Jym980+UzrsLqXl4nQzjjVd05wxtVegOHWQrj+D9iY3hxBL7XOcK25nMEs1bu9AbZLZsgUnjanrEKk0JrcFAp2nT2LdkCSckCE7WIfdzIM6Ar282GpcvnXQ5nKU5m+8bxeRXPuPFl+7nVIkJ9HA4T4PWCLt2GWjVqjlJSen89FMCY8deGjSnUle4gcYtflXfVGTAF9TyOa+4yIlHWe1Vz6OyNotrYlAAubm53H777eX6hj4DBjJx8tQren/5FK2hQxxVtKwkjaxXTyYQS8KXW7Es+ZjAfr2wjhhJ/rw/8caG/udSxCrvp/zzKzlW1Wvn/65bW0zXaCe/zv2Yqb+9xOe8RNiKieyLjKAwP62Cc4TTMtJFtwEObnu0iNX/8mFQ3wK+//knxo8fr6aRNSFqcv4tx49n35IlnNQaGGh14E4HTbIZTRcr/rnFuEPi0YpOWCe0JOTBbGaWbGaZPJE2/V0k/qojeqSWgz+4adeuP0lJ37Jo0W/06uVfY+OrDPUeqGz+1atm0tjSyDZs2MCGDZfsIFfrBrzEgHtqeTdJAgMDL8jTu9o0svIpWpNurLptVX169e2PvW9/HC+9gtdnS2hx/1+4N8lCWOgsbDNmIbWMqLSf8s/L/r1lC7SLdrB7b+C54LWKUsXK/r17b+C59m1IYMdnVl5YEcS6h0azf+5sFr8VxgtBtirPnxKvw+2GI7sCef45mYDAQNLT09U0siZITc0/YPp0NgUEkJOXR7YBmp/R4dXCCkUQnA15zRy04CzJBj0r75nGH1Z9yLLpEzB01sGvICK8gGLS09sAsH59AjqdDxaLqUbGVxXqPVDR/CWqswJvbGlkt9xyC7fccssFr3/22WfV+nWnClPXIJXV6L5W5IAASh54mOz9R3mQN9GePk3Q4P4ETBqL1ycfEU7qFfUTu1H5YRG78crOGxvrZuGKHDb95xC76MuN/r/zzExfxtwTxYtvXVkxEptVQ8IRAyf2GWsu4l2lSaPV6+k8fToAJ9zgTBbIbuCMwJyfgdbppiM5yMgc/FM/HD8Y6fp1EckHJPT+cPCswBKm4ehRQd++PbDZXHz99aH6nVSTRt0Dry6qAa9BaiJF668PVRBdqdGwleEUvf0uWScSKLnnPgxbt3CAGIL69cR3waMYYtcTytkKU9Iu+8PCamUov+D979eJZRyzV/yJL2emMLF3Lm1IZODiuRw+6MOIEcredrtoB1vKKeBfHABX1sZkPu8q12oaoy68SkOj66xZAJzQG5FLnLhSgWQzQnZjySnCn3gM2AkPsfJ598d44RkrLY646HCDUuCk8xhF1jggQJFqXbKktoQiVS5PWR64GoV+rTTeOpjXKQ88XESVcX5GI/bJU7FPnkr4aj2Z7/8Pw5bNeL/7Dof5MwFRGlxduvJfumB+JQgpKIgJAQGsOujN+HwHYqWLh8jF56kzrCKJoJ4H0Kan8SoxaDL78j73sGjfYCZGNOOdX21Y3zBRvuZfy0jXBUpucGkAXFmbBTODOHXQgFGnwdug3moq1aftyJF4h4SQl5VFlgHCk/ToWxVDPoRkQU6Ym9ac5SRGbFMzWXN7EWmTBxHQXHn/WU/p3z17/PHxMfG//6Vw9GgWnTuH1OOsmirV2wNXUQ14o0ZCi6tXb1y9elPy6OOE+RvJ3JGI7vBhjm5LAOksuqNH0eTmMAuBabUbNBra0AIppBlfM5BhX/4dd1Q7Bgb78uVYG4e2adi8u/Igs5R4XZVKblu2gNsNf5sTSEq8cnv5eyn532azuVaug0rTQaPT0WXGDOLee49jbmiWJCE7QSRqMVnOoreb6WDM4CStCb+viAeXPIbw+TufuucQ0BYSEgSRvXTE73ExevQwNm36kY8X7WRa+DEOfPIJuadO4YVElEZmQIgGS1eBfpgGKUZPYScfksNbkmhuyyG6c0aK5KSjI/k5zSg4HoKUZEJ3XEYkSOhznOgy8xGkYzYfpUMHb6ZM6cjUqR3RalXHJ8C4cd3Izl5+ze9vLPvftYlqwK8rBFJYcxxhzXkHE0//3/mV86zVJs4uVZ4/tNrEnIdsLH/WxFudzrcp2ysvv7q+GJtV49kmqDhXPXYjvP6Nsvq2WZUvKo8ePvfcc0+1Z6iiEj13LnHvvccJnYGhdgfOJDAEGhHdSwjKKsHZMgkzXWhBEp8+cguPvfkan/5tNi0HCfISILSvkfg9LqzWDgTzBc537mKTO+tc/0XAPuB4qpupWdDyBPwSPYHV4l5ukt5jfNYaRL7MibDBDO+5lt3/7IRteEcyj4STGxgCFi3s1IJ/COw2Q0krdu78hU8/XU7v3s1ZtuxmOnSoSUXpxklsbGx9D6HRo/4UvAyVSZJej9REEF5ZH2WrbxWVmiZi0CAsbdpQZHeQLIEzyQRFJZAJwZk2kGU6kI6LQnbP6oU5L4cxR38kwQxo4WCuBr0XxO88zp81SwhwZxEkYJoeHmwHfxoNbaIU5+43TshIhVWp9zL+1RhW59+LcY6Tb/r9hamLOnF6/w0s6XoP43ZvILRbKiHt06EzMMoFgVp0w73BxwiMwtd3FLt3pzNgwIfs3p1W9SRVVK4A9Vu2HBXlgX/++VIeePhqPmxVlQq9fJuKy4ZW9p5rzwNft7b4opz1cHr3zOXUwXB690yjqrxwoNJjZX2UpyxXc926dUyePFnNA29C1Nb82918M3FvvMERNLRKsiFZQcSb0IfmYi7ypbVPIns17ehJGu8//SdefutZ+nQeS9sBbhK26xg4NJ1BP07HKFlpq4FJevDpB/ruGizDJWa1MLHsNS1p31gpccFNJ97jmzvuZUb4e+AHs/u8zbK3HqL/glX8au3PSxNfwP21hp9GjsaRCwWEEzrJTcb3WgLGmyhcX0JRUQwdO8Lx4z8xbtxnbNgwk9at/Wrl+jQk1M9AHeaBN2UqygP//POqc7Yr4kra10Y50YufV3Xslx0+l+SsX5wbXtnfil4rnzt+MWV7VdnZ2WoeeBOkNuY/8J57iHvjDU5JArsExlNgtLihN4RkOLH6FRDGWTK1etb9YRTzX/6UKb+v4dDAaeh+sdF391x8yaCFRnCjXsbcHwzRWsQoN84W3nyWeRs/hfyRiUGvo8v5BqfjO/69bw1egeDco2fCjd8yuMdPJA5qTlZxEC8++Twlj/cg9C+FaP8I0q+BZGAiajCc3q6jx20+7FtUzKlTMfTvn89vv+1m/vxN7NjxJwwGbY1fn4aG+hlQ88BVapDqusvLxGHK0snK9s83bVLTxVRqn8B27WgzfDhOt5tjbnAkmpDtTkgAS3YWWpebaDKQcNFFl8FLbz/GjKVr0S2zMST+YwJz9mAWWqboZUydwNgFxGA3UqiBU11C+XnNbcx+uA3rez3KVxoN+50ymRLkxHmjszvJ3RaAn7YI0y4ngZZ8jpwcwaRlXWn1oYTBXsqAAcmATE40+AVI7EvTMO4uL9xuyMsbRuvWgezenc4//7mtvi+lSiNGNeBNlOrmrF8sDmMr1LJgZhA3T1YrjqnUDb3nK2UP9mv1SLk2JSf8lBkiwC9XAAAgAElEQVThdhKUUYQ/pzFTTBSn+W1UX77rcA8vP5uHXh8FaJimd6MP0uPVD9ydgVaQ2D4EzV6Ytf5NPnsqjWOP+XDqr3NAhlidD96FJRQeMRKYm4c7RUOU6wyGYjtjpy5h0X/zGBDwDU89/R9STMUM6phPvgtipisplntKdbTrrOHECRg7Vsln/+c/t3HyZE79XECVRo9qwBsxFYq+1BEXr+C//kzPqYMGJoxVbymVuqHz9On4hIWRY3eSJIHjuAHyrJAKoWnFCEkimmRcFNGGdBLnlPL9rIMMjH+P4TqJQKMO/xFO7IFGtD0hJSgY7XcSraankfeRHtfKVJxj3eya9Sx2L19ySopJl8B+WIvshOxtwWiFBL+WMGTCTzzz5j/JXWpk7Jcb6LbrMMb26Wg1MjuFnl69IKtQ0HeONwBffunFnDl9cTolHn/8CiUSVVQuQv22rWMudj1XB0X0pX64khX8s89e+DwgoPL0NBWVq0VrMND3r38FYBda3CkOXJnAES909iICsgpowQlMlNLLfZx+XzxM/wPjWGBdQ773VBb0Wst6+3RMI+2s2vEH/v7gV+x5eSjv/3Inv48bw2kiMGUFcCKnPUfHK+fZKAIw2UvIO+JNqDsTR4qGzj4ZGK0uvMU+jMEafvznDYy5L47ch0MZXFiIG0HwDcqY1xzUMHqSjuJi8PMbgre3njVrjvO//6XU01VUacyoBryOuVpd8oaOr6nyOMjnnrvQiM+bN6/2B6TSpOj7l79g9PMj2ekm2Q32/XrkjFJIgfDkYjRuOz04BS99R7Otu8nXaSnRw9q29zJhUQzf5dyN+zUNscvvYPJbXXh/3LOcadWGX+mKj93Ivt8j8BeCTe4HcGkMZNvzyZfAdlSHLEHGDj80QiYi3QxIDCSXPX+MYYduGs/8xYHtNxmdkNhUBCMGQ4kdIm5Qiqd89hncdddgAJ5/fms9XkWVxopqwOuYmi54ci08+uTl635fKZeTSH3uuRo7lYrKJXgFBDDo8ccB+FnW4khx4koA9hrQW4sITc1Dt3Y3uud/AWCScKML1jK91Xv88OABbgp6j9xJfvR8dA0fvJtD4IR97PQeAC49+7e3RePWYtsAxXIY8SGzEMhsl/zwsRVSmqShpbMAisAn6yhapwszuzFqJUIn7WL1nSfJnZTF2Ha5SEDgAGXMq/doGH+jjtJS0Ol6YTbriY09xb59Z+vnIqo0WlQDXsfURMGT6lJet7wqqnL3l8mjXg3Ll1+7bKKKSmUMfOQRAiIjyXa52eaC0t+1SCkO2AWOtZnsnp0CMgzSQSsvge8Nbsb0WMsH/zeengt/JeWh5ujnFtLz38vJnORPoeTLmbhWOIq9MO0BexI0dznZmqF4kI65SnHJUHwmCIEMSW0Rkp3wzABk7PSnFO3/FfJ6yR8Z5fwBrzaKyltsMfTvBvlWzhVV+fRTuPPOPgC88cbO+rmAKo0WNQ+8HBUJucDSKxBmKU/dCblc2vZKhVyqan/+77r1ThauKGb1v3wA0wVCLsUFeUCwR6TlfI5n+edlAi5lHDt2jNzcXFXIpQlRV/Mf8fbbfDt1KntdLpxFbqLXwtk98MtJcLtlumignwHMw2VEK9De4CYvxJf0tiEcoTOJtOWgJpp0wsk41IKizEDY68K6R4dPiZP0zaWYvbtxtiScMNI4JWvoeCYHyQocz0HTCfzTUkgO9yNEOgjagex4rB93L1zCrZOGMrBZMTszfYgcVMpvh7xYe1imew+J/fs0BAR0RaP5H19/fYgnn+xFaOj1VTNA/QyoQi51QlMScin/WtlKu0yEpazdpIn6c4VLQgNd7NmniLWkxOvOCTNcLNBQ9vzZZy89ZjabVSGXJkhdzD944kQ0X3zBN7feyiGnk0OFQKFyLEYLI7zAPBJ03YDhkBvhR1K7cA6JbhyRO3NAdCdeiiTjUEvyEkNgN7BXB2lOiveWohEFWEu+YS8xTCCNE8HN6JR7FmdWM4zmTMgIQB+ehl9RCIV+aXRAYs+c9tzw91/ovnM/ft1C2Znpw35vL8IC4WS6jidvMbN/Xynr1/szeXIH1qw5zurVZ3j66WG1fr3qGvUzUIdCLkKI61/fT+UcFQXWPfqk6wJ3/9+ecrPpJw0LV+TQMvK8C/6HHy+t/Q3q3rdK3dN11izujosjes4cmjVvTpSfmenBRsb3MuA7z4hmig7rDDPxPSLY1b4HmzUj+JnhbGU4h7JiSPy1A3m/h8B6YLsE+0thdxFI+5CkL+je3chzK19AaLUkZGRRKoMzXiieuhSlHGlYphcAXchG0mk4+PBIbn3rK4755NLMIHOkGCaMU8Z7RtYREAC7d0uMHTsQgA8+2IPb3XQ8VCrV49wKXAhxEtgDfA3EoxTkUWkklN+vvtr99bLAuvIlQsvvkxt1Gvy99EyfJPPARYVKVq+7sPY3XJo+pqJSV4TGxDD9yy8rPW4GIoE2SIwBbMhoEBT6adD2B3tXcE0FuQRycw0EBkqkpralWbNuREUFIIQgZ8wYTsXGctLXTExSBlI+aE8nQA/wzj6CiGyOQbMHHeP59fYIHnj2e17JPMv4Tq1YesAPe1tlLKt/E/zxFgOL3nFw+HAz2ra1kJCQz6ZN8Ywb164uLpdKI6f8CjxBluXZsix/I8uyarwbGeVX0Veba365wDqLtwEhBBPGKqVEFzx2/ra5aZLukupjla2+Q0NDr3A2Kiq1i8bznzdaTEJDMyMEGSE8CFqFQOvWGnr21NK6tZlBgyJo1y7wXFncrnPmAHDCrGxfOgtag9sJGS0R7kKC8wKRsNILKLEY+HzIQ/jc05z8IyUA/JAHw2LA5oCwnsqP3i+/dHL77b0A+Phj9etX5coob8D3XOmbhBBtanwkKtWifHpaTeSaa4QgwNtQ4bHyBrrMqJfV/q6KW2+99doHpKLSQOg0dSoavZ6U9LOUyOA6Uaq40ZMVlbWgbMUF3gpFnCU2+A/8+ZVAkjfoibFI5DkhZojS18+nNHTvriEvD8LDYxAC1qw5Rn6+rT6mptLIKP+tm11pq0uZUdMDUake5VfR1c01N+m1BJkNGHRqlqGKysWYLBaixoxBliROe/sgpWUi5WsgMR7sYMo9hJAkBPsxIfD7wx6+ufM0EeE7GdReiVlKDQCjHn45DJNnKSmZGzboGTGiLXa7mxUrDtfnFFUaCeW/occIIe66kgcwu74GrHJ5KnKJX4luukYInvq7hL+XHo2m5quKff755zXep4pKfdDpppsAiPdV5IGdxa1AckF2a4RkJSSvGTI2eiHTauxx7rz9a+7Z8xbukHwANmTD+H5KqqqmhR4hYN06F9Ondwdg2bJD9TMxlUZFeQMeBdxzhY9edTtMlepSlW56makO9jHw0gvXvuq+XPBaRkbGNfetotKQ6DhlCgjBmfR0HDK44p3KgVTFjR6Yq7jRW5AKwKnZbRm8YSfHi5LpH+SmxA2RfZW3/HBAw7BhWux20GjaYTRq2bo1kdTUwjqfl0rjovy39SJZlvtcyQN4rb4GrFKzGHUaAs3KXndZkM7VYjIrqWT9B6npLypNA3OzZkQMHIjb6SLRaEJKSEUqAhISwAWm3MMgywj2Y5AFSQH5HB0zgOErfqRPlGKYz5jBbIJdJ2HMFMWNvnatYOLE9sgyLF+uutFVqqa8kIt8Fe/bVdMDaQhcz0psF7fVCgE0x1VSSH4JnFdVu5iLXz///P77NWRnS7SM9GPhihyWv+lH316llYwfrFarqsTWxLie599y9GiSd+zgdEAwHTJSsBeG4+WbhjsrDG3zs1jyupAfmEdHm5ODXjqS5w5m6uvf8t9bbwMCiM2SGRnjYO3vRgrMTjQamY0bnbz2WhtWrz7GsmX7ue229vU9zWpzPd8DV0JdKbFFXembZFleVZ2TNlSaghKbRgh8TTpMei1woUJSZWpJlamtPfNMNsHBwaTESyyYGcRbr+sIDq5cBlJVYmuaXK/z7zl7NjtfeIEzuQVIMkgZJmgB2uxm0PwsIflG8gOhrTaHg4SSfWMLht+TTG7aCXoHtGd3nobOg4ys/R1+PW3mhhusbNnixts7Bi+vTcTFZVBaqiciolrf8Q2C6/UeuFLqQoltthBiZHU6U2m4aIXA30tPiK/xnPGuKWxWJZVswtiq988jIiJq9LwqKvVJcOfOWNq2pbSoiAwE7qNJyE4gIQ1k8M47BbKMSXcYDYJUfSo7505gzNLvGdC2GIA0PzAZYOcxGH1jmRtdZuJEZeW9atXR+pqeSiOg/DfuKCBACHGzasivDwRg0mkJ8Dbg56W/xHDXtWLarFmz6vaEKiq1iBCC9hMnApAY1gKcLlwFgVCYDYX+aGzJ+JT6IWly6YwBCYmc24czdela9MF5APyQBaN7KP1pwhWH6IYNLm68sSugGnCVqjlnwGVZ3ivL8iqPe3x3PY5JpZroPClgwT5G/L31leZzq3rlKirV45wBdylxHa58RROd7DYAhOQq7u8oFE+po4cZq68Z6bcttDfL5Dqhaz/lLVuOahg4sCwavS0Gg5bt25M4e7a47iak0qg4980uhDgphPhaCDEdaFuPY1K5BrQeox1kNhDko9Qaro1c7urw8ccf1/cQVFRqlDYjRqAzmTibmoZVBtfxLI8qmxUAn7xMAMycACBDJPDekAWceTmGzqmK2lpBMGg0sOUgjJ9ctgqHMWMikWVFmU1FpSJULfRGQEXa5gLQazX4mnQE+xgJ9hhtnbbhqqfl5eXV9xBUVGoUvZcXbYYPB+CMrz9yVi5SoQ5ST4NNg7bwEDqnFicnCMcLGzaOyRO5/c3mWH9zAPBjLgzuDE4X+LVT9sHXrXMxZUpnAFavVg24SsWoWuiNgDJt8w0bwdugxeKtBKP5GHV4G3TnVt+1gVpZTEWlatpNmADAGX8lw8JV2kbJQc3rgMBNQHYzALqiGOzO05NYPj+RdsE/E2yA+BIY0F/p6/czGrp21VBQABZLezQawebNCRQUqNroKpeiaqE3YLQagZdBy3RPxa9ZUwz4mvQYddprFl25Wi7eJ6/IoKtGXqUp0278eADOZGQhyeBK9Rw4awLAP08xvsEkAdBy7O+0n7eJ2XveYEwLRcFN9iRorN8Nk6cqbvStW7UMGdIKp1Ni/fqTdTQblcaEqoXegBACDFotfib9Obe4n0nPlAm6K0rTqgsqCnxTg+FUmjKB7dsr6WTFxZ50skRkB5CQCBJ4Fx4HWcbNfrzRk0M2JdNH0nvbXrp4nQZgRwl0iYACK7TsrrjRv/vOxdSpHQFYs+Z4Pc1OpSFzsZDLPVf4vp4Xv+AJfitbFubJsry5mmNrEggBRp0Wk16DUacl216Ml6Fm87QbCpGRkfU9BBWVGkcIQdS4cex+/32Smrek+dlkXMXN0BsyoTAErSULi7Ub+T6ZRKPjN5xE+pSy9cahtFv/JYbI59mZB/f3hyPJcCJfQ3i4ICVFpl27TsCPrF9/ErvdhdGou+x4VJoONamFHlQ+DU0I8eeK9sqFEKPUPXTw8uxlN/M14e+luMXrk7pwg9/kqeCkonK9cc6N7lEJdhV4lMeyFd94UJ4XABGenUotyXx32yS6fb6KG0IkZMBXWWzz/S7BZE80elycFzExoRQVOdi8OaFuJqPSaChvwKurhR7nSUPrLstygSzLHwBRF4vCyLL8EzD6Gsba6NEIgZ9JcY/5merfaJdHdYOrqFw7bUeORKPTkZaSiu2CdDKlcIk5T9kY13MYgSCNJApGDSI0OYPRTiV++IiAEH9IyIDuQ5XviTVrXEybplj2b79Vo9FVLqTGtNBlWd4rhLgb+EAI0RbYBHwNintdluVvyjW/muogdUZtFTMRKCtuo06L1Q6VFw6pSvj/0vdc2rZ8mwvb10ZBgQv7rHxOZSxZsoR58+apxUyaEE1p/s379yd1+3bOePvS8WwW7mI92rTTYNOh4TB6e3ecxiTaOIeRoLfR3Z3Hurnj6f7dxzCkDxsyZW6MsbN8m4mEEjs+PjIHDkg88UQLAFavPsqLLw5ocPoOl6Mp3QMVUVfFTGYLIVZUZ+9aluUCYJbHgM8HNnsGGC+EmA3EAwHAqWqMudao6WImQoDZoMPbcGnUeFXi/ldaVKSi16oqTlIbBQWupBhKGUIItZhJE6SpzL/T5Mmkbt9OcnAzOiYX4ba1Red7Ald2FLqWx2lW0JzUZgl00dtIAMKMuXx32yTenfYoMePf5kCxlk79TLANfj5mZsKEElascJGV1Zo2bSwkJuZz6pSNQYMaX02BpnIPVEZdFDMp00KfXl0tdFmWE2RZflKW5UCgHfAkyoo8AHhCluXXq9N/Q0cAZqOOYLMRs1FXZylfKioq9UfZPnhCdg6yDK5kZVdSTlPc4f55iiSqP/EAZBJPbvcurDNPJvyDLEwHJdL9wMug1AgfNlZ537ffnnejf/ONqo2ucp6KtNC/qckIco8xX+XZE18AXPdLsGAfRWSlsbm6qoOaC67S1AmNicE3PJzivHyyxPl0Mm3SGZDAkH8AZIGDAwThRTHFDBMOVnR8gD8/J9PyiIv1WTDGk+Njs+jQ6WDbNjcjR3YBFFU2Wb6acCWV65krSiwWQvh5VuYvCyHe8/ydLoTwu5qTeYLbEq9ppI0AH0+KR1My3GWoQXAqTR0hxHlVtuYR4HTiKglHYyuCgpYIVz7BhS0AF9EocSCtyODsrS7W/eEI2nZWUm3Qs6/S348HBCNHapEkyMgIJTTUTHx8Hvv2na2nGao0NC5rwIUQLwN5wAqUFfR84AnP8zwhxFdCiNa1OspGgK9Jh1nN0aySTp061fcQVFRqlfaTJgGQICk/4l05AcqBzDAAAnOVr9wwUgCwkoB2usykHm/xUuF/AShsBlpPcZNxN5a50d1Mm6Z8flQ3ukoZVRpwIcRyFKO9CmUfeyaKAf8A2Iuy3TsLJUjtn7U71IaLr0nRJK9vHn+8pMJ/NxQmeb7cVFSuVyJHj0aj15OWnEypDK5DqcgSkKhUJfPKVVTZJPZgQk8mmQwDlt03i1GrFqFxu4nNheHR4HKDrqUOIWDjRhcTJig1wpcvP6K60VWAKgy4EOJmlHzt3rIsz5Jl+bWyvWxZlu/xCLoEoKi37QOeFELE1s2wGw5+Jn2DMN4ATzxRUuG/VVRU6gajry9tR4xAliQSAoOR8/JxFflAVhJYLWhsifiVNkOiiB4oq+vOZHKkdxdywoOYu+c7jhbDoAFKfxsPaRg2TIvDAfn54QQHe3PiRA4HD2bW4yxVGgpVrcCfBP4sy/Leyhp49rQXy7LcG8WQjxVCfFXTg2yo+Jn0163saW3w1ltv1fcQVFRqnQ5TpgCQ4KOk+NrzlTxushSpjZAcM6DsfwM4SUCP4L+P3cJz3/4TZBlbCyUN9ce9cONNiqFfudLN9OmKG/2rrw7V2XxUGi5VGfCoigRbKkOW5cUoKWN9hBB/qvbIGjj+XqrxvlqcTmd9D0FFpdbpNHUqAPFJKThlkI7lKaJQiUUA+OQokqgadqNHRxopDEHw4/SR+LiKGLd3A+vy4IZu4HCBoZUOjQY2bHAxaVIMoBhw1Y2uUpUBj7/azmRZjgfGAq9ebYR6Y0GgGG+TXjXeKioql+LXsiUt+vfHZbeT6OsH6ZlIRSZIOQGlPmisx/ApDcBFNj0wAhBDJrJGw8q/381LXz3HkUKZIYOV/tbt0zB6tBanE1JTw2jRwpeEhHx27Eiux1mqNASqMuDXJHfqMeIfoAS3XVcIwE813ioqKpeh8803A3DSoiiQOYvaKAcyOwAQmq2sb9qSBoCDUxgRvDvnBg6LwQxenMzRJAm9Djbth0k3K270zz93MXduNABLl+6vq+moNFCqMuDV8c+8jBKtft0gAH/vmjHeDTFCXEVFpeboOktZv5xOTsUhg/Owx41+Kh8An6wTIMsIduGFnrOkMxIZWaNhRb/HeORfek7+YmNCb5AkKPLTYzbDjh1ubrihBwBffXWYkhJ1W6opU5UBv2bFNI8m+nWjZlJmvGuqelhTjRCPjo6u7yGoqNQJltataTVkCC67nZM+/shpGbgLfCA9HooC0ZTGYyluhpt8+nq+hjuhuMTzb3Hw8T05TDv+FjGDlXXUF9sEM2cqq/CtW33p378FhYV2li8/XD8TVGkQVBnEJoQYUY2+G2TFsWvB4m1oUKU/Gytjx46t7yGoqNQZMbfdBsBhLx8AnDktlQMpiu5VWIby9duCEwBkcZRQNOSOzaXTB6FMzV9Ohy//QVgAHE2GfhMUA/7JJ07uvLMPAO++e2Fl5x9+lLj7QQc//Nh0Kv41ZapKYLYAmzyFODaiFCPZJMvyvivsu9GFSF5SThS4+nKily+rWVV5veqUEy3/vKpjNcW19KnRaNRyok2Ipjz/sNGj0Xl7k5qSSp4BLHGJmCJAPnwKTXswZu1C07otTl0cLd3tSNFZGWXPZJkxmNwWOdx6/4dseGYMj5lOsiH8dg7f8xHfe62lk9AR8pXETXe7cIYLVj/4C+tO3Mnk9A9Zq72TCcu6sXLuYdq+/jg/TZ/GDxNnYytohi3FjJwt8E6RcWfJaArc+BY5kOylSFI63t6J+PgIYmKCmTQpEovFVCPXoSnfA1C75URFZakIQggJxWj3QTHmcN4oX9agCyF+lGW5US25+vTpI8fFxV3wmhBldcGvjCtpn52dXWl5vcqOVdTvxW3LP6/qWE1xtX2+8cYbPProozU6hoZObVz3xkRTn/93d93F3o8+omerFgzPTMU4rhXG5kkwrB20OkVu21GcCU/DzRhW408o4Xwo9aHdv5cR/fQifG0l7I+O5aYVMSybeYCpB8dz2jyVxKh7mZn4HhNL1jC/eyzjv4gh9qkDjBr5KWuO3kuvId/TRneMbh8exe9kMU8tfo7N0SNJ29saV4kBDgK7AZsMR22Q5ARKgW3AYcxmPU89NZTHHx+MTndFJTMqpanfA1XNXwix2yOKdk1U9X9mjyzLYz0lQaNQgtK+AQrwpIoB/9/efcdHVaWPH//cmUkmvQdCaCGUAFICRKkKSgcBUZo0aYJi2dVVcdWfq+6uinVVBMEVqSKCgPQuvYYAARIgEEJCCem9TLu/P+6oWb6AoGQm5Xm/XvMyc8+d6znP3PDknnvuOUcURbEqirJRUZSXFEWJLPP5SncFLoQQd1P7558H4OS1DIpVMMUUoFqBkwVgA7/LcSg2FQO78UPhWk4ij/R/jtYvfYqupIgQHXS5OIs1o2IZkT6Lh1zgfPjTDPyhFT+EPU2BCwz0mMXGv8YyuMEshrdewudTB9J96BpqDr7KJy+8wlMdN9D/iZ1M+W42DTqdwehTBK2B7mbwUKC1O27tdYA70IvQ0MEUFlp47bXt9Ou3mPz8UidGUNzKrRL40l9+sC8J+rWqqkNvM6FvQrtyF0KIaqtmq1bU694dc2kp0X5BqOlZmC4HQVYqpNRCZ7pKndQgVEqISr6AS6dv0W88gItBzyAX6O9tYEKnn/imfX+G/OsnQjZD/Zpf89XQWFwufkWmFbocW8Msrz70LfiJoovuBGXl4BFnxkspIHt1S8a9UZNv+r3O2P98x+RZ3xAWeQYXUxbUdyHo4QLcjColNbzo+rwnbu5w5UoYffpMITjYky1bEunX7zuKi2W0e0V00wSuquqHtyj7vYTek9+63YUQotpq/+qrAMRk5pBpg9JDxdiKgJgiKIXA5FgKjxRxqOMSlPgM/Fx1jNFZqRfkgk8/C9amRnS9rFib6klpGsqjm07Rp+F46uauYrmXOzaTjauHfFBNoNtrozjHSETueYqO+/DgwEUs+PgKQY8dZeHO4dSakU3tqV7UK0glwGAmI8iL/s8quLvCzgt6Jv3bE3d32LjRnQkTJlKnjg979iQzceJqmfmtAvpzNzfsbpDQK9W9byGEKC81IiNpM2kSNouFDW7elOYXUhzth5qXC4eDObM9j91dkym5YiHUBR5XbPiG6vDrbcZS3xW3PqUUBbuR2KI22SG+HLdEkrlwDJbwIAqzilkfHIRPWh6x8Y1xs5VyZXNtsEGbnJO0iDjEP+aMIdL3CAX1vFnQ9SWeeduLGru9Gf3QNVwU+DEX/vqsVtevdun55+fuAHz0kZ733huNl5crS5acZM6cI06MoriRu5LAr6eq6la05UaF+FW7du2cXQUhnKL3xx8T0Lgx6Xn5LLMZOJeQw+l1OpZ/mc7S18FUqNJEB4/pwKsBePewQX0FY08ThUHuJEWGku3jy0Hu45xrY5JdGrF7nraC89lreaQrBhocTeBqajANTYnE7orEw6UI3Xbtn/hOrbaSf7SQGqNjmf9qNlP2vMNuNYUXIgsBmFcKkwdoS5jOP+7CC39zxWqFd95x58svtcVZXnxxM4mJ2c4JoLihckngdlvL8diiEurWrZuzqyCEUxh9fBi9aRP+DRuSbrawygw/XLFxKgX0wAMG6GcEz3vBvRsorUDXXSU71JvzreuS7hHEDqUbSWoYsbTkiLEdF+95gNgeI9GbTSz0bgzA1QP+qCrUOZ9EXo4nXZocJP5HD9zcrTwamkZYrwRarlxK+8hknvt/szhfJ4X2/ipXS8HlPggPgRMXoUZ7I82b60hIsHH+fENGjGhBUZGZqVPXSVd6BVJuCVxV1VfL69iicjKZTM6ughBO49+gAVNiYuj29tvUiYqiVkhN2gX5Mb6eDx3u9cF7nBeugz2wDXInv3sA51o04ETTSKIN7dhEf2LV1myjO/G5LblytD4XdzVha5f3seiMKGnxnMafRjlnOR7bhABjDod+DAGgm0cGis2Ib80j3E9NVFQOz+jHwMUbcDlwgFHtslGA2cnw6gStru8uV3j/E+058OnTTbz0Um/8/NzYtOk8P/10xkkRFNcrzytwcRdVhfnTv/jiC2dXQQinMvr40PXNN5l4+DCTr6bycHo2tS/m4r4vF93MfJQXCwli3RYAACAASURBVNH3KcK7eSaN/BN5gF2MYxHv8iZzlfHsVTpzxqcpZ5oFUjJQIfGFurR68mkA4sO0md4aFqmoKjzolYpqq0Woz0XqZbUGoD4HccON08GZxH30Aq8/9wFr3C4zrr4NiwpbDNC3HeQXw5aLBh5/3EBpKfznP3reeacbAK+8sgWz2eqM8InrSAKvJKrr/OlCVHc6+z/TruhwRYdegZpu4KqDQB/o849X0BuNJJ05SXZAEOqZBKz6NujNhSiX2wDgf+kIOtWTIg7RFe3KPH5kDfarvbA9bsT3ch5GHSy7AhOHaRNHzd4Iz7zkhosLLF5spkuXtjRpEkhCQhZz58oQp4pAErgQQlRi3rVq0fqJJ0BVOVo3HIDSk/a1pKIPghKCUniSOjnaINJa9qvwi0oym1r/janveHNgvYUn62v3tpcWwvAuYLLA4oM6Jk92QVXhvffM/POf2vIY//rXbkpLLY5vrPgfksCFEKKS6/jiiwCcOh5LsYcn1kMxWN1aQGEmpHcAwP9yPApGCthNZ7QFVRqOOMfiZ6/QoWAJ90Vk46qD5VdgzCDtuHO3woSpRlxdYflyCy1aRNCyZQ0uXcrj229vd1kMUV4kgQshRCUXFBFBo759sZSUEN9cu99tuhSsFR6LB8UDXe4uahZ2BKA+cRgwYOj1MyO/KeAfa/4fP+fEM7qOigqsLobBHaDUDMsO6xg/XrsK/+gjM2+88QAAH3ywF4ul+ixMVBFJAhcO07FjR2dXQYgq695nngHgeFIyqgrmbQdQjSGQegZKegAQfDUHgHw20IZGAOjr5bNldH+6vjeTvk2058Lnp8Bk+1X4Vxvh6eeM6HSwaJGZTp0iaNQogAsXcli+PM7BrRRlSQIXDtOpUydnV0GIKqtRnz74hYWRc+kSKS0jobgYk0kbxMaZYgAM6avxtdyLSgkt0JYnPs4xcl57gUEL1nEy4xh9a0CJDY66QqemkFMIe5J0PPKIAbMZZs2y8NJL2h/jH3+8X54LdyJJ4MJhCgoKnF0FIaosnV5Pm0mTADipNwJg3peIqujg5HbQdwRbIbXStWUqTKwnjHqYMBFW08zyKY/S4l+fMLaRNl/DrCR4boB27M/Xwl9fcAVg9mwzQ4e2IijIg+joK+zZk+zYhopfSQIXDjN79mxnV0GIKq3N+PEoej0JR45QHBiE7cwZrD5dwGaFy9rANbfUHbioNSklmfvQ5j2PJ4akl56FNW589Wwu4Qk2UorBEA51guDsZSjy0BMVpSMzU2XlSnjqKW1U++efH3Jae6s7SeBCCFFFeIeG0rhvX2wWC2ciWgBgTtKSNMcOgT4QpSiWkIL7tP05jCeeXOMaXQI8mHvfazz/toUmZ7RHxGZfhKf6aB+fuV7huee0q/AZM0w89VQUBoOOlSvjSUnJdWxDBQAGZ1egIsnKymLs2LHXbV1ARkbGHRwl6Hf3z829+cl+s7Ibbb9+W9n3tyq7W+70mIWFhWRlZWGzVZ+Rq+UR98qkurcfHB+DhkOGcHbtWo4nJdNaBdPWPbiMDcGQkUhp/iCMHj/hkXwBmuvIZTv3FE/jkMcZMkuOUjy2NWtGxdNxfCA7da3YmqHwYotsXPR+rD0Mrw0sIDBQz7FjNo4fNzNgQENWrkzgk0928/rrHSpE+yuaX9q/adMmNm3adH2x7585tiTwMgICAliwYMH/bFu4EIKCgu7oOLez/632uVnZjbZfv63s+1uV3S13ckxPT08CAgLueh0quvKIe2VS3dsPjo2B/+OPs/Oll8hKTCSjbUuC405gU7oDqzEmq9AUPPI34msbTq5+L608cjgEJLidY8xjXchJ+Ir7dydyodNa5qXAz6o/j3WC73fD2tgAnnyyhPffN/Hddx68+GIXVq5MYPHieN5/vw9G441TSnU/B4KCghg1ahSjRo36n+0LFy78U3/dSBe6EEJUIXpXV1qMHAlAvL/2LLj56DWtMHYrGKPAmk/NTH8ASthAOGGYMdOAy/zw7Aga7DjAuKIYAL5NgQn2bvRvtsKESa4oCixbZqZJk9q0bl2T9PQieaTMCSSBC4fp2rWrs6sgRLUQ+cQTAMQdPYbVzQ3r/oPYAu8FUxGkNQfA49oBDARSQiLt8NH25yj3e4XyzStPUOfjN2jlAxkmyPCHZnXhahacTNPRr58BkwnmzbMwdeq9AHz55WHnNLYakwQuHCYqKsrZVRCiWghp04bge+6hOCuL5DbagDVzRqhWeOos6DxQ8nYTXKzNzeDHMdxx5wpX6YGF758egu/BI/w9LxqA/ybD5F7ax7/aCE8/7QLA7NkmHn+8BT4+Rvbvv8SxY6mObWg1JwlcOEx2drazqyBEtaAoirbACRBn0kaUm3YcR3X1gMQDoNeycVBaPgC5bKQNTQEoIZ4a7j7Mem08D331Ou562J4BnTuA0QU2H4UmkQbq1VNITFTZv1/PuHHa9K0zZ8pVuCNJAhcOM3fuXGdXQYhqo9Xo0Sg6HecOHaakVijqxSSsXvdrhUleABjSVuOptsZGEc3QxlPFcoJH8WH5pMHsiqtNu6/ScTthY1kGDOuiffybrQqTJ2uPlM2aZeLpp7Vu9MWLT5CTU+LYhlZjksCFEKIK8q5Vi4a9e2Mzmzn7yzPh5+zLjMZsB9eGYLpMSLa2BKmFrdSlDiWU0JRrqEZX5tw3jRc+slAnzsK8FJhg70afuxVGj3XBYIDVqy14eQXw0EMNKCoyM3++rFLmKJLAhRCiioocNw6Akxe16U7NW3eh+oRD7hUo1C6nvdNOosODQmKIIgSAMxzjQXxJHq/w09jztHM7RLoJrvhCqzBIz4V9iToefdSAzabdC586VRvjMnNmNDabzI/uCJLAhRCiiooYNAj3gACuxZ8mo1UbKCrCbNHuVxN3DdCjZK0n0KQ9IRLMCYwYuUgy/QBd30L8P0jg6w0j8Sgp5KskeLqv9vEv18Mzz2jd6F9/baZPnwjq1PHh7NlMtmw57/C2VkeSwIUQoooyGI20HD0agJPe9kVM9ieBTg+ntoCxB6gWaqRpsyPmsoZImgFQShz1MLK5ayR5Xe/j7eX/ZncWtGkHPh6wNx68autp2VLHtWsqK1daf70Kl/nRHUMSuHCYnj17OrsKQlQ77Z58EoBTh6Ix+fhiO3oUa0BXbYGTi9pEL66pq/BQ78FKHi3sg9mOcZxh9ufDP/7kRSZv+4b7zhxk7lUY31079hdrFZ5/XrsK/+wzE5MmtcXNzcD69QmcPn0nU1CLP0ISuHCYVq1aObsKQlQ7NVq0oG6nTpjy8znXqi0ApkQ3rTB6B7iEQWkStbKbaGWspT51KaWUBqTgjo6NtdxJ/OxTpnwxj13vFdPQ04aiwHc7oUd/FwIDFaKjbcTHGxkzRvs9//TT/U5obfUiCVw4zLVr15xdBSGqpainnwbgaMoVVBXMG3eg+jaCnEuQax/MdvUgBvwpJp4OaN3txzjEQLQpV1eM7sL6Rk/z77dzWPuzmYH3gckC32xTeOYZbWKXDz808eKLHQGYP/84qakFjm5qtSIJ/BY2bLbRqKWJDZurz+pZ5WnRokXOroIQ1VLzoUPxrFGDtDNnuBLZFoqKMOXfoxUePaPNzJazjZDCBwDw5Gf88CWDTLqSgwKsJouWUxrx/bhkmp9fQKcaVhoVmPhqsY1xE11xd4e1ay2YTP4MHtyU0lIrH3+8z3mNrgYkgd/CynUWpi/LZOU6i7OrIoQQf5jBaCRq6lQAYlQ9AKYNh1CNfnDhMJi1oeWBl8+h4EoeO+hCfQDi2U93fDCh4jEsm1MTGjMyZQ7Hfohn+qxM/Eot/HBYx5NPavfC33mnlNde0yaMmTkzmrS0Ikc3t9qQBH4Lg/sbmDY0kMH9ZdVVIUTldu/UqRjc3Tl36DAZDRujXr2K2arNhU5MGigGdOkrqFmsPVJWg3144sklLvMwWhL+Xsng7wO96f7PbbTOWMmqURcwFlzkgxXw1HOu6PXw448WdLqa3N+tG6ENh/Dqa5ec1eQqTxL4LfTtpePcCVf69pIwCSEqN8/gYNpNmQLAAXdfAErXxqDqveDMbjD1AWzUuHgBBRfy2EQ36gFwht08iDfF2LhQI5WIUB9efesNhg08yfb4B1m+6UFWPLuads0LaeqZwoynY6gX2ozpy7xJvuBL4saDcPUqb/1DJni5m+TSUgghqoku06YRM2cOCdHRXG7WiNoXzmEq6IbRfQfsS4ZuRvSZawjNm8Jln10EsQFfWpFKKj2zz3FpyU6Stx1i+KkrDMguJlUPebUNhBfEsyDLh2nLYP2oK3yePJDtaQ+zcvQU/l/O1wRO2cAac19Saj7JN3u34D65JvsGDKYw3ws13Re1yBVbFringr5QxVpkwy3fij7PhsFgQ1EyMBqTMZmKCQ72pGPHOnTsWBedTnF2SJ1KErhwmH79+jm7CkJUa14hIXR+9VV2vPkm2wpNjFSBnw5gGFYDfUosXO0NIZsIPreTa5G+FOqieaiwJaun7+D4J9NpW2gCoNB+vMvAkqvQwQADCj9n3dinGVJzFh4dS3g4cjkPRyynyMvIxcahfPfacwyd2pBFb/rx6X8G88BbH/POp6+yu2cXsi4Ek54WiuqngwwFDumhRA/ZVjhRDLl+gAcQA2wHrEREBPLBBz0ZODDCGaGsEKRvWDhMs2bNnF0FIaq9zi+/TGBEBJnJyexu0ARKSig+7I1qAXbuhMJ6KMWnaXwhkNS1+Wxv/gb6f+6GQhP1dNDLAMO84dGG0DwUbMA+C5iKf+Lj433ofOgnrOlANJgP6fHIL6VxbDI9us5n2czzdJzwA4v2DOWLMS9g/n/3MHDUNoLDUgmLOoveYILGwEAbeJnBXw9dPDBGlACuQAeCg8cTGlqTM2cyGTToe555Zh0WS/V8UkgSuHCYlJQUZ1dBiGrP4ObGo4sXo3d15djpsxzw9MMad57iE/VQi0tgWzE5CXpWT1rPoQGXKE42UcMAw1xhsJdC2/bQdDg0HWVg6EwY8bEOFy+FM6WwOdgFCiBzt5GSEgMuCVbOx4VhsFkZU3M+419+gx6913Jv4REOJw9iwMJmXL7Ung8H/h1vv1waRJzB5WIq+Oigv4Ww4Iug01HapAbPfuJJgwYK6eneuLuP4u23+2M06pk5M5qRI3/Eaq1+SVwSuHCYH374wdlVEEIAoe3aMWjePFAU9mfmsEI1cCImmeiVrvy4IJ0vJliJ3wEuCnQ1wON6qFcffAapuEQpKF3ApaeF3HAvbE+Fce+G+ihGHfGXzOxsHIRbVinnj9eiFAMNjyWx9XR3XDHT7lwsCaeaEOydTqeOP/L5x0UYpiURnJvBsJ7r8PwilBBLKgHnL4K3O7YR9flYe/qNGT/reXOWJ23b6jh/XmXlyqasW/cEvr5Gli2L469/3ejMkDqFJHAhhKiGWj7+OMNXrsTNz49kk4WNZlifaeJkKthUaKqDJ1whKgA8HwKP7qBrA7oBKqbmBi42qUViizqkeNTmeJf+mL/UxrgculBIQu2a1I1PYd2VHgC0jznIz2e7EuiSTWhsOleSQ3hs7ALubTMf/w4XWLO1L3El3XlzqomApLq8/hcX2vtDcjFs84Pp47Q6Pz9Xx+yFHjRqpOPYMRsLFwaxZs3juLrqmTHjMAsXHndSNJ1DErgQQlRTTQcN4rmEBHp88AHNHnmEiHZt6dK6BZM6tGBg9ybUGNUIzxfDMQwLwzo6nIJ+jUlp1ZoT7boRX+N+Dij92cVjnKYDe8e/zNUB3TEUFPOtVxNsikLXrVtZn9UTb1sB6k4917KD6FJvH8c+b4XVqjB88o+ceS8Li5sBj78lsmrEWZoE7ORb3TXevbcIfxdYnwYB7eCxTpBfDC8u0LFihTvu7jB/vpn09FC+tP/xMHXqepKScpwcVceRUehCCFGNeQQF0fnll393PwPgZX/VvdEOOsibNZ8Z2yMIPr6b1MEDCd2wmk4ZhVDbl4fYjtnzJeAjXh+6iBTzM2S4LeOTj65yki4cHR5LHUMSY6fOpP+kRcyu7cEXLZswOkbhlTg48CTsiYPdcXDgip4PP3Tj2WdLeOaZEuLjI9m48Rw//hjPlClr2bhxFIpS9R8xkytwIYQQd4VP7dp0efVVAHZfvITq44Nu3z4stUcC4LJrHXg9AJZMaqfmYSCYYuUEndHhhhuXHlPJf2YEH45/m5O2ArxqZ9MzGLLN8Okl+ExbGZW/L4DHx7jQqZOe1FSVd94xMXNmf/z93di8+Tw//hjvrBA4lCRw4TCPPPKIs6sghChnHV54Ac8aNbgSE8Pl/oMBKF19FDUgDFLjIVtbTFx39Stqm8YAkMVsuqFN67r8tSY0yi5m4MJ1/DP3Kh+2sKFX4OtkaB0JD9wDmfkwfYXCjBluKAp88YWJvDx33n1XO/bLL2+htLTqr2EhXehlZGVlMXbs2Ou2LiAj404Wpg/63f1zc3PvuOxG26/fVvb9rcruljs9pq+vL1lZWdhs1edxj/KIe2VS3dsP1TMGradOZd9bb7HnxCmG+fnDoQPkD3gWn6wZWLYtwvpYX4yFG3A7dwRDRB1M+hRq51/Ey8OTq4Z0kj94lnZjd7D5UDALe6cxsqYfC6+58dqJUt54tJhep/z4Yq3KmC45DB+u4/vv9bz6ai4zZtTjs88COH06i48+2smUKa2dHYpfv/9NmzaxadOm64t9/8yxFVWVuWl/ERUVpUZHR//PNkWBOwnR7eyfkZFBUFDQHZXdaPv128q+v1XZ3XKnxzx//jwNGza8q3Wo6Moj7pVJdW8/VM8YmAoK+LRuXUpychg6bjR1vl+EoVcfPFrEQVYyPPERWF8CnSdZ937FRcO7GAknl9dZywZCqcWhPp0Y9FkN3p/pxtL3/YjYqmBRIe5B+PsXsOoAvDAI/tLDRuPGBVgscOqUJwkJ5xg06Htq1PDkwoW/4OHh4tRY3Or7VxTliKqqUX/02NKFLhxm1apVzq6CEMIBXL28aPfUUwAcTcsCd3csmzdibazdC2fvT+DXE2yF+KdexIValJJIEwrwxJMrXOXBSe6sGRmPsXUcJ9yzGV8PVOCDc/CPEdphZm8ET38dEye6oKrw3nsmBgxoQlRUKGlphXz99RHnBMBBJIELIYS46+6dOhVFrydx82ZKBz0GgOlIJrj5wLndwCAAlKtfUUPVEnsW39GR9tr2IQfp3WInb6z/B/NI46WGKgqw6BLUDIW+7aCoFL5cB9OmGdHr4bvvzCQnq7zxhrYe+Ucf7cdksjq87Y4iCVwIIcRd51u3LhEDB2KzWDjl5QeAeen3qG3sV+FHj4F7BJguEZhpRIc7BRykNQEYMHCWBLrOmkCbAycwHI4mx6uAwbXArMKsJHjlUe0wX66HkFAdw4cbsFrh889NDBgQQfPmwVy6lMfSpSedEwAHkAQuhBCiXETZu9GPr12LrkMnyM/HnBWqFUYvAf8JAOhT5xPAQAAKWU0rWgJw1COek2++xAt//4LFZPDXcO2js5OgQzOIbADpubB0N/ztb0YA/vtfE4WF8Le/dQTgk08OUFXHekkCF0IIUS7Ce/TAu25dcpKSuNKhCwDmleuh8QNQWgjJrqBzg9ytBBVrCTebNdxHKwBiOErE+MnUSk4ld+sWIgLMtPKBNBOsSIXnB2j/nxnroG1bPV266MnLgwULzIwc2ZLgYA+OHUtlz55kp7S/vEkCFw4zbNgwZ1dBCOFAik5HxPDhAJxMSgZvb6yHDmCt3Vfb4eBSCBwKgHvaLjxohZV8PDhJKLUopph0lyRmP/Ypxf9syXtb8pgapn30qyQYcT/4e0H0OYhOgOeecwVg5kwTRqOeKVPaATBjxmFHNtthJIELh6lb94YTMAohqrCm9gR+evVqbAO1iV3Mx9LBzRsuHABFW/CEtPkEqlp5FiuIQnu66ggxJBb2ZticMI5/V8iIOiqeetidBRdNME6bu4XZm+CRRwyEhCjExdnYs8fKlClR6PUKK1bEk5pa4NiGO4AkcOEw8fHVY3pDIcRvfMPCqNu5M+aiIhJD6gBgXvYDauQQbYdT8WBsAKZL+Od6oOBGAYdphj8GDCRygccetjBrWh79E2cRbyhgRG3to3OTYXJv7eclu6DUojBhgvbc95w5ZurU8WHAgAgsFhtz5x51dNPLnSRw4TDr1693dhWEEE7QavRoAOIOR6PUD0O9fAmrQbvPzeHvIFibUlWfthw/egJQyCbuoTkAgb1iiGtlZFTiXA7F7GRiPe2j81OgYSh0aQ6FJbBsL0yapHWjL19uJjtb/bUb/b//jcFmq1qD2SSBCyGEKFfNhw5FZzCQuG0bpoe157/Ne09BQD1tZraCptqOmSsIsGqX1FmsJhJtKtRjHOfxwBDmvjyWe/79H1r6W2nqpQ1m25gGE+y98HO3QoMGOnr00FNSAkuWmOnZM5x69Xy5cCGH7dsvOLzt5UkSuBBCiHLlERhIw169UK1Wzrl7A2Be9SNqpH1g6/Fd4NUebAV4Z6ViIBgTKdSiAG+8ySKL5/+SxrknR9N67zEOxh1gnH1IzfwUGNoZPN1gbzwkXIEJE7Sr8G+/NaHX65gwIRKAb76pWt3oksCFEEKUuxaPPw5A/K5d6O5pAdnZWIrtfeExyyBAG+ympH+PP/0AyGEdre3PhB8jlj4etZn/wig83p3OqDqgAGuugVkPQ7TFzFj4szaYzdcXoqNtnDxpZdy4SBQFVq6MJzu72KHtLk+SwIUQQpS7iIED0RuNXNy9m5LeWoK27DgEIc2gMBOyQgAd5GwkwKJNhZrDJlrRAoATnKQH3vz49FAiNu/G7dIpugeByQY/XIEnHtL+Pwt/BqNRYcQIbTDb/Plm6tf3o3v3cEpLrXz/fdWZmU0SuHCY0faBLEKI6sfo40Pjvn1BVTmn17q4zWt/Qm2lzZPOsU3g+yCoZtwz4zEShoVMvEgimGCKKCKdi0T5hPLds8M5+ey/GGPvRl90Cbq2gLpBkJQGe+Jg7FgtgS9aZMZiURk3TrufPm/ecYe3vbxIAhcOU7NmTWdXQQjhRM2HapO2xO/YiS6yDeTnYymwT616fBX4aclcyVhapht9w6/d6LGc4GECWPT847Q+uIFHS5Nw18OeLEgugVFdtUMt2gEdO+pp1EhHaqrKtm1WBg9uhre3K4cOXeb06QyHtru8SAIXDhMbG+vsKgghnKjJww+jNxpJ3rOH0h7aaHPLz/ugTmsozoU0X1AMkLMdf7O2KlkuW2lJBABxxNMJNwoP12XM/bvYPG0dg0K0Yy+5DGMe1H7+YS+UmhXGjPntKtzDw4WhQ7XH0hYurBpX4ZLAhcNs2bLF2VUQQjiR0ceHRn362LvRtcVHzGt/Qm1pX1rs+AZtnXCsuGUcwZ0IrOTjQjy1qY0JExc4R811QUz4KJj1Sc2YqLsMwOJL0KwutAmH3EJYFw2jR2sJfMUKMwUFKmPHat3oCxfGVolnwiWBCyGEcJjmQ7QZ2OJ37EDXOlLrRs+voRXGrgZfbTpVMpbihzZnenaZbvQTnOTh/jremelKnYgTdJv3AYGucCofTuTB6G7axxftgPBwHZ076ykqgpUrzdx/f33q1/clJSWPnTuTHNfociIJXAghhMM0GTAAnYsLybt3Y+qpJWjLtr1QJxJK8iDNCxRXyNuJv0mbRS2X7TSnEQoKZzjLpF56il67wvaPm8HChUxwvQbAd5e1BU4UBdZHQ3bBb1fhCxea0ekUxozRZoBbsKDy39KTBC6EEMJh3Hx9adizJ6rNxjlXdwDM61ajtrJ3ox9bB/59ABVj5n48aImNYhSOUZ96WLCQwBnylvmRUSuY2JEDee6nTwHtPnhIADzUCkwWWL4Xhg1zwcUFtm2zcuWKjTFjtG705cvjKCoyOyMEd40kcCGEEA71y2j00zt3at3oeXn/243u84j2c8ZS/OgDQDbraWVfJ/w4sQxx9wdg+iujqLPoawIuZJJcDHuz/rcbPSBA4eGHDdhs2tSqTZoE0r59bQoKTKxaddpRTS4XksCFw0yYMMHZVRBCVAARgwahMxhI2rEDcy/7pC5b90DdNlCSD6luoHODvD34l2r3vnPZSXPqo0PHeRJ5fyLUxpUT9YL44f43iVqajtsJG4svwaMdwc0Vdp2Ci2n8Ohp9wQLtivuXwWzz51fu0eiSwIXD+Pv7O7sKQogKwN3fn/AePbS50Y32bvSyo9GPrQH//gC4ZuzCk7aolGLmAI1oiA0bJ4mjP9q/KfNCRjHlW18aHi/ghytgNMIg7Sk0vtsJ/foZCAhQiI21cfy4lREjWuDqqmfr1kQuX85zePvvFkngwmGio6OdXQUhRAXRfJi2kMnpspO65AZrhSdWg88vo9GX4I+WzLNZR2t7N3ossTxsT+BXHs1j9rBchl/+hmwzbEiDMd20jy/4GVxdFUaMMADa1KoBAe4MHBiBzaaycGHlHcwmCVw4zM6dO51dBSFEBdH0kUfQubhwYft2TL21BG3e9DOEtYfSQrisgt4bCqLxK24EGMhjH42pgQsuXCSZQEpoijvmXrmcqOfJCwkf0PbcERZegl5toIYvnL4E0efgiSe06VsXLzZjNpedWvUYqlo5nwmXBC6EEMLh3P39tSVGbTYS9No9asu61agt7N3oR5ZDgDaYzSV9PT50AqwUsZ1maOuHHyeWgQQAoJ9oxvyvd5kzczLrL1vIs8FI+9Sq87fBvffqaNpUR1qaysaNFnr3bkRIiBdnzmSyb1+KQ9t+t0gCF0II4RS/LDF6atNm9B06QXExlmueoOjg5HrwHKDtmL4Yf/VhALJYTSTa1fNRjtMPP3SAz4Bc1CdHofP359lVnzLm699WKPtulza16vjx2h8Kc+eaMRh0PPGEdpzKuk64JHAhhBBO0XTQIFw8PEjZu5dC+9zopp/WQtMeYDXDuWvgEgIlCfgVeKPDkyJOUAcFL7zIReiQ1wAAEM9JREFUIIMSrtEZH1Q9rFNyuPTZHKatmM7Zy+eIDIfIBtqELj8d1FYo0+th7VoLqak2Jk5sA8DSpafIzS1xZij+EEngQgghnMLVy4umj2jd5HEFxWAwYN26GVvjgdoOBxdD8CgAdGnf429/JjyHNUTaB7PFcJRH7d3oK8ikR/sGfDTyTZb+MJxjqcVM6qUd6uvNEBKi4+GHDVgsMG+emcaNA+nWLYyiIjOLF59wYMvvDkngwmGmTJni7CoIISqYVmPHAhC7bBn63v3AZsN8KhfcvOHCAVDv13ZMX0KAVZt6NYtVtP11NPoJOuOGPwbOUEKCvojiqc+RUKsxRU89w6gHVNxdYdtxSLgCkydr3ejTp5uw2VSeekqbrnXWrOibDmZ7661yDMCfIAlcOIyXl5ezqyCEqGDCe/TAu3Ztss+f52q7+wAwL1qI2m64tkPMLvC6F6y5eGYmYCQcCxkYiaMedSmllHjiGGS/Cv+BDFpcUHmrzrecOuaJ28yPGWH/G2DWBujd20BYmEJOjsqGDRYGD25GzZqenDyZxq5dF4H/m7DfftsRkbhzksCFw+zbt8/ZVRBCVDA6vZ7IceMAiD1xCiUkBNuZ01jdorQdDsyHQK1cSZ1NENo0rBks5T7uBeAghxhqT+DryWb/djPvvpXDZ/f9m9IvZ/JWxucAzN0KxSaFqVO1R8o+/9yEq6ueKVO0q/DPPjsIVNyEfT1J4MJh9u/f7+wqCCEqoLaTJoGiEPfjj5gf0668Tat3QN22UJgJFw2g94X8/QQUhKFgJJ+9NMYTDzy4ylUUrtEZb0pQCehfyLRJgZxv7cUj722n7pLP+fDIS+QXWJm3DSZO1BL45s1WTpyw8vTT9+LqqmfVqtMkJGQ6MRJ3RhK4EEIIp/ILC6NJ//5YTSZO6FxBr8eycjm2e7THzNj5FdQYD4Dhyn8JQBvkls13v16F72E/Y9BmcovudYXz0QZq3qdjh3cY7/79CAeDh/Hh4TdYOe8MPr4K8Nu98JAQL+65pxWqCh98sNexjf8TJIELIYRwuvZ/+QsA0QsXwoBHwGLBtO8yeNeAlKNQEAXoIeN7apRoD3hnsoooGmLAwGlOE0EJjXEjDTO+o7N4oaF27BW73Hl8aV02NX6Spas7823vr2nR0oSnj44lS8ycPWvl6NHO6HQK8+YdB3KcE4Q7JAlcCCGE0zXo3p2QyEgKUlM5HdYYANM3X2NrN1HbYdt/IXgEqBbcLi3Bl4dQKaWQpbSlDSoqO9jFZGoC0PDDazxRz0YNV4hrYmDaxEAiHwmjddfj/FTUhbeXFfFw0wTa2GL4xz9KgUBGjmyJxWIDdvxarw2bbTRqaWLDZptjA3IbqnQCVxSlh6IokxVFma4oymRn10cIIcSNKYpCl9deA2DPwoWoPXpBYSGmQ3ng7gdnd0Dhg4AO0r6hVrE2f3oGS+hME/ToieUErSimIW4U+Zn4yZDBq42hpKWOcyNceft5Hd6NarPNGMHrY3wZ0uwCKxjLs98/xFC+563XO+HlE0Gjls35clYGACvXWZi+LJOV6yxOiszNVdkErihKOICqqnNUVZ0GTFcUpa2Tq1WtPffcc86ughCiAmv+2GOEtGlD/uXLxNQJB8D09RysLexX4Ws+hcCxoFpwvzALP/qhYiKfWXRAWz90Axt4kVoAzCSVYWFmwj2AejArGaY/ASXeOk4HetDl00eZ/Uocn/AMz/IF4T2j6NW6OdOXeTN3cQEWi43B/Q1MGxrI4P6Gm9bbWc+JV9kEDrQFys4cshWIclJdBODq6ursKgghKjBFp6P3p58CsPfbb8kdOBjMZkqWHEINCoerp+BskLZKWfYa6qY3RIc7OWwmCiteeHGRZIzE8wA+5GPlPX0K/2mhTdDyxmmIuAf6tgNcYPKX8Oobnhyu/ygTW27hpfZbGO+1jbVDT9Hy2Nd89O/t9O2l49wJV/r2unm6dNZjZ1U2gauquhx4ssymcCDRSdURwI4dO5xdBSFEBRfWtSuREyZgLS1lTWw8puBgrHt2U5rXWdth06dg067IDedfpW7RCABS+ScD0DpZN7KZySh4oWM7uWSEpMHPUGSFodHw4WTADGsOwYyNCu0fMDB9WS5rEurh+tfn2HzCjUcKf2TQW4+zb9YaZ4ThtlT4BK4oip/9Pvaym5RPVhRliP31StkyVVVz7PuEA1mqqm51QJV/16ZNm+647Ebbr99W9v31ZQsWLLiTKt6WW7XjRo4cOXLX61DRlUfcK5Pq3n6QGPyR9vf5z38IatqU9NOn+cmnBsWKgmnOQkpND6JarbDsGyi6H9WSh7LlUwqWBHL2kwskvDyBqL/spdm/tpP50ZP8d9YsRm9YSOaOd1kY9Q9e8/wvHSxzWJH6DW+OWkWfrmvZmLydgPZXefsZV0yBqYydcZXc/s1ZP+9r3g8bQOPnRzGgwTwOZUH8VdgbBztPQnyKyuUrNo4csQIWEhNt5OeXcuZMBjExVzl/PguLxVau379SkRcyt9+zDre//buqqu2uK5+MlpiX29+HA9NUVZ1y3X6zr992I1FRUWp0dPR1dYA7CdHt7D927Nibfqk3K7vR9uu3lX1/fVlkZCTHjh273Wbcllu140Y+/vhj/va3v93VOlR05RH3yqS6tx8kBn+0/dkXLjDvgQfIu3QJT18f2hbmUUsHurr+5Llnc6kAkjJ15OXfeHR4mA66GCDYBVwagjESdH5AY+AesLjryAzx51qdQEoMRpKpx2maUognadTgDBFkEUjTnWeYMfxFPu33Ah8NeBkuAyeAq4BJhVQznDdBgQ0oBs4BMUAWHh4uuLhcZMWKV3nooQb/p46KohxRVfUP39q9+V35CkBV1Rgg5haDz6aUTeqqqiYqitKj7A72q/Jp5VhNIYQQd5l/gwZM2LePFSNHkrxnD7t/KTifXWYvG25AiA78FNAbINcASUWQZINkEzygQuQZMCeDx/1gMIEtScHQxUZNayZ+6XlcbBKKi28iddUUYmiDTlEJ4RrnrOGc6tqCUdFz+br7VPRuZqb3eQ3qAWdV2K9APVeo4wLnC+CMO6gtgXvw9j5Jfv4OIJTu3RfQv39j5swZQGio912LUYVO4LeiKIofcKPEnqMoSg9VVbfak/nyMl3pPSpKN7oQQohb861bl3E7d3J27Vrili8nMy4ONSMdH3MpIaqVuu42QvxtGGra0DdQ0IXZUGvYyPPW8/MPJmK/L2GHGTIbutH9UgkF2xWKHgkglEzM2wyk9KhDeFASDU5eZmOt7tQJT6YDh1h/Poi8et40cknE/VwJ+63tmbz1SxZ2mUhepsKsR1+CJi7ofHII/TmHSwX1obE3Ub1stCgqZcFcM/n5rWjXrhUpKS9QUtKYdesSaNNmNqtXj6B9+zp3JT4V/h74LYRz4+lysoBw+1X7MmCLoijZiqKo9rI70qrVyj9XSyGEEH+YotMRMXAggxcsYFJ0NE8mXWT45VS6Xkkn/HwmHtHZuK7LRT8jB+WlPHRjC/AbnMvgJcU8ungxOoOBE+dLOH7/Q+isKl6rcjD798XFYiF861VU6xBcsNAnZRsNsgeCqqNfw/VM0Pvio/pQu9EVBrtG07lWBJO2zOSt3Z/wz70LaOZhwxbih8vEMDilUCsAoi/oiPd1Z+sOD+rVUzhyBPLynmPv3qk89FAD0tIK6d59AXv2JN+d2FTke+C/sCfjr8t2l9uvrmerqtrwun2XAYdVVf3gD/x/8vm/f9QUAncSbV8g90/sc7OyG22/flvZ99eXBQEZv1OvO3U7bb1eYyDhLtejIiuPuFcm1b39IDGQ9mvtDwL7ZO2/samq+of71CttF3p5+DOBFEIIIRypMnehA/YFYP+Xn8NrIYQQQjhYZU7g0dw4WQegjeEXQgghqqxKm8DtI8sT7aPRy/KTkeZCCCGqusqSwG/UVQ4wHfh1lTH7YDdJ3kIIIaq8Cj0K3T6z2hCgJ9AD+AA4r6rqnDL7TEab49wPCP8jo8+FEEKIyqZCJ3BRfuy3HoYBDe3LrQohhKhEKksXurj7fpl/V0btCyFEJSQJvJqyD/ST5VWFEKKSkgQuhBBCVEIyE1slVuY+dk9VVYfeoHwyv83/LgP87hL7NL5TgPfQ5uMfAuTcYHDlTWNfmb6bP3ueVfZY3Kr91eVcsNcRoCHabbdpvywSVaa8yp4DcOsYOO08UFVVXpXwhbYS2xD768gNyicDQ8q8D0ebO77sPj2u3yav24r9EOAIoALZwPQ7if3tfDcV5fVnz7PKHovbaH+VPxeAyTdo8/nqcg7cZgycch44PTDy+tMnVtub/MNyo23nr3svCfyPxXzI75TfMva3891UtNcfPc+qSixu0f4qfS7YE8n0G2zP/qXtVf0cuM0YOOU8kHvgVdDvrZXu6PpUJ78X+6r03fzZtlalWNxIFWr/5BtsywICqtE5cNMY/N4HyzMGcg+8arrlWunwP/dswhVFmayWuVcjfp89fn5ocW6r/na/6vdin/U75ZXJn21rlYhFVT4XVFVNBPxvUBSOth5FlT8HbiMGgHPOA0ngVVMAvw2GKCsH+3PfqvYYmUw7+8fEwK+/2CiKkqUoyhZVVXvy+7H/3e+mEvmzba0Ksah254J9sNVWVVVj7Emr2p0DZWNg3+SU80ASuBB36Jdf0jLvYxRFibJP/Suqkep2LtjbNUVV1XbOrouz3CgGzjoP5B541SVrpTtWItqgQPj92Fel7+bPtrUqxeIXVflcmA50v25bdTsHbhSDGyn380ASeNUka6WXE0VRwhVFyb7FLr8X+6r03fzZtlbqWFS3c0FRlOlc9/w31ewcuFEMnHkeSAKvglRZK728vXeDbeFo98RuGfuq9N382bZWkVhUi3PBfs93dtmuYkVRelSnc+BmMbD/6JTzQBJ45SdrpTuQ/Zf3f0aMKooyBPihzC/278W+Mn43f/Q8qyqx+D/try7ngj1JRZcZoOV33eNNVf4cuFUMnHkeyHKilZSsle5cZaZV/GVU/42mRbxp7CvLd3M3zrPKHIs7aD9UwXPB3v7zNyn2V3+bSrSqnwO3GwNw4HkgCVwIIYSohKQLXQghhKiEJIELIYQQlZAkcCGEEKISkgQuhBBCVEKSwIUQQohKSBK4EEIIUQnJYiZCCIezzzo1HW0CjMOqqi53cpWEqHTkClwI4QzbgOmqqk4D7nV2ZYSojCSBCyEcSlGUZWhzRP8yzWRFXnlKiApLutCFEA5TZmpSf/t7P66bR1oIcXvkClwI4UjTgTlllmMcBix1Yn2EqLQkgQshHKLM1feyMpsbqqpa4dZ9FqIykAQuhHCUKUDOL2sc/7K+snOrJETlJQlcCOEok7GvcWy/9+1XZiCbEOIOySA2IUS5s3ef+wFb7Jumq6o6xYlVEqLSkytwIYQj9LD/d6uiKLOBac6sjBBVgSRwIYQjtEN7XGwKMLvMKHQhxB8kCVwI4QgBaF3oS2XUuRB3hyRwIYQjhAPLJXkLcfdIAhdClCtFUV6x/5jl1IoIUcVIAhdClBtFUXoAy9FmWwu4SbkQ4g9QVFV1dh2EEFWc/bnvI6qqNiyzbTba42TyLLgQf4AkcCGEQyiKMgToyW+Ll7wno9GF+OMkgQshhBCVkNwDF0IIISohSeBCCCFEJSQJXAghhKiEJIELIYQQlZAkcCGEEKISkgQuhBBCVEKSwIUQQohKSBK4EEIIUQlJAhdCCCEqIUngQgghRCUkCVwIIYSohP4/1/B+LYfzT4YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = len(param_values)\n",
    "\n",
    "values = range(N)\n",
    "\n",
    "\n",
    "\n",
    "label_size = 18\n",
    "title_size = 25\n",
    "legend_size = 18\n",
    "handle_length = 1.5\n",
    "xscale = 'log'\n",
    "yscale = 'log'\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "fig, (ax1) = plt.subplots(1,1,figsize=(7,5))\n",
    "ax = ax1\n",
    "\n",
    "\n",
    "\n",
    "#verr = err_p_Dl_low_ell\n",
    "#v = Dl_low_ell\n",
    "#verr2 =err_m_Dl_low_ell\n",
    "#verr2[verr>=v] = v[verr>=v]*.999999\n",
    "\n",
    "ax.errorbar(low_ell, Dl_low_ell, \n",
    "            yerr=[verr2,verr],elinewidth=1,capthick=1,capsize=0,\n",
    "             label=r'$\\mathrm{Planck\\,\\,data}$',color = 'b',ls='None',\n",
    "            fmt = 'o',\n",
    "            markerfacecolor='lightblue',markersize=2,linewidth='0.2',markeredgewidth=0.5)\n",
    "\n",
    "\n",
    "ax.set_xscale('log')\n",
    "ax.set_xlim((1.7, 30.))\n",
    "\n",
    "if param == 'A_s': \n",
    "    ax.set_ylim(100,3e4)\n",
    "else: \n",
    "    ax.set_ylim(100,7e3)\n",
    "\n",
    "\n",
    "divider = make_axes_locatable(ax)\n",
    "axLin = divider.append_axes(\"right\",size=4.5, pad=0, sharey=ax)\n",
    "\n",
    "\n",
    "axLin.errorbar(high_ell, Dl_high_ell, \n",
    "            yerr=err_Dl_high_ell,elinewidth=1,capthick=1,capsize=0,\n",
    "            label=r'$\\mathrm{Planck\\,\\,data}$',color = 'b',ls='None',\n",
    "            fmt = 'o',\n",
    "            markerfacecolor='lightblue',markersize=2,linewidth='0.2',markeredgewidth=0.5)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "axLin.set_xlim((30., 2550.))\n",
    "\n",
    "\n",
    "ax.set_yscale(yscale)\n",
    "ax.set_ylabel('$D_{\\ell}^\\mathrm{TT}\\,\\,\\,\\,[\\mu\\mathrm{K}^2]$',size=title_size)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "if yscale == 'log':\n",
    "    axLin.yaxis.set_label_coords(-.5,0.5) \n",
    "if yscale == 'linear':\n",
    "    axLin.yaxis.set_label_coords(-.5,0.5) \n",
    "\n",
    "\n",
    "\n",
    "ax.tick_params(axis = 'x',which='both',length=5,direction='in', pad=10)\n",
    "axLin.tick_params(axis = 'x',which='both',length=5,direction='in', pad=10)\n",
    "\n",
    "ax.plot(Plc_BF[:,0][Plc_BF[:,0]<30],Plc_BF[:,1][Plc_BF[:,0]<30],color = 'r',ls='-',lw=1)\n",
    "axLin.plot(Plc_BF[:,0][Plc_BF[:,0]>=30],Plc_BF[:,1][Plc_BF[:,0]>=30],color = 'r',ls='-',lw=1)\n",
    "\n",
    "\n",
    "\n",
    "cosmiv_variance_low_ell_low = np.asarray(Plc_BF[:,1][Plc_BF[:,0]<30])-np.sqrt(2./(2.*np.asarray(Plc_BF[:,0][Plc_BF[:,0]<30])+1.))*np.asarray(Plc_BF[:,1][Plc_BF[:,0]<30])\n",
    "cosmiv_variance_low_ell_high = np.asarray(Plc_BF[:,1][Plc_BF[:,0]<30])+np.sqrt(2./(2.*np.asarray(Plc_BF[:,0][Plc_BF[:,0]<30])+1.))*np.asarray(Plc_BF[:,1][Plc_BF[:,0]<30])\n",
    "\n",
    "ax.fill_between(Plc_BF[:,0][Plc_BF[:,0]<30],cosmiv_variance_low_ell_low,cosmiv_variance_low_ell_high,alpha=0.1)\n",
    "\n",
    "\n",
    "#unbined cosmic variance\n",
    "#cosmiv_variance_high_ell_low = np.asarray(Plc_BF[:,1][Plc_BF[:,0]>=30])-np.sqrt(2./(2.*np.asarray(Plc_BF[:,0][Plc_BF[:,0]>=30])+1.))*np.asarray(Plc_BF[:,1][Plc_BF[:,0]>=30])\n",
    "#cosmiv_variance_high_ell_high = np.asarray(Plc_BF[:,1][Plc_BF[:,0]>=30])+np.sqrt(2./(2.*np.asarray(Plc_BF[:,0][Plc_BF[:,0]>=30])+1.))*np.asarray(Plc_BF[:,1][Plc_BF[:,0]>=30])\n",
    "#axLin.fill_between(Plc_BF[:,0][Plc_BF[:,0]>=30],cosmiv_variance_high_ell_low,cosmiv_variance_high_ell_high,alpha=0.1)\n",
    "\n",
    "xs = []\n",
    "ys = []\n",
    "\n",
    "for idx in values:\n",
    "    p3 = \"%.4e\" % param_values[idx]\n",
    "    if param == '100*theta_s':\n",
    "        p3 = \"%.4e\" % param_values[idx]\n",
    "    new_root = 'output/bard_'+str_param+'_'+str(p3[0])+'d'+str(p3[2:])+'_'+'cls.dat'\n",
    "    DATA = np.loadtxt(new_root)\n",
    "    xs.append(DATA[:,0][DATA[:,0]<30])\n",
    "    ys.append(DATA[:,1][DATA[:,0]<30])\n",
    "ys = np.asarray(ys)*(2.7255e6)**2\n",
    "\n",
    "if param == 'A_s':\n",
    "    multiline(xs, ys, np.log(1e10*param_values),ax=ax, cmap='jet', lw=2)\n",
    "else:\n",
    "    multiline(xs, ys, param_values,ax=ax, cmap='jet', lw=2)\n",
    "\n",
    "\n",
    "\n",
    "xs = []\n",
    "ys = []\n",
    "\n",
    "for idx in values:\n",
    "    p3 = \"%.4e\" % param_values[idx]\n",
    "    if param == '100*theta_s':\n",
    "        p3 = \"%.4e\" % param_values[idx]\n",
    "    new_root = 'output/bard_'+str_param+'_'+str(p3[0])+'d'+str(p3[2:])+'_'+'cls.dat'\n",
    "    DATA = np.loadtxt(new_root)\n",
    "    xs.append(DATA[:,0][DATA[:,0]>30])\n",
    "    ys.append(DATA[:,1][DATA[:,0]>30])\n",
    "ys = np.asarray(ys)*(2.7255e6)**2\n",
    "\n",
    "\n",
    "if param == 'A_s':\n",
    "    lc = multiline(xs, ys, np.log(1e10*param_values),ax=axLin, cmap='jet', lw=2)\n",
    "else:\n",
    "    lc = multiline(xs, ys, param_values,ax=axLin, cmap='jet', lw=2)\n",
    "\n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "axins = inset_axes(axLin,\n",
    "                   width=\"5%\",  # width = 5% of parent_bbox width\n",
    "                   height=\"35%\",  # height : 50%\n",
    "                   loc='lower left',\n",
    "                   bbox_to_anchor=(0.8, 0.5, 1, 1),\n",
    "                   bbox_transform=axLin.transAxes,\n",
    "                   borderpad=0,\n",
    "                   )\n",
    "axcb = fig.colorbar(lc,cax=axins)\n",
    "\n",
    "\n",
    "if param == 'A_s':\n",
    "    axcb.set_label(col_label,rotation=0,labelpad = -10,y=1.2,size=15)\n",
    "else:\n",
    "    axcb.set_label(col_label,rotation=0,labelpad = -45,y=1.2,size=15)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "axLin.set_xlabel('$\\ell$',size=title_size)\n",
    "axLin.xaxis.set_label_coords(.4,-0.1) \n",
    "axLin.spines['left'].set_visible(False)\n",
    "ax.spines['right'].set_visible(False)\n",
    "\n",
    "ax.xaxis.set_ticks_position('both')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "ax.tick_params(axis = 'y',which='both',length=5,direction='in', pad=5)\n",
    "\n",
    "axLin.yaxis.set_ticks_position('right')\n",
    "axLin.tick_params(axis = 'y',which='both',length=5,direction='in', pad=5)\n",
    "\n",
    "\n",
    "\n",
    "axLin.xaxis.set_ticks_position('both')\n",
    "\n",
    "\n",
    "ax.axvline(30,ls='--',alpha=0.6,c='k')\n",
    "\n",
    "plt.setp(ax.get_yticklabels(), fontsize=label_size)\n",
    "plt.setp(ax.get_xticklabels(), fontsize=label_size)\n",
    "plt.setp(axLin.get_xticklabels(), fontsize=label_size)\n",
    "plt.setp(axLin.get_yticklabels(), visible=False)\n",
    "\n",
    "ax.grid( b=True, which='both', alpha=0.3, linestyle='-')\n",
    "axLin.grid( b=True, which='both', alpha=0.3, linestyle='-')\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "\n",
    "#figure_name = 'bard_varying_' + param + '_cl_TT_day-2.pdf'\n",
    "#plt.savefig(path_to_notebook + figure_name)  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}