\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.19?urlpath=lab/tree/tutorials/api/astro_dark_matter.ipynb)\n",
"- You may download all the notebooks in the documentation as a\n",
"[tar file](../../_downloads/notebooks-0.19.tar).\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 spatial and spectral models\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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:59.116454Z",
"iopub.status.busy": "2021-11-22T21:09:59.115718Z",
"iopub.status.idle": "2021-11-22T21:09:59.781879Z",
"shell.execute_reply": "2021-11-22T21:09:59.782052Z"
}
},
"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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:59.784530Z",
"iopub.status.busy": "2021-11-22T21:09:59.784238Z",
"iopub.status.idle": "2021-11-22T21:09:59.785430Z",
"shell.execute_reply": "2021-11-22T21:09:59.785605Z"
}
},
"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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:59.788121Z",
"iopub.status.busy": "2021-11-22T21:09:59.787790Z",
"iopub.status.idle": "2021-11-22T21:09:59.790345Z",
"shell.execute_reply": "2021-11-22T21:09:59.790528Z"
}
},
"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": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:59.805948Z",
"iopub.status.busy": "2021-11-22T21:09:59.805604Z",
"iopub.status.idle": "2021-11-22T21:10:00.155951Z",
"shell.execute_reply": "2021-11-22T21:10:00.156160Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"LOCAL_DENSITY: 0.3 GeV / cm3\n",
"DISTANCE_GC: 8.33 kpc\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABXDElEQVR4nO3dd1zV1f/A8dfhsmQ6AFFAGW5lCahgImriSHFkmjO1nGV9G1ZW37Ky8UsbX1tmroaVpWWuyhyICm5xD2QpLnCAyOZyfn9cQRwoyLgXOE8f9wH3cz/jfUDe93M/n3PeR0gpURRFUWo+I30HoCiKolQNlfAVRVFqCZXwFUVRagmV8BVFUWoJlfAVRVFqCZXwFUVRagljfQdwL3Z2dtLV1VXfYSiKUkVOnDgBQMuWLfUcSfW2d+/eS1JK+9uXG3TCd3V1Zc+ePfoOQ1GUKhISEgJAeHi4XuOo7oQQiXdbbpAJXwjRH+jfrFkzfYeiKEoVeuONN/QdQo0mDHmkrb+/v1Rn+IqiKGUjhNgrpfS/fbm6aasoisGIjo4mOjpa32HUWAZ5SUdRlIqRl5dHUlIS2dnZ+g6lVK5evQrAsWPH9BxJ9WBubo6zszMmJialWl8lfEWpwZKSkrC2tsbV1RUhhL7DuS8jI91FB9VL5/6klFy+fJmkpCTc3NxKtY1BXtIRQvQXQsxPS0vTdyiKUq1lZ2fToEGDapHslbIRQtCgQYMyfXozyIQvpVwtpZxoa2v7QNufj03j8tnrFRyVolRPKtnXXGX93Rpkwi+v7ctj+GXWLjYsPsq1S1n6DkdRajUhBC+++GLR8zlz5jBz5kwAZs6ciZOTEz4+Pvj4+PDxxx+zceNGBg4cWLT+Bx98QPEu2qtXryYsLAzQjdXx9PTE29ub0NBQLly4UKbYpk+fTtu2bZk+fTrz5s3j+++/B2Ds2LEsX778AVtsuGrkNfx+T3uz759EDoYnEbPnIm27OOHf1xULG1N9h6YotY6ZmRm///47M2bMwM7O7o7Xn3/+eV566SUArl+/TkpKCm+//XbR61FRUdjY2JCcnIyDgwORkZF07ty56PXNmzdjZ2fHa6+9xvvvv8/cuXOLXpNSIqUsujdwu2+++YaUlBTMzMwqqrkGrUae4ZtbmRD0aDNGvRNIq6BGHI44yw//jWLHn7HkZOXrOzxFqVWMjY2ZOHEin3766X3XtbKyws3NDVtbW06dOgXA2bNnefTRR4mMjAQgMjKSoKCgO7YNDg7m1KlTJCQk0Lp1a6ZOnUr79u05c+YM06dPp127dnh6erJs2TIAwsLCyMjIoGPHjixbtoyZM2cyZ86cO/a7d+9eunbtip+fH7169eL8+fPl+XHoVZWe4QshLIEI4C0p5ZrKPp5VPTO6jWyF78NN2Lk6jr1/JXI44iztezXFK8QZY1NNZYegKAbj7dVHOHruWoXus01jG97q3/a+6z399NN4eXnx8ssv3/Hap59+yo8//gjoLvE8/PDDBAUFERkZiVarpXnz5nTq1Il//vmHfv36cfDgQQICAu7Yz5o1a/D09AR0NXkWL17MV199xYoVK4iOjubAgQNcunSJgIAAgoODWbVqFVZWVkX9/gsvMxWXl5fHtGnT+PPPP7G3t2fZsmW8/vrrLFq0qAw/JcNRrjN8IcQiIUSyEOLwbct7CyFOCCFOCSFeLfbSK8Cv5Tnmg6jb0IJeT7Vj6GsBNHS1Jer3WH58cwdHtp6lQFtQ1eEoSq1jY2PDmDFjbrncUuj5558vGnDVunVrzp49S+fOnYmMjCQyMpLAwEA6dOjAzp072b9/Py1btsTc3Lxo+27duuHj48O1a9eYMWMGAE2bNqVTp04AbNu2jeHDh6PRaGjYsCFdu3Zl9+7dpYr7xIkTHD58mJ49e+Lj48OsWbNISkqqgJ+IfpT3DH8J8AXwfeECIYQG+BLoCSQBu4UQq4DGwFHA/M7dVA37Jtb0n+bN2ZNX2bEylvClJ9j/72k6hrnTrL0Dwkj1ZlBqrtKciVem//znP7Rv355x48bdd92goCA+//xztFotEyZMwNramuzsbMLDw2+5fg83r+EXSk1NxdLSsuh5ecrHSClp27YtUVFRD7wPQ1KuM3wpZQRw5bbFHYBTUso4KWUu8AswAOgGdAJGABOEEHq7f+DUoh6Dp/vRd4onGmMj1i84wq8f7CbxyOVy/edQFKVk9evXZ+jQoSxcuPC+67Zp04Zz586xdetWfH19AfDx8WHevHl3vX5/L8HBwSxbtgytVktKSgoRERF06NChVNu2bNmSlJSUooSfl5fHkSNHynR8Q1IZSdcJOFPseRLgJKV8XUr5H+An4Fsp5V2vpQghJgoh9ggh9qSkpFRCeEXHwc3bnmFvdODhsa3JycxnzecHWPnJfi7EqQFfilIZXnzxRS5dunTf9YQQdOzYETs7u6KyAYGBgcTFxZU54Q8aNAgvLy+8vb3p3r07H330EY6OjqXa1tTUlOXLl/PKK6/g7e2Nj49P0c3j6qjc1TKFEK7AGilluxvPHwN6SSmfuvF8NNBBSjmtDPssLI88ISYmplzxlZY2v4AjW8+x568Esq7l4uplR6cB7jRwsqqS4ytKZTh27BitW7fWdxilpiZAKbu7/Y5LqpZZGb10kgCXYs+dgXOVcJwKpTE2wqubM62DGnFg4xn2r0/kl1m7aNGhIR37u2NjV0ffISpKjefi4nL/lZQHVhmXdHYDzYUQbkIIU+BxYFVZdlDe0grlYWKmwb+vK6PfC8K3ZxNi96Ww9K0dRPxyksxruVUej6LUJhYWFlhYWOg7jBqrXGf4QoifgRDATgiRhK5//UIhxDPAP4AGWCSlLNNdDkOY8crc0oSgwc3w6ubC7nXxHI44y7HIc3j3cME3tClmdWrkIGVF0atr13TjBGxsbPQcSc2kZrwqpdSLmexcHcepPcmYWRqrwVtKtaCu4dd8ZbmGb5ClFQyxPPKtg7dsbhm8pVWDtxRFqQYMMuHr8xr+/egGb/kw8AVfrOubEb70BD+/vZOYPReRBYb7aUlRFMUgE355z/B/2nmaRdviycipvEJpavCWopSORqMpKn/s4+PDhx9+CMBTTz3F0aNHK+w4CQkJ/PTTT/ddLzw8HFtbW3x9fWnduvUtlTlLIyUlhY4dO+Lr68vWrVvp27cvqampgK74myEzyDuPUsrVwGp/f/8JD7L9tlMprDt0gf9tjGFUpyY8EeSKg3XFV3QoHLzV1NOOmN0X2bU6jjWfH6Bx87oEDvLA0d3wPqEoSlWrU6fOXScmX7BgQYUepzDhjxgx4r7rdunShTVr1pCRkYGPjw/9+vXDz8+v6PX8/HyMje+eHjdu3EirVq347rvvivZVXRjkGX55fTXSj9+nBhHk0YCvwmN56MPNvLL8IKeS0yvleEZGgpYdHRkxsxPBj7fg6sVMVny0l7VfHVQzbylKCUJCQijslGFlZcXrr7/OY489xujRo7l48SKgm+yk8Gz64YcfLlq+ZcuWok8Mvr6+pKen8+qrr7J161Z8fHz49NNPyc7OZty4cXh6euLr68vmzZvviMHS0hI/Pz9iY2OZOXMmEydOJDQ0lDFjxpCYmEiPHj3w8vKiR48enD59mujoaF5++WXWrVuHj48PWVlZuLq63nX08OzZswkICMDLy4u33nqrEn+SpWeQZ/gV0S2zfZN6fD3Kj4RLGSzYFsdve5JYtucMD7d2YEIXdzq41a/wqd80xkZ4hjjTKrARBzadYf/60/wyaxctOzjSob+bGryl6Ndfr8KFQxW7T0dP6PPhPVfJysrCx8en6PmMGTMYNmzYLetkZGTQqVMn3nvvPV5++WW+/fZb3njjDR566CF27NiBEIIFCxbw0Ucf8fHHHzNnzhy+/PJLOnfuzPXr1zE3N+fDDz9kzpw5rFmjq7z+8ccfA3Do0CGOHz9OaGgoJ0+evOW4ly9fZseOHfz3v//l6NGj7N27l23btlGnTh369+/PmDFjeOKJJ1i0aBHPPvssK1eu5J133mHPnj188cUXJbZ5/fr1xMTEsGvXLqSUhIWFERERQXBwcFl+uhXOIBN+eS/pFOdqZ8msgZ48/3ALvo9K5PuoBDYc24G3S10mBbvTq60jmgqukmlipsG/jyvtgp10M29tVjNvKbVXSZd0ijM1NaVfv36kpqbSqlWrono1SUlJDBs2jPPnz5Obm4ubmxsAnTt35oUXXmDkyJEMHjwYZ2fnO/a5bds2pk3TVXRp1aoVTZs2LUr4hUXZjIyMePXVV2nbti2//fYbYWFh1KmjOzGLiori999/B2D06NF3reVfkvXr17N+/fqiwm/Xr18nJiZGJfyq0sDKjOd7tmByVw+W70tiwdY4pi7dR9MGFjz1kBtD/FyoU8F96tXgLcWg3OdMXJ9MTEwQQnDx4kXS09PJz9d1uJg2bRovvPACYWFhhIeHF01S8uqrr/LII4+wbt06OnXqxIYNG+7Y5706TxRew79d8bLKtyvLFQEpJTNmzGDSpEml3qYqGOQ1/Mrsh1/HVMPoTk3Z9GIIX49sTz0LU/775xGCPtzIJ/+e5PL1nAo/ZuHMWyPe6oiblx17/0rkhzci2bc+kfxcbYUfT1FqirS0NJycnACKbpICxMbG4unpySuvvIK/vz/Hjx/H2tqa9PSb9+mCg4NZunQpACdPnuT06dNlGtAVFBTEL7/8AsDSpUt56KGHSr1tr169WLRoEdev6+7hnT17luTk5FJvX1kM8hSzIi/plERjJOjj2Yje7RzZk3iVb7bEMXdjDN9siWWInzNPdXHHza7kd/sHUbehBaFPtcM3NJ0df8YR9XssBzclEfCIK62DGmGkMcj3X0Upl9uv4ffu3buoa+b9zJw5k8ceewwnJyc6depEfHw8AJ999hmbN29Go9HQpk0b+vTpg5GREcbGxnh7ezN27FimTp3K5MmT8fT0xNjYmCVLlpRpsvK5c+cyfvx4Zs+ejb29PYsXLy71tqGhoRw7dozAwEBAd1P6xx9/xMHBodT7qAyqtEIxp5Kvs2BrHL/vO0teQQGhbRoyMdgDv6b1KuV452KuEvVHHBfi0rB1qKNm3lIqnCqtUPNV+9IK+tLMwYoPH/Vi2yvdmBriQVTsZR79OpIhX0ey/sgFCip4JG3j5vUYPL09fad6qcFbiqJUOnWGfw8ZOfks232GhdviOZuahbu9JRO6uDPI1wlzk4q9wVtQIInZdYGdq+NJv5ytBm8pFaK6neHn5upKkJuaqp5spVWWM3yDTPj6mPHqXvK1Baw7fIH5EbEcPnsNOyszxgY1ZVSnptS1qNj/mNr8Ao5uO8fudWrmLaX8qlvCV8qu2if8Qvo+w7+dlJKo2MvM3xpH+IkU6phoGBbgwpMPueFSv2InbcjL0XJw8xn2/XOa3Ox8NXhLeSDVLeFfuXIF0E14rpSOvqc4rLGEEAQ1syOomR0nLqQzPyKOpTt1g7n6ejZiYrA7Xs51K+RYJmYa/Hq70rbLbYO3HmqMX19XLG1L39tAUaqLlJQUQCX8yqIS/gNq6WjNx0O9md6rJYu3x/PTztOsOXieTu71mRTsQUhL+wop3VB88NaedfEc3nqOY1Hn8e7ugm9oE8wsTCqgNYqi1Aaql045OdqaM6NvayJndOf1vq1JvJzJuCW76fVZBL/tOUNufsVMjmJVz4yQG4O3XL3s2Pt3Ij+8EaUGbykG70FLBr///vtF3yckJNCuXbuKCumBFC/25urqiqenJ97e3oSGhnLhwoUy7Wv69Om0bduW6dOnM2/ePL7//nsAxo4dy/Llyys89kJVlvCFEK2FEPOEEMuFEFOq6rhVxdrchAnB7kS83I1PhnpjJATTlx+ky0eb+Do8lrSsvAo5zi0zb7nZ6mbe+m+UmnlLqXGKJ/zyKizVUJE2b97MgQMH8Pf3vyNWKSUFBSX/PX7zzTfs27eP2bNnM3nyZMaMGVPh8d1NuRK+EGKRECJZCHH4tuW9hRAnhBCnhBCvAkgpj0kpJwNDgTtuJtQUJhojBrd35q/nuvDd+A40c7Di//4+TucPNzFrzVHOpWZVyHF0M295M+hFX6wb1FEzbykG7/z58wQHB+Pj40O7du3YunUrAD///DOenp60a9eOOXPmALpaOYUjdEeOHAmAVqtlwoQJtG3bltDQULKydH9LsbGx9O7dGz8/P7p06cLx48cB3dnyCy+8QLdu3XjllVcYO3YsU6ZMoVu3bri7u7NlyxbGjx9P69atGTt2bFGcU6ZMwd/fn7Zt25aqrHFwcDCnTp0iISGB1q1bM3XqVNq3b8+ZM2eYPn067dq1w9PTk2XLlgEQFhZGRkYGHTt2ZNmyZcycObOo3cXt3buXrl274ufnR69evTh//vyD//BvKO81/CXAF8D3hQuEEBrgS6AnkATsFkKsklIeFUKEAa/e2KZGE0LQtYU9XVvYc/hsGvMj4lgcmcCSyAT6ezdmYrA7rRvZlPs4hYO3Eg5dZsfKWNYvOMI+l0Q6DfSgSZuKLwGtVF//t+v/OH7leIXus1X9VrzS4ZVSrfvTTz/Rq1cvXn/9dbRaLZmZmZw7d45XXnmFvXv3Uq9ePXr27MmhQ4f48MMP+eKLL4qqbCYkJBATE8PPP//Mt99+y9ChQ1mxYgWjRo1i4sSJzJs3j+bNm7Nz506mTp3Kpk2bAF0NnQ0bNqDRaBg7dixXr15l06ZNrFq1iv79+7N9+3YWLFhAQEAA0dHR+Pj48N5771G/fn20Wi09evTg4MGDeHl5ldiuNWvW4OnpCehGCi9evJivvvqKFStWEB0dzYEDB7h06RIBAQEEBwezatUqrKysitpWWBCuuLy8PKZNm8aff/6Jvb09y5Yt4/XXX2fRokWl/+XcRbkSvpQyQgjhetviDsApKWUcgBDiF2AAcFRKuQpYJYRYC9x1LjIhxERgIkCTJk3KE57BaOdky9zhvrzcuyWLtiXwy+7T/LH/LF2a2zEp2IPOzRqUKzELIXDzsqNpuwbE7L7IzlVq5i3F8AQEBDB+/Hjy8vIYOHAgPj4+bNq0iZCQEOzt7QEYNWoUkZGRDBky5I7t3dzcimry+Pn5kZCQwPXr14mMjOSxxx4rWi8n52YBxMceewyN5uYgyf79+yOEwNPTk4YNGxYl6rZt25KQkICPjw+//vor8+fPJz8/n/Pnz3P06NG7Jvxu3bqh0Wjw8vJi1qxZpKam0rRpUzp16gToyjMPHz4cjUZDw4YN6dq1K7t37yYsLOy+P6sTJ05w+PBhevbsCeg+3TRq1Oi+291PZfTScQLOFHueBHQUQoQAgwEzYF1JG0sp5wPzQdcPvxLi0xvneha82b8Nz/VoztJdiSzensCohTtp08iGicHuPOLVCJNyFFArnHmrmZ8DR7aeY89fCaz4aK8avKUAlPpMvLIEBwcTERHB2rVrGT16NNOnT8fG5tZPuenp6UWXam5XvPCZRqMhKyuLgoIC6tatW2K9/dvLHRfuw8jI6Jb9GRkZkZ+fT3x8PHPmzGH37t3Uq1ePsWPHkp2dfdd9b968GTs7u6LnqamptxyvPGOcpJS0bduWqKioB97H3VTGTdu7napKKWW4lPJZKeUkKeWX99xBJZZHNgS2FiZMDWnGtle68dGjXuRqC/jPsmhCZoezYGsc18s5+brG2Aivbs6MeqcTHcPcOXfyKr/M2sWGxUe5dqli7iEoSlklJibi4ODAhAkTePLJJ9m3bx8dO3Zky5YtXLp0Ca1Wy6+//lp01m1iYkJe3r07O9jY2ODm5sZvv/0G6BLlgQMHHjjGa9euYWlpia2tLRcvXuSvv/564H0FBwezbNkytFotKSkpRERE0KFDh1Jt27JlS1JSUooSfl5eHkeOHHngWApVxhl+EuBS7LkzcK4SjlPtmRlrGBrgwhA/ZzafSOabiDhmrT3G3I0xjOzUlHFBrjjYPPjk66bmxvj3LTbzVriaeUvRn/DwcGbPno2JiQlWVlZ8//33NGrUiA8++IBu3bohpaRTp0706NEDgIkTJ+Ll5UX79u157733Stzv0qVLmTJlCrNmzSIvL4/HH38cb2/vB4rR29sbX19f2rZti7u7O507d36g/QAMGjSIqKgovL29EULw0Ucf4ejoWKptTU1NWb58Oc8++yxpaWnk5+fzn//8h7Zt2z5wPFABpRVuXMNfI6Vsd+O5MXAS6AGcBXYDI6SUZX57MrTSClVh/+mrfLs1jr8PX8DYyIiBvrobvM0crMu97+tXc9i9Np5jkefRmBjh3d1ZzbxVw1W30gqqPHLZVVlpBSHEz0AIYCeESALeklIuFEI8A/wDaIBFZU32FTGJeXXl26QeX430I/FyBgu2xvPb3jP8uieJHq0cmBDsTsdyTL5uVc+MbqNa4duzCTtXxbH3r0QOR5zFr5crniFOGFfwFI+KohgWgyyeVu5qmQVaMKoZyetKRi4/3Jh8/XJGLt7OtkwM9qB3u/JPvp5yOp0dK2M5ffQKlnXN1MxbNZA6w6/5VLXMZaPg7D6wbwl2LXVf7VvpvlpUz6JM2XlaVuxLYsHWeOIvZdCkvgVPdXHjsQqYfP3siatErYzlYvw16ja0oGOYOx6+9mrmrRqguiV8rVZXJqR4V0rl3lTC3/sdJG6HlBNw6STkZd58zdJBl/gdWuse9q3BoRXUqZxpDCuatkDy79GLzI+IZd/pVOpZmDA60JUnApvSwOrBK2hKKYk/cImdq+K4ci4D+ybWdBrgjosavFWtVbeEr5RdtU/4FToBSkEBXEvSJf+U45B8HFKO6Z7nXr+5nnVjaNgGHNpAw7a6r/YtwdhwyxDvSbjCvC1xbDh2ETNjIx71c2ZCOSdfLyiQnNx1gV03Zt5yalGXTgPV4K3qqrol/OTkZAC9T/ZdnVT7hF+oUnvpSAmpp2+8CRyFi0ch+RhcOgFa3TRrCA3YNde9ATRsp3s4eoK1IxjQWW9sim7y9RX7zpKnLZx83R2/pg9++UqbV8CRbWfZsy6BrPQ83Lzt6DjAnQaN1eCt6qS6JXx1Db/sqn3C1+sUh9p8uBILF4/ceBzWvRmknb65jkUDXeJ39ARHL93XBs1Bo9/ujSnpOXwXmcAPOxJJy8rDr2k9Jga707N1Q4we8Hp8bnY+BzedYd/60+TlaGnV0ZGAfmrmrerCEBK+RqPB09MTKSUajYYvvviCoKCgu65bUsIPCQlhzpw5+PuXve5iamoqP/30E1OnTgUoKnLWsmVLcnNzCQ4O5quvvsLIqHSdFXJycnjkkUe4dOkSM2bM4N9//+WFF16gTZs2uLq6smfPnltG4Fa2aj/jlZRyNbDa399/QpUfXGN84yZvS2g3+ObyrFTdJ4ELh+HCQbhwCHbOB+2Nuh3G5rpPAo5e0MgbGvvoLgtV4SUhe2szXurVkikhHvy65wwLtsYz6Ye9uNtb8tRD7gxuX/bJ13WDt9xoG+zEvr8TORR+lpO7L9Iu2Am/PmrwlnJ/derUKSp98M8//zBjxgy2bNlS6u0Lb+Q+CK1WS2pqKl999VVRwgfw8PAgOjqa/Px8unfvzsqVKxk8+Obfe35+PsbGd0+P+/fvJy8vr6hNw4YNe+D4qppBJnyDVKcuNA3SPQpp8+FyDJw/qHsTOH8ADv8OexfrXjcy0d0YbuwDjXygsa/uspBx5SZJSzNjxnV2Y3SnpkWTr7/2xyE++fcETwS6Mjqw7JOv17EypfOQ5nh1d2HP2ngOhSdxNPI8Pg+74PtwE0zV4C2lFK5du0a9eroOEuHh4cyZM4c1a9YA8Mwzz+Dk5MTgwYNxdXVl/PjxrF+/nmeeeaZo+4KCAsaNG4eLiwtvv/02r776KuHh4eTk5PD0008zadIkwsPDefvtt2nUqBHR0dF4eXkRGxuLj48PPXv25Omnny7an7GxMUFBQZw6dYolS5awdu1asrOzycjIYPny5YwfP564uDgsLCyYP38+jo6OjBo1ipSUFHx8fFixYgVPPvnkXT99/Pjjj8ydO5fc3Fw6duzIV199pffeR+qvtDw0xjd7+3jfeJeXEq4m6JL/+Wg4Fw3HVsO+GxWkNaa6TwKN24NTe91X+5aVMm7AWGNEmHdj+ns1Iir2Mt9ExPHxvyf5Kjz2gSdft65vTrfRrfHp2YSdq+LZszaBw+Fn8evTlHZdnTAu4ycIpepceP99co5VbHlks9atcHzttXuuU1jXPjs7m/PnzxeVLr4fc3Nztm3bBsC8efPIz89n5MiRtGvXjtdff5358+dja2vL7t27ycnJoXPnzoSGhgKwa9cuDh8+jJubGwkJCRw+fPiWUsuFMjMz2bhxI++88w4XL14kKiqKgwcPUr9+faZNm4avry8rV65k06ZNjBkzhujoaBYsWHDLG9XdHDt2jGXLlrF9+3ZMTEyYOnUqS5curbKJTkpikAm/Wo+0FQLqu+kebQfqlkkJqYlwbr/ucXYfHPoN9izUvW5iqTv7d/YDJz9wDgCbxhUY0s3J149fuMb8iDh+3HFz8vVJwR54OpetF049R0t6T2xHcuI1dqyMZfvyUxzYeIaAfm606uSoBm8pRYpf0omKimLMmDEcPnz4rus2atSo6Pr97ZdKJk2axNChQ3n99dcBWL9+PQcPHiyaEjAtLY2YmBhMTU3p0KEDbm5uJcZUeMYvhGDAgAH06dOHJUuW0LNnz6IJ1Ldt28aKFSsA6N69O5cvX6a0BR03btzI3r17CQgIAHRveobQ88ggE75er+FXBiGgnqvu0XaQbllBAVw+BWf3wrl9kLQHor6CghvVAa0bg7O/Lvm7dNDdFzAp/43SVo42fDLUh+m9WrJke0LR5OuB7g2Y2NWdkBZlm3zdoakNYc/5knT8ClEr49j8w3Gi/z1NxzB33H0rZiJ3pWLc70y8KgQGBnLp0iVSUlIwNja+ZRrA28sQ317aOCgoiM2bN/Piiy9ibm6OlJLPP/+cXr163bJeeHj4HdvervAa/u3uV964tP+fpZQ88cQTfPDBB6Vav6qo0zB9MTIC+xbgMxz6zoaJm+G1s/DURuj9f7p7BecPwL//hUW94AMX+LYH/D0DjvwB18o33Vkj2zrM6Nua7TO681rfVsRfymDc4t30/mwry/cmlXnydedW9Rnyih99JulK2/49/zDLP9zDmWNXyhWnUrMcP34crVZLgwYNaNq0KUePHiUnJ4e0tDQ2btxIWlpaiROCP/nkk/Tt25fHHnuM/Px8evXqxddff11UQvnkyZNkZGTcsZ21tTXp6elljjU4OJilS5cCujcROzu7O+r3l6RHjx4sX768aFzBlStXSExMLHMMFc0gz/BrLWOzG2f1/sBk3bLryZC0G87shDO7Yc8i2PGV7rW6TcClIzQJ1D3sW+neSMrAxtyEicEejA1yY/WBc8yPiOOl3w4w558TjOvsyvCOTbAxNynVvoQQuPva4+ptx4kdF9i1Oo5V/4vGuVU9Og30oKFr+ad0VKqfwmv4oDvz/e6779BoNLi4uDB06FC8vLxo3rw5vr6+ZGZm3vOyyQsvvEBaWhqjR49m6dKlJCQk0L59e6SU2Nvbs3Llyju2adCgAZ07d6Zdu3b06dPnlpu29zJz5kzGjRuHl5cXFhYWfPfdd6Vuc5s2bZg1axahoaEUFBRgYmLCl19+SdOmTUu9j8pgkP3wC9XG8sj3lZ+r6xJ6ZofuTeD0Drh+Ufeaue3N5N+0s653kKZ0ybqQlJItJ1OYHxFHZOxlrMyMGdGxCeM6u9LItmyXlPLztBzecpa9fyeSfT0PD197Og5wp57jg48EVsrGEPrhl4UaeFV2auBVbVLYK+h0FCRG6t4ALt/4mZlY6O4BNO0Mrg/pPjmUYVzAoaQ05m+NY+3BcxgJQZhPYyYFe9DSsWy1+XOz8oneeIbof0+Tn6ulVVAjAh5xw7r+g0/uopSOSvg1X7VP+IXUGf4Dup588w0gYbtutDBSNzjMOQBcu4BbF3DyL9WYgDNXMlm4LZ5lu8+QlaclpKU9E4PdCXQv2+TrWem57P0rkUMRSQgEniFO+PV2xdyqbJ9ClNJTCb/mUwlfuVXWVUiMgoRtkLBVd0kICcZ1oEkncAsG9666wWH3GA+QmpnLjzsSWRKZwKXruXg62TIx2J0+7RwxLkM3zGuXs9i9Op4TOy9gYqbBp2cTvHu4YGqubilVtOqW8E+ePAlAixYt9BxJ9aESvnJvmVd0Z//xEbo3gOSjuuVmtrozf7eu4B6iKxx3lzP47Dwtv+87y4KtccRdysClfh2e7OzG0AAXLExLn7Qvn7vOzj/jiD9wiTrWJvj3daXtQ05oTFTnsYpS3RK+UnYq4Stlcz1Zl/zjt0BcuK6KKICNsy7xe3QD925g2eCWzQoKJP8eu8j8iDj2Jl6lroUJozs15YkgV+zKUJv/QlwaO1bGcvZkKtYNzOnY343mHRwfuOCbcpNK+DWfQSZ8IcRA4BHAAfhSSrn+ftuohK8nV+J0iT92s+5NIDsNELpePx49oFkP3b2AYj2A9iZe4Zstcfx77CImGiOGlLE2v5SSM8euEPVHLJfOXKd+Y0s6DXDH1ctODd4qh+qW8M+dOwdA48YVN9K8pitLwi/XZ2chxCIhRLIQ4vBty3sLIU4IIU4JIV4FkFKulFJOAMYC1ae8XG1U3x38x8OwH+DleN1gsG6vgcYMtn0Ki/vAR+66qST3fgfXzuHXtD7zx/iz4YWuPNreieV7k+j+cTiTftjD3sSr9z2kEIImbRowdEYAoU+1RZtXwLqvD/H77H2ci0mt/DYrlUYIwejRo4ue5+fnY29vT79+/e5YNz09/YEGSd3NkiVLsLe3x8fHhzZt2vDtt9+Wafvjx4/j4+ODr68vsbGxRSWdExISaNeuXYXEWNXKe5dsCfAF8H3hAiGEBvgS6AkkAbuFEKuklDcuFPPGjdeV6sBIc3MwWNeXdWWi4yPg1L9waqOuMBzoqoA274lH8158MCCAF3q2LKrN/8+Ri/jfqM3/8H1q8wsjQXP/hrj72nNs+3l2r43nj4/30bRdAzoNdMfOuWxdQhX9s7S05PDhw2RlZVGnTh3+/fdfnJycKmTf9ypjDLp6PF988QXJycm0bduWsLAwGjZsWKrtV65cyYABA3j77bcBiIyMrJCY9alcZ/hSygjg9rHzHYBTUso4KWUu8AswQOj8H/CXlHJfSfsUQkwUQuwRQuxJSUkpT3hKZahTF9qEQdjn8PwRmBIFPd/RzQkc+Tks7g2zPbBf/zQvNT5E5HPteat/G86nZTPxh708/OkWftl1muy8e9c412iMaBfsxKh3Awkc5MGFuDSWvbeb9QuPkJaSec9tFcPTp08f1q5dC8DPP//M8OHDi167cuUKAwcOxMvLi2HDhhV1zSy+vFOnThw8eBDQjYCdOHEioaGhjBkzhpSUFB599FECAgIICAhg+/btdxzfwcEBDw8PEhMTGTt2LC+88ALdunXjlVdeITo6mk6dOuHl5cWgQYO4evUq69at47PPPmPBggV069YNACurO2d702q1TJ8+nYCAALy8vPjmm28q/GdXkSqjH5wTcKbY8ySgIzANeBiwFUI0k1LOu9vGUsr5QojzQH9TU1O/SohPqShC6OYBbtgGOj+nu9YfuwlO/gMx6+HQb1gKDeOaBjGmSy82Sz8+3afl1d8P8fG/Jxkb5Mqojk2xtSi5H76JqYb2vZrS5qHG7F9/moObzhC7N5k2XRrj39cVS1vDnXPY0Gz99SSXzly//4plYOdiRZeh9+9C+fjjj/POO+/Qr18/Dh48yPjx49m6dSsAb731VlEZ4u+++45XXnmFsLCwW5YXL08MsHfvXrZt20adOnUYMWIEzz//PA899BCnT5+mV69eHDt27Jbjx8XFERcXR2EF3pMnT7JhwwY0Gg1eXl58/vnndO3alTfffJO3336bzz77jMmTJ2NlZcVLL71UYrsWLlx41xLN96rUqU+VkfDv9nldSinnAnNLs4MaVy2ztjC31VUDbTsICrS6MtAn1sHJv9H8+wYPAz3sW5Hk142FKa2Y808WX20+xbCAJjzZxQ2nuiWXbjC3NCFwkAde3ZzZvS6BI1vPcTzqPN49XPANbYqZmoDFoHl5eZGQkMDPP/9M3759b3mteBnizp07k5qaSlpa2j3LE4eFhVGnju7/y4YNGzh69GjR/q5du1Z0H2DZsmVs27YNMzMzvvnmm6LSx4899hgajYa0tDRSU1Pp2rUrAE888QSPPfZYqdtVUonm2pTwkwCXYs+dgXNl2UG1roev6BhpwCVA93j4LV35hxN/I06sxeXofGZKLa/Xb8hO007M39GG7lFt6OPlwsRgD9o0LrnImmVdM0JGtMSnhwu7Vsex969EDkecxa+XK54hThibqglYSlKaM/HKFBYWxksvvUR4eDiXL18uWl68p2CzZs0wNjZGCHHP8sTFyxgXFBQQFRVV9AZQXOE1/Nvdr3xyaZVUotlQVcYIl91AcyGEmxDCFHgcWFUJx1Gqk3qu0GkyPLEapp+CQfMxadqBhzL+5XuTD4iuM5WQY2/xyRefMX7BVrafunTXP/hCdRtaEPpUO4a+FkDDpjZE/n6KpW/t4Oj2cxRoy1baWaka48eP580338TT0/OW5SWVIS5teeLQ0NBbkvrd6tyXxNbWlnr16hVdXvrhhx+KzvZLo7Qlmg1Fuc7whRA/AyGAnRAiCXhLSrlQCPEM8A+gARZJKY+UZb/qkk4NZ1FfNyWk9zDIzYTYTdQ5tpoBJ9YxUGwhI+lLNn7ny/u2IXiHDKG3r3uJpRvsm1jT/1kfkk5cZcfK2JsTsAxwx91HTcBiSJydnXnuuefuWF68DLGxsTEfffTRHcvvVZ547ty5PP3003h5eZGfn09wcDDz5t31FuFdfffdd0yePJnMzEzc3d1ZvHhxqbd96qmnSlWi2VAY5EhbVS2zlsrPhYQI8g//Sf7R1ZjnXiVDmhFlHIBR20F06jUMC8uSu2VKKYmPvsSOP2O5eiETB1cbAge649yqfhU2wrBUt4FXqnha2RnkSNsHoUba1mLafArit3E28ids4v/CVl4jA3NO24XQuMsobNv2KrHSZ4G2gOM7LrB7TTzXr+bg0qY+gQM9sG9S+/rwq4Rf85Ul4Rtk1wZ101ZBY4xRsxBcmoWANp+Tu9ZxMfJnPFO2YPvH32Suskbbsj/W/o/rav0Xq/JppDGiTefGtOjQkEPhZ9n7dwK/vr+bZv4OdAxzp66Dhf7apSh6pM7wlWol9sIVIv76lfrxq+kh9mAlssmzcMDE6zHwGqqb7P226/Y5WfnsX5/IgY1nKMiXtH6oMQGP1I4+/OoMv+ar9mf4ilISD8f6eIybTPK1sczfepxzu/8gNH0r3Xd8g/GOL5F2LRBew3TJv24TAMzqGNNpgAeeIc7sWZfA0a3nOLHjPN7dXfDtpfrwGxITEzUZTmUyyDN8ddNWKa3rOfn8sus0v209SPuMCEaYR+GpvTEIp2ln8B4ObQaA+c3ufKnJmexaHU/M7ouYWRrr+vB3c8LYpOb14a9uZ/hK2ambtkqtk6ctYFX0OeZHxJGRHMtoi508brod26zTupm9WvcHnxG6yV2MdF08U06ns2NlLKePXsGqnhkB/dxo1ckRozLM3mXoVMKv+VTCV2otKSXhJ1OYFx7LzvjLPGQez0sN9+GVuhGjnDTdpC4+w8FnJNTXDX9POnGVqD9iSU64Rj1HCzoN9MDNu2bU4TeEhG9lZcX166Wr4XP6tG7ynSZNmpS4zsyZM+9b46asxy2NVatWcfToUV599VVWrlxJixYtaNOmTYXt/0Gpa/hKrSWEoFtLB7q1dCD6TCrzIxox6LA7lpqBvOYez0A2YxExByJm6yZz9x2Nc5swhrziR9z+FHb8Gcdf8w7h6G5D4CAPGjevp+8m1SpZWVn6DqFEYWFhhIWFAbrSyf369TOIhF8WBvnZVQjRXwgxv7BQkqI8CB+Xunw10o9NL4YQ5ufOzLiWtI2ZyIymP3Pe7yVIOwN/TISPWyL+ehmPxskMf7MD3Ua1Iv1yNn98vJ81Xx7g8tmKrTBZW0kpmT59Ou3atcPT05Nly5YVvfbRRx/h6enJgAED+PjjjwH49ttvCQgIwNvbm0cffZTMzHuXxY6PjycwMJCAgAD++9//3vLa7Nmzi0oYv/XWW4BuIpPWrVszYcIE2rZtS2hoaNEbzty5c2nTpg1eXl48/vjjgG5ClWeeeYbIyEhWrVrF9OnT8fHxITY2lvbt2xcdKyYmBj8/wyz0a5Bn+Kq0glKR3OwseX+QJ88/3IIlkfH8EJXIzyfaE+jWnZe7XcInZRVi7xLYNR8jJ3/a+I2l+X8HcGj7Vfb9k8gvs3bRsoMjHfq7YWNXckXP6iAkJOSOZUOHDmXq1KlkZmbeUckSYOzYsYwdO5ZLly4xZMiQW14LDw8v9bF///13oqOjOXDgAJcuXSIgIIDg4GCio6NZuXIlO3fu5MyZM6SmpgIwePBgJkzQpYA33niDhQsXMm3atBL3/9xzzzFlyhTGjBnDl1/enGNp/fr1xMTEsGvXLqSUhIWFERERQZMmTYiJieHnn3/m22+/ZejQoaxYsYJRo0bx4YcfEh8fj5mZWVE8hYKCgggLC6Nfv35FPw9bW1uio6Px8fFh8eLFjB07ttQ/l6pkkGf4ilIZ7K3NmN6rFZEzevDGI61JuJLNoL+M6ZP0BGtDw9GGvg856bDqGUw+b0v7/M8Z9WwDfB9uwqm9ySyduYNtv8aQdT1X302plrZt28bw4cPRaDQ0bNiQrl27snv3bjZs2MC4ceOwsNANiKtbty4Ahw8fpkuXLnh6erJ06VKOHLl3Sa7t27cXTaxSfErF9evXs379enx9fWnfvj3Hjx+nsPefm5sbPj4+APj5+ZGQkADoyjmPHDmSH3/88Z4zahV66qmnWLx4MVqtlmXLljFixIiy/GiqjEGe4StKZbIyM+apLu6MCXRl1YFzzI+I5ek/EnCq25onO//CiMbnMD/wPez7DvPd3xLUJAjPERPYHdOSg5vPcDTyHO1Dm+LdwwUTs+rVlfNeZ+QWFhb3fN3Ozq5MZ/S3K6mDiJSy6Aa5mdnNwXBjx45l5cqVeHt7s2TJklId+2432qWUzJgxg0mTJt2yPCEh4ZbjaTSaoks6a9euJSIiglWrVvHuu+/e983m0Ucf5e2336Z79+74+fnRoEGD+8aqDwZ5hq+u4StVwdTYiCF+zvz9XDALn/CncV1z3ll7jE5LM/nE+iWuTj4IPd+F9PNYrx9H9+QBPB66D2f3OuxcFceP/43icMRZtKocc6kEBwezbNkytFotKSkpRERE0KFDB0JDQ1m0aBGZmZm4uroWlUBOT0+nUaNG5OXlFZVJvpfOnTvzyy+/ANyyfq9evVi0aFFRj52zZ8+SnJxc4n4KCgo4c+YM3bp146OPPiI1NfWO3j7W1ta3TLZubm5Or169mDJlCuPGjSv9D6WKGWTCl1KullJOtLW11XcoSi1gZCTo0bohv00OYsWUQAJc6zN3YwyBc6OZebkHSaO3wcgV4ORP/QOz6Hu1N4P91mJrk8eWn07wyzu7OLU3+Z71+xUYNGgQXl5eeHt70717dz766CMcHR3p3bs3YWFh+Pv74+Pjw5w5cwB499136dixIz179qRVq1b33f///vc/vvzySwICAih+shgaGsqIESMIDAzE09OTIUOG3JKsb6fVahk1ahSenp74+vry/PPPF11mKvT4448ze/ZsfH19iY2NBWDkyJEIIQgNDX2An07VUP3wFeUuTiWn882WOFZGn6VAQn+vRkwO8aCV2VXYvQD2fY/MSiWhzhCi0oZy9aoJDq42BA32wKmF4XTlNIR++GVReA3d1dVVr3E8iDlz5pCWlsa7775bpcdV/fAVpZyaOVgz+zFvXghtwcKt8fy06zQro8/RraU9U0KeJyDkVcSh5bjtnEfTrMc5Yd+PXRdHsPKTazT1bEDgQA8aOFnpuxnVTk5Ojr5DeCCDBg0iNjaWTZs26TuUe1Jn+IpSCqmZufwQlciSyAQuZ+Ti17Qek7t60KOlPUant8GOeeQf38DBrEfYlzmMnHxTWgU60qG/O9b1zfUWd3U7w1fVMsvOIM/whRDuwOuArZRyyP3WVxRDUtfClGk9mvNUF3d+23uG+RFxTPh+D80drJjc1YOwoT9ikhpP+53zaLP3Gfal9uHgjv7E7DqPZ7cm+PVxxdxSVYJU9KtcN22FEIuEEMlCiMO3Le8thDghhDglhHgVQEoZJ6V8sjzHUxR9q2OqYUygK+EvhfDZMB80RoIXfztAyOxwlhw3IuvhDzF/aTdBg9wZ2fQtmptsJnpDIj/O2MK+v+PIz9NWecyG/CleKZ+y/m7LdUlHCBEMXAe+l1K2u7FMA5wEegJJwG5guJTy6I3Xl5f2DF9d0lEMnZSSzSeS+WpzLHsSr1Lf0pRxQa6MCXTF1lTC4RVc2riMqMTOnM71w8oil46DWtKisytGRpVfnC0+Ph5ra2saNGhQLYrBlaZ4mqIjpeTy5cukp6fj5uZ2y2uVVi1TCOEKrCmW8AOBmVLKXjeez7gR3Ac3nquEr9RIuxOu8HV4LJuOJ2NpqmFUp6Y8+ZAbDtZmcGoDSet+JzLGl5T8ZjSwzSBoaDtc2rtWaiLOy8sjKSmJ7OzsSjuGoj/m5uY4OzvfMXFMVSb8IUBvKeVTN56PBjoCbwHvoTvzX1D4BnCX/U0EJgI0adLELzExsVzxKUpVO3ruGvO2xLLm4DmMNbrBXZODPWjSwAJ5Zg+n/viTHSdac03riLP9ZYKG+2LfRs3frFScqkz4jwG9bkv4HaSUJVc9unOfasYrpdpLvJzBNxFxLN+TRH5BAf29GzMlxINWjjZoz5/g8LK17DnpTnaBDS0an6bjyCBsPFroO2y9GjVqFAA//vijniOp3gzmkk5ZqEs6Sk1w8Vo2C7fFs3RHIhm5Wh5u3ZCnu3ng26QeORfi2ffjBg6cckEi8Goai9+oHpi71M5uiYXVPMtTs0ep2oRvjO6mbQ/gLLqbtiOklPeuPnTrPtUZvlLjpGbm8l1kIosj40nNzCPIowHPdGtGoEcDriclseuHLRw/7YiZyMC/2Uk8R/RH0+j+JQVqEpXwK0ZJCb+83TJ/BqKAlkKIJCHEk1LKfOAZ4B/gGPBrWZK9otRUdS1Mee7h5mx/pTuv9W1FTPJ1RizYyeCvI9mZbkr3GSMZ9kILHOzz2B7jx0+z9hHzxdvI5OP6Dl2pIdRIW0XRk+w8Lcv3JjFvSyxJV7No3ciGp7t50KddI85GJxK57CCX06xwMDlJZ+/TNB44Huxq9s1ddYZfMfQ+0rYsil3S0XcoilJpzE10XTeHBbiwKvocX4Wf4pmf9uNud5Kp3ZoxaNYjxG+PZcefOfyxpwVuh34jqEMqdftOhXpN9R1+pQgMDNR3CDWaOsNXFAOhLZD8ffgCX2w+xbHz13CpX4cpXZsxwNORY+tPsm/DObT50M5yPf6djanT8z9g7ajvsBUDVGk3bSuDummr1GZSSjYeS+bzTTEcSEqjka05k7t60L9lQw6uOsrRnamYiEz8bVbi1c0FTdfnoI7hlGRW9K9aJfxC6gxfqc2klGyNucTnm2LYnXAVe2szJgW707eJHft+P0LiiSxsNBcIrLcCj9DOiE6TwdRC32GXy6OPPgrAihUr9BxJ9VatruEriqKbnzW4hT1dmtuxI+4KczfGMGvtMeZZmTKhizuh3W3Yu0LDP8lP4/jbMTpvHYxj7+HgMxI01fNP+/Lly/oOoUYzyP8V6qatotwkhCDQowGBHg3YnaBL/B/8dZxvLE15KsSNICMz9q+GFWda0/yHCAIj+mLd9wVo0QuqQcE0peqoSzqKUg3tTbzK/zbGEHEyhXoWJkwIdMPzGhzZdBoK8vGxWEn71ucx7fMmOLXXd7ilprplVgx1SUdRahC/pvX4fnwH9p/WJf6PNp7UJf6+TXFNymbv/sc4tj+VTiffpVVHe8TDb4Kts77DVvTMIBO+uqSjKKXj26QeS8bpEv9nG2L4aOsp6luaMjHUCcvjFmw6PY1D4XE8FD2Cxt1CofN/wMxw59rt0aOHvkOo0dQlHUWpQfadvsqn/55ka8wl7CxNmejqiOmhK2Rc09LMfBuBjn9j03saeA0Do3JVVlEMmOqWqSi1yJ6EK3zy70kiYy/T2MqM8Q0akH/kKmjz8LX4g/bN4jDp9x64BOg7VKUSqISvKLVQZOwlPll/kj2JV2luVYdhplbkxF3HyvgKgZaLad7BGRH6tsGM2O3Tpw8Af/31l54jqd4qpVqmoiiGLcjDjt8mB7JkXADmtqbMupJCuLMg17Yh/6a9yB9bvEn5eDBEfg7aPH2HS1ZWFllZWfoOo8YyyIQvhOgvhJiflpam71AUpdoTQhDS0oFVz3Rm3ig/0iyNeD//GocaaUgxas2vF94l/Lcksj7vCXHh+g5XqUTqko6i1DLaAsmf0Wf5dMNJki9lMdjMCpdL+ZiKTDpa/Uhbf0uMer8HNo2qPDbVD79iqH74iqIAoDESDG7vTD+vxvyy+zRzN55CWuUyVFgTcW0SR7fGE3xkOI16PwYdJlXbMg3KnWrkbzI7P5sCWaDvMBQ9ESWUExCIu64jbvy78aTo+6Kv4sYaNaxMgamxEWMCXXm0vTOLtsXzzZY4nLUFPJLnyu8X36Tlss0E7u6H5eBZ4HzHyWKl6NevX5Ucp7aqkZd0Jm+YzPaz2yshIkXRvREYCSPdG4K4+bxwmZEwQgjdV43Q6F7DCCMj3XNjI2M0QoPGSIOxMMbEyARjoxtfNbqvZhozzDRmmGpMMdeYY6Yxw9zYHAtjCyxMLLAwtsDK1AorEyusTK2wNbXF1swWc2PzB27XlYxcvth0ip+jEuiUZYx/tgZTkU1Hy6W069IIo9A3wdy2An+SSmXRe7dMIYQl8BWQC4RLKZfeb5sHTfj/Jv7L2fSzZQ9SMSiSsv/fLGmb4v/Pb1+n8DWJvOX7oq/yxms3Xr/X1wJZQIEsQCLRFmgpkAVopbbouVbeeBRoyS/IJ1/m677eeOQV5JGrzSW3IJdcbS7Z+dnkaHPI1maTX5B/3/aba8ypZ16PBuYNaFCnAfYW9jS0aEhDi4Y0tmqMs7UzDS0aYmxU8of7M1cymbP+BBF7z9Mv1xjHHGMaGMcT0vA3HB+dAq3DVFE2A1cpCV8IsQjoByRLKdsVW94b+B+gARZIKT8UQowGUqWUq4UQy6SUw+63f3XTVlFuyivIIys/i6y8LDLyM8jIzSA9L53ruddJy00jLUf3uJJ9hctZl7mcfZnkzGSuZF+5ZT/Gwhhna2c86nrgbutOi3otaN2gNS7WLhiJmx33Dial8v7ao1w5kUafHA2mWmPa1FlPoM95zAe+D7ZOFd5GddO2YlTWTdslwBfA98UOpAG+BHoCScBuIcQqwBk4dGM1bTmPqyi1jomRCSamJtiY2pRpu1xtLhczL3Lu+jmS0pM4k36GxGuJxKbFEn4mHK3U/TlamljStkFbfBx8aO/QHm8Hb36eGMjGY8l8tPYoTkk5kPUwcTszCDrxIq0GPozwH69KNFQj5Ur4UsoIIYTrbYs7AKeklHEAQohfgAHokr8zEI2B9v9XlJrIVGOKi7ULLtYudGzU8ZbX8rR5xKbFcuzyMY5cPsKhS4dYcGgBBbIAjdDgbe9NUOMgZj/Rmf0xlnz3dwwhV0zYdGUix384Qtfdo6k/7G2wU4UOq4NyX8O/kfDXFF7SEUIMAXpLKZ+68Xw00BF4Bd2ngWxgW0nX8IUQE4GJAE2aNPFLTEwsV3yKopRNZl4mB1IOsOvCLiLPRXL08lEAGlk2okvjbly60JLEbXXokWWEsRT4Wq/Cv38LjDtPLXcXTnVJp2JUZT/8u93NkVLKDGDc/TaWUs4XQpwH+puamvpVeHSKotyThYkFgY0DCWwcyHPtn+NK9hUikiLYmLiRlbHLyS3IxTHQifVXA2l2xIuC9MGcWnaOkN0TcB79Mji01ncTlBJURsJPAlyKPXcGzlXCcUok8/JAn91Nq6IHQ2mOUZ517rFtTeuPrtxbffP6DGw2kIHNBnI99zqbzmxi1alV7MxcTmL7FcSk9CQkPpg/Y0bT6sOldH7EAfOQBzvbHzp0aCW0QClUGZd0jIGTQA/gLLAbGCGlPFLWfT9oL53TEyaSsXVrmbdTyqn4G0EZvhfFlwtx1+9F8ee3rXPX18SNgVKFDyMjXZ95YXTj+2KvGxmBkbj5mpGuHz0azc3lGk3RMqExAiMNaIwQGuOi50KjAWONbpmx5sa6xghjY91yYxOEsTHCxARhUvi12MPU9MbDDGFmipGZGcLMDGFqhpG5GaJOHYzMzXXLDOBG6bnr5/g95nd+O/kbaRnp+J8eiM+FQCyM0nmo6Uaaj3sW4dBS32HWSpXVLfNnIASwAy4Cb0kpFwoh+gKfoeuWuUhK+V4Z91s449WEmJiYMseVtnYteWeSyrxdxaiCTxal+J2V6vda0jr32rT4NrdsX6yfe0nr3LK6LPb15vdF28rb1rmxnpTy5n4KXy8oKPaavLmtLLi5rEAW7UcWFBT7Xntz/QKtbv0CCVotUhbovi/QIguXFRRAQYFuO+2Nr/lapFare12rRebn3/p9Xt4tyyuCMDfHqE4djCwsMLKog5GFJUaWNx5WVhhZW6GxssbIxhqNjS0aWxs0devqHvXqobG11b1BVYBcbS7/JPzDokPfcfXcdbqdGoFdpjMu5nvo0b8ult2eKnVPnszMTAAsLCwqJLbaSu8Drx6E6oev1DSyoODmm0BeHjI/X/c1NxeZm0tBTm7R9zI3h4LsbGROLjInm4Ks7KKvBVlZFGRlIjMzKcjMpCAjA21GBgUZGRRcz6AgPZ2CjIySAxFCl/jr18O4gR3GDg4Y29tj0tABY8dGmDg2xKRxYzR2dqW+hCelZPu57XyyYx4WJ+sTcLovJiKPQJcIfCY9jajvet99qJu2FaNaJfzynuErigJSq6Xg+nW0166hTbuGNi0V7dVUtKmpaK9cIf/KZbSXL5Ofcon8lBTyU1KQubm37EOYm2Pi7ISpswumrq66h5sbZi2aY1yvXonH3pYUxUf/zMfzeCCN05thUecwgwc5YNtl2D3vD6mEXzGqVcIvpM7wFaXqSCnRpqaSf+ECeecvkHfuHHlJSeSdTSI38TS5iYnInJyi9TV2dpg1b4Z5mzbUadsW87ZtMWnSpOgTgZSSXw7/xV9rN9H+dA+MkLR0jST06ecRlg3uGoNK+BWjWiV8dYavKIZHFhSQf+ECObFx5Jw6RU5MDDnHj5MTE6PrGQdo6tenjq8vFr4+WHTsiHmbNuRJLTNXL0BEGdH4Wgu0lscYMsydxh363HEMlfArRrVK+IXUGb6iGD6Zm0vOqVNkHTpM1v79ZO3fT+6NAZNGNjZYdAjAqmtXrni15rOVK/FI9ANRQMNme3n8mZcRJjcrfKqEXzFUwlcUpcrkp6SQsXMXGTuiyIiMJP/ceQDMPT053KIlBzKaUC/bg+s2hxg2tiOubToBsGTJEgDGjh2rp8hrhmqV8NUlHUWpOaSU5Jw8yfXNm0nfuInsQ4eQCI607c0F+17kGWdj5xXDmAkvqrLLFaRaJfxC6gxfUWqe3NOnubbuL9LWrObK+UyOtB5LhpUz2Za76DOuLw6OjbGzs9N3mNWaSviKohgUKSVZ0dFc/HkZB49CktPD/G/Vs5hZ5HAwJsYgRhNXVyUlfIP8iQoh+gsh5qelpek7FEVRKokQAgtfX9w++pB+P71GE+uNCCnIybLi9yEvc+nv9aUbMa6UmkEmfCnlainlRFtbNX+motQGGhsb+s/5gHpNLZHiOhcc+vLP4jiiBw7n+tZt+g6vxjDIhK8oSu1kXscCp2bOaFts50q9puy0H8G+V+dyeuIkcuLi9B1etacSvqIohkUInn3hv3g+lsEV6ysc8pzMngsOnBowmIsffHjvGkHKPVVGPXxFUZQHMmXKlKLvQ3oMwcsrkc8//x7oxkV7Dzr8upD0f//FceZbWAUH6y/Qasoge+mofviKohSSWi1fznuT3KOBGBdAq9OrcIndik3//ji+8Toada/vDqpbpqIoBu/MmTMAuLi43PHaxo0LiVyroX5mE0zz9xC0cymmDerT+MMPsezU8Y71a7Nq1S1TUZTaafTo0YwePfqur/Xo8STjnvckyWEbucb+rO/6POlmdTk9bhwXZ88uKuCmlEwlfEVRqg1nFz9mvvYcKc3/QBQ0ZEuzJ0joMoQrCxeR+MRY8i5e1HeIBk0lfEVRqhUzc1veev5/WHWJINMkm1NGXYjo/QxZx44RP/hRMnbs0HeIBqvKEr4Qwl0IsVAIsbyqjqkoSs0kjIx4avgHPPR4JmfrHSE/uzXruj5Dvk19Tj/5FFd++FGN0r2LUiV8IcQiIUSyEOLwbct7CyFOCCFOCSFevdc+pJRxUsonyxOsoihKcV07jWfyM0GcdFqHebYrfzcfTnpAdy6+9x4XZr6truvfprT98JcAXwDfFy4QQmiAL4GeQBKwWwixCtAAH9y2/XgpZXK5o1UUpUZ78cUXy7yNq1MH3nnBlbe+nUHTmEeJsnyYZqFOuC9bQm5CAs6fz0VjY1MJ0VY/pe6WKYRwBdZIKdvdeB4IzJRS9rrxfAaAlPL2ZH/7fpZLKYfc4/WJwESAJk2a+CXemDlHURTlXvK0uXz409MY7e2BTU59LC3iCfznC8zc3HD5dj4mjo76DrHKVEa3TCfgTLHnSTeWlRRAAyHEPMC38M3hbqSU86WU/lJKf3t7+3KEpyhKdXPixAlOnDjxQNuaaEx5Y9R8nPse55zNKbIymrO+x/PknDtPwuPDyVGDOMuV8O82NU2JHxeklJellJOllB6l+BSgyiMrSi00adIkJk2a9MDbCyF4InQmoSNsONowHOMcV/4JfJZsrREJI0eRuX9/BUZb/ZQn4ScBxYfDOQPnyheOoihK+XXzHsW0iX3Z5/obmnxHNnlNJL2uM6effKpWd9ssT8LfDTQXQrgJIUyBx4FVFRGUqoevKEp5tXLqxHtTZxDd5jsKsCKq2RhSGnlyZuIk0jdv1nd4elHabpk/A1FASyFEkhDiSSllPvAM8A9wDPhVSnmkIoJSl3QURakIjjZNmDt5HnHtl5JhksuBJsM47daFpGnPcm39en2HV+VU8TRFUQxGSEgIAOHh4RW63/yCfN77/WkKogKxz3DGJXkDLU6uxenTT7Dp2bNCj2UISuqlY5D18IuVR9Z3KIqiVKE33nijUvZrbGTMm4/OY0G9tzm54RqCUHI0lsjnX8C5hib9u1Fn+Iqi1Cpr9i/i7z+SaJ38EHVTo/E58h1N/vcp1t276zu0ClOtyiOra/iKUjtFR0cTHR1dqcfo5zue0U904oDzWlLr+rDHewoJz79MRmRkpR7XEKgzfEVRDEZlXcO/m5jkg3yw5Et844ZglZFE+2MLaf7N/7Dw86v0Y1c2dYavKIpSTHMHLz6Y8hYHWy0l3aoxe9pO5fjU6WQdOnz/jaspg0z4qh++oihVoaF1Yz6ZMpcYz1+4ZmnL3jZTOTL1FXLi4/UdWqUwyISvKIpSVazNbPhk0nwu+q/imoUx+1pMYP+El2vk7Fkq4SuKUuuZaEx4f9zX5Adv45q5lmjXJ9g5/mW0qan6Dq1CGeRN22L98CfEqAp3ilJrRN7oKRMUFKS3GBaseY/UvzywzLWl1ZU/6br0M4zMzfUWz4Mo6aatQSb8QqqXjqIo+vBH+DzifrXCIs8Oj6x/6PndHIRGo++wSq1a9dJRFKV2ioyMLDrL16dBIZNpP86cHM1FTln0Ye3k12vEHLkq4SuKYjBee+01XnvtNX2HAUC3gCEET3OhQCZxWvRg+XNv6jukclMJX1EUpQTt23an58s+GOUnkJITzPevztR3SOVikAlfDbxSFMVQtGzenr6vB2GWFU/61Yf4+u239B3SAzPIhK8GXimKYkiaeLTlkdcfwup6LAXnujD7/deq5TV9g0z4iqIohsaxZRv6vNgBm/RTWCR2Z+bHL5NfkK/vsMpEJXxFUQzGZ599xmeffabvMErk0N6Xh8e3wuZaHA4xobzyvxfJys/Sd1ilphK+oigGw8fHBx8fH32HcU+NHg4hpI81NtcScTvej+e/+A9pOdXjfmOVJXwhxEAhxLdCiD+FEKFVdVxFUaqPDRs2sGHDBn2HcV8uY4bTue0VbNNP0+74YJ6d/xwXMi7oO6z7Ku0k5ouEEMlCiMO3Le8thDghhDglhHj1XvuQUq6UUk4AxgLDHjhiRVFqrFmzZjFr1ix9h1Eqbq+9RIBVNNbpF/A/8jjPLnqe2NRYfYd1T6U9w18C9C6+QAihAb4E+gBtgOFCiDZCCE8hxJrbHg7FNn3jxnaKoijVljAyosWnH9L+2hqsMi/z0NExPP/jdKKTo/UdWolKlfCllBHAldsWdwBOSSnjpJS5wC/AACnlISllv9seyULn/4C/pJT7SjqWEGKiEGKPEGJPSkrKg7ZLURSl0hlZWtLy68/wjf0Oy6w0uh8dz8vLX2fLmS36Du2uynMN3wk4U+x50o1lJZkGPAwMEUJMLmklKeV84G1gn6mpaTnCUxRFqXymzk60+N/7tD/0FZY5OfQ5OoE3177LylMr9R3aHcqT8MVdlpU4EkFKOVdK6SelnCylnHevHauBV4qiVCcW/v64vjqN9vvmYpEjGHhkMv+36RMWHlpoUAO0jMuxbRLgUuy5M3CufOHoFKuHXxG7UxSlmvjmm2/0HcIDq/fYYzgdO472z8/Z6/88ww5P4RvxOZezL/OS/0sYCf33gi9PBLuB5kIINyGEKfA4sKpiwlIUpTZq2bIlLVu21HcYD6zhjFdp2KYxPoe+wTS3HqOOTGDZwd+YsXUGedo8fYdX6m6ZPwNRQEshRJIQ4kkpZT7wDPAPcAz4VUp5pCKCUpd0FKV2Wr16NatXr9Z3GA9MmJjg9L/PsDfPoE3cL2iynXny8Dj+OfUv0zZNIzMvU7/xGdL1pUJqikNFqZ1CQkIACA8P12sc5ZV97BgJw0dwqllPEu1CMbOM5gvPH2lj34Yve3xJPfN6lXr8ajXjlTrDVxSlOjNv3ZpGs2bhcXg1dXMPkJPRnpeOhHHy8nHG/DWGc9cr5HZnmRlkwlcURanubPs9Qv3x42kfuQA0Z0hLD+H1o124nJnM6HWjibla9VcvDDLhqwlQFEWpCRxeeB7LoCC6RXzMdbNUzqUO4O2jbaAgjyf+foJ9F0scg1opDDLhq0s6iqLUBMLYmMYfz8HU3p6w6M9JNcvk5OWRvH/ckQYm1kz8dyLhZ8KrLB6DTPjqDF9RaqcffviBH374Qd9hVCjjevVw/uJzROpVBiX/RppJLvsujOf/TpjT3KYp/9n8H/6I+aNKYjHIhK/O8BWldnJxccHFxeX+K1Yz5m3a0OidtxH7d9HbKprrGsn2pCeZfTKbjg7teTPyzSoZlWuQCV9RlNpp2bJlLFu2TN9hVArbAQOoN3Ikdf74Ad9W6VzHjI3xo/m/uEv0aRrKZ/s+Y/ae2RTIgkqLQSV8RVEMxtdff83XX3+t7zAqTcNXXqaOry/Oi/8PB38zrmkbsP5oP95NOMvIlsP54egPvL7tdfIKKmdUrkEmfHUNX1GUmkiYmuL02WcYWVjQ6bf/I8vbiuR8Dzbs92X6mQSe9ZnGmrg1lTYq1yATvrqGryhKTWXS0AHnTz8h98wZhh/5iTNu5iTkBLAt0panLpxhZuBMDqUcqpTBWQaZ8BVFUWoyi4AAHKa/xPUNG/iP7RHiGggOZ/XhwIYEHr2SzLrB62hWr+KrBauEryiKogf1n3gC6z69SZ37P57vruF0nQK2p48lduVKbGM2VsoxDbJ4WiF/f3+5Z88efYehKEoVuXTpEgB2dnZ6jqRqaK9nkDB0KNq0NLLnLmTF4tM4arUMspuJ45R50NjngfZbrYqnqZu2ilI72dnZ1ZpkD6CxssT587kUZGVhM3smgSObkyY0rLzyJmka9wo/nkEmfHXTVlFqpyVLlrBkyRJ9h1GlzDw8aDzrXbKio/HZ/BP1+7oQK0w4eSWnwo9VnikOFUVRKlRhsh87dqxe46hqNn37krk/mqvf/8CwT3xIeTOYZg5WFX4cgzzDVxRFqW0aTn+JOj4+XHjjv7hcT66UY1RZwhdCtBZCzBNCLBdCTKmq4yqKolQHukFZnyLMzDj73HMUZOpp4JUQYpEQIlkIcfi25b2FECeEEKeEEK/eax9SymNSysnAUOCOu8eKoii1nYmjI43nzEabfp28s2crfP+lPcNfAvQuvkAIoQG+BPoAbYDhQog2QghPIcSa2x4ON7YJA7YBldPJVFEUpZqz6twZj3/+xqx58wrfd6lu2kopI4QQrrct7gCcklLGAQghfgEGSCk/APqVsJ9VwCohxFrgpweOWlGUGmndunX6DsEgGJmZVcp+y9NLxwk4U+x5EtCxpJWFECHAYMAMKPG3KoSYCEwEaNKkSTnCUxSlurGwsNB3CDVaeRK+uMuyEoftSinDgfD77VRKOV8IcR7ob2pq6vfA0SmKUu189dVXAEydOlXPkdRM5emlkwQUn5rGGaiQ8m5q4JWi1E6//vorv/76q77DqLHKk/B3A82FEG5CCFPgcWBVRQSlSisoiqJUvNJ2y/wZiAJaCiGShBBPSinzgWeAf4BjwK9SyiOVF6qiKIpSHqXtpTO8hOXruMcN2AclpVwNrPb3959Q0ftWFEWprQyytIK6pKMoilLxDLoevhAiBUgEbIHi2b/485K+twMuVUAYtx+7POuW9Pq92ne/59WxzaVZVl3bXNrf8d2WqTbXvjaX9Fp529tUSml/x1IppcE/gPklPb/H93sq49jlWbek1+/VvjK2v1q0uTTLqmubS/s7Vm1Wbb7XaxXV3tsfBnlJ5y5W3+N5Sd9X1rHLs25Jr9+rffd7Xh3bXJpl1bXNpf0d322ZanPta/P9fh4VyqAv6ZSHEGKPvMsUXzWZanPtoNpc81VWe6vLGf6DmK/vAPRAtbl2UG2u+SqlvTX2DF9RFEW5VU0+w1cURVGKUQlfURSlllAJX1EUpZaolQm/Ns6vK4QYKIT4VgjxpxAiVN/xVAUhhLsQYqEQYrm+Y6ksQghLIcR3N363I/UdT1WoDb/X21XY329ldO6vzAewCEgGDt+2vDdwAjgFvFrKfRkBC/Xdpipuc71a2Obl+m5PZbUdGA30v/H9Mn3HXpW/7+r2e62gNpfr71fvjX6AH1Iw0L74DwnQALGAO2AKHEA3z64nsOa2h8ONbcKASGCEvttUVW2+sd3HQHt9t6mK21ytEkMZ2z4D8Lmxzk/6jr0q2lxdf68V1OZy/f2WZ8YrvZC1cH7dimizEEIAHwJ/SSn3VXLI5VZRv+fqqCxtRzcRkTMQTTW+RFvGNh+t4vAqRVnaLIQ4RgX8/Vbb/yC3udv8uk4lrSyECBFCzBVCfEMllHeuImVqMzANeBgYIoSYXJmBVaKy/p4bCCHmAb5CiBmVHVwlK6ntvwOPCiG+ppKH5evBXdtcw36vtyvp91whf7/V7gy/BJUyv66BK2ub5wJzKy+cKlHWNl8Gquub2+3u2nYpZQYwrqqDqSIltbkm/V5vV1KbK+Tvt6ac4Vfa/LoGTLW5drS5UG1su2pzBbe5piT8Sptf14CpNteONheqjW1Xba7oNuv7TvUD3Nn+GTgP5KF7N3zyxvK+wEl0d7hf13ecqs2qzartqs2G1mZVPE1RFKWWqCmXdBRFUZT7UAlfURSlllAJX1EUpZZQCV9RFKWWUAlfURSlllAJX1EUpZZQCV9RFKWWUAlfURSlllAJX1EUpZb4fzPoej9uMnlNAAAAAElFTkSuQmCC\n",
"text/plain": [
"