\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.13?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/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](..\/astro/darkmatter/index.rst). 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",
" DMAnnihilation,\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": [
"The assumed local density is 0.3 GeV / cm3 at a distance to the GC of 8.33 kpc\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd1xUx/qHn9mlVxWwAdJULFQRaRHBgpoolthiSdREYtpNbhIT80szicn1GlNuTDHGqCnGGGNiLCnGgmBAxd4FaQpYEBWVIrB7fn8srFhQEJAF5/Gzn+XMOTPz7oLvmTPzzvcViqIgkUgkkqaPqqENkEgkEsndQTp8iUQiuUeQDl8ikUjuEaTDl0gkknsE6fAlEonkHkE6fIlEIrlHMGpoA26Fvb294urq2tBmSCSSu8TRo0cB8PT0bGBLGjc7d+48qyiKw/XlBu3wXV1d2bFjR0ObIZFI7hIREREAxMbGNqgdjR0hRObNyg3S4QshBgOD27dv39CmSCSSu8hrr73W0CY0aYQh77Tt3r27Ikf4EolEUjOEEDsVRel+fblctJVIJAbDnj172LNnT0Ob0WQxyCkdiURSN5SWlpKVlUVxcXFDm1Itzp8/D8Dhw4cb2JLGgZmZGU5OThgbG1freunwJZImTFZWFtbW1ri6uiKEaGhzbotKpZt0kFE6t0dRFPLy8sjKysLNza1adQxySkcIMVgIMT8/P7+hTZFIGjXFxcXY2dk1CmcvqRlCCOzs7Gr09GaQDl9RlNWKosTY2treUf2TqfnkZV+uY6skksaJdPZNl5r+bg3S4deWf35O4ceZ21m/6BAXzxY1tDkSyT2NEIIXXnhBfzxnzhxmzJgBwIwZM3B0dMTPzw8/Pz8++OADNmzYwNChQ/XX/+c//6FyiPbq1auJjo4GdHt1vL298fX1JSoqilOnTtXItvj4eLp27Yqfnx/Z2dmMGDEC0O0DGDRo0J1+ZIOlSTr8QU/54t+3Hcd2nWHJm1uJ+zGZwoslDW2WRHJPYmpqyi+//MLZs2dvev7f//63Pjrnww8/ZODAgSQmJurPJyYmYmNjw5kzZwBISEggLCxMf37Tpk3s3buX7t278957793QvkajqdK2JUuW8OKLL7Jnzx4cHR35+eef7/RjNgqapMM3szIm9MH2jH87hE6hbTgQl813ryey9bdUrhSVNbR5Esk9hZGRETExMXz00Ue3vdbKygo3NzdsbW05duwYANnZ2Tz44IMkJCQAOocfGhp6Q93w8HB9HSsrK9544w2CgoJITExkw4YN+Pv74+3tzeTJk7ly5QoLFizgp59+4u2332bcuHFkZGTg5eV1Q7sFBQVMnjyZwMBA/P39+e2332rzdTQodzVKRwhhCcQBbyqKsqa++7NqbkrkuE74923HttVp7PwjkwNx2XTr74JPhBNGJur6NkEiMRjeWn2QQzkX67TNLm1teHNw19te99RTT+Hj48NLL710w7mPPvqI77//HtBN8fTt25fQ0FASEhLQaDR06NCB4OBg/vrrLwYNGsS+ffsIDAy8oZ01a9bg7e0N6Jy0l5cXb7/9NsXFxXTo0IENGzbQsWNHHn74Yb744guee+45tmzZwqBBgxgxYgQZGRk3tf3dd9+ld+/eLFy4kAsXLtCjRw/69u2LpaVlDb4pw6BWI3whxEIhxBkhxIHrygcIIY4KIY4JIaZXOvUy8FNt+rwTmrWyoP9jXoz6v0BaudqS+Esq37+xlYPx2Wg12rttjkRyz2FjY8PDDz/MJ598csO5ylM6nTt3Jjs7m7CwMBISEkhISCAkJIQePXqwbds2du/ejaenJ2ZmZvr6kZGR+Pn5cfHiRV555RUA1Go1Dz74IKATZHNzc6Njx44APPLII8TFxVXb9nXr1jFr1iz8/PyIiIiguLiY48eP1+braDBqO8JfDHwKfFtRIIRQA58B/YAsIEkIsQpoCxwCzG5s5u7g0M6awc/4kp18nq0rU4ldcpTdfx8nKNqd9t1aIlQymkHSdKnOSLw+ee655+jWrRuTJk267bWhoaHMnTsXjUbDlClTsLa2pri4mNjY2Gvm70E3h29vb39NmZmZGWq17gm+tvIxiqKwYsWKJrE3oFYjfEVR4oBz1xX3AI4pipKmKEoJ8CMwBIgEgoGxwBQhRIOtHzh2bM7waQHc/4Q3aiMV6xYc5Kf/JJF5MK/WfxwSieTmtGjRglGjRvH111/f9touXbqQk5NDfHw8/v7+APj5+TFv3rybzt/fik6dOpGRkaGf3//uu+/o1atXtev379+fuXPn6n3D7t27a9S/IVEfTtcROFHpOAtwVBTlVUVRngN+AL5SFOWmcylCiBghxA4hxI7c3Nx6ME/fD26+Dox+rQd9J3bmSmEZa+buZeWHuzmVJjd8SST1wQsvvFBltE5lhBAEBQVhb2+vlw0ICQkhLS2txg7fzMyMRYsWMXLkSLy9vVGpVEydOrXa9V9//XVKS0vx8fHBy8uL119/vUb9GxK1VssUQrgCaxRF8So/Hgn0VxTlsfLjCUAPRVGeqUGbFfLIU1JSUmplX3XRlGk5GJ/Djj8yKLpYgquPPcFD3LFztLor/Usk9cHhw4fp3LlzQ5tRbWQClJpzs99xVWqZ9RGlkwU4Vzp2AnLqoZ86RW2kwifSic6hbdi74QS712Xy48ztdOzRiqDB7tjYmze0iRJJk8fZ2fn2F0numPqY0kkCOggh3IQQJsAYYFVNGqittEJtMDZV0/1+Vya8G4p/v3ak7sqVm7ckkruEhYUFFhYWDW1Gk6VWI3whxFIgArAXQmShi6//WgjxNPAXoAYWKopysIbtNnjGKzNLY0KHt8cn0pmk39M5EJfN4YQcfPs44x/lgqm5FBqVSOqaixd1+wRsbGwa2JKmicx4VU0unC5k2+o0ju04g6mlkdy8JWkUyDn8pk9N5vANUlrBEOWRr928ZXPN5i2N3LwlkUgaAQbp8BtyDv926DZv+TH0eX+sW5gSu+QoS9/aRsqO0yhaw31akkgkEoN0+LUd4f+w7TgLt6RTcKX+hNLk5i2JpHqo1Wq9/LGfnx+zZs0C4LHHHuPQoUN11k9GRgY//PDDba+LjY3F1tYWf39/OnfuzFtvvVXjvqZNm0bXrl2ZNm0a8+bN49tvdWIDEydONGjFTYNceVQUZTWwunv37lPupP6WY7n8vv8U/9uQwvjgdjwS6kpL67pXdKjYvOXibU9K0mm2r05jzdy9tO3QjJBhHrR2N7wnFInkbmNubn7TxOQLFiyo034qHP7YsWNve23Pnj1Zs2YNBQUF+Pn5MWjQIAICAvTny8rKMDKq2j1++eWX5ObmYmpqWie23y0McoRfWz4fF8AvT4YS6mHH57Gp3DdrEy//vI9jZy7VS38qlcAzqDVjZwQTPqYj508XsmL2TtZ+vk9m3pJIqiAiIoKKoAwrKyteffVVRo4cyYQJEzh9+jSgS3YSFBSEv78/ffv21Zdv3rxZ/8Tg7+/PpUuXmD59OvHx8fj5+fHRRx9RXFzMpEmT8Pb2xt/fn02bNt1gg6WlJQEBAaSmprJ48WJGjhzJ4MGDiYqKQlEUpk2bhpeXF97e3ixbtgyA6OhoCgoKCAoKYtmyZcyYMYM5c+bc0PbOnTvp1asXAQEB9O/fn5MnT9bXV1ltDHKEXxdhmd3aNeeL8QFknC1gwZY0lu/IYtmOE/Tt3JIpPd3p4daizlO/qY1UeEc40SmkDXs3nmD3uuP8OHM7nj1a02Owm9y8JWlY/pgOp/bXbZutvWHgrFteUlRUhJ+fn/74lVdeYfTo0ddcU1BQQHBwMO+++y4vvfQSX331Fa+99hr33XcfW7duRQjBggULmD17Nh988AFz5szhs88+IywsjMuXL2NmZsasWbOYM2cOa9bolNc/+OADAPbv38+RI0eIiooiOTn5mn7z8vLYunUrr7/+OklJSSQmJrJv3z5atGjBihUr2LNnD3v37uXs2bMEBgYSHh7OqlWrsLKy0j+1VGTvqkxpaSnPPPMMv/32Gw4ODixbtoxXX32VhQsX1vgrrksM0uHXdkqnMq72lswc6s2/+3bk28RMvk3MYP3hrfg6N+PxcHf6d22Nuo5VMo1N1XQf6IpXuCO7/spk36YsUnacpmtPR7rf74qFjUmd9ieRGDJVTelUxsTEhEGDBnHhwgU6deqkT3aSlZXF6NGjOXnyJCUlJbi5uQEQFhbG888/z7hx4xg+fDhOTk43tLllyxaeeUan6NKpUydcXFz0Dr9ClE2lUjF9+nS6du1KUlIS/fr1o0WLFvr6Dz30EGq1mlatWtGrVy+SkpL06RVvxdGjRzlw4AD9+vUDdFm32rRpU81vrP4wSIdfH9hZmfLvfh2Z2suDn3dlsSA+jSeX7MLFzoLH7nNjRIAz5nUcUy83b0kMituMxBsSY2NjhBCcPn2aS5cuUVamC7h45plneP7554mOjiY2NlY/mp4+fToPPPAAv//+O8HBwaxfv/6GNm8VPFExh389lZOa1Cb4QlEUunbtek2qRkPAIOfw6zMO39xEzYRgFza+EMEX47rR3MKE1387SOisDXz4dzJ5l6/UeZ8VmbfGvhmEm489O//I5LvXEti1LpOykqrzbUok9zr5+fk4OjoC8M033+jLU1NT8fb25uWXX6Z79+4cOXIEa2trLl26uk4XHh7OkiVLAEhOTub48eM12tAVHh7OsmXL0Gg05ObmEhcXR48ePapV19PTk9zcXL3DLy0t5eDBGgkO1AsGOcSsyymdqlCrBAO92zDAqzU7Ms/z5eY0PtmQwpebUxkR4MRjPd1xs6/bFGbNWlkQ9ZgX/lGX2PpbGom/pLJvYxaBD7jSObQNKrVB3n8lklpx/Rz+gAED9KGZt2PGjBmMHDkSR0dHgoODSU9PB+Djjz9m06ZNqNVqunTpwsCBA1GpVBgZGeHr68vEiRN58sknmTp1Kt7e3hgZGbF48eIaRdUMGzaMxMREfH19EUIwe/ZsWrduXa26JiYm/Pzzz/zrX/8iPz+fsrIynnvuObp2bdgkNFJaoRLHzlxmQXwav+zKplSrJapLK2LCPQhwaV4v/eWknCfx1zROpeVj29JcZt6S1DlSWqHp0+ilFRqK9i2tmPWgD1tejuTJCA8SU/N48IsERnyRwLqDp9DW8U7ath2aM3xaN+5/0kdu3pJIJPWOHOHfgoIrZSxLOsHXW9LJvlCEu4MlU3q6M8zfETPjul3g1WoVUrafYtvqdC7lFcvNW5I6obGN8EtKdBLkJiYykq261GSEb5AOvyEyXt2KMo2W3w+cYn5cKgeyL2JvZcrEUBfGB7vQzKJu/zA1ZVoObckh6XeZeUtSexqbw5fUnEbv8Cto6BH+9SiKQmJqHvPj04g9mou5sZrRgc48ep8bzi3qNmlD6RUN+zadYNdfxykpLpObtyR3RGNz+OfOnQPQx8JLbk9DpzhssgghCG1vT2h7e46eusT8uDSWbNNt5rrfuw0x4e74ODWrk76MTdUEDHCla8/rNm/d15aA+12xtG1cGh4SSXXIzc0FpMOvL6TDv0M8W1vzwShfpvX3ZNE/6fyw7Thr9p0k2L0Fj4d7EOHpUCfSDZU3b+34PZ0D8TkcTjyJb29n/KPaYWphXAefRiKR3AvIKJ1a0trWjFfu70zCK7159f7OZOYVMmlxEv0/jmP5jhOUlNVNchSr5qZElG/ecvWxZ+efmXz3WqLcvCUxeKys7mz96b333tP/nJGRgZeXV12ZdEdUlj6OiIjA09MTX19fwsLC9OGk1eXIkSN64bfU1FRCQ0OB+v+cd83hCyE6CyHmCSF+FkI8cbf6vVtYmxkzJdyduJci+XCULyohmPbzPnrO3sgXsankF5XWST/XZN5ys9Vl3no9UWbekjQ5Kjv82lIh1VCXLFmyhL179/LII48wbdq0G85rNFUPxFauXMmQIUPYvXs3Hh4eeu2g+qZWDl8IsVAIcUYIceC68gFCiKNCiGNCiOkAiqIcVhRlKjAKuGExoalgrFYxvJsTfzzbk28m96B9Syv+++cRwmZtZOaaQ+RcKKqTfnSZt3wZ9oI/1nbmMvOWxOA5efIk4eHh+Pn54eXlRXx8PABLly7F29sbLy8vvczw9OnT9Tt0x40bB+gc6JQpU+jatStRUVEUFen+L6WmpjJgwAACAgLo2bMnR44cAXQj8ueff57IyEhefvllZsyYwSOPPEJUVBSurq788ssvvPTSS3h7ezNgwABKS3WDsrfffpvAwEC8vLyIiYm57Z6Y8PBwjh07BoCrqytvv/029913H8uXL2fPnj0EBwfj4+PDsGHDOH/+PL///jsff/wxCxYsIDIyErj5U5BGo2HatGkEBgbi4+PDl19+WdtfQa3n8BcDnwLfVhQIIdTAZ0A/IAtIEkKsUhTlkBAiGpheXqdJI4SgV0cHenV04EB2PvPj0liUkMHihAwG+7YlJtydzm1sat1PxeatjP15bF2ZyroFB9nlnEnwUA/adal7CWhJ4+W/2//LkXNH6rTNTi068XKPl6t17Q8//ED//v159dVX0Wg0FBYWkpOTw8svv8zOnTtp3rw5/fr1Y//+/cyaNYtPP/1Ur7KZkZFBSkoKS5cu5auvvmLUqFGsWLGC8ePHExMTw7x58+jQoQPbtm3jySefZOPGjYBOQ2f9+vWo1WpmzJhBamoqmzZt4tChQ4SEhLBixQpmz57NsGHDWLt2LUOHDuXpp5/mjTfeAGDChAmsWbOGwYMHV/m5Vq9ejbe3t/7YzMyMLVu2AODj48PcuXPp1asXb7zxBm+99RYff/wxU6dOxcrKihdffLHKdr/++mtsbW1JSkriypUrhIWFERUVpVcMvRNq5fAVRYkTQrheV9wDOKYoShqAEOJHYAhwSFGUVcAqIcRa4Ka5yIQQMUAMQLt27WpjnsHg5WjLJw/589IATxZuyeDHpOP8ujubnh3seTzcg7D2drVyzEII3HzscfGyIyXpNNtWycxbEsMjMDCQyZMnU1paytChQ/Hz82Pjxo1ERETg4OAAwPjx40lISGDEiBE31Hdzc9Nr8gQEBJCRkcHly5dJSEhg5MiR+uuuXLkqgDhy5EjU6qubJAcOHIixsTHe3t5oNBoGDBgAgLe3NxkZGQBs2rSJ2bNnU1hYyLlz5+jatetNHf64ceMwNzfH1dWVuXPn6ssrtP7z8/O5cOECvXr1AuCRRx65xs7bsW7dOvbt26dfN8jPzyclJaXhHH4VOAInKh1nAUFCiAhgOGAK/F5VZUVR5gPzQReHXw/2NRhOzS14Y3AXnu3TgSXbM1n0Twbjv95GlzY2xIS784BPG4xrIaBWkXmrfUBLDsbnsOOPDFbM3ik3b0kAqj0Sry/Cw8OJi4tj7dq1TJgwgWnTpmFjc+1T7qVLl/RTNddTWfhMrVZTVFSEVqulWbNmVertV5Y7rtyGSqXSSzJXHJeVlVFcXMyTTz7Jjh07cHZ2ZsaMGRQXF9+07SVLltC9+42z09f3eacoisLcuXPp379/nbQH9bNoe7OhqqIoSqyiKP9SFOVxRVE+u2UD9SiPbAjYWhjzZER7trwcyewHfSjRaHlu2R4i3o9lQXwal2uZfF1tpMIn0onxbwcTFO1OTvJ5fpy5nfWLDnHxbN2sIUgkNSUzM5OWLVsyZcoUHn30UXbt2kVQUBCbN2/m7NmzaDQafvrpJ/30iLGxsX5evSpsbGxwc3Nj+fLlgM5J7t27945trHDu9vb2XL58uVYJyW1tbWnevLl+reK7777Tj/arQ//+/fniiy/030FycjIFBQV3bA/Uzwg/C3CudOwE5NRDP40eUyM1owKdGRHgxKajZ/gyLo2Zaw/zyYYUxgW7MCnUlZY2d5583cTMiO73V8q8FSszb0kajtjYWN5//32MjY2xsrLi22+/pU2bNvznP/8hMjISRVEIDg6mT58+AMTExODj40O3bt149913q2x3yZIlPPHEE8ycOZPS0lLGjBmDr6/vHdnYrFkzpkyZgre3N66urgQGBt5ROxV88803TJ06lcLCQtzd3Vm0aFG16z722GNkZGTQrVs3FEXBwcGBlStX1sqeWksrlM/hr1EUxav82AhIBvoA2UASMFZRlBqr/xuatMLdYPfx83wVn8afB05hpFIx1F+3wNu+pXWt2758/gpJa9M5nHAStbEK395OMvNWE6exSStIeeSac9ekFYQQS4EIwF4IkQW8qSjK10KIp4G/ADWwsKbOvi6SmDdW/Ns15/NxAWTmFbAgPp3lO0/w044s+nRqyZRwd4JqkXzdqrkpkeM74d+vHdtWpbHzj0wOxGUT0N8V7whHjOo4xaNEIjEsDFI8rdZqmVoNqJqG8zpXUMJ35cnX8wpK8HWyJSbcgwFetU++nnv8EltXpnL80Dksm5nKzFtNEDnCb/pItcxl4yF7Fzh4gr2n7t2hk+7donGKMhWXalixK4sF8emkny2gXQsLHuvpxsg6SL6effQ8iStTOZ1+kWatLAiKdsfD30Fm3moCNDaHX7E7tXIopeTWSIe/8xvI/Adyj8LZZCgtvHrOsqXO8bfsrHs5dIaWncC8ftIY1jUarcLfh04zPy6VXccv0NzCmAkhrjwS4oKd1Z0raCqKQvres2xblca5nAIc2lkTPMQdZ7l5q1HT2By+pOY0eodfpwlQtFq4mKVz/rlH4MwRyD2sOy65fPU667bQqgu07AKtuureHTzByHBliHdknGPe5jTWHz6NqZGKBwOcmFLL5OtarULy9lNsL8+85dixGcFD5eatxkpjc/hnzpwBoGXLlg1sSeOh0Tv8Cuo1SkdR4MLx8pvAITh9CM4chrNHQaNLs4ZQg30H3Q2glZfu1dobrFuDAY16U3N1yddX7MqmVFORfN2dAJc7n77SlGo5uCWbHb9nUHSpFDdfe4KGuGPXVm7eakw0Nocv5/BrTqN3+A2a4lBTBudS4fTB8tcB3c0g//jVayzsdI6/tTe09tG923UAdcOGN+ZeusI3CRl8tzWT/KJSAlyaExPuTr/OrVDd4Xx8SXEZ+zaeYNe645Re0dApqDWBg2TmrcaCITh8tVqNt7c3iqKgVqv59NNP9XLA11OVw4+IiGDOnDk33dlaHRYvXkxUVBRt27bVt3fy5EnMzMywsrJi4cKFNbrJHDlyhDFjxiCE4Oeff2bChAkkJCSQkZHBoEGDOHDgwO0bqSMafcYrRVFWA6u7d+8+5a53rjYqX+T1BK/hV8uLLuieBE4dgFP74NR+2DYfNOW6HUZmuieB1j7Qxhfa+ummhe7ilJCDtSkv9vfkiQgPftpxggXx6Tz+3U7cHSx57D53hnerefJ13eYtN7qGO7Lrz0z2x2aTnHQar3BHAgbKzVuS22Nubq6XPvjrr7945ZVX2Lx5c7Xr30pmuLr1Fy9ejJeXl97hw1VphPnz5zNt2jRWrVp1Q72qFo8r5I3feustgLsmb1xbDNLhGyTmzcAlVPeqQFMGeSlwcp/uJnByLxz4BXaW76ZTGesWhtv6QRs/aOuvmxYyql8naWlqxKQwNyYEu+iTr//fr/v58O+jPBLiyoSQmidfN7cyIWxEB3x6O7NjbTr7Y7M4lHASv77O+Pdth4ncvCWpBhcvXqR5c12ARGxsLHPmzGHNmjUAPP300zg6OjJ8+HBcXV2ZPHky69at4+mnn9bX12q1TJo0CWdnZ2bOnMm6det48803uXLlCh4eHixatAgrK6tr6k+dOpUdO3boxc4SExOvsSk8PJyPP/4Y4IZ+O3XqpN8p6+HhwcKFC0lMTOTjjz9GrVYTFxfHpk2bsLKy4vLly9e0q9FomD59OrGxsVy5coWnnnqKxx9/vD6/3tsi/5fWBrXR1WgfX51CHooC5zN0zv/kHsjZA4dXw65yBWm1ie5JoG03cOyme3fwrJd9A0ZqFdG+bRns04bE1Dy+jEvjg7+T+Tw29Y6Tr1u3MCNyQmf8+rVj26p0dqzN4EBsNgEDXfDq5YhRDZ8gJHePU++9x5XDdSuPbNq5E63/7/9ueU2Frn1xcTEnT57USxffjsoyw/PmzaOsrIxx48bh5eXFq6++ytmzZ5k5cybr16/H0tKS//73v3z44Yd6aePK9RcsWFDllFBjkjeuLQbp8Bv1TlshoIWb7tV1qK5MUeBCJuTs1r2yd8H+5bDja915Y0vd6N8pABwDwCkQbNpW3UeNTbqafP3IqYvMj0vj+61Xk68/Hu6Bt1PNonCat7ZkQIwXZzIvsnVlKv/8fIy9G04QOMiNTsGt5eYtiZ7KUzqJiYk8/PDDVc5xt2nTRj+XXiEzXMHjjz/OqFGjePXVVwHYunUrhw4dIiwsDICSkhJCQkL0119f/3oao7xxbTFIh9+gc/j1gRDQ3FX36jpMV6bVQt4xyN4JObsgawckfg7acnVA67bg1F3n/J176NYFjGu/UNqptQ0fjvJjWn9PFv+ToU++HuJuR0wvdyI61iz5eksXG6Kf9SfryDkSV6ax6bsj7Pn7OEHR7rj7100id0ndcLuR+N0gJCSEs2fPkpubi5GREVrt1bSc18sQXy8zHBoayqZNm3jhhRcwMzNDURT69evH0qVLb9rX7WSKG6O8cW0xSId/T6BSgUNH3cvvIV1Z2RXdYnDWDshK0r0Oly8kqYx1Tt+5R/krGGza3HH3bWzNeeX+zjzVuz0/bj/Owi0ZTFqUhGcra6aEuxPt2xYTo+qP0p06tWDEy81J33OWrb+l8uf8A7R0sSZ4qAfOnRvn7mZJ3XPkyBE0Gg12dna4uLhw6NAhrly5QnFxMRs2bMDb25tTp07dtO6jjz5KXFwcI0eO5NdffyU4OJinnnqKY8eO0b59ewoLC8nKyqJjx4431LW2tubSpUs1srWyvHHPnj3vWN64d+/eGBsbk5ycjKOjY53dUO4E6fANCSPT8lF9d2CqruzyGZ3jP7ENTiTBjoWw9XPduWbtwDkI2oXoXg6ddDeSGmBjZkxMuAcTQ91YvTeH+XFpvLh8L3P+OsqkMFceCmqHjZlxtdoSQuDu74Crrz1Ht55i++o0Vv1vD06dmhM81INWrrVP6ShpfFTM4YNu1PvNN9+gVqtxdnZm1GWKGaIAACAASURBVKhR+Pj40KFDB/z9/SksLORWeTCef/558vPzmTBhAkuWLGHx4sU89NBD+ixXM2fOvKnDnzhxIlOnTr3pou2tMDR549pikHH4FdyL8si3paxE9xRwYqvuJnB8K1w+rTtnZnvV+buE6aKD1NVz1hUoisLm5Fzmx6WRkJqHlakRY4PaMSnMlTa2NZtSKivVcGBzNjv/zKT4cike/g4EDXGneeuGG+HcaxhCHH5NkBuvao7ceHUvUREVdDwRMhN0N4C88u/M2EK3BuASBq736Z4carAvYH9WPvPj01i7LweVEET7teXxcA88W9dMm7+kqIw9G06w5+/jlJVo6BTahsAH3LBucefJXSTVQzr8pk+jd/gVyBH+HXL5zNUbQMY/ut3CKLrNYU6B4NoT3HqCY/dq7Qk4ca6Qr7eksyzpBEWlGiI8HYgJdyfEvWbJ14sulbDzj0z2x2UhEHhHOBIwwBUzq5o9hUiqj3T4TR/p8CXXUnQeMhMhYwtkxOumhFDAyBzaBYNbOLj30m0Ou8V+gAuFJXy/NZPFCRmcvVyCt6MtMeHuDPRqjVENwjAv5hWRtDqdo9tOYWyqxq9fO3z7OGNiJpeU6prG5vCTk5MBbjoPL7k50uFLbk3hOd3oPz1OdwM4c0hXbmqrG/m79QL3CJ1w3E1G8MWlGn7Zlc2C+DTSzhbg3MKcR8PcGBXojIVJ9Z12Xs5ltv2WRvres5hbG9P9fle63ueI2ljG8NcVjc3hS2qOdPiSmnH5jM75p2+GtFidiiiAjZPO8XtEgnskWNpdU02rVfj78Gnmx6WxM/M8zSyMmRDswiOhrtjXQJv/VFo+W1emkp18AWs7M4IGu9GhR+s7FnyTXEU6/KaPQTp8IcRQ4AGgJfCZoijrbldHOvwG4lyazvGnbtLdBIrzAaGL+vHoA+376NYCKkUA7cw8x5eb0/j78GmM1SpG1FCbX1EUThw+R+KvqZw9cZkWbS0JHuKOq4+93LxVCxqbw8/JyQG4RuRMcmvumsMXQiwEBgFnFEXxqlQ+APgfuiTmCxRFmVXpXHNgjqIoj96ufenwDQCtRicHkboRjm3Q7QlQNGBqo5v3b98POvTTS0HcXJvfgwCX6mUUU7QKx3adYdtvaeTnFtHa3ZaQYR607dCsPj9lk8UQHL4QgvHjx/Pdd98BUFZWRps2bQgKCtILp1VQl4u2GRkZdO7cGU9PT0pKSggPD+fzzz9HVYO9Kp988glffPEF3bp1Y/To0Rw6dIjp06czY8aM22rp3C3upjzyYuBT4NtKHamBz4B+QBaQJIRYpShK+UQxr5WflzQGVOqrm8F6vaSTiU6Pg2N/624Ah1frrmvlBR364dGhP/8ZEsjz/Tz12vx/HTxN93Jt/r630eYXKkGH7q1w93fg8D8nSVqbzq8f7MLFy47goe7YO9UsJFTS8FhaWnLgwAGKioowNzfn77//xtHRsU7aLisrw8ioajfm4eHBnj17KCsro3fv3qxcuZLhw6/Knt9KAhng888/548//tDr30RHR9eJ3Q1FrVbHFEWJA85dV9wDOKYoSpqiKCXAj8AQoeO/wB+Kouyqqk0hRIwQYocQYkdubm5tzJPUB+bNoEs0RM+Ffx+EJxKh39u6nMAJc2HRAHjfA4d1T/Fi2/0kPNuNNwd34WR+MTHf7aTvR5v5cftxiktvrXGuVqvwCndk/DshhAzz4FRaPsveTWLd1wfJzy28ZV2J4TFw4EDWrl0LwNKlS3nooYf0586dO8fQoUPx8fFh9OjR+lF+5fLg4GD27dsHwIwZM4iJiSEqKoqHH34YjUbDtGnTCAwMxMfHhy+//PKG/o2MjAgNDeXYsWPExsYSGRnJ2LFj9SqZH374IV5eXnh5eemlkqdOnUpaWhrR0dF89NFHLF68+Bqp5gpSU1MZMGAAAQEB9OzZkyNH6laRtC6pjzg4R+BEpeMsIAh4BugL2Aoh2iuKMu9mlRVFmS+EOAkMNjExCagH+yR1hRC6PMCtukDYs7q5/tSNkPwXpKyD/cuxFGomuYTycM/+bFIC+GiXhum/7OeDv5OZGOrK+CAXbC2qjsM3NlHTrb8LXe5ry+51x9m38QSpO8/QpWdbut/viqWt4eYcNjTif0rm7InLt7+wBtg7W9Fz1O1DKMeMGcPbb7/NoEGD2LdvH5MnTyY+Ph6AN998E39/f1auXMk333zDyy+/THR09DXlGzdu5OGHH9arbu7cuZMtW7Zgbm7O/PnzbypDXHntp7CwkA0bNvD2228DsH37dg4cOICbmxs7d+5k0aJFbNu2DUVRCAoKolevXsybN48///yTTZs2YW9vz+LFi2/62WJiYpg3bx4dOnRg27ZtPPnkk9WWgL7b1IfDv9nzuqIoyifAJ9VpoMmpZd4rmNnq1EC7DtPN/WfvgqO/Q/KfqP9+jb5AH4dOZAVE8nVuJ+b8VcTnm44xOrAdj/Z0w7FZ1dINZpbGhAzzwCfSiaTfMzgYn8ORxJP49nHGP8oFU5mAxaDx8fEhIyODpUuXcv/9919zbsuWLaxYsQKAsLAwLly4QH5+/jXlvXv3Ji8vT6+zEx0djbm57u+lKhnijh07kpqaip+fH0IIhgwZwsCBA4mNjaVHjx76aZotW7YwbNgwvajZ8OHDiY+Px9/f/7af6/LlyyQkJFwjm1yh62OI1Mf/kizAudKxE5BTkwYatR6+RIdKDc6BulffN3XyD0f/RBxdi/Oh+cxQNLzaohXbTIKZv7ULvRO7MNDHmZhwD7q0rVpkzbKZKRFjPfHr48z21Wns/COTA3HZBPR3xTvCESMTmYClKqozEq9PoqOjefHFF4mNjSUvL09fXjlwpH379hgZGSGE4GYBJRWj9sqKk1XJEGdkZOjn8K/n+vp3ilarpVmzZjftwxCpjx0uSUAHIYSbEMIEGAOsuk0dSVOnuSsET4VHVsO0YzBsPsYuPbiv4G++Nf4Pe8yfJOLwm3z46cdMXhDPP8fO3vI/YrNWFkQ95sWo/wuklYsNCb8cY8mbWzn0Tw5ajbbKepKGY/LkybzxxhvXZJcCXYrBJUuWALq0h/b29tjY2FRZfj0VMsSlpbpcEsnJyRQUFFTbrvDwcFauXElhYSEFBQX8+uuv9OzZs1p1bWxscHNzY/ny5YDu5rF3795q9323qdUIXwixFIgA7IUQWcCbiqJ8LYR4GvgLXVjmQkVRDtakXTml08SxaKFLCek7GkoKIXUj5odXM+To7wwVmynI+owN3/jznm0EvhEjGODvXqV0g0M7awb/y4+so+fZujL1agKWIe64+8kELIaEk5MTzz777A3lM2bMYNKkSfj4+GBkZMTs2bNvKLewsOCbb765abu1lSHu1q0bEydOpEePHvr2qjOdU8GSJUt44oknmDlzJqWlpYwZMwZfX99q17+bGOROW6mWeY9SVgIZcZQd+I2yQ6sxKzlPgWJKolEgqq7DCO4/GgvLqsMyFUXRJ2A5f6qQlq42hAx1x6nTvZuAxRDi8GuCFE+rOQa50/ZOkBuv7mE0ZWjTt5Cd8AM26X9gq1ykADOO20fQtud4bLv2r1LpU6vRcmTrKZLWpHP5/BWcu7QgZKgHDu3uvRh+6fCbPndz41W9IBdtJaiNULWPwLl9BGjKSN7+O6cTluKduxnbX/+kcJU1Gs/BWHcfo9P6r6TyqVKr6BLWlo49WrE/Npudf2bw03tJtO/ekqBod5q1tGi4zyWRNCByhC9pVKSeOkfcHz/RIn01fcQOrEQxpRYtMfYZCT6jdHl/r5u3v1JUxu51mezdcAJtmULn+9oS+MC9EcMvR/hNn0Y/wpdIqsKjdQs8Jk3lzMWJzI8/Qk7Sr0Rdiqf31i8x2voZin1HhM9onfNv1g4AU3Mjgod44B3hxI7fMzgUn8PRrSfx7e2Mf38Zw29IGBvLZDj1iUGO8OWiraS6XL5Sxo/bj7M8fh/dCuIYa5aIt6ZctsklDHwfgi5DwOxqON+FM4VsX51OStJpTC2NdDH8kY4YGTe9GP7GNsKX1By5aCu55yjVaFm1J4f5cWkUnEllgsU2xpj8g23RcV1mr86DwW+sLrlLuVpi7vFLbF2ZyvFD57BqbkrgIDc6BbdGVYPsXYaOdPhNH+nwJfcsiqIQm5zLvNhUtqXncZ9ZOi+22oXPhQ2oruTrkrr4PQR+46CFbmt91tHzJP6aypmMizRvbUHwUA/cfJuGDr8hOHwrKysuX66ehs/x47rkO+3atavymupKE9ek3+qwatUqvTzyypUr6dixI126dKmz9u8UOYcvuWcRQhDp2ZJIz5bsOXGB+XFtGHbAHUv1UP7PPZ2hbMIibg7Eva9L5u4/Aacu0Yx4OYC03bls/S2NP+btp7W7TbkOf/V0/CV1Q1FRUUObUCXR0dF6eeSVK1cyaNAgg3D4NcEgn12FEIOFEPMrhJIkkjvBz7kZn48LYOMLEUQHuDMjzZOuKTG84rKUkwEvQv4J+DUGPvBE/PESHm3P8NAbPYgc34lLecX8+sFu1ny2l7zsulWYvFdRFIVp06bh5eWFt7c3y5Yt05+bPXs23t7eDBkyhA8++ACAr776isDAQHx9fXnwwQcpLLy1LHZ6ejohISEEBgby+uuvX3Pu/fff18snv/nmm8DVBClTpkyha9euREVF6W84n3zyCV26dMHHx4cxY8YA6OWRExISWLVqFdOmTcPPz4/U1FS6deum7yslJYWAAMMU+jXIEb6UVpDUJW72lrw3zJt/9+3I4oR0vkvMZOnRboS49ealyLP45a5C7FwM2+ejcuxOl4CJdHh9CPv/Oc+uvzL5ceZ2PHu0psdgN2zsq1b0bAxERETcUDZq1CiefPJJCgsLb1CyBJg4cSITJ07k7NmzjBgx4ppzsbGx1e77l19+Yc+ePezdu5ezZ88SGBhIeHg4e/bsYeXKlWzbto0TJ05w4cIFQKdaOWWKzgW89tprfP311zzzzDNVtv/ss8/yxBNP8PDDD/PZZ1dzLK1bt46UlBS2b9+OoihER0cTFxdHu3btSElJYenSpXz11VeMGjWKFStWMH78eGbNmkV6ejqmpqZ6eyoIDQ0lOjqaQYMG6b8PW1tb9uzZg5+fH4sWLWLixInV/l7uJgY5wpdI6gMHa1Om9e9Ewit9eO2BzmScK2bYH0YMzHqEtVGxaKLegyuXYNXTGM/tSreyuYz/lx3+fdtxbOcZlszYypafUii6XNLQH6VRsmXLFh566CHUajWtWrWiV69eJCUlsX79eiZNmoSFhW5DXLNmunSWBw4coGfPnnh7e7NkyRIOHry1JNc///yjT6wyYcIEffm6detYt24d/v7+dOvWjSNHjlAR/efm5oafnx8AAQEBZGRkADo553HjxvH999/fMqNWBY899hiLFi1Co9GwbNkyxo4dW7Mv5y5hkCN8iaQ+sTI14rGe7jwc4sqqvTnMj0vlqV8zcGzWmUfDfmRs2xzM9n4Lu77BLOkrQtuF4j12CkkpnuzbdIJDCTl0i3LBt48zxqaNK5TzViNyCwuLW563t7ev0Yj+eqoKEFEURb9Abmp6dTPcxIkTWblyJb6+vixevLhafd9soV1RFF555RUef/zxa8ozMjKu6U+tVuundNauXUtcXByrVq3inXfeue3N5sEHH+Stt96id+/eBAQEYGdnd1tbGwKDHOHLOXzJ3cDESMWIACf+fDacrx/pTttmZry99jDBSwr50PpFzk/dB/3egUsnsV43id5nhjAmahdO7uZsW5XG968nciAuG42UY64W4eHhLFu2DI1GQ25uLnFxcfTo0YOoqCgWLlxIYWEhrq6uegnkS5cu0aZNG0pLS/UyybciLCyMH3/8EeCa6/v378/ChQv1ETvZ2dmcOXOmyna0Wi0nTpwgMjKS2bNnc+HChRuifaytrbl06ZL+2MzMjP79+/PEE08wadKk6n8pdxmDdPiKoqxWFCXG1ta2oU2R3AOoVII+nVuxfGooK54IIdC1BZ9sSCHkkz3MyOtD1oQtMG4FOHanxd6Z3H9+AMMD1mJrU8rmH47y49vbObbzTK0SadwLDBs2DB8fH3x9fenduzezZ8+mdevWDBgwgOjoaLp3746fnx9z5swB4J133iEoKIh+/frRqVOn27b/v//9j88++4zAwEAqDxajoqIYO3YsISEheHt7M2LEiGuc9fVoNBrGjx+Pt7c3/v7+/Pvf/9ZPM1UwZswY3n//ffz9/UlNTQVg3LhxCCGIioq6k6/nriDj8CWSm3DszCW+3JzGyj3ZaBUY7NOGqREedDI9D0kLYNe3KEUXyDAfQWL+KM6fN6alqw2hwz1w7Gg4oZyGEIdfEyrm0F1dXRvUjjthzpw55Ofn884779zVfmUcvkRSS9q3tOb9kb48H9WRr+PT+WH7cVbuySHS04EnIv5NYMR0xP6fcds2D5eiMRx1GMT202NZ+eFFXLztCBnqgZ2jVUN/jEaHIeeDvRXDhg0jNTXVYJOXVyBH+BJJNbhQWMJ3iZksTsggr6CEAJfmTO3lQR9PB1THt8DWeZQdWc++ogfYVTiaK2UmdAppTY/B7li3MGswuxvbCF+qZdYcgxzhCyHcgVcBW0VRRtzueonEkGhmYcIzfTrwWE93lu88wfy4NKZ8u4MOLa2Y2suD6FHfY3whnW7b5tFl59PsujCQfVsHk7L9JN6R7QgY6IqZpVSClDQstVq0FUIsFEKcEUIcuK58gBDiqBDimBBiOoCiKGmKojxam/4kkobG3ETNwyGuxL4Ywcej/VCrBC8s30vE+7EsPqKiqO8szF5MInSYO+Nc3qSD8Sb2rM/k+1c2s+vPNMpKNXfdZkN+ipfUjpr+bms1pSOECAcuA98qiuJVXqYGkoF+QBaQBDykKMqh8vM/V3eEL6d0JIaOoihsOnqGzzelsiPzPC0sTZgU6srDIa7YmihwYAVnNywjMTOM4yUBWFmUEDTMk45hrqhU9S/Olp6ejrW1NXZ2do1CDK464mkSHYqikJeXx6VLl3Bzc7vmXL2pZQohXIE1lRx+CDBDUZT+5cevlBv3n/Jj6fAlTZKkjHN8EZvKxiNnsDRRMz7YhUfvc6OltSkcW0/W77+QkOJPbll77GwLCB3lhXM313p1xKWlpWRlZVFcXFxvfUgaDjMzM5ycnG5IHHM3Hf4IYICiKI+VH08AgoA3gXfRjfwXVNwAbtJeDBAD0K5du4DMzMxa2SeR3G0O5Vxk3uZU1uzLwUit29w1NdyDdnYWKCd2cOzX39h6tDMXNa1xcsgj9CF/HLrI/M2SuuNuOvyRQP/rHH4PRVGqVj26sU2Z8UrS6MnMK+DLuDR+3pFFmVbLYN+2PBHhQafWNmhOHuXAsrXsSHanWGtDx7bHCRoXio1Hx4Y2u0EZP348AN9//30DW9K4MZgpnZogp3QkTYHTF4v5eks6S7ZmUlCioW/nVjwV6YF/u+ZcOZXOru/Xs/eYMwoCH5dUAsb3wcz53gxLrFDzrI1mj+TuOnwjdIu2fYBsdIu2YxVFubX60LVtyhG+pMlxobCEbxIyWZSQzoXCUkI97Hg6sj0hHnZczspi+3ebOXK8NaaigO7tk/EeOxh1m9tLCjQlpMOvG6py+LUNy1wKJAKeQogsIcSjiqKUAU8DfwGHgZ9q4uwlkqZKMwsTnu3bgX9e7s3/3d+JlDOXGbtgG8O/SGDbJRN6vzKO0c93pKVDKf+kBPDDzF2kfPoWypkjDW26pIkgd9pKJA1EcamGn3dmMW9zKlnni+jcxoanIj0Y6NWG7D2ZJCzbR16+FS2NkwnzPU7boZPBvmkv7soRft3Q4Dtta0KlKZ2GNkUiqTfMjHWhm6MDnVm1J4fPY4/x9A+7cbdP5snI9gyb+QDp/6Sy9bcr/LqjI277lxPa4wLN7n8Smrs0tPn1QkhISEOb0KSRI3yJxEDQaBX+PHCKTzcd4/DJizi3MOeJXu0Z4t2aw+uS2bU+B00ZeFmuo3uYEeb9ngPr1g1ttsQAqbdF2/pALtpK7mUURWHD4TPM3ZjC3qx82tiaMbWXB4M9W7Fv1SEObbuAsSiku81KfCKdUfd6FswNR5JZ0vA0KodfgRzhS+5lFEUhPuUsczemkJRxHgdrUx4Pd+f+dvbs+uUgmUeLsFGfIqT5CjyiwhDBU8HEoqHNrhUPPvggACtWrGhgSxo3jWoOXyKR6PKzhnd0oGcHe7amneOTDSnMXHuYeVYmTOnpTlRvG3auUPPXmadovfwwYfHDaT3gIfAbB+rG+V87Ly+voU1o0hjkX4VctJVIriKEIMTDjhAPO5IydI7/P38c4UtLEx6LcCNUZcru1bDiRGc6fBdHSNz9WN//PHTsD41AME1y95BTOhJJI2Rn5nn+tyGFuORcmlsYMyXEDe+LcHDjcdCW4Wexkm6dT2Iy8A1w7NbQ5lYbGZZZN8gpHYmkCRHg0pxvJ/dg93Gd45+9IVnn+O93wTWrmJ27R3J49wWCk9+hU5ADou8bYOvU0GZLGhiDdPhySkciqR7+7ZqzeJLO8X+8PoXZ8cdoYWlCTJQjlkcs2Hj8GfbHpnHfnrG0jYyCsOfA1HBz7fbp06ehTWjSyCkdiaQJsev4eT76O5n4lLPYW5oQ49oak/3nKLioob3ZFkJa/4nNgGfAZzSoaqWsIjFgZFimRHIPsSPjHB/+nUxCah5trUyZbGdH2cHzoCnF3+JXurVPw3jQu+Ac2NCmSuoB6fAlknuQhNSzfLgumR2Z5+lgZc5oEyuupF3GyugcIZaL6NDDCRH1lsHs2B04cCAAf/zxRwNb0ripF7VMiURi2IR62LN8agiLJwViZmvCzHO5xDoJSmxb8Xf+C/y62ZfcD4ZDwlzQlDa0uRQVFVFUVNTQZjRZDNLhCyEGCyHm5+fnN7QpEkmjRwhBhGdLVj0dxrzxAeRbqniv7CL726jJVXXmp1PvELs8i6K5/SAttqHNldQjckpHIrnH0GgVftuTzUfrkzlztojhplY4ny3DRBQSZPU9XbtbohrwLti0ueu2yTj8ukHG4UskEgDUKsHwbk4M8mnLj0nH+WTDMRSrEkYJa+IuPs6h+HTCDz5EmwEjocfjjVamQXIjTfI3WVxWjFbRNrQZkgZCVCEnIBA3vUaU/ys/0P+sfxflVzQxmQITIxUPh7jyYDcnFm5J58vNaThptDxQ6sovp9/Ac9kmQpIGYTl8JjjdMFisFwYNGnRX+rlXaZJTOlPXT+Wf7H/qwSKJRHcjUAmV7oYgrh5XlKmECiF072qh1p1DhUqlOzZSGaEWatQqNUbCCGOVMUaq8ne17t1UbYqp2hQTtQlmajNM1aaYGZlhYWSBhbEFFkYWWJlYYWVshZWJFbYmttia2mJmZHbHn+tcQQmfbjzG0sQMgouM6F6sxkQUE2S5BK+ebVBFvQFmtnX3RUrqjQYPyxRCWAKfAyVArKIoS25X504d/t+Zf5N9KbvmRkoMCoWa/21WVafy3/n111ScU1Cu+Vn/rpSfKz9/q3etokWraFFQ0Gg1aBUtGkWjP9Yo5S+thjJtGWVKme69/FWqLaVEU0KJtoQSTQnFZcVc0VyhWFNMmbbstp/fTG1Gc7Pm2JnZYWduh4OFA60sWtHKohVtrdriZO1EK4tWGKmqfrg/ca6QOeuOErfzJINKjGh9xQg7o3QiWi2n9YNPQOdoKcpm4NSLwxdCLAQGAWcURfGqVD4A+B+gBhYoijJLCDEBuKAoymohxDJFUUbfrn25aCuRXKVUW0pRWRFFpUUUlBVQUFLApdJLXC65TH5JPvlXdK9zxefIK8ojrziPM4VnOFd87pp2jIQRTtZOeDTzwN3WnY7NO9LZrjPO1s6oxNXAvX1ZF3hv7SHOHc1n4BU1JhojupivI8TvJGZD3wNbxzr/jHLRtm6or0XbxcCnwLeVOlIDnwH9gCwgSQixCnAC9pdfpqllvxLJPYexyhhjE2NsTGxqVK9EU8LpwtPkXM4h61IWJy6dIPNiJqn5qcSeiEWj6P47Whpb0tWuK34t/ejWshu+LX1ZGhPChsNnmL32EI5ZV6CoL2nbCgg9+gKdhvZFdJ8sJRoaEbVy+IqixAkhXK8r7gEcUxQlDUAI8SMwBJ3zdwL2YKDx/xJJU8REbYKztTPO1s4EtQm65lypppTU/FQO5x3mYN5B9p/dz4L9C9AqWtRCja+DL6FtQ3n/kTB2p1jyzZ8pRJwzZuO5GI58d5BeSRNoMfotsJdCh42BWs/hlzv8NRVTOkKIEcAARVEeKz+eAAQBL6N7GigGtlQ1hy+EiAFiANq1axeQmZlZK/skEknNKCwtZG/uXraf2k5CTgKH8g4B0MayDT3bRnL2lCeZW8zpU6TCSBH4W6+i++COGIU9WesQTjmlUzfczTj8m63mKIqiFACTbldZUZT5QoiTwGATE5OAOrdOIpHcEgtjC0LahhDSNoRnuz3LueJzxGXFsSFzAytTf6ZEW0LrEEfWnQ+h/UEftJeGc2xZDhFJU3Ca8BK07NzQH0FSBfXh8LMA50rHTkBOPfRTJUppKTRkuOndiGCoTh+1ueYWdZtaPLrk1rQwa8HQ9kMZ2n4ol0sus/HERlYdW8W2wp/J7LaClNx+RKSH81vKBDrNWkLYAy0xi7iz0f6oUaPq4RNIKqiPKR0jIBnoA2QDScBYRVEO1rTtO43SOT4lhoL4+BrXk9SSyjeCGvwsKpcLcdOfReXj66656TlRvlGq4qVS6WLmhar850rnVSpQiavnVLo4etTqq+Vqtb5MqFWgUoNahVAb6Y+FWg1Gal2Zkbr8WiOEkZGu3MgYYWSEMDZGGFe8V3qZmJS/TBGmJqhMTRGmpggTU1Rmpghzc1RmZroyA1gozbmcwy8pv7A8eTn5BZfofnwofqdCsFBd4j6XDXSY9C9ES8+GNvOepL7CMpcCEYA9cBp4U1GUr4UQ9wMfOTdI0gAAFcRJREFUowvLXKgoyrs1bLci49WUlJSUGtuVv3YtpSeyalyvbrgLTxbV+J1V6/da1TW3qlq5zjX1K8W5V3XNNZcrld6v/qyvq1x3Tfl1iqJcbafivFZb6Zxyta6ivVqmVfTtKFptpZ81V6/XanTXaxXQaFAUre5nrQalokyrBa1WV09T/l6mQdFodOc1GpSysmt/Li29prwuEGZmqMzNUVlYoLIwR2Vhicqy/GVlhcraCrWVNSoba9Q2tqhtbVA3a6Z7NW+O2tZWd4OqA0o0JfyV8RcL93/D+ZzLRB4bi32hE85mO+gzuBmWkY9VO5KnsLAQAAsLizqx7V6lwTde3QkyDl/S1FC02qs3gdJSlLIy3XtJCUpJCdorJfqflZIraIuLUa6UoFwpRltUrH/XFhWhLSpEKSxEW1iItqAATUEB2oICtJcL0F66hLagoGpDhNA5/hbNMbL7//buPK6qOn3g+OdhFxDDEEmWNEUKxwYIK50Wy3LUac+ZsjK3XMvKFseln9mMbWOLS+ZSNpbTovWbTM2yrJymtCZTyl2QFMgQklBAZP3+/hD7MchF5C7nLs/79bqvF+fce855Ho734XjOd4kiIDqagHbtCGwfTUDMWQTGtCewQwf8o6KafQvPGMOXB77kua8WELqnLT1yBhAoVfSM/5yU0XcjbTuech/60NYxPKrg23uFr5QCU1NDbWkpNUeOUHP4CDWHi6n5pZia4mJqioqoLjpEzaFDVBf+THVhIdWFhZjKyv/ah4SEEBgXS1BcPEEdOx5/depEcNdEAiIjbR77i7yN/G3tIrrv6kmHki6EttrGTTdG0+bSW5p8PqQF3zE8quCfoFf4SrmOMYaa4mKq8/Op+imfqgMHqMrLo+rHPCr351C5fz+mouLXz/tHRRGc2IWQ5GRadetGSLduBCYk/Po/AmMMb237gA/e/5S0nD74YUjquIG+d09Aws5sNAYt+I7hUQVfr/CVcj+mtpbq/Hwq9mZTkZVFRWYmFbt2UZGZebxlHODfti2tUlMJTU0h9KKLCElOpsrUMH3Vy8hGPzoc6UpN2E4G3nIOHS7sf9IxtOA7hkcV/BP0Cl8p92cqK6nIyqJ86zbKt2yhfMsWKus6TPpFRBB6YQ/CL7+covPPY9aKFXTefwFILe27fMut90xEAv9/hE8t+I6hBV8p5TLVhYWUff0fyr7aSNmGDVQf+AmAkO7d2dY1ie/KEog81pnSiK3cMvQiOiZfDMCSJUsAGDp0qEWRewePKvh6S0cp72GMoWLPHko/+4ySTz7l2NatGITt3fqR3+73VAUcI+r8TO4c+aAOu+wgHlXwT9ArfKW8T2VODkfWfMDh1aso+uko288bSll4HMfC/kP/YQOIjulAVFSU1WF6NC34Sim3YoyhPCODg28u4/sdkBd7FbNX3ktwaAXfZ2a6RW9iT2Wr4Lvlb1RErhWRRYcPH7Y6FKWUk4gIoampdPrbU1zzxhQSWn+CGKGiPJx/DpzIzx9+1Lwe46rZ3LLgG2NWGWNGtWmj82cq5Qv8IyK49pkniTw7DCOl5EcPYO3fs8m4YRCl//7C6vC8hlsWfKWUbwppFUpslzhqun5JUeTZfN3uNjZPmkPOqNFUZGdbHZ7H04KvlHIvItz7wP/Q/Y9lFLUuYmv3MWzKjybr+ps4+ORTTY8RpJrkjPHwlVKqRcaOHfvrz737DOT88/czd+5rwBUcbNeZC5cvpuTjj4mZ/ijhl11mXaAeyi1b6Wg7fKXUCaamhnkLplG5oycBtXBuzkri9/6biGuvJeaRqfjrs76TaLNMpZTby83NBSA+Pv6k9z75ZDEb3ven7dEEgqo30evr1wk6sy0dnnqKsIsvOunzvsyjmmUqpXzT4MGDGTx4cKPv9ekzgmETupMX/QWVAel8dPkESoLPIGfYMA7OnPnrAG7KNi34SimPERd/AdOn3Edh4rtIbXv+1WUI+y4dSNHiV9g/ZChVBw9aHaJb04KvlPIowSFteHTCbMIv/ZyjgcfI8ruUz/vdQ/nOnfxw082UffWV1SG6LZcVfBE5R0QWi8g7rjqmUso7iZ8fdw16kktuPcqPkdupPnYeay6/h+qItuSMuIuipf/QXrqNaFbBF5FXRKRARLY1WN9PRHaLSJaITGpqH8aYbGPMCHuCVUqp+i6/eDhj7unFntg1hBzryIeJgyjpcSUHH3+c/OmP6X39BprbDn8J8ALw2okVIuIPzAOuBvKAb0RkJeAPPNlg++HGmAK7o1VKebUHH3zwtLfpGHshf3mgI4++NJmzM29mY9hVdOkbyznLllC5bx9xc+fgHxHhhGg9T7ObZYpIR2C1MeY3dcs9genGmN/XLU8GMMY0LPYN9/OOMWZgE++PAkYBJCQkXLC/buYcpZRqSlVNJU+9cTd+3/YhoqItYaE/0HPtCwR36kT8S4sIjImxOkSXcUazzFggt95yXt06WwGcKSILgNQTfxwaY4xZZIxJN8akt2vXzo7wlFKeZvfu3ezevbtF2wb6B/HIHYuIG7CLAxFZlJcl8lGfCVQc+Il9tw6iQjtx2lXwG5uaxuZ/F4wxh4wxY4wxnZvxvwAdHlkpHzR69GhGjx7d4u1FhCF9p9P3tgh2tF9PQEVH1va8l2M1fuy7/Q6ObtniwGg9jz0FPw+o3x0uDjhgXzhKKWW/K357B+NHDWBzx7fxr47h0/NHUXJGHDkj7vLpZpv2FPxvgEQR6SQiQcCtwEpHBKXj4Sul7HVu7MU8Pm4yGcmvUks4G7vcSeFZ3ckdNZqSzz6zOjxLNLdZ5pvARiBJRPJEZIQxphq4B1gL7ASWG2O2OyIovaWjlHKEmIgE5oxZQHba65QFVvJdwi3kdLqUvPH3cuSjj6wOz+V08DSllNvo3bs3AOvXr3fofqtrq3n8n3dTu7En7criiC9YR9c97xP7/HNEXH21Q4/lDmy10nHL8fDrDY9sdShKKRd65JFHnLLfAL8Apt28gJcjH2PPuiMIfanwD8NMeIA4Ly36jdErfKWUT1m95RU+fDeP8wou4YziDFK2v0rC7OdpfeWVVofmMB41PLLew1fKN2VkZJCRkeHUY1yTOpzBQy7mu7j3KT4jhU2/Hcu+CRMp27DBqcd1B3qFr5RyG866h9+YzILveXLJPFKzBxJelkfazsUkLpxN6AUXOP3YzqZX+EopVU9i9Pk8OfZRvj/3dUrCO7Cp2zh2jXuY8q3bTr2xh3LLgq/t8JVSrtC+dQeeGzuHzO5vcSSsDd8mj2P7uD9T8cMPVofmFG5Z8JVSylVaB0fw3OhFHExfyZHQADZ3HcmWkRO9cvYsLfhKKZ8X6B/IE8PmU33ZFxwJqSGj4xC+Hj6RmuJiq0NzKLd8aFuvHf7ITB3hTimfsaGupUyvXr0si+Hl1Y9T/EFnwirbcG7Re1z++iz8QkIsi6clbD20dcuCf4K20lFKWeHd9QvIXh5OaFUUncvXcvWrzyD+/laH1Wwe1UpHKeWbNmzY8OtVvpVu7D2GtGEhVPgfJCu0P++PmeoVc+RqwVdKuY0pU6YwZcoUq8MA4IoeA7lsfDy1Jo8c6cM7902zOiS7acFXSikb0rpdydUTU/Cr3kdhxWW8Nmm61SHZxS0Lvna8Ukq5i6TENAZM7UVw+Q+U/HIJ8x971OqQWswtC752vFJKuZOEzt34w9RLCC/dS+2BS5n5xBSPvKfvlgVfKaXcTUxSMv0fvJCIkixC91/J9GcnUl1bbXVYp0ULvlLKbcyaNYtZs2ZZHYZN0WmpXDX8XCKOZBOd2Zc/z36Q8upyq8NqNi34Sim3kZKSQkpKitVhNOmsq3rTu39rIo7sp9Oua5jwwv0crvCM540uK/gicoOIvCQi74lIX1cdVynlOdatW8e6deusDuOU4u8cxO+6FdGmJIff7LqJexfdR35ZvtVhnVJzJzF/RUQKRGRbg/X9RGS3iGSJyKSm9mGMWWGMGQkMBW5pccRKKa81Y8YMZsyYYXUYzdJpykP0CM+gdUk+6dtv5d5XJrC3eK/VYTWpuVf4S4B+9VeIiD8wD+gPJAODRCRZRLqLyOoGr+h6mz5St51SSnks8fOj6/NPkXZkNeFHD3HJjjuZ8I+HyShw7oxd9mhWwTfGfA4UNVh9IZBljMk2xlQCbwHXG2O2GmOuafAqkOOeBj4wxmy2dSwRGSUim0RkU2FhYUvzUkopp/MLCyNp/ixS975KWPlhrtwxnInvTOVfuf+yOrRG2XMPPxbIrbecV7fOlvHAVcBAERlj60PGmEXAY8DmoKAgO8JTSinnC4qLpevsJ0jb+iJhFRX03zGSae//lRVZK6wO7ST2FHxpZJ3NngjGmDnGmAuMMWOMMQua2rF2vFJKeZLQ9HQ6ThpP2uY5hFYIN2wfw9OfPsfirYvdqoNWgB3b5gHx9ZbjgAP2hXNcvfHwHbE7pZSHWLhwodUhtFjkH/9I7M5d1Lw3l2/TJ3DLtrEslLkcOnaIh9Ifwk+sbwVvTwTfAIki0klEgoBbgZWOCUsp5YuSkpJISkqyOowWaz95Eu2TO5CydSFBlZHcsX0ky75/m8n/nkxVTZXV4TW7WeabwEYgSUTyRGSEMaYauAdYC+wElhtjtjsiKL2lo5RvWrVqFatWrbI6jBaTwEBiZ8+iXUgZydlv4X8sjhHbhrE262PGfzqeo1VHrY3Pne4vnaBTHCrlm3r37g3A+vXrLY3DXsd27mTfoNvI6nI1+6P6EhyWwQvd/0Fyu2Tm9ZlHZEikU4/vUTNe6RW+UsqThZx3HmfNmEHnbas4o/I7KsrSeGj7dew5tIs7P7iTA6UOedx52tyy4CullKdrc80faDt8OGkbXgb/XA6X9Gbqjks5dLSAwWsGk/mL6+9euGXB1wlQlFLeIPqBCYT16sUVnz9LaXAxB4qv57EdyVBbxZAPh7D5oM0+qE7hlgVfb+kopbyBBATQ4dlnCGrXjusy5lIcfJQ9h27niV0xnBnYmlEfj2J97nqXxeOWBV+v8JXyTUuXLmXp0qVWh+FQAZGRxL0wFyn+hRsL3uZwYCWb84fz9O4QEiPO5v7P7ufdzHddEotbFny9wlfKN8XHxxMfH3/qD3qYkORkzvrLY8iW/9AvPINSf8OXeSOYuecYF0WnMW3DNJf0ynXLgq+U8k3Lli1j2bJlVofhFG2uv57I22+n1btLST23hFKC+eSHwTyd/TP9z+7LrM2zmLlpJrWm1mkxaMFXSrmN+fPnM3/+fKvDcJr2f55Iq9RU4v7+NNHpwRypOZOPdlzDX/f9yO1Jg1i6YylTv5hKVa1zeuW6ZcHXe/hKKW8kQUHEzpqFX2goF7/9NOW/DaegujPrtqTycO4+7k0Zz+rs1U7rleuWBV/v4SulvFVg+2jinn+OytxcBm1/g9xOIeyr6MEXG9pwV34u03tOZ2vhVqd0znLLgq+UUt4stEcPoh9+iNJ167i/zXayzxS2lffnu3X7uLmogDU3raFLpONHC9aCr5RSFmg7ZAit+/ejeM5sJlzpT06rWr4sGcreFStok/mJU47ploOnnZCenm42bdpkdRhKKRf5+eefAYiKirI4EteoKS1j35/+RM3hwxybs5j//XsOMTU13Bg1nZixC6BDSov261GDp+lDW6V8U1RUlM8UewD/8DDi5s6htryciJnT6Xl7IofFnxVF0zjsf47Dj+eWBV8f2irlm5YsWcKSJUusDsOlgjt3psOMv1KekUHKZ2/QdkA8eyWQPUUVDj+WPVMcKqWUQ50o9kOHDrU0DleLGDCAo1sy+OW1pdzyXAqF0y6jS3S4w4/jllf4Sinla9o//BCtUlLIf+R/iC8tcMoxXFbwReQ8EVkgIu+IyFhXHVcppTzB8U5ZzyPBwfx4333UHrWo45WIvCIiBSKyrcH6fiKyW0SyRGRSU/swxuw0xowB/gSc9PRYKaV8XWBMDB2emUlNSSlVP/7o8P039wp/CdCv/goR8QfmAf2BZGCQiCSLSHcRWd3gFV23zXXAF4BzGpkqpZSHC//d7+i89kOCExMdvu9mPbQ1xnwuIh0brL4QyDLGZAOIyFvA9caYJ4FrbOxnJbBSRN4H3mhp0Eop77RmzRqrQ3ALfsHBTtmvPa10YoHcest5wEW2PiwivYGbgGDA5lkVkVHAKICEhAQ7wlNKeZrQ0FCrQ/Bq9hR8aWSdzW67xpj1wPpT7dQYs0hEfgKuDQoKuqDF0SmlPM6LL74IwLhx4yyOxDvZ00onD6g/NU0c4JDh3bTjlVK+afny5SxfvtzqMLyWPQX/GyBRRDqJSBBwK7DSEUHp0ApKKeV4zW2W+SawEUgSkTwRGWGMqQbuAdYCO4HlxpjtzgtVKaWUPZrbSmeQjfVraOIBbEsZY1YBq9LT00c6et9KKeWr3HJoBb2lo5RSjufW4+GLSCGwH2gD1K/+9Zdt/RwF/OyAMBoe257P2nq/qfxOteyJOTdnnafm3Nxz3Ng6zdn3crb1nr35nm2MaXfSWmOM27+ARbaWm/h5kzOObc9nbb3fVH6nmb9H5NycdZ6ac3PPseasOTf1nqPybfhyy1s6jVjVxLKtn511bHs+a+v9pvI71bIn5tycdZ6ac3PPcWPrNGffy/lUvw+HcutbOvYQkU2mkSm+vJnm7Bs0Z+/nrHw95Qq/JRZZHYAFNGffoDl7P6fk67VX+Eoppf6bN1/hK6WUqkcLvlJK+Qgt+Eop5SN8suD74vy6InKDiLwkIu+JSF+r43EFETlHRBaLyDtWx+IsIhImIq/WndvbrY7HFXzhvDbksO+vMxr3O/MFvAIUANsarO8H7AaygEnN3JcfsNjqnFycc6QP5vyO1fk4K3dgMHBt3c/LrI7dlefb086rg3K26/tredIt+CVdBqTV/yUB/sBe4BwgCPiO4/PsdgdWN3hF121zHbABuM3qnFyVc912zwJpVufk4pw9qjCcZu6TgZS6z7xhdeyuyNlTz6uDcrbr+2vPjFeWMD44v64jchYRAZ4CPjDGbHZuxPZz1Hn2RKeTO8cnIooDMvDgW7SnmfMO10bnHKeTs4jsxAHfX4/9B9JAY/Prxtr6sIj0FpE5IrIQJwzv7CKnlTMwHrgKGCgiY5wZmBOd7nk+U0QWAKkiMtnZwTmZrdz/CdwsIvNxcrd8CzSas5ed14ZsnWeHfH897grfBqfMr+vmTjfnOcAc54XjEqeb8yHAU/+4NdRo7saYMmCYq4NxEVs5e9N5bchWzg75/nrLFb7T5td1Y5qzb+R8gi/mrjk7OGdvKfhOm1/XjWnOvpHzCb6Yu+bs6JytflLdgifbbwI/AVUc/2s4om79AGAPx59wT7U6Ts1Zc9bcNWd3y1kHT1NKKR/hLbd0lFJKnYIWfKWU8hFa8JVSykdowVdKKR+hBV8ppXyEFnyllPIRWvCVUspHaMFXSikfoQVfKaV8xP8BYJBaiP7nEzQAAAAASUVORK5CYII=\n",
"text/plain": [
"