{ "cells": [ { "cell_type": "markdown", "id": "cbf4ce41", "metadata": {}, "source": [ "# Source model interface\n", "\n", "`icecube_tools` has a simple source modelling interface built in that we demonstrate here." ] }, { "cell_type": "code", "execution_count": 1, "id": "17d9c882", "metadata": { "execution": { "iopub.execute_input": "2024-03-22T09:57:01.076068Z", "iopub.status.busy": "2024-03-22T09:57:01.075685Z", "iopub.status.idle": "2024-03-22T09:57:01.755046Z", "shell.execute_reply": "2024-03-22T09:57:01.754346Z" } }, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "from icecube_tools.source.flux_model import PowerLawFlux, BrokenPowerLawFlux\n", "from icecube_tools.source.source_model import PointSource, DiffuseSource" ] }, { "cell_type": "markdown", "id": "97e7a920", "metadata": {}, "source": [ "## Spectral shape" ] }, { "cell_type": "markdown", "id": "ed90839a", "metadata": {}, "source": [ "We start by defining a spectral shape, such as a power law or broken power law. Let's start with the definition of a simple power law flux." ] }, { "cell_type": "code", "execution_count": 2, "id": "20db9489", "metadata": { "execution": { "iopub.execute_input": "2024-03-22T09:57:01.758170Z", "iopub.status.busy": "2024-03-22T09:57:01.757906Z", "iopub.status.idle": "2024-03-22T09:57:01.761478Z", "shell.execute_reply": "2024-03-22T09:57:01.760901Z" } }, "outputs": [], "source": [ "# Parameters of power law flux\n", "flux_norm = 1e-18 # Flux normalisation in units of GeV^-1 cm^-2 s^-1 (sr^-1)\n", "norm_energy = 1e5 # Energy of normalisation in units of GeV\n", "spectral_index = 2.0 # Assumed negative slope\n", "min_energy = 1e4 # GeV\n", "max_energy = 1e8 # GeV\n", "\n", "# Instantiate\n", "power_law = PowerLawFlux(flux_norm, norm_energy, spectral_index, min_energy, max_energy)" ] }, { "cell_type": "code", "execution_count": 3, "id": "03647354", "metadata": { "execution": { "iopub.execute_input": "2024-03-22T09:57:01.763772Z", "iopub.status.busy": "2024-03-22T09:57:01.763578Z", "iopub.status.idle": "2024-03-22T09:57:02.472037Z", "shell.execute_reply": "2024-03-22T09:57:02.471348Z" } }, "outputs": [ { "data": { "text/plain": [ "Text(0, 0.5, 'F $[GeV^-1 cm^-2 s^-1 (sr^-1)]$')" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAG1CAYAAADpzbD2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABwGUlEQVR4nO3deVxU9f4/8NcsDPuioiAK4oImsgwioia54b6kZZlaoZW3BZci62r3lnpbLO2WC1OW5ZJb7mZmXa9cSXAJZXNBxQUFF0BCtgEGmDm/P/w13whBZmA4M/B6Ph7zSM45c+ZFB4b3fD6f8/lIBEEQQEREREQPJBU7ABEREZE5Y7FEREREVAcWS0RERER1YLFEREREVAcWS0RERER1YLFEREREVAcWS0RERER1YLFEREREVAe52AEsnU6nw+3bt+Ho6AiJRCJ2HCIiIqoHQRBQXFwMDw8PSKV1tx2xWGqg27dvw9PTU+wYREREZISsrCx07NixzmNYLDWQo6MjgPv/s52cnEROQ0RERPVRVFQET09P/d/xurBYaqA/ut6cnJxYLFE1FRUVWLlyJQBg3rx5UCgUIiciIqK/qs8QGgkX0m2YoqIiODs7o7CwkMUSVaNWq+Hg4AAAKCkpgb29vciJiIjoD4b8/WbLkpFUKhVUKhW0Wq3YUchMyeVyRERE6P9NRESWiS1LDcSWJSIiIstjyN9vzrNEREREVAcWS0RERER1YLFEZCJqtRouLi5wcXGBWq0WOw4RERmJo06JTKiwsFDsCERE1EAslohMxNbWFunp6fp/ExGRZWpR3XCTJk1Cq1atMHny5Br7MjIyMGTIEPj6+sLf35/dJtRgUqkUPj4+8PHxeei6Q0REZL5a1Dv4vHnz8N133z1w34wZM/Cvf/0LaWlp+PXXX2Ftbd3E6YiIiMgctahiafDgwQ9cA+b8+fOwsrJCWFgYAKB169acRJAarLKyUj95aWVlpdhxiIjISGZTLB09ehTjx4+Hh4cHJBIJ9u3bV+MYlUoFb29v2NjYIDQ0FAkJCY3y2pcvX4aDgwPGjx+P3r1746OPPmqU8zZUpVYndgRqgIqKCsyePRuzZ89GRUWF2HGIiMhIZtN8olarERgYiBdeeAFPPPFEjf3bt29HVFQU1qxZg9DQUKxYsQIjR47EpUuX0K5dOwCAUqlEVVVVjeceOnQIHh4etb52VVUV4uLikJKSgnbt2mHUqFEICQnB8OHDaxyr0Wig0Wj0XxcVFRnz7T5UcXklJn1xHM+EeOLFgZ3rtdAfmReZTKYfHyeTyUROQ0RExjKbYmn06NEYPXp0rfs/++wzzJo1CzNnzgQArFmzBj/99BPWrVuHBQsWAABSUlKMeu0OHTqgT58+8PT0BACMGTMGKSkpDyyWli5diiVLlhj1OobYlXgTV3JL8MFPF3Dy2u/49KlAuNhx1XpLYmNjg507d4odg4iIGshsuuHqUlFRgcTERISHh+u3SaVShIeH48SJEw0+f0hICHJzc3Hv3j3odDocPXoUPXv2fOCxCxcuRGFhof6RlZXV4Nd/kBkDvPH+472gkElx+EIuxqyMQ+KNfJO8FhEREdXOIoqlvLw8aLVauLm5Vdvu5uaG7Ozsep8nPDwcTz31FA4ePIiOHTvqCy25XI6PPvoIjz32GAICAuDj44Nx48Y98BzW1tZwcnLCpk2b0K9fPwwbNsz4b6wOEokEz/X3xp7XBsC7jR1uF5bj6a9OYs2vV6HTce1jIiKipmI23XBN4fDhw7Xue1g34F9FRkYiMjJSv2qxqfh1cMaPcwbinb3n8GPqbXz880WcvPY7/v1UINo4cHoDc1ZaWgofHx8A928isLOzEzkREREZwyJallxdXSGTyZCTk1Nte05ODtzd3UXJpFKp4Ovri5CQEJO/lqONFVY9o8TSJ/xhLZci9tJdjFkVh9+u/W7y1ybjCYKA27dv4/bt2xAEtgYSEVkqiyiWFAoFgoODERMTo9+m0+kQExOD/v37i5IpMjISaWlpOHXqVJO8nkQiwdS+XtgX+Si6tLVHTpEGU9eeRPT/LrNbzkzZ2NggOTkZycnJsLGxETsOEREZyWy64UpKSnDlyhX91xkZGUhJSUHr1q3h5eWFqKgoREREoE+fPujbty9WrFgBtVqtvzuupejZ3gk/zh6Id/edw57kW/j0UDpOXsvH51OUaOvIbjlzIpPJoFQqxY5BREQNJBHMpH8gNjYWQ4YMqbE9IiICGzZsAABER0dj+fLlyM7OhlKpxKpVqxAaGtrESe/7Y2ZmrVaL9PR0FBYWwsnJqUkz7DydhXd/OIfySh3aOlpj5RQlBnRzbdIMREREluiPMcf1+fttNsWSpTLkf7YpXM4pRuTWJKTnlEAiAeYO9cHcYT6QSTmJpdgqKyuxZcsWAMD06dNhZWUlciIiIvoDi6UmJHaxBABlFVos2n8OO07fBAD069Iaq54JQjsnjpMRk1qthoODA4D73cz29vYiJyIioj8Y8vfbIgZ4m6OmvBvuYWwVMiybHIjPpwTCTiHDyWv5GL0yDkfT74odrUWTyWQYM2YMxowZw+VOiIgsGFuWGsgcWpb+7OrdEkRuScLF7GJIJMBrg7vijfDukMtYFxMREf2BLUstWNe2DtgX+Simh3pBEADVkauYtvY33CksEzsaERGRRWKxZCRz6ob7KxsrGT6c5I/VU4PgYC1HwvV8jFkZhyMXc8WORkREZHHYDddA5tYN91fX89SYvS0J524VAQBefqwL5o/sASt2y5lcaWkpAgMDAQCpqalc7oSIyIywG470vF3tsfvVAZgxwBsA8NXRa5jy1QncKmC3nKkJgoArV67gypUrXO6EiMiCsVhqAazlMiye0Atrnu0NRxs5kjILMGZlHP6blvPwJ5PRbGxsEB8fj/j4eC53QkRkwdgNZyRzmMHbGFn5pZi9NQmpNwsBAC882hkLRj8ChZx1MxERtRyclLIJmfuYpQepqNJh2S8X8U18BgAgsKMzoqf1hmdrjqkhIqKWgWOWqE4KuRT/HOeLb57vA2dbK6TeLMSYVXH4+ewdsaM1K1VVVdi5cyd27tyJqqoqseMQEZGR2LLUQJbYsvRntwrKMGdrEpIyCwAAEf07YeGYnrCx4ozTDcXlToiIzBdblqjeOrjYYvvL/fHyoC4AgI0nbuDJL4/jep5a5GSWTyqVYtCgQRg0aBCkUv6qERFZKrYsGclSB3jX5cilXLy5IxX56go4WMux9Al/jA/0EDsWERFRo+MA7yZk6d1wf5VdWI6525KRcD0fADC1rxcWjfdltxwRETUr7IYjo7k722DrrFDMGdoNEgmwLSETE1XHcPVuidjRiIiIRMFiiWqQy6R4c0QPfPdCX7g6KHAxuxjjV8djb/JNsaNZlLKyMiiVSiiVSpSVccZ0IiJLxWKJahXm0xYH54ahf5c2KK3Q4o3tqXhrZypKK3gbfH3odDqkpqYiNTUVOp1O7DhERGQkFktUp3ZONtj8UiheD/eBRALsTLyJx6OPIT2nWOxoZs/GxgaHDh3CoUOHuNwJEZEF4wDvBmpuA7zrcvxqHl7/PgW5xRrYWEnxr8f98FRwR0gkErGjERERGYQDvMkkBnR1xcF5YQjzcUV5pQ5v7zqDqB2pUGvYLUdERM0XiyUjqVQq+Pr6IiQkROwoTcrVwRobZ/bFWyN7QCoB9ibfwvjoeFy4UyR2NLNTVVWFn376CT/99BOXOyEismDshmugltQN91cJGfmYuy0Z2UXlUMilWDy+F6b29WS33P/H5U6IiMwXu+GoSfTt3BoH54VhSI+2qKjS4Z29ZzH3+xQUl1eKHc0sSKVS9OnTB3369OFyJ0REFowtSw3UkluW/qDTCfgm/hqW/XIJVToB3m3sED2tN/w6OIsdjYiI6IHYskRNSiqV4G+PdcX2l/ujg4strv9eiie+OI7vTlwHa3EiIrJ0LJao0QR3aoWf5g7EcF83VGh1eO+H83htSxIKy9gtR0RElovFEjUqFzsFvn4uGO+N84WVTIKfz2Vj3Oo4pGYViB2tyZWVleHRRx/Fo48+yuVOiIgsWIsqliZNmoRWrVph8uTJ1bZfunRJv4aXUqmEra0t9u3bJ07IZkAikeCFgZ2x65UB8Gxti6z8Mkxecxzfxme0qG45nU6H48eP4/jx41zuhIjIgrWoAd6xsbEoLi7Gxo0bsWvXrgceU1JSAm9vb9y4caNet3pzgHfdCssqsWD3Gfx8LhsAEN7TDZ8+FQAXO4XIyUyvqqoKBw4cAACMGzcOcrlc5ERERPQHDvCuxeDBg+Ho6FjnMfv378ewYcM4J04jcba1whfTe+P9x3tBIZPi8IUcjF0Vj8Qb98SOZnJyuRwTJ07ExIkTWSgREVkwsymWjh49ivHjx8PDwwMSieSB3WAqlQre3t6wsbFBaGgoEhISGj3Hjh07MGXKlEY/b0smkUjwXH9v7HltALzb2OFWQRmmfHUCX/16FTpdi2nYJCIiC2U2xZJarUZgYCBUKtUD92/fvh1RUVFYtGgRkpKSEBgYiJEjRyI3N1d/jFKphJ+fX43H7du365WhqKgIx48fx5gxY2o9RqPRoKioqNqD6sevgzN+nDMQ4wM9UKUTsPTni3hx4ynkqyvEjmYSWq0WsbGxiI2NhVarFTsOEREZySzHLEkkEuzduxcTJ07UbwsNDUVISAiio6MB3B886+npiTlz5mDBggX1PndsbCyio6MfOGZp06ZN+M9//oPNmzfX+vzFixdjyZIlNbZzzFL9CYKA709lYfH+89BU6eDuZINVU4PQt3NrsaM1Ki53QkRkvprdmKWKigokJiYiPDxcv00qlSI8PBwnTpxotNepTxfcwoULUVhYqH9kZWU12uu3FBKJBFP7emFf5KPo0tYe2UXlmLr2JFRHrjSrbjmJRAJfX1/4+vpyvTwiIgtmEcVSXl4etFot3Nzcqm13c3NDdnZ2vc8THh6Op556CgcPHkTHjh2rFVqFhYVISEjAyJEj6zyHtbU1nJycsGnTJvTr1w/Dhg0z7JshvZ7tnfDj7IF4IqgDtDoBy/9zCRHrE5BXohE7WqOws7PD+fPncf78edjZ2Ykdh4iIjGQRxVJjOXz4MO7evYvS0lLcvHkT/fv31+9zdnZGTk4OFIr63dIeGRmJtLQ0nDp1ylRxWwR7azn+/XQglk0OgI2VFHGX8zB6ZRyOX80TOxoREREACymWXF1dIZPJkJOTU217Tk4O3N3dRcmkUqng6+uLkJAQUV6/OZFIJHi6jyd+nD0QPu0ccLdYg2e/+Q0rDqdD24y65YiIyDJZRLGkUCgQHByMmJgY/TadToeYmJhqrUNNiS1Ljc/HzRH7Zw/E0306QicAKw5fxnPf/obconKxoxmlrKwMw4cPx/Dhw7ncCRGRBTObmfJKSkpw5coV/dcZGRlISUlB69at4eXlhaioKERERKBPnz7o27cvVqxYAbVajZkzZ4qYmhqbrUKGZZMD0b9rG/xj7zkcv/o7xqyKw+dTlAjzaSt2PIPodDocPnxY/28iIrJMZlMsnT59GkOGDNF/HRUVBQCIiIjAhg0bMGXKFNy9exfvvfcesrOzoVQq8csvv9QY9N1UVCoVVCoV588xkUlBHRHQ0QWRW5JwMbsYz69LQOTgbng93AdymUU0iMLa2lo/DYW1tbXIaYiIyFhmOc+SJeHacKZVXqnFvw6kYetvmQCAvt6tsXKqEu2dbUVORkRElqzZzbNELZeNlQwfTfLHqqlBcLCWI+F6PsasjMORS7kPfzIREVEjYLFkJN4N17QmBHrgwJyB6OXhhHullZi5/hSW/nwBlVrzHQuk1Wpx6tQpnDp1it21REQWjN1wDcRuuKZVXqnFRwcv4LsTNwAAvb1csHpab3RwMb9uOS53QkRkvtgNR82WjZUM/3rcD19O7w1HGzmSMgswZmUc/puW8/AnNzGJRIJOnTqhU6dOXO6EiMiCsWXJSH++Gy49PZ0tSyLI/L0Uc7YlIfVmIQDgxYGd8fdRj0Ah52cAIiKqmyEtSyyWGojdcOKqqNLh458vYt2xDABAoKcLoqcGwbM112IjIqLasRuOWgyFXIr3xvvi6+eC4WQjR2pWAcasisMv5+6IHY2IiJoJFkvULIzo5Y6D88IQ5OWC4vIqvLI5CYt+OAdNlXh3oZWXl2PixImYOHEiysstc8kWIiJiN5zROGbJPFVqdfj00CV89es1AIBfBydET+0Nb9emvxONd8MREZkvjllqQhyzZJ6OXMxF1I4U3CuthIO1HB8/6Y9xAR5NmqGyshIbNmwAAMyYMQNWVlZN+vpERFQ7FktNiMWS+bpTWIa525Jx6vo9AMC0UC+8N84XNlYykZMREZHYOMCbCEB7Z1tsm9UPkUO6QiIBtv6WiYmqY7h6t0TsaEREZEFYLFGzJpdJ8dbIR/DdC33Rxl6Bi9nFGL86HnuTb5r8tXU6Hc6fP4/z589DpzPfZVmIiKhuLJaMxLXhLEuYT1v8PC8M/bu0QWmFFm9sT8Xbu1JRVmG6u+XKysrg5+cHPz8/lJWVmex1iIjItDhmqYE4ZsmyaHUCVsVcxqr/XYYgAD7tHPDF9N7wcXNs9NdSq9Xw9vYGAFy/fp13wxERmREO8G5CLJYs0/EreZi3PQV3izWwsZLiX4/74angjlzDjYioheAAb6KHGNDNFQfnhiHMxxXllTq8vesM3tyRCrWmSuxoRERkZlgsUYvV1tEaG2f2xfwR3SGVAHuSb2F8dDwu3CkSOxoREZkRFkvUokmlEswe6oNts/rBzcka1+6qMVF1DFt/y0RDe6jLy8sxffp0TJ8+ncudEBFZMI5ZaiCOWWo+8tUViNqRgthLdwEA4wM98NEkPzjaGDfzNpc7ISIyX4b8/ZY3USYis9faXoF1ESFYG3cNy/5zCT+m3sbZmwWIntYbfh2cDT6fQqHA559/rv83ERFZJrYsGYkL6TZviTfuYc7WJNwuLIdCJsU/x/XEc/068W45IqJmglMHNCF2wzVfBaUVmL/zDA5fyAEAjPF3x9InAuBsywVxiYgsnUmKpf379xscZPjw4bC1tTX4eZaExVLzJggCvo3PwCe/XESlVoBna1tET+2NQE+Xhz5Xp9MhMzMTAODl5QWplPdTEBGZC5MUS4a+0UskEly+fBldunQx6HmWhsVSy5CSVYDZW5Nw814ZrGQSLBjdEy886l1ntxwHeBMRmS+TTUqZnZ0NnU5Xr4ednV2Dvgkic6L0dMFPc8Mwqpc7KrUC3j+Qhr9tSkRBaUWdz7Ozs+PvAhGRhat3sRQREWFQl9qzzz7LlhZqVpxtrfDls72xZEIvKGRS/DctB2NXxSMp894Dj7e3t4darYZarWarEhGRBeMA7wZiN1zLdO5WISK3JuHG76WQSyV4e1QPvDSwC6RS3i1HRGQJuDZcLSZNmoRWrVph8uTJNfZ9/vnn6NWrF3x9fTF37twGz95MzZtfB2ccmDMQ4wLao0on4KODF/HSd6eRr667W46IiCyP0cVSZWUlsrKycOnSJeTn5zdmJpOZN28evvvuuxrb7969i+joaCQmJuLs2bNITEzEyZMnRUhIlsTRxgqrpwbho0n+UMil+N/FXIxdFYdT1+//Pmg0GsyaNQuzZs2CRqMROS0RERnLoGKpuLgYX375JQYNGgQnJyd4e3ujZ8+eaNu2LTp16oRZs2bh1KlTpsraYIMHD4ajo+MD91VVVaG8vByVlZWorKxEu3btmjgdWSKJRIJpoV74IfJRdGlrjzuF5Xjm65NQHbmCiopKfPPNN/jmm29QVVUldlQiIjJSvYulzz77DN7e3li/fj3Cw8Oxb98+pKSkID09HSdOnMCiRYtQVVWFESNGYNSoUbh8+bJBQY4ePYrx48fDw8MDEokE+/btq3GMSqWCt7c3bGxsEBoaioSEBINeozZt27bF/Pnz4eXlBQ8PD4SHh6Nr164GnUOtVlfruquoqIBara7RovDHgF+dTqffVllZCbVaXWOxVUOOLS0thVqthlar1W+rqqqCWq1GWVmZ0ceWlZVBrVZX+2Ov1WoNPra0tLTaseXl5VCr1aisrDTqWJ1Op///82cajQZqtRoVFRVGHSsIgv7YB13P2o59xN0RP84eiCeCOqCqsgKf/JiKlzYmYOF7i/HBBx/AysrqgdezMX5OHnQ9G+Pn5I/r2dCfk79ez4b+nNR2PRv6c/Ln62nIsYb83vM9ouW+Rxh77fkeUfexDX2PqDehnp555hnh3LlzDz2urKxM+PLLL4Vvv/22vqcWBEEQDh48KPzjH/8Q9uzZIwAQ9u7dW23/999/LygUCmHdunXC+fPnhVmzZgkuLi5CTk6O/pjAwEChV69eNR63bt3SH3PkyBHhySefrHbu/Px8YcSIEcLvv/8ulJaWCoMGDRJ+/fXXB+YsLy8XCgsL9Y+srCwBgABAyM3N1R/3wQcfCACEl156qdrz7ezsBABCRkaGftvnn38uABCmTZtW7VhXV1cBQLX/719//bUAQHj88cerHdupUycBgJCQkKDftnnzZgGAEB4eXu1YX19fAYBw5MgR/ba9e/cKAIQBAwZUO7ZPnz4CAOHAgQP6bYcOHRIACIGBgdWOHTRokABA2LFjh35bfHy8AEDo1q1btWPHjBkjABDWr1+v35acnCwAEDw8PKodO3nyZAGAEB0drd+Wnp4uABCcnZ2rHRsRESEAEJYtW6bfdvPmTQGAIJfLqx372muvCQCERYsW6bfdu3dPfz0rKir02+fPny8AEObPn6/fVlFRoT/23r17giAIgk6nEybPel0AIDgEjRVCPvivcPxKniAIgiCXywUAws2bN/XnWLZsmQBAiIiIqJbN2dlZACCkp6frt0VHRwsAhMmTJ1c71sPDQwAgJCcn67etX79eACCMGTOm2rHdunUTAAjx8fH6bTt27BAACIMGDap2bGBgoABAOHTokH7bgQMHBABCnz59qh07YMCAGr+3R44cEQAIvr6+1Y4NDw8XAAibN2/Wb0tISBAACJ06dap27OOPPy4AEL7++mv9tnPnzgkABFdX12rHTps2TQAgfP755/ptGRkZAgDBzs6u2rEvvfSSAED44IMP9Ntyc3P11/PP5s2bJwAQ3nnnHf22kpIS/bElJSX67e+8844AQJg3b161c/A94j6+R9y3aNEiAYDw2muvVXs9vkfc11TvEYWFhQIAobCwUHiYei+ku23btnodZ2Njg1deeaW+p9UbPXo0Ro8eXev+zz77DLNmzcLMmTMBAGvWrMFPP/2EdevWYcGCBQCAlJQUg18XAA4fPoxu3bqhdevWAICxY8fi5MmTeOyxx2ocu3TpUixZssSo16HmTyKRoJeHM3YBcLG1Qm6xBtO/OYl5w7qLHY2IiIzUqFMHZGVlYdGiRVi3bl2DziORSLB3715MnDgRwP1mSDs7O+zatUu/Dbg/91NBQQF++OGHep87NjYW0dHR2LVrl37byZMn8eqrr+LEiROwsrLChAkT8Le//Q2PP/54jedrNJpqzaFFRUXw9PTE7du34e7urp/RuaKiApWVlZDL5bC2ttYf/0ezoK2trX5W9MrKSlRUVEAmk8HGxsaoY0tLSyEIAmxsbCCTyQDcbw7VaDSQSqXV5sgy5NiysjLodDpYW1tDLr9fW2u1WpSXlxt0rEQiqTY5Y3l5ObRaLRQKBaysrAw+VqfT6Zt5/zyHkUajQVVVFaysrKBQKAw+VhAEfTOvnZ1djetpyLEVWuDdnQn4IfUOpLZOCPW0xyeTA9CpXSv99WyMn5MHXc/G+Dn543o29Ofkr9ezoT8ntV3Phv6c/Pl6NvTnpLbracixfI9o/u8R9bn2fI8w3XuEaAvppqamonfv3tX6L43x12Lp9u3b6NChA44fP47+/fvrj3v77bfx66+/4rfffqvXecPDw5Gamgq1Wo3WrVtj586d+vP94x//wJ49eyCVSjFs2DCsXLmyzqUsVCoVVCoVtFot0tPTOc8S1fDn5U56/H0vymEFVwcFVkwJwkAfV5HTERG1bIYUS/XuhgMevpjutWvXDDldkzt8+HCt+z788EN8+OGH9T5XZGQkIiMj9f+zieqy4+X+ePuHS7iYXYzn1v2GyMHd8Hq4D+SyFjXVGRGRRTKoWJo4cSIkEkmdEzbW1RpjLFdXV8hkMuTk5FTbnpOTA3d390Z/vfr4c8sS0YPY29tX+13ZF+mKJT+mYVtCJqKPXEFCRj5WTQ2Cu7NNHWchIiKxGfSxtn379tizZ0+ti+cmJSWZJKRCoUBwcDBiYmL023Q6HWJiYqp1yzWlyMhIpKWlmfW8UmRebKxkWPqEP1ZNDYK9QoaE6/kYsyoORy7lih2NiIjqYFCxFBwcjMTExFr3P6zVqS4lJSVISUnR39GWkZGBlJQUZGZmAgCioqKwdu1abNy4ERcuXMCrr74KtVqtvzuOyFJMCPTAgblh6OXhhHx1BWauP4WlP19ApVb38CcTEVGTM2iAd1xcHNRqNUaNGvXA/Wq1GqdPn8agQYMMDhIbG4shQ4bU2B4REYENGzYAAKKjo7F8+XJkZ2dDqVRi1apVCA0NNfi1GgMHeNPDaDQa/P3vfwcAfPLJJ9XuZAGA8kotPjp4Ad+duAEACO7UCqumBqGDi22NcxERUeMS7W64lsiQ/9nUsvz5briSkpJqt7D+2cGzd/D3XWdQrKmCi50VPp0ciHBft6aMSkTU4pjsbjgiqj8rKyu88847+n/XZox/e/h5OGP2tiScuVmIl747jZcGdsbbox6BQs675YiIxMaWJSOxG44aW0WVDh//fBHrjmUAAAI9XRA9NQiere0e8kwiIjIUu+GaELvhqLEdOp+N+TtTUVReBUcbOZZPDsQoP3GmyCAiaq4M+fttUBt/dnZ2g4IRtSRCLSuOP8yIXu44OC8MQV4uKC6vwiubE7F4/3loqjinFxGRGAwqlkaMGGGqHBZHpVLB19cXISEhYkchM1VaWgoHBwc4ODjo14eqr46t7LDj5f54+bEuAIANx69j8pcncON3tSmiEhFRHQzqhvP398fZs2dNmcfisBuOalPfu+Ee5n8Xc/DmjlTcK62Eg7UcHz/pj3EBHo0ZlYioxTFZN5wpljIhaq7s7OxQUlKCkpKSaqtkG2roI244OC8MId6tUKKpwuytyfjH3rMor2S3HBFRU+B9yUQmIpFIYG9vD3t7+wZ/0GjvbItts/ohckhXSCTAlt8yMemL47h2t6SR0hIRUW1YLBmJY5aoqcllUrw18hFsnNkXbewVuHCnCONWx2Nf8i2xoxERNWsGjVkKCgpCcnKyKfNYHI5ZotpUVFRgyZIlAIBFixZBoVA02rlzi8ox9/tknLyWDwCY0scTiyf0gq1C1mivQUTUnHGepSbEYolq01gDvGuj1QlYFXMZq/53GYIA9HBzhGp6ELq1c2zU1yEiao643AmRGZDL5Zg3b57+341NJpXgjeHdEdq5NeZtT8GlnGKMX30M70/0w+Tgjo3+ekRELZXBLUtarRY//vgjhg0bBkdHfoJlyxKZg7vFGkTtSEHc5TwAwBO9O+CDiX6wU/DzEBHRg5hs6gAAkMlkmDp1Ku7evWt0wOaAA7zJnLR1tMbGmX0xf0R3SCXAnqRbGL86Hhezi8SORkRk8Yy6Gy4kJAQZGRmNncWiREZGIi0tDadOnRI7ChEAQCqVYPZQH2yb1Q9uTta4eleNx6OP4fuETIOWWyEiouqMKpbmzJmDd955B1lZWY2dh6jZUKvVkEgkkEgkUKubbpmS0C5tcHBuGAZ1bwtNlQ4L9pzFvO9TUKKparIMRETNiVF3w0ml92ssBwcHTJgwAYMHD0ZQUBD8/f0b9fZoS8AxS1QbU98N9zA6nYCv465h+X8uQasT0NnVHtHTgtDLw7lJcxARmSOTTx1w/fp1nDlzBikpKUhNTUVKSgquX78OuVyOHj164MyZM0aHtzQslqg2giAgL+/+gGtXV1fRlgtKvJGPOVuTcbuwHAq5FO+O88WzoV5cvoiIWjRR5lkqLi5GSkoKzpw5g8jIyMY4pUVgsUSW4J66Am/tSsXhC7kAgLH+7bH0SX842ViJnIyISBwmvRsOAF544QVs2LBB//WNGzcQHx+PgICAFlUoEVmKVvYKrH2+D/45tifkUgl+OnsH41bF48zNArGjERGZPaOKpYMHD+KRRx4BABQUFCA4OBgTJ06Er68v0tPTGzUgkaWqqKjAhx9+iA8//BAVFRVix4FEIsFLYV2w85X+6OBii8z8Ujz55XGsi8/g3XJERHUwqlgqLCxEhw4dAAC7d++Gu7s7ioqKMGXKFCxYsKBRA5orzrNED1NZWYl//vOf+Oc//4nKykqx4+gFebXCwblhGNnLDZVaAf86kIaXNyWisNR8MhIRmROjiiVPT0/9PEs7d+7EjBkzYG1tjVdeeQXHjh1r1IDmivMs0cPI5XK89NJLeOmll0yy3ElDONtZYc2zwVg83hcKmRSH0nIwZlUckjLviR2NiMjsGDXA+6OPPsKOHTswfvx4fPzxx7hw4QK6deuGixcvIjg4uEnnlBEbB3iTpTt7sxCRW5OQmV8KuVSCt0f1wEsDu0Aq5d1yRNR8mXyA98KFC/HUU0/h6NGj+Pjjj9GtWzcAwKlTp+Dl5WXMKYlIJP4dnXFg7kCMDWiPKp2Ajw5exEvfncY9tfjjrIiIzEGjTR0AAMuXL0d5eTnefffdxjql2WPLEjUXgiBga0ImlvyYhooqHdo722DV1CCEeLcWOxoRUaMTZZ6llorFEtVGrVajXbt2AIDc3Nwmn8HbWGm3izB7axKu5akhk0oQNbw7Xh3Uld1yRNSsmKQbLjMz06AQt27dMuh4ouaotLQUpaWlYscwiK+HE/bPGYiJSg9odQKW/+cSItYnIK9EI3Y0IiJR1LtYCgkJwcsvv1zn3V+FhYVYu3Yt/Pz8sHv37kYJ2JgmTZqEVq1aYfLkyTX2ffrpp+jVqxf8/PywefNmEdJRc2Nra4uMjAxkZGTA1tZW7DgGcbCW4/MpSix7MgA2VlLEXc7DmJVxOHH1d7GjERE1uXp3w/3+++/48MMPsW7dOtjY2CA4OBgeHh6wsbHBvXv3kJaWhvPnz6N379549913MWbMGFNnN1hsbCyKi4uxceNG7Nq1S7/97NmziIiIwPHjxyEIAoYMGYJffvkFLi4uDz0nu+GoubuUXYzIrUm4klsCqQSYN6w7Zg/tBhm75YjIgpmkG65Nmzb47LPPcOfOHURHR8PHxwd5eXm4fPkyAGD69OlITEzEiRMnzLJQAoDBgwfD0dGxxvYLFy6gf//+sLGxga2tLQIDA/HLL7+IkJDI/PRwd8T+2Y9icnBH6ATg88PpeO7b35BbXC52NCKiJmHw1AG2traYPHkyVqxYgb179+KXX37B5s2b8eabb8LPz8/oIEePHsX48ePh4eEBiUSCffv21ThGpVLB29sbNjY2CA0NRUJCgtGv92d+fn6IjY1FQUEB7t27h9jYWI65ogarrKzEihUrsGLFCrOawdsYdgo5Pn0qEP9+KhC2VjIcv/o7xqyMR/zlPLGjERGZnNlMK6xWqxEYGIgXXngBTzzxRI3927dvR1RUFNasWYPQ0FCsWLECI0eOxKVLl/R3HCmVSlRVVdV47qFDh+Dh4VHra/v6+mLu3LkYOnQonJ2d0a9fP8hksgceq9FooNH830DXoqIiQ79VaiEqKirwxhtvAABmzZoFKysrkRM13JPBHRHo6YLZW5NwMbsYz637DbOHdMO8YT6Qy4yato2IyOw16rtbVlYWXnjhBaOeO3r0aHzwwQeYNGnSA/d/9tlnmDVrFmbOnAlfX1+sWbMGdnZ2WLdunf6YlJQUnDt3rsajrkLpDy+//DKSkpJw5MgRWFlZwcfH54HHLV26FM7OzvqHp6enUd8vNX8ymQzTpk3DtGnTai2+LVG3dg7YF/kopvb1hCAAq/93BdO++Q3ZheyWI6LmqVGLpfz8fGzcuLExTwng/if0xMREhIeH67dJpVKEh4fjxIkTjfIaubm5AIBLly4hISEBI0eOfOBxCxcuRGFhof6RlZXVKK9PzY+NjQ22bNmCLVu2wMbGRuw4jcrGSoalTwRg5TNK2CtkSMjIx5hVcYi9lCt2NCKiRmdQN9z+/fvr3H/t2rUGhalNXl4etFot3Nzcqm13c3PDxYsX632e8PBwpKamQq1Wo2PHjti5cyf69+8PAHj88cdRWFgIe3t7rF+/vtaFT62trWFtbQ2VSgWVSgWtVmv8N0Zk4R5XdoB/B2fM3pqMtDtFmLH+FF4d3BVRw7vDit1yRNRMGFQsTZw4ERKJBHXNNiCRmO/txIcPH651n6EtVJGRkYiMjNTfekjUUnVp64A9rw3Ahz9dwKaTN/Bl7FUkZORj9dQgeLhY1vxSREQPYtBHv/bt22PPnj3Q6XQPfCQlJZkkpKurK2QyGXJycqptz8nJgbu7u0le82FUKhV8fX0REhIiyuuT+VOr1Wjbti3atm0LtVotdhyTsrGS4f2Jfvhiem84WsuReOMexqyKQ8yFnIc/mYjIzBlULAUHByMxMbHW/Q9rdTKWQqFAcHAwYmJi9Nt0Oh1iYmL03WhNLTIyEmlpaXXOaE6Ul5eHvLyWc3v9GP/2+GluGAI6OqOgtBIvbjyNDw7cX5iXiMhSGdQN99Zbb9X5Cblbt244cuSIUUFKSkpw5coV/dcZGRlISUlB69at4eXlhaioKERERKBPnz7o27cvVqxYAbVajZkzZxr1ekSmZmtri3Pnzun/3VJ4tbHDzlf645OfL2HdsQx8E5+B0zfuYfXUIHi2thM7HhGRweq93ImpxcbGYsiQITW2R0REYMOGDQCA6OhoLF++HNnZ2VAqlVi1ahVCQ0ObOOl9fx7gnZ6ezuVOiB7g0PlszN+ZiqLyKjjZyLH8qUCM7CVO1zkR0Z8ZstyJ2RRLloprwxHV7ea9UszZlozkzAIAwIwB3lg45hFYy5vP3FNEZHlMsjYcERmmsrISa9euxdq1ay1+uZOG6NjKDjte7o+/PdYFALDh+HVM/vIEbvzevAe9E1HzYVDLUnZ2tmh3n5kbdsPRw6jVajg4OAC4PybP3t5e5ETi+9/FHLy5IxX3SivhaC3Hx08GYGxAe7FjEVELZLJuuICAAJw5c6bBAZsTdsNRbcrLy/HMM88AAL7//vtmN4u3se4UlmHutmScun4PAPBsPy/8c6wvbKzYLUdETcdkxZK/vz/Onj3b4IDNCYslIsNVaXX4/HA6VEeuAgB6tneCaloQurR1EDkZEbUUJhuzZM6zczc1TkpJZDy5TIq3Rj6CjS/0RRt7BS7cKcL41fH4IeWW2NGIiGpgN1wDsWWJqGFyisox7/tknLyWDwB4JsQTi8b3gq2C3XJEZDq8G47IDJSWlsLb2xve3t4oLS0VO47ZcnOywZaX+mHuMB9IJMD3p7IwUXUMV3KLxY5GRATAwGJJJuMnPaL6EgQBN27cwI0bN0yyDFBzIpNKEDW8Oza/GApXB2tcyinG+NXHsCvxptjRiIgMK5aSk5NNlcPicMwSPYyNjQ0SEhKQkJDAO+Hq6dFurjg4byAe7dYGZZVazN+Zijd3pKK0okrsaETUgnEG7wbimCWixqfVCfjiyBV8fjgdOgHo2tYeX0wPRg93R7GjEVEzwTFLRGTRZFIJ5gzzwbZZ/eDmZI2rd9WYEB2P7xMy2aVJRE2OxRKRiVRVVWHLli3YsmULqqrYjWSM0C5tcHBuGAZ1bwtNlQ4L9pzF69tTUKLh/08iajpGdcO98MILeOyxxzBjxgwAwI0bN5CWloYBAwbA2dm5sTOaNXbDUW243Enj0ekEfHX0Gj49dAlanYDOrvaInhaEXh4t6/2GiBqPybvhDh48iEceeQQAUFBQgODgYEycOBG+vr64dOmSMae0OBzgTQ8jlUoRHh6O8PBwSKVsxG0IqVSCVwd3xfa/9UN7Zxtk5Kkx6Yvj2HSSdxoSkekZ1bJka2uL9PR0eHp64ttvv8Xnn3+OxMRELFy4ENevX8eePXtMkdUssWWJqGndU1dg/s5UxFzMBQCM9W+PpU/6w8nGSuRkRGRJTN6y5OnpiYyMDADAzp07MWPGDFhbW+OVV17BsWPHjDklEVG9tLJX4JuIPvjn2J6QSyX46ewdjFsVjzM3C8SORkTNlFHF0owZMzB37ly8++67iImJwcSJEwEAOp0OJSUljZmPiKgGiUSCl8K6YOcr/dHBxRaZ+aV48svjWH8sg91yRNTojCqWFi5ciKeeegpHjx7Fxx9/jG7dugEATp06BS8vr0YNSGSpSktL0atXL/Tq1YvLnZhIkFcrHJwbhhG+bqjUCljyYxpe2ZyIwtJKsaMRUTPSqJNSLl++HOXl5Xj33Xcb65Rmj2OWqDa8G67pCIKADcev46ODF1CpFdDBxRbR04IQ5NVK7GhEZKYM+fvNGbwbiMUS1Uar1SIuLg4AEBYWxrUVm8CZmwWYvTUZmfmlkEsl+PuoR/BSWGdIJBKxoxGRmTHpAO+ysjLEx8cjLS2txr7y8nJ89913hp6SqFmSyWQYPHgwBg8ezEKpiQR0dMGBuQMx1r89qnQCPjx4AS9tPI176gqxoxGRBTOoWEpPT0fPnj3x2GOPwd/fH4MGDcKdO3f0+wsLCzFz5sxGD2mOOM8SkXlysrFC9LQgvD/RDwq5FDEXczF2VRxOX88XOxoRWSiDiqW///3v8PPzQ25uLi5dugRHR0c8+uijyMzMNFU+sxUZGYm0tDScOnVK7ChkpqqqqrBv3z7s27ePy500MYlEguf6dcLe1wags6s9bheWY8rXJ/FF7BXodBx5QESGMWjMkpubGw4fPgx/f38A9wdVvvbaazh48CCOHDkCe3t7eHh4QKvVmiywueGYJaoNB3ibhxJNFf6x9yx+SLkNABjUvS0+ezoQbRysRU5GRGIy2ZilsrIyyOVy/dcSiQRffvklxo8fj0GDBiE9Pd24xETNkFQqxYABAzBgwAAudyIiB2s5VkxR4pMn/WEtl+LX9LsYsyoOv137XexoRGQh5A8/5P888sgjOH36NHr27Flte3R0NABgwoQJjZeMyMLZ2tpyRnszIZFIMCXEC0rPVnhtSyKu3lVj6tqTeCO8O14b0g0yKe+WI6LaGfRxd9KkSdi2bdsD90VHR2Pq1KmcPZeIzFYPd0f8OGcgnuzdEToB+Pd/0/H8ut9wt1gjdjQiMmOcZ6mBOGaJyDLtSryJd/edQ1mlFq4O1lj5jBKPdnMVOxYRNRGTL6RribKysjB48GD4+voiICAAO3furLb/wIED6NGjB3x8fPDNN9+IlJKak7KyMoSEhCAkJARlZWVix6G/mBzcET/OeRQ93ByRV6LBs9/+hs/+mw4t75Yjor9oMS1Ld+7cQU5ODpRKJbKzsxEcHIz09HTY29ujqqoKvr6+OHLkCJydnREcHIzjx4+jTZs2Dz0vW5aoNrwbzjKUVWix5Mfz+P5UFgAgtHNrrJoaBDcnG5GTEZEpmaxlKTs7u0HBxNS+fXsolUoAgLu7O1xdXZGff3+SuoSEBPTq1QsdOnSAg4MDRo8ejUOHDomYlpoDa2trHDhwAAcOHIC1NW9TN1e2Chk+fjIAK59Rwl4hw28Z+RizMg6/pt8VOxoRmQmDiqURI0aYKgeOHj2K8ePHw8PDAxKJBPv27atxjEqlgre3N2xsbBAaGoqEhASjXisxMRFarRaenp4AgNu3b6NDhw76/R06dMCtW7eMOjfRH+RyOcaOHYuxY8dWm3KDzNPjyg74cc5A9GzvhN/VFYhYl4BPfrmIKq1O7GhEJDKDiiVT9tip1WoEBgZCpVI9cP/27dsRFRWFRYsWISkpCYGBgRg5ciRyc3P1xyiVSvj5+dV43L59W39Mfn4+nn/+eXz99ddG5dRoNCgqKqr2IKLmoUtbB+x9bQCe7ecFAPgy9iqe+fokbhdwzBlRS2bQx11Trtw9evRojB49utb9n332GWbNmqVfe27NmjX46aefsG7dOixYsAAAkJKSUudraDQaTJw4EQsWLMCAAQP02z08PKq1JN26dQt9+/Z94DmWLl2KJUuW1PfbohZMq9Xif//7HwBg6NChXEzXQthYyfDBRH/069IGC3efxekb9zBmVRw+ezoQQx9xEzseEYnAIu6Gq6ioQGJiIsLDw/XbpFIpwsPDceLEiXqdQxAEzJgxA0OHDsVzzz1XbV/fvn1x7tw53Lp1CyUlJfj5558xcuTIB55n4cKFKCws1D+ysrKM/8aoWSsvL8eIESMwYsQIlJeXix2HDDQuwAMH5g6EfwdnFJRW4oUNp/HhT2moqGK3HFFLYxHFUl5eHrRaLdzcqn+qc3Nzq/eg82PHjmH79u3Yt28flEollEolzp49C+D+2JJ///vfGDJkCJRKJd58881a74SztraGk5MTNm3ahH79+mHYsGEN++ao2ZJKpQgMDERgYCCXO7FQndrYY9er/TFjgDcAYG1cBp7+6gSy8kvFDUZETcqgbjhL7kYYOHAgdLraPxFOmDDBoOVaIiMjERkZqb/1kOivbG1tH9o1TObPWi7D4gm90K9LG7y9KxUpWQUYuyoOy58KxMhe7mLHI6ImYNDH3eTkZFPlqJOrqytkMhlycnKqbc/JyYG7uzhvViqVCr6+vggJCRHl9YmoaY3yc8dPc8MQ6OmCovIqvLwpEYv3n4emSit2NCIyMYvoG1AoFAgODkZMTIx+m06nQ0xMDPr37y9KpsjISKSlpeHUqVOivD4RNT3P1nbY+XJ/zArrDADYcPw6Jn95Ajd+V4ucjIhMyWyKpZKSEqSkpOi7LTIyMpCSkoLMzEwAQFRUFNauXYuNGzfiwoULePXVV6FWq/V3xxGZm7KyMgwePBiDBw/mcifNiEIuxT/G+uLbiD5wsbPC2VuFGLcqHj+duSN2NCIyEbNZ7iQ2NhZDhgypsT0iIgIbNmwAAERHR2P58uXIzs6GUqnEqlWrEBoa2sRJ71OpVFCpVNBqtUhPT+dyJ1QDlztp/m4XlGHOtmQk3rgHAHi2nxf+OdYXNlaWO76TqKUwZLkTo4qloqIirF+/HtnZ2ejcuTMCAwPh7+8POzs7o0NbKq4NR7WpqqrC3r17AQCTJk3iLN7NVKVWh8/+m44vY68CAHq2d4JqWhC6tHUQORkR1cXkxVJ4eDhSU1MREhKCzMxMXLp0CQDQtWtXBAYGYvv27cYlt0AslogIAGIv5SJqRyry1RWwV8jw0RP+eFzZ4eFPJCJRGPL326iPuidOnEBsbKz+TjCNRoOzZ88iJSUFqampxpzS4vy5G46IaHCPdjg4Nwxzv09GQkY+5n2fghNXf8ei8b1gq2C3HJElM6plqX///vjiiy8QFBRkikwWhS1LVButVouTJ08CAPr162fR85RR/VVpdVgVcxmrj1yBIAA93Byhmh6Ebu0cxY5GRH9i8m64uLg4LFu2DLt27YK1tbXRQZsDFktUGw7wbtniL+fh9e0pyCvRwNZKhvcn+mFycEexYxHR/2fI32+jpg7w9vZGUVERfH198c4772D//v0tbo00TkpJDyORSNCtWzd069bNpItQk3ka6OOKg/MGYkDXNiir1GL+zlS8uSMVpRVVYkcjIgMZ1bLUt29f5OTkYNCgQcjMzERqaiqKiorQunVrBAUF4dChQ6bIapbYskREddHqBKiOXMGKw+nQCUC3dg5QTeuNHu7sliMSk8kHeJ87dw4nTpxAYGCgftv169eRnJyMM2fOGHNKIqJmSSaVYO4wH/Tt3BpztyXjSm4JJkTHY8mEXpgS4slWRyILYFTL0qBBg7B06VIMGDDAFJksCluWiKi+8ko0iNqRiqPpdwEAjys98OEkfzhYcw4uoqZm8jFL8+bNw+LFi1FQUGDM05sFjlmihykvL8fYsWMxduxYlJeXix2HzICrgzU2zAjB26N6QCaV4IeU25iwOh5pt4vEjkZEdTCqZUkqvV9jtWnTBpMmTUJoaCiCgoLg5+cHhULR6CHNGVuWqDa8G47qcvp6PuZsS8adwnIo5FK8N84X00O92C1H1ERMPnXAjRs3kJqaqp+EMiUlBdevX4dcLkePHj1a1LglFktUm8rKSmzZsgUAMH36dFhZWYmciMzNPXUF3tyZiv9dzAUAjA1oj6VP+MPJhj8rRKZm8mLpQYqLi5GSkoIzZ84gMjKyMU5pEVgsEVFDCIKAb+Iy8MkvF1GlE9CpjR2ip/aGf0dnsaMRNWsmL5aOHz8OJycn+Pn5GR2yuWCxRESNISnzHuZsTcatgjIoZFK8M+YRRAzwZrcckYmYfIB3ZGQkfvvttxrbr169iuLiYmNOaXE4wJseRqvVIiUlBSkpKVxDkB6qt1crHJwbhhG+bqjQ6rD4xzS8sjkRhaWVYkcjavGMalmys7PD2bNn0bVr12rbv/rqK/z44484cOBAowU0d2xZotpwgDcZQxAEbDh+HR8dvIBKrYCOrWwRPa03lJ4uYkcjalZM3rLk5OSEe/fu1dgeFhamXziUqKWTSCTw8PCAh4cHu1Ko3iQSCWY+2hm7Xx0Ar9Z2uHmvDJO/PI5v4q6hkYaYEpGBjCqWRo0ahU8//bTmyaRSVFRUNDgUUXNgZ2eHW7du4datW7CzsxM7DlmYgI4uODB3IMb6t0eVTsAHP13ArO9Oo6CU77FETc2oYun999/Hr7/+iieffBJnz54FcH8Cvk8++QQBAQGNGpCIqKVysrFC9LQgvD/RDwq5FIcv5GLMyjgk3sgXOxpRi2JUseTp6YmTJ0+irKwMgYGBsLW1haOjI3788UcsX768sTMSEbVYEokEz/XrhL2vDUBnV3vcLizH01+dxJexV6HTsVuOqCk0eJ6lzMxMpKSkwMrKCqGhoWjdunVjZbMIHOBNtSkvL8dzzz0HANi0aRNsbGxETkSWrkRThX/sPYsfUm4DAAZ1b4vPng5EGwdrkZMRWR5RJqVsqVgsUW14NxyZgiAI2H4qC4v2n4emSgc3J2useiYIoV3aiB2NyKKY/G642mRlZeGFF15ozFOaLc6zRA+jUCgQHR2N6OjoFrdmIpmORCLBM3298MPsR9G1rT1yijSYuvYkVsdchpbdckQm0agtS6mpqejdu3eLmoCPLUtEJJbSiiq8u+88difdBAAM7OaKz6co0daR3XJED2PI32+5ISfev39/nfuvXbtmyOmIiKgB7BRy/PvpQPTv2gbv7juH+Ct5GL0yDqueUWJAN1ex4xE1Gwa1LEmlUkgkkjonRpNIJGxZIgKg0+lw9epVAEDXrl0hlTZqrzdRNZdzihG5NQnpOSWQSIA5Q30wb5gPZFJOiEr0ICYbs9S+fXvs2bMHOp3ugY+kpKQGBSdqTsrKytC9e3d0794dZWVlYsehZs7HzRE/RA7EMyGeEARgVcxlTP/mJHKKysWORmTxDCqWgoODkZiYWOv+h7U6EbU0zs7OcHZ2FjsGtRC2Chk+fjIAK59Rwl4hw8lr+RizMg5H0++KHY3IohnUDRcXFwe1Wo1Ro0Y9cL9arcbp06cxaNCgRgto7tgNR0Tm6NrdEkRuTcaFO0UAgNcGd0XU8O6Qy9gdTASYsBsuLCys1kIJAOzt7c22UMrKysLgwYPh6+uLgIAA7Ny5s9r+SZMmoVWrVpg8ebJICYmIGk+Xtg7Y+9oATA/1AgB8EXsVU9eexO0CdgkTGarFTEp5584d5OTkQKlUIjs7G8HBwUhPT9dPFBgbG4vi4mJs3LgRu3btqvd52bJERObuwJnbWLD7LEo0VWhlZ4V/Px2IoY+4iR2LSFSiTUppztq3bw+lUgkAcHd3h6urK/Lz/28xysGDB8PR0VGkdNQcaTQazJgxAzNmzIBGoxE7DrVg4wI88NPcgfDr4IR7pZV4YcNpfHTwAiq1OrGjEVkEg4ql7OxsU+XA0aNHMX78eHh4eEAikWDfvn01jlGpVPD29oaNjQ1CQ0ORkJBg1GslJiZCq9XC09OzgamJaldVVYWNGzdi48aNqKqqEjsOtXCd2thj96sDMGOANwDg66PX8PRXJ3DzXqm4wYgsgEHF0ogRI0yVA2q1GoGBgVCpVA/cv337dkRFRWHRokVISkpCYGAgRo4cidzcXP0xSqUSfn5+NR63b9/WH5Ofn4/nn38eX3/9tVE5NRoNioqKqj2IHsTKygrLli3DsmXLYGVlJXYcIljLZVg8oRfWPNsbjjZyJGcWYMzKOBw6b7oPwkTNgUFjlvz9/XH27FlT5gFwfwqCvXv3YuLEifptoaGhCAkJQXR0NID7E/55enpizpw5WLBgQb3Oq9FoMHz4cMyaNUu/GvyfxcbGIjo6us4xS4sXL8aSJUtqbOeYJSKyJFn5pZi9LRmpWQUAgJmPemPh6J5QyFvM6Axq4Uw2ZkkiEWcm2IqKCiQmJiI8PFy/TSqVIjw8HCdOnKjXOQRBwIwZMzB06NAHFkr1tXDhQhQWFuofWVlZRp+LiEgsnq3tsPPl/nhpYGcAwPpj1zF5zXFk/s5uOaK/soiPEHl5edBqtXBzq373hpubW73HUR07dgzbt2/Hvn37oFQqoVQqq7WShYeH46mnnsLBgwfRsWPHWoswa2trODk5YdOmTejXrx+GDRtm/DdGzZpOp8OtW7dw69Yt6HQcSEvmRyGX4p/jfPHN833gbGuFMzcLMXZVHA6evSN2NCKzYtBCupZs4MCBdf7BOnz4sEHni4yMRGRkpL4Zj+ivysrK0LFjRwBASUmJfpoKInMT7uuGg/PCMHdbMhJv3MNrW5LwXL9O+MfYnrCxkokdj0h0BrUsyWTi/NK4urpCJpMhJyen2vacnBy4u7uLkkmlUsHX1xchISGivD5ZBrlcDrm8xXwmIQvWwcUW3/+tH14Z1BUAsOnkDTz55XFk5KlFTkYkPoOKpeTkZFPlqJNCoUBwcDBiYmL023Q6HWJiYtC/f39RMkVGRiItLQ2nTp0S5fXJ/Nnb26OyshKVlZVsVSKLYCWTYsHoR7BhZgha2ytw/nYRxq2Kw/7U2w9/MlEzZjZjlkpKSpCSkoKUlBQAQEZGBlJSUpCZmQkAiIqKwtq1a7Fx40ZcuHABr776KtRqNWbOnCliaiKi5mdwj3Y4ODcMfTu3hrpCi7nbkrFwz1mUV2rFjkYkCoOmDrh37x5atWplkiCxsbEYMmRIje0RERHYsGEDACA6OhrLly9HdnY2lEolVq1ahdDQUJPkeRiVSgWVSgWtVov09HROHUBEzU6VVoeVMZcRfeQKBAF4xN0R0dN6o1s7B7GjETWYIVMHGFQstW3bFu+//z5efvll0aYRMDdcG45qo9FoEBUVBQD47LPPYG1tLXIiIuPEX87D69tTkFeigZ1Chg8m+uGJ3h3FjkXUICYrlpYuXYqlS5eiS5cuWL16NcLCwhoc1tKxWKLaqNVqODjc/wTOu+HI0uUWl+P171Nw/OrvAICngjtiyeO9YKfgDQxkmUw2KeXChQtx6dIlBAUFYciQIZg6dSpu3brVoLCWinfD0cNYWVlh0aJFWLRoEZc7IYvXztEGm14MxRvh3SGVADsTb+Lx6GNIzykWOxqRyRnUsvRniYmJeP3115GcnIwFCxbgrbfeapHdDGxZIqKW5sTV3zHv+2TkFmtgYyXFvyb44ak+HTk8gyyKyVqW/iw4OBhxcXH49ttv8e2336Jnz57Yu3evsacjIiIL0b9rGxycF4YwH1eUV+rw9u4ziNqRCrWmSuxoRCbR4KkDpkyZgosXL+LFF19EREQEhg8f3hi5zB674ehhBEFAQUEBCgoKYGQDLpHZcnWwxsaZffHWyB6QSSXYm3wL41fH48KdIrGjETU6o7vhKioqcPHiRZw7d07/+O233/TruLUU7Iaj2nCAN7UUp67nY87WZGQXlUMhl2LReF9M6+vFbjkya4b8/TboNoYlS5boC6OrV6+iqqoKzs7O8PPzQ0BAAMaMGYOAgIAGhSciIssS4t0aB+eFYf7OVPzvYi7+sfccTlz9HUuf8IejDW9uIMtnUMuSn58f/P39ERAQoP+vl5eXKfOZPbYsUW0EQUBV1f0xHHK5nJ+yqdnT6QR8E38Ny365hCqdgE5t7KCa1ht+HbjYOJkfk82zRP+HM3gTET1YUuY9zNmajFsFZVDIpPjH2J54vn8nfmAgs8JiqQmxZYmIqKaC0gq8tesM/puWAwAY1csdn0wOgLMtu+XIPJh86oDCwkL87W9/Q7du3dCzZ0/cuXPHqKBEzVlFRQXeeustvPXWW6ioqBA7DlGTcrFT4OvngvHeOF9YyST45Xw2xq6KQ0pWgdjRiAxmVLEUGRmJs2fPYtmyZbhx4wbKysoAAG+88Qaio6MbNSCRpaqsrMSnn36KTz/9FJWVlWLHIWpyEokELwzsjF2vDIBna1vcvFeGp9Ycxzdx1zidBlkUo4qln3/+GV988QWeeOIJyGQy/faRI0di48aNjRaOyJJZWVlh/vz5mD9/Ppc7oRYt0NMFB+aEYbSfOyq1Aj746QJmfXcaBaVscSXLYFSxJAgCHB0da2z38fHB5cuXGxzKEnBSSnoYhUKB5cuXY/ny5VAoFGLHIRKVs60VvpjeG+8/3gsKmRSHL+RizMo4JN7IFzsa0UMZVSyNHj0aW7ZsqbFdrVa3mLsdIiMjkZaWhlOnTokdhYjIIkgkEjzX3xt7XhsA7zZ2uF1Yjqe/Ook1v16FTsduOTJfBk1K+YelS5eiT58+AO63MkkkEpSXl+P9999H7969GzUgkaXiPEtED+bXwRk/zhmId/aew4+pt/Hxzxdx8trv+PdTgWjj0PIWZCfzZ1TLkpeXF44fP47jx4+jtLQUffv2hYuLC3799Vd88sknjZ2RyCKVlpZCoVBAoVCgtLRU7DhEZsXRxgqrnlFi6RP+sJZLEXvpLsasikNCBrvlyPw0eJ6lzMxMpKamwsrKCqGhoWjVqlVjZbMInGeJasO14Yjq58KdIkRuTcK1u2pIJUDU8O54bXA3SKVsjSXT4aSUTYjFEtVGEAQUFhYCAJydndkNR1QHtaYK7+47hz3JtwAAYT6u+OxpJdo6sluOTMNkk1Lu3r0bSqVS//WCBQuwbt06JCYmQqPRGBWWqLmSSCRwcXGBi4sLCyWih7C3luOzKUosnxwAGysp4i7nYcyqOBy/kid2NCLDWpbGjRuH8PBwvP766wAAR0dHaLValJeXQyaToWfPnjh69ChcXFxMFNf8sGWJiKhxXc4pRuTWJKTnlEAiAeYO9cHcYT6QsVuOGpHJWpbOnz+PESNGVNt29uxZXLt2DXv27IGVlRXWrFljeGILxHmW6GEqKiqwePFiLF68mMudEBnAx80RP0QOxJQ+nhAEYGXMZTz7zW/ILSoXOxq1UAa1LNnY2ODq1avo0KEDAMDFxQVJSUno0qULAOD777/H6tWrcezYMdOkNUNsWaLacIA3UcPtS76Fd/aeRWmFFm3sFfh8ihKPdW8rdixqBgz5+23QPEuurq64fv26vljKzs6utoyDUqlEWlqaEZGJmh+5XI7XXntN/28iMtzEoA7w7+iMyC1JuJhdjIj1CXhtcFe8Ed4dcplRs98QGcyglqXnn38ecrkc69ate+D+9PR09O7dGyUlJY0W0NyxZYmIyPTKK7V4/0AatvyWCQDo690aK6cq0d7ZVuRkZKlMNmbprbfewtatW7Fy5coH7j927Ji+S46IiKix2FjJ8OEkf6yeGgQHazkSrudjzMo4HLmYK3Y0agEMKpb8/f2xefNmvPXWWwgPD8fu3buRmZmJ27dvY8eOHVi4cCGmT59uqqxERNTCjQ/0wIE5A+HXwQn3Sisxc8MpLD14AZVandjRqBkzalLK5ORkvPHGGzh69Kh+/hhBEDB+/Hjs2rWr2jim5o7dcFQbtVqtn0ajoKCAA7yJGpGmSoulBy9iw/HrAIAgLxesnhqEjq3sxA1GFsNk3XB/CAoKQmxsLK5fv479+/djy5YtSE1NxQ8//GC2hVJWVhYGDx4MX19fBAQEYOfOnfXaR9QQVVVV+sV0iajxWMtlWDyhF9Y82xuONnIkZxZg7Kp4HDqfLXY0aobq3bJ05swZ+Pn5QSqtX311/vx59OjRw2zuArpz5w5ycnKgVCqRnZ2N4OBgpKenw97evs59D8OWJaqNTqfDnTt3AADt27ev9+8OERkmK78Us7cmIfXm/eWFXni0MxaMfgQKOX/nqHYmaVkKCgrC77//Xu8Q/fv3R2ZmZr2PN7X27dvrl2pxd3eHq6sr8vPzH7qPyFhSqRQdOnRAhw4dWCgRmZBnazvsfGUAXhzYGQCw7lgGnlpzHFn5pSIno+ai3s0+giDg3XffhZ1d/fqDDZ2x+OjRo1i+fDkSExNx584d7N27FxMnTqx2jEqlwvLly5GdnY3AwECsXr0affv2Neh1ACAxMRFarRaenp4G7SMiIvOkkEvx7jhf9OvSBvN3piL1ZiHGrIrD8skBGOXXXux4ZOHqXSw99thjuHTpUr1P3L9/f9ja1n/+C7VajcDAQLzwwgt44oknauzfvn07oqKisGbNGoSGhmLFihUYOXIkLl26hHbt2gG4Pynmg8aHHDp0CB4eHgCA/Px8PP/881i7dm2N4+ra9weNRlNt0eCioqJ6f4/UslRUVOin2Zg3bx4UCoXIiYiav+G+bjg4LwxztiYhKbMAr2xOQkT/TnhnbE9Yy2VixyMLZdTdcKYmkUhqtCyFhoYiJCQE0dHRAO6PB/H09MScOXOwYMGCep1Xo9Fg+PDhmDVrFp577rl67/uzxYsXY8mSJTW2c8wS/RWXOyEST6VWh08PXcJXv14DAPh1cEL01N7wduXvId1n8rvhmlpFRQUSExMRHh6u3yaVShEeHo4TJ07U6xyCIGDGjBkYOnRojWKorn1/tXDhQhQWFuofWVlZhn9D1CLI5XJEREQgIiLCbG50IGoprGRSLBzdE+tnhqC1vQLnbhVh3Op4/Jh6W+xoZIEsoljKy8uDVquFm5tbte1ubm7Izq7fbaLHjh3D9u3bsW/fPiiVSiiVSpw9e/ah+/7K2toaTk5O2LRpE/r164dhw4Y17JujZsva2hobNmzAhg0bYG1tLXYcohZpSI92ODg3DH07t0aJpgpztiXjnb1nUV6pFTsaWZAW83F34MCB0OkePMNrXftqExkZicjISH0zHhERmSd3ZxtsfSkUK2MuI/rIFWz9LRNJN+5BNb03urZ1EDseWQCLaFlydXWFTCZDTk5Ote05OTlwd3cXJZNKpYKvry9CQkJEeX0iIqo/uUyKN0f0wHcv9IWrgwIXs4sxfnU89ibfFDsaWQCDiqVr165BjPHgCoUCwcHBiImJ0W/T6XSIiYlB//79mzwPcL9lKS0tDadOnRLl9cn8/bHciYuLC9RqtdhxiAhAmE9bHJwbhv5d2qC0Qos3tqfi7V2pKKtgtxzVzqBiycfHB3fv3tV/PWXKlBqtPcYqKSlBSkoKUlJSAAAZGRlISUnRT2wZFRWFtWvXYuPGjbhw4QJeffVVqNVqzJw5s1Fen8gU/rgRgIjMRzsnG2x+KRRvhHeHVALsOH0TE6LjkZ5TLHY0MlMGTR0glUqRnZ2tn9fI0dERqamp6NKlS4ODxMbGYsiQITW2R0REYMOGDQCA6Oho/aSUSqUSq1atQmhoaINf2xgqlQoqlQparRbp6emcOoBq0Ol0uHr1KgCga9eunMWbyAyduPo75n6fjLvFGthYSfGvx/3wVHBH/SLx1HwZMnWA2RRLloprwxERWba8Eg3e2J6CuMt5AIBJQR3wwUQ/2Fu3mHugWiSTzbMkkUhqVNusvomIyJK5Olhj48y+eGtkD0glwN7kWxgfHY8Ld7hCA91ncMvS6NGj9XPG/Pjjjxg6dGiNmYn37NnTuCnNELvh6GEqKyvx9ddfAwD+9re/wcrKSuRERPQwCRn5mLstGdlF5VDIpVg8vhem9vVkw0AzZLJuuPoOpl6/fn19T2nx2A1HteFyJ0SWKV9dgTd3pODIpfs3NI0P9MBHk/zgaMMPPM2JyYolqonFEtWmvLxcv3zOpk2bYGNjI3IiIqovnU7A2rhrWP6fS6jSCfBuY4foab3h14GTEDcXLJaaALvhiIiav8Qb9zB3WzJuFZRBIZPiH2N74vn+ndgt1wywWGpCbFkiImreCkorMH/nGRy+cH9ewdF+7vj4yQA427JbzpKZ7G44IiKilsbFToG1zwfj3XG+sJJJ8PO5bIxbHYfUrAKxo1ETYbFEZCKlpaXo0KEDOnTogNLSUrHjEFEDSCQSvDiwM3a9MgCerW2RlV+GyWuO49v4DFGWAaOmxWLJSFxIlx5GEATcvn0bt2/f5pspUTMR6OmCA3PCMNrPHZVaAe8fSMOs7xJRUFohdjQyIY5ZaiCOWaLaaLVanD17FgDg7+8PmUwmciIiaiyCIGDTyRv44MAFVGh16OBii1VTgxDcqZXY0aieOMC7CbFYIiJquc7dKsTsrUm4/nsp5FIJ3hrZA7PCukAq5d1y5o4DvImIiJqAXwdn/DhnIMYHeqBKJ2Dpzxfx4sZTyFezW645YbFEZCKVlZXYsGEDNmzYgMrKSrHjEJGJONpYYdUzSnw0yR/WcimOXLqLMSvjkJCRL3Y0aiTshjMSJ6Wkh+FyJ0Qtz4U7RYjcmoRrd9WQSSWIGt4drw7qym45M8QxS02IY5aoNuXl5XjyyScBALt37+ZyJ0QthFpThXf3ncOe5FsAgDAfV3z2tBJtHa1FTkZ/xmKpCbFYIiKivxIEATsTb+K9H86hvFKHto7WWPmMEgO6uoodjf4/DvAmIiISkUQiwdN9PPHj7IHwaeeAu8UaPPvNb1hxOB1aHdsoLA2LJSIiIhPxcXPE/tkD8XSfjtAJwIrDl/HsN78ht6hc7GhkABZLRCZSWloKHx8f+Pj4cLkTohbMViHDssmB+HxKIOwUMpy49jvGrIpD3OW7YkejemKxRGQigiDgypUruHLlCpc7ISJMCuqI/bMH4hF3R+SVVOD5dQlY/p+LqNLqxI5GD8FiichEbGxsEB8fj/j4eN4JR0QAgG7tHLAv8lFMC/WCIACqI1cxbe1vuFNYJnY0qgPvhjMS51kiIqKG+DH1NhbuOYsSTRVa2Vnhs6eVGPJIO7FjtRicOqAJceoAIiIy1vU8NSK3JuH87SIAwMuPdcH8kT1gJWPHj6lx6gAiM1BVVYWdO3di586dqKqqEjsOEZkhb1d77H51ACL6dwIAfHX0GqZ8dQK3CtgtZ07YstRAbFmi2nC5EyIyxM9n7+Dt3WdQXF4FZ1srfPpUIIb7uokdq9liyxKRGZBKpRg0aBAGDRoEqZS/akRUt9H+7XFwbhgCOzqjsKwSs747jfcPpKGiinfLiY0tSw3EliUiImpMFVU6fPLLRXwbnwEACPR0QfTUIHi2thM5WfPCliUiIiILpZBL8e44X6x9vg+cba2QmlWAMavi8Mu5O2JHa7FaTLGUlZWFwYMHw9fXFwEBAdi5c6d+X0FBAfr06QOlUgk/Pz+sXbtWxKRERETAcF83HJwXht5eLigur8Irm5Ow6Idz0FRpxY7W4rSYbrg7d+4gJycHSqUS2dnZCA4ORnp6Ouzt7aHVaqHRaGBnZwe1Wg0/Pz+cPn0abdq0eeh52Q1HtSkrK0P//v0BACdOnICtra3IiYjIElVqdfj00CV89es1AIBfBydET+0Nb1feNNIQ7IZ7gPbt20OpVAIA3N3d4erqivz8fACATCaDnd39vmCNRgNBELg8BTWYTqdDamoqUlNTodNxgCYRGcdKJsXC0T2xfkYIWtlZ4dytIoxbHY8DZ26LHa3FMJti6ejRoxg/fjw8PDwgkUiwb9++GseoVCp4e3vDxsYGoaGhSEhIMOq1EhMTodVq4enpqd9WUFCAwMBAdOzYEW+99RZcXV2N/VaIANxf7uTQoUM4dOgQlzshogYb8kg7HJwXhhDvVijRVGH21mS8s/csyivZLWdqZlMsqdVqBAYGQqVSPXD/9u3bERUVhUWLFiEpKQmBgYEYOXIkcnNz9cf8Mebor4/bt/+v+s7Pz8fzzz+Pr7/+utr5XVxckJqaioyMDGzduhU5OTkPzKHRaFBUVFTtQfQgMpkMw4cPx/DhwyGTycSOQ0TNQHtnW2yb1Q+RQ7pCIgG2/paJiapjuHq3ROxozZpZjlmSSCTYu3cvJk6cqN8WGhqKkJAQREdHA7jfxeHp6Yk5c+ZgwYIF9TqvRqPB8OHDMWvWLDz33HO1Hvfaa69h6NChmDx5co19ixcvxpIlS2ps55glIiJqSkfT7+KN7Sn4XV0BO4UMH07yw6SgjmLHshjNbsxSRUUFEhMTER4ert8mlUoRHh6OEydO1OscgiBgxowZGDp0aI1CKScnB8XFxQDuFz1Hjx5Fjx49HniehQsXorCwUP/Iysoy8rui5q6qqgo//fQTfvrpJy53QkSN7rHubfHzvDD069IapRVavLE9FW/vSkVZBbvlGptFFEt5eXnQarVwc6s+7bubmxuys7PrdY5jx45h+/bt2LdvH5RKJZRKJc6ePQsAuHHjBsLCwhAYGIiwsDDMmTMH/v7+DzyPtbU1nJycsGnTJvTr1w/Dhg1r2DdHzZZGo8G4ceMwbtw4aDQaseMQUTPUzskGW17qh3nDfCCRADtO38TjqnhczikWO1qzIhc7QFMZOHBgrXck9e3bFykpKQadLzIyEpGRkfpmPKK/kkql6NOnj/7fRESmIJNK8Mbw7gjt3BrztqcgPacE46Pj8f7jfniqj+fDT0APZRHv4K6urpDJZDUGXefk5MDd3V2UTCqVCr6+vggJCRHl9cn82dra4tSpUzh16hTnWCIikxvQzRUH54ZhYDdXlFfq8NauM4jakQK1hsMAGsoiiiWFQoHg4GDExMTot+l0OsTExOgn/WtqkZGRSEtLw6lTp0R5fSIior9q62iN717oi/kjukMqAfYk3cKE6HhczOad2w1hNsVSSUkJUlJS9N1hGRkZSElJQWZmJgAgKioKa9euxcaNG3HhwgW8+uqrUKvVmDlzpoipiYiIzItUKsHsoT7YNqsf3JyscfWuGo9HH8O2hExOuGwks5k6IDY2FkOGDKmxPSIiAhs2bAAAREdHY/ny5cjOzoZSqcSqVasQGhraxEnvU6lUUKlU0Gq1SE9P59QBVENZWZn+Ds7Dhw+zK46ImtzvJRq8uTMVsZfuAgDGB3rgo0l+cLSxEjmZ+AyZOsBsiiVLxbXhqDZqtRoODg4A7rec2ttzHSciano6nYC1cdew7D+XoNUJ8G5jh+hpveHXoWXfnNTs5lkiskTW1tbYu3cv9u7dC2tra7HjEFELJZVK8PKgrtjxcn94ONvg+u+leOKL49h04jq75eqJLUtGYjccERFZmoLSCszfeQaHL9y/u3yMvzs+fjIATi2wW47dcE2I3XBERGRJBEHAt/EZ+OSXi6jUCvBsbYvoqb0R6OkidrQmxW44IjOg1WoRGxuL2NhYaLVcfoCIzINEIsFLYV2w85UB6NjKFln5ZZi85ji+jc9gt1wt2LJkJHbD0cNwgDcRmbvCskr8fdcZ/HL+/tJhw33dsHxyAFzsFCInMz12wzUhdsNRbUpLS/UzvJ86dQp2dnYiJyIiqkkQBHx34gY+/OkCKrQ6dHCxxaqpQQju1ErsaCbFYqkJsVgiIqLm4NytQkRuTcKN30shl0rw1sgemBXWBVKpROxoJsExS0RERGQQvw7OODBnIMYFtEeVTsDSny/ipe9OI19dIXY00bFYMhIX0iUioubG0cYKq6cG4cNJflDIpfjfxVyMWRmHU9fzxY4mKnbDNRC74ag2ZWVlmDBhAgBg//79XO6EiCxK2u0izN6ahGt5asikEkQN745XB3VtNt1yHLPUhFgsUW14NxwRWTq1pgr/3HcOe5NvAQDCfFzx+RQlXB0sf1UCjlkiMgPW1tbYvHkzNm/ezOVOiMgi2VvL8dnTgVj2ZABsrKSIu5yH0SvjcPxqntjRmhRblhqILUtERNQSpOcUI3JLEi7nlkAqAeYN647ZQ7tBZqHdcmxZagIc4E1ERC1JdzdH/DD7UTwV3BE6Afj8cDqe+/Y35BaXix3N5Niy1EBsWaLaaLVaJCUlAQB69+4NmUwmciIiosaxJ+km/rnvHEortHB1UGDFlCAM9HEVO5ZBOMC7CbFYotpwgDcRNWdXckswe2sSLmYXQyIBZg/phnnDfCCXWUanFbvhiMyARCJBp06d0KlTJ0gkltmnT0RUm27tHLAv8lFM7esFQQBW/+8Kpn3zG7ILm1+3HFuWGogtS0RE1NLtT72NhbvPQF2hRWt7BT57OhCDe7QTO1ad2LJERERETWZCoAcOzA1DLw8n5KsrMGP9KXz880VUanViR2sULJaIiIiowTq72mP3qwPwfP9OAIA1v17FM1+fxO2CMpGTNRyLJSITKS8vx8SJEzFx4kSUlze/Pnwior+ysZLhX4/74cvpveFoLUfijXsYsyoOh9NyxI7WIByzZCSVSgWVSgWtVov09HSOWaIaeDccEbVkmb+XYva2JJy5WQgAeGlgZ7w96hEo5ObRTsOpA5oQB3hTbSorK7FhwwYAwIwZM2BlZSVuICKiJlZRpcPHP1/EumMZAIBATxdETw2CZ2s7kZOxWGpSLJaIiIjqduh8NubvTEVReRUcbeRYPjkAo/zai5qJd8MRERGR2RjRyx0H54UhyMsFxeVVeGVzEhb9cA6aKq3Y0eqFxRKRieh0Opw/fx7nz5+HTtc8bp8lIjJWx1Z22PFyf7z8WBcAwMYTN/Dkl8dxPU8tcrKHYzdcA7EbjmrDAd5ERA/2v4s5eHNHKu6VVsLBWo6Pn/THuACPJs3AbjgiM+Hq6gpXV8taXJKIyNSGPuKGg/PCEOLdCiWaKszemox39p5FeaV5dsu1mGIpKysLgwcPhq+vLwICArBz584ax5SWlqJTp06YP3++CAmpubG3t8fdu3dx9+5dtioREf1Fe2dbbJvVD68N7goA2PpbJiaqjuHq3RKRk9XUYooluVyOFStWIC0tDYcOHcLrr78Otbp6P+mHH36Ifv36iZSQiIioZZHLpHh71CPY+EJftLFX4GJ2Mcavjse+5FtiR6umxRRL7du3h1KpBAC4u7vD1dUV+fn5+v2XL1/GxYsXMXr0aJESEhERtUyDurfFwXlh6NelNUortHh9ewr+vusMyirMo1vObIqlo0ePYvz48fDw8IBEIsG+fftqHKNSqeDt7Q0bGxuEhoYiISHBqNdKTEyEVquFp6enftv8+fOxdOlSY+MT1VBeXo7p06dj+vTpXO6EiOgh3JxssOWlfpg7zAcSCbD9dBYeV8Xjck6x2NHMp1hSq9UIDAyESqV64P7t27cjKioKixYtQlJSEgIDAzFy5Ejk5ubqj1EqlfDz86vxuH37tv6Y/Px8PP/88/j666/123744Qd0794d3bt3f2hOjUaDoqKiag+iB9Fqtdi6dSu2bt0KrdY8Ph0REZkzmVSCqOHdseXFULR1tEZ6TgkmRB/DztNZouYyy6kDJBIJ9u7di4kTJ+q3hYaGIiQkBNHR0QDuz2Hj6emJOXPmYMGCBfU6r0ajwfDhwzFr1iw899xz+u0LFy7E5s2bIZPJUFJSgsrKSrz55pt47733apxj8eLFWLJkSY3tnDqA/qqyslJf/EdGRnK5EyIiA9wt1uCN7SmIv5KHHm6O+HHOwEZdV87ilzv5a7FUUVEBOzs77Nq1q1oBFRERgYKCAvzwww8PPacgCJg2bRp69OiBxYsX13rchg0bcO7cOXz66acP3K/RaKDRaPRfFxUVwdPTk8USERFRI9PqBKz59SpG9nJDt3aOjXruZjfPUl5eHrRaLdzc3Kptd3NzQ3Z2dr3OcezYMWzfvh379u2DUqmEUqnE2bNnDc5ibW0NJycnbNq0Cf369cOwYcMMPgcRERE9nEwqQeSQbo1eKBlKLuqrN6GBAwfWa8mJGTNm1Ot8kZGRiIyM1FemRH+l0+mQmZkJAPDy8oJUahGfTYiI6C8solhydXWFTCZDTk5Ote05OTlwd3cXJZNKpYJKpeLAXapVWVkZOnfuDIDLnRARWTKL+KirUCgQHByMmJgY/TadToeYmBj0799flEyRkZFIS0vDqVOnRHl9sgx2dnaws7MTOwYRETWA2bQslZSU4MqVK/qvMzIykJKSgtatW8PLywtRUVGIiIhAnz590LdvX6xYsQJqtRozZ84UMTVR7ezt7WvMEk9ERJbHbIql06dPY8iQIfqvo6KiANy/423Dhg2YMmUK7t69i/feew/Z2dlQKpX45Zdfagz6birshiMiImoZzHLqAEtiyK2HREREZB6a3dQBRJZIo9Fg1qxZmDVrVrW5uYiIyLKwZclIf+6GS09PZ8sS1aBWq+Hg4ACAd8MREZkbi5/B25KwG45qU1FRgeXLlwMA3nrrLSgUCpETERHRH1gsNSEWS0RERJaHY5aagEqlgq+vL0JCQsSOQkRERCbElqUGYssS1UYQBOTl5QG4Pwu9RCIROREREf3BkL/fZjPPElFzU1painbt2gHgAG8iIkvGYqmB/miYKyoqEjkJmZs/z95dVFTECUyJiMzIH3+369PBxmLJSH9MHVBRUQEA8PT0FDkRmTMPDw+xIxAR0QMUFxfD2dm5zmM4ZqmBdDodbt++DUdHx2pjUkJCQmpdZLe2fX/dXlRUBE9PT2RlZYk+Hqqu76cpz2fI8x52bEP2P2jfg7aZyzVsjtevPsfU93ettu3mcv2AlnkN+T5qmvPxGt4nCAKKi4vh4eEBqbTu+93YstRAUqkUHTt2rLFdJpPVelFr21fbdicnJ9F/yev6fpryfIY872HHNmT/g/bVdbzY17A5Xr/6HGPo7xp/Bxv3eQ29hnwfNc35eA3/z8NalP7AqQNMJDIy0uB9dT1HbI2dzdjzGfK8hx3bkP0P2sfr17jPq8+xxl5D/g5axjXk+6hpzsdraDh2w5kxTktg+XgNLRuvn+XjNbR85nAN2bJkxqytrbFo0SJYW1uLHYWMxGto2Xj9LB+voeUzh2vIliUiIiKiOrBliYiIiKgOLJaIiIiI6sBiiYiIiKgOLJaIiIiI6sBiiYiIiKgOLJYsWGlpKTp16oT58+eLHYWM4O3tjYCAACiVSgwZMkTsOGSEjIwMDBkyBL6+vvD396+2eDKZt0uXLkGpVOoftra22Ldvn9ixyECff/45evXqBV9fX8ydO7dei+Iag8udWLAPP/wQ/fr1EzsGNcDx48fh4OAgdgwy0owZM/DBBx8gLCwM+fn5nMvHgvTo0QMpKSkAgJKSEnh7e2P48OHihiKD3L17F9HR0Th//jysrKzw2GOP4eTJk+jfv3+jvxZblizU5cuXcfHiRYwePVrsKEQt0h9v0GFhYQCA1q1bQy7n509LtH//fgwbNgz29vZiRyEDVVVVoby8HJWVlaisrES7du1M8joslkRw9OhRjB8/Hh4eHpBIJA9s+lWpVPD29oaNjQ1CQ0ORkJBQbf/8+fOxdOnSJkpMf9UY11AikWDQoEEICQnBli1bmig5/aGh1/Dy5ctwcHDA+PHj0bt3b3z00UdNmJ4a43fwDzt27MCUKVNMnJj+qqHXsG3btpg/fz68vLzg4eGB8PBwdO3a1SRZWSyJQK1WIzAwECqV6oH7t2/fjqioKCxatAhJSUkIDAzEyJEjkZubCwD44Ycf0L17d3Tv3r0pY9OfNPQaAkB8fDwSExOxf/9+fPTRRzhz5kxTxSc0/BpWVVUhLi4OX3zxBU6cOIH//ve/+O9//9uU30KL1hi/g8D9dceOHz+OMWPGNEVs+pOGXsN79+7hwIEDuH79Om7duoXjx4/j6NGjpgkrkKgACHv37q22rW/fvkJkZKT+a61WK3h4eAhLly4VBEEQFixYIHTs2FHo1KmT0KZNG8HJyUlYsmRJU8amPzHmGv7V/PnzhfXr15swJdXFmGt4/PhxYcSIEfr9y5YtE5YtW9Ykeam6hvwOfvfdd8L06dObIibVwZhruGPHDuG1117T71+2bJnwySefmCQfW5bMTEVFBRITExEeHq7fJpVKER4ejhMnTgAAli5diqysLFy/fh2ffvopZs2ahffee0+syPQX9bmGarUaxcXFAO4PLv3f//6HXr16iZKXaqrPNQwJCUFubi7u3bsHnU6Ho0ePomfPnmJFpj+pz/X7A7vgzFN9rqGnpyeOHz+O8vJyaLVaxMbGokePHibJw9GIZiYvLw9arRZubm7Vtru5ueHixYsipSJD1Oca5uTkYNKkSQAArVaLWbNmISQkpMmz0oPV5xrK5XJ89NFHeOyxxyAIAkaMGIFx48aJEZf+or7vo4WFhUhISMDu3bubOiI9RH2uYb9+/TBmzBgEBQVBKpVi2LBhmDBhgknysFiycDNmzBA7AhmhS5cuSE1NFTsGNdDo0aN5R6oFc3Z2Rk5OjtgxqAE+/PBDfPjhhyZ/HXbDmRlXV1fIZLIav8A5OTlwd3cXKRUZgtfQ8vEaWjZeP8tnbteQxZKZUSgUCA4ORkxMjH6bTqdDTEyMSSbaosbHa2j5eA0tG6+f5TO3a8huOBGUlJTgypUr+q8zMjKQkpKC1q1bw8vLC1FRUYiIiECfPn3Qt29frFixAmq1GjNnzhQxNf0Zr6Hl4zW0bLx+ls+irqFJ7rGjOh05ckQAUOMRERGhP2b16tWCl5eXoFAohL59+wonT54ULzDVwGto+XgNLRuvn+WzpGsoEQQTrTpHRERE1AxwzBIRERFRHVgsEREREdWBxRIRERFRHVgsEREREdWBxRIRERFRHVgsEREREdWBxRIRERFRHVgsEREREdWBxRIRERFRHVgsEVGLMmPGDEgkEkgkEuzbt0+0HLGxsfocEydOFC0HET0ciyUislh/Lnz+/Bg1alSdzxs1ahTu3LmD0aNHV9t+5MgRjBs3Dm3btoWNjQ26du2KKVOm4OjRo/XO5O/vj1deeeWB+zZt2gRra2vk5eVhwIABuHPnDp5++ul6n5uIxMFiiYgs2h+Fz58f27Ztq/M51tbWcHd3h7W1tX7bF198gWHDhqFNmzbYvn07Ll26hL1792LAgAF444036p3nxRdfxPfff4+ysrIa+9avX48JEybA1dUVCoUC7u7usLW1rf83S0SiYLFERBbtj8Lnz49WrVoZdI7MzEy8/vrreP3117Fx40YMHToUnTp1QkBAAObNm4fTp09XOz4+Ph5hYWGwtbWFp6cn5s6dC7VaDQB49tlnUVZWht27d1d7TkZGBmJjY/Hiiy827BsmoibHYomIWrzdu3ejsrISb7/99gP3SyQS/b+vXr2KUaNG4cknn8SZM2ewfft2xMfHY/bs2QAAV1dXPP7441i3bl21c2zYsAEdO3bEiBEjTPeNEJFJsFgiIot24MABODg4VHt89NFHBp0jPT0dTk5OcHd312/bvXt3tXOePXsWALB06VJMnz4dr7/+Onx8fDBgwACsWrUK3333HcrLywHc74qLjY1FRkYGAEAQBGzcuBERERGQSvm2S2Rp5GIHICJqiCFDhuDLL7+stq1169YGn+fPrUcAMHLkSKSkpODWrVsYPHgwtFotACA1NRVnzpzBli1b9McKggCdToeMjAz07NkTw4cPR8eOHbF+/Xr861//QkxMDDIzMzFz5kwjvkMiEhuLJSKyaPb29ujWrVuDzuHj44PCwkJkZ2frW5ccHBzQrVs3yOXV3yZLSkrw8ssvY+7cuTXO4+XlBQCQSqWYMWMGNm7ciMWLF2P9+vUYMmQIunTp0qCcRCQOtgcTUYs3efJkWFlZ4ZNPPnnosb1790ZaWhq6detW46FQKPTHzZw5E1lZWdizZw/27t3Lgd1EFowtS0Rk0TQaDbKzs6ttk8vlcHV1rfc5vLy88O9//xvz5s1Dfn4+ZsyYgc6dOyM/Px+bN28GAMhkMgDA3//+d/Tr1w+zZ8/GSy+9BHt7e6SlpeG///0voqOj9efs3Lkzhg4dir/97W+wtrbGE0880QjfLRGJgS1LRGTRfvnlF7Rv377aY+DAgQafZ86cOTh06BDu3r2LyZMnw8fHB2PGjEFGRgZ++eUX+Pv7AwACAgLw66+/Ij09HWFhYQgKCsJ7770HDw+PGud88cUXce/ePUybNg02NjYN/l6JSBwSQRAEsUMQETWVGTNmoKCgQNSlTv7M3PIQUU1sWSKiFueP6QYOHDggWoa4uDg4ODhUu6uOiMwTW5aIqEXJzc1FUVERAKB9+/awt7cXJUdZWRlu3boF4P6dd3+e44mIzAuLJSIiIqI6sBuOiIiIqA4sloiIiIjqwGKJiIiIqA4sloiIiIjqwGKJiIiIqA4sloiIiIjqwGKJiIiIqA4sloiIiIjq8P8AS7PM4bpuMFAAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "energies = np.geomspace(min_energy, max_energy)\n", "\n", "fig, ax = plt.subplots()\n", "ax.plot(energies, [power_law.spectrum(e) for e in energies])\n", "ax.axhline(flux_norm, color=\"k\", linestyle=\":\")\n", "ax.axvline(norm_energy, color=\"k\", linestyle=\":\")\n", "ax.set_xscale(\"log\")\n", "ax.set_yscale(\"log\")\n", "ax.set_xlabel(\"E [GeV]\")\n", "ax.set_ylabel(\"F $[GeV^-1 cm^-2 s^-1 (sr^-1)]$\")" ] }, { "cell_type": "markdown", "id": "308c4bb5", "metadata": {}, "source": [ "We can also use the `PowerLawFlux` class to perform some simple calculations, such as integration of the flux." ] }, { "cell_type": "code", "execution_count": 4, "id": "bd2ee5dc", "metadata": { "execution": { "iopub.execute_input": "2024-03-22T09:57:02.474743Z", "iopub.status.busy": "2024-03-22T09:57:02.474420Z", "iopub.status.idle": "2024-03-22T09:57:02.479677Z", "shell.execute_reply": "2024-03-22T09:57:02.479142Z" } }, "outputs": [ { "data": { "text/plain": [ "array([9.999e-13])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "total_flux = power_law.integrated_spectrum(min_energy, max_energy) # cm^-2 s^-1 (sr^-1)\n", "total_flux" ] }, { "cell_type": "code", "execution_count": 5, "id": "967dcd06", "metadata": { "execution": { "iopub.execute_input": "2024-03-22T09:57:02.482005Z", "iopub.status.busy": "2024-03-22T09:57:02.481630Z", "iopub.status.idle": "2024-03-22T09:57:02.485870Z", "shell.execute_reply": "2024-03-22T09:57:02.485335Z" } }, "outputs": [ { "data": { "text/plain": [ "9.210340371976183e-08" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "total_energy_flux = power_law.total_flux_density() # GeV cm^-2 s^-1 (sr^-1)\n", "total_energy_flux" ] }, { "cell_type": "markdown", "id": "b2cd56a4", "metadata": {}, "source": [ "Sampling from the power law shape is also possible:" ] }, { "cell_type": "code", "execution_count": 6, "id": "e064de49", "metadata": { "execution": { "iopub.execute_input": "2024-03-22T09:57:02.488270Z", "iopub.status.busy": "2024-03-22T09:57:02.487910Z", "iopub.status.idle": "2024-03-22T09:57:03.294714Z", "shell.execute_reply": "2024-03-22T09:57:03.294033Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAG1CAYAAAABTQXdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABWlUlEQVR4nO3deVxU9f7H8dcMq4qghmIoRmlWuIAaopULiuG+pGa7W6aGWZdrpd2bZou0WplTmKmYbVamlXpNIxXX6xa2uGSGSRqoqSCoLDPz+2N+cSM3wIHDDO/n43Ee955lznmPJ+Tj+S7HZLfb7YiIiIi4OLPRAUREREScQUWNiIiIuAUVNSIiIuIWVNSIiIiIW1BRIyIiIm5BRY2IiIi4BRU1IiIi4hZU1IiIiIhb8DQ6QEWx2WwcPnyYmjVrYjKZjI4jIiIiJWC32zl16hTBwcGYzRd/FlNliprDhw8TEhJidAwREREpg/T0dBo2bHjRY6pMUVOzZk3A8Yfi7+9vcBoREREpiezsbEJCQop+j19MlSlq/mxy8vf3V1EjIiLiYkrSdUQdhUVERMQtuH1RY7FYCAsLIzIy0ugoIiIiUo5MdrvdbnSIipCdnU1AQABZWVlqfhIRqSKsVisFBQVGx5BL8Pb2vuDIptL8/q4yfWpERKTqsNvtZGRkcPLkSaOjSAmYzWauvvpqvL29L+s8KmpERMTt/FnQ1KtXj+rVq2t+skrsz3nkfv/9dxo1anRZ90pFjYiIuBWr1VpU0FxxxRVGx5ESqFu3LocPH6awsBAvL68yn8ftOwqLiEjV8mcfmurVqxucRErqz2Ynq9V6WedRUSMiIm5JTU6uw1n3yuWKmtOnT3PVVVcxYcIEo6OIiIhIJeJyRc1zzz1Hu3btjI4hIiIilYxLdRTet28fe/bsoU+fPvzwww9GxxERERcTOnFZhV3rwPO9KuxaFcFkMrF48WL69+9vdJQLctqTmpSUFPr06UNwcDAmk4klS5acc4zFYiE0NBRfX1+ioqLYsmVLqa4xYcIEEhISnJTYiaya2ElERC7f0aNHGTt2LI0aNcLHx4f69esTGxvLhg0bjI7mEpz2pCY3N5fw8HBGjBjBbbfdds7+hQsXEh8fT2JiIlFRUbz22mvExsayd+9e6tWrB0BERASFhYXnfHblypVs3bqVpk2b0rRpUzZu3HjJPHl5eeTl5RWtZ2dnX8a3u4jCfJjXA67pDJ0ngYdLPfwSEZFKZODAgeTn5zN//nyuueYaMjMzSU5O5o8//jA6mktw2pOaHj168OyzzzJgwIDz7p8+fTqjRo1i+PDhhIWFkZiYSPXq1Zk7d27RMampqfzwww/nLMHBwWzevJmPPvqI0NBQJkyYwOzZs3n66acvmCchIYGAgICiJSQkxFlftbi9y+HQNlj3MszvA1mHyuc6IiLi1k6ePMm6det44YUXiI6O5qqrrqJt27ZMmjSJvn37Ao7fpS1atKBGjRqEhITw4IMPkpOTU3SOpKQkatWqxdKlS7nuuuuoXr06gwYN4vTp08yfP5/Q0FBq167N+PHjiw2fDg0N5ZlnnuHOO++kRo0aNGjQAIvFctG86enp3H777dSqVYs6derQr18/Dhw4ULR/zZo1tG3blho1alCrVi1uvvlmfv31V+f+of1NhXQUzs/PZ/v27cTExPzvwmYzMTExbNq0qUTnSEhIID09nQMHDvDyyy8zatQoJk+efMHjJ02aRFZWVtGSnp5+2d/jvJr1h4FzwLsmHNwIibfATyvL51oiIuK2/Pz88PPzY8mSJcVaGv7KbDYzY8YMfvzxR+bPn88333zDY489VuyY06dPM2PGDD766CNWrFjBmjVrGDBgAMuXL2f58uUsWLCAWbNm8emnnxb73EsvvUR4eDjffvstEydO5OGHH2bVqlXnzVFQUEBsbCw1a9Zk3bp1bNiwAT8/P7p3705+fj6FhYX079+fTp068d1337Fp0yYeeOCBch9mXyFtJceOHcNqtRIUFFRse1BQEHv27CmXa/r4+ODj44PFYsFisVz2hD4X1WIQBLeCT4ZBxnfwwWC4aTx0nQweZZ8ZUUREqg5PT0+SkpIYNWoUiYmJtG7dmk6dOnHHHXfQsmVLAB555JGi40NDQ3n22WcZM2YMb775ZtH2goIC3nrrLRo3bgzAoEGDWLBgAZmZmfj5+REWFkZ0dDSrV69myJAhRZ+7+eabmThxIgBNmzZlw4YNvPrqq3Tr1u2crAsXLsRms/HOO+8UFSrz5s2jVq1arFmzhhtvvJGsrCx69+5dlOOGG25w7h/YebjckG6AYcOG8fLLL5fo2Li4OHbt2sXWrVvLN9QVjWHkKmj7gGN94wyY1xNOltMTIhERcTsDBw7k8OHDfPHFF3Tv3p01a9bQunVrkpKSAPj666/p2rUrDRo0oGbNmtx777388ccfnD59uugc1atXLyokwPEAITQ0FD8/v2Lbjhw5Uuza7du3P2d99+7d5825c+dOfv75Z2rWrFn0hKlOnTqcPXuW/fv3U6dOHYYNG0ZsbCx9+vTh9ddf5/fff7/cP55LqpCiJjAwEA8PDzIzM4ttz8zMpH79+uV6bYvFQlhYGJGRkeV6ndCJywh9MpnQlM6MyX+EbHt1+G0LJ1+N4v4nnq7QYYQiIuK6fH196datG08++SQbN25k2LBhTJkyhQMHDtC7d29atmzJokWL2L59e1G/l/z8/KLP//3dSSaT6bzbbDZbmTPm5OTQpk0bUlNTiy0//fQTd911F+B4crNp0yZuuukmFi5cSNOmTdm8eXOZr1kSFVLUeHt706ZNG5KTk4u22Ww2kpOTz6kMna3CntT8xQpbW3rlP0eq7RpqmXJ5x/sVnvRc4BgpJSIiUgphYWHk5uayfft2bDYbr7zyCu3ataNp06YcPnzYadf5e8GxefPmCzYZtW7dmn379lGvXj2aNGlSbAkICCg6rlWrVkyaNImNGzfSvHlzPvjgA6flPR+nFTU5OTlFlRpAWloaqampHDx4EID4+Hhmz57N/Pnz2b17N2PHjiU3N5fhw4c7K0Klkm4PYnD+U7xT2AOAkZ7/gbmxcOKAscFERKRS+uOPP+jSpQvvvfce3333HWlpaXzyySe8+OKL9OvXjyZNmlBQUMAbb7zBL7/8woIFC0hMTHTa9Tds2MCLL77ITz/9hMVi4ZNPPuHhhx8+77F33303gYGB9OvXj3Xr1pGWlsaaNWsYP348v/32G2lpaUyaNIlNmzbx66+/snLlSvbt21fu/Wqc1lF427ZtREdHF63Hx8cDMHToUJKSkhgyZAhHjx5l8uTJZGRkEBERwYoVK87pPOxsFdJR+AIK8OTZwnvZbAvjZa9Eah3eAYkdod9MCOtb4XlERKq6yjzLr5+fH1FRUbz66qvs37+fgoICQkJCGDVqFE888QTVqlVj+vTpvPDCC0yaNImOHTuSkJDAfffd55Tr//Of/2Tbtm1MnToVf39/pk+fTmxs7HmPrV69OikpKTz++OPcdtttnDp1igYNGtC1a1f8/f05c+YMe/bsYf78+fzxxx9ceeWVxMXFMXr0aKdkvRCT3W63l+sVKons7GwCAgLIysrC39/f6ee/VJ+ZYI6xscl78Nv/z6Lc9gG49Vnw9HF6FhGRquzs2bOkpaVx9dVX4+vra3QclxAaGsojjzxSbHRVRbrYPSvN72+XHP3kig4TCMOXw83//yhvy9swpxv8sd/YYCIiIm7C7Yuaihr9VCIeXtDtabjrE6hWB37fCbM6wQ+LjE4mIiLi8tT85CSlHbJdnz943dtClNkx+eD7hV15uvBe9j5//tdMiIhIyaj5yfWo+cnFZXAFd+X/izcK+2Ozm7jbM5kl3pPh2D6jo4mIiLgkty9qKlXz099Y8eCVwtu5r2AiR+3+3GA+6GiO+u5jo6OJiIi4HLcvaoyYfK+01tta0DMvgY3WMCjIhc9GwefjIP/0pT8sIiIiQBUoalzFUWpzT8ET0HkSYIJvF8DsLnCkfF74KSIi4m4q5C3dUjI2zISuaEF78xO87mWh3tHdnLF04MnC4Xxq7VR0XGWePEpERMQobv+kpjL3qbmQTbZm9MxLIMXagmqmfF72msUrXm9RnbNGRxMRERe2Zs0aTCYTJ0+eLPFnQkNDee2118otkzO5fVHjCn1qzucYAQwteJwXC27Hajcx0GMdX3j/m+tNB42OJiIi5WTYsGGYTCbGjBlzzr64uDhMJhPDhg2r+GAuwu2LGldmx8yb1v7cmf9vMuy1aWI+zBLvJ2HbPKga0wuJiFQ5ISEhfPTRR5w5c6Zo29mzZ/nggw9o1KiRgckqPxU1LmCL/QZ65iWw2hqOr6kAlj4Ci0bC2Wyjo4mIiJO1bt2akJAQPvvss6Jtn332GY0aNaJVq1ZF2/Ly8hg/fjz16tXD19eXW2655ZxWieXLl9O0aVOqVatGdHQ0Bw4cOOd669evp0OHDlSrVo2QkBDGjx9Pbm5uuX2/8qSixkUcx58RBY+SUHAnmDwcr1Z4u5PjVQsiInJxdjvk51b8Usan6iNGjGDevHlF63PnzmX48OHFjnnsscdYtGgR8+fPZ8eOHTRp0oTY2FiOHz8OQHp6Orfddht9+vQhNTWV+++/n4kTJxY7x/79++nevTsDBw7ku+++Y+HChaxfv55x48aVKbfR3H70k8ViwWKxYLVajY5y2eyYmWXtw6Thw+HTEXD8F3gnBmKnQeT9YDIZHVFEpHIqOA3Tgiv+uk8cBu8apf7YPffcw6RJk/j1118B2LBhAx999BFr1qwBIDc3l7feeoukpCR69OgBwOzZs1m1ahVz5szh0Ucf5a233qJx48a88sorAFx33XV8//33vPDCC0XXSUhI4O677y56O/e1117LjBkz6NSpE2+99ZbLvWbC7Z/UuGpH4YtqFAVj1kHTHmDNh+UT4JOhcOak0clERMQJ6tatS69evUhKSmLevHn06tWLwMDAov379++noKCAm2++uWibl5cXbdu2Zffu3QDs3r2bqKioYudt3759sfWdO3eSlJSEn59f0RIbG4vNZiMtLa0cv2H5cPsnNe7ofy/PvIeRHlfwuOeHeO/6nIM/bGRcwXi+szcGNJ+NiEgRr+qOpyZGXLeMRowYUdQMZLFYnJWomJycHEaPHs348ePP2eeKnZJV1Lg0E3OsPdlma8pMrzdoZD7Kp95P8XzhXcy1djc6nIhI5WEylakZyEjdu3cnPz8fk8lEbGxssX2NGzfG29ubDRs2cNVVVwFQUFDA1q1bi5qSbrjhBr744otin9u8eXOx9datW7Nr1y6aNGlSfl+kArl981NVsNPehF750/iPNRJvk5XJXguY7TUdTh83OpqIiJSRh4cHu3fvZteuXXh4eBTbV6NGDcaOHcujjz7KihUr2LVrF6NGjeL06dOMHDkSgDFjxrBv3z4effRR9u7dywcffEBSUlKx8zz++ONs3LiRcePGkZqayr59+/j8889dtqOwiho3kU0NxhY8wuSCoeTZPenmsR1mdYT0LUZHExGRMvL398ff3/+8+55//nkGDhzIvffeS+vWrfn555/56quvqF27NuBoPlq0aBFLliwhPDycxMREpk2bVuwcLVu2ZO3atfz000906NCBVq1aMXnyZIKDDehU7QQmu71qzOKWnZ1NQEAAWVlZF/wP5HL8r5+L8ZqZ0rB4zSDUnAlmT+g6Gdo/BGbVsCLi/s6ePUtaWhpXX321y43eqaouds9K8/tbv+Xc0I/2q+md/xw0uw1shbBqMnw4BHL/MDqaiIhIuXH7osYVX2jpDDlUh0Fzofdr4OED+1ZC4i3w60ajo4mIiJQLty9q3HKempIymeDG4TDqG7jiWjh1GJJ6Q8rLYLMZnU5ERMSp3L6oEaB+c3hgDbQcAnYrfPMMvD8Qco4anUxERMRpVNRUFT5+MGAW9LOAZzXY/w0k3gxpKUYnExERcQpNvufGzj8iqzbXmqZi8XqdpjmH4N1+0Olx6PgomD3Oc7yIiGuqIoN73YKz7pWKmipon70hffOf5WnPJG73XAtrEtiQ/AWPFMRxlFrnHK/XLYiIK/Hy8gLg9OnTVKtWzeA0UhL5+fkA50wyWFoqaqqos/jwWOFoNtnCeNZrLjd7/Mhy80QeKYhjg62F0fFERMrMw8ODWrVqceTIEQCqV6+OyWQyOJVciM1m4+jRo1SvXh1Pz8srS1yqqAkNDcXf3x+z2Uzt2rVZvXq10ZFc3mJbB77Lv4aZXjO4wZzOAq/nsVj78VrhQKyoOUpEXFP9+vUBigobqdzMZjONGjW67OLTpYoagI0bN+Ln52d0DLey396A/vnPMNlzAXd7JvOQ5xLamvcwPn8cmdQxOp6ISKmZTCauvPJK6tWrR0FBgdFx5BK8vb0xO2HWe5craqR85OHNvwpHstl2A9O85hBl3sNyn0nEFzwIqE+NiLgmDw+Py+6nIa7DaUO6U1JS6NOnD8HBwZhMJpYsWXLOMRaLhdDQUHx9fYmKimLLltK9bNFkMtGpUyciIyN5//33nZRc/upL2030yX+WH2yhXGE6xXzvF2DVFLDqXzoiIlK5Oe1JTW5uLuHh4YwYMYLbbrvtnP0LFy4kPj6exMREoqKieO2114iNjWXv3r3Uq1cPgIiICAoLC8/57MqVKwkODmb9+vU0aNCA33//nZiYGFq0aEHLli3PmycvL4+8vLyi9ezsbCd9U/d3wH4lA/Of4gnP9xnquQo2vAYHN8HAOVArxOh4IiIi51Uub+k2mUwsXryY/v37F22LiooiMjKSmTNnAo7eziEhITz00ENMnDix1Nd49NFHadasGcOGDTvv/qeeeoqpU6ees70qvKXbmbqbt5BYcx7kZUG12tD/Lbiuh9GxRESkiqh0b+nOz89n+/btxMTE/O/CZjMxMTFs2rSpROfIzc3l1KlTAOTk5PDNN9/QrFmzCx4/adIksrKyipb09PTL+xJV1ApbWxi9FoJbw5kT8OEd8NW/oDDf6GgiIiLFVEhH4WPHjmG1WgkKCiq2PSgoiD179pToHJmZmQwYMAAAq9XKqFGjLvrmbR8fH3x8fLBYLFgsFqxWa9m/QBUX+uIuvHiEiZ4fMtLzP7BpJqkb/sO4gof4ze5oOtQEfSIiYjSXGf10zTXXsHPnzlJ/Li4ujri4uKLHV1I2BXjyTOG9bLbdwMteiUSY97Pc+wkeLRjNV7YLF5ciIiIVpUKanwIDA/Hw8CAzM7PY9szMzKIJksqLxWIhLCzsok91pORW2W6kZ14CO2xN8DedZpb3q0zxnA+FeZf+sIiISDmqkKLG29ubNm3akJycXLTNZrORnJxM+/bty/XacXFx7Nq1i61bt5brdaqSQ9Tl9vzJJBb2BmC451cwpxsc/8XgZCIiUpU5rfkpJyeHn3/+uWg9LS2N1NRU6tSpQ6NGjYiPj2fo0KHceOONtG3bltdee43c3FyGDx/urAhSgQrx5PnCu/iv7QZe8XqLOr/v5NTrNzGxYBTLbO3O+xn1uxERkfLktKJm27ZtREdHF63Hx8cDMHToUJKSkhgyZAhHjx5l8uTJZGRkEBERwYoVK87pPOxs6ihcvlbbWtEzL4EZ3jNpa96LxXsG7Qt/5JnCe8nD2+h4IiJShZTLPDWVUWnGuZeFu85TU1IeWPmH56c86PEFZpOd3bZGxBWM5xd7cNExelIjIiKlVenmqRH3Z8WDlwuHMLTgcY7Z/bnBfJAvvf9FP/N6o6OJiEgV4fZFjUY/Vax1tpb0zEtgkzWMGqY8Xvd+k+c938YXjY4SEZHypeYnJ6nqzU9/Z8bGeM/PGO+xGLPJzl5bQ64btwjqXW90NBERcSFqfhLD2TDzWuEg7i54giP2Wlxn/g1mR8O3eru6iIiUD7cvatT8ZKxNtmb0zEtgnbU5FJyGzx+ExWMgL8foaCIi4mbU/OQkan66OBM20nrugdXTwG6DwKYwOAmCLvxSUhERkdL8/naZdz+Ja7NjJnR5GG1N/2KG90zqH/uJs2924qnCoXxkjQZMRcdq6LeIiJSF2zc/SeWyxX4DPfMSWGMNx9dUwPNe7/C6l4UanDE6moiIuDi3L2rUp6byOY4/wwseJaHgTgrtZvp5bGSp9xM0Mx0wOpqIiLgwty9q9ELLysmOmVnWPtyeP5lD9iu42pzJZ95TuMdjFVSNbl4iIuJkbl/USOW2w96UXnnTWGVtjY+pgGe95sEnQ+FsltHRRETExaioEcOdpCajCv7JMwX3kG/3gF2fw6yOcGiH0dFERMSFqKiRSsLEHGtPBudPgVqN4MQBmHMrbH5LzVEiIlIibl/UqKOwa9lpbwKj18ENfcBWACsmwsJ74MwJo6OJiEglp8n3nEST7zmbnfs8VvIvz/fxMRXymz2QcfnjSbU3ATSXjYhIVaF3P4kbMPGuNZaB+U9xwBZEQ9MxPvGeyv0eyzBhMzqciIhUQipqpFL7wX4NvfOfY6m1HV4mK//2ep93vF6B08eNjiYiIpWMihqp9HKozriCh3iiYCR5di+6enwLibfAwc1GRxMRkUpERY24CBMfWLvSP/9p9tuuhOxDMK8nrHsFbGqOEhERFTXiYnbbr6Jv/rPQ4nawWyH5aXh/EOQcNTqaiIgYTEWNuJxcqsFtb0PfmeBZDfYnO5qjDqw3OpqIiBjI0+gA5c1isWCxWLBarUZHEScKnbQcqENT01NYvGZwbc4hrPN683rhQGZa+2PDrGHfIiJVjNs/qdELLd3bT/YQ+uY/w8eFnfAw2Yn3+pQFXgnU5aTR0UREpIK5fVEj7u8MvjxWOJp/5I/ltN2Hmz1+ZLnPRNi/2uhoIiJSgVTUiNtYbOtAn/xn2W0Loa4pGxYMgG+eBWuh0dFERKQCqKgRt7Lf3oD++c/wQWEXwA4pL8G7fSH7sNHRRESknKmoEbeThzdPFN4PA+eAtx/8usExOmrf10ZHExGRcqSiRtxXi0EwOgXqt4DTf8D7A+Hrp8BaYHQyEREpBy5V1KSlpREdHU1YWBgtWrQgNzfX6EhS2V3RGEZ+DZGjHOvrX4WkXpD1m7G5RETE6VyqqBk2bBhPP/00u3btYu3atfj4+BgdSVyBly/0ehkGzwcff0j/r6M5au8Ko5OJiIgTuUxR8+OPP+Ll5UWHDh0AqFOnDp6ebj93oDhTs/6O5qjgVnDmBHw4BL76FxTmG51MREScwGlFTUpKCn369CE4OBiTycSSJUvOOcZisRAaGoqvry9RUVFs2bKlxOfft28ffn5+9OnTh9atWzNt2jRnRZeqpM7VMOIraPegY33TTJjXA078amwuERG5bE4ranJzcwkPD8disZx3/8KFC4mPj2fKlCns2LGD8PBwYmNjOXLkSNExERERNG/e/Jzl8OHDFBYWsm7dOt588002bdrEqlWrWLVq1QXz5OXlkZ2dXWwRAcDTB7onwB0fgG8AHNoGszrA7i+NTiYiIpfBae03PXr0oEePHhfcP336dEaNGsXw4cMBSExMZNmyZcydO5eJEycCkJqaesHPN2jQgBtvvJGQkBAAevbsSWpqKt26dTvv8QkJCUydOrWM30aqhOt7wZj18OkI+G0rLLwH2o6GW59xFD4iIuJSKqRTSn5+Ptu3b2fSpElF28xmMzExMWzatKlE54iMjOTIkSOcOHGCgIAAUlJSGD169AWPnzRpEvHx8UXr2dnZRQWRVA2hE5dd8pgDz/eC4f+B5Kmw8Q3YMsvRkXjwPKhzTQWkFBERZ6mQjsLHjh3DarUSFBRUbHtQUBAZGRklOoenpyfTpk2jY8eOtGzZkmuvvZbevXtf8HgfHx/8/f1ZsGAB7dq1o2vXrpf1HcSNeXjBrc/CXR9DtTrweyokdoQfPjM6mYiIlIJLDR+6VBPX+cTFxREXF0d2djYBAQHllEzcQtPY/zVHpW+GT4fDgXUQm+AYFi4iIpVahRQ1gYGBeHh4kJmZWWx7ZmYm9evXL9drWywWLBYLVqu1XK8jrul8TVQexBHvWY84zy9g21x2bUkmbPxnENjEgIQiIlJSFdL85O3tTZs2bUhOTi7aZrPZSE5Opn379uV67bi4OHbt2sXWrVvL9TriPqx48FLhHdyX/zjH7P6EmX+FWR3hu4+NjiYiIhfhtKImJyeH1NTUohFMaWlppKamcvDgQQDi4+OZPXs28+fPZ/fu3YwdO5bc3Nyi0VAilU2KLZyeeQlssoZBQS58Ngo+Hwf5p42OJiIi5+G05qdt27YRHR1dtP7nyKOhQ4eSlJTEkCFDOHr0KJMnTyYjI4OIiAhWrFhxTudhZ1Pzk1yOI9Tm7oIn+KXL97D2Rfh2ARzaDoOToO51RscTEZG/MNntdrvRISrCnx2Fs7Ky8Pf3d/r5SzJ8WFzXged7wS9rYNEoyD0CXtWh1ysQcZfR0URE3Fppfn+7zLufRAx3TWcYu8HxvwWnYclYWDwG8vW2eBGRysDtixqLxUJYWBiRkZFGRxF34FcP7vkMov8NJjPs/BDe7gyZPxqdTESkynP7okajn8TpzB7Q6VEY+iXUvBKO/QSzu8D2+VA1WnNFRColty9qRMpN6C2OyfqaxEDhWfhyvGOEVN4po5OJiFRJLjWjcFlo9JOUqxqBcNcnsPF1Clc9jef3n/DLznWMKxjPLnvoOYcfeL5XxWcUEaki3P5JjZqfpNyZzXDLP7g9fzKH7FdwjTmDxd5TuMdjFaDmKBGRiuL2RY1IRdlhb0qvvGmssrbGx1TAs17zmOk1g5posj4RkYqgokbEiU5Sk1EF/+SZgrspsHvQ2+O/LPV+ghamX4yOJiLi9ty+qNGQbql4JuZYezE4fwq/2QO5ynyERd5TGOaxQqOjRETKkWYUdhLNKCzn408OL3rNprvH//fpur439JsJ1WobG0xExEVoRmGRSiIbP8YUPMKUgqHk2T1hz1JI7Ai/bTM6moiI21FRI1LuTMy3xjIw/ymofTVkHYS5sbDxDTVHiYg4kYoakQryg/0aGL0Wmg0AWyGs/Dd8eAecPm50NBERt+D2RY06Ckul4hsAg+ZBr+ng4QM/rYDEDnBws9HJRERcntsXNZp8TyodkwkiR8L9X0OdxpD9G8zrCeumg81mdDoREZfl9kWNSKV1ZUtHc1SLwWC3QvJU+GAw5B4zOpmIiEtSUSNiJJ+acNts6PsGePrCz19D4i1wYIPRyUREXI6KGhGjmUzQ+j4YtRoCr4NTv8P83rD2JbDpRawiIiWlokaksggKgwdWQ/hdYLfB6mdhwQA4lWl0MhERl6CiRqQy8a4BA96C/m+BV3VIW+tojvpljdHJREQqPRU1IpVRxF3wwBqoFwa5R+Dd/rB6mpqjREQuwtPoAOXNYrFgsViwWvXLQFxM3evg/mRY8TjseBfWvuDoQDzwHUKn7SjRKQ4836ucQ4qIVB5u/6RG89SIS/Ou7hgZdds74O0Hv66HxJvpaN5pdDIRkUrH7Z/UiFQmJXmb+3mfrrQcDMGt4JNhkPk973q/wJuFfXmlcDBWPJwfVETEBbn9kxoRtxHYxDEL8Y0jAXjQ8ws+8n6GK/nD4GAiIpWDihoRV+LlC72nE5c/nmx7NSLNP7HcZxJdzCXrYyMi4s5U1Ii4oGW2dvTOn8Z3tqupbcphrvfLPOH5Pl4UGh1NRMQw6lMjUsmUpN8NwEF7EIPyn2KS5wcM9/yKBzyXEWney0MFD/GbvW45pxQRqXxc5knN3r17iYiIKFqqVavGkiVLjI4lYqh8vJhaOJQH8v9Blr06rcw/s8x7EreaNdpPRKoel3lSc91115GamgpATk4OoaGhdOvWzdhQIpXESlsku/JDecPrDVqZf+Zt71eZVxgLhTHg6WN0PBGRCuEyT2r+6osvvqBr167UqFHD6CgilcZv9rrcnj+ZtwsdQ8KHe34Fc26F478YnExEpGI4rahJSUmhT58+BAcHYzKZzts0ZLFYCA0NxdfXl6ioKLZs2VKma3388ccMGTLkMhOLuJ8CPJlWeDcj8idwwu4Hv6fCrE7w42Kjo4mIlDunFTW5ubmEh4djsVjOu3/hwoXEx8czZcoUduzYQXh4OLGxsRw5cqTomIiICJo3b37Ocvjw4aJjsrOz2bhxIz179rxonry8PLKzs4stIlXFN7bW9MxLgJB2kJftmLRvaTwUnDU6mohIuTHZ7Xa7009qMrF48WL69+9ftC0qKorIyEhmzpwJgM1mIyQkhIceeoiJEyeW+NwLFizgq6++4r333rvocU899RRTp049Z3tWVhb+/v4lvl5JlXTEikhFOvDcrY4XYa6f7thQvwUMng9XNDY2mIhICWVnZxMQEFCi398V0qcmPz+f7du3ExMT878Lm83ExMSwadOmUp2rpE1PkyZNIisrq2hJT08vdW4RVxf6r5WEfn0j9+U/zjG7P2R8T86Mmxj/xL8InbhMxbiIuJUKKWqOHTuG1WolKCio2PagoCAyMjJKfJ6srCy2bNlCbGzsJY/18fHB39+fBQsW0K5dO7p27Vrq3CLuIsUWTs+8BDbbbsDPdJYZ3jNJ8JyND/lGRxMRcRqXGv0UEBBAZmYm3t7eJf6M3tIt4nCE2tyd/wSvF96GzW7iTs/VfO79JBzda3Q0ERGnqJCiJjAwEA8PDzIzM4ttz8zMpH79+uV6bYvFQlhYGJGRkeV6HRFXYMWDVwsHcU/BJI7aA7jenA5vd4bUD42OJiJy2SqkqPH29qZNmzYkJycXbbPZbCQnJ9O+fftyvbae1Iica6OtOT3ynme9tRkUnIYlY2DJg5Cfa3Q0EZEyc1pRk5OTQ2pqatGsv2lpaaSmpnLw4EEA4uPjmT17NvPnz2f37t2MHTuW3Nxchg8f7qwIIlIKxwjgvoJJEP1vMJkh9X14OxoydxkdTUSkTJz2moRt27YRHR1dtB4fHw/A0KFDSUpKYsiQIRw9epTJkyeTkZFBREQEK1asOKfzsLNZLBYsFgtWq7VcryPiimyYodOjcFV7WHQ/HNsLs7tAzxeh1b1gMhkdUUSkxMplnprKqDTj3MtCQ2PFVR143vFaBXKPwWcPwP7/byZucTv0ng4+NY0LJyJVXqWbp0ZEXECNQLj7U+g6BUwe8P3Hjk7EGd8bnUxEpETcvqjR6CeRUjCboUM8DF8O/g3gj59hdlfYOgeqxkNdEXFhan5yEjU/ibupxSle8Uqkq8e3jg3NBkCf18E3wNhgIlKlqPlJRC7bSWpyf8E/ebbgbjB7Ot70PasTHP7W6GgiIufl9kWNmp9Eys6OmXesvWD4CghoBCfSYM6t8N9Zao4SkUrH7YsaTb4n4gQhkTAmBa7vDdZ8+M9jsPAeOHPC6GQiIkWcNk+NiLiv//UZu5NhHnV4wvN9vPcsJX3XfxlX8BA77U3+NzRcRMQgbv+kRkScyUSStTsD86fyq60eIeajfOo9lZEey9QcJSKGc/uiRn1qRJzve/s19M6fxjJrW7xMVp70eh8+vBNOHzc6mohUYRrS7SQa0i1Vk517PL7mSc/38DEVgH9DGDQXGkUZHUxE3ISGdItIBTHxnrUbA/KnQp3GkP0bzOsB618Fm83ocCJSxaioEZHLtsseCqPXQovBYLfC10/BB4Md75MSEakgKmpExDl8asJts6HPDPD0hZ+/hsRb4MAGo5OJSBXh9kWNOgqLVCCTCdoMhVHfQGBTOPU7zO8Na18Cm9XodCLi5ty+qNHkeyIGCGoGD6yB8LvAboPVz8KCAZBzxOhkIuLG3L6oERGDeNeAAW9B/7fAqzqkrYW3boZf1hidTETclIoaESlfEXfBqNVQ9wbIPQLv9ofV09QcJSJOp6JGRMpfvesd/Wxa3wfYYe0L8G4/yP7d6GQi4kZU1IhIxfCuDn3fgNveAW8/OLDOMTrq56+NTiYibkJFjYhUrJaD4YG1ENQCTh+D9wbC11PBWmh0MhFxcSpqRKTiBTaB+7+GG0c41tdPdwz9zjpkbC4RcWmeRgcobxaLBYvFgtWqTokilYqXL/R+FUI7wBfj4eAmR3PUgERC5176qc2B53tVQEgRcSV6oaWT6IWWIpd2wULk+C/wyTD4fScAswp78VLhEAov8u8uFTUiVYNeaCkirqXONTByFbQdDcBoz2V87P00DThqcDARcSUqakSkcvD0gZ4vwpD3yLJXp7X5Z5b7TKKbeZvRyUTERaj5yUnU/CTiPA1NR5npNYMI834A5hZ2J6HwLgpK2Q1QTVQirk/NTyLi0n6z12Vw/hRmF/YEYITnCj71fooQU6bByUSkMlNRIyKVUgGePFd4DyPz/8kJux/h5l9Y5v0EPcz/NTqaiFRSLlXUvPrqqzRr1oywsDDGjx9PFWk5E6nSkm1t6JU3jW22pvibzvCW9+s87TkPH/KNjiYilYzLFDVHjx5l5syZbN++ne+//57t27ezefNmo2OJSAU4TCB35P+bNwv7AnCf5yo+855CqEnvjhKR/3GZogagsLCQs2fPUlBQQEFBAfXq1TM6kohUkEI8ebHwDobmP84f9po0M//KUu9/0de80ehoIlJJOK2oSUlJoU+fPgQHB2MymViyZMk5x1gsFkJDQ/H19SUqKootW7aU+Px169ZlwoQJNGrUiODgYGJiYmjcuLGz4ouIi1hrC6dnXgL/tV2Pn+ksM7xnMs1ztpqjRMR5RU1ubi7h4eFYLJbz7l+4cCHx8fFMmTKFHTt2EB4eTmxsLEeOHCk6JiIigubNm5+zHD58mBMnTrB06VIOHDjAoUOH2LhxIykpKRfMk5eXR3Z2drFFRNxDJnW4K/9fzCjsj81u4i7P1Xzu/SSNTXp3lEhVVi7z1JhMJhYvXkz//v2LtkVFRREZGcnMmTMBsNlshISE8NBDDzFx4sRLnvOTTz5hzZo1RUXTSy+9hN1u57HHHjvv8U899RRTp049Z7vmqRFxLzebv+c1rzepa8ritN2HfxcM5zNbR0Dz1Ii4g0o3T01+fj7bt28nJibmfxc2m4mJiWHTpk0lOkdISAgbN27k7NmzWK1W1qxZw3XXXXfB4ydNmkRWVlbRkp6eftnfQ0Qqnw22FvTMS2CDtRnVTXlM907kJc9EqnHW6GgiUsEqpKg5duwYVquVoKCgYtuDgoLIyMgo0TnatWtHz549adWqFS1btqRx48b07dv3gsf7+Pjg7+/PggULaNeuHV27dr2s7yAilddRanFvwSSmFwzCajcx2DOFL7yfhMxdRkcTkQrkUqOfnnvuOXbv3s2PP/7IjBkzMJlMl/xMXFwcu3btYuvWrRWQUESMYsPMDOtt3F3wLzLttbjWfAhmd4Ed74LmtBKpEiqkqAkMDMTDw4PMzOJTnGdmZlK/fv1yvbbFYiEsLIzIyMhyvY6IVA6bbWH0zEsgxdoCCs/AFw/BZw9A3imjo4lIOauQosbb25s2bdqQnJxctM1ms5GcnEz79u3L9dp6UiNS9fxBAEMLHoeuk8HkAd9/DG93hozvjY4mIuXIaUVNTk4OqamppKamApCWlkZqaioHDx4EID4+ntmzZzN//nx2797N2LFjyc3NZfjw4c6KICJSxI4ZOvwThi0D/wbwx88wuytsnaPmKBE35emsE23bto3o6Oii9fj4eACGDh1KUlISQ4YM4ejRo0yePJmMjAwiIiJYsWLFOZ2Hnc1isWCxWLBareV6HRGppK5qD6PXwZKxsO8rWBYPB9ZBnxng6/zpHUTEOOUyT01lVJpx7mWheWpEKp9i89TYbLBpJiRPBVsh1L4aBidBcIRR8USkBCrdPDUiIoYzm+Hm8TB8BQQ0ghNpMKcb/HeWmqNE3ITbFzUa/SQixYREwpgUuL43WPPhP4/Bx/fCmZNGJxORy+T2RY1GP4nIOarVhiHvQfcXwOwFu7+EWR3gt+1GJxORy+D2RY2IyHmZTNBuDIxcCbVD4eRBmHsrbJyp5igRF+X2RY2an0Tkohq0htEpENbP0YF45b/gwzvh9HGjk4lIKbl9UaPmJxG5JN8AGDwfer0CHj7w038gsQMc/K/RyUSkFNy+qBERKRGTCSLvh/u/hjqNIfs3mNcD1r/qGA4uIpWeihoRkb+6siWMXgvNB4HdCl8/BR/cDrnHjE4mIpfg9kWN+tSISKn51ISB70Cf18HTF35eBYm3wIENRicTkYtw+6JGfWpEpExMJmgzDEZ9A4FN4dTvML83rH0JbHrtikhl5LR3P4mIuKKSvOLkwNTVsHwC7PwQVj8Lv66H22aDX70KSCgiJeX2T2pERC6bjx8MSIR+b4JXdfhljaM56pe1RicTkb9QUSMiUlKt7oZRq6HuDZCTCe/2g9UJao4SqSTcvqhRR2ERcap61zv62bS6F7DD2ucdxc2pDKOTiVR5bl/UqKOwiDidd3XoN9PRr8arBhxYB2/dDD8nG51MpEpz+6JGRKTctLzd8YqFoBZw+hi8NxCSnwZrodHJRKokFTUiIpcjsIljFuIbRwJ2WPeKY+h31iGjk4lUOSpqREQul5cv9J4Og+aBd004uMkxOuqnlUYnE6lSVNSIiDhL89tgTApcGQ5njsMHg2Hlk2AtMDqZSJWgokZExJnqXAMjV0Hb0Y71jTMcL8Y8edDYXCJVgIoaERFn8/SBni/C7QvAJwB+2wqJHWDPpWcvFpGyc/uiRvPUiIhhwvo6mqMatIGzJ+Gju+A/E6Ew3+hkIm7J7YsazVMjIoaqHQrDV0D7cY71/74Fc2+F42mGxhJxR25f1IiIGM7TG2Kfgzs/At9acPhbmNURflxidDIRt6KiRkSkolzXA8ash5AoyMuGT4bCsn9CwVmjk4m4BRU1IiIVqVYIDFsGNz/iWN/6DszpBn/sNzSWiDtQUSMiUtE8vKDbVLj7U6h+BWR852iO+v5To5OJuDRPowOIiFRZ13ZzNEd9OhIOboRFIx0vx+z+PHhVu+hHQydeenj4ged7OSupiEtwqSc1L7/8Ms2aNaN58+a89957RscREbl8/sEw9Evo+Chggu1JMLsrHP3J6GQiLsdliprvv/+eDz74gO3bt7N161ZmzpzJyZMnjY4lInL5PDyhy7/h3s+gRl048iO83Rl2fmR0MhGX4jJFze7du2nfvj2+vr5Uq1aN8PBwVqxYYXQsERHnadzF0RwV2gEKcmHxaFgSB/m5RicTcQlOK2pSUlLo06cPwcHBmEwmlixZcs4xFouF0NBQfH19iYqKYsuWLSU+f/PmzVmzZg0nT57kxIkTrFmzhkOHDjkrvohI5VCzPtz3OXR+AkxmSH0PZneBI7uNTiZS6Tmto3Bubi7h4eGMGDGC22677Zz9CxcuJD4+nsTERKKionjttdeIjY1l79691KtXD4CIiAgKCwvP+ezKlSsJCwtj/PjxdOnShYCAANq1a4eHh8cF8+Tl5ZGXl1e0np2d7YRvKSJSAcwe0PlxuOomR+fho3vg7Wjo9TJE3A0mk9EJRSolpxU1PXr0oEePHhfcP336dEaNGsXw4cMBSExMZNmyZcydO5eJEycCkJqaetFrjB49mtGjHW++vf/++7n22msveGxCQgJTp04t5bcQEalEru4AYzbA4gdg/zfweRykrYNerxidTKRSqpA+Nfn5+Wzfvp2YmJj/XdhsJiYmhk2bNpX4PEeOHAFg7969bNmyhdjY2AseO2nSJLKysoqW9PT0sn8BERGj+NWFuxdBlycdzVHffQRvd+Z600Gjk4lUOhUyT82xY8ewWq0EBQUV2x4UFMSePXtKfJ5+/fqRlZVFjRo1mDdvHp6eF47v4+ODj48PFosFi8WC1Wotc34REUOZzdBxgqM56tOR8Mc+Pvd+kqcK7+NDaxdAzVEi4GKT75Xmqc6f4uLiiIuLIzs7m4CAgHJIJSJSQa66yTE6askYfPatJMFrDu3Nu3iiYCQ5VDc6nYjhKqT5KTAwEA8PDzIzM4ttz8zMpH79+uV6bYvFQlhYGJGRkeV6HRGRClHjCrhzIdMK7qTA7kFfj0186f0vmpnSjE4mYrgKKWq8vb1p06YNycnJRdtsNhvJycm0b9++XK8dFxfHrl272Lp1a7leR0SkwpjNvG3tw5D8J/nNHsjV5kw+857CfR5fAXaj04kYxmlFTU5ODqmpqUUjmNLS0khNTeXgQUdntvj4eGbPns38+fPZvXs3Y8eOJTc3t2g0lIiIlM4Oe1N65U1jlbUNPqZCnvaaz5ter+OPJuuTqslpfWq2bdtGdHR00Xp8fDwAQ4cOJSkpiSFDhnD06FEmT55MRkYGERERrFix4pzOw86mjsIi4s6y8GNUQTwjbCuY6PkBPT220NyUxriC8UZHE6lwJrvdXiWeVf7ZUTgrKwt/f3+nn78kb8wVkYpVkrdUu+rbrs+Xu6VpPzO9ZtDIfJR8uwfe3Z+FdmM1WZ+4tNL8/naZdz+JiMjFfWdvTO/8aSy3tsXbZIWvJsFHd8Hp40ZHE6kQbl/UaPSTiFQl2dTgwYKHebJgGHh4w97lMKsjpJf8XXsirsrtixqNfhKRqsfEAuutcP/XUOcayEqHeT1gw+tgsxkdTqTcuH1RIyJSZV0ZDg+sheYDwVYIqybDh0Mg9w+jk4mUC7cvatT8JCJVmq8/DJwDfV4HT1/YtxISb4FfNxqdTMTp3L6oUfOTiFR5JhO0GQb3J8MV18Kpw5DUG1JeVnOUuBW3L2pEROT/1W8OD6yBlneA3QrfPAPv3QY5R4xOJuIUKmpERKoSHz8YkAj9LOBZDX5Z7WiOSksxOpnIZXP7okZ9akRE/sZkglb3OJ7a1L0ecjLh3X6w5nmwafZ1cV1uX9SoT42IyAXUux5GrXYUOHYbrElwFDenMoxOJlImbl/UiIjIRXhXdzRFDXgbvGrAgXWO5qj93xidTKTUVNSIiAiED3E0RwU1h9yjsOA2SH4GrIVGJxMpMRU1IiLiULepYxbiNsMBO6x7Geb3gaxDRicTKRG3L2rUUVhEpBS8qkGf1xwT9nnXhIMbHc1RP600OpnIJbl9UaOOwiIiZdBiEIxeC/Vbwpnj8MFgWPkkWAuMTiZyQW5f1IiISBld0djRHNX2Acf6xhkwryecTDc2l8gFqKgREZEL8/SBni/B7e+CTwD8tsXRHLVnudHJRM6hokZERC4trJ+jOSq4NZw9CR/dCSuegMJ8o5OJFFFRIyIiJVPnahjxFbSLc6xvtsDcWDhxwNBYIn9SUSMiIiXn6Q3dp8EdH4JvLTi8AxI7wq4vjE4moqJGRETK4PqeMGYdNGwLeVnw8b2w/FEozDM6mVRhbl/UaJ4aEZFyUqsRDF8ONz/sWN/yNszpBn/sNzaXVFluX9RonhoRkXLk4QXdnoa7P4VqdeD3nTCrE/ywyOhkUgW5fVEjIiIV4NpuMGY9NLoJ8k/BpyPgy0eg4IzRyaQKUVEjIiLOEdAAhn4JHSYAJtg+D96JgWP7jE4mVYSKGhERcR4PT+j6JNz7GdSoC5k/OJqjdi40OplUASpqRETE+Rp3cTRHhXaAglxY/AB8Hgf5p41OJm5MRY2IiJSPmvXhvs+h8yTABN++B7O7wJE9RicTN1Upi5oBAwZQu3ZtBg0adM6+pUuXct1113HttdfyzjvvGJBORERKzOwBnSfC0C/ALwiO7oa3O8O37xudTNxQpSxqHn74Yd59991zthcWFhIfH88333zDt99+y0svvcQff/xhQEIRESmVqzvCmA1wTTQUnoHPH4TFYyAvx+hk4kYqZVHTuXNnatasec72LVu20KxZMxo0aICfnx89evRg5cqVBiQUEZFS86sL93wGXZ4Ekxl2fgizoyHzR6OTiZsodVGTkpJCnz59CA4OxmQysWTJknOOsVgshIaG4uvrS1RUFFu2bHFGVg4fPkyDBg2K1hs0aMChQ4eccm4REakAZjN0nADDlkHNYDj2k6OfzfYksNuNTicurtRFTW5uLuHh4VgslvPuX7hwIfHx8UyZMoUdO3YQHh5ObGwsR44cKTomIiKC5s2bn7McPny47N/kb/Ly8sjOzi62iIhIJXHVTY7RUU26QeFZ+PJhWDQSzurvaik7z9J+oEePHvTo0eOC+6dPn86oUaMYPnw4AImJiSxbtoy5c+cyceJEAFJTU8sUNjg4uNiTmUOHDtG2bdvzHpuQkMDUqVPLdB0REakANa6Auz6GTW/A11Mdr1Y4/C0MToIrw41OJy7IqX1q8vPz2b59OzExMf+7gNlMTEwMmzZtuuzzt23blh9++IFDhw6Rk5PDf/7zH2JjY8977KRJk8jKyipa0tPTL/v6IiLiZGaz44WYI1ZAQAgc/8UxC/GW2WqOklJzalFz7NgxrFYrQUFBxbYHBQWRkZFR4vPExMQwePBgli9fTsOGDYsKIk9PT1555RWio6OJiIjgn//8J1dcccV5z+Hj44O/vz8LFiygXbt2dO3atexfTEREyldIWxidAtf1Ams+LJ8AnwyFs1lGJxMXUurmp4rw9ddfX3Bf37596du3b4nPFRcXR1xcHNnZ2QQEBDgjnoiIlIfqdeCO92HzW7BqMuz6HA6nwuB50KCN0enEBTj1SU1gYCAeHh5kZmYW256ZmUn9+vWdeakSs1gshIWFERkZacj1RUSkFEwmaP8gjPwKajWCk7/CnFhHoaPmKLkEpxY13t7etGnThuTk5KJtNpuN5ORk2rdv78xLlVhcXBy7du1i69athlxfRETKoEEbGL0ObugDtgJYMRE+uhtOHzc6mVRipS5qcnJySE1NLRrBlJaWRmpqKgcPHgQgPj6e2bNnM3/+fHbv3s3YsWPJzc0tGg0lIiJSItVqwe0LoOfL4OENe5fBrI6Qrn+kyvmVuk/Ntm3biI6OLlqPj48HYOjQoSQlJTFkyBCOHj3K5MmTycjIICIighUrVpzTebiiWCwWLBYLVqvVkOuLiMhlMJmg7ShoGAmfDIMTaTCvO3SdDO0fcoyeEvl/Jru9ajRS/tlROCsrC39/f6efP3TiMqefU0Quz4Hne13ymJL87JbkPBXNVXNflrPZjkn6fvzMsX7trdA/0THfjbit0vz+VokrIiKuwdcfBs2F3q+Chw/sWwmJt8CvG41OJpWE2xc1Gv0kIuJGTCa4cQSMSoYrmsCpw5DUG1JeBpvN6HRiMLcvajT6SUTEDdVvAQ+shZZDwG6Fb56B9wdCzlGjk4mB3L6oERERN+XjBwNmQd+Z4FkN9n8DiTdDWorRycQgbl/UqPlJRMSNmUzQ+l54YDXUvR5yMuHdfrDmebBp1GtV4/ZFjZqfRESqgHo3wKhvIOIesNtgTYKjuDlV8vcOiutz+6JGRESqCO8a0N/iaJLyqgEH1jlGR+3/xuhkUkFU1IiIiHsJvwMeWAP1mkHuUVhwGyQ/A9ZCo5NJOXP7okZ9akREqqC6TR3DvtsMA+yw7mWY3weyDxudTMqR2xc16lMjIlJFeVWDPq/DwDng7QcHNzqao/atMjqZlBO3L2pERKSKazEIRqdA/ZZw+g94fxCsmgzWAqOTiZOpqBEREfd3RWMYuQoiRznWN7wOSb3gZLqxucSpVNSIiEjV4OULvV6G298FnwBI/6+jOWrPcqOTiZO4fVGjjsIiIlJMWD8YvRaCW8PZk/DRnfDVv6Aw3+hkcpncvqhRR2ERETlHnathxFfQ7kHH+qaZMK87nDhgaCy5PG5f1IiIiJyXpzd0T4A7PgTfWnBoOyR2hF1fGJ1MykhFjYiIVG3X94Qx66BhJORlwcf3wvJHoTDP6GRSSipqREREajWC4f+Bmx92rG95G+Z0gz/2G5tLSkVFjYiICICHF3R7Gu76BKrVgd93wqxO8MNnRieTElJRIyIi8ldNb4Ux66FRe8g/BZ8Oh6X/gIIzRieTS1BRIyIi8ncBDWDoUujwT8AE2+bCOzFwbJ/RyeQi3L6o0Tw1IiJSJh6e0HUy3LMIqgdC5g+O5qjvPjY6mVyA2xc1mqdGREQuS5OuMHYDhHaAglz4bBR8Pg7yTxudTP7G7YsaERGRy1azPtz3OXSaCJjg2wUwuwsc2WN0MvkLFTUiIiIlYfaA6EmO4sYvCI7uhtnR8O37RieT/6eiRkREpDSu6eQYHXVNZyg4DZ8/CIvHQF6O0cmqPBU1IiIipeVXD+5ZDF3+DSYz7PzQ8dQm80ejk1VpKmpERETKwmyGjo86hn7XvBKO/eToZ7N9PtjtRqerkiplUTNgwABq167NoEGDSrVPRESkwoXe7GiOatINCs/Cl+Nh0f2Qd8roZFVOpSxqHn74Yd59991S7xMRETFEjUC462OImQomD/jhU5jV0fGqBakwlbKo6dy5MzVr1iz1PhEREcOYzXDLI44XY/o3hOO/wDvdYMtsNUdVkFIXNSkpKfTp04fg4GBMJhNLliw55xiLxUJoaCi+vr5ERUWxZcsWZ2QVERGp/BpFwZh10LQHWPNg+QT4ZCiczTI6mdsrdVGTm5tLeHg4FovlvPsXLlxIfHw8U6ZMYceOHYSHhxMbG8uRI0eKjomIiKB58+bnLIcPHy77N/mbvLw8srOziy0iIiIVonoduPNDiJ0GZk/Y9TkkdoBDO4xO5tY8S/uBHj160KNHjwvunz59OqNGjWL48OEAJCYmsmzZMubOncvEiRMBSE1NLVvaUkhISGDq1Knlfh0REZHzMpmgfRyEtINPh8HJX2HOrXDrMxA1xrFfnMqpfWry8/PZvn07MTEx/7uA2UxMTAybNm1y5qUuadKkSWRlZRUt6enpFXp9ERERABq2gdHr4PreYCuAFRNh4T1w5oTRydyOU4uaY8eOYbVaCQoKKrY9KCiIjIyMEp8nJiaGwYMHs3z5cho2bFisILrYvr/y8fHB39+fBQsW0K5dO7p27Vq2LyUiInK5qtWCIe9Bj5fAwxv2LIXEjpCuly07U6mbnyrC119/XaZ95xMXF0dcXBzZ2dkEBARcbjQREZGyMZkg6gEIiYRPhsOJNJjXHbpOgfbjHKOn5LI49U8wMDAQDw8PMjMzi23PzMykfv36zrxUiVksFsLCwoiMjDTk+iIiIsUEt4LRa6HZALAVwqon4cM74PRxo5O5PKcWNd7e3rRp04bk5OSibTabjeTkZNq3b+/MS5VYXFwcu3btYutWPeITEZFKwjcABs2D3q+Chw/s+woSb4FfK7b/qbspdVGTk5NDampq0QimtLQ0UlNTOXjwIADx8fHMnj2b+fPns3v3bsaOHUtubm7RaCgRERHB0Rx14wgYlQxXNIHsQ5DUC9a9Ajab0elcUqn71Gzbto3o6Oii9fj4eACGDh1KUlISQ4YM4ejRo0yePJmMjAwiIiJYsWLFOZ2HK4rFYsFisWC1Wg25voiIyEXVbwEPrIGl8fD9x5D8NBzYAANmgV9do9O5FJPdXjXmbv6zo3BWVhb+/v5OP3/oxGVOP6eIXJ4Dz/e65DEl+dktyXkqmqvmlouw2+Hb92D5o1B4Bvzqw6A5EHqL0ckMVZrf3+pqLSIiUhmYTND6XnhgNQReBzkZML8PrHkBbGptKAm3L2o0+klERFxKvRschU3EPWC3wZppsGAAnMq89GerOLcvajT6SUREXI53DehvcfSr8aoOaWsdo6N+WWN0skrN7YsaERERlxV+BzywFuo1g9wj8G5/+OY5sBYanaxScvuiRs1PIiLi0uo2dQz7bj0UsEPKi/BuX8j+3ehklY7bFzVqfhIREZfnVQ36zoCBc8DbD37dAIk3w77SvTrI3bl9USMiIuI2WgyC0SmOuW1O/wHvD4SvnwJrgdHJKgUVNSIiIq7kisYw8muIvN+xvv5Vx0zEWb8Zm6sScPuiRn1qRETE7Xj5Qq9XYPB88PGH9P86RkftXWF0MkO5fVGjPjUiIuK2mvV3NEcFt4IzJ+DDIfDVv6Aw3+hkhnD7okZERMSt1bkaRnwFUWMd65tmwrzucOJXY3MZQEWNiIiIq/P0gR7Pw5D3wTcADm2HWR1g95dGJ6tQKmpERETcxQ29Ycx6aHAjnM2ChffA8segMM/oZBXC7YsadRQWEZEqpVYjGLECbnrIsb5lFsy5FY7/YmyuCuD2RY06CouISJXj4QW3Pgt3fQzVasPvqTCrE/y42Ohk5crtixoREZEqq2msozkqpB3kZcMnw2BpPBScNTpZuVBRIyIi4s4CGsKwZXBLvGN92xx4JwaO/WxsrnKgokZERMTdeXhCzBS4ZxFUD4TM7+HtTvDdJ0YncyoVNSIiIlVFkxhHc1RoB8jPgc/uh8/HQf5po5M5hYoaERGRqsT/Srjvc+j0OGCCbxfAO13h6F6jk102FTUiIiJVjdkDop9wFDc16sGRXfB2Z0j9wOhkl8XtixrNUyMiInIB13SCsRvgms5QcBqWjIXFYyA/1+hkZeL2RY3mqREREbkIv3pwz2cQ/W8wmWHnh46nNpk/Gp2s1Ny+qBEREZFLMHtAp0dh6JdQ80o49hPM7gLb54PdbnS6ElNRIyIiIg6htzhGRzWJgcKz8OV4WHQ/5J0yOlmJqKgRERGR/6kRCHd9AjFPgckDfvjU8YqF378zOtklqagRERGR4sxmuOUfMHw5+DeA4/sdsxBvfadSN0epqBEREZHza9TO0RzVtAdY82DZPx3vjzqbZXSy86qURc2AAQOoXbs2gwYNKrY9PT2dzp07ExYWRsuWLfnkE/ea3llERKTSqV4H7vwQbn0OzJ6wawnM6giHdhid7ByVsqh5+OGHeffdd8/Z7unpyWuvvcauXbtYuXIljzzyCLm5rjmWXkRExGWYTHDTOBjxFQQ0ghMHYM6tsDmxUjVHVcqipnPnztSsWfOc7VdeeSUREREA1K9fn8DAQI4fP17B6URERKqohjfCmBS4vjfYCmDF47DwHjhzwuhkQBmKmpSUFPr06UNwcDAmk4klS5acc4zFYiE0NBRfX1+ioqLYsmWLM7IWs337dqxWKyEhIU4/t4iIiFxAtdow5D3o8SJ4eMOepZDYEX7bZnSy0hc1ubm5hIeHY7FYzrt/4cKFxMfHM2XKFHbs2EF4eDixsbEcOXKk6JiIiAiaN29+znL48OESZTh+/Dj33Xcfb7/99gWPycvLIzs7u9giIiIiTmAyQdRoGLkSaodC1kGYGwsb3zC0OcqztB/o0aMHPXr0uOD+6dOnM2rUKIYPHw5AYmIiy5YtY+7cuUycOBGA1NTUsqXFUaz079+fiRMnctNNN13wuISEBKZOnVrm64iIiMglBLeC0Snw5cPw42JHH5tW90K1WobEcWqfmvz8fLZv305MTMz/LmA2ExMTw6ZNmy77/Ha7nWHDhtGlSxfuvffeix47adIksrKyipb09PTLvr6IiIj8jW8ADJoHvabDoDmGFTTg5KLm2LFjWK1WgoKCim0PCgoiIyOjxOeJiYlh8ODBLF++nIYNGxYVRBs2bGDhwoUsWbKEiIgIIiIi+P777897Dh8fH/z9/VmwYAHt2rWja9euZf9iIiIicmEmE0SOdMxrY6BSNz9VhK+//vq822+55RZsNlupzhUXF0dcXBzZ2dkEBAQ4I56IiIhUQk59UhMYGIiHhweZmZnFtmdmZlK/fn1nXqrELBYLYWFhREZGGnJ9ERERqRhOLWq8vb1p06YNycnJRdtsNhvJycm0b9/emZcqsbi4OHbt2sXWrVsNub6IiIhUjFI3P+Xk5PDzzz8XraelpZGamkqdOnVo1KgR8fHxDB06lBtvvJG2bdvy2muvkZubWzQaSkRERKQ8lLqo2bZtG9HR0UXr8fHxAAwdOpSkpCSGDBnC0aNHmTx5MhkZGURERLBixYpzOg9XFIvFgsViwWq1GnJ9ERERqRilLmo6d+6M/RIT64wbN45x48aVOZQzqaOwiIhI1VAp3/0kIiIiUlpuX9Ro9JOIiEjV4PZFjUY/iYiIVA1uX9SIiIhI1eD2RY2an0RERKoGty9q1PwkIiJSNbh9USMiIiJVQ6V8oWV5+HNunezs7HI5vy3vdLmcV0TKriQ/7yX52S2vvzcuh6vmFimtP/87vtQceQAme0mOcmF/ziicn5/P/v37jY4jIiIiZZCenk7Dhg0veozbFzV/stlsHD58mJo1a2IymYq2R0ZGXrC/zYX2/X17dnY2ISEhpKen4+/v7/zwpXCx71OR5yvN50pybFnu04X2nW+bu95DV7h/F9uvn0HdQyNUxXtYmX8X2u12Tp06RXBwMGbzxXvNVJnmJ7PZfN4Kz8PD44J/+Bfad6Ht/v7+hv8wXuz7VOT5SvO5khxblvt0oX0XO97d7qEr3L+L7dfPoO6hEariPazsvwtL+pqjKt9ROC4urtT7LvYZozk7W1nPV5rPleTYstynC+2rzPcPnJvPFe7fxfbrZ1D30AhV8R66y+/CKtP8VJ7+fFlmVlaW4f/CkLLRPXRtun+uT/fQ9VWGe1jln9Q4g4+PD1OmTMHHx8foKFJGuoeuTffP9ekeur7KcA/1pEZERETcgp7UiIiIiFtQUSMiIiJuQUWNiIiIuAUVNSIiIuIWVNSIiIiIW1BRUwFOnz7NVVddxYQJE4yOImUQGhpKy5YtiYiIIDo62ug4UgZpaWlER0cTFhZGixYtyM3NNTqSlNDevXuJiIgoWqpVq8aSJUuMjiWl9Oqrr9KsWTPCwsIYP358iV5OWRZV5jUJRnruuedo166d0THkMmzcuBE/Pz+jY0gZDRs2jGeffZYOHTpw/PhxzYXiQq677jpSU1MByMnJITQ0lG7duhkbSkrl6NGjzJw5kx9//BEvLy86duzI5s2bad++vdOvpSc15Wzfvn3s2bOHHj16GB1FpEr68y/SDh06AFCnTh08PfXvOVf0xRdf0LVrV2rUqGF0FCmlwsJCzp49S0FBAQUFBdSrV69crqOi5iJSUlLo06cPwcHBmEym8z7ytFgshIaG4uvrS1RUFFu2bCm2f8KECSQkJFRQYvk7Z9xDk8lEp06diIyM5P3336+g5PKny72H+/btw8/Pjz59+tC6dWumTZtWgenFGT+Df/r4448ZMmRIOSeWv7vce1i3bl0mTJhAo0aNCA4OJiYmhsaNG5dLVhU1F5Gbm0t4eDgWi+W8+xcuXEh8fDxTpkxhx44dhIeHExsby5EjRwD4/PPPadq0KU2bNq3I2PIXl3sPAdavX8/27dv54osvmDZtGt99911FxRcu/x4WFhaybt063nzzTTZt2sSqVatYtWpVRX6FKs0ZP4PgeK/Qxo0b6dmzZ0XElr+43Ht44sQJli5dyoEDBzh06BAbN24kJSWlfMLapUQA++LFi4tta9u2rT0uLq5o3Wq12oODg+0JCQl2u91unzhxor1hw4b2q666yn7FFVfY/f397VOnTq3I2PIXZbmHfzdhwgT7vHnzyjGlXExZ7uHGjRvtt956a9H+F1980f7iiy9WSF4p7nJ+Bt9991373XffXREx5SLKcg8//vhj+4MPPli0/8UXX7S/8MIL5ZJPT2rKKD8/n+3btxMTE1O0zWw2ExMTw6ZNmwBISEggPT2dAwcO8PLLLzNq1CgmT55sVGT5m5Lcw9zcXE6dOgU4Oil+8803NGvWzJC8cq6S3MPIyEiOHDnCiRMnsNlspKSkcMMNNxgVWf6iJPfvT2p6qpxKcg9DQkLYuHEjZ8+exWq1smbNGq677rpyyaPecmV07NgxrFYrQUFBxbYHBQWxZ88eg1JJaZTkHmZmZjJgwAAArFYro0aNIjIyssKzyvmV5B56enoybdo0OnbsiN1u59Zbb6V3795GxJW/Kenfo1lZWWzZsoVFixZVdES5hJLcw3bt2tGzZ09atWqF2Wyma9eu9O3bt1zyqKipIMOGDTM6gpTBNddcw86dO42OIZepR48eGoHowgICAsjMzDQ6hlyG5557jueee67cr6PmpzIKDAzEw8PjnB+0zMxM6tevb1AqKQ3dQ9ene+jadP9cX2W7hypqysjb25s2bdqQnJxctM1ms5GcnFwuEwqJ8+keuj7dQ9em++f6Kts9VPPTReTk5PDzzz8XraelpZGamkqdOnVo1KgR8fHxDB06lBtvvJG2bdvy2muvkZuby/Dhww1MLX+le+j6dA9dm+6f63Ope1guY6rcxOrVq+3AOcvQoUOLjnnjjTfsjRo1snt7e9vbtm1r37x5s3GB5Ry6h65P99C16f65Ple6hya7vZzeKiUiIiJSgdSnRkRERNyCihoRERFxCypqRERExC2oqBERERG3oKJGRERE3IKKGhEREXELKmpERETELaioEREREbegokZERETcgooaEamUhg0bhslkwmQysWTJEsNyrFmzpihH//79DcshIpemokZEyt1fC5S/Lt27d7/o57p3787vv/9Ojx49im1fvXo1vXv3pm7duvj6+tK4cWOGDBlCSkpKiTO1aNGCMWPGnHffggUL8PHx4dixY9x00038/vvv3H777SU+t4gYQ0WNiFSIPwuUvy4ffvjhRT/j4+ND/fr18fHxKdr25ptv0rVrV6644goWLlzI3r17Wbx4MTfddBP/+Mc/Spxn5MiRfPTRR5w5c+acffPmzaNv374EBgbi7e1N/fr1qVatWsm/rIgYQkWNiFSIPwuUvy61a9cu1TkOHjzII488wiOPPML8+fPp0qULV111FS1btuThhx9m27ZtxY5fv349HTp0oFq1aoSEhDB+/Hhyc3MBuOeeezhz5gyLFi0q9pm0tDTWrFnDyJEjL+8Li0iFU1EjIi5j0aJFFBQU8Nhjj513v8lkKvr/+/fvp3v37gwcOJDvvvuOhQsXsn79esaNGwdAYGAg/fr1Y+7cucXOkZSURMOGDbn11lvL74uISLlQUSMiFWLp0qX4+fkVW6ZNm1aqc/z000/4+/tTv379om2LFi0qds7vv/8egISEBO6++24eeeQRrr32Wm666SZmzJjBu+++y9mzZwFHE9SaNWtIS0sDwG63M3/+fIYOHYrZrL8eRVyNp9EBRKRqiI6O5q233iq2rU6dOqU+z1+fxgDExsaSmprKoUOH6Ny5M1arFYCdO3fy3Xff8f777xcda7fbsdlspKWlccMNN9CtWzcaNmzIvHnzePrpp0lOTubgwYMMHz68DN9QRIymokZEKkSNGjVo0qTJZZ3j2muvJSsri4yMjKKnNX5+fjRp0gRPz+J/neXk5DB69GjGjx9/znkaNWoEgNlsZtiwYcyfP5+nnnqKefPmER0dzTXXXHNZOUXEGHq+KiIuY9CgQXh5efHCCy9c8tjWrVuza9cumjRpcs7i7e1ddNzw4cNJT0/ns88+Y/HixeogLOLC9KRGRCpEXl4eGRkZxbZ5enoSGBhY4nM0atSIV155hYcffpjjx48zbNgwrr76ao4fP857770HgIeHBwCPP/447dq1Y9y4cdx///3UqFGDXbt2sWrVKmbOnFl0zquvvpouXbrwwAMP4OPjw2233eaEbysiRtCTGhGpECtWrODKK68sttxyyy2lPs9DDz3EypUrOXr0KIMGDeLaa6+lZ8+epKWlsWLFClq0aAFAy5YtWbt2LT/99BMdOnSgVatWTJ48meDg4HPOOXLkSE6cOMFdd92Fr6/vZX9XETGGyW63240OISLyd8OGDePkyZOGviLhrypbHhE5l57UiEil9ecw8KVLlxqWYd26dfj5+RUbRSUilZOe1IhIpXTkyBGys7MBuPLKK6lRo4YhOc6cOcOhQ4cAx0irv86RIyKVi4oaERERcQtqfhIRERG3oKJGRERE3IKKGhEREXELKmpERETELaioEREREbegokZERETcgooaERERcQsqakRERMQt/B/AiIRVDX0WmgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "samples = power_law.sample(1000)\n", "\n", "fig, ax = plt.subplots()\n", "ax.hist(samples, bins=energies, density=True, label=\"Samples\")\n", "ax.plot(energies, [power_law.spectrum(e) / total_flux for e in energies], label=\"Model\")\n", "ax.set_xscale(\"log\")\n", "ax.set_yscale(\"log\")\n", "ax.set_xlabel(\"E [GeV]\")\n", "ax.legend()" ] }, { "cell_type": "markdown", "id": "3a750ec9", "metadata": {}, "source": [ "The `BrokenPowerLaw` class is also available and behaves in a very similar way." ] }, { "cell_type": "markdown", "id": "984c6a57", "metadata": {}, "source": [ "## Diffuse and point sources\n", "\n", "Once the spectral shape is defined, we can specify either a `DiffuseSource` or a `PointSource`. It is assumed that diffuse sources are isotropic and the flux model describes the per-steradian flux over the entire $4\\pi$ sky. We also specify a redshift of the source such that adiabatic neutrino energy losses can be accounted for. Naturally, `PointSource` objects also have a direction specified in (ra, dec) coordinates." ] }, { "cell_type": "code", "execution_count": 7, "id": "184ea965", "metadata": { "execution": { "iopub.execute_input": "2024-03-22T09:57:03.297545Z", "iopub.status.busy": "2024-03-22T09:57:03.297028Z", "iopub.status.idle": "2024-03-22T09:57:03.300908Z", "shell.execute_reply": "2024-03-22T09:57:03.300349Z" } }, "outputs": [], "source": [ "diffuse_source = DiffuseSource(power_law, z=0.0)\n", "\n", "ra = np.deg2rad(50)\n", "dec = np.deg2rad(-10)\n", "point_source = PointSource(power_law, z=0.5, coord=(ra, dec))" ] }, { "cell_type": "markdown", "id": "be206106", "metadata": {}, "source": [ "The original flux model can now be accessed from within the source along with its other properties:" ] }, { "cell_type": "code", "execution_count": 8, "id": "4bf3ae48", "metadata": { "execution": { "iopub.execute_input": "2024-03-22T09:57:03.303422Z", "iopub.status.busy": "2024-03-22T09:57:03.303036Z", "iopub.status.idle": "2024-03-22T09:57:03.307224Z", "shell.execute_reply": "2024-03-22T09:57:03.306541Z" } }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "diffuse_source.flux_model" ] }, { "cell_type": "markdown", "id": "f36b54d2", "metadata": {}, "source": [ "Sources and lists of sources can be used as input to simulations, as demonstrated in the simulation notebook." ] }, { "cell_type": "code", "execution_count": null, "id": "fb8ae56b", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "icecube_dev", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.18" } }, "nbformat": 4, "nbformat_minor": 5 }