\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.18.2?urlpath=lab/tree/astro_dark_matter.ipynb)\n",
"- You can contribute with your own notebooks in this\n",
"[GitHub repository](https://github.com/gammapy/gammapy/tree/master/docs/tutorials).\n",
"- **Source files:**\n",
"[astro_dark_matter.ipynb](../_static/notebooks/astro_dark_matter.ipynb) |\n",
"[astro_dark_matter.py](../_static/notebooks/astro_dark_matter.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Dark matter utilities\n",
"\n",
"## Introduction \n",
"\n",
"Gammapy has some convenience methods for dark matter analyses in `~gammapy.astro.darkmatter`. These include J-Factor computation and calculation the expected gamma flux for a number of annihilation channels. They are presented in this notebook. \n",
"\n",
"The basic concepts of indirect dark matter searches, however, are not explained. So this is aimed at people who already know what the want to do. A good introduction to indirect dark matter searches is given for example in https://arxiv.org/pdf/1012.4515.pdf (Chapter 1 and 5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"As always, we start with some setup for the notebook, and with imports."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from gammapy.astro.darkmatter import (\n",
" profiles,\n",
" JFactory,\n",
" PrimaryFlux,\n",
" DarkMatterAnnihilationSpectralModel,\n",
")\n",
"\n",
"from gammapy.maps import WcsGeom, WcsNDMap\n",
"from astropy.coordinates import SkyCoord\n",
"from matplotlib.colors import LogNorm\n",
"from regions import CircleSkyRegion\n",
"import astropy.units as u\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Profiles\n",
"\n",
"The following dark matter profiles are currently implemented. Each model can be scaled to a given density at a certain distance. These parameters are controlled by ``profiles.DMProfile.LOCAL_DENSITY`` and ``profiles.DMProfile.DISTANCE_GC``"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[gammapy.astro.darkmatter.profiles.NFWProfile,\n",
" gammapy.astro.darkmatter.profiles.EinastoProfile,\n",
" gammapy.astro.darkmatter.profiles.IsothermalProfile,\n",
" gammapy.astro.darkmatter.profiles.BurkertProfile,\n",
" gammapy.astro.darkmatter.profiles.MooreProfile]"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"profiles.DMProfile.__subclasses__()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"LOCAL_DENSITY: 0.3 GeV / cm3\n",
"DISTANCE_GC: 8.33 kpc\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd1zWVfvH34ebJUNUQFFAlorIXrISwQFqiqPUcqXmrp6eLDObVjZ+ak89Vj5mpjbMKC1zNJyICipiuAeyZKgMBRVB4Ob7++MW3AoCMjxvX/fr5rvOub63cH3Pfc51fS6hKAoSiUQiafpo1bcBEolEInk4SIcvkUgkjwjS4UskEskjgnT4EolE8oggHb5EIpE8IkiHL5FIJI8I2vVtwL0wMzNTbG1t69sMiUTykDhx4gQAjo6O9WxJ4yY+Pj5XURTzW/c3aIdva2vLvn376tsMiUTykAgJCQEgKiqqXu1o7Agh0u60v0E6fCHEAGBAhw4d6tsUiUTyEHnzzTfr24QmjWjImbY+Pj6KHOFLJBJJ9RBCxCuK4nPrfrloK5FIGgwJCQkkJCTUtxlNlgY5pSORSGqH0tJSMjIyKC4urm9TqsSFCxcAOHbsWD1b0jjQ19fHysoKHR2dKp0vHb5E0oTJyMjA2NgYW1tbhBD1bc590dLSTDrIKJ37oygKeXl5ZGRkYGdnV6VrGuSUjhBigBBicUFBQX2bIpE0aoqLizE1NW0Uzl5SPYQQmJqaVuvbW4N0+IqirFMUZZKJickDXX8mqYC8zMu1bJVE0jiRzr7pUt3/2wbp8GvKrlWJ/DRnL5uXHeViblF9myORPNIIIXj55Zcrt+fPn8/s2bMBmD17NpaWlnh4eODh4cEnn3zCli1bGDRoUOX5H330ETeGaK9bt46IiAhAk6vj6uqKu7s7YWFhnD17tlq2zZgxA2dnZ2bMmMGiRYv47rvvABg7diyrVq160FtusDTJOfz+z7mz/+80DkZlkLjvHM7dLPHpZ4tBc936Nk0ieeTQ09Pj119/ZdasWZiZmd12/KWXXuKVV14B4PLly+Tk5PDuu+9WHo+NjaV58+ZkZ2fTunVrYmJiCAoKqjy+bds2zMzMeP311/nwww9ZsGBB5TFFUVAUpXJt4Fa++uorcnJy0NPTq63bbdA0yRG+vpEOgU90YNR7AXQObMvh6Ey+fyuW3b8ncbWorL7Nk0geKbS1tZk0aRKffvrpfc81MjLCzs4OExMTTp06BUBmZiZPPPEEMTExAMTExBAYGHjbtcHBwZw6dYrU1FScnJyYNm0aXl5epKenM2PGDFxcXHB1dSUyMhKAiIgICgsL8fPzIzIyktmzZzN//vzb2o2Pj6d79+54e3sTHh7OmTNnavJx1CsPdYQvhDAEooF3FEVZX9f9GbXUI3RkZzx7tWfPumTi/0zjcHQmXuE2uIVYoa2rqmsTJJIGw7vrjnA062KtttmlXXPeGeB83/Oee+453NzcePXVV2879umnn/LDDz8AmimeXr16ERgYSExMDGq1mo4dO+Lv78/ff/9N//79OXjwIL6+vre1s379elxdXQGNJs+yZctYuHAhq1evJiEhgQMHDpCbm4uvry/BwcGsXbsWIyOjyrj/immmGyktLeWFF17g999/x9zcnMjISN544w2WLl1anY+pwVCjEb4QYqkQIlsIcfiW/X2EECeEEKeEEK/dcGgm8HNN+nwQWrQxIHyCC8Ne96WNrQmxvybxw9u7ObIjk3J1+cM2RyJ55GjevDljxoy5abqlgpdeeqky4crJyYnMzEyCgoKIiYkhJiaGgIAAunbtyp49e/jnn39wdHREX1+/8vrQ0FA8PDy4ePEis2bNAsDGxgZ/f38Adu7cydNPP41KpaJNmzZ0796duLi4Ktl94sQJDh8+TO/evfHw8GDOnDlkZGTUwidSP9R0hL8c+AL4rmKHEEIFfAn0BjKAOCHEWqAdcBTQv72Zh4N5e2MGvOBO5skL7F6TRNSKE/yz6TR+EfZ08GqN0JLRDJKmS1VG4nXJv//9b7y8vBg3btx9zw0MDOTzzz9HrVYzceJEjI2NKS4uJioq6qb5e7g+h19Bfn4+hoaGlds1kY9RFAVnZ2diY2MfuI2GRI1G+IqiRAPnb9ndFTilKEqyoiglwE/AQCAU8AdGABOFEPW2fmDZqSVDZnjTb6orKm0tNi45ws8fxZF2JK9GvxwSieTutGrVimHDhvHNN9/c99wuXbqQlZXFjh078PT0BMDDw4NFixbdcf7+XgQHBxMZGYlarSYnJ4fo6Gi6du1apWsdHR3JycmpdPilpaUcOXKkWv03JOrC6VoC6TdsZwCWiqK8oSjKv4Efga8VRbnjXIoQYpIQYp8QYl9OTk4dmFfZD3bu5gx/syu9xjpx9UoZ6z8/wJr//MPZZJnwJZHUBS+//DK5ubn3PU8IgZ+fH2ZmZpWyAQEBASQnJ1fb4Q8ePBg3Nzfc3d3p0aMHc+fOxcLCokrX6urqsmrVKmbOnIm7uzseHh6Vi8eNkRqrZQohbIH1iqK4XNseCoQrijLh2vZooKuiKC9Uo80KeeSJiYmJNbKvqqjLyjmyI4t9f6ZSdLEEWzcz/AfaY2pp9FD6l0jqgmPHjuHk5FTfZlQZWQCl+tzp//huapl1EaWTAVjfsG0FZNVBP7WKSlsLt1ArnALbcmBLOv9sTOOnOXvp1LUNfgPsaW7WrL5NlEiaPNbW1vc/SfLA1MWUThzQUQhhJ4TQBZ4C1langZpKK9QEHT0VPv1sGf1BIJ6925O0P4cV7+wm+qeTXLlY8tDtkUgeJQwMDDAwMKhvM5osNRrhCyFWAiGAmRAiA018/TdCiOeBvwEVsFRRlGqtcjSEilf6hjoEDumAW6g1cX+kcDg6k2MxWbj3tMYzzAa9Zk0ySVkiqVcuXtTkCTRv3ryeLWmayIpXVST/3BX2rEvm1L5s9Ay1ZfKWpFEg5/CbPtWZw2+Q0goNUR755uSt5jclb6ll8pZEImkENEiHX59z+PdDk7zlwaDpnhi30iNqxQlWvruHxH3nUMob7rcliUQiaZAOv6Yj/B/3nGbpzhQKr9adUJpM3pJIqoZKpaqUP/bw8ODjjz8GYMKECRw9erTW+klNTeXHH3+873lRUVGYmJjg6emJk5PTTcqcVSEnJwc/Pz88PT3ZsWMH/fr1Iz8/H9CIvzVkGuTKo6Io64B1Pj4+Ex/k+p2ncvjj0Fn+uyWRUf7teSbQltbGta/oUJG8ZeNqRmLcOfauS2b95wdo17EFAYMdsLBveN9QJJKHTbNmze5YmHzJkiW12k+Fwx8xYsR9z+3WrRvr16+nsLAQDw8P+vfvj7e3d+XxsrIytLXv7B63bNlC586d+fbbbyvbaiw0yBF+TVk40ptfpwUS6GDKwqgkHvt4GzNXHeRU9qU66U9LS+DoZ8GI2f4EP9WJC+eusHpuPBsWHpSVtySSuxASEkJFUIaRkRFvvPEGQ4cOZfTo0Zw7dw7QFDupGE336tWrcv/27dsrvzF4enpy6dIlXnvtNXbs2IGHhweffvopxcXFjBs3DldXVzw9Pdm2bdttNhgaGuLt7U1SUhKzZ89m0qRJhIWFMWbMGNLS0ujZsydubm707NmT06dPk5CQwKuvvsoff/yBh4cHRUVF2Nra3jF7eN68efj6+uLm5sY777xTh59k1WmQI/zaCMv0at+S/43yJjW3kCU7k/llXwaR+9Lp5dSaid3s6WrXqtZLv6m0tXANsaJzQFsObE3nn42n+WnOXhy7WtB1gJ1M3pLUL3++BmcP1W6bFq7Q9+N7nlJUVISHh0fl9qxZsxg+fPhN5xQWFuLv788HH3zAq6++ytdff82bb77JY489xu7duxFCsGTJEubOncsnn3zC/Pnz+fLLLwkKCuLy5cvo6+vz8ccfM3/+fNav1yivf/LJJwAcOnSI48ePExYWxsmTJ2/qNy8vj927d/PWW29x9OhR4uPj2blzJ82aNWPAgAGMGTOGZ555hqVLl/Kvf/2LNWvW8N5777Fv3z6++OKLu97zxo0bSUxMZO/evSiKQkREBNHR0QQHB1fr461tGqTDr+mUzo3YmhkyZ5ArL/XqxHexaXwXm8rmY7txt27B5GB7wp0tUNWySqaOngqfvra4BFtqKm9tk5W3JI8ud5vSuRFdXV369+9Pfn4+nTt3rtSrycjIYPjw4Zw5c4aSkhLs7OwACAoKYvr06YwcOZIhQ4ZgZWV1W5s7d+7khRc0ii6dO3fGxsam0uFXiLJpaWnx2muv4ezszC+//EJERATNmmkGZrGxsfz6668AjB49+o5a/ndj48aNbNy4sVL47fLlyyQmJkqH/7AwNdLjpd6dmNLdgVX7M1iyI5lpK/ZjY2rAhMfseNLbmma1HFMvk7ckDYr7jMTrEx0dHYQQnDt3jkuXLlFWpgm4eOGFF5g+fToRERFERUVVFil57bXXePzxx/njjz/w9/dn8+bNt7V5r+CJijn8W7lRVvlWqjMjoCgKs2bNYvLkyVW+5mHQIOfw6zIOv5muitH+Nmx9OYT/jfSipYEub/1+hMCPt/CfTSfJu3y11vusqLw14h0/7NzMiP8zje/fjGH/xjTKStS13p9E0lQoKCjA0tISoHKRFCApKQlXV1dmzpyJj48Px48fx9jYmEuXrq/TBQcHs2LFCgBOnjzJ6dOnq5XQFRgYyE8//QTAihUreOyxx6p8bXh4OEuXLuXyZc0aXmZmJtnZ2VW+vq5okEPM2pzSuRsqLUFf17b0cbFgX9oFvtqezIItiXy1PYknva2Y0M0eO7O7P+0fhBZtDAib4IJn2CV2/55M7K9JHNyage/jtjgFtkVL1SCfvxJJjbh1Dr9Pnz6VoZn3Y/bs2QwdOhRLS0v8/f1JSUkB4LPPPmPbtm2oVCq6dOlC37590dLSQltbG3d3d8aOHcu0adOYMmUKrq6uaGtrs3z58moVK1+wYAHjx49n3rx5mJubs2zZsipfGxYWxrFjxwgICAA0i9I//PADrVu3rnIbdYGUVriBU9mXWbIjmV/3Z1JaXk5YlzZMCnbA26ZlnfSXlXiB2N+SOZtcgEnrZrLylqTWkdIKTZ9GL61QX3RobcTHT7ixc2Yo00IciE3K44n/xfDk/2LYeOQs5bWcSduuY0uGzPCi3zQ3mbwlkUjqHDnCvweFV8uIjEvnm50pZOYXYW9uyMRu9gz2tERfp3YXeMvLFRL3nmXPuhQu5RXL5C1JrdDYRvglJRoJcl1dGclWVaozwm+QDr8+Kl7dizJ1OX8cPsvi6CQOZ17EzEiPsYE2jPK3oYVB7f5iqsvKObozi7g/ZOUtSc1pbA5fUn0avcOvoL5H+LeiKAqxSXks3pFM1IkcmumoGO5rzbOP2WHdqnaLNpReVXNwWzr7/z5NSXGZTN6SPBCNzeGfP38e0BQ8l1SN+i5x2GQRQhDYwYzADmacOHuJxdHJrNijSebq59qWScH2uFm1qJW+dPRUePexxbnbLclbj7XDu58thiZVjzaQSBoLOTk5gHT4dYV0+A+Io4UxnwxzZ0a4I8t2pfDjntOsP3gGf/tWTA52IMTRvFakG25M3tr3RwqHd2RxLPYM7j2s8Qxrj56BTi3cjUQieRSQUTo1xMJEn1n9nIiZ1YM3+jmRlneFccvjCP8sml/2pVNSVjvFUYxa6hFyLXnL1s2M+L/S+P7NWJm8JWnwPKhk8Icfflj5c2pqKi4uLrVl0gNxo9ibra0trq6uuLu7ExYWxtmzZ6vV1owZM3B2dmbGjBksWrSI7777DoCxY8eyatWqWre9gofm8IUQTkKIRUKIVUKIqQ+r34eFsb4OE4PtiX41lP8Mc0dLCGasOki3uVv5X1QSBUWltdLPTZW37Ew0lbfeipWVtyRNjhsdfk2pkGqoTbZt28aBAwfw8fG5zVZFUSgvv/vf41dffcX+/fuZN28eU6ZMYcyYMbVu352okcMXQiwVQmQLIQ7fsr+PEOKEEOKUEOI1AEVRjimKMgUYBty2mNBU0FFpMcTLij9f7Ma347vSobUR//fXcYI+3sqc9UfJyi+qlX40lbfcGfyyJ8amzWTlLUmD58yZMwQHB+Ph4YGLiws7duwAYOXKlbi6uuLi4sL8+fMBjVZORYbuyJEjAVCr1UycOBFnZ2fCwsIoKtL8LSUlJdGnTx+8vb3p1q0bx48fBzSj5enTpxMaGsrMmTMZO3YsU6dOJTQ0FHt7e7Zv38748eNxcnJi7NixlXZOnToVHx8fnJ2dqyRrHBwczKlTp0hNTcXJyYlp06bh5eVFeno6M2bMwMXFBVdXVyIjIwGIiIigsLAQPz8/IiMjmT17duV930h8fDzdu3fH29ub8PBwzpw58+Af/jVqOoe/HPgC+K5ihxBCBXwJ9AYygDghxFpFUY4KISKA165d06QRQtC9kzndO5lzOLOAxdHJLItJZXlMKgPc2zEp2B6nts1r3E9F8lbqoTx2r0li45Ij7LdOw3+QA+271L4EtKTx8n97/4/j54/XapudW3VmZteZVTr3xx9/JDw8nDfeeAO1Ws2VK1fIyspi5syZxMfH07JlS3r37s2hQ4f4+OOP+eKLLypVNlNTU0lMTGTlypV8/fXXDBs2jNWrVzNq1CgmTZrEokWL6NixI3v27GHatGls3boV0GjobN68GZVKxdixY7lw4QJbt25l7dq1DBgwgF27drFkyRJ8fX1JSEjAw8ODDz74gFatWqFWq+nZsycHDx7Ezc3trve1fv16XF1dAU2m8LJly1i4cCGrV68mISGBAwcOkJubi6+vL8HBwaxduxYjI6PKe6sQhLuR0tJSXnjhBX7//XfMzc2JjIzkjTfeYOnSpdX577mNGjl8RVGihRC2t+zuCpxSFCUZQAjxEzAQOKooylpgrRBiA3DHWmRCiEnAJID27dvXxLwGg4ulCQue9uTVPo4s3ZnKT3Gn+e2fTLp1NGNysANBHUxr5JiFENi5mWHjYkpi3Dn2rJWVtyQND19fX8aPH09paSmDBg3Cw8ODrVu3EhISgrm5OQCjRo0iJiaGJ5988rbr7ezsKjV5vL29SU1N5fLly8TExDB06NDK865evS6AOHToUFSq60mSAwYMQAiBq6srbdq0qXTUzs7OpKam4uHhwc8//8zixYspKyvjzJkzHD169I4OPzQ0FJVKhZubG3PmzCE/Px8bGxv8/f0BjTzz008/jUqlok2bNnTv3p24uDgiIiLu+1mdOHGCw4cP07t3b0Dz7aZt27b3ve5+1EWUjiWQfsN2BuAnhAgBhgB6wB93u1hRlMXAYtDE4deBffWGVUsD3h7QhRd7dmTF3jSW7Upl1Dd76NK2OZOC7XncrS06NRBQq6i81cG7NUd2ZLHvz1RWz42XyVsSgCqPxOuK4OBgoqOj2bBhA6NHj2bGjBk0b37zt9xLly5VTtXcyo3CZyqViqKiIsrLy2nRosVd9fZvlTuuaENLS+um9rS0tCgrKyMlJYX58+cTFxdHy5YtGTt2LMXFxXdse9u2bZiZmVVu5+fn39RfTXKcFEXB2dmZ2NjYB27jTtTFou2dhqqKoihRiqL8S1GUyYqifHnPBupQHrkhYGKgw7SQDuycGcrcJ9woUZfz78gEQuZFsWRHMpdrWHxdpa2FW6gVo97zxy/CnqyTF/hpzl42LzvKxdzaWUOQSKpLWloarVu3ZuLEiTz77LPs378fPz8/tm/fTm5uLmq1mp9//rly1K2jo0Np6b2DHZo3b46dnR2//PILoHGUBw4ceGAbL168iKGhISYmJpw7d44///zzgdsKDg4mMjIStVpNTk4O0dHRdO3atUrXOjo6kpOTU+nwS0tLOXLkyAPbUkFdjPAzAOsbtq2ArDrop9Gjp61imK81T3pbse1ENl9FJzNnwzEWbElkpL8N4wJtad38wYuv6+pr49PvhspbUbLylqT+iIqKYt68eejo6GBkZMR3331H27Zt+eijjwgNDUVRFPz9/enZsycAkyZNws3NDS8vLz744IO7trtixQqmTp3KnDlzKC0t5amnnsLd3f2BbHR3d8fT0xNnZ2fs7e0JCgp6oHYABg8eTGxsLO7u7gghmDt3LhYWFlW6VldXl1WrVvGvf/2LgoICysrK+Pe//42zs/MD2wO1IK1wbQ5/vaIoLte2tYGTQE8gE4gDRiiKUu3HU0OTVngY/HP6Al/vSOavw2fR1tJikKdmgbdDa+Mat335wlXiNqRwLOYMKh0t3HtYycpbTZzGJq0g5ZGrz0OTVhBCrARCADMhRAbwjqIo3wghngf+BlTA0uo6+9ooYt5Y8WzfkoUjvUnLK2TJjhR+iU/n530Z9OzcmonB9vjVoPi6UUs9Qkd1xrN3e/asTSb+zzQOR2fiHW6La4gl2rVc4lEikTQsGqR4Wo3VMsvVoNU0nNf5whK+v1Z8Pa+wBHcrEyYFO9DHpebF13NOX2L3miROHz2PYQs9WXmrCSJH+E0fqZYZOQoy94O5I5g5at7NO2veDRqnKFNxqZrV+zNYsiOFlNxC2rcyYEI3O4bWQvH1zBMXiF2TxLmUi7RoY4BfhD0Onuay8lYToLE5fLVaIxNyYyil5N5Ihx//LaTtgpwTkHsSSq9cP2bYWuP4WztpXuZO0LozNKubMoa1jbpcYdPRcyyOTmL/6XxaGugwOsCWZwJsMDV6cAVNRVFIOZDLnrXJnM8qxLy9Mf4D7bGWyVuNmsbm8CXVp9E7/FotgFJeDhczNM4/5zhkH4ecY5rtksvXzzNuB226QOsu0MZZ827uCNoNV4Z4X+p5Fm1PZvOxc+hpa/GEtxUTa1h8vbxc4eTes+y9VnnLslML/AfJ5K3GSmNz+NnZ2QD1Xuy7MdHoHX4FdRqloyiQf/raQ+AonDsK2ccg9wSoNWXWECow66h5ALRx0bwsXMHYAhrQqDcpR1N8ffX+TErVFcXX7fG2efDpK3VpOUd2ZrLvj1SKLpVi526G30B7TNvJ5K3GRGNz+HIOv/o0eodfryUO1WVwPgnOHbn2Oqx5GBScvn6OganG8Vu4goWb5t20I6jqN7wx59JVvo1J5fvdaRQUleJt05JJwfb0dmqD1gPOx5cUl3Fwazr7N56m9Kqazn4W+PaXlbcaCw3B4atUKlxdXVEUBZVKxRdffEFgYOAdz72bww8JCWH+/Pn4+FRfdzE/P58ff/yRadOmAVSKnDk6OlJSUkJwcDALFy5ES6tqwQpXr17l8ccfJzc3l1mzZrFp0yamT59Oly5dsLW1Zd++fTdl4NY1jb7ilaIo64B1Pj4+Ex965yrta4u8juAy5Pr+onzNN4Gzh+HsQTh7CPYsBvU13Q5tfc03AQs3aOsO7Tw000IPcUrI3FiPV8IdmRriwM/70lmyI4XJ38djb27IhMfsGeJV/eLrmuQtO5yDLdn/VxqHojI5GXcOl2BLvPvK5C3J/WnWrFml9MHff//NrFmz2L59e5Wvr1jIfRDUajX5+fksXLiw0uEDODg4kJCQQFlZGT169GDNmjUMGXL9772srAxt7Tu7x3/++YfS0tLKexo+fPgD2/ewaZAOv0HSrAXYBGpeFajLIC8RzhzUPATOHIDDv0L8Ms1xLR3NwnA7D2jrAe08NdNC2nXrJA31tBkXZMdof5vK4uuv/3aI/2w6wTMBtowOqH7x9WZGugQ92RG3Htbs25DCoagMjsacwaOXNZ692qMrk7ckVeDixYu0bKkJkIiKimL+/PmsX78egOeffx5LS0uGDBmCra0t48ePZ+PGjTz//POV15eXlzNu3Disra159913ee2114iKiuLq1as899xzTJ48maioKN59913atm1LQkICbm5uJCUl4eHhQe/evXnuuecq29PW1iYwMJBTp06xfPlyNmzYQHFxMYWFhaxatYrx48eTnJyMgYEBixcvxsLCglGjRpGTk4OHhwerV6/m2WefveO3jx9++IEFCxZQUlKCn58fCxcurPfoI/lXWhNU2tejfdyvPeUVBS6kapz/mQTISoBj62D/NQVpla7mm0A7L7D00rybO9ZJ3oC2SosI93YMcGtLbFIeX0Un88mmkyyMSnrg4uvGrfQJHe2ER+/27Fmbwr4NqRyOysS7rw0u3S3RruY3CMnD4+yHH3L1WO3KI+s5dcbi9dfveU6Frn1xcTFnzpyplC6+H/r6+uzcuROARYsWUVZWxsiRI3FxceGNN95g8eLFmJiYEBcXx9WrVwkKCiIsLAyAvXv3cvjwYezs7EhNTeXw4cM3SS1XcOXKFbZs2cJ7773HuXPniI2N5eDBg7Rq1YoXXngBT09P1qxZw9atWxkzZgwJCQksWbLkpgfVnTh27BiRkZHs2rULHR0dpk2bxooVKx5aoZO70SAdfqPOtBUCWtlpXs6DNPsUBfLTIOsfzStzPxz6BfZ9ozmuY6gZ/Vt5g6U3WPlC83a1aNL14uvHz15kcXQyP+y+Xnx9crADrlbVi8JpaWFIn0kuZKddZPeaJHatOsWBLen49rejs7+FTN6SVHLjlE5sbCxjxozh8OHDdzy3bdu2lfP3t06VTJ48mWHDhvHGG28AsHHjRg4ePFhZErCgoIDExER0dXXp2rUrdnZ2d7WpYsQvhGDgwIH07duX5cuX07t378oC6jt37mT16tUA9OjRg7y8PKoq6Lhlyxbi4+Px9fUFNA+9hhB51CAdfr3O4dcFQkBLW83LebBmX3k55J2CzHjI2g8Z+yB2IZRfUwc0bgdWPhrnb91Vsy6gU/OF0s4WzfnPMA9mhDuyfFdqZfH1AHtTJnW3J6RT9Yqvt7ZpTsSLnmQcP0/smmS2fX+chE2n8Yuwx96zdgq5S2qH+43EHwYBAQHk5uaSk5ODtrb2TWUAb5UhvlXaODAwkG3btvHyyy+jr6+Poih8/vnnhIeH33ReVFTUbdfeSsUc/q3cT964qr/PiqLwzDPP8NFHH1Xp/IeFHIbVF1paYN4JPJ6GfvNg0jZ4PRMmbIE+/6dZKzhzADa9BUvD4SNr+Lon/DULjvwGF2tW7qytSTNm9XNi16wevN6vMym5hYxbFkefz3awKj6j2sXXrTq34smZ3vSdrJG2/WvxYRZXcC8AACAASURBVFZ9vI/0Y+drZKekaXH8+HHUajWmpqbY2Nhw9OhRrl69SkFBAVu2bKGgoOCuBcGfffZZ+vXrx9ChQykrKyM8PJz//e9/lRLKJ0+epLCw8LbrjI2NuXTpUrVtDQ4OZsWKFYDmIWJmZnabfv/d6NmzJ6tWrarMKzh//jxpaWnVtqG2aZAj/EcWbb1ro3ofYIpm3+VsyIiD9D2QHgf7lsLuhZpjLdqDtR+0D9C8zDtrHiTVoLm+DpOCHRgbaMe6A1ksjk7mlV8OMP/vE4wLsuVpv/Y019epUltCCOw9zbF1N+PE7rPsXZfM2v8mYNW5Jf6DHGhjW/OSjpLGR8UcPmhGvt9++y0qlQpra2uGDRuGm5sbHTt2xNPTkytXrtxz2mT69OkUFBQwevRoVqxYQWpqKl5eXiiKgrm5OWvWrLntGlNTU4KCgnBxcaFv3743Ldrei9mzZzNu3Djc3NwwMDDg22+/rfI9d+nShTlz5hAWFkZ5eTk6Ojp8+eWX2NjYVLmNuqBBxuFX8CjKI9+XshJNSGj6bs1D4PRuuHxOc0zf5LrztwnSRAepquasK1AUhe0nc1gcnUxMUh5GetqM8GvPuCBb2ppUb0qprFTN4e2ZxP+VRvHlUhw8zfEbaE9LiwfPBJZUj4YQh18dZOJV9ZGJV48SFVFBp2MhLUbzAMi79pnpGGjWAGyCwPYxzTeHauQFHMooYPGOZDYczEJLCCI82jE52AFHi+pp85cUlZGwJZ2ETacpK1HTObAtvo/bYdzqwYu7SKqGdPhNn0bv8CuQI/wH5HL29QdA6i5NtjCKJjnMyhdsu4FdN7D0qVJOQPr5K3yzM4XIuHSKStWEOJozKdieAPvqFV8vulRC/J9pHIrOQCBwDbHEu48t+kbV+xYiqTrS4Td9pMOX3EzRBUiLhdSdkLpDMyWEAtrNoL0/2AWDfXdNctg98gHyr5Tww+40lsekknu5BFdLEyYF29PXxQLtaoRhXswrIm5dCif2nEVHT4VH7/a497RGV18uKdU2jc3hnzx5EoBOnTrVsyWNB+nwJffmynnN6D8lWvMAyD6q2a9nohn523UH+xCNcNwdRvDFpWp+3Z/Jkh3JJOcWYt2qGc8G2THM1xoD3ao77bysy+z5PZmUA7k0M9bBp58tzo9ZotKRwWO1RWNz+JLqIx2+pHpcztY4/5TtkBylUREFaG6lcfwOoWAfCoamN11WXq6w6dg5FkcnE592gRYGOoz2t+GZQFvMqqHNfza5gN1rksg8mY+xqT5+A+zo2NXigQXfJNeRDr/p0yAdvhBiEPA40Br4UlGUjfe7Rjr8euJ8ssbxJ23TPASKCwChifpx6AkdemrWAm6IAIpPO89X25PZdOwcOiotnqymNr+iKKQfO0/sb0nkpl+mVTtD/AfaY+tmJpO3akBjc/hZWVkAtGtXe5nmTZ3qOPwafXcWQiwVQmQLIQ7fsr+PEOKEEOKUEOI1AEVR1iiKMhEYCzQeeblHkVb24DMehn8Pr6ZoksFCXweVHuz8FJb1hbn2mlKS8d/CxSy8bVqxeIwPm6d35wkvS1bFZ9Djkygmf7+P+LQL9+1SCEH7LqYMm+VL2ARn1KXl/PG/Q/w6bz9ZifkP4aYldYUQgtGjR1dul5WVYW5uTv/+/W8799KlSw+UJHUnli9fjrm5OR4eHnTp0oWvv/66WtcfP34cDw8PPD09SUpKqpR0Tk1NxcXFpVZsfNjUdJVsOfAF8F3FDiGECvgS6A1kAHFCiLWKolybKObNa8cljQEt1fVksO6vamSiU6Lh1CY4tUUjDAcaFdCOvXHoGM5HA32Z3tuxUpv/7yPn8Lmmzd/rPtr8QkvQ0acN9p7mHNt1hrgNKfz2yX5sXEzxH2SPmVX1QkIl9Y+hoSGHDx+mqKiIZs2asWnTJiwtLWul7XvJGINGj+eLL74gOzsbZ2dnIiIiaNOmTZWuX7NmDQMHDuTdd98FICYmplZsrk9qNMJXFCUauDV3vitwSlGUZEVRSoCfgIFCw/8BfyqKsv9ubQohJgkh9gkh9uXk5NTEPEld0KwFdImAiM/hpSMwNRZ6v6epCRzzOSzrA/McMN/4HK+0O0TMi168M6ALZwqKmfR9PL0+3c5Pe09TXHpvjXOVSguXYEtGvR9AwGAHziYXEPlBHBu/OUJBzpV7XitpePTt25cNGzYAsHLlSp5++unKY+fPn2fQoEG4ubkxfPjwytDMG/f7+/tz8OBBQJMBO2nSJMLCwhgzZgw5OTk88cQT+Pr64uvry65du27rv3Xr1jg4OJCWlsbYsWOZPn06oaGhzJw5k4SEBPz9/XFzc2Pw4MFcuHCBP/74g88++4wlS5YQGhoKgJHR7dXe1Go1M2bMwNfXFzc3N7766qta/+xqk7qIg7ME0m/YzgD8gBeAXoCJEKKDoiiL7nSxoiiLhRBngAG6urredWCfpLYQQlMHuE0XCHpRM9eftBVO/g2JG+HQLxgKFeNsAhnTLZxtijef7lfz2q+H+GTTScYG2jLKzwYTg7vH4evoqvAKt6HLY+34Z+NpDm5NJyk+my7d2uHTzxZDk4Zbc7ihsePnk+SmX77/idXAzNqIbsPuH0L51FNP8d5779G/f38OHjzI+PHj2bFjBwDvvPNOpQzxt99+y8yZM4mIiLhp/43yxADx8fHs3LmTZs2aMWLECF566SUee+wxTp8+TXh4OMeOHbup/+TkZJKTk6lQ4D158iSbN29GpVLh5ubG559/Tvfu3Xn77bd59913+eyzz5gyZQpGRka88sord72vb7755o4SzfdS6qxP6sLh3+n7uqIoygJgQVUaaHJqmY8K+iYaNVDnwVCu1shAn/gDTv6FatOb9AJ6mncmwzuUb3I6M//vIhZuO8Vw3/Y8280OyxZ3l27QN9QhYLADbqFWxP2RypEdWRyPPYN7T2s8w2zQkwVYGjRubm6kpqaycuVK+vXrd9OxG2WIg4KCyM/Pp6Cg4J7yxBERETRrpvl92bx5M0ePHq1s7+LFi5XrAJGRkezcuRM9PT2++uqrSunjoUOHolKpKCgoID8/n+7duwPwzDPPMHTo0Crf190kmh8lh58BWN+wbQVkVaeBRq2HL9GgpQJrX82r1zsa+YcTfyFObMD66GJmK2reaNWGPbr+LN7dhR6xXejrZs2kYAe6tLu7yJphCz1CRjji0dOaveuSif8zjcPRmXiH2+IaYom2rizAcjeqMhKvSyIiInjllVeIiooiLy+vcv+NkYIdOnRAW1sbIcQ95YlvlDEuLy8nNja28gFwIxVz+LdyP/nkqnI3ieaGSl1kuMQBHYUQdkIIXeApYG0d9CNpTLS0Bf8p8Mw6mHEKBi9Gx6YrjxVu4judj0hoNo2QY+/wny8+Y/ySHew6lXvHP/gKWrQxIGyCC8Ne96WNTXNifj3Find2c3RXFuXq6kk7Sx4O48eP5+2338bV1fWm/XeTIa6qPHFYWNhNTv1OOvd3w8TEhJYtW1ZOL33//feVo/2qUFWJ5oZCjUb4QoiVQAhgJoTIAN5RFOUbIcTzwN+ACliqKMqR6rQrp3SaOAatNCUh3YdDyRVI2kqzY+sYeOIPBontFGZ8yZZvPfnQJAT3kCfp42l/V+kG8/bGDPiXBxknLrB7TdL1AiwD7bH3kAVYGhJWVla8+OKLt+2/UYZYW1ubuXPn3rb/XvLECxYs4LnnnsPNzY2ysjKCg4NZtOiOS4R35Ntvv2XKlClcuXIFe3t7li1bVuVrJ0yYUCWJ5oZCg8y0lWqZjyhlJZAaTdnh3yk7ug79kgsUKnrEavui5TwY//DhGBjePSxTURRSEnLZ/XsSF85eobVtcwIG2WPVudVDvImGRWNLvJLiadWnQWbaPggy0/YRRl1GecpOMmN+pHnKn5goFylEn9NmIbTrNgoT5/C7Kn2Wq8s5vvsscetTuHzhKtZdWhEwyAHz9o9eDL90+E2f6jj8BhnaIBdtJai00eoQgnWHEFCXcXLvH5yLWYlrznZMfvuLK2uNUTsOwNjnKY3W/w0qn1oqLboEtaNT1zYcisok/q9Ufv4wjg4+rfGLsKdFa4P6uy+JpB6RI3xJoyLp7Hmi//yZVinr6Cn2YSSKKTVojY7bUHAbpin2fsu8/dWiMv7ZmMaBLemUlyk4PdYO38cfjRh+OcJv+jT6Eb5EcjccLFrhMG4K2RfHsnjHcbLifiPs0g567P4K7d1foph1QrgN1zj/Fu0B0Gumjf9AB1xDrNj3RypHd2RxYvcZ3HtY4xkuY/gbEjo6shhOXdIgR/hy0VZSVS5fLeOnvaf5ZcdBvAqjGaEfi6v6WhKOTRC4Pw1dBoL+9XC+/Owr7F2XQmLcOfQMtTUx/KGWaOs0vRj+xjbCl1QfuWgreeQoVZezNiGLxdHJFGYnMdpgD0/p7sKk6LSmspfTAPAYoSnuoqUJ8cw5fYnda5I4ffQ8Ri318O1vR2d/C7SqUb2roSMdftNHOnzJI4uiKESdzGFRVBJ7UvJ4TD+FV9rsxy1/C1pXCzRFXTyeBo+R0EqT/p5x4gKxvyWRnXqRlhYG+A9ywM69aejwNwSHb2RkxOXLVdPwOX1aU3ynffv2dz1n9uzZ99W4qW6/VWHt2rUcPXqU1157jTVr1tCpUye6dOlSa+0/KHIOX/LIIoQg1LE1oY6tSUjPZ3F0WwYftsdQNYjX7VMYxDYMoudD9DxNMXfP0Vh1ieDJmd4k/5PD7t+T+XPRISzsmxMw2IF2HVvW9y09UhQVFdW3CXclIiKCiIgIQCOd3L9//wbh8KtDg/zuKoQYIIRYXCGUJJE8CB7WLVg40putL4cQ4W3P7GRHnBMnMctmJWe8X4GCdPhtEnziiPjzVRzaZfP0210JHdWZS3nF/PbJP6z/8gB5mbWrMPmooigKM2bMwMXFBVdXVyIjIyuPzZ07F1dXVwYOHMgnn3wCwNdff42vry/u7u488cQTXLlyb1nslJQUAgIC8PX15a233rrp2Lx58yoljN955x1AU8jEycmJiRMn4uzsTFhYWOUDZ8GCBXTp0gU3NzeeeuopQFNQ5fnnnycmJoa1a9cyY8YMPDw8SEpKwsvLq7KvxMREvL0bptBvgxzhS2kFSW1iZ2bIh4NdealXJ5bHpPB9bBorT3gRYNeDV0Nz8chZi4hfDnsXo2XpQxfvsXR8ayCHdl1g/99p/DRnL45dLeg6wI7mZndX9GwMhISE3LZv2LBhTJs2jStXrtymZAkwduxYxo4dS25uLk8++eRNx6Kioqrc96+//kpCQgIHDhwgNzcXX19fgoODSUhIYM2aNezZs4f09HTy8zUVzoYMGcLEiRoX8Oabb/LNN9/wwgsv3LX9F198kalTpzJmzBi+/PJ6jaWNGzeSmJjI3r17URSFiIgIoqOjad++PYmJiaxcuZKvv/6aYcOGsXr1akaNGsXHH39MSkoKenp6lfZUEBgYSEREBP3796/8PExMTEhISMDDw4Nly5YxduzYKn8uD5MGOcKXSOoCc2M9ZoR3JmZWT9583InU88UM/lObvhnPsCEsCnXYh3D1Eqx9Hp3PnfEq+5xR/zLFs1d7TsVns2L2bnb+nEjR5ZL6vpVGyc6dO3n66adRqVS0adOG7t27ExcXx+bNmxk3bhwGBpqEuBYtWgBw+PBhunXrhqurKytWrODIkXtLcu3atauysMqNJRU3btzIxo0b8fT0xMvLi+PHj1MR/WdnZ4eHhwcA3t7epKamAho555EjR/LDDz/cs6JWBRMmTGDZsmWo1WoiIyMZMWJE9T6ch0SDHOFLJHWJkZ42E7rZMybAlrUHslgcncRzv6Vi2cKJZ4N+YkS7LPQPfAf7v0U/7msC2wfiOmIicYmOHNyWztGYLLzCbHDvaY2OXuMK5bzXiNzAwOCex83MzKo1or+VuwWIKIpSuUCup3c9GW7s2LGsWbMGd3d3li9fXqW+77TQrigKs2bNYvLkyTftT01Nvak/lUpVOaWzYcMGoqOjWbt2Le+///59HzZPPPEE7777Lj169MDb2xtTU9P72lofNMgRvpzDlzwMdLW1eNLbir9eDOabZ3xo10Kf9zYcw3/FFf5j/AoXphyE3u/DpTMYbxxHj+yBPBW2Hyv7ZuxZm8wPb8VyODoTtZRjrhLBwcFERkaiVqvJyckhOjqarl27EhYWxtKlS7ly5Qq2traVEsiXLl2ibdu2lJaWVsok34ugoCB++ukngJvODw8PZ+nSpZURO5mZmWRnZ9+1nfLyctLT0wkNDWXu3Lnk5+ffFu1jbGx8U7F1fX19wsPDmTp1KuPGjav6h/KQaZAOX1GUdYqiTDIxMalvUySPAFpagp5ObfhlSiCrpwbga9uKBVsSCViQwOy8nmSM3gkjV4OlD60OzKHfhT4M8d6ASfNStv94gp/e28up+Ox76vdLYPDgwbi5ueHu7k6PHj2YO3cuFhYW9OnTh4iICHx8fPDw8GD+/PkAvP/++/j5+dG7d286d+583/b/+9//8uWXX+Lr68uNg8WwsDBGjBhBQEAArq6uPPnkkzc561tRq9WMGjUKV1dXPD09eemllyqnmSp46qmnmDdvHp6eniQlJQEwcuRIhBCEhYU9yMfzUJBx+BLJHTiVfYmvtiezJiGTcgUGuLVlSogDnfUuQNwS2P8dSlE+qc2eJLZgGBcu6NDatjmBQxyw7NRwQjkbQhx+daiYQ7e1ta1XOx6E+fPnU1BQwPvvv/9Q+5Vx+BJJDenQ2ph5Q92ZHtaJb3ak8OPe06xJyCLU0ZypIS/hG/Ia4tAq7PYswqboKU6Y92fvuRGs+c9FbFxNCRjkgKmlUX3fRqPj6tWr9W3CAzF48GCSkpLYunVrfZtyT+QIXyKpAvlXSvg+No3lMankFZbgbdOSKd0d6OlojtbpnbB7EWXHN3Ow6HH2XxnO1TJdOgdY0HWAPcat9OvN7sY2wpdqmdWnQY7whRD2wBuAiaIoT97vfImkIdHCQJcXenZkQjd7folPZ3F0MhO/20fH1kZM6e5AxLAf0MlPwWvPIrrEP8/+/L4c3D2AxL1ncA1tj3dfW/QNpRKkpH6p0aKtEGKpECJbCHH4lv19hBAnhBCnhBCvASiKkqwoyrM16U8iqW+a6aoYE2BL1CshfDbcA5WW4OVfDhAyL4rlx7Uo6vUx+q/EETjYnpE279BRZxsJm9P4YdZ29v+VTFmp+qHb3JC/xUtqRnX/b2s0pSOECAYuA98piuJybZ8KOAn0BjKAOOBpRVGOXju+qqojfDmlI2noKIrCthPZLNyWxL60C7Qy1GVcoC1jAmwx0VXg8Gpyt0QSmxbE6RJvjAxK8BvsSKcgW7S06l6cLSUlBWNjY0xNTRuFGFxVxNMkGhRFIS8vj0uXLmFnZ3fTsTpTyxRC2ALrb3D4AcBsRVHCr23PumbcR9e2pcOXNEniUs/zv6gkth7PxlBXxSh/G559zI7WxnpwajMZf/xKTKInOWUdMDUpJHCYC9ZetnXqiEtLS8nIyKC4uLjO+pDUH/r6+lhZWd1WOOZhOvwngT6Koky4tj0a8APeAT5AM/JfUvEAuEN7k4BJAO3bt/dOS0urkX0SycPmaNZFFm1PYv3BLLRVmuSuKcEOtDc1QEnfx6nffmf3CScuqi2wMs8j8GlPzLvI+s2S2uNhOvyhQPgtDr+roih3Vz26vU1Z8UrS6EnLK+Sr6GRW7cugrLycAe7tmBriQGeL5qjPnOBw5Ab2nbSnuLw5ndqdxm9kIM0dOtW32fXKqFGjAPjhhx/q2ZLGTYOZ0qkOckpH0hQ4d7GYb3amsGJ3GoUlano5teG5UAc827fk6tkU9v+wmQOnrFEQuNkk4T2qJ/rWj2ZYYoWaZ000eyQP1+Fro1m07Qlkolm0HaEoyr3Vh25uU47wJU2O/CslfBuTxrKYFPKvlBLoYMrzoR0IcDDlckYGe7/fzvHTFuiJQnw6nMR1xABUbe8vKdCUkA6/dribw69pWOZKIBZwFEJkCCGeVRSlDHge+Bs4BvxcHWcvkTRVWhjo8mKvjuya2YPX+3UmMfsyI5bsYcj/YthzSZces0YyfHonWpuXsivRmx/n7Cfxi3dRso/Xt+mSJoLMtJVI6oniUjWr4jNYtD2JjAtFOLVtznOhDvR1aUtmQhoxkQfJKzCitc5JgtxP027QeDBr2ou7coRfO9R7pm11uGFKp75NkUjqDH0dTejmcF9r1iZksTDqFM//+A/2ZieZFtqBwXMeJ2VXErt/v8pv+zphd+gXArvm06LfNGhpU9/m1wkBAQH1bUKTRo7wJZIGgrpc4a/DZ/li2ymOnbmIdatmTO3egYGuFhzbeJL9m7NQl4GL4UZ8grRp1vvfYGxR32ZLGiB1tmhbF8hFW8mjjKIobDmWzedbEzmQUUBbE32mdHdggGMbDq49ytE9+eiIK/g0X4NbqDWq7i9Cs4YjySypfxqVw69AjvAljzKKorAjMZfPtyYSl3oBc2M9Jgfb06+9Gft/PULaiSKaq84S0HI1DmFBCP8poGtQ32bXiCeeeAKA1atX17MljZtGNYcvkUg09VmDO5nTraMZu5PPs2BLInM2HGORkS4Tu9kT1qM58atV/J39HBa/HCNoxxAs+jwNHiNB1Tj/tPPy8urbhCZNg/ytkIu2Esl1hBAEOJgS4GBKXKrG8X/053G+MtRlQogdgVp6/LMOVqc70fH7aAKi+2Hcbzp0CodGIJgmeXjIKR2JpBESn3aB/25JJPpkDi0NdJgYYIfrRTiy9TSUl+FhsAYvpzPo9n0bLL3q29wqI8Myawc5pSORNCG8bVry3fiu/HNa4/jnbjmpcfz9bLDNKCb+n6Ec+ycf/5Pv09nPHNHrbTCxqm+zJfVMg3T4ckpHIqkanu1bsnycxvF/tjmRuTtO0cpQl0lhlhgeN2Dr6Rc4FJXMYwkjaBcaBkH/Br2GW2u3Z8+e9W1Ck0ZO6UgkTYj9py/w6aaT7EjMxcxQl0m2FugeOk/hRTUd9HcSYPEXzfu8AG7DQatGyiqSBowMy5RIHiH2pZ7nP5tOEpOURzsjPcabmlJ25AKoS/E0+A2vDsno9P8ArH3r21RJHSAdvkTyCBKTlMt/Np5kX9oFOho1Y7iuEVeTL2OkfZ4Aw2V07GqFCHu3wWTs9u3bF4A///yzni1p3NSJWqZEImnYBDqY8cuUAJaP80XfRJc553OIshKUmLRhU8HL/LbdnZxPhkDM56AurW9zKSoqoqioqL7NaLI0SIcvhBgghFhcUFBQ36ZIJI0eIQQhjq1Z+3wQi0Z5U2CoxYdlFznUVkWOlhM/n32fqF8yKPq8NyRH1be5kjpETulIJI8Y6nKF3xMy+XTzSbJzixiiZ4R1bhm64gp+Rj/g7GOIVp8PoHnbh26bjMOvHWQcvkQiAUClJRjiZUV/t3b8FHeaBVtOoRiVMEwYE31xMkd3pBB85Gna9hkKXSc3WpkGye00yf/J4rJiypXy+jZDUk+Iu8gJCMQdzxHX/l3bqPy58l1cO6OJyRToamsxJsCWJ7ysWLozha+2J2OlLufxUlt+Pfc2jpHbCIjrj+GQOWB122CxTujfv/9D6edRpUlO6UzZPIVdmbvqwCKJRPMg0BJamgeCuL5dsU9LaCGE5l0lVJpjaKGlpdnW1tJGJVSotFRoC210tHTQ1rr2rtK866n00FPpoavSRV+lj55KD31tfQy0DTDQMcBA2wAjXSOMdIww0jXCRNcEEz0T9LX1H/i+zheW8MXWU6yMTcW/SBufYhW6ohg/wxW4dGuLVtjboG9Sex+kpM6o97BMIYQhsBAoAaIURVlxv2se1OFvSttE5qXM6hspaVAoVP93827X3Ph7fus5FccUlJt+rnxXrh27dvxe7+VKOeVKOQoK6nI15Uo5akVdua1Wrr3K1ZSVl1GmlGner71Ky0spUZdQUl5CibqE4rJirqqvUqwupqy87L73r6/Sp6V+S0z1TTFtZoq5gTltDNrQxqAN7YzaYWVsRRuDNmhr3f3Lffr5K8zfeILo+DP0L9HG4qo2ptophLT5BYsnpoJThBRla+DUicMXQiwF+gPZiqK43LC/D/BfQAUsURTlYyHEaCBfUZR1QohIRVGG3699uWgrkVyntLyUorIiikqLKCwrpLCkkEull7hccpmCkgIKrmpe54vPk1eUR15xHtlXsjlffP6mdrSFNlbGVji0cMDexJ5OLTvhZOqEtbE1WuJ64N7BjHw+3HCU8ycK6HtVha5amy7NNhLgcQb9QR+CiWWt36NctK0d6mrRdjnwBfDdDR2pgC+B3kAGECeEWAtYAYeunaauYb8SySOHjpYOOro6NNdtXq3rStQlnLtyjqzLWWRcyiD9UjppF9NIKkgiKj0KtaL5czTUMcTZ1BmP1h54tfbCvbU7KycFsOVYNnM3HMUy4yoU9SJ5TyGBJ16m86BeCJ/xUqKhEVEjh68oSrQQwvaW3V2BU4qiJAMIIX4CBqJx/lZAAg00/l8iaYroqnSxNrbG2tgav7Z+Nx0rVZeSVJDEsbxjHMk7wqHcQyw5tIRypRyVUOFu7k5gu0DmPRPEP4mGfPtXIiHnddh6fhLHvz9C97jRtBr+LphJocPGQI3n8K85/PUVUzpCiCeBPoqiTLi2PRrwA2ai+TZQDOy82xy+EGISMAmgffv23mlpaTWyTyKRVI8rpVc4kHOAvWf3EpMVw9G8owC0NWxLt3ah5J51JG1nM3oWaaGtCDyN1+IzoBPaQdNqHMIpp3Rqh4cZh3+n1RxFUZRCYNz9LlYUZbEQ4gwwQFdX17vWrZNIJPfEQMeAgHYBBLQL4EWvFzlffJ7ojGi2pG1hTdIqSspLsAiwZOOFADoccaP80hBORWYREjcRq9GvQmun+r4FyV2oC4efAVjfsG0FZNVBP3dFKS2F+gw3fRgRDFXpoybn3OPaRraCsQAAFuxJREFUphaPLrk3rfRbMajDIAZ1GMTlkstsTd/K2lNr2XNlFWleq0nM6U1ISjC/J46m88crCHq8NfohDzbaHzZsWB3cgaSCupjS0QZOAj2BTCAOGKEoypHqtv2gUTqnJ06icMeOal8nqSE3Pgiq8bO4cb8Qd/xZ3Lh9yzl3PCauJUpVvLS0NDHzQuvazzcc19ICLXH9mJYmjh6V6vp+lapyn1BpgZYKVFoIlXbltlCpQFul2aetunauNkJbW7NfW4f/b+/Ow6MqrweOf092kpAQCAFJgiCEaCg2YFChLihKgbpLq6jIJquC4kJZrGKLW3FhUUEUilIX0J8im6Ko1CpoRYiyk4BAIgYQDGRfZt7fH0kwxkwImeXOZM7nee6TmXfucg6XObm59733laAgJDgYCa76WW0KCamcQpHQEAJCQ5HQUCQklICwUKRJEwLCwiravOBC6cH8g7yT8Q5v7X6L4wV5pB24jtScHoQH5HHRmR+TNHQ8EpdsdZh+yV3dMt8AegGxwCHgYWPMAhHpD8ykolvmQmPMo6e53qoRr0ZkZGScdlzHV62iLCv7tJdzDQ/8ZVGPfVav/eponroWrb7Mr5av1s/d0Ty/mt1U+/nL65PLmhrzVM5njPllPVWf2+3VPjO/LGvsv7TZzcn1GLu92mvbL/PbbRXz2w3YbBhjr3htt2Gq2ux2sNsrlrNV/iy3YWy2is9tNkx5+a9fl5X9qt0VJCyMgCZNCAgPJyC8CQHhEQREVE6RkQQ0jSQwsikBUU0JjIomMDqKwGbNKqaYGAKjoyt+QblAqa2UNfvWsHDLK/x8MJ/LMm8htjCBxLCN9L66GRGX3VHvnjyFhYUAhIeHuyQ2f2X5jVcNof3wVWNj7PZffgmUlWHKyyt+lpZiSkuxl5SefG1KS7AXF2NKSjElxdiLik/+tBcVYS8qxBQWYi8sxF5QgK2gAHtBAfb8Aux5edgLChwHIlJR+JvHENQilqC4OIJatiS4VRxBrc8guHUrgtu0ITA2tt6n8IwxfHHwC575ch7hu5vT/UB/gqWMHomfkTrqTqR5u1OuQy/auoZPFXxnj/CVUmBsNuz5+dhOnMB2/AS247nYfs7FlpuL7dgxyo8dxXb0KOVHfqL8yBHKjxzBlJb+ah0SFkZwQjwhCYmEtGtXMbVvT2inJIJiYhxu+/PsDfxzzXy67OxBm7yOhDfZyg3XxxF98U11Xh/Sgu8aPlXwq+gRvlKeY4zBlptLeU4OZT/mUHbwIGXZ2ZT9kE3p/gOU7t+PKSk5OX9gbCyhSR0JS0mhSefOhHXuTHDbtif/IjDG8ObW93l/1Sd0O9CbAAzJ7dbT584JSESLWmPQgu8aPlXw9QhfKe9j7HbKc3Io2bOXksxMSjIyKNm5k5KMjIqecUBg8+Y06dqV8K6phF9wAWEpKZQZG9NWvIxsCKDNiU7YInYw4KazaHN+v99sQwu+a/hUwa+iR/hKeT9TWkpJZiZFW7ZStHkzRZs3U1p5w2RAVBTh53cn8tJLOXbuOcxctowO+88DsdOq4zfcfNdEJPiXJ3xqwXcNLfhKKY8pP3KEgq/+R8GXGyhYv57ygz8CENalC1s7JfNtQVtiijuQH7WFm4ZcQLuUCwFYtGgRAEOGDLEo8sbBpwq+ntJRqvEwxlCyezf5n35K3sefULxlCwZhW+e+5LT8I2VBxcSem8HtI+7Txy67iE8V/Cp6hK9U41N64AAnVr/P8ZUrOPZjIdvOGUJBZALFEf+j39D+xLVuQ2xsrNVh+jQt+Eopr2KMoSg9nUNvLOG77ZAdfwWzlo8nNLyE7zIyvOJuYl/lqOB75b+oiFwtIvOPHz9udShKKTcREcK7dqX9P5/gqten0Lbpx4gRSooieWfARH764MP63TGu6s0rC74xZoUxZmR0tI6fqZQ/CIyK4uqnHifmzAiM5JMT1581/9pL+nUDyf/v51aH12h4ZcFXSvmnsCbhxHdMwNbpC47FnMlXLW9h06TZHBg5ipK9e60Oz+dpwVdKeRcRxt/7N7r8uYBjTY+xpctoNubEkXntDRx6/Im6nxGk6uSO5+ErpVSDjBkz5uTrXr0HcO65+5kz51XgMg617MD5SxeQ99FHtJ72MJGXXGJdoD7KK3vpaD98pVQVY7Px/LyHKN3egyA7nH1gOYl7/kvU1VfT+sGpBOq1vt/QbplKKa+XlZUFQGJi4m8++/jjBaxfFUjzwraElG+k51evEdKiOW2eeIKICy/4zfz+zKe6ZSql/NOgQYMYNGhQrZ/17j2coRO6kB33OaVBaXx46QTyQptxYOhQDs2YcfIBbsoxLfhKKZ+RkHge06bczZGkdxF7K/7TcTD7Lh7AsQUL2T94CGWHDlkdolfTgq+U8imhYdE8PGEWkRd/RmFwMZkBF/NZ37so2rGD72+4kYIvv7Q6RK/lsYIvImeJyAIRedtT21RKNU4SEMAdAx/nopsL+SFmG+XF57D60rsoj2rOgeF3cGzxv/Uu3VrUq+CLyEIROSwiW2u09xWRXSKSKSKT6lqHMWavMWa4M8EqpVR1l144jNF39WR3/GrCitvxQdJA8rpfzqFHHyVn2iN6Xr+G+vbDXwQ8B7xa1SAigcDzwJVANvC1iCwHAoHHayw/zBhz2OlolVKN2n333Xfay7SLP5+/39uOh1+azJkZN7Ih4go69onnrCWLKN23j4Q5swmMinJDtL6n3t0yRaQdsNIY87vK9z2AacaYP1a+nwxgjKlZ7Guu521jzIA6Ph8JjARo27btefsrR85RSqm6lNlKeeL1Own4pjdRJc2JCP+eHmueI7R9exJfmk9w69ZWh+gx7uiWGQ9kVXufXdnmKIAWIjIP6Fr1y6E2xpj5xpg0Y0xay5YtnQhPKeVrdu3axa5duxq0bHBgCA/eNp+E/js5GJVJUUESH/aeQMnBH9l380BK9CZOpwp+bUPTOPxzwRhz1Bgz2hjToR5/BejjkZXyQ6NGjWLUqFENXl5EGNxnGn1uiWJ7q3UElbRjTY/xFNsC2HfrbRRu3uzCaH2PMwU/G6h+O1wCcNC5cJRSynmX/f42xo3sz6Z2bxFY3ppPzh1JXrMEDgy/w6+7bTpT8L8GkkSkvYiEADcDy10RlD4PXynlrLPjL+TRsZNJT3kFO5Fs6Hg7R87oQtbIUeR9+qnV4Vmivt0y3wA2AMkiki0iw40x5cBdwBpgB7DUGLPNFUHpKR2llCu0jmrL7NHz2NvtNQqCS/m27U0caH8x2ePGc+LDD60Oz+P04WlKKa/Rq1cvANatW+fS9Zbby3n0nTuxb+hBy4IEEg+vpdPuVcQ/+wxRV17p0m15A0e9dLzyefjVHo9sdShKKQ968MEH3bLeoIAgHrpxHi/HPMLutScQ+lASGIGZcC8JjbTo10aP8JVSfmXl5oV88G425xy+iGa56aRue4W2s56l6eWXWx2ay/jU45H1HL5S/ik9PZ309HS3buOqrsMYNPhCvk1YRW6zVDb+fgz7JkykYP16t27XG+gRvlLKa7jrHH5tMg5/x+OLnqfr3gFEFmTTbccCkl6cRfh557l92+6mR/hKKVVNUty5PD7mYb47+zXyItuwsfNYdo59gKItW0+9sI/yyoKv/fCVUp7Qqmkbnhkzm4wub3IiIppvUsaybexfKfn+e6tDcwuvLPhKKeUpTUOjeGbUfA6lLedEeBCbOo1g84iJjXL0LC34Sim/FxwYzGND51J+yeecCLOR3m4wXw2biC031+rQXMorL9pW64c/IkOfcKeU31hf2VOmZ8+elsXw8spHyX2/AxGl0Zx97D0ufW0mAWFhlsXTEI4u2nplwa+ivXSUUlZ4d9089i6NJLwslg5Fa7jylaeQwECrw6o3n+qlo5TyT+vXrz95lG+l63uNptvQMEoCD5EZ3o9Vo6c2ijFyteArpbzGlClTmDJlitVhAHBZ9wFcMi4Ru8nmgPTm7bsfsjokp2nBV0opB7p1vpwrJ6YSUL6PIyWX8OqkaVaH5BSvLPh645VSylskJ3Wj/9SehBZ9T97PFzH3kYetDqnBvLLg641XSilv0rZDZ/409SIi8/dgP3gxMx6b4pPn9L2y4CullLdpnZxCv/vOJyovk/D9lzPt6YmU28utDuu0aMFXSnmNmTNnMnPmTKvDcCiuW1euGHY2USf2EpfRh7/Ouo+i8iKrw6o3LfhKKa+RmppKamqq1WHU6YwretGrX1OiTuyn/c6rmPDcPRwv8Y3rjR4r+CJynYi8JCLviUgfT21XKeU71q5dy9q1a60O45QSbx/IHzofIzrvAL/beQPj599NTkGO1WGdUn0HMV8oIodFZGuN9r4isktEMkVkUl3rMMYsM8aMAIYANzU4YqVUozV9+nSmT59udRj10n7K/XSPTKdpXg5p225m/MIJ7MndY3VYdarvEf4ioG/1BhEJBJ4H+gEpwEARSRGRLiKyssYUV23RByuXU0opnyUBAXR69gm6nVhJZOFRLtp+OxP+/QDph907Ypcz6lXwjTGfAcdqNJ8PZBpj9hpjSoE3gWuNMVuMMVfVmA5LhSeB940xmxxtS0RGishGEdl45MiRhuallFJuFxARQfLcmXTd8woRRce5fPswJr49lf9k/cfq0GrlzDn8eCCr2vvsyjZHxgFXAANEZLSjmYwx84FHgE0hISFOhKeUUu4XkhBPp1mP0W3LC0SUlNBv+wgeWvUPlmUuszq033Cm4EstbQ7vRDDGzDbGnGeMGW2MmVfXivXGK6WULwlPS6PdpHF02zSb8BLhum2jefKTZ1iwZYFX3aAV5MSy2UBitfcJwEHnwqlQ7Xn4rlidUspHvPjii1aH0GAxf/4z8Tt2YntvDt+kTeCmrWN4UeZwtPgo96fdT4BY3wvemQi+BpJEpL2IhAA3A8tdE5ZSyh8lJyeTnJxsdRgN1mryJFqltCF1y4uElMZw27YRLPnuLSb/dzJltjKrw6t3t8w3gA1Asohki8hwY0w5cBewBtgBLDXGbHNFUHpKRyn/tGLFClasWGF1GA0mwcHEz5pJy7ACUva+SWBxAsO3DmVN5keM+2QchWWF1sbnTeeXqugQh0r5p169egGwbt06S+NwVvGOHewbeAuZHa9kf2wfQiPSea7Lv0lpmcLzvZ8nJizGrdv3qRGv9AhfKeXLws45hzOmT6fD1hU0K/2WkoJu3L/tGnYf3cnt79/OwXyXXO48bV5Z8JVSytdFX/Unmg8bRrf1L0NgFsfzejF1+8UcLTzMoNWDyPjZ82cvvLLg6wAoSqnGIO7eCUT07Mllnz1NfmguB3Ov5ZHtKWAvY/AHg9l0yOE9qG7hlQVfT+kopRoDCQqizdNPEdKyJdekzyE3tJDdR2/lsZ2taRHclJEfjWRd1jqPxeOVBV+P8JXyT4sXL2bx4sVWh+FSQTExJDw3B8n9mesPv8Xx4FI25QzjyV1hJEWdyT2f3sO7Ge96JBavLPh6hK+Uf0pMTCQxMfHUM/qYsJQUzvj7I8jm/9E3Mp38QMMX2cOZsbuYC+K68dD6hzxyV65XFnyllH9asmQJS5YssToMt4i+9lpibr2VJu8upuvZeeQTysffD+LJvT/R78w+zNw0kxkbZ2A3drfFoAVfKeU15s6dy9y5c60Ow21a/XUiTbp2JeFfTxKXFsoJWws+3H4V/9j3A7cmD2Tx9sVM/XwqZXb33JXrlQVfz+ErpRojCQkhfuZMAsLDufCtJyn6fSSHyzuwdnNXHsjax/jUcazcu9Jtd+V6ZcHXc/hKqcYquFUcCc8+Q2lWFgO3vU5W+zD2lXTn8/XR3JGTxbQe09hyZItbbs7yyoKvlFKNWXj37sQ9cD/5a9dyT/Q29rYQthb149u1+7jx2GFW37CajjGuf1qwFnyllLJA88GDadqvL7mzZzHh8kAONLHzRd4Q9ixbRnTGx27Zplc+PK1KWlqa2bhxo9VhKKU85KeffgIgNjbW4kg8w5ZfwL6//AXb8eMUz17A//3rAK1tNq6PnUbrMfOgTWqD1utTD0/Ti7ZK+afY2Fi/KfYAgZERJMyZjb2oiKgZ0+hxaxLHJZBlxx7ieOBZLt+eVxZ8vWirlH9atGgRixYtsjoMjwrt0IE20/9BUXo6qZ++TvP+ieyRYHYfK3H5tpwZ4lAppVyqqtgPGTLE0jg8Lap/fwo3p/Pzq4u56ZlUjjx0CR3jIl2+Ha88wldKKX/T6oH7aZKaSs6DfyMx/7BbtuGxgi8i54jIPBF5W0TGeGq7SinlCypuynoWCQ3lh7vvxl5o0Y1XIrJQRA6LyNYa7X1FZJeIZIrIpLrWYYzZYYwZDfwF+M3VY6WU8nfBrVvT5qkZ2PLyKfvhB5evv75H+IuAvtUbRCQQeB7oB6QAA0UkRUS6iMjKGlNc5TLXAJ8D7ulkqpRSPi7yD3+gw5oPCE1Kcvm663XR1hjzmYi0q9F8PpBpjNkLICJvAtcaYx4HrnKwnuXAchFZBbze0KCVUo3T6tWrrQ7BKwSEhrplvc700okHsqq9zwYucDSziPQCbgBCAYd7VURGAiMB2rZt60R4SilfEx4ebnUIjZozBV9qaXN4264xZh2w7lQrNcbMF5EfgatDQkLOa3B0Simf88ILLwAwduxYiyNpnJzppZMNVB+aJgFwyePd9MYrpfzT0qVLWbp0qdVhNFrOFPyvgSQRaS8iIcDNwHJXBKWPVlBKKderb7fMN4ANQLKIZIvIcGNMOXAXsAbYASw1xmxzX6hKKaWcUd9eOgMdtK+mjguwDWWMWQGsSEtLG+HqdSullL/yykcr6CkdpZRyPa9+Hr6IHAH2A9FA9epf/b2j17HATy4Io+a2nZnX0ed15Xeq976Yc33afDXn+u7j2to0Z//L2dFnzuZ7pjGm5W9ajTFePwHzHb2v4/VGd2zbmXkdfV5XfqeZv0/kXJ82X825vvtYc9ac6/rMVfnWnLzylE4tVtTx3tFrd23bmXkdfV5Xfqd674s516fNV3Ou7z6urU1z9r+cT/Xv4VJefUrHGSKy0dQyxFdjpjn7B8258XNXvr5yhN8Q860OwAKas3/QnBs/t+TbaI/wlVJK/VpjPsJXSilVjRZ8pZTyE1rwlVLKT/hlwffH8XVF5DoReUlE3hORPlbH4wkicpaILBCRt62OxV1EJEJEXqnct7daHY8n+MN+rcll3193dO535wQsBA4DW2u09wV2AZnApHquKwBYYHVOHs45xg9zftvqfNyVOzAIuLry9RKrY/fk/va1/eqinJ36/lqedAP+kS4BulX/RwICgT3AWUAI8C0V4+x2AVbWmOIql7kGWA/cYnVOnsq5crmngW5W5+ThnH2qMJxm7pOB1Mp5Xrc6dk/k7Kv71UU5O/X9dWbEK0sYPxxf1xU5i4gATwDvG2M2uTdi57lqP/ui08mdioGIEoB0fPgU7WnmvN2z0bnH6eQsIjtwwffXZ/+D1FDb+LrxjmYWkV4iMltEXsQNj3f2kNPKGRgHXAEMEJHR7gzMjU53P7cQkXlAVxGZ7O7g3MxR7u8AN4rIXNx8W74Fas25ke3XmhztZ5d8f33uCN8Bt4yv6+VON+fZwGz3heMRp5vzUcBXf7nVVGvuxpgCYKing/EQRzk3pv1ak6OcXfL9bSxH+G4bX9eLac7+kXMVf8xdc3Zxzo2l4LttfF0vpjn7R85V/DF3zdnVOVt9pboBV7bfAH4Eyqj4bTi8sr0/sJuKK9xTrY5Tc9acNXfN2dty1oenKaWUn2gsp3SUUkqdghZ8pZTyE1rwlVLKT2jBV0opP6EFXyml/IQWfKWU8hNa8JVSyk9owVdKKT+hBV8ppfzE/wMz6Ho/7WxIjQAAAABJRU5ErkJggg==\n",
"text/plain": [
"