\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.8?urlpath=lab/tree/source_population_model.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",
"[source_population_model.ipynb](../_static/notebooks/source_population_model.ipynb) |\n",
"[source_population_model.py](../_static/notebooks/source_population_model.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Astrophysical source population modeling with Gammapy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction\n",
"\n",
"The [gammapy.astro.population](..\/astro/population/index.rst) package contains some simple Galactic source population models.\n",
"\n",
"Here we provide some Python code to compute observable parameter distributions for Galactic gamma-ray source populations.\n",
"\n",
"* Observables: Flux, GLON, GLAT\n",
"* Source classes: Pulsar (PSR), Supernova remnant (SNR), pulsar wind nebula (PWN)\n",
"\n",
"References:\n",
"\n",
"* Section 6.2 in the Fermi-LAT collaboration paper [\"The First Fermi-LAT Catalog of Sources Above 10 GeV\"](http://adsabs.harvard.edu/abs/2013arXiv1306.6772T)\n",
"* Axel Donath's bachelor thesis [\"Modelling Galactic gamma-ray source populations\"](http://pubman.mpdl.mpg.de/pubman/item/escidoc:912132:1/component/escidoc:912131/BScThesis_ddonath.pdf), specifically Chapter 4.\n",
"* Casanova & Dingus (2008), [\"Constraints on the TeV source population and its contribution to the galactic diffuse TeV emission\"](http://adsabs.harvard.edu/abs/2008APh....29...63C)\n",
"* Strong (2007), [\"Source population synthesis and the Galactic diffuse gamma-ray emission\"](http://adsabs.harvard.edu/abs/2007Ap%26SS.309...35S)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import astropy.units as u\n",
"from gammapy.utils.random import sample_powerlaw\n",
"from gammapy.astro import population"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulate positions"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Spatial distribution using Lorimer (2006) model\n",
"n_sources = int(1e5)\n",
"\n",
"table = population.make_base_catalog_galactic(\n",
" n_sources=n_sources,\n",
" rad_dis=\"L06\",\n",
" vel_dis=\"F06B\",\n",
" max_age=1e6 * u.yr,\n",
" spiralarms=True,\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulate luminosities"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Several source population models, e.g. the 1FHL paper or Strong (2007), use power-law luminosity functions.\n",
"\n",
"Here we implement the \"reference model\" from the 1FHL catalog paper section 6.2."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Source luminosity (ph s^-1)\n",
"\n",
"luminosity = sample_powerlaw(x_min=1e34, x_max=1e37, gamma=1.5, size=n_sources)\n",
"table[\"luminosity\"] = luminosity"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute observable parameters"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"
\n",
" name dtype unit description \n",
"---------- ------- --------- --------------------------------------\n",
" age float64 yr Age of the source\n",
" n_ISM float64 1 / cm3 Interstellar medium density\n",
" spiralarm str18 Which spiralarm?\n",
" x_birth float64 kpc Galactocentric x coordinate at birth\n",
" y_birth float64 kpc Galactocentric y coordinate at birth\n",
" z_birth float64 kpc Galactocentric z coordinate at birth\n",
" x float64 kpc Galactocentric x coordinate\n",
" y float64 kpc Galactocentric y coordinate\n",
" z float64 kpc Galactocentric z coordinate\n",
" vx float64 km / s Galactocentric velocity in x direction\n",
" vy float64 km / s Galactocentric velocity in y direction\n",
" vz float64 km / s Galactocentric velocity in z direction\n",
" v_abs float64 km / s Galactocentric velocity (absolute)\n",
"luminosity float64 \n",
" distance float64 pc Distance observer to source center\n",
" GLON float64 deg Galactic longitude\n",
" GLAT float64 deg Galactic latitude\n",
" VGLON float64 deg / Myr Velocity in Galactic longitude\n",
" VGLAT float64 deg / Myr Velocity in Galactic latitude\n",
" RA float64 deg Right ascension\n",
" DEC float64 deg Declination\n",
" flux float64 1 / kpc2 Source flux\n"
]
}
],
"source": [
"table = population.add_observed_parameters(table)\n",
"table.info()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Check output"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The simulation is done, you could save the simulated catalog to a file.\n",
"\n",
"Here we just plot a few distributions to check if the results look OK."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnX+QXNV157+ne96IHjlLS0F2pDaDFIqIWBHSBMWootqUURzkWAYGCJZZSHkrqVVcFdcuCjW1Q8xGI4yDahUCW/m5ZOOKUyZY2MBYoGwk22iXjVKyI2VGiImlMhgQtFRBiTTYaFpST8/dP7rv6PXre9+77/d73edTpdLMm+737vt17rnnJwkhwDAMw3Q/hbQHwDAMwyQDC3yGYZgegQU+wzBMj8ACn2EYpkdggc8wDNMjsMBnGIbpEVjgMwzD9Ags8BmGYXoEFvgMwzA9Ql/aA7Bz1VVXieXLl6c9DIZhmFxx5MiRfxVCLPH6XKYE/vLly3H48OG0h8EwDJMriOgtk8+FNukQ0dVEdICIvk9EU0T0X1rbFxPRt4joB63/F4U9FsMwDBOcKGz4swAeEEL8LID1AH6biD4CYBTAd4QQ1wH4Tut3hmEYJiVCC3whxGkhxD+1fv4xgO8DqAC4HcBXWh/7CoDhsMdiGIZhghNplA4RLQcwBOC7AD4khDgNNCcFAB+M8lgMwzCMPyIT+ET0AQDPArhfCPEjH9/bSkSHiejwmTNnohoOwzAM4yCSKB0istAU9k8JIZ5rbf4XIloqhDhNREsBvKv6rhDiSQBPAsC6deu4GwvDJMD4RBW79p3AqekalpVLGNm0EsNDlbSHxcRMFFE6BOAvAXxfCPGHtj/tAfDZ1s+fBfDNsMdiGCY84xNVPPjcMVSnaxAAqtM1PPjcMYxPVNMeGhMzUZh0NgD4dQAbiWiy9e+TAHYC+BUi+gGAX2n9zjBMyuzadwK1eqNtW63ewK59J1IaEZMUoU06Qoi/B0CaP/9y2P0zDBMtp6ZrvrYz3QPX0mGYHmNZueRrO9M9sMBnmB5jZNNKlKxi27aSVcTIppXzv49PVLFh50tYMboXG3a+xPb9LiFTtXQYhokfGY2ji9KRTl1p55dOXft3mXzCAp9hepDhoYpWeLs5dVng5xs26TAM0wY7dbsXFvgMw7TBTt3uhQU+wzBtmDh1mXzCNnyGYdrwcuoy+YUFPsMwHbg5dZn8wiYdhmGYHoEFPsMwTI/AAp9hGKZHYBs+07VwzXeGaYcFPtOVcHkAhumETTpMV8I13xmmExb4TFfC5QEYphM26TCpE4etfVm5hKpCuHN5AKaXYQ2fSZW4+qtyeQCG6SQSgU9EXyaid4noVdu2MSKqOvrcMkwbcdnah4cqePTO1aiUSyAAlXIJj965mh22TE8TlUnnrwD8MYC/dmx/XAjxBxEdg+lC4rS1c3kAhmknEg1fCPEygLNR7IvpLbgUL8MkR9w2/M8T0Sstk8+imI/FBCDt3qVsa88PaT8rTHjijNL5MwBfBCBa/z8G4DecHyKirQC2AsDg4GCMw8kHSWaHZiE5iUvx5oMsPCtMeEgIEc2OiJYDeFEI8XN+/mZn3bp14vDhw5GMJ484Xyqgqe3G5WzcsPMlZehipVzCwdGNkR+PyS/8rGQbIjoihFjn9bnYTDpEtNT26x0AXtV9lmmSdHYoJycxpvCz0h1EYtIhoqcBfAzAVUT0DoDtAD5GRGvRNOm8CeC3ojhW3nEz2ST9UnFyEmMKPyvdQVRROvcIIZYKISwhxIeFEH8phPh1IcRqIcQNQojbhBCnozhWnvFKMko6YoUdpowp/Kx0B5xpmyBeJpukXypOTmJM4WelO+BaOgniZbJJI2KFk5MYU/hZyT8s8BPExA7KLxXDMHHBJp0EYTsowzBpwhp+gnCSEcMwacICP2HYZMMwTFqwwGcYJjK4cXy2YYHPMEwkcL2d7MNOW4ZhIoEbx2cfFvgMw0QC19vJPizwGYaJBG5mk31Y4DMMEwlueSbcPCUbsNOWYRjfuEXj7Np3AtXpGopEqNUbGNszhfOXZlFvNHtvsDM3PVjgMwzjC5NoHPvfp2v1jn1IZy4L/GRhgc8wCZP3WHW3aJzhoYry7yrYmZs8LPCZVJHCT5oAGkKgkrAQjEsAq/YLIPex6l7ROKaCnJ25ycMCn0kNp2mgIeK38TqF8M3XL8GzR6odAvjwW2dx4PiZwBPRQ+PH8NShk5Ado+V+F/QVXLXjPOBV9VX3dzt+iwbmfVWUFSJrYh4Fvd7EPG+EfQl1jbElRSLMCeF737pxqZrEEwDVG6DbbtJUfnyiim27J5Xf10EA3ti52cc30kN1He3XRfV3q0D4wBV9mJ6pB7qfbsdjzJuYR9XT9ssAPgXgXSHEz7W2LQawG8ByNHvafloIcS6K4zHpE0UavdfSP4jG7zYulW1ZJ5R120208V37TvgS9gBQIML4RDUXAsyr6mvUVWG9fAaMOVGZdP4KwB8D+GvbtlEA3xFC7CSi0dbv/zWi4zEpY/ISeq0ArixZyggOFaYvuNu4onISeu0nyHEaQrSZkrwEZdomDq+qr1FWheUM3uiIROALIV4mouWOzbcD+Fjr568A+D9ggd8VjE9UtaYY+RK6adpAUzCbCnvnvoN85tR0zXiC0ZlzJF7ORhMbtopavaG0+wPtKxvdtTWdLPKGSac4xow4M20/JIQ4DQCt/z8Y47GYhJDCRod8CXWa9tieKTz43LFAAtHkBXdL7ydSf2dhf7GtOfe96wdhFdUftork6WxUZZya4pxoavUG7t892Zadqru2Tx06iep0DQKXJ4FuyGi9+folvrYzelKP0iGirQC2AsDg4GDKo2G8cIuxtkde6DRtEw17wCpgpj7Xts1E0AJNYaty8I1sWoltuyeV35m51MDUwxvbtr149LRyrAv7+zy1ZpUNe+bSLM7N+FvR2LFr+7prq5osusHOvfeV08rtB46fSXgk+SdOgf8vRLRUCHGaiJYCeFf1ISHEkwCeBJpROjGOh4kAN7OKPWoiqFnDKhIuzM51bG80zB4NN4ehjPd3oloVvKeZmHTbVeNwM8Oo8DIlSQFeHrCMJ4+827nHJ6rac837uaVBnCadPQA+2/r5swC+GeOxmITQmUwq5VKbgFOZNayCxqZi28fC/j7MKaTeHGBcV314qIKDoxvxxs7NODi6cX5cfprIR135cXiogkfvXN1mOrpv/WCHKcmLU9M1+Imkzrud2+2e5/3c0iCqsMyn0XTQXkVE7wDYDmAngGeI6DcBnARwdxTHYtLFzWRix49ZgwA8vmUthocqWDG6V3vssBqdn3BB0/P0e3wv88reV067au/LyiXj60BAoPGmHQFkx+1cw9yLXiWqKJ17NH/65Sj2z2QHP0LTKeB0wlzY9utmCopCozMNF4w6ltyU7beu0pp+5ISjM005sV9XU7LWplD3PJRLVu59E2mQutOWyR9OYSiX3V4voO7lrdgE+cimlRj5xtH5UroSq2DmtI2SKGPJ/RwTaC8xrCrr4OUPAIBFA5bv44/tmcpUkpNupTV22yrX72VplZIlWOAzvgmqBZqYSeT3d7wwNW/aKJcsjN22qmdeWJOkJuDyhAuC0q7vt2rK+ERVG0WVloM0yEora6uULMG1dBjf6GrgVMolHBzdqPjGZVjzip4Vo3u19YD81OcZeni/1n+waMDCxO/dEmyACRPm+cwridbSYXqLMKnuaZhJup0oMlHdwh8B4P0Ls7mp9cOlGPRwT1vGN93erDpv/Vf9hJvq8Ap5rc8JPPDM0cxfC6D7n88wsMBnfBOFgImbIEJ7fKKKoYf34/7dk20lCrbtnsRD4/pyEmmjivH3WzrYRPuVBd6yLvTz8HymBdvwmUBk2RY/PlHFyNePom7L4LIKhF13r9GO0SsT1p4r0I149SZwUmk1j8lqsbYsP59xYGrDZ4HPpI7Jy6lqhVguWSBCR1ONtTv2K6NNyiULk9tv6TimaamCRQMWBvr7fHfB8nN+aQkok9IPXnBTkvRgpy2TOYL2eNW1QrQLdfv3dKGFcruz/aBpXZpzM/X5z5o2ZzEJEcxCGKEz/LE80Cwl7Ucf7JZibd0Ma/hMIvjVIO3tDc9fnDWunV/xKNr2xJa1vtsPmh734OjGjklNN3a5MnDLmk0qFNJPS0gv8tSqUZL26ioK2KTDZAq/NuIwFAjKAmwy8zRMmWIdBODe9YNtKwcvSlbRSJj6aZ7uF5P+tKalHORYVbHuWRGqXk3sgXyapkwFPkfpMImQZAy0SthbRcLmG5bGIuyBZrtGP8K+SGSsOcfZzMStJSRwufLoE1vWapvCOMfqjIqSk0razVlU43jq0EnX8+82WOAziZB2DPSWX7g6toYZJasIIvda9s7PN3yurOMSQqZJSsNDFSzsN3P5OQW6blJxdvKKGz9N7Ls1SYsFPhMYP7HuI5tWGmmIOoIUArPz7BF9H94wDFgFPHrnakx7rBzs7RWvsAqBzicOIeQnScm0+QvQPkG5jTtJbd/P9UtbQYkLFvhMIGSsu315fP/uSazdsV/58vrREJ2UrCK237qqraqmX8KEG+rYcO1i/P6dN2DHC1Oe2r1dobdH+/ghDiEURVMYHVLAen0vbhOKVEz8rMC6NUmLwzKZQIztmWpLbJJM1+oY+cZRjO2Zwnu19vh4PxqixKsscMkqYkFfwTiKJwrKJQu1egMHXz+Lg6+fTeSYfoWQqZM0bFMYNwpEWDG6F1eWLK0jXRKXCWV8oqost+1EtpeM00GeBThKh/FElfTkBxnBcuD4GaVZxdm0XH7+keHV2rG4xfLHSblk4ccXZ9Fwk14xsLC/CKtY6JhEVXhF3oTBb9SOKX4rWZpOaG4VQHVjyEpEkR848YqJBKeG5FfYA03N6alDJ3Hv+kFlCNzv+xBEbtU25Ut6ZcnC+UuznlqdxCqS0WdlGGXSwh4Azl9qADBLzNI5SR945qj2O6bI67/cpRWlX4KsXkwT1fyYzk5N1zKRBBcnsQt8InoTwI/RfFpnTWYhJj7GJ6oY2zM1bwJZNGBh+63tzUXsGo6uuYZfBIADx8/g0TtXa7WnMJqV0zRRsgrGAr/eEJ4mBwCYEwIXZ+fcP5QQUoBv2z3Zca105hFZ/AwIL7zKJSsSM1oQE4pbKGmY81pWLsW276yQlIZ/sxDiXxM6FqNgfKLa1kVKcm6maXMHoM6ujFCZPTVd02roYTQr1bnZTUQmzAlvTT8rwl7iLO9w+K2zOHD8jOsti0p4jd22qqNAnUTaw70I2pDET71704lJrjK27Z70dcy8wVE6PYAUprrlbb0h5qMkVBpOVLhFa+g0qx0vTLnu0+vc/LCwv68tfDJP1OoNPHXopJFtPQrhNTxUwa6718xHThVbF65SLuHe9YMdkT9OwkTC+AklHbttFaxC+021CoQN1y6eH3ORCHfd2FREur2WfhIavgCwn4gEgP8phHgygWNGTh4qHuowEeJSUIQVBuWShU+tWdqRder1guuOe26mjofGjykduEC0E1SQKCIdQZzbYTE9WnnAwoadL4V+Tt38KeuuWdxmOrSjMiP6waQ3sn2MQHsUkiynIO9PQwg8e6SKddcs9rXvPJKEwN8ghDhFRB8E8C0iOi6EeFn+kYi2AtgKAIODgwkMxz95qXiow0SIE5rnoGuXZ4JVoPlm4+uuWexr8nM77lcPncS6axYrvx/lUvsKq4DFCxdEEoEiyzcnGS5qQoGa7Qrliiiu51TuS9WXIIywt+/b9PlyTkwbdr6ktdNLE1MWFbcoSDQsk4jGALwvhPgD1d+zGpZp0hQ5y42TTQuXSQeaU8OxCoT+vkIrUkSNX61NVcTqq4dOaj+va0ASdVG2+xSRREGQoaW7v/e20s6dNeJ4TrP6TkTV9D1LZKJ4GhEtJKKfkD8DuAXAq3EeMw50AsW+PcuNk0c2rYSJabracqo62+XtunsNygP9yu8UifDElrWY+L1btM5YZ/kFVRGrZ4+4p9YLQJmCr8oUDYOMJDJFVy5CRiXtunsNyqVwZSF0FCN0OMTxnGb1neh2O70bcTttPwTg74noKIDvAdgrhPi7mI8ZOboXy749yw/R8FAF9643M5eNT1TnKyS+sXMzDo5uxPBQRfuSzgnhmgA08o328gsyC1e1pPZCOnHtEwiAjgnqvvWDgcswyEgit+/Lu14pl7Dr19Z47mty+y24z/D6m2IVCY99ek3oGkOSOJ7TrL4TvdzzNlaBL4T4oRBiTevfKiHEl+I8XlzonG/27Vl/iHROTye6qBi/L+/4RBXbdk92hDnWGyKUXfvcTL2jzC6AtgnqkeHVODi6MZDQl+ejWzlsuHYx3ti5GW/aJkPdcezXZu8rp32PxY16Q+D+3ZORRCcR1GWNwxL0nQjSgN4P9lUscLlU9a59JzLfoD0sHJZpgO6Ftm9XmULy1kQB0Gcm+nl5xyeq+J0YukqpcCu8pRPa5ZKFDdcu7jBz2c9H3k+7Bl0uWbh7XaembnJt3ARzmpGg9ph5r8qVfgWx7p0AoN1PUrXzh4cq8/fNmc/QzUKfSysYYBqq5RamlgXCRI3oIiMAdIT4je2Zgt8UJdNkHRU6c9PwUAWH3zqLp7/7NhpCoEiEe266en6146wR5Gz8AQAXbAlc07W6MqLFb9SIk8e3rG3rJStEM0Q0iQnTeQxZp37XvhMdWdBBotCc74SzVIc088nPJpnp2u1ZtSq4eJohWY2x94MsaewWNVIuWZjcbtZHVVeky2+Ey6IBC5tvWNoRHWMVCQv7++YLhrn1h9W11fMqIqbr20pohmnWFBm7RWraz00ygL1KVLhdb100SVLYr1VUETe6Ymayf2+SETTdFK3DxdMC4CbUs6S9PzR+TKu16pDnVp8T2qQgGUdvik5D8osQTR+DV+y+SjhbRcLZ8xfni3nZw0NNNDhd4pYAlMIeaPputu2exP27J5W1YJz1iuZ3qOBTa5Z2bJP3Km1VzH6tooq40Zm25HZdPkZcTmXVsQpE88ELQHcoexIW+C2ynDhl56HxY23x6g0h8NVDJ/HVQye1haic59YQAiWriLturODA8TOBH+SowuukYPSaVJ2mk/KAhfdm6m2OYXttIN34pINyZNPKwOfgtHvLOjbV6Zov85Sz7aJuxZEW1ekaVozuRUGjJEQtiJPMdNXV97cXmQOQC7lgCpt0WmQ1ScTJigf3ulavVNU993NuOm1GtV1XF72/SBCAcbVKORY/E874RBUPPHNUG0ElHepuSVlRNk8J6oNwmg+iTiSLkyA19n/2v/1v5crJbtry8wxGUd9f9xy5PUNZkwts0vFJVpNE7IxPVD1LFaucTqbnplvlHH7rbJt9XW6/68YKnv7e2x314S81BK774EK8dub8/HgHWvZw3fDdNCdVVq69Fopuf17U6g1cYRWM6+G7EfTbUkOOq7FI1BSJMCdEIIE7PlHFrMJ/VEDTtGW375dLVltmdVwr8OGhSqAKmVHKhSRNRizwWyRpOwyKad9P58Noem46m7f0Fzi3P3fkHcxpHMA/ePd82+8m5YpVk5UqqsOtBINfoohjD0N1uoahh/fj/QuzgUswhIlw8ktDCDyhKHFhwq59J5QTa6m/iN3/+Hbb36ZrdYx8PZ7oHaeALQ9YyudgmYuGH5VcSNqUzHH4LdJInPIb12yqVTgfRtNzc2ucoWLGRWMPinMMO16YCq19Z51zM/VQ9XYEgpdZuG/9oO/vBo1V1z1f5y81lPe4Pne5bHdUK3BVnP/7F2Y7SmTI9yNuueA2kcUBC/wWSSdOBUkwMdEqdPkBJuem23+UNVu8cI4hbQ1cxaIBK3DphrgIUop5YX8RjwyvxmOfXuMr+SuoQCoHKAMhBXpUZRpUArY+J7Cwv0/5fsQtF5I2JbNJx0aSoZdufUdVbesAfVSBbM/n5vg0OTddhMRdN1YiM6OUSxYWLuhTRrNkqRSFs7G6nXMzdQz05//VsYpNfU8mqDl7GLgRRLN+/8KszxG2l7qIInpHN+73anVtPkScciFpUzJr+CnhZj7RafxS23BWX5wTlx/+sA4slTbzyPBq105Q9u5BbljFZpz/wdGNeHPnZjy+Za2r5pRWiruJE9fEuZr15lnTtfq8KfGR4dVt98OLIJp1ELOVs9RFWE07awXdkjYl519NySkmjUZUTinpwHKGEkaVEq7TZu69aVCr5f/Tyffw2KebVSN1MeSqevlemlNcdkwvFvb3hQ7VlCGLQHveQBjnbBzYnYSmRKlZe+HneTEhax2twpbl8AsL/JTQmWecqF6UpOx+9miGPpe1YFTdgpzRE2mEKFZCHFdmMDtNa85Vi71+T9JtEFXIstMX6nNGCV+y/6sfgtzPOPwkQQVsnKGTSZqSWeCnhPPBM81kHJ+oJpL16AwX84qqlJNN0IdXFZ6WNAR4dt5y++7rj34SwGXhsG33ZFsxNCko7Ak74xNV3K+JA08SP87xva+cNi63LTFVcCQExKZ1+31G85KFbwIL/BSxP3i6Ql/yoVfWZ7HhtSz1q6H4bQ4eJAIjzPH8QAAG+ouuLRqBpgP5xaPB6tbbE6js99EuSFWCYnio4npfs0iQyCmngnOlR+VWgewI026qqskCPyO4LTW96qsUiVwdWF4aimoy8Gseev/CbFvBKR26icfP8UpWET8/eCUOvn7W6PMC8BT2QPDy0fbJ1mviUgmKsdtWZap+TlSo7rV9dbN2x37X6z308H5Mz9Rjt2t7kYcsfFNY4GcI3VLTS4i4tRnUfd8eS62aDLw0MCcyScZtHKqJ5/7dkxjbMwWrSLjkEhlTKZfa6sWbCntTrCK5RiK5jcsujEyEgP0zUijW6o15m/6iAQvTM+p6+PJ4cZRhiCpjd3yiit997pW2sFbV6mbstlWu5iy5kqhO17Bt9yQOv3XWtynJZKxeK988ZOGbErvAJ6JPAPgfAIoA/pcQYmfcx8wr9ofvypIFImhffDteD56bhqKbDK6wCr5r23sJIN3E5TWxFIlw8/VLsPeV07ElYlkFwrTPfasKaJk4J3XmH1nFdPutzRLVKhPf8p8sYVtM3cRkxm5DCFfhX3CZGN16LjhXN8NDFWP/hQDw1KGTWHfN4sg0fVPbfNYie8IQaxw+ERUB/AmAXwXwEQD3ENFH4jxmXnFm3k7X6jhnIOzdHjxZukG3j2UtrVnFuZk67rrR34tFcI+dD7oEliWg48y6nanP+dLYrAIpr7uuraJEOoYBb9uwM+78rhsr+IfXz8ZaN8dL2APNvA8dXvH2zmfAmVPihkC0obqmZQ2SzsKPk7g1/I8CeE0I8UMAIKKvAbgdwD/HfNzcEcRpqYptl3jZ/eVE4WYa2P29t33ZteULqXsR0gq1NMUt09i+uiiXLIzdpr7uTl+M01ksADx7pIp11yz2tA07TXxuk3eUeB3DLVzSa1K3T6rjE1X86IK/STxKu7lJv4QsNkAKQ9wCvwLgbdvv7wC4KeZj5hLTB5mA0FE2TruzzkRQnxMg6mxb6GbqkS+LyiYaNOQxSkjTarBkFVwd535sx3bhsGHnSzh/qf3eSi3Sr204C05CL1OG26TuDLXcte+E62pBt/+ocBtrnkMv3Yhb4KusfW23mIi2AtgKAIODgzEPJ1vYbfa62Ho7fpou6IQDAW378LKjTs/Uce/6wbaWirJTluplIVy25TtfGmd3pzTQZQzPzon5KKOgeQTOmv26awQ078/jW9b6sg1nYYXklXQ1smml1oZ/7/rBtjBkv+cStd3cKzcgr6GXbsQt8N8BcLXt9w8DOGX/gBDiSQBPAs2OVzGPJzOoHHZu+H3YTbXH8Ymqa8bnQH+xrdlIQwg8e6SKu26sdDQdV9l+TfqiJkVfgbQrjHpDuBauc0Pl/PNaySwrl7QrCgDzq6QrSxYuzTaM+gkkgdekLc/pwedeme9sJRvCP3XoJA4cPzPfwMYPXqHHQbBff7eJuZuIu3jaPwK4johWEFE/gM8A2BPzMXOBzuRSJAKhaSdeNGAFdhKZFGWSgsptspm51FA6tg4cP9PhyNLtxavEbVKoui3ZcStc54Zf/4v9PgwPVXBwdCPe2Lm5pXG+gvt3T7Y577Mi7AE/AvDy4l42hJfX9qlDJ337q7xCj4Mir/8iTeJg2s9s1MSq4QshZono8wD2oRmW+WUhxFScx8wLuhdnToi2HqdBMakZYiKo3IS4yqnotqoY2bQytpBCN4LEl9tXJl6x2n61wAV9BRx+62zbPgf6Cx1dwvwSpg6QKSYC0Ou5CnL/r/QRzQP4yyzXlW62iupIrDwTexy+EOJvAfxt3MfJG0kkc3jZo8MsV1Xj9IpXDlJ3PSxhhOCp6ZpRrLZf2/p0rd5m8olKSHvtp1Iu4dz5i6FWDCYCMA4ziJ+kOL+1b3a8MKX0OSzs7+sq+z3A9fBTQ8Zim26Pg6CTi86f4BWvPD5RxYHjZ0K15DPFKhLKJQunWlUpg1AgwtieKc9Yba/Y+yywaMDCyKaVSsFmuWVSOfZhIgDjMIP4SYrTxdfveKHTuDA+UdXmd7yXo/pGpnBphZTQOb+SjGTxU8FQljbwWh7rVhW6rNI46scUWjYcmT8QtARxQwhtDoJdi3U6/5JsKm7K+xdm8YXnjykbu3zgij4M9Pe5rhDsGcBujE9Ucf6i/85WXviZRNySCZ31ntwSubrNfg+wwI8dv8XC4owKUI3lrhsr8yGXOopExuGgOnRaV9Q14QlNG3lNYbYoEmFOCKPKmV44hYGz8qm96UkW+vLW5wTqmnOenqlj+62dBdzkxOXWOtOOV7KfKWFbX7qZ2Jxhlm7vW7fZ7wEW+LHiZktMuiCTaiwj3ziKRkPAy6J7z01Xe3zCG90LGHUDEBkRokI6xK99MJxLyUsAOROvsiDw3XALEfVjw46ixLXMbD5w/EzgcYxsWqnNLXEKeN17WC6Zma/yBgv8GHGr1ZF0QSbVWLz6tgLAwv4i1l2zOPTxs9DdSU6mYcehC5FVraDSTpQywe5UDyPkwq5OTVcSXgwPVdryAOw4o3107+HYbd7mqzzCAj9G3Mw2SfeyDPoynr/UiCTFPG1hbxUIM5dmsWJ0b6j9ODU/e8t84H/1AAAZxElEQVRCuylCrqCiIEid/gLguXKT+w5zX70a8wBNZ6+ufeKiAQubb1g6r9FLm7rfMTkn29mG+uwvzbaPIen3MG1Y4MeIl9kmyYJMYdLyo0gxTyJG3M6iAQsD/X3z2arnL81GYlr50YXLjr/xiSpGvnF0fqXknNJMVlAmLFxg3lSd0CxhYFJK2kSTdYtndyuFbOf9C3Vli0yp0YdtH+inPaYqJLVbCqOZwGGZMWKS7ZoUYcM9wy7XRzatNA7/C4uMKBnZtBLLyiVM1+qRCd85AYztaYb37XhhKrL9uuFnovzFaxfjwPEzrsLeNHtbCnSZ9VudrmHk60fnM5C9SiFLdGH/bv0Y/JRBjrM9ZrfBAj9GnHHp5ZKFK6wCtu2exIadLxmn7kdB2HDPsM7k4aEK+vuSedwevbNZ2VL2F/DC7zw0Xatj7Y79mXTGHnz9rOs5V8olvLFzMw6ObvTUasf2dCYk1efE/IQXdsXm1o/Bj4Lh57MJ6RyZhQV+zMhaHY9vWYuLs3PzTU2kjXftjv1YMbo39gnA9KW4b/2gMono/MXZUOMbn6iGDoX0gx+tb040VwV+mnHkqem4xO/qUneO07WmWctUdhZaJbZVY9EpEn4UDD+f9VuOudvoKoEvOzwlIUD9oouSma7VAxXs8ovJS1Epl/DI8Go8eufqjmJS07V6qPFF2anIC1l8zA+1eiOXQtwUnQlH98543ecvPH/MOLnsP9w0qM3AjsLsqdqH22SUNdmQJF3jtPVbPyNpTDTsWr2BsT1TsYzXK6vWWfNm174THSaLMM7bbiszmyd0fRTc3hmvCdrPak02jzHpEBYkSkZ+dscLU/PP7BVWAbNzQuljyZpsSJKuEfhe/UHTxjRKRi6Xox6z88WyN0n3UwEyqODOQvOOoBSoaQqolEs4f3E2sZVAFCUa3LRl3Ttj2ljcBLd2iJKoomQu2LzDtfocrAJhkSbTOUuyIUm6RuCnUarAD37q1sT1IPp5saLOBPZz/nEghXYQpI1fCk6TUMQoEACe2LK2TXOVyIzUF4+enp+AnDHtXtpy3O9GXBFp4xPVtmtSbikvHSbTOYGB/j5Mt/xmTrIiG5KkawR+0qUK/OLUsN3ERRYexKgzgVVLd9kGUP6+/CdL+IfXz0ZaeGzAKkCAOs7j5wevxMHXzxrvR2qE0jTilWwUFXKS1sXDe/XalTZ61aruygAJXaZElTXrxJn7ALg70OX1yqJs8FOzPypIpJwBaWfdunXi8OHDgb6rKtxUsoqRt0WLirU79msfVD+9a+PEnkUqSyPE8SLbj5MUQUs92JvIA+i4PlFSLlmY3H5LoO+aZMDGRZzPr67JjttYdMpLmrIhanlFREeEEOu8Ptc1Gn7eUqTdSrRnpUrf5f6k8TnDo6qw6JegwtmegAS6nE3bEAJWgTAHoOEw9xQIEMKfPd4qUOB6LqYZsHEQd2Khn9WvHEsWZUNaPseuEfhAvlKk3Ro6ZOkc4n4ws5ol6eUwVQnT+pyYtyfb7ctjt63y7PQ1YBXQ31fEezW1E92O3RRQsgqozc5BiOaq5Z6brsaLR0+nIuzjMuPYcXP+l0sWFi7oUwr1rMmGtHyOsQl8IhoD8J8AyBTP3221O2Sgf3BNohqSwMvMEtWDmQV/hQqBYFEy79XqHT2JxyeqePZI1XVfM/U51BsCj29Z6yqYHho/1jZx2GvDNIRoa52YJEmZIUc2reyw4QOXV0RZEupupOVXiFvDf1wI8QcxHyOXJF0e2Q8mZhZZZtaruJbXMtpNY6s4bOVRNvw2IYiOXB6w5p2k8pxNVzH1OYHffe4V1+uZZD9gU7yeW+dzIJ31QXxDqph7uYrKi7AH0nv/Y3PatjT89/0I/DBO2zyShpfeBBPH2KIBS9klSTqeABg5pUydV7oyxFnCKhIg2s09Yds42q/F0MPJ1e9xXmNCszDbm/9WUwpqQG0j9+OjcTY/KQ9YEAJGZq48EuX7b+q0jVvg/0cAPwJwGMADQohzbt/pFYGfVUEvWTG611OgymgV3cSgi1pRLf29rkdajl0/FInwE1eoyxiHjeCRQjXKhCgdVpGw69fWADB3cjrNTEDTxPKBK/p8T1Buk3nakTVZJpEoHSL6NoCfUvzpCwD+DMAX0bx/XwTwGIDfUOxjK4CtADA4OBhmOLkg6yUgALOsWLdKh4A+Ckb1HS+HWlYdu3bmhMB7mhDIsOGa9sYgQTBdEanMKlLo6xqT6MxM9TkRaDXiNs5ezY6NklACXwjxcZPPEdFfAHhRs48nATwJNDX8MOPJA1kvAQG49wQFmgJkZNPKQHHeAk2TkZ9VTVYcu26Cs+CixYc1QXlNrl779zr2mw4nM2CumOzadyJR81pWnoW8Elu1TCJaavv1DgCvxnWsPJH1EhBA84V2Vsu0c+/6QQwPVVxzCdzwWxk0jYxI2b9g0YA1X+HxF6/V9/Z10+LDCsSbr1+ivQblkoU3dm5WCu0wmDYmSfq5TTs71i9Zq+AbZ3nk/05Ex4joFQA3A9gW47EAZO/iqoii/ncSbL91lbLk7H3rB+fT+d1yCSRFzazgp6vRyKaVxrXXo2Thgj5sv3UV3ti5GSObVuKfTr7X8ZkkxnXg+BltGeFPrVk6/8zrJuBFA5Z2AtdtN1VMon5u3a5nVqLYTJGrJHvHsDhLoJsQm8AXQvy6EGK1EOIGIcRtQojTcR0LyObFVZGltoduqLp1lQcsPHXo5Pxk6vWyl6wiHvv0Gu1LbKodDg9VcO/6wY79xCls5TN0/+5JLB/di/t3Tyr9CFGYMwju5yKb3jtryt91YwXPHqnOP/O6RYYQwOYbljajiGxYRcL2W9XZvF6KiVSutE77AK2lZCSQPMdFA1YzkQ1mLRmzRhTtG6Oma2rp6B6+rNSlsaOLS85q1I6qYJVVJGz5havx7JGqUhDaw/Z0JYX93hvVddMd34SF/cVEu3DpGLAKuPPGD2uTpnQ1dXTPPFGn8HeGPJpk8/oJue0Y27WLMXXqx759PN0UiaOLdiOgIzkvLD1XSycPtnGJPSolD1E7qmbd9YbA3ldO464bK8ooDWnTrk7XUCwQrAJ1xKf7XdWoonnWXbM4cHz+zKUG7ls/iKe/+3bgSBq3ssum45mpz+EplwzZ85dmlT0SdM+26lRq9Qae/u7bmBOirfibM0nMXooA6Kxualrk7tAPz+H1Rz+Jh8aP+cr+zVoAQxiyWKWza1oc5sU27iSLyz4nuvC6czN1HDh+xlOoNeYE6nNi3p6/aMDCgr5omrnLnsFv7tyMx7es9fXdAhGeOnQSV1jBXwOnsJd29Eq5hMe3rDUuleFat6chlM+D33E3hJg3Vf3OM5MY+frRNhPo/bsnsXbH/vn7Ia+t9GFI85HpsYCm/0GFlwkrbaLwB2bRfNs1Gn6WSxW4kaeViQo/42wIAatIeP/C7Ly2L1c0h986G9qsJVszKs0c6BSqUihFadZZdmWnmUpV+8Uv9ussSx/XbHV0/DInmrkDTmTvYqAz/NKP6UxO7tpVCJqTYtY0YCC6VXcWq3R2jYavcmrlwRaYh5VJuaSO5CiXLN/jrDdERyXHWr2Bpw6djMThrtOq7l0/OP9s6CKHosAp4IaHKrACODCdyBo90oEcZ537KMIv77npagD651gmeWVNAwaiXXXbV0kHRzemLo+6RsMHslcC1YQ8rEzGblvVUV/dXq89irIHTl0zqC3XRKtaMbo31FjdWFYudTiXZ0Jo4pL3L8z6ylwNm+ylCr/UmXMGFCWaZeiu2/Odlgbs1dgn76tuN7pC4Ge9No0bXg+9SZ2ZuM/d5MWUf6MQvWOdBH3BvCb+uBqql6wibr5+SYc5IAr81rcPewucmnnQrlFez07SSprTXGMPLhj5+lEA2XS2RkXuwzLz1trQD17nlsVzV4VwBiWukFrVdSsWCHNzIpCgtLc9TLpVYxzonqE8K1YSr0qwstRy1t4rL3omLDMPtWmC4mVLfOCZox3hhGmfu12j071YbqGMkjjNWl5ap0m1UIlzUkqioqWKRa1SwkFs+wv7iygP9HsK8jyaTJ14rRqna/VMOlujIvcCv5vtbbpzkE5NPxUpk0QKBp3gNLFOLOiLN57ATXiZmnxUk1IUzcz9JIQRmrWNHhle7eqb0Nn0rSLhS3dkV3ONGtN72w2Tm4rcR+nkIcolKLpzKBK5OkkLRJmoJ+Q2fhX2rTI8MI3xu9XuKRK5RoGZCnu5H2eBtie2rMXUw59wjd+X10/G+ksHqVtEjD0nwP79Xb+2xlWw5aE+lR9UkUF23IoGdgO51/DzEOUSFN25eUXE2B1RaWbt6sYva8DYt6s0ULv5yi2qImqGhyrKpuMmdlxdbLmTOSFc0+uDOEm9ImL8XqssZYFH5T+Q31GV9narLdQt5F7Dz2v8vQm6c/PT6DzNrF3d+B8ZXt2xXacXSyEjhahzMrNrnOMTVazdsR/LR/di+eheDD28P7BG+sjw6nmt2M9z5aVBSrxWoEGe66jfhaxkgUdVGFGuVrbtnsTCBX24z5abYbLa6QZyH6XTi+iic3SafxzFmqJGFz3hZROXTtPxiWpHrgBwuWVf0qF/UhsdUNjjsx7xIQlb/CsqrTyKwohZjGiLEtMondxr+L2IX80/D/4MXdall01cOqh37TuhjFXX1aGJE3t25dTDn8ATAVYKWSCMfyzKcuVRBGZkZbWSNrm34fcqKpuszu6cB3+GLhTOK65dCh+3l9/tb+MTVex4YWo+i1XGYUcpkPMa8RHGPxZluHQUiVDdHM3nBxb4XcL4RBXPHql2LMEJAtt2T2LXvhPaJXVWEmp0k5hbed2br18CALiyZGlj0MuayAtV6d7pWn0+4zKKa5CVaxuEMPHopgLW5PpEEZjRzdmzfmCB3yXoqhnKOi66CIssRWKo0JXXtf99fKKK85dmtZ85N1PH0MP7sfmGpfMVOd0miPqciCR5LevX1gTT1YlTcJcHLGXtH7uANb0+USRCdXM0nx/YadslmGaHOh1dWe8U5nVesqxB1OUMgji6O4qmXVIXPMvKtfWLThtXOUStAgGEthIbMvS24mGui7OkRl5XW14kUlqBiO4GMAbgZwF8VAhx2Pa3BwH8JoAGgP8shNgX5liMO6ZCz7mkzrpt0+u8lpVLsYzV71Jfpa3q0I03ywLJTRtXrS7rcwLlkoWFC/o6upHJ7+qiyuJ69vLqS4mSsFE6rwK4E8DL9o1E9BEAnwGwCsAnAPwpEXkHJzOBCRr/nfVMZbfzkkvyqMdqFcj3Ut9PgxDVeFVRLdtaDdTtGa4mma9xZMe6OWF1Avq9Wh0HRzcq8yxq9YY24zorz143EkrgCyG+L4RQxTXdDuBrQoiLQog3ALwG4KNhjsW44wzVXDRgdTTeUNkss9qEQmI/L6C9LIAMb3QrheCXAauAXXf7j9s31Up111YlUJ0a8UPjxzxDHaMMh7TjthL0Uhp0320IkelnrxuJy2lbAXDI9vs7rW0dENFWAFsBYHBwMKbhhCPLS207ziWrybjzUBnQayk+PFQJVKXSKhA+cEUfpmfqoc9bZ3qSZg2va+s1Ycgm5F7VUaMKh/TjhPVyiOqujd2Wn9Vnr9vwFPhE9G0AP6X40xeEEN/UfU2xTel7E0I8CeBJoOm09RpP0uQ50sLUZtkNtk1dDZsiEeaEwLJyCTdfvyRQ39wwoYO6mH5TgWrHpDpqFD4Z1TNvFQhWkdqcsKadq6Ku8cMEx1PgCyE+HmC/7wC42vb7hwGcCrCf1Onmevs68rKisRO0I5MXD40fa0tmiyJ0ULVPE3RlJuwmlSjizb2csH47V+VhFdkrxGXS2QPgb4joDwEsA3AdgO/FdKxYyXoUiymmQjyvK5o4hMr4RLUjcxnQT/huQs/eR9UNGc3irB6qqzLqtHmrJj5C8z5u2PmS0TVxc8JObr/F9bs6WJPPBmHDMu8A8EcAlgDYS0STQohNQogpInoGwD8DmAXw20KIcF2uU6IbMvT8CHHdimbHC1OZ19CiFiq79p3Q5gCEMZG4YY9Tt5d8WNBXwLprFmPdNYtd74N94tOFQ9o/Zx+j3G/BYCXB5JNQAl8I8TyA5zV/+xKAL4XZfxbohgw9P2YpnSA7N1OfFz52wSH3n+WJwA23lY+bUA9rInFDHvdCK0sauNwQ5tE7V3smJcmJT5VUp7rvusbedvL2zDNquLSCB91gf/RjljJN4KrVGxjbM4WLs3O5M/9IvFY+umtBgC/h59f8t6xcisR3ZHrfdROS3eGdt2eeUcPlkQ2wl7s9OLoxdw++n+Qq0wQuoKl15rnkrFfJXNW1kD1k/TwDflYDUpOOwndket91+5RdufL4zDNqWOD3AH6Sq1S19sslf30+8+LQ9hKqqmth7yFrim7iAJoJcuWS1VErP4oMaNP7nvVsayY62KTTA/g1S6kSuFR+jCusgmdFxCxj4pCPwhEcxCwYhe/I9Ljd4KdizOBqmYwRKucmgFy3jct627sk8yHymHvBXMa0WiYLfCYUeRcUeR8/wwAs8BmGYXoGbmLOMAzDtMECn2EYpkdggc8wDNMjcFgmkzrsOGWYZGCBz6SKW3kDIN8lLRgma7DAZ1LFrTrnhXp+6/QwTBZhGz6TKm7VOfNcp4dhsggLfCZV/JZhyEudHobJIizwmVTxU50TyE+dHobJImzDZ1LFWeBL120JSK+gF0cRMd0CC3wmdewVKVeM7tV+LkxRs6BCO689fhlGRSiTDhHdTURTRDRHROts25cTUY2IJlv//jz8UJleQGeyqZRLoYT9g88dQ3W6BoHLQnt8our5XV0U0QPPHMWK0b3YsPMlo/0wTBYIa8N/FcCdAF5W/O11IcTa1r/PhTwO0yP4adZiildnKzd0TuKGEL4nj25mfKKKDTtf4kkw44QS+EKI7wshOE6OiQxVl6mw9enDtAs0cRL3erhomBUUkyxx2vBXENEEgB8BeEgI8f9UHyKirQC2AsDg4GCMw2HyQhRdpuyYdLbSoeoGpaKXw0WjaLjOJIOnhk9E3yaiVxX/bnf52mkAg0KIIQC/A+BviOjfqT4ohHhSCLFOCLFuyZIlwc6CYVwIYyZyrjiKRMrP9XK4aBQN15lk8NTwhRAf97tTIcRFABdbPx8hotcB/AwA7m7CJE6QnrLO78vP6toimvoYujHEM8wKikmWWEw6RLQEwFkhRIOIfhrAdQB+GMexGMaEqMxEYSaPbg3x5Cbo+SGUwCeiOwD8EYAlAPYS0aQQYhOAXwLwMBHNAmgA+JwQ4mzo0TJMBgg6eXSrrTvsCopJjlACXwjxPIDnFdufBfBsmH0zTLfRzbbuqB3tTDxwLR2GSQidTZtt3UxSsMBnmISII6mMYfzAtXQYJiHY1s2kDQt8hkkQtnUzacICn+l5ujE2nmFUsMBneppujY1nGBXstGV6mjCVNBkmb7DAZ3qabo6NZxgnLPCZnoZj45leggU+09NwbDzTS7DTlulpODae6SVY4DM9D8fGM70Cm3QYhmF6BBb4DMMwPQILfIZhmB6BBT7DMEyPwAKfYRimRyAhRNpjmIeIzgB4K4JdXQXgXyPYT5bgc8oHfE75oNvO6RohxBKvD2VK4EcFER0WQqxLexxRwueUD/ic8kE3npMJbNJhGIbpEVjgMwzD9AjdKvCfTHsAMcDnlA/4nPJBN56TJ11pw2cYhmE66VYNn2EYhnHQNQKfiO4moikimiOidbbty4moRkSTrX9/nuY4/aA7p9bfHiSi14joBBFtSmuMYSGiMSKq2u7PJ9MeU1CI6BOt+/EaEY2mPZ4oIKI3iehY694cTns8QSCiLxPRu0T0qm3bYiL6FhH9oPX/ojTHmBRdI/ABvArgTgAvK/72uhBibevf5xIeVxiU50REHwHwGQCrAHwCwJ8SUbHz67nhcdv9+du0BxOE1vX/EwC/CuAjAO5p3adu4ObWvclrGONfofme2BkF8B0hxHUAvtP6vevpGoEvhPi+EKKrGpG6nNPtAL4mhLgohHgDwGsAPprs6BgHHwXwmhDih0KISwC+huZ9YlJGCPEygLOOzbcD+Err568AGE50UCnRNQLfgxVENEFE/5eI/n3ag4mACoC3bb+/09qWVz5PRK+0lt55XVp32z2RCAD7iegIEW1NezAR8iEhxGkAaP3/wZTHkwi5aoBCRN8G8FOKP31BCPFNzddOAxgUQvwbEd0IYJyIVgkhfhTbQH0Q8JxIsS2z4VZu5wjgzwB8Ec3xfxHAYwB+I7nRRUau7okPNgghThHRBwF8i4iOtzRmJofkSuALIT4e4DsXAVxs/XyEiF4H8DMAMuGACnJOaGqPV9t+/zCAU9GMKHpMz5GI/gLAizEPJy5ydU9MEUKcav3/LhE9j6bpqhsE/r8Q0VIhxGkiWgrg3bQHlARdb9IhoiXSoUlEPw3gOgA/THdUodkD4DNEtICIVqB5Tt9LeUyBaL1skjvQdFTnkX8EcB0RrSCifjSd6ntSHlMoiGghEf2E/BnALcjv/XGyB8BnWz9/FoBuNd1V5ErDd4OI7gDwRwCWANhLRJNCiE0AfgnAw0Q0C6AB4HNCCKcDJ5PozkkIMUVEzwD4ZwCzAH5bCNFIc6wh+O9EtBZN88ebAH4r3eEEQwgxS0SfB7APQBHAl4UQUykPKywfAvA8EQFNWfE3Qoi/S3dI/iGipwF8DMBVRPQOgO0AdgJ4hoh+E8BJAHenN8Lk4ExbhmGYHqHrTToMwzBMExb4DMMwPQILfIZhmB6BBT7DMEyPwAKfYRimR2CBzzAM0yOwwGcYhukRWOAzDMP0CP8fm/DGKxXor38AAAAASUVORK5CYII=\n",
"text/plain": [
"