{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6e68d5e0",
   "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": "66ae85a1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-08T10:36:23.902432Z",
     "iopub.status.busy": "2024-11-08T10:36:23.902234Z",
     "iopub.status.idle": "2024-11-08T10:36:24.604132Z",
     "shell.execute_reply": "2024-11-08T10:36:24.603544Z"
    }
   },
   "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": "b1ffdd62",
   "metadata": {},
   "source": [
    "## Spectral shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eabff136",
   "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": "1e044f9c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-08T10:36:24.606503Z",
     "iopub.status.busy": "2024-11-08T10:36:24.606225Z",
     "iopub.status.idle": "2024-11-08T10:36:24.609941Z",
     "shell.execute_reply": "2024-11-08T10:36:24.609365Z"
    }
   },
   "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": "8da1435b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-08T10:36:24.611832Z",
     "iopub.status.busy": "2024-11-08T10:36:24.611478Z",
     "iopub.status.idle": "2024-11-08T10:36:25.263893Z",
     "shell.execute_reply": "2024-11-08T10:36:25.263209Z"
    }
   },
   "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": "iVBORw0KGgoAAAANSUhEUgAAAksAAAG1CAYAAADpzbD2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABwGUlEQVR4nO3deVxU9f4/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": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "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": "1373b1e5",
   "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": "85223ef8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-08T10:36:25.266084Z",
     "iopub.status.busy": "2024-11-08T10:36:25.265684Z",
     "iopub.status.idle": "2024-11-08T10:36:25.270097Z",
     "shell.execute_reply": "2024-11-08T10:36:25.269554Z"
    }
   },
   "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": "2827bf3a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-08T10:36:25.271874Z",
     "iopub.status.busy": "2024-11-08T10:36:25.271682Z",
     "iopub.status.idle": "2024-11-08T10:36:25.275886Z",
     "shell.execute_reply": "2024-11-08T10:36:25.275368Z"
    }
   },
   "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": "30eeeb49",
   "metadata": {},
   "source": [
    "Sampling from the power law shape is also possible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2b8f262c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-08T10:36:25.277713Z",
     "iopub.status.busy": "2024-11-08T10:36:25.277519Z",
     "iopub.status.idle": "2024-11-08T10:36:25.965773Z",
     "shell.execute_reply": "2024-11-08T10:36:25.965051Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fe3857b5670>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAG1CAYAAAABTQXdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXGElEQVR4nO3deVxU9f7H8dfMsKmISyqGYZRmhSmkIVq5JYWYmKblvS0X0SwNsy63RbulWaat5vU6RZlbtlmZ5i0zjSLM5ecWba4ZJrngDoHKDDPz+2Nu3MgNcODMDO/n4zGPR3POmXPe4wn5eL6byeVyuRARERHxcWajA4iIiIh4gooaERER8QsqakRERMQvqKgRERERv6CiRkRERPyCihoRERHxCypqRERExC+oqBERERG/EGB0gJridDrZs2cP9evXx2QyGR1HREREKsDlcvHbb78RERGB2XzmZzG1pqjZs2cPkZGRRscQERGRKsjLy+OCCy444zG1pqipX78+4P5DCQsLMziNiIiIVERhYSGRkZFlv8fPpNYUNb83OYWFhamoERER8TEV6TqijsIiIiLiF/y+qLFarURHRxMXF2d0FBEREalGJpfL5TI6RE0oLCykQYMGFBQUqPlJRKSWcDgc2O12o2PIWQQFBZ12ZFNlfn/Xmj41IiJSe7hcLvbt28fRo0eNjiIVYDabueiiiwgKCjqn86ioERERv/N7QdOsWTPq1q2r+cm82O/zyO3du5eWLVue071SUSMiIn7F4XCUFTTnnXee0XGkApo2bcqePXsoLS0lMDCwyufx+47CIiJSu/zeh6Zu3boGJ5GK+r3ZyeFwnNN5VNSIiIhfUpOT7/DUvfK5oubYsWNceOGFPPjgg0ZHERERES/ic0XN008/TefOnY2OISIiIl7GpzoKb9++nS1btpCcnMwPP/xgdBwREfExUWM+qbFr7Xzmxhq7Vk0wmUwsXLiQ/v37Gx3ltDz2pCY7O5vk5GQiIiIwmUwsWrTopGOsVitRUVGEhIQQHx/P2rVrK3WNBx98kMmTJ3sosQc5NLGTiIicuwMHDjBy5EhatmxJcHAwzZs3JzExkZUrVxodzSd47ElNcXExMTExDB06lJtvvvmk/fPnzyc9PZ2MjAzi4+OZOnUqiYmJbN26lWbNmgEQGxtLaWnpSZ9dtmwZ69ato02bNrRp04ZVq1adNU9JSQklJSVl7wsLC8/h251BqQ1mJ8HFPaDHWLD41MMvERHxIgMHDsRmszF37lwuvvhi8vPzyczM5NChQ0ZH8wkee1KTlJTExIkTGTBgwCn3T5kyheHDh5Oamkp0dDQZGRnUrVuXWbNmlR2Tk5PDDz/8cNIrIiKCNWvW8O677xIVFcWDDz7IjBkzePLJJ0+bZ/LkyTRo0KDsFRkZ6amvWt7WJbB7Pax4AeYmQ8Hu6rmOiIj4taNHj7JixQqeffZZevbsyYUXXkinTp0YO3Ys/fr1A9y/S9u1a0e9evWIjIzk3nvvpaioqOwcc+bMoWHDhnz88cdceuml1K1bl0GDBnHs2DHmzp1LVFQUjRo1YvTo0eWGT0dFRfHUU0/x17/+lXr16tGiRQusVusZ8+bl5XHrrbfSsGFDGjduzE033cTOnTvL9mdlZdGpUyfq1atHw4YNueaaa/jll188+4f2JzXSUdhms7FhwwYSEhL+d2GzmYSEBFavXl2hc0yePJm8vDx27tzJCy+8wPDhwxk3btxpjx87diwFBQVlr7y8vHP+HqfUtj8MnAlB9WHXKsi4FrYtq55riYiI3woNDSU0NJRFixaVa2n4I7PZzLRp0/jxxx+ZO3cuX3zxBQ8//HC5Y44dO8a0adN49913Wbp0KVlZWQwYMIAlS5awZMkS5s2bx6uvvsoHH3xQ7nPPP/88MTExfPPNN4wZM4b777+f5cuXnzKH3W4nMTGR+vXrs2LFClauXEloaCi9e/fGZrNRWlpK//796d69O9999x2rV6/m7rvvrvZh9jXSVnLw4EEcDgfh4eHltoeHh7Nly5ZquWZwcDDBwcFYrVasVus5T+hzRu0GQcSV8P4Q2PcdvH0LXD0aeo0DS9VnRhQRkdojICCAOXPmMHz4cDIyMujQoQPdu3fnL3/5C+3btwfggQceKDs+KiqKiRMnMmLECF5++eWy7Xa7nVdeeYVWrVoBMGjQIObNm0d+fj6hoaFER0fTs2dPvvzySwYPHlz2uWuuuYYxY8YA0KZNG1auXMlLL73E9ddff1LW+fPn43Q6ef3118sKldmzZ9OwYUOysrK46qqrKCgooG/fvmU5Lr/8cs/+gZ2Czw3pBhgyZAgvvPBChY5NS0tj06ZNrFu3rnpDndcKhi2HTne736+aBrP7wNFqekIkIiJ+Z+DAgezZs4fFixfTu3dvsrKy6NChA3PmzAHg888/p1evXrRo0YL69etz5513cujQIY4dO1Z2jrp165YVEuB+gBAVFUVoaGi5bfv37y937S5dupz0fvPmzafM+e233/LTTz9Rv379sidMjRs35sSJE+zYsYPGjRszZMgQEhMTSU5O5l//+hd79+491z+es6qRoqZJkyZYLBby8/PLbc/Pz6d58+bVem2r1Up0dDRxcXHVep2oMZ8Q9XgmUdk9GGF7gEJXXfh1LUdfiueuR5+s0WGEIiLiu0JCQrj++ut5/PHHWbVqFUOGDGH8+PHs3LmTvn370r59exYsWMCGDRvK+r3YbLayz/957SSTyXTKbU6ns8oZi4qK6NixIzk5OeVe27Zt47bbbgPcT25Wr17N1Vdfzfz582nTpg1r1qyp8jUrokaKmqCgIDp27EhmZmbZNqfTSWZm5kmVoafV2JOaP1jq7MSNtqfJcV5MQ1Mxrwe9yOMB89wjpURERCohOjqa4uJiNmzYgNPp5MUXX6Rz5860adOGPXv2eOw6fy441qxZc9omow4dOrB9+3aaNWtG69aty70aNGhQdtyVV17J2LFjWbVqFVdccQVvv/22x/KeiseKmqKiorJKDSA3N5ecnBx27doFQHp6OjNmzGDu3Lls3ryZkSNHUlxcTGpqqqcieJU8Vzi32J7g9dIkAIYFfAqzEuHITmODiYiIVzp06BDXXXcdb775Jt999x25ubm8//77PPfcc9x00020bt0au93Ov//9b37++WfmzZtHRkaGx66/cuVKnnvuObZt24bVauX999/n/vvvP+Wxt99+O02aNOGmm25ixYoV5ObmkpWVxejRo/n111/Jzc1l7NixrF69ml9++YVly5axffv2au9X47GOwuvXr6dnz55l79PT0wFISUlhzpw5DB48mAMHDjBu3Dj27dtHbGwsS5cuPanzsKfVSEfh07ATwMTSO1njjOaFwAwa7tkIGd3gpukQ3a/G84iI1HbePMtvaGgo8fHxvPTSS+zYsQO73U5kZCTDhw/n0UcfpU6dOkyZMoVnn32WsWPH0q1bNyZPnszf/vY3j1z/H//4B+vXr2fChAmEhYUxZcoUEhMTT3ls3bp1yc7O5pFHHuHmm2/mt99+o0WLFvTq1YuwsDCOHz/Oli1bmDt3LocOHeL8888nLS2Ne+65xyNZT8fkcrlc1XoFL1FYWEiDBg0oKCggLCzM4+c/W5+ZCA6yqvWb8Ot/Z1HudDfcMBECgj2eRUSkNjtx4gS5ublcdNFFhISEGB3HJ0RFRfHAAw+UG11Vk850zyrz+9snRz/5oj00gdQlcM1/H+WtfQ1mXg+HdhgbTERExE/4/Zz+RjY//VnUP5cB8fQwP8SUwFdovPdbfpt2DWPtd/Gx091h2psfjYqIiHgzv39SY8Top7PJcl5Jn5LJ/J/zMuqbjjM96N88HTCTYDQ6SkREat7OnTsNa3ryJL8varzVPs7jNts/+Xdpf5wuE7cHZLIoaBwc3G50NBEREZ/k90VNTU2+VxUOLLxYeit/s4/hgCuMy8274NXu8N17RkcTERHxOX5f1Hhj89Offe1sR5+SyaxyRIO9GD4cDh+NAtuxs39YREREgFpQ1PiKAzTiDvuj0GMsYIJv5sGM62B/9Sz4KSIi4m9U1HgRJ2boMQZSFkNoOBzYDK/1gG/eMjqaiIiI1/P7osab+9Sc1kXdYMTXcHFPKD0OH90LC0dASZHRyURExIdlZWVhMpk4evRohT8TFRXF1KlTqy2TJ/l9UeMLfWr+KGrMJ+7XxHVctGkYz9lvxeEywbfv8NPTcfQe+4rREUVEpJoMGTIEk8nEiBEjTtqXlpaGyWRiyJAhNR/MR/h9UePLXJh52dGfv9oeY5+rEa3Ne1gU9Disnw21Y3ULEZFaJzIyknfffZfjx4+XbTtx4gRvv/02LVu2NDCZ91NR4wPWui6nT8lkvnTEEGKyw8cPwIJhcKLQ6GgiIuJhHTp0IDIykg8//LBs24cffkjLli258sory7aVlJQwevRomjVrRkhICNdee+1JrRJLliyhTZs21KlTh549e7Jz586Trvf111/TtWtX6tSpQ2RkJKNHj6a4uLjavl91UlHjIw4TxlD7Q0y2/xVMFvhhAbzWHfZ+a3Q0ERHv53KBrbjmX1V8qj506FBmz55d9n7WrFmkpqaWO+bhhx9mwYIFzJ07l40bN9K6dWsSExM5fPgwAHl5edx8880kJyeTk5PDXXfdxZgxY8qdY8eOHfTu3ZuBAwfy3XffMX/+fL7++mtGjRpVpdxG09pPPsSFmVcdyYxNTYUPhsLhn+H1BEicBHF3gclkdEQREe9kPwaTImr+uo/ugaB6lf7YHXfcwdixY/nll18AWLlyJe+++y5ZWVkAFBcX88orrzBnzhySkpIAmDFjBsuXL2fmzJk89NBDvPLKK7Rq1YoXX3wRgEsvvZTvv/+eZ599tuw6kydP5vbbby9bIuGSSy5h2rRpdO/enVdeecXnVjn3+yc1vtZRuEJaxsOIFdAmCRw2WPIgvJ8Cx48anUxERDygadOm3HjjjcyZM4fZs2dz44030qRJk7L9O3bswG63c80115RtCwwMpFOnTmzevBmAzZs3Ex8fX+68Xbp0Kff+22+/Zc6cOYSGhpa9EhMTcTqd5ObmVuM3rB5+/6TGH0WN+eS//3UHwyzn8UjAOwRt+ohdP6xilH0037laAVrxW0SkTGBd91MTI65bRUOHDi1rBrJarZ5KVE5RURH33HMPo0ePPmmfL3ZKVlHj00zMdPRhvbMN0wP/TUvzAT4IeoJnSm9jlqO30eFERLyHyVSlZiAj9e7dG5vNhslkIjExsdy+Vq1aERQUxMqVK7nwwgsBsNvtrFu3rqwp6fLLL2fx4sXlPrdmzZpy7zt06MCmTZto3bp19X2RGuT3zU+1wbeu1txom8SnjjiCTA7GBc5jRuAUOHbY6GgiIlJFFouFzZs3s2nTJiwWS7l99erVY+TIkTz00EMsXbqUTZs2MXz4cI4dO8awYcMAGDFiBNu3b+ehhx5i69atvP3228yZM6fceR555BFWrVrFqFGjyMnJYfv27Xz00Uc+21FYRY2fKKQeI+0PMM6eQokrgOstG+DVbpC31uhoIiJSRWFhYYSFhZ1y3zPPPMPAgQO588476dChAz/99BOfffYZjRo1AtzNRwsWLGDRokXExMSQkZHBpEmTyp2jffv2fPXVV2zbto2uXbty5ZVXMm7cOCIiDOhU7QEml6t2zOJWWFhIgwYNKCgoOO3/IOfif/1cjNfWlIs1cBpR5nwwB0CvcdDlPjCrhhUR/3fixAlyc3O56KKLfG70Tm11pntWmd/f+i3nh350XURf29PQ9mZwlsLycfDOYCg+ZHQ0ERGRauP3RY1PLmjpAUXUhUGzoO9UsATD9mWQcS38ssroaCIiItVCzU8e4k3NT392mWkX1sB/0cq8l1KXmSmlg3jF0Q8XZg37FhG/o+Yn36PmJ6mwLa6WJNue5kPHtQSYnDwc+B5zA5/lPAqMjiYiIuIxKmpqiWOEkG4fyUP2uznuCqKb5Xs+DR4LudlGRxMREfEIFTW1ion3HT3oZ5vINmcLmpmOwhs3QdYz4PT9tbFERP6olvSu8AueuleaUbgW2u66gH62iTwZMIdbA76CrMmszFzMA/Y0DtDwpOPV70ZEfElgYCAAx44do06dOgankYqw2WwAJ00yWFkqamqpEwTzcOk9rHZGMzFwFtdYfmSJeQwP2NNY6WxndDwRkSqzWCw0bNiQ/fv3A1C3bl1MJpPBqeR0nE4nBw4coG7dugQEnFtZ4lNFTVRUFGFhYZjNZho1asSXX35pdCSft9DZle9sFzM9cBqXm/OYF/gMVsdNTC0diINzq5hFRIzSvHlzgLLCRryb2WymZcuW51x8+lRRA7Bq1SpCQ0ONjuFXdrha0N/2FOMC5nF7QCb3BSyik3kLo22jyKex0fFERCrNZDJx/vnn06xZM+x2u9Fx5CyCgoIwe2DWe58raqR6lBDEP0uHscZ5OZMCZxJv3sKS4LGk2+8F1KdGRHyTxWI5534a4js8NvopOzub5ORkIiIiMJlMLFq06KRjrFYrUVFRhISEEB8fz9q1lVts0WQy0b17d+Li4njrrbc8lFz+6D/Oq0m2TeQHZxTnmX5jbtCzsHw8OPQvHRER8W4ee1JTXFxMTEwMQ4cO5eabbz5p//z580lPTycjI4P4+HimTp1KYmIiW7dupVmzZgDExsZSWlp60meXLVtGREQEX3/9NS1atGDv3r0kJCTQrl072rdvf8o8JSUllJSUlL0vLCz00Df1fztd5zPQ9gSPBrxFSsByWDkVdq2GgTOhYaTR8URERE6pWpZJMJlMLFy4kP79+5dti4+PJy4ujunTpwPu3s6RkZHcd999jBkzptLXeOihh2jbti1Dhgw55f4nnniCCRMmnLS9Ni6TcC56m9eSUX82lBRAnUbQ/xW4NMnoWCIiUkt43TIJNpuNDRs2kJCQ8L8Lm80kJCSwevXqCp2juLiY3377DYCioiK++OIL2rZte9rjx44dS0FBQdkrLy/v3L5ELbXU2Qnu+QoiOsDxI/DOX+Czf0KpzehoIiIi5dRIR+GDBw/icDgIDw8vtz08PJwtW7ZU6Bz5+fkMGDAAAIfDwfDhw8+48nZwcDDBwcFYrVasVisOh2bMraqo5zYRyAOMCXiHYQGfwurp5Kz8lFH2+/jV5W461AR9IiJiNJ8Z/XTxxRfz7bffVvpzaWlppKWllT2+kqqxE8BTpXeyxnk5LwRmEGvewZKgR3nIfg+fOU9fXIqIiNSUGml+atKkCRaLhfz8/HLb8/PzyyZIqi5Wq5Xo6OgzPtWRilvuvIo+JZPZ6GxNmOkYrwa9xPiAuVBacvYPi4iIVKMaKWqCgoLo2LEjmZmZZducTieZmZl06dKlWq+dlpbGpk2bWLduXbVepzbZTVNutY0jo7QvAKkBn8HM6+HwzwYnExGR2sxjRU1RURE5OTnk5OQAkJubS05ODrt27QIgPT2dGTNmMHfuXDZv3szIkSMpLi4mNTXVUxGkBpUSwDOlt5Fqe4jDrlDY+y1kdIMfPjQ6moiI1FIe61Ozfv16evbsWfY+PT0dgJSUFObMmcPgwYM5cOAA48aNY9++fcTGxrJ06dKTOg97mjoKV68vnVfSp2Qya9q87Z7L5oNU2LkCEidBoFbHFRGRmlMt89R4o8qMc68Kf52npqJ2Pp0IWZNgxRTABeFXwC1zoMklRkcTEREf5nXz1EgtYAmAXuPgjgVQtwnk/wCvdofv3jM6mYiI1BJ+X9Ro9FMNa90LRnwNUV3BXgwfDoePRoHtmNHJRETEz/l9UaPRTwYIOx/+9hF0HwOY4Jt5MOM62F+xiRZFRESqwu+LGjGI2QI9x7qLm9BwOLAZZvSEb7S6uoiIVA+/L2rU/GSwi7u7m6Mu7gH2Y/DRvbBwBJQUGZ1MRET8jN8XNWp+8gKhzeCOhXDdY2Ayw7fvuJ/a5P9odDIREfEjPrP2k3i3igxp3/nMQ9DyalgwDA5uc/ezSXoWOqSAyVQDKUVExJ9pnhoPqe3z1FRGYwqZEvgKPSzuBUo/clzNo/ZhFOOerE8rfouIyO80T80fqE+N9zlMGKn2h5hs/yulLjM3WVbxcdCjtDXtNDqaiIj4ML8vatSnxju5MPOqI5lbbePY7TqPi8z5fBg0njssy6F2PDwUEREP8/uiRrzbRlcbbiyZxHJHB4JNdiYGzob3U+BEgdHRRETEx6ioEcMdpT7D7f/gKfsd2FwW2PQRvNoNdm80OpqIiPgQFTXiJUzMdPThFtt4aNgSjuyEmTfAmlfUHCUiIhXi90WNOgr7lm9dreGeFXB5MjjtsHQMzL8Djh8xOpqIiHg5vy9q1FHYB9VpCLfOg6TnwRIEWz6GjG7w63qjk4mIiBfT5Hvidf43508LrjCNY3rgv4kq2IV9xg08W/oXZjqSyH0m2dCMIiLiffz+SY34th9cF9PX9jQfOzoTaHLwWOBbvB74Ihw7bHQ0ERHxMipqxOsVUZdR9vt41D6MElcgvSzfQMa1sGuN0dFERMSLqKgRH2HibUcv+tueZIfzfCjcDbP7wIoXwek0OpyIiHgB9akRn7LZdSH9bBOZGDiLAZaVkPkkXy1bRLp9JIdoUO5YrSElIlK76EmN+Jxi6vB3+708ZL+b464gulu+Y0nwWOJNm42OJiIiBvL7okbz1PgrE+87enCT7Sm2O1sQbjrK20ETGW35EDNqjhIRqY38vqjRPDX+bZsrkn62p3ivtDsWk4v0wA+YFziZphw1OpqIiNQwvy9qxP8dJ4SHS+/h77aRHHMFc43lR5YEj4EdXxodTUREapCKGvEbC51dSbZNZLMzkqamQpg3AL6YCI5So6OJiEgNUFEjfmWHqwX9bU/xdul1gAuyn4c3+kHhHqOjiYhINVNRI36nhCAeLb0LBs6EoFD4ZaV7sr7tnxsdTUREqpGKGvFf7QbBPdnQvB0cOwRvDYTPnwCH3ehkIiJSDUwul8tldIiKys3NZejQoeTn52OxWFizZg316tWr0GcLCwtp0KABBQUFhIWFeTzb/xZhFG8TjI1HA94iJWA5AOudbbjPdh97OU8T9ImIeLnK/P72qSc1Q4YM4cknn2TTpk189dVXBAcHGx1JfEAJQYwvTWWk7X4KXXW4yryNJcFjuc680ehoIiLiQT5T1Pz4448EBgbStWtXABo3bkxAgFZ5kIr71BnPjbZJfOu8mEamImYFvQCf/RNKbUZHExERD/BYUZOdnU1ycjIRERGYTCYWLVp00jFWq5WoqChCQkKIj49n7dq1FT7/9u3bCQ0NJTk5mQ4dOjBp0iRPRZdaJM8Vzi228cwsTXJvWD0dZifBkV+MDSYiIufMY0VNcXExMTExWK3WU+6fP38+6enpjB8/no0bNxITE0NiYiL79+8vOyY2NpYrrrjipNeePXsoLS1lxYoVvPzyy6xevZrly5ezfPny0+YpKSmhsLCw3EsEwEYgT5XeyXBbOoQ0gN3r4dWusPk/RkcTEZFzUC0dhU0mEwsXLqR///5l2+Lj44mLi2P69OkAOJ1OIiMjue+++xgzZsxZz7l69WqeeOIJPvvsMwCef/55AB566KFTHv/EE08wYcKEk7aro7D80c4x7eCDofDrf5fR6HQP3PAUBKi/loiIN/C6jsI2m40NGzaQkJDwvwubzSQkJLB69eoKnSMuLo79+/dz5MgRnE4n2dnZXH755ac9fuzYsRQUFJS98vLyzvl7iB9q2BJSP4Wr73O/X/sqzLwBDv9sbC4REam0GilqDh48iMPhIDw8vNz28PBw9u3bV6FzBAQEMGnSJLp160b79u255JJL6Nu372mPDw4OJiwsjHnz5tG5c2d69ep1Tt9B/JglEG6YCLe9B3Uaw94cyOgGP3xodDIREakEnxo+lJSURFJSUqU+k5aWRlpaWtnjK5HTapMII752N0flrYEPUmHnCkicDIEhRqcTEZGzqJEnNU2aNMFisZCfn19ue35+Ps2bN6/Wa1utVqKjo4mLi6vW64ifaNAChnwC16a736+fBa8nwMGfjM0lIiJnVSNFTVBQEB07diQzM7Nsm9PpJDMzky5dulTrtdPS0ti0aRPr1q2r1uuIH7EEQMJ4uGMB1G0C+d/Dq93gu/eMTiYiImfgsaKmqKiInJwccnJyAPeSBjk5OezatQuA9PR0ZsyYwdy5c9m8eTMjR46kuLiY1NRUT0UQ8azWCe7mqKiuYC+GD4fDR6PAdszoZCIicgoeG9KdlZVFz549T9qekpLCnDlzAJg+fTrPP/88+/btIzY2lmnTphEfH++Jy5+W1WrFarXicDjYtm2bhnRLpZlx8nPi9/DVc4ALmkXDLXOg6aVGRxMR8XuVGdLtUwtangstaCnnYuczN8LPWbBgOBTvh8C6cOOLEHub0dFERPyaippTUFEjntCEAl4KtNLV8gMACxxdecyeynFCtOK3iEg18LrJ94yk0U/iSQdpQIp9DC/Yb8HhMjHQsoL/BD3GpaZdRkcTEan19KTGQ/SkpvaJN23mX0HTaW46wglXIONLhzDf0QMwnfYzepojIlI5elIjUgP+z3U5fUomk+WIIcRk59nAGUwNtFKP40ZHExGplfy+qFHzk1Snw4SRan+IZ+x/odRlpr9lFf8J+ifRpp1GRxMRqXX8vqjR5HtS3VyYyXD041bbOHa7zuNi8z4WBo3nDstyoFa07oqIeAW/L2pEaspGVxtuLJnEckcHgk12JgbOZnrgNOqjyfpERGqCihoRDzpKfYbb/8FT9tuxuyz0tfwfHwc9SjvTz0ZHExHxe35f1KhPjdQ8EzMdN3KLbTy/uppwoXk/C4LGM8SyFGrHYEMREUNoSLeHaEi3nEoYRTwXOIPelv/26bqsL9w0Heo0MjaYiIiP0JBuES9RSCgj7A8w3p5CiSsAtnwMGd3g1/VGRxMR8TsqakSqnYm5jkQG2p6ARhdBwS6YlQir/q3mKBERDwowOoBIbfGD62La7f0nkwNfpy9rYNljfP7pAh60j+Ao9cuO06zDIiJV4/dPatRRWLzJb9RllP0+/mkfSokrkATLN3wS/CgdTVuNjiYi4vP8vqjR5HvifUy85UhggG0CPzub08J0iPlBTzHSshgTTqPDiYj4LL8vakS81SZXFMm2p1nkuJoAk5NHAt9lduDzUHzQ6GgiIj5JfWpEDFRMHR6wp7HK2ZYnA+bQw/It+567itG2Uax1XX7az6nfjYjIyfSkRsRwJt5z9KSfbSLbnS1objrCO0ETGWVZiFnNUSIiFaaiRsRLbHNF0s/2FB84umExuXgw8H3eCJxMU44aHU1ExCeoqBHxIscJ4UH7CP5hG8ExVzDXWn5kSfBYrjb/YHQ0ERGvp6JGxAstcHYj2TaRLc5ImpoKeDNwMn8P+EDNUSIiZ+D3az9ZrVasVisOh4Nt27Zp7SfxKSGUMD7gDf4a8CUAa5yXM9o2iv1UbO0odSgWEV+ntZ/+QPPUiC87QTBjS4cz2pZGkSuEzubNfBo8hm7mb42OJiLidfy+qBHxB4ud15Bse5pNzgs5z/QbbwQ9y8MB72LBYXQ0ERGvoaJGxEfkus5ngG0C80oTALg3YDHvBj3F+RwyOJmIiHdQUSPiQ0oI4vHSoaTZRlPoqkOceRtLgsdynXmj0dFERAynokbEB33i7Exf2yS+c15EI1MRs4Je4NGAtwik1OhoIiKGUVEj4qN2ucIZZHuC2aWJANwd8AnvBT3JBaYDBicTETGGzxQ1W7duJTY2tuxVp04dFi1aZHQsEUPZCGRCaQp32/5OgasuV5p/4pOgsdxg1mg/Eal9fGZBy0svvZScnBwAioqKiIqK4vrrrzc2lIiXWOaMY5Mtin8H/psrzT/xWtBL7ic4pQkQEGx0PBGRGuEzT2r+aPHixfTq1Yt69eoZHUXEa/zqasqttnG8VuqecC814DOYeQMc/tngZCIiNcNjRU12djbJyclERERgMplO2TRktVqJiooiJCSE+Ph41q5dW6VrvffeewwePPgcE4v4HzsBTCq9naG2BzniCoW9OfBqd/hxodHRRESqnceKmuLiYmJiYrBarafcP3/+fNLT0xk/fjwbN24kJiaGxMRE9u/fX3ZMbGwsV1xxxUmvPXv2lB1TWFjIqlWr6NOnzxnzlJSUUFhYWO4lUlt84exAn5LJENkZSgrh/SHwcTrYTxgdTUSk2lTL2k8mk4mFCxfSv3//sm3x8fHExcUxffp0AJxOJ5GRkdx3332MGTOmwueeN28en332GW+++eYZj3viiSeYMGHCSdu19pPUJjufvgG+nARfT3FvaN4ObpkL57UyNpiISAV53dpPNpuNDRs2kJCQ8L8Lm80kJCSwevXqSp2rok1PY8eOpaCgoOyVl5dX6dwiPs8SCAnj4Y4FULcJ7PseXu0G339gdDIREY+rkaLm4MGDOBwOwsPDy20PDw9n3759FT5PQUEBa9euJTEx8azHBgcHExYWxrx58+jcuTO9evWqdG4Rv9E6AUZ8DRdeC7YiWDAMFt8H9uNGJxMR8RifGv3UoEED8vPzCQoKqvBntEq3yH+FnQ9/+wi6PwKYYOMbMOM6OLDV6GQiIh5RI0VNkyZNsFgs5Ofnl9uen59P8+bNq/XaVquV6Oho4uLiqvU6Ij7BEgA9H4W/LYJ6zWD/JnitB+S8Y3QyEZFzViNFTVBQEB07diQzM7Nsm9PpJDMzky5dulTrtfWkRuQULu4BI1fCRd3BfgwWjYBF94Kt2OhkIiJV5rGipqioiJycnLJZf3Nzc8nJyWHXrl0ApKenM2PGDObOncvmzZsZOXIkxcXFpKameiqCiFRGaDO4cyH0fAxMZsh5C17rCfmbjE4mIlIlHlsmYf369fTs2bPsfXp6OgApKSnMmTOHwYMHc+DAAcaNG8e+ffuIjY1l6dKlJ3Ue9jSr1YrVasXhcFTrdUR8ktkC3R+CC7vAgrvg4FZ3P5s+z8GVd4LJZHRCEZEKq5Z5arxRZca5V4XmqRFvtPOZGyt+cPFB+PBu2PHfZuJ2t0LfKRBcv3rCiYhUgNfNUyMiPqBeE7j9A+g1HkwW+P49dyfifd8bnUxEpEL8vqjR6CeRSjCboWs6pC6BsBZw6CeY0QvWzYTa8VBXRHyY3xc1Gv0kUgUtO7sn67skERwl8Ek6fJAKJwqMTiYiclp+X9SISBXVbQx/fRdumAjmAPdK3692hz3fGJ1MROSUPDb6yVtp9JPUZhXpwH7GzsRmM1x9n3u17w+GwpFcmHmDu9DpdLdGR4mIV/H7JzVqfhLxgMg4GJENl/UFhw0+fRjm3wHHjxidTESkjN8XNSLiIXUaweA3ofezYA6ELR+7V/z+dYPRyUREABU1IlIZJhN0HgHDlkGjKDi6C2bdAKuma3SUiBjO7yff+2Ofmm3btmnyPZEqOGW/mxMFsPg+2PSR+32bJOj/sruDsYiIh2jyvT9QnxqRahLSAG6ZCze+CJZg2PYpZHSFXf9ndDIRqaX8vqgRkWpkMkHcXXDX59C4FRT+CrOT4OuXwOk0Op2I1DIqakTk3J3fHu75CtrdAi4HfP4EvH2Lez0pEZEaoqJGRDwjuD7cPAOSp0FACPz0OWRcCztXGp1MRGoJvy9qtPaTSA0ymaBjCgz/Apq0gd/2wty+8NXz4NQEmCJSvfy+qFFHYREDhLeFu7Mg5jZwOeHLiTBvABTtNzqZiPgxvy9qRMQgQfVgwCvQ/xUIrAu5X8Er18DPWUYnExE/paJGRKpX7G0w/EtoejkU74c3+sOXk9QcJSIep6JGRKpfs8vc/Ww6/A1wwVfPwhs3QeFeo5OJiB/x+1W6RcRLBNUlalVv+pnrMylwJqE7V3DwxU6k20eS7Ywpd+gZVw4XETkNFTUiclYVWQakooXIYuc1fG+7GGvgNKLNv/BG0LNYS/sxpfQWHFjONaqI1GJqfhKRGpfrOp8Btgm8WdoLgLSAxbwTNJHmHDI4mYj4Mr8vajRPjYh3KiGIx0qHkWYbzW+uOnQyb2VJ8Fh6mr8xOpqI+Ci/X6X7d5VZ5bMqtEq3SNW1NOVjDfwX7cw73Ruuvg96jQdLoKG5RMR4WqVbRHzKLlc4A20TmF2a6N6w6t/uhTGP7jI2mIj4FBU1IuIVbAQyoTSFe2x/h5AG8Os699pRW/QUVEQqRkWNiHiVz5xxcM8KaNERThTAu7fBp2Og1GZ0NBHxcipqRMT7NLoQUpdCl1Hu9//3Csy6AQ7nGptLRLyaihoR8U4BQZD4NPz1XajTCPZ8A692gx8XGZ1MRLyUTxU1L730Em3btiU6OprRo0dTSwZuidRulya5m6Mi46GkEN5PgU/+AfYTRicTES/jMzMKHzhwgOnTp/Pjjz8SGBhIt27dWLNmDV26dDE6moh42KmmSAggjfSAcO4NWAzrXoe8/4Nb5sJ5rQxIKCLeyKee1JSWlnLixAnsdjt2u51mzZoZHUlEakgpATxX+hdSbI9wyFUf9n3vbo76/gOjo4mIl/BYUZOdnU1ycjIRERGYTCYWLVp00jFWq5WoqChCQkKIj49n7dq1FT5/06ZNefDBB2nZsiUREREkJCTQqpX+hSZS23zljKFPyWS48BqwFcGCYbB4NNiPGx1NRAzmsean4uJiYmJiGDp0KDfffPNJ++fPn096ejoZGRnEx8czdepUEhMT2bp1a9kTl9jYWEpLS0/67LJly6hTpw4ff/wxO3fupE6dOiQlJZGdnU23bt1OmaekpISSkpKy94WFhR76piJitHwa02rrCO4PaMooy0eYN85ly/ovSLOPZoerRdlxWu1bpHaplmUSTCYTCxcupH///mXb4uPjiYuLY/r06QA4nU4iIyO57777GDNmzFnP+f7775OVlYXVagXg+eefx+Vy8fDDD5/y+CeeeIIJEyactF3LJIj4l2vM3zM18GWamgo45grmMXsqHzrd/9hRUSPi+7xumQSbzcaGDRtISEj434XNZhISEli9enWFzhEZGcmqVas4ceIEDoeDrKwsLr300tMeP3bsWAoKCspeeXl55/w9RMT7rHS2o0/JZFY62lLXVMKUoAyeD8igDhodJVLb1EhRc/DgQRwOB+Hh4eW2h4eHs2/fvgqdo3PnzvTp04crr7yS9u3b06pVK/r163fa44ODgwkLC2PevHl07tyZXr16ndN3EBHvdYCG3GkfyxT7IBwuE7cEZLM46HHI32R0NBGpQT41+unpp59m8+bN/Pjjj0ybNg2TyXTWz6SlpbFp0ybWrVtXAwlFxChOzExz3Mzt9n+S72rIJebdMOM62PgGaE4rkVqhRoqaJk2aYLFYyM/PL7c9Pz+f5s2bV+u1rVYr0dHRxMXFVet1RMQ7rHFG06dkMtmOdlB6HBbfBx/eDSW/GR1NRKpZjRQ1QUFBdOzYkczMzLJtTqeTzMzMap88T09qRGqfQzQgxf4I9BoHJgt8/x681sM9t42I+C2PDekuKirip59+Knufm5tLTk4OjRs3pmXLlqSnp5OSksJVV11Fp06dmDp1KsXFxaSmpnoqgohIGRdm6PoPaHm1ey6bQz/BjF7QezJcNRT+23xdkZGLGkUl4hs8VtSsX7+enj17lr1PT08HICUlhTlz5jB48GAOHDjAuHHj2LdvH7GxsSxduvSkzsOeZrVasVqtOByOar2OiHipC7u4145aNBK2fwafpMPOFZA8DUI8P72DiBinWuap8UaVGedeFZqnRsT7lHvC4nTC6umQOQGcpdDoIrhlDlHTdlfuPCJSo7xunhoREcOZzXDNaEhdCg1awpFcmHk9KZbPgFrxbzsRv+f3RY1GP4lIOZFxMCIbLusLDhsTAufySuBUwig2OpmInCO/L2o0+klETlKnEQx+E3o/i81lIcmyjk+CHiXG9NPZPysiXsvvixoRkVMymaDzCAbaJvCLsxmR5gN8EDSBYZZPUHOUiG/y+6JGzU8icibfuy6mr20Snzg6EWhy8HjgW8wIfJEGFBkdTUQqye+LGjU/icjZ/EZd0uz385g9lRJXINdbNrIkeCwdTNuMjiYileD3RY2ISMWYeNNxPQNsE/jZ2ZwWpkO8F/QkIyyL3cPBRcTrqagREfmDTa4okm1P85HjagJMTsYEvgtv3wrFB42OJiJn4fdFjfrUiEhlFVOH++1pjLHfxQlXIPy0HDKuhZ0rjY4mImfg90WN+tSISNWYeNdxHTfZnoImbeC3vTC3L3z1PDi17IqIN/L7okZE5FxsdbWE4V9CzF/B5YQvJ8KbN0PRfqOjicifqKgRETmb4FAYkAE3vQyBdeHnLHdz1M9fGZ1MRP5ARY2ISEVdebv7qU3Ty6EoH964Cb6crOYoES8RYHSA6ma1WrFarTgc+ktHRKomaswn5d6H8DBPBMzlLwFZ8NUzrP5iMaPtaax75g5jAooIUAue1KijsIh42gmCGVN6N/fb7qXYFUwXyyY+DR4LP2UaHU2kVvP7okZEpLp85LyWvrZJbHJeSBNTIbw5EDKfBEep0dFEaiUVNSIi5yDXdT4DbBOYV5oAuGDFi+6h3wW7jY4mUuuoqBEROUclBPF46VAYNBuC6sOu1e7RUduWGR1NpFZRUSMi4ilX3AwjsuH8GDh+GN6+BZY9Dg670clEagUVNSIintT4Yhi2HDrd436/ahrMToKju4zNJVILqKgREfG0gGDo8xzcOg+CG8Cv6yCjK2z55OyfFZEq8/uiRgtaiohhovu5m6NadIQTR+Hd2+DTMVBqMzqZiF/y+6JG89SIiKEaRUHqUugyyv3+/16BWTfA4VxDY4n4I7+fUVhExHABQZD4NERdCwtHwJ5v4NVu0O/f0Lb/SYf/eQbjU9n5zI3VEFTEt/n9kxoREa9xaRKM+Boi46GkEN5PgU/+AfYTRicT8QsqakREalLDSBjyCVzzgPv9utdh5vVwaIehsUT8gYoaEZGaZgmE6yfA7R9A3fNg33fu5qjvPzA6mYhPU58aERGjXHK9uznqg2GwaxUsGAY7VxBMD0oIOuNH1e9G5GQ+9aTmhRdeoG3btlxxxRW8+eabRscRETl3YRGQ8h/o9hBggg1zWBT0OK1MWjtKpLJ8pqj5/vvvefvtt9mwYQPr1q1j+vTpHD161OhYIiLnzhIA1z0Gd34I9ZpyuTmPxUGPMcC8wuhkIj7FZ4qazZs306VLF0JCQqhTpw4xMTEsXbrU6FgiIp7T6joY8TWrHNHUM5XwUtArPBfwKnXQ6CiRivBYUZOdnU1ycjIRERGYTCYWLVp00jFWq5WoqChCQkKIj49n7dq1FT7/FVdcQVZWFkePHuXIkSNkZWWxe7cez4qIn6nfnDvsjzLFPgiHy8StAV/xUdDjXGL61ehkIl7PYx2Fi4uLiYmJYejQodx8880n7Z8/fz7p6elkZGQQHx/P1KlTSUxMZOvWrTRr1gyA2NhYSktLT/rssmXLiI6OZvTo0Vx33XU0aNCAzp07Y7FYTpunpKSEkpKSsveFhYUe+JYiItXPiZlpjptZ67qMaYHTaWPezeKgxxhXOoT3Hd0Bk9ERRbySx57UJCUlMXHiRAYMGHDK/VOmTGH48OGkpqYSHR1NRkYGdevWZdasWWXH5OTk8MMPP5z0ioiIAOCee+5h48aNfPnllwQGBnLJJZecNs/kyZNp0KBB2SsyMtJTX1VEpEascUaTVDKZbEc76phsPB/4Gi8GvkJdNUeJnFKN9Kmx2Wxs2LCBhISE/13YbCYhIYHVq1dX+Dz79+8HYOvWraxdu5bExMTTHjt27FgKCgrKXnl5eVX/AiIiBjlEA1Lsj/Cc/VYcLhMDLV/zn6B/cplpl9HRRLxOjcxTc/DgQRwOB+Hh4eW2h4eHs2XLlgqf56abbqKgoIB69eoxe/ZsAgJOHz84OJjg4GCsVitWqxWHw1Hl/CIiRnJh5mVHf9Y5L2Na0HRamffyUdDjPFH6N95xXIeao0TcfGb0E8Dq1avLVtzu2LFjhT6jVbpFxF+sc11Gn5JJfOGIJdhkZ3LgTKYFTieUY0ZHE/EKNVLUNGnSBIvFQn5+frnt+fn5NG/evFqvbbVaiY6OJi4urlqvIyJSE44QxjD7g0yy/xW7y0I/y2r+E/RP2ppyjY4mYrgaKWqCgoLo2LEjmZmZZducTieZmZl06dKlWq+tJzUi4m9cmHnNkcxg2+P86mrCReZ8Pgwaz98snwEuo+OJGMZjRU1RURE5OTnk5OQAkJubS05ODrt2uTuzpaenM2PGDObOncvmzZsZOXIkxcXFpKameiqCiEitstHVhhtLJrHc0ZFgUylPBs7l5cB/EUax0dFEDOGxjsLr16+nZ8+eZe/T09MBSElJYc6cOQwePJgDBw4wbtw49u3bR2xsLEuXLj2p87CnqaOwiPizAkIZbk9nqHMpYwLepo9lLVeYchllH210NJEaZ3K5XLXiWWVhYSENGjSgoKCAsLAwj5+/IivmikjNqsgq1Z782fXUqthVzdTetIPpgdNoaT6AzWUhqPdE6DwSTBodJb6rMr+/fWr0k4iInN53rlb0tU1iiaMTQSYHfDYW3r0Njh02OppIjfD7okajn0SkNimkHvfa7+dx+xCwBMHWJfBqN8ir+Fp7Ir6qRibfM1JaWhppaWllj69ERPyfiXmOG3jqrqHw/hA4/DPMToJe46DLfWA27t+zFWla81QzntQ+fv+kRkSk1jo/Bu7+Cq4YCM5SWD4O3hkMxYeMTiZSLfy+qFHzk4jUaiFhMHAmJP8LAkJg+zLIuBZ+WWV0MhGP8/uiRpPviUitZzJBxyFwVyacdwn8tgfm9IXsF8DpNDqdiMf4fVEjIiL/1fwKuDsL2v8FXA744il482Yo2m90MhGPUFEjIlKbBIfCgAy4yQoBdeDnL93NUbnZRicTOWd+X9SoT42IyJ+YTHDlHe6nNk0vg6J8eOMmyHoGnJp9XXyX3xc16lMjInIazS6D4V+6CxyXE7Imu4ub3/YZnUykSvy+qBERkTMIqutuihrwGgTWg50r3M1RO74wOplIpamoERERiBnsbo4KvwKKD8C8myHzKXCUGp1MpMJU1IiIiFvTNnDX59AxFXDBihdgbjIU7DY6mUiF+H1Ro47CIiKVEFgHkqe6J+wLqg+7Vrmbo7YtMzqZyFn5fVGjjsIiIlXQbhDc8xU0bw/HD8Pbt8Cyx8FhNzqZyGn5fVEjIiJVdF4rd3NUp7vd71dNg9l94GiesblETkNFjYiInF5AMPR5Hm59A4IbwK9r3c1RW5YYnUzkJCpqRETk7KJvcjdHRXSAE0fh3b/C0keh1GZ0MpEyKmpERKRiGl8EQz+Dzmnu92usMCsRjuw0NJbI71TUiIhIxQUEQe9J8Jd3IKQh7NkIGd1g02Kjk4moqBERkSq4rA+MWAEXdIKSAnjvTljyEJSWGJ1MajG/L2o0T42ISDVp2BJSl8A197vfr30NZl4Ph3YYm0tqLb8vajRPjYhINbIEwvVPwu0fQJ3GsPdbeLU7/LDA6GRSC/l9USMiIjXgkuthxNfQ8mqw/QYfDIX/PAD240Ynk1pERY2IiHhGgxaQ8h/o+iBggg2z4fUEOLjd6GRSS6ioERERz7EEQK/H4c4PoV5TyP/B3Rz17Xyjk0ktoKJGREQ8r9V17uaoqK5gL4aFd8NHaYSg0VFSfVTUiIhI9ajfHP72EfQYC5jgmzf5KOhxWpt+NTqZ+CmvLGoGDBhAo0aNGDRo0En7Pv74Yy699FIuueQSXn/9dQPSiYhIhZkt0GMMpCyG0HAuNf/Kf4IeY5DlK6OTiR/yyqLm/vvv54033jhpe2lpKenp6XzxxRd88803PP/88xw6dMiAhCIiUikXdYMRK8l2tKOOycYLga/yYuAr1OWE0cnEj3hlUdOjRw/q169/0va1a9fStm1bWrRoQWhoKElJSSxbtsyAhCIiUmmhTUmxP8Jz9ltxuEwMtKxgcdBjXGraZXQy8ROVLmqys7NJTk4mIiICk8nEokWLTjrGarUSFRVFSEgI8fHxrF271hNZ2bNnDy1atCh736JFC3bv3u2Rc4uISPVzYeZlR3/+Ynucva7GtDbv4aOgx/mL5QvAZXQ88XGVLmqKi4uJiYnBarWecv/8+fNJT09n/PjxbNy4kZiYGBITE9m/f3/ZMbGxsVxxxRUnvfbs2VP1b/InJSUlFBYWlnuJiIh3WOe6jD4lk/jSEUOIyc4zga8zLXA6oRwzOpr4sIDKfiApKYmkpKTT7p8yZQrDhw8nNTUVgIyMDD755BNmzZrFmDFjAMjJyalS2IiIiHJPZnbv3k2nTp1OeezkyZOZMGFCla4jIiLV7whhDLU/xHDnJzwcMJ9+ltW0M/0Me9vA+TFGxxMf5NE+NTabjQ0bNpCQkPC/C5jNJCQksHr16nM+f6dOnfjhhx/YvXs3RUVFfPrppyQmJp7y2LFjx1JQUFD2ysvLO+fri4iIZ7kw85ojmVtt4/jV1YSLzPnuWYjXzgCXmqOkcjxa1Bw8eBCHw0F4eHi57eHh4ezbt6/C50lISOCWW25hyZIlXHDBBWUFUUBAAC+++CI9e/YkNjaWf/zjH5x33nmnPEdwcDBhYWHMmzePzp0706tXr6p/MRERqVYbXW24sWQSyxwdwWGDJQ/C+ylwosDoaOJDKt38VBM+//zz0+7r168f/fr1q/C50tLSSEtLo7CwkAYNGnginoiIVIMCQrnbns7OPrtg+TjY9BHsyYFbZkOLjkbHEx/g0Sc1TZo0wWKxkJ+fX257fn4+zZs39+SlKsxqtRIdHU1cXJwh1xcRkcowQZd7Ydhn0LAlHP0FZibCmlfUHCVn5dGiJigoiI4dO5KZmVm2zel0kpmZSZcuXTx5qQpLS0tj06ZNrFu3zpDri4hIFbToCPesgMuTwWmHpWPg3dvh2GGjk4kXq3RRU1RURE5OTtkIptzcXHJycti1yz15Unp6OjNmzGDu3Lls3ryZkSNHUlxcXDYaSkREpELqNIRb50GfF8ASBFs/gVe7QZ7+kSqnVuk+NevXr6dnz55l79PT0wFISUlhzpw5DB48mAMHDjBu3Dj27dtHbGwsS5cuPanzcE2xWq1YrVYcDoch1xcRkXNgMkGn4XBBHLw/BI7kwuze0GscdLkPzF45Mb4YpNJFTY8ePXCdpV1z1KhRjBo1qsqhPEkdhUVE/EBELNyTDf+5H3780N2ReOfX0D8D6p16FKzUPipxRUTEN4SEwaBZ0PclsATD9mWQcS38ssroZOIl/L6o0egnERE/YjLBVUNheCac1xp+2wNz+kL2C+B0Gp1ODOb3RY1GP4mI+KHm7eDur6D9YHA54Iun4K2BUHTA6GRiIL8vakRExE8Fh8KAV6HfdAioAzu+gIxrIDfb6GRiEL8vatT8JCLix0wm6HAn3P0lNL0MivLhjZsg6xlwatRrbeP3RY2an0REaoFml8PwLyD2DnA5IWuyu7j5reLrDorv8/uiRkREaomgetDf6m6SCqwHO1e4R0ft+MLoZFJDVNSIiIh/ifkL3J0FzdpC8QGYdzNkPgWOUqOTSTXz+6JGfWpERGqhpm3cw747DgFcsOIFmJsMhXuMTibVyO+LGvWpERGppQLrQPK/YOBMCAqFXavczVHblxudTKqJ3xc1IiJSy7Ub5F5ioXl7OHYI3hrkXmbBYTc6mXiYihoREfF/57WCYcshbrj7/cp/wZwb4WiesbnEo1TUiIhI7RAYAje+ALe+AcENIO//3M1RW5YYnUw8xO+LGnUUFhGRcqJvgnu+gogOcOIovPtX+OyfUGozOpmcI78vatRRWERETtL4Ihj6GXS+1/1+9XSY3RuO7DQ0lpwbvy9qRERETikgCHpPhr+8AyENYfcGyOgGmxYbnUyqSEWNiIjUbpf1gREr4II4KCmA9+6EJQ9BaYnRyaSSVNSIiIg0bAmpn8I197vfr30NZl4Ph3YYm0sqRUWNiIgIgCUQrn8Sbnsf6jSGvd/Cq93hhw+NTiYVpKJGRETkj9rcACO+hpZdwPYbfJAKH/8d7MeNTiZnoaJGRETkzxq0gJSPoes/ABOsnwWvJ8DB7UYnkzPw+6JG89SIiEiVWAKg1zi4YwHUbQL5P7ibo757z+hkchp+X9RonhoRETknrXvByJUQ1RXsxfDhcPhoFNiOGZ1M/sTvixoREZFzVr85/O0j6D4GMME382DGdbB/i9HJ5A9U1IiIiFSE2QI9x7qLm9BwOLAZZvSEb94yOpn8l4oaERGRyri4u3t01MU9wH4MProXFo6AkiKjk9V6KmpEREQqK7QZ3LEQrnsMTGb49h33U5v8H41OVqupqBEREakKsxm6PeQe+l3/fDi4zd3PZsNccLmMTlcreWVRM2DAABo1asSgQYMqtU9ERKTGRV3jbo5qfT2UnoD/jIYFd0HJb0Ynq3W8sqi5//77eeONNyq9T0RExBD1msBt70HCBDBZ4IcP4NVu7qUWpMZ4ZVHTo0cP6tevX+l9IiIihjGb4doH3Atjhl0Ah3+G16+HtTPUHFVDKl3UZGdnk5ycTEREBCaTiUWLFp10jNVqJSoqipCQEOLj41m7dq0nsoqIiHi/lvEwYgW0SQJHCSx5EN5PgRMFRifze5UuaoqLi4mJicFqtZ5y//z580lPT2f8+PFs3LiRmJgYEhMT2b9/f9kxsbGxXHHFFSe99uzZU/Vv8iclJSUUFhaWe4mIiNSIuo3hr+9A4iQwB8CmjyCjK+zeaHQyvxZQ2Q8kJSWRlJR02v1Tpkxh+PDhpKamApCRkcEnn3zCrFmzGDNmDAA5OTlVS1sJkydPZsKECdV+HRERkVMymaBLGkR2hg+GwNFfYOYNcMNTED/CvV88yqN9amw2Gxs2bCAhIeF/FzCbSUhIYPXq1Z681FmNHTuWgoKCsldeXl6NXl9ERASACzrCPSvgsr7gtMPSMTD/Djh+xOhkfsejRc3BgwdxOByEh4eX2x4eHs6+ffsqfJ6EhARuueUWlixZwgUXXFCuIDrTvj8KDg4mLCyMefPm0blzZ3r16lW1LyUiInKu6jSEwW9C0vNgCYItH0NGN8jTYsueVOnmp5rw+eefV2nfqaSlpZGWlkZhYSENGjQ412giIiJVYzJB/N0QGQfvp8KRXJjdG3qNhy6j3KOn5Jx49E+wSZMmWCwW8vPzy23Pz8+nefPmnrxUhVmtVqKjo4mLizPk+iIiIuVEXAn3fAVtB4CzFJY/Du/8BY4dNjqZz/NoURMUFETHjh3JzMws2+Z0OsnMzKRLly6evFSFpaWlsWnTJtat0yM+ERHxEiENYNBs6PsSWIJh+2eQcS38UrP9T/1NpYuaoqIicnJyykYw5ebmkpOTw65duwBIT09nxowZzJ07l82bNzNy5EiKi4vLRkOJiIgI7uaoq4bC8Ew4rzUU7oY5N8KKF8HpNDqdT6p0n5r169fTs2fPsvfp6ekApKSkMGfOHAYPHsyBAwcYN24c+/btIzY2lqVLl57UebimWK1WrFYrDofDkOuLiIicUfN2cHcWfJwO378HmU/CzpUw4FUIbWp0Op9S6aKmR48euM4y3fOoUaMYNWpUlUN5kjoKi4iI1wuuDze/Bhd1gyUPwY5Md3PUoJkQda3R6XyGulqLiIh4A5MJOtwJd38JTS6Fon0wNxmyngWnWhsqwu+LGo1+EhERn9LscndhE3sHuJyQNQnmDYDf8s/+2VrO74sajX4SERGfE1QP+lvd/WoC60LuV+7mqJ+zjE7m1fy+qBEREfFZMX+Bu7+CZm2heD+80R++eBocpUYn80p+X9So+UlERHxa0zbuYd8dUgAXZD8Hb/SDwr1GJ/M6fl/UqPlJRER8XmAd6DcNBs6EoFD4ZSVkXAPbK7d0kL/z+6JGRETEb7QbBPdku+e2OXYI3hoInz8BDrvRybyCihoRERFfcl4rGPY5xN3lfv/1S+6ZiAt+NTaXF/D7okZ9akRExO8EhsCNL8ItcyE4DPL+zz06autSo5MZyu+LGvWpERERv9W2v7s5KuJKOH4E3hkMn/0TSm1GJzOE3xc1IiIifq3xRTD0M4gf6X6/ejrM7g1HfjE2lwFU1IiIiPi6gGBIegYGvwUhDWD3Bni1K2z+j9HJapSKGhEREX9xeV8Y8TW0uApOFMD8O2DJw1BaYnSyGuH3RY06CouISK3SsCUMXQpX3+d+v/ZVmHkDHP7Z2Fw1wO+LGnUUFhGRWscSCDdMhNvegzqNYG8OvNodflxodLJq5fdFjYiISK3VJtHdHBXZGUoK4f0h8HE62E8YnaxaqKgRERHxZw0ugCGfwLXp7vfrZ8LrCXDwJ2NzVQMVNSIiIv7OEgAJ4+GOBVC3CeR/D691h+/eNzqZR6moERERqS1aJ7ibo6K6gq0IPrwLPhoFtmNGJ/MIFTUiIiK1Sdj58LePoPsjgAm+mQev94IDW41Ods5U1IiIiNQ2Zgv0fNRd3NRrBvs3wWs9IOdto5OdE78vajRPjYiIyGlc3B1GroSLe4D9GCwaCQtHgK3Y6GRV4vdFjeapEREROYPQZnDHh9DzMTCZ4dt33E9t8n80Olml+X1RIyIiImdhtkD3hyDlP1D/fDi4DWZcBxvmgstldLoKU1EjIiIiblHXukdHtU6A0hPwn9Gw4C4o+c3oZBWiokZERET+p14TuO19SHgCTBb44QP3Egt7vzM62VmpqBEREZHyzGa49u+QugTCWsDhHe5ZiNe97tXNUSpqRERE5NRadnY3R7VJAkcJfPIP9/pRJwqMTnZKXlnUDBgwgEaNGjFo0KBy2/Py8ujRowfR0dG0b9+e99/3r+mdRUREvE7dxvDXd+CGp8EcAJsWwavdYPdGo5OdxCuLmvvvv5833njjpO0BAQFMnTqVTZs2sWzZMh544AGKi31zLL2IiIjPMJng6lEw9DNo0BKO7ISZN8CaDK9qjvLKoqZHjx7Ur1//pO3nn38+sbGxADRv3pwmTZpw+PDhGk4nIiJSS11wFYzIhsv6gtMOSx+B+XfA8SNGJwOqUNRkZ2eTnJxMREQEJpOJRYsWnXSM1WolKiqKkJAQ4uPjWbt2rSeylrNhwwYcDgeRkZEeP7eIiIicRp1GMPhNSHoOLEGw5WPI6Aa/rjc6WeWLmuLiYmJiYrBarafcP3/+fNLT0xk/fjwbN24kJiaGxMRE9u/fX3ZMbGwsV1xxxUmvPXv2VCjD4cOH+dvf/sZrr7122mNKSkooLCws9xIREREPMJkg/h4YtgwaRUHBLpiVCKv+bWhzVEBlP5CUlERSUtJp90+ZMoXhw4eTmpoKQEZGBp988gmzZs1izJgxAOTk5FQtLe5ipX///owZM4arr776tMdNnjyZCRMmVPk6IiIichYRV8I92fCf++HHhe4+NlfeCXUaGhLHo31qbDYbGzZsICEh4X8XMJtJSEhg9erV53x+l8vFkCFDuO6667jzzjvPeOzYsWMpKCgoe+Xl5Z3z9UVERORPQhrAoNlw4xQYNNOwggY8XNQcPHgQh8NBeHh4ue3h4eHs27evwudJSEjglltuYcmSJVxwwQVlBdHKlSuZP38+ixYtIjY2ltjYWL7//vtTniM4OJiwsDDmzZtH586d6dWrV9W/mIiIiJyeyQRxw9zz2hio0s1PNeHzzz8/5fZrr70Wp9NZqXOlpaWRlpZGYWEhDRo08EQ8ERER8UIefVLTpEkTLBYL+fn55bbn5+fTvHlzT16qwqxWK9HR0cTFxRlyfREREakZHi1qgoKC6NixI5mZmWXbnE4nmZmZdOnSxZOXqrC0tDQ2bdrEunXrDLm+iIiI1IxKNz8VFRXx008/lb3Pzc0lJyeHxo0b07JlS9LT00lJSeGqq66iU6dOTJ06leLi4rLRUCIiIiLVodJFzfr16+nZs2fZ+/T0dABSUlKYM2cOgwcP5sCBA4wbN459+/YRGxvL0qVLT+o8XFOsVitWqxWHw2HI9UVERKRmVLqo6dGjB66zTKwzatQoRo0aVeVQnqSOwiIiIrWDV679JCIiIlJZfl/UaPSTiIhI7eD3RY1GP4mIiNQOfl/UiIiISO3g90WNmp9ERERqB78vatT8JCIiUjv4fVEjIiIitYNXLmhZHX6fW6ewsLBazu8sOVYt5xWRqqvIz7snf3Y99feLpzJV199356Ii380bc4txfv//4Wxz5AGYXBU5yof9PqOwzWZjx44dRscRERGRKsjLy+OCCy444zF+X9T8zul0smfPHurXr4/JZCrbHhcXd9r+Nqfb9+fthYWFREZGkpeXR1hYmOfDV8KZvk9Nnq8yn6vIsVW5T6fbd6pt/noPfeH+nWm/fgZ1D41QG++hN/8udLlc/Pbbb0RERGA2n7nXTK1pfjKbzaes8CwWy2n/8E+373Tbw8LCDP9hPNP3qcnzVeZzFTm2KvfpdPvOdLy/3UNfuH9n2q+fQd1DI9TGe+jtvwsrusxRre8onJaWVul9Z/qM0Tydrarnq8znKnJsVe7T6fZ58/0Dz+bzhft3pv36GdQ9NEJtvIf+8ruw1jQ/VaffF8ssKCgw/F8YUjW6h75N98/36R76Pm+4h7X+SY0nBAcHM378eIKDg42OIlWke+jbdP98n+6h7/OGe6gnNSIiIuIX9KRGRERE/IKKGhEREfELKmpERETEL6ioEREREb+gokZERET8goqaGnDs2DEuvPBCHnzwQaOjSBVERUXRvn17YmNj6dmzp9FxpApyc3Pp2bMn0dHRtGvXjuLiYqMjSQVt3bqV2NjYsledOnVYtGiR0bGkkl566SXatm1LdHQ0o0ePrtDilFVRa5ZJMNLTTz9N586djY4h52DVqlWEhoYaHUOqaMiQIUycOJGuXbty+PBhzYXiQy699FJycnIAKCoqIioqiuuvv97YUFIpBw4cYPr06fz4448EBgbSrVs31qxZQ5cuXTx+LT2pqWbbt29ny5YtJCUlGR1FpFb6/S/Srl27AtC4cWMCAvTvOV+0ePFievXqRb169YyOIpVUWlrKiRMnsNvt2O12mjVrVi3XUVFzBtnZ2SQnJxMREYHJZDrlI0+r1UpUVBQhISHEx8ezdu3acvsffPBBJk+eXEOJ5c88cQ9NJhPdu3cnLi6Ot956q4aSy+/O9R5u376d0NBQkpOT6dChA5MmTarB9OKJn8HfvffeewwePLiaE8ufnes9bNq0KQ8++CAtW7YkIiKChIQEWrVqVS1ZVdScQXFxMTExMVit1lPunz9/Punp6YwfP56NGzcSExNDYmIi+/fvB+Cjjz6iTZs2tGnTpiZjyx+c6z0E+Prrr9mwYQOLFy9m0qRJfPfddzUVXzj3e1haWsqKFSt4+eWXWb16NcuXL2f58uU1+RVqNU/8DIJ7XaFVq1bRp0+fmogtf3Cu9/DIkSN8/PHH7Ny5k927d7Nq1Sqys7OrJ6xLKgRwLVy4sNy2Tp06udLS0sreOxwOV0REhGvy5Mkul8vlGjNmjOuCCy5wXXjhha7zzjvPFRYW5powYUJNxpY/qMo9/LMHH3zQNXv27GpMKWdSlXu4atUq1w033FC2/7nnnnM999xzNZJXyjuXn8E33njDdfvtt9dETDmDqtzD9957z3XvvfeW7X/uuedczz77bLXk05OaKrLZbGzYsIGEhISybWazmYSEBFavXg3A5MmTycvLY+fOnbzwwgsMHz6ccePGGRVZ/qQi97C4uJjffvsNcHdS/OKLL2jbtq0heeVkFbmHcXFx7N+/nyNHjuB0OsnOzubyyy83KrL8QUXu3+/U9OSdKnIPIyMjWbVqFSdOnMDhcJCVlcWll15aLXnUW66KDh48iMPhIDw8vNz28PBwtmzZYlAqqYyK3MP8/HwGDBgAgMPhYPjw4cTFxdV4Vjm1itzDgIAAJk2aRLdu3XC5XNxwww307dvXiLjyJxX9e7SgoIC1a9eyYMGCmo4oZ1GRe9i5c2f69OnDlVdeidlsplevXvTr169a8qioqSFDhgwxOoJUwcUXX8y3335rdAw5R0lJSRqB6MMaNGhAfn6+0THkHDz99NM8/fTT1X4dNT9VUZMmTbBYLCf9oOXn59O8eXODUkll6B76Pt1D36b75/u87R6qqKmioKAgOnbsSGZmZtk2p9NJZmZmtUwoJJ6ne+j7dA99m+6f7/O2e6jmpzMoKirip59+Knufm5tLTk4OjRs3pmXLlqSnp5OSksJVV11Fp06dmDp1KsXFxaSmphqYWv5I99D36R76Nt0/3+dT97BaxlT5iS+//NIFnPRKSUkpO+bf//63q2XLlq6goCBXp06dXGvWrDEusJxE99D36R76Nt0/3+dL99DkclXTqlIiIiIiNUh9akRERMQvqKgRERERv6CiRkRERPyCihoRERHxCypqRERExC+oqBERERG/oKJGRERE/IKKGhEREfELKmpERETEL6ioERGvNGTIEEwmEyaTiUWLFhmWIysrqyxH//79DcshImenokZEqt0fC5Q/vnr37n3Gz/Xu3Zu9e/eSlJRUbvuXX35J3759adq0KSEhIbRq1YrBgweTnZ1d4Uzt2rVjxIgRp9w3b948goODOXjwIFdffTV79+7l1ltvrfC5RcQYKmpEpEb8XqD88fXOO++c8TPBwcE0b96c4ODgsm0vv/wyvXr14rzzzmP+/Pls3bqVhQsXcvXVV/P3v/+9wnmGDRvGu+++y/Hjx0/aN3v2bPr160eTJk0ICgqiefPm1KlTp+JfVkQMoaJGRGrE7wXKH1+NGjWq1Dl27drFAw88wAMPPMDcuXO57rrruPDCC2nfvj33338/69evL3f8119/TdeuXalTpw6RkZGMHj2a4uJiAO644w6OHz/OggULyn0mNzeXrKwshg0bdm5fWERqnIoaEfEZCxYswG638/DDD59yv8lkKvvvHTt20Lt3bwYOHMh3333H/Pnz+frrrxk1ahQATZo04aabbmLWrFnlzjFnzhwuuOACbrjhhur7IiJSLVTUiEiN+PjjjwkNDS33mjRpUqXOsW3bNsLCwmjevHnZtgULFpQ75/fffw/A5MmTuf3223nggQe45JJLuPrqq5k2bRpvvPEGJ06cANxNUFlZWeTm5gLgcrmYO3cuKSkpmM3661HE1wQYHUBEaoeePXvyyiuvlNvWuHHjSp/nj09jABITE8nJyWH37t306NEDh8MBwLfffst3333HW2+9VXasy+XC6XSSm5vL5ZdfzvXXX88FF1zA7NmzefLJJ8nMzGTXrl2kpqZW4RuKiNFU1IhIjahXrx6tW7c+p3NccsklFBQUsG/fvrKnNaGhobRu3ZqAgPJ/nRUVFXHPPfcwevTok87TsmVLAMxmM0OGDGHu3Lk88cQTzJ49m549e3LxxRefU04RMYaer4qIzxg0aBCBgYE8++yzZz22Q4cObNq0idatW5/0CgoKKjsuNTWVvLw8PvzwQxYuXKgOwiI+TE9qRKRGlJSUsG/fvnLbAgICaNKkSYXP0bJlS1588UXuv/9+Dh8+zJAhQ7jooos4fPgwb775JgAWiwWARx55hM6dOzNq1Cjuuusu6tWrx6ZNm1i+fDnTp08vO+dFF13Eddddx913301wcDA333yzB76tiBhBT2pEpEYsXbqU888/v9zr2muvrfR57rvvPpYtW8aBAwcYNGgQl1xyCX369CE3N5elS5fSrl07ANq3b89XX33Ftm3b6Nq1K1deeSXjxo0jIiLipHMOGzaMI0eOcNtttxESEnLO31VEjGFyuVwuo0OIiPzZkCFDOHr0qKFLJPyRt+URkZPpSY2IeK3fh4F//PHHhmVYsWIFoaGh5UZRiYh30pMaEfFK+/fvp7CwEIDzzz+fevXqGZLj+PHj7N69G3CPtPrjHDki4l1U1IiIiIhfUPOTiIiI+AUVNSIiIuIXVNSIiIiIX1BRIyIiIn5BRY2IiIj4BRU1IiIi4hdU1IiIiIhfUFEjIiIifuH/AVfL/u2dJoAaAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "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": "e686bbee",
   "metadata": {},
   "source": [
    "The `BrokenPowerLaw` class is also available and behaves in a very similar way."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc6f180b",
   "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": "e8f36727",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-08T10:36:25.967997Z",
     "iopub.status.busy": "2024-11-08T10:36:25.967771Z",
     "iopub.status.idle": "2024-11-08T10:36:25.971325Z",
     "shell.execute_reply": "2024-11-08T10:36:25.970800Z"
    }
   },
   "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": "e105c6d2",
   "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": "be5a3fdc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-11-08T10:36:25.973398Z",
     "iopub.status.busy": "2024-11-08T10:36:25.973012Z",
     "iopub.status.idle": "2024-11-08T10:36:25.976947Z",
     "shell.execute_reply": "2024-11-08T10:36:25.976446Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<icecube_tools.source.flux_model.PowerLawFlux at 0x7fe3c4f9f190>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diffuse_source.flux_model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "516b232a",
   "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": "b65dc4ef",
   "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.20"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}