{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Brían Ó Fearraigh, Rodrigo G. Ruiz\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual, Output\n",
    "import ipywidgets as widgets\n",
    "import IPython\n",
    "from IPython.display import clear_output, display\n",
    "from array import array\n",
    "from ctypes import string_at\n",
    "\n",
    "\n",
    "import matplotlib.mlab as mlab\n",
    "from matplotlib.colors import LogNorm\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lateral Spread of Muon Bundles\n",
    "From APP 25 (2006) 1-13, the muon radial distance $R$ from the shower axis is taken into account in order to correctly parameterize the energy of the muons in a bundle.\n",
    "\n",
    "The muon lateral distrubution in a plane perpendicular to the shower axis can be described as\n",
    "\n",
    "$ \\cfrac{dN}{dR} = C \\cfrac{R}{(R + R_{0})^{\\alpha}} $.\n",
    "\n",
    "The average value of the of the radial distribtution is given by: \n",
    "\n",
    "$\\langle R \\rangle = 2 R_{0} / (\\alpha - 3 ) $.\n",
    "\n",
    "$C$ represents a normalization given by\n",
    "\n",
    "$C = (\\alpha - 1 )(\\alpha - 2 ) \\cdot R_{0}^{\\alpha-2} $.\n",
    "\n",
    "## The parameter $\\langle R \\rangle$\n",
    "\n",
    "In MUPAGE, for a given depth $h$, multiplicity $M$ and zenith angle $\\theta$, the average radial distance $\\langle R \\rangle$ is parameterised as\n",
    "\n",
    "$\\langle R \\rangle = \\rho(h, \\theta, M) = \\rho_{0}(M) \\cdot h^{\\rho_{1}} \\cdot F(\\theta)$\n",
    "\n",
    "where \n",
    "\n",
    "$\\rho_{0}(M) = \\rho_{0a} \\cdot M + \\rho_{0b} $\n",
    "\n",
    "and\n",
    "\n",
    "$F(\\theta) = \\cfrac{1}{e^{\\theta - \\theta_{0} \\cdot f} +1 } $.\n",
    "\n",
    "## The parameter $\\alpha$\n",
    "\n",
    "This parameter is given as\n",
    "\n",
    "$\\alpha = \\alpha(h,M) =  \\alpha_{0}(M) \\cdot e^{\\alpha_{1}(M) \\cdot h} $,\n",
    "\n",
    "where\n",
    "\n",
    "$\\alpha_{0}(M) = \\alpha_{0a} \\cdot M + \\alpha_{0b} $\n",
    "\n",
    "and\n",
    "\n",
    "$\\alpha_{1}(M) = \\alpha_{1a} \\cdot M + \\alpha_{1b} $.\n",
    "\n",
    "The lateral spread is plotted below using different values of the above parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define functions and parameters #########################\n",
    "\n",
    "theta = np.arange(0, np.pi/2, 0.001)\n",
    "h = 2.785\n",
    "multiplicity=2\n",
    "\n",
    "rho_0_a = -1.786\n",
    "rho_0_b = 28.26\n",
    "rho_1 = -1.06\n",
    "theta_0 = 1.3\n",
    "f = 10.4\n",
    "\n",
    "alpha_0_a = -0.448\n",
    "alpha_0_b = 4.969\n",
    "alpha_1_a = 0.0194\n",
    "alpha_1_b = 0.276\n",
    "\n",
    "def alpha_0(M, alpha_0_a, alpha_0_b):\n",
    "    alpha_0 = alpha_0_a * M + alpha_0_b\n",
    "    return alpha_0\n",
    "\n",
    "def alpha_1(M, alpha_1_a, alpha_1_b):\n",
    "    alpha_1 = alpha_1_a * M + alpha_1_b\n",
    "    return alpha_1\n",
    "\n",
    "def alpha(h , M, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b):\n",
    "    alpha = float(alpha_0(M, alpha_0_a, alpha_0_b) * np.exp(alpha_1(M, alpha_1_a, alpha_1_b) * h))\n",
    "    return alpha\n",
    "\n",
    "def F(theta, f, theta_0):\n",
    "    F = np.exp(f * (theta - theta_0)) + 1\n",
    "    return 1./F\n",
    "\n",
    "def rho_0(M, rho_0_a, rho_0_b):\n",
    "    rho_0 = rho_0_a * M + rho_0_b\n",
    "    return rho_0\n",
    "\n",
    "def R(M, theta, h, rho_0_a, rho_0_b, rho_1, f, theta_0):\n",
    "    R = rho_0(M, rho_0_a, rho_0_b) * F(theta, f, theta_0) * np.power(h , rho_1)\n",
    "    return R\n",
    "\n",
    "def R_0(M , theta , h, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0):\n",
    "    R_0 = 0.5 * R(M, theta, h, rho_0_a, rho_0_b, rho_1, f, theta_0) * (alpha(h , M, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b) - 3)\n",
    "    return R_0\n",
    "\n",
    "def C(M , theta , h, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b):\n",
    "    C = (alpha(h , M, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b) - 1) * (alpha(h , M, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b) - 2) * np.power(R_0(M, theta, h, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0) , (alpha(h , M, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b) - 2))  \n",
    "    return C\n",
    "\n",
    "def dN_dR(M , theta , h, R, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ):\n",
    "    dN_dR = C(M , theta , h, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b) * R /np.power((R + R_0(M , theta , h, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0)) , alpha(h , M, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b)) \n",
    "    return dN_dR\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzQAAAIpCAYAAAB0TX2GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmczWX/x/HXNWMPYw2l8CNkUFFRyJKolKVIGLLXjZTcJW26K3dISqvuIsowVJJkKWWnFJWsyc7Y1ymMMXP9/vjOnHOGWc6ZOTNnzsz7+XjM43yv73J9PyP3bT5zXdfnMtZaREREREREglFIoAMQERERERHJKCU0IiIiIiIStJTQiIiIiIhI0FJCIyIiIiIiQUsJjYiIiIiIBC0lNCIiIiIiErSU0IiIiIiISNBSQiMiIiIiIkFLCY2IiIiIiAQtJTQiIiIiIhK0lNCIiIiIiEjQUkIjIiIiIiJBSwmNiIiIiIgELSU0IiIiIiIStII6oTHGFDbGvGSM+dMYc84YE22MmWSMudLHfpoaY0YYY74xxhwxxlhjzK407s9vjGlljHnHGLPBGHPGGHPWGLPZGDPWGFM209+ciIiIiIiky1hrAx1DhhhjCgGLgYbAAWA5UBm4GTgCNLTW7vCyr9+A6y46vdtaWzmV+1sC3yU2dwHrgPzALUAZ4CDQzFq71etvSEREREREfBbMIzTP4SQzq4Hq1trO1toGwFCgLDDJh76+TeyvNRDuxf0JwEyggbW2irX2fmttW6AasBAoD3zsw/tFRERERCQDgnKExhhTADgMhAH1rLW/XnT9d6AucKO1dq2PfZfHGfFJdYQmneevAPYnNitba3f72oeIiIiIiHgnWEdoGuEkM9svTmYSfZ74eW/2heSw1kbjTHkDuCK73y8iIiIikpcEa0KTtN5lXSrXk87XzYZYkjHGlABKJjYPZvf7RURERETyknyBDiCDrk783JfK9aTzlbIhlosNxPlz/cNau9ObB4wxG1O5VB04C+z1U2wiIiIiIjnBVcAZa235zHYUrAlN0cTPM6lc/yfxs1g2xOJijLkBp7gAwDA/dBlSsGDBYlWrVq3lh75ERERERHKE7du3Exsb65e+gjWhyXGMMeWAWUAh4E1r7Xxvn7XWplhZzRizsWrVqrU2bkxtAEdEREREJPiEh4ezadMmv8xCCtY1NH8nfhZJ5fpliZ8x2RALxphiwDycfXA+wykdLSIiIiIiWSxYE5o9iZ8VU7medD7LSyYnbvA5B6iHs59NhLU2IavfKyIiIiIiwZvQ/J74WS+V60nn12dlEMaYfMAMoBmwCrjPWns+K98pIiIiIiJuwZrQrAROAVWNMdencL1j4ufXWRWAMcYAHwNtgd+ANtbaf9J+SkRERERE/CkoE5rEUZB3EpvvGmOS1sxgjHkCZ/+ZpdbatR7nBxljthhjXvVTGG8CEcAWoJW19qSf+hURERERES8Fc5WzV4CWwK3ANmPMcpx9ZxoAR4DeF91fBqgBVLi4I2NMX6BvYjN/4mcFY8yPHrcNsNauS7y/HTA48fxe4DVnwOYSo6y1W3z8vkRERERExEtBm9BYa88ZY5oDw4GuQHvgODAZeN5am9qmmympiJMIeSpw0bniHsclPY7vSKPfyTgjOCIiIiIikgWMtTbQMUgqjDEba9WqpX1oREREJEey1qKfJfM2YwypzFRKU+I+NJtS24/RF0E7QiMiIiIi2S8+Pp5jx44RExPD+fMq7ioQGhpKkSJFKF68OMWKFctQgpMZSmhERERExCvx8fHs2bOHc+fOBToUyUHi4+OJiYkhJiaGEiVKUK5cOUJCsq/2mBIaEREREfHKsWPHOHfuHKGhoZQrV47LLrssW39wlZzHWktsbCwxMTEcP36ckydPUqhQIUqWLJn+w36ihEZEREREvBITEwNAuXLlCAsLC3A0klMUKVKEIkWKkC9fPg4fPsyJEyeyNaFRSi0iIiIi6bLWutbMXHbZZencLXlR8eJOUeDY2NhsLRahhEZERERE0uX5A6qmmUlKQkNDXcdKaERERERERLyghEZERERERIKWEhoREREREQlaSmhE/Cw+Hlavhpdfhp494fbboUYNqFoVbrgBmjWDzp3htddg6VL4++9ARywiIiIZZYxxfa1evTrV+2bOnOm6r3LlytkXYDoOHTrExIkT6dChAxUrVqRAgQKUKFGCpk2bMmXKlGxdC5NRKtss4icrV8L//gfz5sHRo+nfP3Om8xkaCi1bQrdu0L49FCuWtXGKiIhI1oiMjOSWW25J8drUqVOzORrvDB06lMjISPLly8eNN95I48aN2b9/PytWrGDZsmXMnTuXqKioZAv+cxqN0Ihk0vLlTkLSuDF88ol3yYyn+HhYuBB69IBy5aBfP9i+PWtiFREREf8LDQ2lTp06zJgxgwsXLlxy/dixYyxYsIB69eoFILq0lS5dmpEjRxIdHc3q1auJiopi+fLl/PjjjxQvXpzPP/+ciRMnBjrMNCmhEcmggwfh3nvhttvg+++TX8uXz5la9vTT8N57MGcOfPstfPYZfPSRc75Fi0tHY86eda5Xrw4REbB5c7Z9OyIiIpIJ3bp14+jRoyxcuPCSazNmzCAuLo6IiIgARJa28ePH88wzz1C2bNlk52+66SaefvppAKZPnx6I0LymKWciGTB3LvTuDUeOJD/fqBE8+ii0bg0lSqTfT3w8rFgB06Y5yc6JE875hASIjISoKHj8cRgxQlPRREREcrKuXbsyfPhwpk6dSps2bZJdmzp1KkWLFqVdu3Y88cQTAYrQd9dddx0A0dHRAY4kbRqhEfHBhQsweLAzMuOZzDRuDIsWOdPPOnf2LpkBZ/1M06bwwQdw4ABMmgTVqrmvx8fD669DzZowYwYEwbo8ERGRPOmqq67itttuY86cOfztUfFnx44drF69mg4dOlCkSJEARui7HTt2AFC+fPkAR5I2JTQiXoqNhU6d4O233ecKFXKmlC1b5lQzMybj/RcsCL16OdPMIiOhShX3tehoePBB6NLFPYojIiIiOUtERARnzpxh1qxZrnORkZGuaxmxZMmSZJXUvPnyRxW1uLg43nvvPQDatWuX6f6ykqaciXjhzBno0MFZB5Pkuutg+nS49lr/vitfPuja1Xnf6NEwapSTTIEzSrNqFXz6qTOyIyIiklNYazl16lSgw/BZWFgYJjO/kfTQsWNHBg0aRGRkJD169ACchKZChQrcfvvtHLl4rroXypcvz0MPPeTTM2XKlPH5PRd7/vnn2bx5M1WqVOGRRx7JdH9ZSQmNSDpiYqBNG2c6WZJ77nHKLhcunHXvLVwYXnwRuneHPn2cPWsA9u6F5s3hlVdg+PDMjQqJiIj4y6lTpyhZsmSgw/DZiRMnKOHtXPF0lChRgjZt2vDVV19x8OBB9u7dy9atWxkyZEiGyx7XrFmTyZMn+yU+b0VFRTFmzBgKFSrEtGnTcvxUOU05E0lDfLwzWuKZzHTuDLNmZW0y46lqVaeK2qhRzugNOGtpnn3WieWff7InDhEREUlfREQE8fHxREVFufaeyYnVzVLzww8/0LNnT0JCQpg+fToNGzYMdEjp0giNSBqeecapaJakd29n88zs3lsqNBSGDXP2u+nc2b1PzWefwdat8NVXkIM2HRYREcmz7r77bkqUKMEnn3xCdHQ01157bab2n9myZQujRo3y6ZkyZcowduxYn9/1888/065dO86fP8/EiRNp3769z30EghIakVRMmQJjxrjbbdvChx9CSADHNevXhzVrnAIB333nnFu/Hm65BebPh+uvD1xsIiKSt4WFhXEiCCvXhIWF+bW/ggUL0qlTJz788EMABg8enKn+Dh48yJQpU3x6plKlSj4nNJs2beKuu+7i77//5o033qBXr14+PR9ISmhEUvDjj9C/v7tdpw5MnRrYZCZJqVIwb56zOefrrzvnDh50NvicPdvZsFNERCS7GWP8thYl2HXv3p1Zs2ZhjKFbt26Z6qtZs2bYLN63YdeuXbRq1Ypjx47x4osv8vjjj2fp+/wtB/x4JpKz/P23s27m/HmnXbYszJmTsza2zJcPxo6FiRPd099iYuCuu5xiBSIiIhI4TZo04ejRoxw5coRKlSoFOpw0HT58mFatWrF//36GDh3KiBEjAh2SzzRCI3KRp56CnTud49BQ+OKLnLs+pXdvuPxyeOABOHvWScK6dHE+g2j9oYiIiATIww8/zLZt2yhSpAhHjx6lZ8+el9yT0TU52UUJjYiHRYvg/ffd7eHDoUmTwMXjjXvugR9+cEpLHz8OCQnQowfExTkbdYqIiIikJmnd05kzZ1Jdq5ORNTnZSQmNSKJTp5wRjyTXXQfPPx+4eHzRsKGzT83tt8Phw05Z59694cIF6Ncv0NGJiIjkXr6sbylfvnyWr4fx1ZIlSwIdQqZpDY1IoiefdDatBMif36lyVqBAYGPyRe3asGQJlC/vPte/v/N9iIiIiORWSmhEgF9/hY8+crdfeMEZoQk2117rjNRceaX7XO/ezjogERERkdxICY3kedbC0KHOJ0B4uFMSOVhVrw6LF0O5ck47IcEpFDB/fmDjEhEREckKSmgkz/v6aycBSPL6605Z5GB2zTXOxpslSzrtuDi47z5YuTKwcYmIiIj4mxIaydPOn4d//9vdbt3a+coN6tSBBQugaFGnfe4c3HsvbN4c2LhERERE/EkJjeRpEybAtm3OcUiIs1llbnLzzc4IVFJxgxMn4M47ITo6sHGJiIiI+IsSGsmzYmLgP/9xt/v1cyqF5TbNmsHUqWCM096zx0lqTp0KaFgiIiIifqGERvKsDz5wNqIEZ1qWZ3KT23TqBG++6W7/8Qd07uzsUyMiIiISzJTQSJ507pyz+D/JwIHuqmC51eDB8NRT7vbChcnbIiIiIsFICY3kSZMnw8GDznGhQjBkSEDDyTavvgr33+9uv/EGTJwYuHhEREREMksJjeQ5Fy7AmDHudu/euX90JklICEyZAtdf7z73r3/BsmWBi0lEREQkM5TQSJ4zYwbs3Okch4bCk08GNp7sdtll8NVX7iQuLs4ZtUn6MxEREREJJkpoJE9JSIBRo9ztbt2gcuWAhRMwV18NX37pLud89Ci0betUfhMREREJJkpoJE/59lvYsMHdHjYscLEE2i23wEcfudsbNjgJXnx84GISERER8ZUSGslTJkxwH997L9SqFbhYcoLu3ZMndV9/Dc89F7h4REREgo0xxvW1evXqVO+bOXOm677KOWh6yIkTJxg+fDgtW7akUqVKFClShCJFihAeHs5TTz3F0aNHAx1iupTQSJ6xd6/zA3uSgQMDF0tOMnKkk9wlGTUKZs8OXDwiIiLBKjIyMtVrU6dOzcZIvLd//35GjRrFr7/+SqVKlWjbti3Nmzfn6NGjvPbaa1x33XXszOELbZXQSJ7x0UfOGhqAKlXgjjsCG09OERoKkZEQHu4+99BD8NdfgYtJREQkmISGhlKnTh1mzJjBhRR2rT527BgLFiygXr16AYgubVdddRW//PILR44cYdmyZURFRfHNN9+we/duunfvTnR0NE/m8ApKSmgkT4iLgw8/dLcfftgpYSyOYsXgiy+gaFGnffo0dOwIZ88GNi4REZFg0a1bN44ePcrChQsvuTZjxgzi4uKIiIgIQGRpCwsLo379+oRc9INRoUKF+O9//wvADz/8EIjQvKYf6SRP+PprOHDAOc6f39l7RpKrUQMmTXK3f/9d0/JERES81bVrV4wxKU4tmzp1KkWLFqVdu3YBiCzj8ufPD0CBpLKoOZQSGskTPIsBdOwIZcsGLpYUnT0Lu3bB2rWwfj38+Sfs2QPnz2drGJ06wWOPudsffwwTJ2ZrCCIiIkHpqquu4rbbbmPOnDn8/fffrvM7duxg9erVdOjQgSJFigQwQt/ExcXx4osvAtCmTZvABpOOfIEOQCSr/fUXfPedu/3II4GLBYDjx2HZMli1yvnasAFOnUr53pAQZ8FP9epw3XXQvDk0auTsjplFxoyBNWsgqVDLwIFQrx7ccEOWvVJERHIDa1P/9ywnCwsDY/zSVUREBEuXLmXWrFn06NEDcBcKyOh0syVLltC8eXOfnqlUqRK7du3y+V19+vQhPj6eEydOsHbtWvbv30+jRo0YM2aMz31lJyU0kut9/LH7uFYtaNIkAEGcPw/z5sGUKTB3LqSwYDBFCQmwfbvzNX++U4Isf35nE5nOneGBB6BMGb+GWqAAzJzpJDBHj0JsrDOqtXYtlCjh11eJiEhucuoUlCwZ6Ch8d+KE3/6B69ixI4MGDSIyMjJZQlOhQgVuv/12jhw54nOf5cuX56GHHvLpmTIZ/NlgypQpxHtsSNesWTM+/vhjSpcunaH+sosSGsnVrIVp09ztPn389ksY75w758x3e/VVOHw4/fvDwpwk5uzZ1JOeuDhnhGfZMmd+WOvWMGAA3Hmn3yodVKwI06dDq1bOn+GOHU7ls9mzs/nPT0REJIiUKFGCNm3a8NVXX3Hw4EH27t3L1q1bGTJkCKGhoRnqs2bNmkyePNm/gaYiqULbgQMHWLlyJcOHD6dOnTp8/vnntG7dOltiyAitoZFc7ccfnaUp4Pwg/uCD2fTi+HhnhX316jBkyKXJzOWXOwtW3ngDVq50NsmJjYWTJ50SY3FxTjK0ZQvMmQOjRzubxRQvnryfCxfgm2+gTRuoXdsp5XbunF++hZYt4aWX3O05c2D8eL90LSIikmtFREQQHx9PVFSUq0BATqxulpYKFSrQsWNHFi1ahDGGnj178s8//wQ6rFQZa22gY5BUGGM21qpVq9bGjRsDHUrQevRReOcd57hFC/j++2x46c6d0L27k6h4KlgQOnRwhjpatoR8GRggvXDBmfs1c6YzhJJUus1TxYowYgT07Jmxd3hISHDyqHnznHb+/E6SmAPL6IuISBZLSEhg69atANSoUeOSMr95dQ2NMYbQ0FDX6EZsbCzly5enSpUqREdHU6pUKTZt2gTAwYMHqVChgk9rXLZs2cKoUaN8iqlMmTKMHTvWp2dS06ZNG+bNm8f3339PixYt0rw33b8jHsLDw9m0adMma214qjd5SVPOJNeKi4MZM9ztrl2z+IXWwuTJMHgweFQ3IX9+pxLBM89A+fKZe0e+fNCggfM1ZgwsXuxkbHPmOO8H2LcP+vWDsWOdqW7t22f4/6hDQpxlP3XrOrlTXJwzyrVunXvPGhEREcD5t0aLLSlYsCCdOnXiw8QN8AYPHpyp/g4ePMiUKVN8eqZSpUp+S2iS1uNkZP1PdtGUM8m1vv8ekv63V6AA3H9/Fr4sNhYiIpwNbjyTmc6dnRLMb72V+WTmYqGhzkjP7NmwdauzjsazTvzWrXDffc50tO3bM/yaMmUgMtKdE23b5ox8iYiISMq6d+9O6dKlKVOmDN26dctUX82aNcNa69NXRiqcpSQ+Pp4VK1YAULVqVb/0mRWU0Eiu5VkMoE2bLPyl0fHjzup5zxeWKAFRUc5X5cpZ9GIP11wD777rJE+9eiUvDjB/PoSHOwtiMrivTfPmzgBTksmTk3+7IiIi4takSROOHj3KkSNHqFSpUqDDSVNUVBR//PHHJeePHz9O//792bFjB3Xq1KF+/foBiM47SmgkVzpzBr780t3Osulmu3Y5+8IsW+Y+16wZ/PGHMzqT3SpVcooR/PGHM3qTJDbWWVfToIFzLQNGjHCqRSd55JFMDfyIiIhIDrBgwQLq1q1L1apVad++PV27dqVp06ZUqlSJSZMmceWVVzJjxgxMDi5zGtQJjTGmsDHmJWPMn8aYc8aYaGPMJGPMlT7209QYM8IY840x5ogxxhpjdnnxXKgxZogx5g9jzNnEZ2caY67N8DclfvH11+6ZX8WKOSM0frdjh5PMbNniPterF3z7rbMwP5Bq1XLiiIqCChXc53/7DerXd9bWeNSZ90b+/M6oTFiY046JgS5dMjzoIyIiIjlA3759GTBgAMWKFWPlypV89tlnrF+/ntq1azNy5Eg2btzItdfm7B9tg7YogDGmEPAD0BA4AHwFVAZ6AfcYYxpaa3d42d144Dof3x8CfAZ0AE4C3wBlgI5AG2NMc2vtGl/6FP+ZOdN9fP/9ULiwn19w4ADccQdER7vPvfgivPBCztmoxRhnlOiuu2DYMGc/HHBW9j/zDHz3nbM4xjPhSUflyk5l6AcecNo//+x8yz4WXxEREck1fKkYXL58eZ/uzw6NGzemcePGgQ4jU4J5hOY5nGRmNVDdWtvZWtsAGAqUBSb50Ne3if21BrwtHdcbJ5nZBtS01na01jYDOgFFgEhjTNAmjMHs7FlYsMDd7tTJzy9IWjOzwyNfnjDBmZOVU5IZT8WLw/vvO38oV3oMXi5eDNdf7yQ2PujUCfr2dbdHj/a5CxERERG/CcqExhhTABiU2BxorXWVlbLWjgPWA02NMV6tXrLWPmWtHWmt/RY47mUYTyR+PmWtPeTR1xfAHKAa0M7LvsSPvv/eWUMDcNllzv4zfnPmDNx9N2zY4D43diw8/LAfX5JFWrd21s8kDa+As+Fn69ZOwYCEBK+7Gj8ePEefe/aEY8f8F6qIiIiIt4IyoQEaAWHAdmvtrylc/zzx896seLkxpgpwLXAWZ6pZtr5f0vbVV+7jO++EQoX81LG1TuLy00/uc88+C0OH+ukF2aBkSWddzfvvOxt9gvN9jRjhzM2LifGqmyJFnG6SqkRHRztFAnLYKLqIiIjkAcGa0CStd1mXyvWk83Wz+P0brLVxAXi/pCIhwSkIkKSdP8fI3n8fpk51tx95BF5+2Y8vyCbGOLH/+CNUq+Y+P3s2NGzobDTjhbp14b//dbc//xw+/dTPsYqIiIikI1gTmqsTP/elcj3pfFYV/vbr+40xG1P6AnLuDkY51Jo1cChxAmBoqB+rm61eDY8/7m7fdpuzWWZOXDPjreuvd1b133WX+9ymTU5Ss3y5V10MGeLsUZNk0CCnkrWIiIhIdgnWhKZo4ueZVK7/k/hZLJe+X1LhOd3sttugVCk/dHr4sLMSPi5xMK5CBZgxw6ljHOxKlHCGtDx3zTx+3NnDJjIy3cdDQmDKFPempTEx0L27zxWhRURERDIsWBOaXMVaG57SF6BtC33kmdD4ZbpZ0rqZ/fuddr588NlnUL68HzrPIUJDYeRImD7dva7m/HmIiIBXXkl3YcxVV8F777nbK1bAa69lYbwiIiIiHoI1oUmqalYkleuXJX56t8I5+N4vKdi2DTZvdrf9ktBERTlrS5K89pqzmWZu9OCDTom40qXd555/3plql04FtC5dnK8kL7wAv6ZUrkNERETEz4I1odmT+JnaduxJ53fn0vdLCjxHZ+rWdTaBzJRDh5xFIUmaN4fBgzPZaQ7XqJFTLKB6dfe5t96C3r3hwoU0H333XWe0BpzZed26OXsCiYiIiGSlYE1ofk/8rJfK9aTz67P4/bWNMSktpMjq90sK5sxxH2d6dMZa+Ne/nPUk4GxoM3Gis2gkt6tWDVauhPoe2zhNmeKsIzp3LtXHSpZ0bkuqk7B5MwwblsWxioiISJ4XrD+drQROAVWNMdencL1j4ufXKVzLNGvtTmAzUBhIqY5Wlr5fLnX6NKxa5W7fm9kdgGbMgC+/dLdfew2qVMlkp0GkTBn44Qdo2tR9bvZsp2xcGnvVNG8OTzzhbr/9tjOLTURERCSrBGVCY609D7yT2HzXGJO0ZgVjzBM4+78stdau9Tg/yBizxRjzqp/CGJf4OcYYc7nHe+4D2gJ/AV+l9KD435Il7spapUpBvdTG7rzx99/Jfypv3twpDJDXFC8O8+fDPfe4z/3wg1MBLWnkKgUjR0KdOu52795OwikiIiKSFYIyoUn0CvATcCuwzRgzwxjzI/A6cATofdH9ZYAaQIWLOzLG9DXG/Jj4/DeJpysknUv8uvhH5EnAl8A1wBZjzGfGmMXA58BZIMJam/aiA/Gbb791H7ds6RTuyrBRo+DAAee4YEH46KO8MdUsJYULw6xZzoKYJGvWODWxDx5M8ZGCBeGTT5yCcAB79iTPD0VERET8KWh/SrPWngOaAy/j7AfTHmcjy8lAPWvtDh+6qwg0SPxKSlwKeJxrABS/6P0JQCdgKBAN3APUAb4AbrTW/pSR70sy5rvv3Md33JGJjnbtgrFj3e2hQ+H//i8THeYC+fM7GcrAge5zGzdCixbOHj0puP56p9JZkokTYd68LI5TREQkAIwxrq/Vq1enet/MmTNd91XOdOWirLVs2TJCQkIwxtC3b99Ah5OuoE1oAKy1Z621L1hrq1lrC1prK1hre1lr96Vw74vWWmOt7ZnGtbS+lqTwXLy1dpy1tra1trC1toy1tpO1dlPWfMeSkt274c8/3e1MJTRPPQWxsc5xhQowfHimYss1QkKcBTHPPus+t3kz3H47HDmS4iNPP528rkDfvmnOVBMREQl6kWlsSj116tRsjCTjYmNj6d+/f6DD8ElQJzQikHx0pnp1qFQpgx0tW+Zsmpnk1VehaNFMxZarGONstDlihPvchg3OHL9jxy65PX9+p+pZ0l6dBw7k/qrXIiKSN4WGhlKnTh1mzJjBhRS2OTh27BgLFiygXqYW+WaPV155hT///JM+ffoEOhSvKaGRoOe5fqZVqwx2Yq0zvSzJTTdB9+6ZiivXGjEi+UjN+vWpFgoID4eXX3a3IyOdJTkiIiK5Tbdu3Th69CgLFy685NqMGTOIi4sjIiIiAJF5b+PGjYwZM4Y+ffrQKIg2EldCI0EtPj55WeAMTzebOxd++cXdfvPNvFsIID3GOFmK5yYzv/3mZJMnT15y+xNPwK23utuPPJLq0hsREZGg1bVrV4wxKU4tmzp1KkWLFqVdpjfKyzrWWvr3709YWBijR48OdDg+0U9sEtTWrXMPDOTLB82aZaATa+HFF93te+9N/hN4DpGQkMCJEyfYsWMHW7ZsYefOnURHR3Pq1CmstdkbjDHOlDzPUa21a6F1azh1KtmtoaEwebJTMA2cJTf/+pfzxy4iIpJbXHXVVdx2223MmTOHv//+23V+x44drF69mg4dOlCkSJEARpi2998+Dh2JAAAgAElEQVR/n1WrVvH6669TqlSpQIfjk3yBDkAkMzzXzzRs6Gyd4rO5c53MKInnGpEA+euvv1i8eDF//PEHGzduZOPGjRw+fDjVxKVYsWJcffXVVKpUibp161K/fn3q169P5cqVMcZkTZDGOBuOXrgA48c759ascTbf/PZb8Pg/7WuugTFj4NFHnfasWTB9OnTtmjWhiYhI9rP2kt9pBYWwMOefNH+IiIhg6dKlzJo1ix49egDuQgEZnW62ZMkSmjdv7tMzlSpVYteuXV7fv3//foYPH07z5s3pHoRT7pXQSFDL9PqZlEZnPEtzZRNrLcuXLycqKoqFCxeyY4cvVcchJibGlfjM86iPXLFiRVq1akXr1q254447KFmypH8DNwbeeMOZ+/dO4l63K1dCp04we7ZTGSDRgAFOIrN4sdMeONAZUbviCv+GJCIigXHqFPj7n5nscOIElCjhn746duzIoEGDiIyMTJbQVKhQgdtvv50jqVQGTUv58uV56KGHfHqmTJkyPt0/aNAgzp07x/vvv+/TczmFEhoJWmfOwKpV7naG1s9cPDrjmdxkg6NHj/Lxxx/z0Ucf8adn7el0FChQgPPnz6d73759+5g0aRKTJk0if/78tG7dmq5du9K2bVsuu+yyzITuZgy89RacO+dsQgrOpjM9e8Knn7rWIoWEwKRJULcuxMQ4y2369XP+E2TVIJKIiEh2KlGiBG3atOGrr77i4MGD7N27l61btzJkyBBCM7jrd82aNZk8ebJ/A/Uwa9YsZs+ezQsvvECNGjWy7D1ZSQmNBK0ff4S4OOe4WDG48UYfO7h4dKZtW8imcoonTpxg7NixjB8/nn/++SfFe8LDw7n11lupXbs24eHhVK5cmZIlSxIWFkZoaCjWWs6fP8/JkyfZu3cve/bsYdu2baxbt461a9eyffv2ZP3FxcUxd+5c5s6dy2WXXUb37t0ZMGAAderUyfw3ZAxMmOD8muuLL5xz06ZBqVJOspOYsVSuDOPGOYkMOHnPpEkQRJUhRURE0hQREcGsWbOIiopi586drnM50enTp3n00Ue55ppreOaZZwIdToYpoZGgtXy5+/jWW52iAD754YdsXztz/vx5Xn/9dUaPHs2piyYaG2O44447eOCBB2jdujUVK1ZMsy9jDAULFqRcuXKUK1eOGy/K6A4fPsx3333HwoULWbBgQbJh7n/++YcJEyYwYcIEmjRpwpNPPsk999yTufU2oaFOXeaTJ92l5955B8qUSfZn26ePM/Vs/nynPWSIU/U5w/sHiYhIjhAW5vxeK9iEhfm3v7vvvpsSJUrwySefEB0dzbXXXpup/We2bNnCqFGjfHqmTJkyjB07Nt371q1bR3R0NJUrV6Z169bJrh08eBCAb775hmbNmlG+fHmioqJ8iiPbWGv1lUO/gI21atWykrIWLax1hlmsHTkyAx3cc4+7g9at/R7fxX766ScbHh5ugWRfpUuXts8995zdsWNHlr07Li7Ofvfdd7Z37942LCzskhgAW7duXTt9+nR74cKFzL3s9Glrb7rJ/WcL1r71VrJb9u2ztkQJ9+WWLa1NSMjca0VEJGvFx8fbTZs22U2bNtn4+PhAh5NjADY0NDTZuX79+rn+fR3p8UPKgQMHLGArVarkdf+LFy9O8d/ttL687d+Xvr3p05e/I7Vq1bLARuuHn5lVtlmC0vnzsHq1u92kiY8d/Pmns3gjyZAhfokrJefOnWPo0KHccsstbNy40XU+LCyMl19+mZ07d/Lyyy9TpUqVLIshX758tGzZkokTJxIdHc3EiRMv+W3R+vXr6dKlC3Xq1GH27NlJSbXvihVz5pJde6373ODBzuhNoiuvhLffdl9etMi9/EZERCTYde/endKlS1OmTBm6deuWqb6aNWvm8w/43lY4S6vvjz/+GIA+ffr41GcgKKGRoLRuHZw96xwXLAg33eRjB0llhgFq1cpgibT07dmzhyZNmjBu3DgSEhJc5wcMGMDOnTt57rnnKFasWJa8OzVFihShd+/e/PLLL6xateqSTb42b95Mhw4daNSoEStXrszYS8qUcUrQXX21+1zPnsnK0nXr5ixbSjJ0KOzZk7HXiYiI5CRNmjTh6NGjHDlyhEqaU53llNBIUPJcP3PzzVCokA8Pnzjh7PSY5PHHs6TM1uLFi6lfvz6//PKL69w111zDsmXLePfdd/1fQtlHxhhuueUWZs+ezR9//EGXLl2SraFZvXo1jRs3JiIigv379/v+gooVnY2CypZ12hcuQMeO8Ntvie936ggklcqMiYH+/bXhpoiIiPhGCY0EpWXL3Me33ebjwx9+6NR8BihdGrKg8sgHH3zAHXfcwdGjR13n/vWvf/H777/TxOf5cVmvdu3aTJs2jd9++4177rkn2bXIyEhq1KjBmDFjiEsqK+et6tWd6WdJm2zGxMDdd7uGYipUSD5YtnAhJI5wi4iIiHhFCY0EnYQEWLHC3fYpP4iLS75445FHoHBhv8UGMHr0aB555BHi4+MBKFiwIJMmTeK9996jsJ/f5W9169bl66+/ZtmyZcnW2Pzzzz8MGzaMm266iZ9//tm3Tm+8ET77zKmCBnDggJPUnDwJQPfuTjPJE0/Avn2Z/U5ERESyh7WWCxcueHVv+fLlc/x6lCQ9e/bEWstHQbDIVQmNBJ0NG1w/CxMS4pRs9trXX7t/Ws6f39mu3k+stTzzzDM8/fTTrnNXXnklK1asoFevXn57T3Zo0qQJa9as4YMPPqB06dKu87///jsNGzZkyJAhnEka5fLG3XeD5+7DGzdChw4QG4sx8L//uctmnjoFDz+sqWciIiLiHSU0EnQ8p5vVq+cU1fKa528ZOnZ05jz5gbWWJ554gldffdV1rlq1aqxcufKS/WGCRWhoKP3792fbtm3079/fdT4hIYE333yTG2+8kXWe+/ikp18/ePZZd3vJEujVCxISuPJKZ8PNJPPmwaefZv57EBERkdxPCY0EHc+CAD5NN9u711mkkSRpu3o/ePXVV3nzzTdd7dq1a7N8+fJcUdmkZMmSfPDBByxZsoTq1au7zm/evJmGDRsyevRo1/S6dL38MvTo4W5Pnw6JOxP36gWee3o99hhER/vjOxAREZHcTAmNBBVrM1EQYPJkZwEOQNWq0LSpX2L66KOPeNZj5KF+/fosXbqU8uXL+6X/nKJp06b8/vvvDB061FUNLS4ujqeffpoWLVqwe/fu9DsxxinK0LKl+9zo0fDBB65LSSNuJ086S5w09UxERETSooRGgsr27XDwoLvduLGXDyYkwMSJ7nafPs4CnEyaM2cODz/8sKtdvXp15s+fT6lSpTLdd05UqFAhxo4dy6JFi6hYsaLr/LJly6hbty6RkZHpb8hZoAB8/jnUqeM+N3AgLFrEVVfB66+7T3/9NUyb5udvQkRERHIVJTQSVH780X1cs6azf6NXvv8ekkYQQkLgoYcyHcvvv//Ogw8+6Nows0KFCixcuJCySfuu5GItWrRg/fr1dO7c2XXu9OnTRERE8NBDD6VfMCAszFkoc8UVTjs+3lnTtGULffsmH8AZPDh5EisiIiLiSQmNBJWffnIfN2jgw4OexQDatHH/IJ1BJ06c4L777uPs2bMAhIWFsWDBAipXrpypfoNJyZIlmT59OlOnTqV48eKu859++im33HIL27dvT7uDihWdIZikPWpOnYJ77sEcO8qHH0LRos7p48dhwABNPRMREZGUKaGRoJKhhOboUfjyS3e7b99MxZCQkEBERAQ7duwAwBjDzJkzqVu3bqb6DUbGGLp168b69etp7DH/b/369dSvX5+5c+em3UG9ejB1qru9fTvcdx+VK8Ty2mvu019+CTNn+jl4ERHxSdL6ScA1O0HEk2eRIM+/L1lNCY0EjdhY+O03d9vrhGbaNGdDTXDKNHvu4pgBL7/8MvPmzXO1X3nlFVq1apWpPoNdpUqV+OGHHxgyZIjr3KlTp7j33nsZMWJE2lXQOnSAUaPc7eXL4eGH6d/P0ry5+/TAgXD4cBYELyIiXjHGUKBAAcDZcFnkYqdPnwacTcWzM6HJl21vEsmk335z5yWFCiVfU54mz1Xl3btDvoz/tf/hhx/4z3/+42q3a9cu2UaaeVn+/PkZN24cDRo0oE+fPq5/7F566SXWrFlDZGRk6sUSnnoKtmxxKtEBTJlCyLXXMnHiMOrUgX/+gWPHYNAgjdSIiARSsWLFOHbsGIcOHQLgsssuI8QPRXYkeFlriY2NJSYmhuPHjwPOtPTsZNKtSCQBY4zZWKtWrVobN24MdCg5wltvOXuTANx6K6xc6cVDO3Y4JZqT/PYbXHddht5/6tQp6tSpw969ewG45ppr+PnnnwlL2uJeXDZu3Mh9993Hn3/+6TpXuXJl5s6dS3h4eMoPnT8Pd9yRvC73F1/wTvR9PPqo+9Rnnzn1A0REJPvFx8ezZ88ezp07F+hQJIcqUaIE5cqVSzfRDQ8PZ9OmTZustan8YOA9pdQSNDK0fiYqyn187bWQiXUujz32mCuZyZcvHzNmzFAyk4rw8HB+/vlnOnTo4Dq3a9cubr31VhZ6bm7qqUAB+OKL5AloRAQDGqxNtt/QgAHOsigREcl+oaGhXH311ZQuXdo1/UwkNDSU4sWLc+WVV1K+fPlsH7XTlDMJGhlKaKZPdx936eJs7JgBs2fPZsqUKa72iy++yA033JChvvKK4sWL88UXXzBmzBiGDx+OtZbTp0/Tpk0b3nrrLQYMGHDpQ2XKwNy50LChU/Xs7FlC2rdl4oy11G1VnrNn4cgRePTR5P9pRUQk+4SGhnL55Zdz+eWXY61Nf/8xydWMMdm6XibFGPSXMOfSlDO3Y8eS7zmzcyekWyF5w4bkC222bYNq1Xx+9+HDh6lduzZHjhwBoEGDBqxYsYJ8mViLk9fMmTOHrl27JltE+uijjzJu3LiU/xwXLYI773T2pwGoV4/xnVfx+LCCrltmzXLqCYiIiEjw0ZQzyXPWrHEfX345VKrkxUOe081uvDFDyQzA0KFDXclM4cKF+eSTT5TM+Kht27YsX76cK6+80nXu7bffpm3btq6KKMm0bAnvvONur1vHo6u70qiR+xcw//qXk+iKiIhI3qaERoKC53Szm2/2YuaYtZdON8uApUuXMtVjn5RRo0ZRvXr1DPWV191www2sWbOG+vXru87Nnz+fRo0asWfPnksfeOQRePxxVzNk9iwmXfcWhQo57UOHkl0WERGRPEoJjQQFn9fP/PyzU+EMnOync2ef3xkXF8fAgQNd7fr16ydri++uuOIKli1bxn333ec6t2HDBm655Rb++OOPSx8YOzbZvkHV33ucVzr97mpPnQpff52lIYuIiEgOp4RGcjxrk0858yqh8Zxu1qQJeEx18tb48eNJWr9kjOG9994jNDTU534kuSJFivDZZ58xfPhw17no6GiaNGnC0qVLk98cGursI1SzpuvU4583pmEd91qchx+GEyeyPGwRERHJoZTQSI7311+QuE8TADfdlM4D1sKXX7rbDz7o8zv37dvHiy++6Gr369ePm2++2ed+JGUhISH897//5YMPPnCVdjx16hStW7fmiy++SH5zWBjMmQMlSgAQevZvPj7WloIFnfU0Bw5o6pmIiEhepoRGcjzP0ZmaNV0/16Zu/XrYtcvdbtfO53cOGzbMVZGrdOnS/Pe///W5D0lf//79mTVrFoUSF8bExsbSqVMn3nvvveQ3XnMNzJgBiclPzegfeOmKD1yXP/kEvvkm28IWERGRHEQJjeR469a5j2+80YsHvvrKfdygAVxxhY/vW8e0adNc7VGjRlG6dGmf+hDvtWvXjkWLFlGyZEkArLUMHDiQ5557LvneBq1awWuvuZpP7BzEzZfvdLX794eTJ7MtbBEREckhlNBIjvfrr+5jr/aynD3bfdy+vc/v81zbUbduXXr16uVzH+KbRo0asXz5cipWrOg6N3LkSPr370980l40AEOGQI8eAOQjno8Pt6FAPud6dDQ88US2hi0iIiI5gBIaydGs9TGh2b07+QM+TjdbtGgR3377ras9atQoFQLIJuHh4axevZrwcPf+Wh999BFdu3bl/Pnzzglj4IMPXJUharGZ/yQ877r/449h/vxsDVtEREQCTAmN5Gi7dyefRnT99ek8MGeO+7h69WTVsdKTkJDA008/7Wo3a9aMO++80+vnJfMqVqzI8uXLady4sevczJkz6dChA2fPnnVOFCoEs2ZBhQoA/DthDDfmcyex/frBqVPZGraIiIgEkBIaydE8B1sqV4bEZRap85xu1q6dFztwun322WesXbvW1R49ejTGh+fFP0qWLMnChQtp3bq169y8efO46667OH36tHPiiiuc/9YFCzpTzy50p4BxRnH274ehQwMRuYiIiASCEhrJ0Xyabnb8OHjuY+LD+pn4+Hief949daljx44q0xxARYoU4auvvuL+++93nVu6dCktW7bk2LFjzombb4b//Q+A2mxkhH3Rde/EibBwYXZGLCIiIoGihEZyNJ8SmnnzIGkBeblyXu7A6Zg5cybbtm0DnD1SXnnlFR8jFX8rWLAgUVFRPPTQQ65zP//8M82aNePAgQPOiR49XJUAnmIM9fnFdW/fvpp6JiIikhcooZEczaeExnO62b33OrvMeyEhISHZPjNdunShRo0aPkQpWSVfvnxMmjSJQYMGuc5t2LCB2267jT179jgnRo+GVq2cqWf0Ij/O1LN9++DJJwMRtYiIiGQnJTSSYx054qyHSJJmQhMXBx7VyXypbjZnzhw2bNjganuWbZbACwkJ4a233uKZZ55xnfvrr79o1qwZu3fvhnz5ICoKqlWjDht4gZdc9334YfK/FiIiIpL7KKGRHMtzdKZs2XT2x1y9GmJinOOCBaF5c6/eYa1l5MiRrnaHDh2SlQ2WnMEYw8iRIxk1apTr3M6dO2natCk7d+50qkXMng1FizKM0dyAezfWvn0hqZaAiIiI5D5KaCTHuni6WZoFxzxXgN92G1x2mVfv+O677/jlF/e6i2effdbHKCU7DRs2jHHjxrnau3fvpmnTpmzfvh3Cw2HKFPJzgcn0JB9xAOzdq6lnIiIiuZkSGsmxfFo/s2CB+9ij3G96PBf/33nnndSvX9/rZyUwhgwZwvjx413tvXv30qxZM6eow333wTPPUJc/eJ6XXff873+waFEgohUREZGspoRGciyvE5pDh2Cde4oRXm6GuWbNGpYvX+5qa3QmeAwePJh33nnH1d63bx/NmjVj69at8NJLcNddDOdVrsf9l6hPH/esRBEREck9lNBIjvT335BYRRlIJ6H57jv3ccWKUKuWV+/w/C3/rbfemmx3esn5Bg4cyIQJE1zt6OhomjVrxuY//4Rp08hfrXKyqWd79sBTTwUqWhEREckqSmgkR/r9d7DWOS5aFKpVS+Pmi6ebpbnYxhEdHc3MmTNd7cceeyyDkUogPfzww3z44YeYxP/mBw8epHnz5mzcvx9mz+a6ojt4Dve0wgkT4PvvAxWtiIiIZAUlNJIjeU43u+46CEntb2pCQvK6vF5ON5swYQIXLlwAoGLFinTo0CGDkUqg9e3bl4kTJ7qSmkOHDtG8eXP+SEiAKVMYzqtcx2+u+/v0sZp6JiIikosooZEc6fff3cdpTjf79VdnwxpwNtJs2TLdvs+dO5dsqtLAgQPJnz9/BiOVnKBXr15MmTKFkMTM98iRI7Ro0YLfq1alwLNPJZt6tnu34emnAxmtiIiI+JMSGsmRPPa5pG7dNG70LNfcoAGUKJFu31FRURxJTIIKFy5Mv379Mhil5CTdu3fn008/dSU1R48e5fbbb2dDp05cf/eVPMN/Xfe+9x4sXhyoSEVERMSflNBIjmNt8oSmdu00bvZcP+PFdDNrbbJiABEREZQuXToDUUpO1LVrV6ZNm0ZoaCgAx44d4/ZWrdj6wgs8W20mdVjvurdPj/P8/XegIhURERF/UUIjOc6ePST7QTM8PJUbY2Jg1Sp324uEZtWqVfz2m3s9xeDBgzMYpeRUnTt3JjIy0jVSc/jwYZq1b0/0+NeYXHgAoThrp3buK8DwoecDGaqIiIj4QVAnNMaYwsaYl4wxfxpjzhljoo0xk4wxV2agr5LGmPHGmN3GmNjEzzeNManOYTLGVDfGfJx473ljTIwx5mdjzBBjTIHMfXd5l+fozNVXQ/Hiqdy4YgXExzvHJUpAvXrp9v3hhx+6jlu0aEHtNId/JFh17tyZKVOmJKt+1rh/f64Y25XhvOq6753/FWDp4oRAhSkiIiJ+ELQJjTGmEPAD8DxQFPgK2Av0An41xvyfD32VAdYAg4ELwGwgBngM+MkYUyqFZ24FfgV6Av8kPrMKCAfGAd8aY/Jl8NvL07yebrZkifu4aVOnKEAaTp06laxUc//+/TMWoASFiIgIJk2a5Grv37+fBqNHM+iRY9TmD9f53vef5J9/AhGhiIiI+EPQJjTAc0BDYDVQ3Vrb2VrbABgKlAUmpfXwRd4EqgGzgBqJfdUG3gaq4yQoF3sHKAIMt9bWstY+YK1tndjPDqAp0D1j31re5nVC47mqu1mzdPudPn06Z8+eBaBUqVK0b98+YwFK0OjZsycffPCBq71nzx5unT+b9298zzX1bMeJUgx7YGegQhQREZFMCsqEJnE616DE5kBrrWvFhbV2HLAeaGqMqe9FXxWALsB5YIC19oLH5SeBI0CEMeZyj2eKAjcAZ4Axnv1Za6Nxkh2Am3z81gQvE5pTp2DtWne7efN0+/3oo49cxz169KBgwYIZjFCCSf/+/XnnnXdc7R27d/Po0QU8GeZOdN6dV4XvJ+8NRHgiIiKSSUGZ0ACNgDBgu7X21xSuf574ea8Xfd2J8+ew3Fp7yPOCtTYW+BoIBe72uBQHeDPx/pgX94iHCxdg82Z3O9WEZsUKZ1NNgFKloE6dNPv99ddfWeuRAPXp0yeTkUowGThwIG+88Yar/duuXawLG08d486ee/cL5dTe04EIT0RERDIhWBOa6xI/16VyPel8WjuYZLivxERnGc6Us6c8bzbGXAEMxEl6PvXi/eJh+3aIjXWOQ0KgZs1UbvScbta0qXNzGiZOnOg6btiwoYoB5EGPP/44Y8a4B1S/3bON8LJDXRtu7rlwBU80/smdKIuIiEhQCNaE5urEz32pXE86XykL+3oEpwjBq8aYTcaYmcaYBcBfgAHaWGv/9OL9GGM2pvQFVPXm+dzEc7pZtWpQuHAqN3omNOlMNzt79ixTp051tbWRZt715JNP8sorr7jaUYe/5a4ir7vak/bcwTc9ZgQiNBEREcmgYE1oiiZ+nknlelLNomJZ1Ze1divQGGcE51qgE9AaKAQsBjZ68W65iFfrZ06cgF89Zhqmk9DMmjWLU6dOAVC0aFEeeOCBTEYpwezZZ59lxIgRrvbXZ56nZoh7gLZvZDOORS5I6VERERHJgYI1oQk4Y0wL4HcgP9ACKA5UAUbilI5eaYwp601f1trwlL6A7VkUfo7lVUKzfDlY6xyXKQO1aqXZ56efumf+PfjggxQtWjSNuyUvGDFiBMOHD09sXWBfQgQFOQfAQSrwaM8Y2Lo1cAGKiIiI14I1oUmqalYkleuXJX7GZEVfifvSfIaTzNxlrV1srY2x1u6y1j4PvAtUBv7txfvFw0aPca1UE5qLyzWnsX7m0KFDfPfdd652RERE5gKUXMEYw8iRI3n88ccB+JvNhPGc6/r0C534vMV7cFpFAkRERHK6YE1o9iR+VkzletL53VnUVxugFPCjtXZ/Cs98lvh5mxfvl0SxsfCnx6ojrxKadKabzZw5k4TERd5XXXUVTZo0yWSUklsYYxg3bhx9+/YF4DBvUIoVruv/in6OQ50GqUiAiIhIDhesCc3viZ/1UrmedH59FvWVlOScSuWZpPMlvXi/JNq6FeLjneMCBZyiAJc4fhzWe/ynSGdDzcjISNdxly5dCEmnGprkLcYYJkyYQNeuXYEEjtOT/InL5o5Slke+7YD9z0uBDVJERETSFKw/3a3ESRqqGmOuT+F6x8TPr73oawHOnjJNPDfPBDDGFMTZyyYemOdx6WDi5w3GmNAU+kzaUHOXF++XRJ7rZ2rWhPz5U7hp5Ur3+pnLL4drr021v7/++ouffvrJ1e7WrZufIpXcJDQ0lMmTJ9O+fXtgO3E86bo2mw5MfWk7zJ4duABFREQkTUGZ0FhrzwNJW3+/a4xJWueCMeYJnD1jllpr13qcH2SM2WKMefWivg4A04ECwHvGmHwel8cAZYGp1trDHucXALE4RQBeNsa4/hyNMTWApF/pfo54zauCACtXuo8bNQJjUu1v2rRpHv3Vpm5db7Ylkrwof/78REVF0apVK2ACoSxyXXuUt9nXbRhs2hS4AEVERCRVQZnQJHoF+Am4FdhmjJlhjPkReB04AvS+6P4yQA2gQgp9PY5TUex+YIsxJsoY8wcwGNgGPOF5c2IS9G/AAsMT3/+5MWYxzhS2K3BGdCb74fvMMzKU0KTCWptsupkzpUgkdQULFuTLL7+kSZPGxNObEJyCAKcoQd8z47Ft28HJkwGOUkRERC4WtAmNtfYc0Bx4GWcPmfY4m19OBupZa3f40NdR4GbgbZyRmg5AGPAWcLO19ngKz7yDU655Nk6FtHY4621+BQYCba21FzL47eVJmze7j1OsxHz+PPz8s7udRkKzdu1a/vSoMKCERrxRpEgR5s6dy003lSeBx1znF3InH25vDl27uhd6iYiISI4QtAkNgLX2rLX2BWttNWttQWttBWttL2vtvhTufdFaa6y1PVPp67i1drC19urEvq621j5mrU31V7LW2iXW2g6J781vrQ2z1t5irX3PWqufenwQGws7PFLQFJfGrFvn3AhQsCDccEOq/XlON2vcuDGVKlXyU6SS21VunUsAACAASURBVBUvXpwFCxZQu/ZaPJfhPcE4ts/fCs8/H7jgRERE5BJBndBI7vHXX+7quPnyQZUqKdzkOd3sppucpCYF1lq++OILV1ujM+KrUqVKsWjRd/zf/40GjgHwD0XpwSfEvzoaZs4MbIAiIiLiooRGcgTPTdmrVUulwtmqVe7jdKab7dnjbC9kjKFDhw5+ilLyknLlyrF0aRRly77gOreKRozhKejVK3n5cBEREQkYJTSSI2zZ4j6uWTOFG6xNPkJz662p9uU5OtOoUSPKly/vhwglL6pYsSKrVz9B4cLuv1Mv8BLrztSA9u3h2LEARiciIiKghEZyCM+EpkaNFG7YuRMOHXK3U0loLp5udv/99/spQsmrqlatyuLFtQkJcZbmXSA/EUzl7M4D8OCDcEG1P0RERAJJCY3kCJ5TzlIcofEcnaleHcqUSbGfjRs3sm3bNlf7vvvu81OEkpc1aFCDCRPOudqbqcVwXoVFi2DYsABGJiIiIkpoJOCs9WLKmZf7z3iOztx4441cffXVfohQBPr1q0aXLtGu9ngeZxG3w7hxMHVqACMTERHJ25TQSMAdOgSnT7vbKU4587IgwKxZs1zHmm4m/jZp0hVUrvyPq92TyZygBAl9+8LatQGMTEREJO9SQiMB5zk6c/nlULLkRTecPAkbNrjbqayf+euvv1jvUXlKCY34W6FC8OWXl5Evn1NjfD8VGci7hMTGEnfvvXD4cIAjFBERyXuU0EjApVsQ4McfnXlpAKVKpXJT8ulmderU4ZprrvFjlCKO66+HV15x/1/ndLoynQfJf+AA5+69F+LiAhidiIhI3qOERgIu3YIAP/3kPr7lFghJ+a/tl19+6TpWMQDJSv/+NzRu7G4P4D32cSWF1qzh74cfDlxgIiIieZASGgm4dAsC/Pyz+/jmm1Ps49ChQ6xZs8bV1maakpVCQ+GTT6BoUad9kpL0ZDIJGIp+/DEx48cHNkAREZE8RAmNBFyaU86sTZ7Q3HRTin3Mnz8fmzgt7aqrrqJu3bp+jlIkuSpV4K233O3vacnbPApAwSFDiFm0KECRiYiI5C1KaCSgzp6F3bvd7UtGaPbsSb7QOpWE5ptvvnEdt2nTBmOMH6MUSVnPntC+vbs9jNFspBYFrOVcmzac2b49YLGJiIjkFUpoJKC2bXOv9y9QACpXvugGz9GZKlVS3FDz/PnzLFy40NVu06aN/wMVSYEx8L//OdX5AGIpRFciiaUAZc+fZ1e9epw7eTKwQYqIiORySmgkoDynm11zjbM2IRkvpputWLGCmJgYAAoVKkSLFi38HKVI6sqWhUmT3O31XM9wXgWg1unTrKxTh7jz5wMUnYiISO6nhEYCKt0KZ14kNHPnznUdt2jRgiJFivgpOhHvtGkDAwe622/wBAtoDcDt+/Yxo0ED4uPjAxSdiIhI7qaERgIqzYIACQnwyy/uthfrZ+655x4/Rifivddeg/BwdzvCTOEwZQHo8ttvvH333SQkJAQoOhERkdxLCY0EVJojNFu3QuJUMkJCoH79S57ftm0bf/75p6ut9TMSKIULw/TpULCg0z5my9E1ZAoWCAV6ffsto3v2dFXjE/l/9u48zsby/+P46xpmDFkqskbZQqGQNfsSaaNVy5dKvxSifFMqpdLGNyEl7fpK+yJJfCtrtoSsWRKSkCWyzpiZ6/fHdc657zObMfvyfj4e53Hu67qXc98ezTSfc13X5yMiIplDAY3kGGvDA5okIzT+6WZ16nhFP3z8ozP16tWjSpUqmXyXImlXr54bqQn6PuEyRpv+AJQCrps0iecGD86ZmxMREcmnFNBIjvnrLzh82GvXrJnogFNcP6PRGckN+veHrl299iOFx7CaugDUBBqPGsV/nnsuZ25OREQkH1JAIznm11+97dKl4YwzEh3w44/edjIBzeHDh5k3b16orYBGcgNj4J13oFw51445UYgeZ33LMaIB6AQUeeQRxo8fn3M3KSIiko8ooJEc4w9oqldPtDM2Fn7+2Ws3aZLk/Llz53LixAkASpUqRbNmzbLgLkVOXdmyMHGi1163pzyDqn4Sag8AVvTrx7vvvpvt9yYiIpLfKKCRHOMPaGrUSLRz9WoX1ICruFm/fpLzv/3229B2+/btKVy4cBbcpUj6dOkC993ntSdsuYIp5/QNtccD79x+O59++mn235yIiEg+ooBGcszmzd52koDGv37mwgtdUJOIP6Dp1KlTJt+dSMY9/7z7zzforsMv8Udpl9s5EvjEWh6+6SamT5+eMzcoIiKSDyigkRyT6gjNSerP7Nixg3Xr1oXaCmgkNypSBN5/36V0BtizrxB31viBuEiX2/ks4LO4OHpecw1z5szJsfsUERHJyxTQSI5JdQ3NihXedjL1Z7777rvQ9rnnnkv1JBcQyR3OPx9efNFrz1xyOi9fPz/Urg+8ERPDVVdcwZIlS7L/BkVERPI4BTSSI/bvh7//9tphIzQnTsCaNV77oouSnJ94upkxJgvuUiRz9OkDV1/ttR/6tDHLe44JtbsDg48coUuXLqxcuTL7b1BERCQPU0AjOcK/fqZECTjrLN/OX37xEgIULgwXXBB2rrU2bIRG080ktzMG3nwTKlZ07dhYuHHhAP659LrQMY8Blx44QKdOndjgrzgrIiIiqVJAIzki8fqZsAEWf7rm8893CxF8Vq9eze7duwEwxtC+ffssvFORzFGmDEyeDBGB37q//mroU+J9bO06oWMmAlX27KFjx45s3bo1J25TREQkz1FAIzkizetnGjRIcq5/ulmjRo0oXbp0Jt+dSNZo2xYef9xrf/hZJG/9a06oqmxRYCpg//iDDh068Oeff+bAXYqIiOQtCmgkR6Sa4cw/QpOG9TMiecnQoS6wCbp3eFnW/OcbN70SqAh8Bez+7Tc6derEnj17cuI2RURE8gwFNJIjUgxorA0PaBKN0Bw/fpx58+aF2gpoJK8pVMhNPQuuGzt+HG4Y1ZQjo18PHdMAmAysX7eOzp07c+DAgRy5VxERkbxAAY3kiBSLam7dCv4/3hKN0CxevJhjx44BUKxYMVq0aJF1NymSRSpWhP/+12v/8gvcu/x2+Pe/Q31XAyOAFStWcPnll3PkyJFsv08REZG8QAGNZLtDhyCwph9ItIbGPzpTtSqUKhV27ty5c0Pbl1xyCUUSJQwQySu6dIGHHvLa77wDk+qNhKuuCvU9ANwJLFy4kKuuuioUzIuIiIhHAY1kO//oTHS0l8oWOGlCAH819bb+hQgiedDw4dC8ude+p18E6x9/P2xkcjzQHpg1axbXXHMNMTEx2X6fIiIiuZkCGsl2iTOcRfj/K0wlIUBMTAyLFy8Otdu0aZNFdyiSPSIj4YMPQknOOHIEbrzjNI59/BVUqOCOAT4DagEzZszg+uuvJzZYp0lEREQU0Ej2S3H9DKQ6QvPjjz9y/PhxAIoWLUrjxo2z6A5Fss8557jpZkGrVsGgF8+GqVOhaFEATgemAaWBr776iptvvpm4uLicuF0REZFcRwGNZLsUa9Ds3Qt//OG1E43Q+KebtWjRgqioqKy5QZFsdvXVMHCg154wAT7+7WJ4771QXw3gcyAK+Oyzz+jZsyfx8fHZfasiIiK5jgIayXYppmz2TzcrUwYqVQo7z58QQNPNJL8ZMQIaNfLa//d/sPnCa+C550J9rYFgcucPPviAO++8k4SEhGy9TxERkdxGAY1kuzQFNA0agDGhZmxsLAsXLgy1lRBA8psiReCjj6BECdf+5x+44QY4PvAhuO220HG9gIcD2xMnTuSee+7BWpvdtysiIpJrKKCRbHXsWPissrCAxr9+JtF0s6VLl4ZS1kZHR9OkSZMsvEuRnFG9Orz5ptdevhwGDDTw2mvQunWo/1mgR2D79ddfZ8CAAQpqRESkwFJAI9lq61Zvu3BhqFzZt3PVKm87UUDjn27WrFkz1Z+RfOuGG6B/f6/9xhsw8f0o+PzzsG8A3jWGloHtl19+mcGDByuoERGRAkkBjWSrLVu87SpVXFADQGwsrF/v7axfP+w8f0Cj6WaS340aBU2beu177oGVf5SG6dOhdGkAoqxlWqFCnBc6ZxRDhw5VUCMiIgWOAhrJVv4RmnPP9e3YuBGCaWgLF4bzzgvtOnHiBAsWLAi1lRBA8ruoKPjkk1DswvHjcO21cOCsmi6dc2CEslR8PN9HRXFW4Lxnn32W4cOH58xNi4iI5BAFNJKt/CM0Vav6dqxZ423XquX+ogtYtmwZR44cASAqKoqm/q+uRfKpypVd0c1gbozNm11uANu8BUyaFDru7NhYvi9alKKB9rBhwxgxYkS236+IiEhOybGAxhhT1hjzfE59vuSMFEdoVq/2tuvVCztn/vz5oe2mTZtStGhRRAqCTp3gySe99pdfwn/+A1x/PYwcGeqvd+wYU4sXD/1CHzJkCKNHj87WexUREckp2R7QGGMqG2PGAVuAwdn9+ZKz0jRCU7du2Dn+dM2tWrXKojsTyZ0efRQuu8xrP/wwzJkDPPCAW1wT0PHwYd4sVSrUHjRoEK+88kr23aiIiEgOyZSAxhgTYYy52RjzrjHmm8B7D2NMhO+YysaYN4Ffgb5ANPBFZny+5B1pGqHxBTTW2rCApkWLFll2byK5UUSEm2F2zjmunZAAPXrAzl0GXnoJLr88dOztBw/yxJlnhtr9+/dXUCMiIvlehgMaY0xhYAYwCfgX0DnwPhn4JHBML2AdcDsQCXwJNLDWXpfRz5e849Ah2LfPa4dGaA4fDh+68U0527x5M3/99Veo3axZsyy+S5Hcp3Rp+PRTb2nZ7t0uvfMJWxg+/BAaNgwd+/jff9M7mE0AF9SMHz8+u29ZREQk22TGCE0/oCMQA0wA7gUeAX4AuhljJgBvA6cB3wINrbXXWGtXpXA9yaf8ozNFikD58oHGunXejtNOCxu68Y/O1K5dm9K+P9RECpKLL3YDMkE//OCmn1G8OEyb5vKgA8ZaXj9yhCvKlg0d269fPwU1IiKSb2VGQNMDiAfaWGv7WmtfsdY+b61tA7wO3BU4brC1tou19udM+EzJg/wBzTnnuKk0QPh0swsu8O0ID2guueSSrL1BkVzurrugZ0+vPWqUq7dJhQquRk1gDU3E8eNMiYvjkgoVQscqqBERkfwqMwKaOsBCa+3SZPb9J/C+3lo7KhM+S/Iw/6yysPUzqSQE8Nef0foZKeiMgVdfDU8EeNttgZq0F1zgopvISAAK7d/PrCJFuLBixdCx/fr149VXX83emxYREclimRHQlAC2prAv+Cfsykz4HMnj/CM0aclwduDAAdauXRtqK6ARgWLF4LPPoGRJ1z50CLp1g4MHgfbt4c03Q8dGbd3K4jJlqFWpUqivb9++CmpERCRfyYyAxuCmnCVhrbWBzeOZ8DmSx6U4QpNCDZolS5YQ/E/ozDPP5LzzzsvaGxTJI2rWhHff9dobNsCtt7oMaPTsCc8+G9oXvWoVy6tXp+rZZ4f6+vbty4QJE7LxjkVERLJOjhXWlIIn2ZTNe/a4lE1BvhGaxOmaIyL0n6tIULdu8PjjXnvaNHjiiUBjyBDo3z+0r9i8eaxs3JizfSM199xzj4IaERHJFzLrL8Rexpj45F6ATWV/XEY+1BhT1BjzlDFmozHmuDHmT2PM28aYSic/O8m1zjDGjDXGbDPGxATexxhjTj/JecWNMcOMMauMMYeNMQeNMWuMMa8YY4qn/+nyn2SLavqnm5UuDeXKhZpaPyOSumHD4Morvfbw4YEkAcbAmDFwnZcZv8QXX7Dqiis42zdSo6BGRETyg8wKaEw6X+n+fGNMNDALeAwojqttsx1X62aFMabaKVyrDPAjMACIA6YAh4CBwBJjzJkpnFcVWAU8gUtL/Q0wF1drpy+QajBUkBw4EJjjHxAaofEHNPXquT/EgLi4OJYsWRLapYBGJKmICHjvPahd2+vr2TPwY1WokKvI2aZNaN8Zr73GittvTxLUvPbaa9l41yIiIpkrwwGNtTYiI68MfPRQoBmwCDjPWnujtbYp8G/gLFztm7QaA9QAPgdqBa5VFxgHnAe8mPgEY0wRXABTBbjbWlvdWnu9tfYqa20toB6wP/2Pl7/4p5sVLQqhEhkpJARYs2YNhw8fBqBQoUI0btw4629SJA8qWRKmTPGSBBw54qaj/f03EB0NX34J9euHji/z9NMsHTw4LKi5++67FdSIiEielScXJRhjooDgBPF+1trDwX3W2hdxoyZtjDGN0nCtCsBNQCzQ11rrnwY3GNgD3GqMKZvo1IFALeBFa22SvwSstWustUdP4bHytcQJAQIDMSkmBPBPN2vQoAHFihXL2hsUycNq1YL33/d+rjZvhptugvh4XG2ab74JFd7EWsoPHszi555TUCMiIvlCngxogEuAUsBma+2KZPZ/Gni/Mpl9iXXB/TvMt9bu9u+w1sYAXwGFgK6Jzvu/wPu4tN50QZZsymZrwZeWmQsuCG0uXrw4tK3pZiInd/nlbg1N0MyZ8OijgUbFiq7jzMDs2dhYKvXrx4Lx45MENSq+KSIieU2GAxpjTJWMvNL5sRcG3pensD/YXz+F/Rm6ljGmMm6K2h/W2u3GmEuMMSOMMROMMQ8ZY2qk4XMLlGRTNu/cCf/84+3wBTQ//vhjaLtp06ZZe3Mi+cQjj8C113rtESPgo48Cjdq1XSq0okVd+59/qNKnD/MnTQoLavr168eYMWOy76ZFREQyqHAmXGMrLpNZeth03kMwEPojhf3B/nOy6FrnB97/NMa8gksA4Pe0MWaItXZUGj4fY8zaFHZVT8v5eUGyIzS//OJ1li8Pp7scCgcOHGDjxo2hXU2aNMn6GxTJB4yBiRNdXZrg8rTbb3exzIUXAs2bw8cfu0U28fGwcyfn9unDvM8+o81117F9+3YA7r//fmJiYnjooYdy7FlERETSKjOmnM1L5rUEL5PZAWBl4PV3oI/AMfPT+ZnBdMgprVE5EngvkUXXOiPw3hC4G5flrDJQAQj+BfCCMebyNHx+gZDsCI0/oPGlafrpp59C22eccQbVq+ebuE4kyxUv7pIEBL4f4NgxF7/s2xc44IorwL9WZuNGqt5zD/OnTaNq6NsGGDJkCE899VSouK2IiEhulRlZztpaa9sFX8BVuLTFa4DLrLWlrbUNA68yuDUrqwPHXJXRz88hwX+3wsBr1tonrbV/WGt3WWtHAqMD+x9Jy8WstRck9wI2Z8G9ZztrUyiquX6911mnTmjTP92scePGmFAGARFJi+rV4cMPXVpncD9/118PJ04EDujdG55+2jth+XLOufde5s2cSc2aNUPdw4YNY+jQoQpqREQkV8uKpADDcVOl2ltrZybeaa39H9AxcMwz6fyMYFazlFJfnRZ4P5RF1zrs234nmXOCfU0D9XIKtH374LDvXyzZKWe+EZqlS5eGtjXdTCR9OneG55/32rNnQ79+7gsGwC24uf9+74B58zh70CDmfvcdtX0/j88++yyDBw9WUCMiIrlWVgQ01wCzrLV7UzogsG8W0D2dn/F74P3sFPYH+7dl0bX821uTOSfYVwhItihnQeIfnSle3Eu0lJYRGgU0Iun3wANwyy1e+403ILTe3xh44QW47TbvgGnTqPDII8yZNYu6vrpQo0aNYuDAgQpqREQkV8qKgOYs0rbQvzBQJp2fsTLw3jCF/cH+VVl0rfXA8cD2GSTlD2IOJ7O/QEk83cwY4OBB+PNPb0cgoNmxYwd/+vpVUFMk/YyBN98Ef+bzf//bJTsD3Jy0N96A7r7vliZPptyzzzJ71iwuuuiiUPe4ceO4++67SUhIyJ6bFxERSaOsCGi2Ah0CqY2TFdjXgeRHN9JiAXAQqG6MuSiZ/dcF3r9Kw7VmAAlAq8TFM40xRXC1bOKB6cH+QH2a4HS6tslcs03g/Tdr7T/J7C9Qfv/d2z4nmCtuwwavs3hxqFQJCJ9uVrlyZcqXL58NdyiSf0VHwxdfeD971rqim6uCX9EULuyqcnbo4J308suUefllZs2aFfalwuuvv07v3r2Jj4/PvgcQERE5iawIaN7CrTuZa4zp6V9DYowpYoz5FzAXt2blrfR8gLU2Fng50HzFGBNc54IxZhCuZsxca+0yX39/Y8x6Y8xzia61E/gAiALGG2P8o0sjcSNO71lr/0p0GyMD748ZY87zfU5V3DoigAnpeb78JpAJFoDKwTA38fqZwMJ/TTcTyXxly7pRmRKBXI2HD8OVV8LuYCnhYNTj/5l76inOePddvv3227DithMnTuRf//oXJ0IZBkRERHJWVgQ0LwIfAefiFscfMcbsMsbswqVGnhjY92ng2PR6Gpf6uQWwyRjzkTFmMTAK2APckej4MkAtXGrlxO7DZRS7FlhvjPnQGLMaGABsAgYlPsFauxB4CrfGZoUx5n/GmG9wU9iqAd9k8PnyDX9AUyVY9SeNGc5EJHPUrRue+ez331065+PBybMlSsD06XD++d5J999PqS++YMaMGbRu3TrU/cEHH3DddddxPHSyiIhIzsn0gMZam2CtvQm4CfgBOAGUDbxOBPputtbeaK1N92Rsa+1xoB1uNOQo0A1X/HIi0NBa+9spXGsv0AQYhxup6Q6UAl4Cmlhr96dw3jBcELQMaIabarYZuB+4ylqreRmETzlLcYQGSEhICKtBoxEakczVtSuM8pX7XbwY7rjDl/msdGn43/98udWB3r0p8f33fPPNN3Ts2DHUPXXqVK644goOHy7wywRFRCSHmazOWhOYwlU60NxnrY3L0g/MR4wxa88///zz165dm9O3kiEVKsCuXW577lxo3RqoVQs2bnSdn38O3buzYcOGULpYYwwHDhygZMmSOXPTIvmUtXDPPeG1NZ98Eh5/3HfQr79Cy5benLSoKJg+neOXXEKPHj348ssvQ4c2b96c6dOnc3qwkqeIiEgaXHDBBaxbt25doPZihmR4hMYY844xppt/HYuftTbOWrs78FIwU8DExHjBDARGaGJjYbOvZmggiPFPN6tdu7aCGZEsYAyMGxeeA2DYMPj4Y99BNWq4kZpgkBIbC1ddRfRPP/HJJ59w8803hw5dtGgR7dq146+/Ei8zFBERyR6ZMeWsF/AZsNcY87Ux5m5jTEo1XaSA2bHD2zYmkMzs118hmCWpcGH3xxMqqCmSXSIj4ZNP4LzzvL5evcD3nQLUrw9ffw3FAjWHjx6Frl2JXLGCSZMm0adPn9ChP//8M61bt+aPP/7IngcQERHxyYyA5mygHzAbaA+MB7YZY5YbY4YZYxplwmdIHuVfP1O+vJu5EpYQoEYN99cV4QGNEgKIZK0zznCZz84IVNI6fhyuuir8Z5YWLWDqVChSxLUPHYLOnYlYuZJXX32VBx54IHTohg0baNmyJZv9o68iIiLZIMMBjbX2T2vtBGttV9xamWuBd4GKwDDgR2PMH8aYV40xXQO1XaSASFPKZiA+Pp6VK1eGui+++OJsuDuRgq1mTbeErXAgWf3u3S5xwIEDvoM6dHApnQNfPHDgAHTqhFm3jpEjRzJ8+PDQodu2baNVq1bk9XV/IiKSt2RqljNr7VFr7RfW2jtw6ZFb4Oq1/A30wRW63GuM+dwYc3viQpaS/6Q1ZfOGDRs4duwYABEREdSrVy+b7lCkYGvbFib4KmatXQvXXOOWzYRcdhl89BEUKuTa+/ZBhw6YTZsYOnQoY8aMCR26c+dOWrduHZaxUEREJCtlRR0aAKyz2Fr7sLW2HlAdV8/lR+ByXFHNP40xC40xnbLqPiRnpTVl8/Lly0NdderUoVhw3r6IZLneveHRR7327Nlw552+dM4A3bvD5MleIZvdu93ozZYtDBw4kLfeegsTKJC7f/9+2rdvz9y5c7PvIUREpMDKsoAmMWvtFmvtWGttB+As4GZcAc5aQPPsug/JXklGaBISkh2hWbFiRairQYMG2XR3IhI0fDjceqvXnjTJZT8Lc+ON8PbbXvuPP6B9e9i+nTvuuIMPPviAwoH5a4cOHaJz585MmTIl629eREQKtGwLaPystf9Yaz+01t6CK7j5ak7ch2S9JCM0O3bAkSNeZzIjNA0bNsymuxORIGPgrbegXTuvb/hw1xemV6/wOWpbt7qgZudObrzxRqZMmUJ0dDQAMTExXHvttbztD4JEREQyWY4ENH7W2nhr7Z6cvg/JGkmSAgSLaYKruFmiBNbasBEaBTQiOSMqyiUJOP98r69PH5g5M9GBffqAb90Mv/4KHTvCnj1cfvnlfPvtt5QqVQqAhIQEevfuzciRI8nqQs4iIlIwFc7oBYwxszJwug1MQZN86J9/4OBBr12lCjBlk9cRKIKxZcsWDvoOvOiii7LpDkUksdNPh+nToVkzVxQ3Ph6uuw7mz4ewH82BA+HYMXj4Yddetw46dYLvv6dly5bMmzePzp07sytQWfehhx7ir7/+YuTIkURE5Ph3aSIiko9kOKAB2qbjHAuYwLvkU/7RmchIKFsW2OQLaGrWBMLXz1SvXj30za6I5IxzznE1NVu3djNEDx92ic4WLoSqVX0HDhniCtg8+aRrr1zpRmq++4769euzcOFCOnXqFKpNM2rUKPbu3csbb7xBZDANtIiISAZlxtdkzU/h1QIYBRzPhM+VXC7xdLOICJINaPzrZ5QQQCR3aNgQPvnEy9S8axd07gx7Ek8QHjYMHnrIa//8s8t+tncvVatWZcGCBWGjru+++y7XXHNNKE27iIhIRmVGYc0laXkB5YHxuNTNRYH5gKab5WPJpmz2r6FJZoRG62dEco/LLoM33vDamzbB5Ze7EZsQY+C55+DBB72+lStDQU25cuWYM2cObdq0Ce2eNm0al156KQfCKniKiIikT5ZPZDbGXG6MWQp8DlwELAI6WWvbWGvnZPXnS85JkrI5Lg5++83rDKyh0QiNSO51++3wzDNez1hnmwAAIABJREFUe+lSt6YmrPCmMfD88+EjNatWuexne/ZQqlQpZsyYQbdu3UK7f/jhB1q3bs3OnTuz/iFERCRfy7KAxhjTxRizBJgKNAKWAJ2ttS2ttd9n1edK7pFkhOb33+HECddhDFSvzs6dO9m9e3foOAU0IrnPww/Dvfd67ZkzXTHOhATfQcGRmiFDvL7Vq91IzZ49REdH88knn9C7d2/f7tW0aNGCDRs2ZP1DiIhIvpXpAY0xppMxZgHwNdAY+Anoaq1tYa39NrM/T3KvVFM2V64M0dFhozOVKlWiXLly2XeDIpImxrgszTfc4PW99174gEzowGefhUce8fpWr3YjNX/9ReHChXnjjTcY4gt6tm7dSosWLVi4cGHWPoSIiORbmRbQGGPaGWPmATNwCQCWA1daa5taa2dk1udI3pFkytlJMpxpdEYk94qIgP/+18UmQS+8AKNGJTrQGHj6aXj0Ua9vzZpQUGOM4bnnnuPFF18M7d6/fz8dOnRgypQpWfsQIiKSL2U4oDHGtDbGzAa+A1oCK4GrrbWNrbVfZ/T6kjdZm8wIzaakNWj8IzRKCCCSuxUpAl98EV6P5oEHYOLERAcaA8OHw9ChXt/atdCuHQSmmN5///18+OGHREVFAXD8+HGuvfZaxo8fn7UPISIi+U5m1KGZg6sncwx4FZgCWGNMi7ScbK3VPIN8aM8eiInx2hqhEckfSpaEb76BFi1gyxbX17s3lCoF3bv7DjQGnnrKC27AFd9s1w5mzYLy5bnxxhspV64c3bp14+DBgyQkJNCvXz9+//13nn32WRXgFBGRNMnM/1sUw6VknodLyZyW17xM/HzJRfwJAUqUcH/sJE7ZfPDgQbZu3RrquiisDLmI5Fbly8P//gfBJW8JCdCjB3z3XaIDg0HNsGFe3y+/QNu28McfALRt25YFCxZw9tlnhw4ZMWIEPXv2JDYslZqIiEjyMmOEZh5uhEYkJMn6mdhY8AUv1KzJ6tWrQ82SJUtyzjnnZNv9iUjG1Kjhgpo2beDAAfcj3q2bC2qaNUt08BNPuPcnn3TvGzZAq1bw/fdQrRoXXHABixYtomvXrqHfC5MnT2bXrl18/vnnlCxZMtueS0RE8p4MBzTW2raZcB+SzyRJ2bxli5fjtVAhqFqV1bNmhY6pX78+xpjsvUkRyZD69WH6dOjYEY4ehSNHoGtXmDsX6tVLdPATT0BkpLeuZutWaN3aRUC1a3P22Wczf/58unfvzuzZswH4/vvvadWqFd988w0VK1bMzkcTEZE8RBOUJUsEZpMAySQEOPdciIpi1apVoa769etn272JSOZp3hymTIHA2n7+/hsuvRQ2b07m4EcfhdGjvfaOHS6oCfwuKFWqFN988w033XRT6JBVq1bRrFkz1q5dm4VPISIieVmGR2iMMT0zcr619r8ZvQfJfXbs8LYrVSLJ+hlAAY1IPtGpE3zwAVx/vRuI3bXL9c2fH/j597vvPihWDO6+26VD3LPHramZMQOaNKFIkSK89957VK5cmZEjRwKwfft2WrRowaeffkqnTp2y/flERCR3y4w1NBNJ3xoaEzhPAU0+lCSgWRae4SwhISFsDU29JPNTRCQvueYaePNNuOMO196yxU1FmzPHSx4QctddcNpp0KsXxMe7YZ0OHeDrr6F1ayIiIhgxYgSVK1dmwIABWGv5559/uOyyyxg/fjx33XVXdj+eiIjkYpkR0DxF0oCmOnArcBT4H7A10H8OcClwGvAekNykBMkH/vzT265YEfgwvAbNtm3bOHToUKirbt262XdzIpIlbr/dJQgYNMi11693IzWzZ0Pp0okOvuUWN1Jz441w4gQcPgxdurhCN507A9C/f38qV67MzTffzNGjR4mPj6dPnz5s2rSJESNGKK2ziIgAmbCGxlr7hLX2yeALeB+4EhewVLHWXmOtHRR4XQtUASYBVwAfZPTzJfex9uRTzvyjM1WrVlUWI5F84v77vbIzAKtXuzU1Bw4kc3D37jB1KkRHu/axY3DVVW5RTsDVV1/NvHnzqFChQqjvhRde4LrrruPo0aNZ9BQiIpKXZMXXW88BfwO3W2v3J95prf0b6A0cCBwr+cyBA+7vkqBKZx4Lz+Ncs6bWz4jkY0OHwiOPeO3ly+Gyy8A3KOvp0sWtnyle3LVjY+G669yinIBGjRqxZMmSsN8VX3zxBW3btmXXrl1Z9BQiIpJXZEVA0xZYbK2NT+kAa20csBhokwWfLznMPzpTpAiUPuCbWRgZCVWqKKARyeeeftqt/w9avBiuuMKld06iTRuXvvn00107Pt5NSXv11dAhlStX5ocffqBr166hvqVLl9K0aVPWrFmTRU8hIiJ5QVYENEWBCic9CsoD0Vnw+ZLD/AFNxYpgfvWtn6lWDQoXVkAjks8ZAy++6JKZBc2b54pvHj+ezAlNm7oMAmed5drWQt++8NRTbhsoUaIEX375Jf379w+d9vvvv9OiRQtmzpyZdQ8jIiK5WlYENKuAVsaYjikdYIzpALQOHCv5TJKEAL/95nXUqMHRo0fZ5KtLowxnIvmTMfDKK3DbbV7ft9/CtddCTEwyJ1x4oYt6zj7b6xs2DO69143aAIULF2bcuHGMHTs2lBTg0KFDXH755bzqG9EREZGCI6vW0EQA04wxbxtjOhtjagdenY0xbwFf49I2P58Fny85LElCAH9AU60a69atIyEhAYDo6Ghq1KiRvTcoItkmIsKlc+7Rw+ubPj2VoKZ2bViwwL0HvfIK3Hxz2AkDBgzgyy+/5LTTTgMgPj6evn370r9/f06cOJFFTyMiIrlRpgc01tqpQF8gAbgNmA6sDbymA7fj0jzfGzhW8pkkAY2/ZHj16mEZzurWrUuhQoWy7+ZEJNsVKgT//a+rVRP09depBDVVqriqnE2ben0ff+wW4fgyC1xxxRXMnz+fSr7qna+88gpdunRh//4kOWlERCSfypIk/tbaCUBtYDgwB9gYeM0N9NWx1o7Pis+WnHeyERqtnxEpeCIj4cMPXabmoK+/dgnNkg1qypRxiQICNWkA1+7QAfbsCXU1aNCAH3/8kcaNG4f6Zs2aRZMmTVi3bl0WPImIiOQ2WVaVzFr7u7V2mLW2g7W2TuDVPtC3Nas+V3Je2Bqa8gmwdavXoYBGpMCKjISPPgoPaqZNSyWoKV7c1am56Savb+lSaNkStm0LdVWsWJG5c+dy8803h/o2b95Ms2bN+Prrr7PgSUREJDdRmWXJdGEjNFF7XBXwAHvuuaxcuTLUVkIAkYIlpaDm+utTCGqiouC991xigKCNG6FFC1i7NtRVtGhR3nvvPZ577jmMMYBLFnDllVfyn//8BxvIlCYiIvmPAhrJVCdOwO7dXrtS7BavUaECu/75h3379oW6FNCIFDzJTT/76qtURmoiImDsWFfcJujPP6FVK1i4MNRljGHIkCFMmTKF4oFCndZaHnzwQXr16sXxZPNFi4hIXqeARjLVrl2hkhEAVPxnvdeoVo21vm9Uy5Urx1nBmhMiUqBERSUNaqZNc3Vqjh1L5gRj4NFH4bXXXIAD8Pffbk3N55+HHXrVVVexaNEizj333FDfpEmTaNu2LTt37sz8hxERkRylgEYylX/9zBlnQNE/wotq+hfpXnDBBdl4ZyKS2yQX1MyYAZdfDocPp3DSXXe5jGdRUa59/Lgb2hkzJuywunXrsnTpUtq0aRPqW7JkCY0bN2bp0qWZ/CQiIpKTFNBIpjpZhjP/CI0CGhGJinJram64weubPRu6dIGDB1M46dprYeZMOP1017YW7r8f7rsvVIAToEyZMvzvf/+jT58+ob4dO3bQqlUr3n777Sx4GhERyQkKaCRTpaWoZtD555+ffTcmIrlWZCRMngz/+pfXt2ABdOrkZpUlq21bd9A553h9Y8e67AJHj4a6oqKimDBhAq+88kqo5lVMTAy9e/fmnnvuITY2NvMfSEREspUCGslUqQU0NtEIjQIaEQkqXBgmToQ77/T6li6F9u3Dys6EO/98WLwYGjb0+r74ItmT+vbty3fffRe2bm/ChAm0bduWP/1zZUVEJM9RQCOZyh/QVCwdA3v3htp7SpTgb9/XrZpyJiJ+ERFuzX///l7fzz+7wZgU1/KXLw9z57qFN0FLlkDz5i69s0/btm1ZtmxZWBHORYsW0ahRI3744YfMexAREclWCmgkU/m/6KwU+ZfXiI5mjS+4KVu2LKVLl87GOxORvCAiAl56CQYP9vrWrXMZmrdsSeGk4sVhyhS4+26vb/NmV6vGl9YZoHLlysybN4/evXuH+nbt2kW7du0YP3686tWIiORBCmgkU4VNOUvY7jWqVWOtMpyJSBoYAyNGwOOPe32bN0PLli64SVbhwjB+vDsxaN8+N/3s00/DDo2OjuaNN95gwoQJREZGAhAXF0e/fv244447OJZs3mgREcmtFNBIpgoLaI796jWUEEBEToEx8OSTMHKk1/fnn9C6tVtbk+JJDz4IH3zgpXWOiXEp1EaMCCuSZYyhT58+zJ07lwoVKoT6J06cSKtWrdi2bVsWPJWIiGQFBTSSaf75J7x2RKW/13iN6tWVsllETtngwfD66y5WAW/QZc6cVE7q0QO+/dYVwwIXyAwZArfd5gIcn+bNm7Ns2TIuueSSUN+yZcto1KgRM2bMyNRnERGRrKGARjKNf/1M4cJw1s5VobatWlUZzkQkXf7v/1wBzsDsMA4fdnVqvvoqlZNat3brZ6pV8/r++18XDf31V9ihFSpUYNasWfT3ZSPYt28fXbt25fHHHyfeV9tGRERyHwU0kmn8080qVICILZtD7QNnnsn+/ftDbQU0InIqbrgBpk6FokVdOyYGuneHSZNSOal2bZfxrE0br2/hQmjcGFatCjs0KiqKcePGMXHiRKKjowGw1jJ8+HA6d+7M7t27M/mJREQksyigkUwTtn6mooWtW0PtjXFxoe2zzjorrBaEiEhadOkC//sflCzp2vHx0LMnjBqVykllyriT/AVufv/dZUCbOjXJ4b169WLJkiXUrFkz1Pf999/ToEED5s+fn0lPIiIimUkBjWSasIDmzKPgC2JWHDgQ2tbojIikV8uWbv2M/zuRBx6Af/8bEhJSOCkqyi3EGT3a5YUGOHIEunVzWQcSpWquX78+P/30E9dff32ob+fOnbRr146RI0eSkOIHiYhITlBAI5nGv4amYvR+X6MiqzZtCjWVEEBEMqJBA1iwAKpW9fpefNGN1sTGpnCSMXDffTBtGpQo4fqshYcegttvT5IsoGTJknz00UeMGzculNo5Pj6ehx56iG7duoUVCRYRkZylgEYyTdgIjfGV9VbKZhHJZDVruuUwF13k9U2eDFdcAYcOpXLiZZfBokXh0dC770KHDkmSBRhj6N+/Pz/88ANVqlQJ9X/11Vc0bNiQpSnmjxYRkeykgEYyzU5fDFPxhK+GgzKciUgWKF8e5s51icuCvv0W2rWDVNfwX3AB/PgjtGrl9S1Y4JIFLF+e5PAmTZqwYsUKLr/88lDf1q1badmyJePGjcMmmrImIiLZSwGNZJpdu7zt8oe9oppHypZl7969obamnIlIZilZEqZPd6VngpYtg0suAd9M16TKlIHvvoM77vD6fv/dnZhM6rQzzzyTqVOn8txzzxERWIcTGxvLgAED6NatG/v27cukJxIRkVOVpwMaY0xRY8xTxpiNxpjjxpg/jTFvG2MqpeNaZxhjxhpjthljYgLvY4wxp6fx/ChjzDpjjDXGxJ38jPzF2vARmvL7vSlm2wsVCm2XLl1aGc5EJFMVKeKmm913n9e3eTM0b+4GXlIUFQVvvukW4ASTBRw/7hbjDBwIJ06EHR4REcGQIUOYNWsWFSpUCPVPnTqVCy+8kLlz52biU4mISFrl2YDGGBMNzAIeA4oDXwLbgduBFcaYaqmcnvhaZYAfgQFAHDAFOAQMBJYYY85Mw2UeAWqfyjPkJwcPhq+prbD759D2Rt+O888/HxMs+S0ikkkiIlxcMmKE17dvn1sa88knqZxoDNx/v0vtXLq01//SS9CxY7Jz19q0acPPP//MZZddFurbsWMH7du354knniAursB9pyUikqPybEADDAWaAYuA86y1N1prmwL/Bs4C3j6Fa40BagCfA7UC16oLjAPOA15M7WRjTB3gYeCNU36KfMI/3axwYcuZu7wRmpUHD4a2a9WqlZ23JSIFiDHw4IPwwQdu8AXcFy033AD/+U+S7MzhOnRwc9UaNPD65s2DRo3ceptEypYty7Rp0xg1alQoC1pCQgJPPvkk7du3Z/v27Zn4ZCIikpo8GdAYY6KA/oFmP2vt4eA+a+2LwCqgjTGmURquVQG4CYgF+lpr/V+tDQb2ALcaY8qmcL4BXgcOAEPS8Tj5Qth0szLxRBCo02AMP/ryOdeuXWAHsUQkm/To4ZbHnOkbW3/wQejbN6w8VlLnnOPmqPXs6fXt2OGSB7z1VpLDIyIiGDRoEAsXLqR69eqh/vnz53PhhRcyZcqUTHgaERE5mTwZ0ACXAKWAzdbaFcns/zTwfmUartUF9+8w31obNrfAWhsDfAUUArqmcH4foCXwb2ttgS1MEJYQoORRr1GpEmt9K3M1QiMi2aFVK5fWuZpv8vGECXD11XD4cMrnUbQoTJzoppwVLuz6YmPhzjvhnnuSLXRz8cUXs3z5cm655ZZQ399//0337t3p168fx48fz5yHEhGRZOXVgObCwHvS/Jrh/fWz8lqB0Z3nge+tte+l4bPyrbARmqIHQtsJlSuzdevWUFsBjYhkl1q1XMmZpk29vunToWVLl9AsRcbAvffC999DWd/g/IQJ0LZteBXhgJIlSzJp0iQmTpzIaaedFuofP348TZo0Yc2aNRl/IBERSVZeDWiCFc7+SGF/sP+cLL7Wy0A00DcNn5MiY8za5F5A9ZOenEv4R2gqRHjF6Q6eeWaoRkNkZCRV/cXsRESyWNmyMHs2XHON17dyJTRpAosXn+Tk1q3dupomTby+RYtcNc9vv01yuDGGXr16sXz5chr41uKsXr2aiy++mNGjR5OQkJDBJxIRkcTyakBTPPB+NIX9RwLvJbLqWsaYq4FrgOettRvT8Dn5WtiUs3gvNtwVXJkL1KhRg8LBKRwiItmkaFGX6WzwYK9v92432PL++yc5+eyzXfXO3r29vj17oHNneOIJiI9Pcsp5553HokWLGDhwYKgvJiaGQYMGcemll/LHHyl9fyYiIumRVwOaHGWMKYEbndkIPJfR61lrL0juBWzO6LWzi3/KWYWjv4W2f/N9G6npZiKSUyIiYORIt7Y/kJSMmBi45RZ4/HFIdeAkOhreeMO9ihRxfdbCk09Cly7w119JTilSpAhjxozhm2++oXz58qH+77//nnr16vHhhx9m4tOJiBRseTWgCS7pLJbC/uAE5kNZdK1ngbNxWdFikp5S8ISN0BzcENpec8j7Z1OGMxHJaXfc4TKg+UvODB8ON94IR1Mapwe3rubOO908tRo1vP7vvnNT0ObNS/a0Ll26sHr1aq7xzXk7cOAAN910E7feeisHDhxI9jwREUm7vBrQBJdznp3C/mD/tiy61pXAceAxY8wc/yuwv5Cv76I03EOeF5YUYN/a0PbSPXtC2xqhEZHcoHVrWLIE6tTx+j791GVGSzVZALjgZdkyuO46r2/nTmjf3lX1TGaop0yZMnz66ae88847FC9ePNQ/efJk6tevz5w5czL2QCIiBVxeDWhWBt4bprA/2L8qC68VDbRJ5hUUbJ+ehnvI02JjXUXuoAoJ3vzwecpwJiK5UPXqbn1/585e3/Llro7mSeOLkiXh449daufg/LX4eBgyBK66KvwXYoAxhttuu42VK1dyySWXhPq3b99O+/btGTx4MDExGvAXEUmPvBrQLAAOAtVTGAEJfnX2VRquNQNIAFolLp5pjCmCG42JB6YH+62151prTXKvwCHxvr45p/ZoeU/i6ePlcfPP4suWZY9vypkCGhHJTUqVgmnTYMAAr2/vXujYEcaOdctkUhRM7fzDD64gZ9DXX0PDhm4IKBnVqlVj7ty5PPPMM6EkKdZaXnjhhVA9GxEROTV5MqCx1sbiFuUDvGKMCSX9N8YMwtWMmWutXebr72+MWW+MeS7RtXYCHwBRwHhjjD8N10jgLOA9a23SVZ8ChE83O71YDNG4bxmP+Capn3XWWZzpL9stIpILFC7sgpd33vHW+8fHw333Qa9ecOzYSS7QpIkb2rniCq/v999dsZsUpqAVKlSIRx55hMWLF4etLVyzZg1NmjRh2LBhxCZTwFNERJKXJwOagKeBJUALYJMx5iNjzGJgFLAHuCPR8WWAWkCFZK51Hy6j2LXAemPMh8aY1cAAYBMwKGseIX8ISwhQ7J/Q9p5iXp4Fjc6ISG52220wf77L0hw0aVIainACnHkmfPmlS6NWqJDri4tzU9A6dUq2ECdAo0aNWLZsGf379w/1xcfH89RTT9GkSRN+/vnnjD2UiEgBkWcDGmvtcaAdMBxXQ6YbrvjlRKChtfa3lM9Ocq29QBNgHG6kpjtQCngJaGKt3Z+pN5/PhCUEKOzNHd/qO0YZzkQkt2vcGH76ySUNCErzupqICFfoZs4cqFzZ6581C+rXh6+SnwFdrFgxxo0bx6xZszj33HND/StXrqRx48Y89dRTnDhxIr2PJCJSIOTZgAbAWnvMWvu4tbaGtbaItbaCtfZ2a22SqmXW2icCa1puS+Fa+621A6y1VQLXqmKtHWitPaWcmoHPKFDVI/0jNBXsjtD2et9cDY3QiEheUK6cy8R8771eX5rX1YAb0lm5Eq691uvbt88lC+jfP8U5bO3atWPVqlXcfffdob64uDiGDRtG06ZNWb16dQaeSkQkf8vTAY3kDmFTzo572a2X7/cGthTQiEheERnpEphNnJh0XU3PniepVwNwxhnwySeuEKdv6i2vvOLW3KxZk+xpJUqU4NVXX+Xbb7+lSpUqof4VK1bQqFEjnnnmGeLi4jL2cCIi+ZACGskw/5SzCoc2hraX+CIdTTkTkbymVy+XxMy/rua996BZM9iwIeXzAK8Q57JlrnZN0Jo1bm7b+PEpDvd07NiR1atXc9ddd4X6Tpw4wdChQ2nWrBmrVqWlIoGISMGhgEYyLGyEJsGbcrY18B4ZGUnVqlWz9Z5ERDLDxRe7mMS/rmb1atf/4YdpuEDt2rB4sRveCTp+HPr1g27d3Hy2ZJQsWZLXXnuNmTNnUtm3JmfZsmU0atSIRx99lOPHj6fzqURE8hcFNJJhYUkBAjVoYkqU4Eigr3r16qF6CyIieU3Zsm5dzf33e32HD8NNN0Hfvi4+SVWRIjB6NEyf7i4WNHUq1K3riuGk4NJLL2X16tX07t071BcXF8ezzz7LhRdeyLx589L5VCIi+YcCGskQaxMlBcBFN/tLlAj1abqZiOR1kZHw4ovw+eeuIGfQq6/CJZfA5s1puMhll8GqVdC5s9e3ezdceaWbnvbPP8meVqpUKd58801mzpwZlglt48aNtGnThrvvvpuDBw+m78FERPIBBTSSIQcPQkyM1w6O0OyIjAz11axZM7tvS0QkS3Tv7qVyDlq+HBo2dMHOSZUr50ZqRo/2Mg4AvPUWXHghpDLicumll7JmzRoGDRpERIT3v+/XXnuNOnXqMGXKlHQ8kYhI3qeARjLEP90sMiKOM3GZzX711U1QQCMi+Um1arBggVsGE/TPPy5T8/33Q2zsSS4QEeHW1CSOjLZuhbZt4YEHUpzHdtpppzFq1CgWLVpEvXr1Qv07d+6ke/fuXHfddez0/2IWESkAFNBIhvinm5WL+psIXNae1YcOhfoV0IhIflOkCLz8Mnz0Efhm2DJmjEsgsG1byueGnH8+LFoETzwBhQq5Pmth1CgX6CxfnuKpTZo0YdmyZTzzzDMU8Y30fPbZZ9SpU4c33niDhISE9D2ciEgeo4BGMiQ8IcDu0PY6X0BTo0aN7LwlEZFsc8MN8NNPUL++17dkicvU/OmnabhAZCQMG+YyofnXG65bB02bwvDhkELtmcjISB555BFWrlxJa18atoMHD3LXXXfRsmVLVq5cmc4nExHJOxTQSIaEJQQ48Xtoe3vgvWjRolSsWDF7b0pEJBudd56LR/7v/7y+Awfg+uuhd2+XEe2kLr7Yjcj40zvHxcHjj0OLFrB2bYqn1qpVi9mzZ/Paa69RsmTJUP+iRYto1KgRgwYN4pDvSyYRkfxGAY1kSNgITfwfoe1gQFOjRo2wxasiIvlR0aLw+uswaRIUL+71v/02NGgAS5em8SKjR8OsWVClite/dKm7yFNPpbhAJyIigrvuuotffvmFG264IdQfHx/P6NGjqVOnDp9++ik2hWKeIiJ5mf7SlAxJLmVzXKFC7An0abqZiBQkt94KP//sZosF/fqrG2R5/nmIj0/DRdq1c+mdb7vN6ztxwk1Na9zYVfpMQcWKFfnoo4+YMWMG1atXD/Xv2LGD66+/nq5du7I5TTmmRUTyDgU0kiH+gCaYsnlfsWIEvwNUQgARKWiqV4f58+HRR8EY1xcXBw8/DB07wvbtqZ8PuGI377wDX38NZ5/t9a9a5aKlIUPg2LEUT+/cuTNr1qxh2LBhREVFhfpnzJhB3bp1GT58ODH+nPsiInmYAhrJkPCkAIEaNL4pZhqhEZGCKDISnn4a5syBypW9/jlzXLmZNCUMAOja1a2f6dPH64uPhxEjXOaBH35I8dTo6GieeOIJ1qxZQ6dOnUL9x48f5/HHH6d+/frMnDnzlJ5LRCQ3UkAjGbLbS2wWCmg2+7710wiNiBRkrVvDypUuG1rQ33+fYsKAkiVhwgSYPdsN/wRt3Og+4N57U71QzZo1mTlzJh9++CEVKlTwnb6RLl26cPXVV2samojkaQpoJN3i4mDfPq9dLpC2eaOvIJyluCwSAAAgAElEQVRGaESkoDvjDPjwQ5g4MWnCgPr1Ye7cNF6obVs35WzQIFecE1zdmpdfhrp1IZXRFmMMN954I+vXr2fAgAFhyVqmTp3K+eefz6OPPsrhNEVYIiK5iwIaSbe9e93/S4PK8heglM0iIokZA716wYoV0KSJ179li4tT7rsPjh5Nw4WKFXOFNxcudIU5g7Ztgy5doEeP8LnAiZQsWZKxY8eybNkyWrVqFeqPjY3l2WefpVatWrz//vvKhiYieYoCGkm3v/7ytouaY5zGEcALaKpXr66UzSIiPjVquGUvjz0GhQp5/WPHuiUxixal8UJNm7q6NY89BoULe/0ffeQKdL7ySqop1S666CLmzp3LBx98QKVKlUL9f/75J7fccgutWrVixYoVp/h0IiI5Q39tSrr5A5py7CaQzCcU0Gj9jIhIUpGRrqTM4sXhgyybNkHLlvDQQ+CbuZuyIkXchZYtg2bNvP5//oH+/aF5czcklAJjDD169GDDhg0MHTqUIkWKhPYtWLCARo0a0adPH/bs2ZPiNUREcgMFNJJu/oQAZa3X8BfVFBGR5F18sYtFHnzQWxKTkAAjR0KjRvDTT2m8UP36sGCBSxxw+ule/9Kl7kPuvx8OHUrx9NNOO43hw4ezbt06unXrFuq31vL6669Ts2ZNXnjhBaV5FpFcSwGNpJt/hCa4fuZoRAQHAn0aoRERSV10tMvA/MMP4P+VuW6dG3R57DGIjU3DhSIiXGrn9evh5pu9/oQEGDPGDQV98UX4wsdEqlWrxhdffMHMmTOpXbt2qP/gwYMMHjyYOnXq8PHHH2t9jYjkOgpoJN0STzkD2BGsIodGaERE0qp5c/j5Z5ccIPhrND7e1bJp3NgNtqRJuXIweTJ8+61bsBP0xx9wzTVw9dXw22+pXuLSSy9l1apVjBo1ilKlSoX6t2zZwo033sgll1zC4sWLT/EJRUSyjgIaSbfkRmi2+BahaoRGRCTtihWD0aNd8c2qVb3+VavcaM3996exbg1Ax46wejU8/jhERXn9X33lRmseeyzVtGqRkZEMGjSIX3/9lXvvvZfCvsQDixYtonnz5vTo0YMtW7ac2kOKiGQBBTSSbmFraBKlbI6OjlbKZhGRdGjd2gUx99zj9QVnjtWtCzNmpPFC0dHw5JPuYm3bev0xMW7op3Zt+OSTVKehlSlThpdeeok1a9Zw1VVXhe376KOPqF27Ng899BAHDx5M+wOKiGQyBTSSbsmN0PgTAihls4hI+hQvDuPHu9Ga887z+rdtg8sug1tuCf8dnKpatWDWLDcVzf9F0/btcMMN0KEDrFlzkkvU4ssvv2T27Nk0bNgw1B8bG8vIkSOpUaMGY8eOVeIAEckR+otT0i25NTTKcCYiknnatIGVK+HRR8PLzbz/PtSpA+++m+oAi8cYlyxg/XqXFzoy0ts3e7YrgjNwIBw4kPI1gLZt27J06VLefffdsPo1e/fu5b777qNWrVr897//JT6VGjgiIplNAY2ki7WpTznT+hkRkcwRHe1miC1f7uppBu3fD7fdBpdeCps3p/FiJUrA88+7EZnLLvP64+PhpZfccNBbb7k5bimIiIigZ8+ebNy4keHDh3PaaaeF9m3bto1evXpx0UUX8dVXXykjmohkCwU0ki6HD4cXfksc0FSvXj37b0pEJB+rV8+Vm3npJTclLei779y+559PY4pncIHL11+7JAH+39d79sCdd0KTJjB3bqqXKFasGEOHDmXz5s3079+fSN+oT3DNTatWrfjhhx9O4SlFRE6dAhpJF/90M0MCZdgLKKAREclKhQrBvffC2rVw+eVe/7Fj8PDDrsbmd9+l8WLGwBVXuNGaZ55xadaCli1ziQS6dYONG1O9TLly5Rg3bhzr16/nZn8NHGDBggW0atWKK6+8ktWrV6fxxkRETo0CGkkXf0BThr0UIoG/gSOBvmrVquXEbYmIFAhVqrjBlQ8/hLJlvf4NG6BTJ7fW/48/0nix6Gh45BF3co8e4fu+/BIuuMCtr9m3L9XLVKtWjcmTJ7NixQou809nA6ZNm8aFF17Irbfeyq+//prGGxMRSRsFNJIuqa2fKVSoEJUrV87+mxIRKUCMgRtvdOv8+/cHf2LJTz5xWZlHjjyFaWhnnw0ffODmtTVr5vXHxbl5bjVqwIsvurTPqbjooouYPn06c+bMoZnvOtZaJk+eTO3atbnjjjtUw0ZEMo0CGkmX1DKcValSJWwutYiIZJ0zzoBx4+Cnn8LjkCNHXEKzCy+E778/hQu2aAELF7rhn3PP9foPHIB//9sV5vz005OmV2vTpg0LFy7kiy++oE6dOqH++Ph43nnnHc477zz69OnD77//fgo3JyKSlAIaSZfUatBoupmISPZr0MANrrz1FpQp4/WvXw8dO7rRnDRPQwsO//zyixvmKVXK2/fbb3D99dCqlfvAVC9j6NatG6tXr+bdd98N+/9DXFwcr7/+OjVq1KBfv37s2LHjFJ5WRMSjgEbSJbUpZwpoRERyRkQE3HGHWw7Tt6+LS4I+/thNQ3v6aZdEIE2io2HwYPj1VzevrVAhb9+CBdCyJVx5JaxaleplChUqRM+ePVm/fj1vvfUW5/pGfk6cOMH48eOpXr069913H7t27Ur7A4uIoIBG0im1KWcKaEREctaZZ8Irr8DSpeG1a44cgcceg1q13HKZNJeJKVPGzWtbs8YFMH7TprnCnP/6lxu9SUVkZCR33HEHGzZs4LXXXgtbbxkTE8PYsWOpVq0agwYN4s8//0zjzYlIQaeARtJFU85ERHK/Ro3ccpg33wyfhrZ9O9x8s1sus3jxKVywdm2YOhXmzAlfsGMtvPee23/vveHD+MmIiorirrvuYtOmTbz88stUqFAhtO/YsWOMHj2aqlWr0rdvX7Zt23YKNygiBZECGkmX5AKa4OxnBTQiIrlHRAT07g2bNsEDD4A/Z8vixdC8OdxyC5zS2vw2bVykNGWKSxIQdOIEvPwyVKsGQ4fCwYOpXqZIkSL069ePzZs3M2bMGMqVKxfaFxsby6uvvkqNGjXo3bs3mzZtOoUbFJGCRAGNpEtya2gU0IiI5F6nnw7/+Q+sWwfdu4fve/99Nw3tscfg8OE0XtAYuPpqt35m4kRXHCfo6FFXrLNaNZdU4MiRFC8DULRoUQYOHMhvv/3G2LFjqVSpUmhfXFwcb7/9NrVr1+aWW25h7dq1abxBESkoFNDIKYuLC6+vVo7dHACOAiVLluSMM87IqVsTEZGTqFEDPv8cZs92S1+Cjh93CQPOOw/eftv9rk+TQoWgVy/YuBHGjIGzzvL27d/vckdXrQovvHDSwKZYsWIMGDCAzZs3M2HChLDkAQkJCbz//vvUrVuX6667jhUrVqT9oUUkX1NAI6ds797wdln+ChudMf60OiIikiu1betq17z1FpQv7/Xv3OmmqNWv72aUpTlxQJEiMHAgbN4MTz4JxYt7+/bscdnSqlVzxTmPHv3/9u47TKry7v/4+7vLAkux0PtSpC5FOggIEkRFpdhiiUaNiYkao5g8+viLxpZoTGwp6hOjYi9YKBEQQlFQpPciUqWD9LqUvX9/3DO7s8vO7syyMDM7n9d1nWtm58w55z7H4zCfuVsRuyrH7bffzooVK3j99ddp2rRpnvUff/wxHTp04KKLLmLixIm4iAspIqWRAo1ELbS5WQUOUJGDam4mIpKAUlP9MM8rVsCDD/pMErRsmW+a1qMHfPllFDutXBkeftiPePa730GFCrnrtm3zk3M2bgzPPVfk+NFpaWncfPPNLFu2jPfee4/MzMw868ePH0+/fv3o1KkTH3zwAccirlYSkdJEgUaiVtCQzcHBNRVoREQST+XKvsvLt9/CTTflnb9m+nQ/BsCllxY53Uxe1av7/jNr1vgQk56eu27rVhg61AebF14oMtikpqZy7bXXsnDhQj799FM6duyYZ/3cuXO59tprad68OS+++CIHi6gBEpHSRYFGoqYRzkRESqeMDHjjDViwAC67LO+6MWNyp5tZuzaKndao4fvPrF4N997rJ+sM2rIF7rkHmjTxNTZF9LFJSUlh8ODBzJo1i4kTJ3LRRRflWb969WruvPNOMjIyePzxx9kR2uFTREotBRqJmgKNiEjp1qYNjB7tm5p17577enC6mWbNfHeZLVui2GmtWr7/zOrVfuPQ9m2bN/sam4wMePxx2LWr0F2ZGX379mXcuHHMnz+f66+/ntTU1Jz1P/zwAw8//DANGjTg17/+NStXroyioCKSaBRoJGqhfWiCTc4UaERESp9eveCrrwqebuZvf/Mtxn7727w/dBWpdm0/Gtrq1XD33XmDzY4dvv9NRgY88ECRE3QCtGvXjnfeeYeVK1dy9913UyGkz87Bgwf5xz/+QbNmzRg4cCCTJ0/WAAIipZACjUQtXA2NmZGRkRGbQomIyCkROt3Ma69B/fq56w4dgmee8aMy33//iaNgFqpOHd9/ZvVq38emYsXcdfv2wZ//DA0bwl13wbp1Re6uYcOGvPDCC3z//fc8+uijVKtWLWedc47Ro0fTt29fOnTowBtvvEFWVlYUhRWReKZAI1ErKNBsAurXr0/ZsmVjUygRETmlUlPhllv8iGjPPOO7xgQdPOj7/zdq5EdLi6rrSp06vo/NunXwhz9A6Fxmhw/DP//pJ8+55RZYvrzI3VWtWpWHH36YdevW8fLLL9OiRYs86+fPn8/NN99MRkYGjz32GNu3b4+isCISjxRoJGqhLQBqsI1jwFbU3ExEJBmUL++7u6xe7UNMSEUI+/fDk0/6YPP73/t5NSNWtSo88ogPNk8/DTVr5q47dgyGDYOWLWHgQN+5p4imYxUqVOD2229nyZIljBkzhv79++dZv3XrVv7whz9Qv359brvtNubPnx9FYUUknijQSNTyD9u8BchGgUZEJJlUrOinmVmzBp56yueRoH37/DDQjRrBQw9F2RStcmW/47Vr4cUXfX+aUKNH+3Gku3aFDz/0YacQKSkpXHLJJXz++ecsXryY2267jXIh/XaysrJ49dVXad++PT179uS9997jyJEjURRYRGJNgUai4tyJTc40IICISPKqVMn3n1mzxoeY0BZje/fCE0/4TDJ0KGzcGH4/JyhfHn71K/juO3jzzbyjEgDMmgU//jE0bepHKNi/v8hdZmZm8sorr7B+/Xoee+wxaobWAgFfffUV119/PQ0aNODhhx9mY1QFFpFYUaCRqOzf75s0BynQiIgI+IqVBx/0FSuPPw5nnZW77uBBP81Mo0bwi19AVKMop6X5yW8WL/aT4fTtm3f92rV+GOj69X0BNm8ucpfVq1fnoYceYt26dbz55pt06dIlz/qtW7fy+OOPk5GRwdVXX80XX3yh0dFE4pgCjUQltO+kkU1VduQEmkaNGsWkTCIiEj/OOMP3n1m71tfYhPaxOXoUXnkFmjeH66/3I6dFzAwuuQQmToQ5c/wOQuaeYfdu34EnIwNuuglmzy5yl+XKlePGG29kxowZzJw5k5/+9Kd5mqMdP36cjz76iD59+tC2bVtefPFF9uzZE0WhReR0UKCRqIQGmirsJJVsNgX+btiwYSyKJCIicejMM32Fybp1fnTmevVy12Vnw3vvQbt2vo//9OlR7rxDB3jnndwhnytXzl139Ci89RZ07gznnQfvv+9fK0Lnzp0ZNmwYGzZs4KmnnjphGoLFixdz5513UqdOHX72s58xY8YM1dqIxAkFGolKaMfOavg/NuJ/5aoROoaniIgIUKGCnz9z1Sp49VXf5SXU6NE+d/Tq5SfwPH48ip03aOCHfF6/Hv7yF6hbN+/66dPhuuv8fDZPPBHRDKDVqlXj/vvvZ9WqVYwcOfKE0dEOHjzIa6+9Rrdu3Tj33HNVayMSBxRoJCqhgaY6vrpmI5CRkUFKim4nEREpWNmycOutsGwZfPCBr50JNW0aDBkCLVr4qWcOHIhi52eeCb/9rR+Z4P33fUIKtWmTH26tfn24+WbfZK0IqampDBw4kM8//5zly5czdOhQqlSpkuc9Cxcu5M4776R27drceuutfPPNN6q1EYmBhP4GambpZvaYma0ws8NmtsnMXjOzukVvfcK+zjazF8xsnZllBR6fN7OzCnhvmpn1N7N/mNliMztoZofMbJmZ/dXMqpfMGcaf0CZnoTU0am4mIiKRSE2Fa66BefPgs8+gR4+861euhLvu8pUvv/89bNkSxc7T0vzIZ1995UdBu+kmn6SCjhyBN96ATp186Hnrrbwj3YTRvHlznnnmGTZu3Mg777xD796986w/dOgQr7/+Ot27d6ddu3Y8//zzmrBT5DRK2EBjZuWBScBDQCVgJLAeuAWYZ2YRD7llZtWAmcDdwDFgBLAP+A0ww8yq5NukN/A5cCdQERgLTACqAfcBC82sebFPLo6Fa3KWv62xiIhIYcxgwABfM/PNN3D11RBa0b9zpx9UICPD1+wsXhzlATp18uHl++/hscegdu2866dP94Gnbl1fu/Pdd0Xusnz58lx//fVMmTKFZcuWMXToUKqGTsADLFq0iHvvvZc6deowZMgQRo0axdEI+vCISPElbKABfg90A6YDzZxzP3bOdcUHiurAa1Hs63ngHOAToHlgX62BvwPNgGfzvT8b+BDo6pxr5Jy70jk3MLCPz4FawOvFP7X4lb/J2T5gP6qhERGR4gvOkblype9vU7Fi7rojR+D116FNG+jf3/e5iaqfTc2avrnZ2rV+JILu3fOu37kTnnkGmjWDfv3go48iGkSgRYsWPPPMM2zYsIF3332XPn365Fl/7NgxRowYwaBBg6hXrx733Xcfi6NOZSISCUvEtp5mVhbYBpwJdHDOzcu3fgHQFujknCu0oayZ1QY24GtmGjjntoasK4ev9akC1HHOFdmb0MzqQM5Ixg2dc+siPrET97WkVatWrZYsWVLcXZS4QYNg1Cj//BmGMoDnaAm8/fbb3HDDDTEtm4iIlA67dsG//uXny9y06cT1jRrBHXf4mpsq+dtQRGLOHHj5ZXj3XT9JTn61asFtt8HPf+7bvkXo22+/ZdiwYbz55ptsKqjgQMeOHbnlllu47rrrTuiTI5JMMjMzWbp06VLnXObJ7itRa2h64MPMqvxhJuCjwOPlEezrYvx1mBoaZgCcc1nAaCAVGBBJwZxzm4Bgw9k6kWyTSPI3OctJbqqhERGREnL22XD//b6P/5tvnjiAwJo18Lvf+aGgf/7zKOezAejY0U+Is3Ej/P3vkJnv+9SWLX5UtEaNfLu4jz/2VUVFaN68OU8++STff/89Y8eO5ZprrqFsaB8eYM6cOdx1113Url2ba665hrFjx3Ls2LEoT0BEQiVqoAl+tM0Nsz74etvTvC8CgwicHfgzmq6MCSG0j2N1tucEGvWhERGRkla2LNx4ox9AYNIkuOKKvP1sDh2Cf//bB57zz4fhwyNqLZbrrLP8CASLFsHUqXDDDXkHEcjOhrFj4aqrfF+boUMhglYTqampXHzxxXzwwQds3ryZf/7zn3Tq1CnPe44cOcLw4cMZMGAA9evX595772X27NkaJU2kGBI10ATrfzeEWR98PZJv2SW5L/ADBZQBFjnn1kSygZktKWgBmkR4zNOmoBqatLQ0aufvbCkiIlJCzOCCC3xFyZo18L//C9Wq5X3P1Kl+9LRGjeCRR/zUNFEdoGdPePtt2LABnn4amuT7J/iHH+C556B1a9/p51//gr17i9x1lSpVuOOOO5g1axaLFi3ivvvuO2Heti1btvD888/TuXNnWrRowaOPPsp3EQxSICJeogaaSoHHAhq+AhAcvb5ymPWnZF9m1h4/WAHA/REcO6EcO+bbNQcFa2gaNGhAampqzMolIiLJo0ED+NOffGAZNsy3Hgu1cSM8+qifS/Pyy/0gAlG16Kpe3bdnW7ECJk70tTbly+d9z8yZcPvtvq/NT38KX37pa3OK0Lp1a/7617+yYcMGRo4cyZAhQyhTpkye96xYsYJHHnmEZs2a0aVLF1544QW2RDV2tUjySdRAE3fMrCZ+lLTywPPOubGRbuucyyxoAVadqvIWx44def+uxg9sQs3NRETk9Ctf3meJWbP8CMw33OCnoQnKzob//AcGDixmrU1KCvTt62ttNm3ys3126JD3PYcO+U4+vXv7Gp2HH45o+Oe0tDQGDhzIJ598wubNm3nppZfo2bPnCe+bNWsW99xzD3Xr1qV///4MGzaMvRHUCokkm0QNNPsDjxXCrA8O+LjvdOzLzCoDY4CGwHD80NGlTmhzs3IcpiIHNKmmiIjElBl06+Zzx/r18OST0DjfTHQbNpxkrc3ZZ/th1ebM8R16fv1r/1qotWvh8cf98M/du8NLL/khoYtQrVo1fvnLXzJ16lTWrl3Lk08+SevWrfO8Jzs7mwkTJnDLLbdQo0YNrr76aj7++GMOFjRCm0gSStRA833gsV6Y9cHXIxky+aT2FZjgcxTQARgP/MQ5V3S9cwLKPyCAoUk1RUQkftSsCQ884CtJxo/3fflDW3SF1to0bAgPPuhblkXl3HNzx5N+/30/OU5Kvq9T33zjA1CtWn4kgxEjIholLSMjgwceeIBFixaxcOFCHnjgARrkGzY6KyuLjz76iKuuuooaNWrw4x//mOHDh3PgwIEwexUp/RI10CwIPHYIsz74eiQDORZ7X2ZWBvgA6AN8DVzhnCv6EytB5R8QIBvYimpoREQkvqSkwIUX+lHPwtXabNzoX2/eHHr08KM4R9Waq3x5+PGP4fPP/UH+8hc/YECoo0fh009hyBCoXRt+9Sv44ouI+tu0adOGJ598kjVr1jB16lR++ctfnjBvzYEDB/jwww+55pprcmpuPvzwQ/bv3x9mryKlU2mYWLO9c25+vvXFnVizfujkmYVNrGlmBrwJ/ASYD1zgnNtdAqcXWra4mljz5Zf9ZzFAPybwNv2pBUyZMoXevXvHtGwiIiKFyc72ffz/9S9fYVJQk7P0dF+hcsstflS1/BUvRXIOFiyAt96Cd96BrVsLfl/duj4MXXstdOrk281F4MiRI4wfP57333+fUaNGsW9fwa3h09PTueSSS7j66qu57LLLqFSpUoHvE4mlpJ9YM1AL8o/An/80s2A/F8xsKD7MfBEaZszsLjNbbmZP5tvXZuA9oCzwYqDWJehpoDrwdmiYCXgeH2aWA/1LOszEo/xNzoJjrqiGRkRE4l1orc3GjfDss9CmTd73HDrkc0i/fn4ggQj7+Ocy803SnnnGd9wZOxauu84npVDBAnTp4vvcPPRQRPPblC1blssuu4y3336bbdu2MXLkSG688UbOOOOMfOdxiE8++YTrrruO6tWrM2TIEN5991327NkTxcmIJI6ErKGBnL4rU4CuwGZgKn6umK7AdqCbc251yPsfAf4AvOGcuznfvqoB3+DnfVkFzAYygdbAd4F97Qx5/yBgRODPCcCmMMV8yjm3/CTOMa5qaH7zG99sGODX/I0B/IbLUlM5fPjwCcNOioiIxDvnfB//YcN8kAnXh79zZ/jJT3ylSs2axTjQ3r2+6dn778OECXD8eMHva9PG19pcdZUPOhHKyspiwoQJfPTRR4wYMSJscElLS6NPnz4MHjyYgQMHUq9euO7DIqdeSdbQJGygATCzdOB/geuB+sBOYBzwkHNuQ773PkKYQBNYXwV4BBgM1MR3D/kU+EP+2hczuxl4PYIiXuCcmxLFKeUvU1wFmhtugHff9c8f4yHq8QSPZmSwdu3amJZLRETkZGVl+QEDXn/dV6wU1M0lNdXX3txwAwweDJUjme0uv+3b4aOP4L33/Gyg4bRu7YPNlVdCZmZUzdL++9//Mnz4cEaMGMHu3eEbkHTs2JHBgwczaNAgWrdujUV4DJGSoECTJOIt0PTv739YAniRX7GLlxnfuzdTpkyJablERERK0ubNfhjot96CRYsKfk96Ogwa5MPNRRflnQMnYuvXw4cf+nAzp5Auv82b+2Bz5ZXQvn1U4WbSpEkMHz6ckSNHsiP/hHIhGjduzKBBgxg0aBA9evRQyws55RRokkS8BZoOHXzVPMCHXM1UPmLPTTfxxhtvxLZgIiIip8iiRb452rvvhp+Ys0oVP5jANdf4wQSKlQVWrIAPPvC1NwsLGaS1UaPccNOlS8QjFxw/fpyvv/6aESNGMHLkSFatCj93d9WqVbnssssYNGgQ/fr1o3KxqqJECqdAkyTiLdDUr+/7OAJMpg//5AtaPfwwjz76aGwLJiIicoplZ/sWYu+84wcWCNeSq2rV3HDTp08xw81338HHH/tl9uzw76tXzx/siiv82NMRHsw5x9KlS3PCzaxZs8K+Ny0tjfPPP58BAwZw6aWX0qxZMzVNkxKhQJMk4inQOAcVKsDhw/7vRbTmlyzh1ldf5dZbb41t4URERE6jrCzfz+add2D0aP93QapV8xUpV18NvXsXM9ysXQuffOLDzddfh3/f2WfDpZfC5ZfDxRdDvpHPCrNx40ZGjRrFyJEjmTRpEkePHg373saNGzNgwAAGDBhAnz59SM8/gptIhBRokkQ8BZoDByB0GPvN1KIHW3ll4kT69u0bu4KJiIjE0N69PtR8+CGMGwdHwkyvXb26H0hgyBDo2xfKlSvGwTZu9KOlffwxfPll+Ak609J89dDll/sliukV9u7dy9ixYxk5ciTjxo1j165dYd+bnp5O3759cwKOpnGQaCjQJIl4CjRr1/pmu0FHSOMsjrFo1Soa559+WUREJAnt2ZM33ISr6DjjDF+ZMmQIXHJJ3h8MI7Ztm58h9JNPYNKk8AcDaNvWB5uBA/1EnhH2uzl27BgzZsxgzJgxjBkzhvnz5xf6/pYtW3LJJZfQv39/evXqRYUKFaI5I0kyCjRJIp4CzezZfhx+gLPYxTqqcJYZhw8fpmzZsrEtnIiISJzZvRtGjQ9R5OUAACAASURBVPLhZvz48HmjXDk/4eeQIT5vVKtWjIPt2+cPMmoUfPYZFDKaGbVqwWWX+UTVt2/UTdPGjh3LmDFjmDBhAvv37w/73nLlytGzZ08uvPBC+vfvT7t27UiJMEhJclCgSRLxFGjGjfO/IgE0ZQWjac4FtWuzaVO4OUVFREQEYNcuX3Pz6afw+edw6FDB70tJgV69fGXKpZf60Zqj7n9//DhMn+7DzejRsLyQ+b3LlIGePX2fm4sv9jU5UQwJPW3aND777DPGjBnD8sKOA1SvXp1+/fpx4YUXcuGFF2pST1GgSRbxFGjeegtuusk/787X/Ike3N+lCzNmzIhtwURERBLIwYM+1Hz6qc8bhcx7SZMmuZUp559fzH43333nDzRqFEyb5gNPOLVr+0l1Lr7YVxtVqRLxYVavXs3YsWMZP348kydPZt++fYW+v2XLljnhpnfv3hoaOgkp0CSJeAo0zz0HQ4f655cziusZxMdXXcXw4cNjWzAREZEEdfQofPGFDzcjRkBhjR4qVfITXF92GQwYADVrFuOAO3f64dn+8x/fRG3nzvDvTUmBrl1za286doTU1IgOc/ToUWbOnMn48eOZMGECM2bMIDvcAAZAamoqnTt35oILLuCCCy6gR48e6n+TBBRokkQ8BZoHH4Qnn/TPb+VVWnMb6++9l2effTa2BRMRESkFsrNh1iyfNf7zHyii/z2dO+fW3rRvH3E//1zHj/sOsuPG+WXGDD9HQzhVq/o+Nz/6EfTrB40bR9w8bffu3UyePDkn4BQ2qSf4uW+6du2aE3C6d+9O+fLlozk7SQAKNEkingLNL34Br7zin/8Pf8bxAHWee4577rkntgUTEREphTZsgDFjfLj573/D97sBX1vTv79fLrywmLU3O3bAhAm5AWfr1sLfn5Hhw01wieKgq1evZsKECYwfP55Jkyaxu7B2d/gBBrp3754TcLp06UK5YrW/k3iiQJMk4inQXHGFrxIHeJrfsZC/Mvijj7jyyitjWzAREZFS7tAhmDIlt/bm++8Lf3+7drkBp2dPiLpyIzsbFizIDTdffVV43xuA1q1zw03v3hGPnnb8+HEWLFjA5MmTmTx5Ml9++WWR/W/S09Pp2rUrvXr1olevXnTv3p1KxRr7WmJJgSZJxFOg6dXL9yUEeJ2beYc3+NPMmXQOjuUsIiIip5xzsHixH535P//xA5oV0j2F9HSfL4IBp1WrYoyctmePT1QTJ/pl6dLC35+aCl26+CZqvXvDeedBxYoRHerYsWPMnTs3J+BMmzaNAwcOFHG4VNq3b58TcHr27En16tUjPDmJFQWaJBFPgaZlS8fy5f4TcDSX8QCf8d/Nm6lVq1aMSyYiIpK8du70GWP8eD962vr1hb+/Th0fbH70I7jgAqhbtxgH3bTJT+YZDDhFHbRMGT+h5/nn+4DToweceWZEhzp69CizZs3KCThfffUVhw8fLnK7Fi1a5AScXr16kZGRgUWd5ORUUqBJEvEUaKpVzWbHTt/jcDrdGJI2l42HD2uSLBERkTjhHKxYkRtuJk/2w0QXpmlTH2z69PGPUf9O6ZwfGjoYbiZN8hPvFCYlBc4914eb3r19u7iqVSM6XFZWFrNmzWLq1KlMnTqVr776ir179xa5Xd26denevXvO0qFDB/XDiTEFmiQRL4Hm+HFIS3M453/ZWEoTLmvkWLV6dUzLJSIiIuFlZfkmaePH+2XOnKK3adHCB5tgyIm65dbx436ItokT/ZjU06ZBBIGDNm18Dc755/smahFOvHn8+HEWL16cE3CmTp3K5s2bi9yubNmydOjQIU/I0WSfp5cCTZKIl0Dzww95P9AWcSZ39W7PlClTYlYmERERic727T5nTJjga2/WrCl6m9atcwNOz54nEXC+/NIHnC+/LLoGB6BBAx9sgku7dr7pWhGcc6xevZpp06blBJwVK1ZEVNR69erRrVs31eKcJgo0SSJeAs2yZb4TIUAaR/iScrx44428+eabMS2XiIiIFN+6dT7YBJeiusKAr8Hp2TN3iWI6Gi87G5Ys8eEmGHC2bSt6uwoV/ESf553n++B06wZnnx3RIbdu3cr06dNzltmzZ3OosHGwA9LS0mjbti2dO3emU6dOdO7cmVatWlEmgmAlRVOgSRLxEmimTvU1wAC12MxL1GHWgw/yxz/+MablEhERkZLhHKxenTfgRNByi1q18gacCCtS8h74229za3C+/hrWro1s28xMH3C6dfOjqrVs6UdYK8LRo0dZsGBBnpCzNsJjVqhQgfbt29O5c+ec5ZxzztGAA8WgQJMk4iXQjBgBQ4b4561ZxB20JeXll7n99ttjWi4RERE5NYJ9/YPh5osvYMuWorerWBG6d/fhpkcPnzMinJIm16ZNPtgEl7lz4ejRorerVMmPptali1+6dvXDuEUQNrZs2cI333yTE3BmzZoV0WhqAGeddRYdO3bME3Lq1aunkFMEBZokES+B5tVX4bbb/PPeTKE3F9D1s88YMGBATMslIiIip4dzvs/NtGm+5ca0abB8edHbmflm6926+aVrV/93BBUpuQ4dgtmz84acH36IbNvatXPDTZcuPvBEMGT00aNHWbJkCbNmzcpZFi1axPGiJhgNqFmzJh07dqR9+/ace+65tG/fnsaNGyvkhFCgSRLxEmiefhruv98/v4KPqcpV3L1oEa1bt45puURERCR2tm/32WLaNL/Mng3HjhW9XaVKPlsEA07XrlCzZhQHDlYfBcPNzJl+ttFIwoYZNG8OHTtChw5+ad8+opBz6NAh5s+fnxNwZs+ezfJIUl3AGWecwbnnnpsTcNq3b0+rVq1IS0uLeB+liQJNkoiXQPPAA/DnP/vnt/EKW/gFb+/ezZkRToolIiIipd/BgzBrVm4NzowZsHt3ZNs2apS3EqV9ex98InbgAMyb5w86c6ZfIu2LA3DOObkBJ7hEMDfO3r17mTNnTp6anHXr1kV82LJly5KZmZkTcNq3b0/btm2pXLly5GVPUAo0SSJeAs3Pfw7//rd/fj9P8U3FJ5iyf39MyyQiIiLxLTvbT/T5zTc+Z3zzDSxc6F8vipkfUa1TJ1+Z0qmTn4uzYsUoCrB1q09YM2fmBp1IExZARsaJISeCmUe3bdvG7NmzmT9/PvPmzWPevHmsWrUq4sOaGeeccw7t2rWjdevWtGnThjZt2tC4cWNSo2qrF98UaJJEvASaKwYd49NRfsiSp/kd45qNYuK338a0TCIiIpJ4DhzwE3x+803uEsloagApKX4gs9CQ066dH9E5IsGmanPm+GXuXL/s2RP5CdSu7YNNu3bQtq1fmjYtcmi3vXv3smDBgpyAM3/+fJYsWcLRSAY7CEhPTyczMzNPyGnTpg01a9ZMyL45CjRJIl4CTZ9uh/hiRjoA/+JnjOy/gf98/nlMyyQiIiKJzznYsCG3Fmf2bJ81Im0Ikprqu8Sce67PGMHHiPvkBEc7CIabuXN9ASIddACgXDk/0kEw4ASXGjUK3SwrK4ulS5fmqcmZP38++6NsBVOtWrUTQk5mZmbcN1tToEkS8RJo2jQ5wOLVvo73dQYz/Rc1+b//+7+YlklERERKp+xsX5Eye3ZuwJk719fuRKpWrbwBp107aNYswjlygikrNOTMneuHk45GzZonhpyWLX0ACiM7O5vVq1czb948Fi9ezKJFi1i0aBGrVq0i2u/sGRkZtGrVKs/SsmXLuOkDrUCTJOIl0NSpcojNu3wNzev0YsPjF/H73/8+pmUSERGR5HH8uO+PExpy5s3zAxFEqnx5aN06b21O27ZRzJOzZYsPNgsW+M5ACxf6SUEjHMoZyK1Sysz0tTrBpWnTQoPOgQMHWLp0KYsWLcoTdLZu3Rr5sQPq1q17QtDp0aPHaW+2pkCTJOIh0DgH6WWPkXXM/6TxLzIp98b93HTTTTErk4iIiMjx434unAUL/DJ/vn+M9jt+gwY+X2Rm+sCTmekrUiIagODwYVi2LDfgLFzoC7F9e3SFSE31I62FhpxWrXz4SU8Pu9n27dtzwk0w6CxevJgDUVRn1a5dm03R1j6VAAWaJBEPgebAgbzDJv6dWmROeo8LLrggZmUSERERCWfLlhNDTrQVKWbQsGFu0AkuLVsWmi9ybd2aG3AWLfKPS5bAkSPRnYwZNG6cN+S0bOmDTpiqpezsbNauXcuyZctYunRpnqWg/jk/+tGP+O9//xtduUqAAk2SiIdA8/33ftTCoD9TlitXLqNJkyYxK5OIiIhINA4d8nkif9DZuze6/ZhBkyZ5Q06LFr5/TpHz5hw96tvNLVkCS5fmLitW+HXRqlXLH7h5c78EnzdqBAVM1umcY8OGDSeEnD59+vDHP/4x+uOfJAWaJBEPgWb+fD+5FUAl9nEHZ/DY4cOUK6Sdp4iIiEi8cw7WrfP5InRZutQHoGjVq+fzRIsWuRmjRQv/ekpKIRsePQorV+YNOUuX+vZ00dbogB/5oHHjgsNOzZo+lcUBBZokEQ+BZuJE6NfPP89gLQMrt+Vv0f6cISIiIpIgsrNh7VpYvDhv0Fm2DLKyot9ferrPE/mDTpG1OseO+SGlly7NTVpLlvganWhGQwh1xhl+AIKmTX2fnaZNfUG6dSve/k5CSQaaSAavkyS2Y0fu86rsiGJgdxEREZHEk5LiKzgaN4aBA3NfP34cVq06sTZnxQo/NkA4hw7lNnXLr25dnymaNPH5Irg0aQKVK5fJDR+DBuVulJ3th5D+9lu/rFiR+7h2rV8fzt69uROLBtWrB+vXR3x94pECjRRq587c51XYSWqdOrErjIiIiEiMpKb6yoxmzWDIkNzXs7N9n+Ngvvj2W99a7NtvYePGwve5caNfpkw5cV2NGnlDTm7YSaFKvXo+iPzoR3k3ysryqSs06ASfh5sstGnTqK5DPFKgkULt2HIEKAv4Gpr0hg1jWh4RERGReJKS4kdEa9gQLroo77p9+/LmimDQWbGi6H4627b55euvT1x39tm5NTmNGuUuDRuWo8E5rSjbqtWJG+3c6Q++cqWfuTT4GOwsncAUaKRQOzYcIjTQnHHOObEtkIiIiEiCqFwZOnb0S6jsbN/Ka8UKX6GycmXusmpV4U3YAHbtglmz/JJfSopvytawYWjQgUaNqtCoUXfqdulOampJnWF8UKCRQu3cnDu6Rjo7qRM6hrOIiIiIRC0lxU+LkZEBF16Yd112NmzenDfkhC4FTCVzwvbr1/tl6tQT15cp4ycTDQadFi3gt78tsVOLCQUaKdSO7bmzUJVhB3XrnhfD0oiIiIiUbsEalrp1oXfvvOucg+3b89bmrFnjxwJYs8b3xylqAONjx2D1ar+AH3VNgUZKtdBRzowd1K1bN3aFEREREUliZn6wgBo14LwCfmPOyvI1M2vW5A06wWXbthO3adTolBf7lFOgkUL9sDt3Jqjj7FSgEREREYlT5crljoZWkIMHc0NO8LF589NZwlNDgUYKteNA+dw/yuylcuXKsSuMiIiIiBRbhQrQqpVfSpOUot8iySo7G3ZnVcz5u2zlI4W8W0RERETk9FOgkbD27oVscsf1K1dNt4uIiIiIxBd9Q5Ww8g4IkE2l+mpuJiIiIiLxRYFGwgoNNGezi/SMerErjIiIiIhIARRoJKydW4/mPK/CTio1aRLD0oiIiIiInEiBRsLasXZvzvOq7KBKy5YxLI2IiIiIyIkUaCSsH9YdyHl+JjuoXRpmXhIRERGRUkWBRsLaum5fzvPy7KRePfWhEREREZH4okAjYW1efzDneTl2Uq1atRiWRkRERETkRAo0EtYP24/nPC+ftgczi2FpREREREROpEAjYe3akzupZvnyh2JYEhERERGRginQSFi7D6bnPK9w5rEYlkREREREpGAKNBLWniOVc56fUV23ioiIiIjEn4T+lmpm6Wb2mJmtMLPDZrbJzF4zs7rF2NfZZvaCma0zs6zA4/NmdlYh26Sa2b1mtsjMDpnZdjP70MxKxYQte46fmfP8rLrphbxTRERERCQ2EjbQmFl5YBLwEFAJGAmsB24B5plZ4yj2VQ2YCdwNHANGAPuA3wAzzKxKAdukAMOBZ4F6wGfAEuAqYLaZdSn2ycWBY8dgH7lZrlazqjEsjYiIiIhIwRI20AC/B7oB04FmzrkfO+e6AvcB1YHXotjX88A5wCdA88C+WgN/B5rhQ0t+twJDgO+AFs65q5xzfYCrgQrAO2ZWplhnFgd2bsvbZ6Ze26grvURERERETrmEDDRmVha4K/Dnnc65/cF1zrlngYVAbzPrGMG+agPXAUeAO5xzod/kfwdsB35iZjXybTo08Pg/zrmtIcf/GBiFD0iDojqxOLJ5Wc4pUYajNOxcKlrRiYiIiEgpk5CBBugBnAmscs7NK2D9R4HHyyPY18X46zA1NJgAOOeygNFAKjAg+LqZNQJaAofwTc1O5vhxafXcdTnPq7CD2k3PiWFpREREREQKlqiBpl3gcW6Y9cHX256ifQW3WeycO3qSx49L65dty3l+JjtIS0uLYWlERERERAqWqH08GgQeN4RZH3w94xTtqySPj5ktCbOqSSTbnwrb1u7NeV4pZU+siiEiIiIiUqhEraGpFHg8GGb9gcBj5TDrT3ZfJXn8uLRzS27FU+W0vYW8U0REREQkdhK1hqZUcc5lFvR6oOam1WkuDgCHDlrO88rpBwp5p4iIiIhI7CRqoAmOalYhzPqKgcd9p2hfJXn8uPT62lt5aX8WK2d+h6W0iHVxREREREQKlKiB5vvAY70w64Ovrwuz/mT3VZLHj1vlK5Wjdd/WsS6GiIiIiEhYidqHZkHgsUOY9cHXF56ifQW3aW1mBQ3/Fc3xRURERESkmBI10HwF7AGamNm5Bay/KvA4OoJ9jQOygV75J880s3L4uWSOA2OCrzvn1gDLgHTg0pM8voiIiIiIFFNCBhrn3BHgH4E//2lmwT4rmNlQ/PwvXzjn5oS8fpeZLTezJ/PtazPwHlAWeNHMQpvhPQ1UB952zm0jr2eD7wkNQmZ2BTAQWAmMPInTFBERERGRIiRqHxqAJ4B+wHnAd2Y2FT/vS1dgO3BrvvdXA5oDtQvY1z1AN+BKYLmZzQYygdbAd8DQArZ5DRgADAlsMzFwjN7AIeAnzrljJ3OCIiIiIiJSuISsoQFwzh0GLgAex88HMxgfaIYBHZxzq6PY1w9AF+Dv+JqaIcCZwN+ALs65nQVskw1cDdwHbAIuA9oAHwOdnHMzintuIiIiIiISGXPOxboMEoaZLWnVqlWrJUuWxLooIiIiIiIlJjMzk6VLly4NNx9jNBK2hkZERERERESBRkREREREEpYCjYiIiIiIJCwFGhERERERSVgKNCIiIiIikrAUaEREREREJGEp0IiIiIiISMJSoBERERERkYSlQCMiIiIiIglLgUZERERERBKWOediXQYJw8z2litXrnKTJk1iXRQRERERkRKzatUqsrKy9jnnzjjZfSnQxDEz2wJUANbHqAjBJLUqRsdPdLp+xadrd3J0/U6Orl/x6dqdHF2/k6PrV3yxuHb1gYPOuVonuyMFGgnLzJYAOOcyY12WRKTrV3y6didH1+/k6PoVn67dydH1Ozm6fsWX6NdOfWhERERERCRhKdCIiIiIiEjCUqAREREREZGEpUAjIiIiIiIJS4FGREREREQSlkY5ExERERGRhKUaGhERERERSVgKNCIiIiIikrAUaEREREREJGEp0IiIiIiISMJSoBERERERkYSlQCMiIiIiIglLgUZERERERBKWAk0SMbN0M3vMzFaY2WEz22Rmr5lZ3WLs62wze8HM1plZVuDxeTM761SUPZbMrIKZDTazV83s28C1O2BmC8zsYTOrFOX+1pqZK2RpcarOJVbMbEoR53xxlPtLivvPzPoUcd2Cy8MR7q9U3ntm1tHMHjCzT8xsQ/B8ItjuZjObaWb7zWynmY0xs/OKWYZUM7vXzBaZ2SEz225mH5pZy+Ls73SK5vqZWYqZ9TKzp81sjpntC/w/uMrMXjazRsU4/rAi7stfnvxZnhrR3ntm9kgR5/pUMcqQFPde4P2RfB5OiuL4iXzvFfu7SWn87CtzOg8msWNm5YFJQDdgMzASaAjcAlxmZt2cc6sj3Fc1YDpwDrAaGAFkAr8BLjGz7s65nSV+ErFzPfBK4PkyYBRwBnAe8ChwnZn1ds5ti3K/b4R5fU+xSpkYPgb2F/D6xkh3kGT33xbC3yepwE8Cz6dGud/Sdu89BAyKZgMzex5/zxwCxgPlgQuB/mZ2lXNuRBT7SgGGA0OA3cBnQDXgKuBSM7vAOTczmvKdZtFcv8bAl4HnW/D/rhwHugC3A9eb2QDn3LRilOPzwD7z+7YY+zpdor73Ar4CVhbw+pxodpJk9x6E/+wCuBR/7tF+HkJi3nvF+m5Saj/7nHNakmABngAc8DVQKeT1oYHXp0Sxr7cD23wMlAl5/W+B14fF+nxL+Nr9FPg/oGW+12sDcwPn/G4U+1vr/9eL/bmdxms4JXCdGpbAvpLq/ivkOlwSON/vAYtwm1J57wH3A48BlwO1gMOFnSfQL3DtfgCahrzeHcgCdgFnRXH82wL7WwHUDHn9ysDr34Xeq/G2RHP9gCb4L0F9Q+87oBzweuB81wFpURx/WGC7PrG+Fqfh3nskcK43l9Dxk+beK2I/ZwW3Df1/OoLtEvnei/q7SWn+7Iv5fxAtp+E/MpTFJ2cHtC9g/YLAuo4R7Ks2/te4rNCbN7CuHLANOAbUiPV5n6Zr2z1w7Q4DZSPcZm1xPrATeaGEAo3uvzzn+07gmj4ZxTZJce9F8KVyTODa3VPAuhcC6+6L4nhLA9sMLmDdyMC6K2N9XUrq+hWyXXrIvzW9o9guYb9URnvtKPlAo3vPb/fzwLlOj3K7UnPv5TuvAr+blObPPvWhSQ49gDOBVc65eQWs/yjweHkE+7oY3/dqqnNua+gK51wWMBrfFGZA8YubUBYEHssBVWNZkCSh+w8ws4rkNtN4K5ZlSTRmlo6vXYDcz75Q0XweEugz0hLffOOzk91fInPOHcL/UgtQJ5ZlSQa69/IINr/V56F3wneT0v7Zpz40yaFd4HFumPXB19uW0L5ujXBfpUHjwONRIKp+G2b2O3zzjSxgCfCpc257yRYv7vzMzKoC2fgvPiOcc99Hsb3uP+8KoCIwzzm3NNqNk/TeC2qO/0d+u3NuQwHro/k8hNx7crFz7mgJ7C9hBdrTZwT+LKg/QlGuMLMr8T9KrAFGO+eWl1T54kxfMzsX339hAzDWORdV/xl07wFgZg2AXvh/hz8o5m5K271X0HeTUv3Zp0CTHBoEHgu6gUNfzwiz/lTtqzT4TeBxXKCGIBpP5/v7OTP7tXPutRIoV7z6fb6//2pmjzvnHo9we91/3sn+GpmM915QofeQc+6Ame0Gzjazys65fSezP5LnngS4DqgBbMf314zWr/P9/Wczewn4jXPu2MkWLs7cmO/vx83sY3xTtIIGTimI7j3vBsDwoXBHMfdR2u69gr6blOrPPjU5Sw7BofsOhll/IPBY+TTvK6GZ2QDgZ/hfQB6KYtNR+F/YM4AKQGvgWfwvJ/82s+KMmBPvvsT/A94Ef87Ngf+H7+/ymJn9ppBtQyX9/WdmtYEf4fsSvRfl5sl47+VX1D0E+kyMmpnVB54P/PlwlD/wzAN+CTTD35eNgTvx/XHuAP5SgkWNtZXAb/EjM1YC6uO/kG/Ed6SO5kcK3XveyfzAU+ruvUK+m5Tqzz7V0IgUg/n5Ot7G/yr0O+fcgiI2yeGcuzvfS0uA+8xsOfAv4M/4znSlhnMu/zwpK4A/mdls/HCZj5jZvwJt8KVw1+GbRYxzzkXVrCcZ7z059QJ9uj7BD9c6wjn3cjTbO+deyPfSGuBFM/sC32zlLjN71jm3vkQKHEPOubfzvXQAeNfMJgOLgMGBaRS+Of2lSzxm1gFohQ8go6PdvrTdeyfz3STRqYYmOQSrryuEWV8x8FhU9WJJ7yshmZ+IdBxwNvBsAR+IxfUqfpSu5mbWsIT2Gdecc+OB2fghN7tGsEnS33+cms6vyXTvFXUPgT4TI2Zmafh5KDoB0/BzY5QI59wSfK1iGXytZKnlnNuMH/Ya/OAnkUjqey8g+Hk4vBjNvsNKxHsvgu8mpfqzT4EmOQQ7XdcLsz74+rrTvK+EY2ZV8HMwZOD/8fltSe3bOZcNrAr8Wbuk9psAvgs8RnLOyX7/tQTa4/8hiXjys6Ik2b1X6D0UqG04C9gVQRvyIvdHKb4nA4MAvIGfE2k+cPkpqGWN5vMh0UV7rkl774GfoR64NvBn/pqvkpAw916E301K9WefAk1yCFY5dgizPvj6wtO8r4RiZpWAsfjq7U+An7vAYOsl6OzA44FC31W6RHPOSXv/BQQ7En/inCusHXRxJMu99y1+dLfqgV8084v2Hgrek60DtRUnu79E8nd8E8gVwEXOud2n4BjJcl9C9OeazPce+JqT2vgvzFNPwf4T4t6L4rtJqf7sU6BJDl8Be4AmgWEi87sq8BhJ+9Nx+CF3e5lZjdAVZlYOP974cfzkTaVG4NxGAl3wfT6uc84dL+FjZOI7yx8EEnm4yIiZWXX8cJsQfijmUEl5/wGYmZHbnKdE51pIpnsvUIMwKfDn1QW8JZrPQ5xza4Bl+EklLz3Z/SUKM3sC32n6e+BC59y2U3CMcuRe00g+HxJW4P/vIYE/IzrXZL33QgSbm71d0j8uJsq9F813k1L/2Xc6Zu/UEvsFeAI/Y+tXQMWQ14cGXp+S7/134b/YnDALOb5q1+EnTSoT8npwltlhsT7fEr52qfhfPRx+tK4KEWxT4PXDT/jYt4D3tyV3xt0XYn3OJXz9zgMGA6n5Xm+Ib3PvgJGRXL/AuqS6/0LO7/zA+W0AUgp5X9LfexQ9W3u/wPn+ADQN0AdZVQAABzxJREFUeb17YNtdwFn5tukSuK4TC9jfbYH9rQBqhLx+ReD170Lv1XhfIrh+9wbOa3Po9StinwVeP6AFvuaxXL7XqwOfBo4zH7BYX5eTvXaBc7oTqJzv9UrAyyHXtEK+9br3TnxfBXzfDAc0T8Z7j+J9Nym1n30a5Sx5PIG/kc8DvjOzqfi2ll3xcwbcmu/91fC/2BbUdvQeoBt+iMnlgZGqMvFDwH6HD0mlyV3k/nL2A34ElILe91vn3A+B5+GuXxfgD2a2Dl9dexA/TGQHfOfDKcADJVn4ONAM36Z3i5nNxY9GkwF0xE8qtwT4eb5tdP+dKPhr5LvO93kJJ+nuPTO7lLzDk5YNvB46UtTjzrnPAJxz/zWzF/BzNcw3swmBbS7Ejw50izux+VRwuPHyBRThNXxgHIK/Jyfi/zv0xs+i/RMXx3NZRHP9ArX8zwReWwP8vzCfh/92zk0L+Tvc9asFvAm8EPh/eTtQB//5UBkf4K9xgW9J8SbKe68i8A/gKTObhQ8v1fH/D1bFfzZe5U5sTqp7L/D/bojB+CA4yzn3bRGHKZX3HsX4blKqP/tinTC1nL4FXy34GH4c/Cz8h+nrQL0C3vsIhfzaDVQB/oZvbpAVeHyBfMm+NCwh16KopWFR1w//K8ir+DalP+DHid8BTMb/0pF6Os/tNF2/lsCLwBz8SFpH8f9wT8eHj3Tdf0Vew3L42Z4d0DbC+zVp7j3g5gj+/7w5zHaz8W3kd+HboZ8X5hh9AvtZG2Z9auB+Xoz/h/wH/OhfrWJ9fUry+oVch6iud7jrh/8C+Vzg82AzcAT/y/ucwL18dqyvTwleu8rAU/gfDzbgfxE/ELhn/grU1b0X8f+7YwLr7o7gGKX13nskwv8XG4a57qXqs88ChREREREREUk4GhRAREREREQSlgKNiIiIiIgkLAUaERERERFJWAo0IiIiIiKSsBRoREREREQkYSnQiIiIiIhIwlKgERERERGRhKVAIyIiIiIiCUuBRkREREREEpYCjYiIiIiIJCwFGhERERERSVgKNCIiIiIikrAUaEREJKGYmcu3ZJvZHjP7xszuMbO0IrY/x8yOmNlfTnE5BwfKd82pPI6ISLIz51ysyyAiIhIxMwv+w/VG4DEVaAich/+hbiJwsXPuWJjthwOXAI2cc9tPYTkNmAdUAlo6546eqmOJiCQzBRoREUkowUDjnLN8r3cFpgDlgRudc28XsG0HYA7wN+fcb05DWa8F3gPucM69dKqPJyKSjNTkTERESgXn3AxgWODPi8K87VeBxzdPeYG8kcA+4Jen6XgiIklHgUZEREqTJYHHGvlXmFkl4FrgO+fcnALWDwv0eeljZv3M7Esz22dm28zsFTM7M/C+Gmb2f2a20cwOm9lMM+tTUGGcc4eAEUDbQA2SiIiUMAUaEREpTSoHHrcVsK43vj/LlCL2MQQYB1jgMQu4DRhpZtWA6fgaoKnAfKAzMM7M2oTZX/B4l0Z0BiIiEhUFGhERKU0uDjyOK2Bdr8DjrCL2cRcwyDnXyzl3NdAKWIwPRF8AXwNNnXPXOue6AU8A5YDfhtnfzMBj78hOQUREoqFAIyIiCc3MUsysiZm9BJyP77fyQQFvbRt4/LaIXb7rnPss+Idzbh/wSuDPesDd+UYs+yvgCB9Ylgcezy3iuCIiUgxlYl0AERGR4ggZvjnUK8DtruAhPIP9anYVsevxBby2OvA42zmXZ3vn3B4z2wnULmhnzrljZrYPOMPMyjrnjhRxfBERiYICjYiIJKrgPDTlgXZAC+Dn+CZhwwp4/5mBx31F7HdjAa/tL2RdcH3VQva5F9+/5ywK7t8jIiLFpEAjIiIJyTl3c+jfZvY74Gngn2Y22Tm3Lt8mewKPlSlcdjHXFSYYpnYXc3sREQlDfWhERKRUcM79Bd9crALwhwLeEqwZqXLaCgWYWRp+dLW9am4mIlLyFGhERKQ0eSDweKOZZeRbtyDw2Pw0lgd8UzjwQzyLiEgJU6AREZFSwzk3Dz+RZRngf/Ktnhp47HxaCwVdAo9fnObjiogkBQUaEREpbR7BD6N8q5nVCnn9S3zn/T6nuTzB431W2JtERKR4FGhERKRUcc4tAD7Fj342NOT1/cB7wDlmdlpqacwsHRgMLHTOzTgdxxQRSTZW8FD9IiIipY+ZnQvMA/7hnPv1aTjedcC7wB3OuZdO9fFERJKRAo2IiCQVM/sQuBRo5Jw7ZXPCmJnhw1MloJVGOBMROTXU5ExERJLN/wJpnDhoQEkbhJ/w80GFGRGRU0c1NCIiIiIikrBUQyMiIiIiIglLgUZERERERBKWAo2IiIiIiCQsBRoREREREUlYCjQiIiIiIpKwFGhERERERCRhKdCIiIiIiEjCUqAREREREZGEpUAjIiIiIiIJS4FGREREREQSlgKNiIiIiIgkLAUaERERERFJWAo0IiIiIiKSsBRoREREREQkYf1/qvmM2BIFHD8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 900x600 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot lateral spread for different multiplicities #########################\n",
    "\n",
    "def plot_lateral_spread_M():\n",
    "    Rmin = 0.0\n",
    "    Rmax = 20\n",
    "    npointsR = 200\n",
    "    deltaR = (Rmax-Rmin)/npointsR\n",
    "    Rs = []\n",
    "    Spreads = []\n",
    "    theta = 0.0\n",
    "    ms = [2.0 , 3.0 , 4.0]\n",
    "    for i in range(len(ms)):\n",
    "        spread = []\n",
    "        rs = []\n",
    "        for j in range(npointsR):\n",
    "            r = Rmin + deltaR * j\n",
    "            rs.append(r)\n",
    "            spread.append(dN_dR(ms[i] , theta , h , r, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ))\n",
    "        Spreads.append(spread)\n",
    "        Rs.append(rs)   \n",
    "    fig = plt.figure()\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.plot(Rs[0], Spreads[0],'black', label=r'M = 2')\n",
    "    ax.plot(Rs[1], Spreads[1],'red', label=r'M = 3')\n",
    "    ax.plot(Rs[2], Spreads[2],'blue', label=r'M = 4')\n",
    "    plt.xlabel(r'R(m)', fontsize = 10)\n",
    "    plt.ylabel(r'dN/dR', fontsize = 10)\n",
    "    plt.rcParams['figure.dpi'] = 150\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    \n",
    "plot_lateral_spread_M()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interactive look at the lateral distribution when one varies the depth $h$, zenith angle and multiplicity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lateral_distance = np.arange(0,20, 0.01)\n",
    "zenith_angle = 0 #radians\n",
    "\n",
    "opts = dict(continuous_update=False, readout=True,readout_format='.3f')\n",
    "\n",
    "def plot_basic_parameters( multiplicity=widgets.FloatSlider(min=2, max= 4, value=multiplicity,**opts),\n",
    "       h=widgets.FloatSlider(min=2.0,max=4.0,value=h,**opts),\n",
    "       zenith_angle=widgets.FloatSlider(min=0,max=np.pi/2,value=zenith_angle,**opts)):\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(8,4))\n",
    "    ax.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b,rho_0_a, rho_0_b, rho_1, f, theta_0 ))\n",
    "    ax.set_xlabel(r'R(m)', fontsize = 10)\n",
    "    ax.set_ylabel(r'dN/dR', fontsize = 10)\n",
    "    ax.set_ylim(0.0,0.12)\n",
    "\n",
    "interactive_plot = interactive(plot_basic_parameters)\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout = {'height': '600px'}\n",
    "interactive_plot\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\alpha$ parameters\n",
    "\n",
    "Change these parameters interactively & statically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#array = np.arange(0,20, 0.01)\n",
    "#mult=2\n",
    "#angle = 0\n",
    "\n",
    "#def plot( alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b):\n",
    "#    plt.figure(3)\n",
    "#    plt.plot(array, dN_dR(mult, angle, h, array, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b))\n",
    "#    plt.xlabel(r'R(m)', fontsize = 10)\n",
    "#    plt.ylabel(r'dN/dR', fontsize = 10)\n",
    " #   plt.ylim(0,0.15)\n",
    "#    plt.show()\n",
    "\n",
    "#interactive_plot = interactive(plot, alpha_0_a=(alpha_0_a-0.01,alpha_0_a+0.01), alpha_0_b=alpha_0_b, alpha_1_a=alpha_1_a, alpha_1_b=alpha_1_b)\n",
    "\n",
    "#interactive_plot\n",
    "\n",
    "lateral_distance = np.arange(0,20, 0.01)\n",
    "zenith_angle = 0 #radians\n",
    "\n",
    "opts = dict(continuous_update=False, readout=True,readout_format='.4f')\n",
    "\n",
    "def plot_alpha_parameters( alpha_0a=widgets.FloatSlider(min=alpha_0_a - 0.05, max= alpha_0_a + 0.05, value=alpha_0_a,step=0.001,**opts),\n",
    "                          alpha_0b=widgets.FloatSlider(min=alpha_0_b - 0.05,max=alpha_0_b + 0.05, value=alpha_0_b,step=0.01,**opts),\n",
    "                          alpha_1a=widgets.FloatLogSlider(min=-1.8,max=-1.5, value=alpha_1_a,step=0.0001,**opts),\n",
    "                          alpha_1b=widgets.FloatSlider(min=alpha_1_b - 0.005,max=alpha_1_b + 0.005, value=alpha_1_b,step=0.001,**opts)):\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(8,4))\n",
    "    ax.plot(lateral_distance, dN_dR(2, 0, h, lateral_distance, alpha_0a, alpha_0b, alpha_1a, alpha_1b, rho_0_a, rho_0_b, rho_1, f, theta_0 ))\n",
    "    ax.set_xlabel(r'R(m)', fontsize = 10)\n",
    "    ax.set_ylabel(r'dN/dR', fontsize = 10)\n",
    "    ax.set_ylim(0,0.12)\n",
    "\n",
    "interactive_plot = interactive(plot_alpha_parameters)\n",
    "output = interactive_plot.children[-1]\n",
    "output.layout = {'height': '600px'}\n",
    "interactive_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ), label = 'alpha_0_a: '+ str(alpha_0_a) + ' (nominal)')\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a + 0.5, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ), label = 'alpha_0_a: '+ str(alpha_0_a+ 0.5)  )\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a - 0.5, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ), label = 'alpha_0_a: '+ str(alpha_0_a- 0.5) )\n",
    "plt.xlabel(r'R(m)', fontsize = 10)\n",
    "plt.ylabel(r'dN/dR', fontsize = 10)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ), label = 'alpha_0_b: '+ str(alpha_0_b) + ' (nominal)')\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b + 0.5, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ), label = 'alpha_0_b: '+ str(alpha_0_b+ 0.5)  )\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b - 0.5, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ), label = 'alpha_0_b: '+ str(alpha_0_b- 0.5) )\n",
    "plt.xlabel(r'R(m)', fontsize = 10)\n",
    "plt.ylabel(r'dN/dR', fontsize = 10)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ), label = 'alpha_1_a: '+ str(alpha_1_a) + ' (nominal)')\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a + 0.05, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ), label = 'alpha_1_a: '+ str(alpha_1_a + 0.05)  )\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a - 0.05, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 ), label = 'alpha_1_a: '+ str(alpha_1_a - 0.05) )\n",
    "plt.xlabel(r'R(m)', fontsize = 10)\n",
    "plt.ylabel(r'dN/dR', fontsize = 10)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0), label = 'rho_0_a: '+ str(rho_0_a) + ' (nominal)')\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a + 1.0, rho_0_b, rho_1, f, theta_0), label = 'rho_0_a: '+ str(rho_0_a + 1.0)  )\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a - 1.0, rho_0_b, rho_1, f, theta_0), label = 'rho_0_a: '+ str(rho_0_a - 1.0) )\n",
    "plt.xlabel(r'R(m)', fontsize = 10)\n",
    "plt.ylabel(r'dN/dR', fontsize = 10)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0), label = 'rho_0_b: '+ str(rho_0_b) + ' (nominal)')\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b + 5, rho_1, f, theta_0), label = 'rho_0_b: '+ str(rho_0_b + 5.0)  )\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b - 5, rho_1, f, theta_0), label = 'rho_0_b: '+ str(rho_0_b - 5.0) )\n",
    "plt.xlabel(r'R(m)', fontsize = 10)\n",
    "plt.ylabel(r'dN/dR', fontsize = 10)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0), label = 'rho_1: '+ str(rho_1) + ' (nominal)')\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1 + 0.2, f, theta_0), label = 'rho_1: '+ str(rho_1 + 0.2)  )\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1 - 0.2, f, theta_0), label = 'rho_1: '+ str(rho_1 - 0.2) )\n",
    "plt.xlabel(r'R(m)', fontsize = 10)\n",
    "plt.ylabel(r'dN/dR', fontsize = 10)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0), label = 'f: '+ str(f) + ' (nominal)')\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f + 5.0, theta_0), label = 'f: '+ str(f + 5.0)  )\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f - 5.0, theta_0), label = 'f: '+ str(f - 5.0) )\n",
    "plt.xlabel(r'R(m)', fontsize = 10)\n",
    "plt.ylabel(r'dN/dR', fontsize = 10)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0), label = 'theta_0: '+ str(theta_0) + ' (nominal)')\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 + 0.7), label = 'theta_0: '+ str(theta_0 + 0.7)  )\n",
    "plt.plot(lateral_distance, dN_dR(multiplicity, zenith_angle, h, lateral_distance, alpha_0_a, alpha_0_b, alpha_1_a, alpha_1_b, rho_0_a, rho_0_b, rho_1, f, theta_0 - 0.7), label = 'theta_0: '+ str(theta_0 - 0.7) )\n",
    "plt.xlabel(r'R(m)', fontsize = 10)\n",
    "plt.ylabel(r'dN/dR', fontsize = 10)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}