diff --git a/docs/tutorials/absorption_in_BLR_with_cubepy.ipynb b/docs/tutorials/absorption_in_BLR_with_cubepy.ipynb deleted file mode 100644 index cabf7f7..0000000 --- a/docs/tutorials/absorption_in_BLR_with_cubepy.ipynb +++ /dev/null @@ -1,423 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/home/pgliwny/anaconda3/envs/agn/lib/python3.8/site-packages/scipy/__init__.py:138: UserWarning: A NumPy version >=1.16.5 and <1.23.0 is required for this version of SciPy (detected version 1.23.3)\n", - " warnings.warn(f\"A NumPy version >={np_minversion} and <{np_maxversion} is required for this version of \"\n" - ] - } - ], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2\n", - "\n", - "import numpy as np\n", - "import astropy.units as u\n", - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "from agnpy.absorption import Absorption\n", - "from agnpy.targets import SphericalShellBLR, lines_dictionary" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This notebook presents how to use the multi-layer BLR model, which considers 27 BLR lines." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": {}, - "outputs": [], - "source": [ - "SUM_RATIO_LINES_BLR = 30.652\n", - "\n", - "def get_absor_blrs(E, L_disk, R_line, r_blob_bh, z):\n", - " \n", - " \"\"\"\"\n", - " Function to calculate absorption in BLR shells.\n", - " Radius and luminosity are normalized to Hbeta\n", - " Parameters\n", - " ----------\n", - " E: numpy array \n", - " Energy array with unit\n", - " L_disk : astropy.units.Quantity\n", - " Luminosity of the disk whose radiation is being reprocessed by the BLR\n", - " R_line astropy.units.Quantity\n", - " radius of the BLR spherical shell\n", - " r_blob_bh : astropy.units.Quantity\n", - " distance between the Broad Line Region and the blob\n", - " z : float\n", - " redshift of the source\n", - " \"\"\"\n", - " nu = E.to(\"Hz\", equivalencies=u.spectral())\n", - " XI_LINE = 0.1\n", - " blrs=dict()\n", - " ec_blrs=dict()\n", - " tau_z_all = np.zeros(nu.shape)\n", - " for line in list(lines_dictionary.keys()):\n", - " xi_line_this=lines_dictionary[line]['L_Hbeta_ratio']\n", - " xi_line = (xi_line_this*XI_LINE)/SUM_RATIO_LINES_BLR\n", - " R_line_this=R_line*lines_dictionary[line]['R_Hbeta_ratio']\n", - " blrs[line] = SphericalShellBLR(L_disk, xi_line, line, R_line_this)\n", - " ec_blrs[line] = Absorption(blrs[line], r=r_blob_bh, z=z)\n", - " tau_z_this = ec_blrs[line].tau_blr_cubepy(nu, eps_abs=1e-6)\n", - " tau_z_all+=tau_z_this\n", - " A_blr = get_absorbtion_tau(np.array(tau_z_all))\n", - " return A_blr\n", - "\n", - "def get_absorbtion_tau(tau):\n", - " \"Function return absorption\"\n", - " return np.exp(-tau)\n", - "\n", - "def log_parabola(energy, amplitude, reference, alpha, beta):\n", - " \"Log Parabola model from gammapy\"\n", - " xx = energy / reference\n", - " exponent = -alpha - beta * np.log(xx)\n", - " return amplitude * np.power(xx, exponent)" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": {}, - "outputs": [], - "source": [ - "E = np.logspace(-1, 1, 20)*u.TeV\n", - "lp = log_parabola(\n", - " E,\n", - " amplitude=1e-12/((u.cm**2)*(u.s)*(u.TeV)),\n", - " reference=1 * u.TeV,\n", - " alpha=2.3,\n", - " beta=0.5\n", - " )\n", - "sed = E*E*lp" - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/home/pgliwny/anaconda3/envs/agn/lib/python3.8/site-packages/astropy/units/quantity.py:611: RuntimeWarning: invalid value encountered in sqrt\n", - " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "CPU times: user 20.9 s, sys: 904 ms, total: 21.8 s\n", - "Wall time: 20.7 s\n" - ] - } - ], - "source": [ - "%%time\n", - "# Set BLR parameters\n", - "L_disk=6e45*u.Unit(\"erg s-1\")\n", - "R_line = 2.45*1e17*u.cm\n", - "r_blob_bh = 1e18*u.cm\n", - "z = 1\n", - "A = get_absor_blrs(E, L_disk, R_line, r_blob_bh, z)" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAGtCAYAAAB5r18AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABTXElEQVR4nO3deZxN9R/H8dd3Zuz7EsoSESlKSdkq9Csq28i+Za1EiBK/pI2oLKmGIWQpe2SJUPGzpLK0kUiWkCJLjJ35/v44Y7LMmHvHvffce+f9fDzuY+Z+z/d8z+dMp/txzj3n8zXWWkRERCR5EW4HICIiEuyULEVERFKgZCkiIpICJUsREZEUKFmKiIikIMrtANySN29ee80115AlS5Yr9jt27NgV+6S0PBS5vU/+2r4vxk3tGN6u52l/T/qlxWMY3N2vcDyGvV3X130DdRyvW7fub2vtNZctsNamyVf58uXt0qVLbUpS6uPJGKHG7X3y1/Z9MW5qx/B2PU/76xhOnpv7FY7HsLfr+rpvoI5jYK1NImfoMqyIiEgKlCxFRERSoGQpIiKSAiVLERGRFChZioiIpEDJUkREJAVKliIiIilQshQREUmBkqWIiEgKlCxFRERSoGQpIiKSAiVLERGRFChZioiIpEDJUkREJAVKliIiIilIs5M/iwSLP/74g1WrVrFy5UpWrlzJ5s2byZQpE9deey25c+cmT5485MmTJ8nfd+zYwZ9//knu3LlJnz6927siEraULEUCKD4+nl9++SUxMa5atYpt27YBkDlzZipWrEiHDh3YsmULGTJk4MCBA/zyyy8cPHiQAwcOcObMmWTHzpYt22UJtWzZspQoUSJQuycStkIyWRpjbgBeAHJYaxsmtNUHHgHyATHW2sXuRSjiOHXqFOvWrWPKlCkMHTqUVatWcfDgQQDy5ctH1apV6dKlC1WrVqVcuXKkS5cOgGXLllGtWrWLxrLWEhcXl5g4Dxw4wIoVKyhQoAAHDhy4qP3gwYNs376dadOmYYxh7NixtG3blvr165MpU6ZA/xlEQl7Ak6UxZhxQG9hnrS1zQXstYDgQCYyx1g5Kbgxr7TagvTFm5gVtnwCfGGNyAYMBJUsJuEOHDvHVV18lnjmuWbOGU6dOAVCqVCmio6OpWrUqVapUoUSJEhhjPB7bGEO2bNnIli0b119/PQDp0qW7LKleaPv27bzyyissW7aM5s2bkzNnTpo1a0bbtm258847vdq+SFrmxpnleOA9YOL5BmNMJBADPADsBtYYY+biJM6Bl6zfzlq77wrj900YSyQgdu7cydChQ/nyyy/ZsGEDAFFRUZQvXz7xrNFaS3R0dMBjK1asGG3atGHcuHEsXbqUDz74gA8++ICRI0dyyy230K5dO4oVKxbwuERCjbHWBn6jxhQF5p8/szTGVAJettbWTHjfB8Bae2mivHScmRdchjXAIGCJtfbzZPo/DjwOkD9//vJjxowha9asV4w1Li7uin1SWh6K3N4nf23fF+NeOMY///zD5MmTmT17NgDlypWjbNmylC1blptuuomMGTOmetue9vek36V94uLi+PLLL/nss8/YtGkTkZGRVKxYkYceeoi7776bqKiQ/HbmMm4ex6FyDPtzXV/3DdRncfXq1ddZa++8bIG1NuAvoCiw4YL3DXEuvZ5/3wp47wrr5wFigd+APgltXYF1Ce1PphRD+fLl7dKlS21KUurjyRihxu198tf2fTHu0qVL7fHjx+2gQYNsjhw5rDHGtmnTxv7+++8+3ban/a/2GN64caNt0qSJzZ8/vwVsvnz5bM+ePe2GDRs8CzSIuXkcB/sxHIh1fd03UJ/FwFqbRM4Ilucsk/riJNlTXmvtAWvtk9ba4jbh7NNa+461tnxCe6zfIpU069y5cyxYsIAbb7yR3r17U7VqVX744Qc++OADChcu7HZ4qXLzzTfz5JNPsmvXLubOnUvlypUZPnw4ZcqU4e677yY2NpbDhw+7HaaI64IlWe4GLvy0KQT84VIsIhex1jJv3jxuvfVW3nrrLQoVKsSyZcuYP38+ZcuWdTs8n0iXLh116tRh9uzZ7Nmzh6FDh3L8+HE6derEtddeS/Pmzdm0aZPbYYq4JliS5RrgRmNMMWNMeqApMNflmERYvXo19913H3Xr1uXMmTO8/PLLiW3hKl++fDzzzDP8+OOPrFmzhnbt2rFgwQJuv/12hg0bRnx8vNshigRcwJOlMWYKsBooZYzZbYxpb609C3QBFgGbgOnW2o2Bjk3kvM2bN/Poo49SuXJltmzZwsiRI9m4cSP33XdfmnncwhjDnXfeSUxMDJs3b6ZmzZr06NGDGjVqsGPHDrfDEwmogCdLa20za+211tp01tpC1tqxCe0LrLUlE76HHBDouEQA9u7dy5NPPsktt9zC4sWLefXVV9m6dStPPvlkYsGAtCh//vx88sknjBs3jvXr13Prrbcybty48zfciYS9YLkMK+KqI0eO8OKLL1KiRAnGjh3LU089xW+//caLL74Ydo8GpZYxhrZt2/Ljjz9yxx130L59e+rVq8dff/3ldmgifqdkKWna6dOneeeddyhevDj9+/enbt26/PLLL7zzzjvky5fP7fCCUtGiRfnyyy8ZNmwYixcvpkyZMnz88cduhyXiV0qWkmbNmTOHm266iW7dunHrrbeyZs0apkyZQvHixd0OLehFRETQvXt31q9fz/XXX0/Dhg1p3bq1HjORsOVKBR83GWPqAHUKFizYMTY2NmiqRgQTt/fJ39VPrLVMnDiR8ePHc8MNN/DEE09QoUIFj27cSW1swVTB52pju9TZs2f58MMPmTRpEnny5OH555+nfPnyqR7PV1TBx/djqIJPGnypgk/y3N4nf1Y/OXHihG3RooUFbOvWre3JkycDEtvSL7+0Nj7eefPHH9Z+8421n39u7axZ1o4fb+2IEf92njHD7mje3NoBA6wdPtzasWOtnTHj3+U7dli7ebO1e/bY5fPnW3v27FXF7Ku/97fffmtLlSplAdulSxd77Ngxn4ybWqrg4/sx0nIFn/AoAinigcOHD3P//ffz1VdfMWDAAPr06eObx0D274d334WePSFHDpgwAWJj4cgR53X0KPed/z1rVhg8GIYOvXycxx+HyEj48kuKTJ0KFz7PmDMnNGzo/P7cczBjBgD3nF9evDhs3er83qULrF/vbKtUKbKXKgX33Qd+fuSlQoUKfPfdd/Tp04fhw4ezePFiJk6cyN133+3X7YoEgpKlpAk///wzTz31FIcOHWLatGk0btz46gc9cwZiYuDll51E2K6dkyyjoiBbNihYELJnh2zZ2HnoEEXPJ6s2baB6dWfZha+IhFsIRozgf40aUa1yZYiLc14nT/673WeegXr1IC6OrT/8QIn8+eHCOSpz5XIS5dGjMGYMd5w8CZ9+CgsXXv0+pyBTpky8/fbb1K1bl7Zt21K5cmX++9//8uKLL5I+fXq/b1/EX5QsJewtWbKERo0aERERwbJly3xzprN4MXTvDps2Qc2aMGwYFC3qLGvRwnldYMeyZRTNksV5U7as87oSYyBDBueVJ8/FyypVcl7A7mXLKHHpfJavvfbv70ePsun11yldrpzz/sQJJ966daFZMyeh+0GNGjX48ccf6d69O/379+fTTz9l4sSJlClTJuWVRYKQ7oaVsBYbG8tDDz1EkSJFGDFihG8S5Zkz0LkznD4Nc+c6Z2ylS1/9uP6QLRt/1awJTZo47/fscc5Sn3sOCheGGjUo8Omnzpmxj+XIkYMPPviA2bNns3v3bsqXL8/gwYM5d+6cz7cl4m9KlhKWzp07xzPPPEOnTp2oWbMmK1eupECBAqkf8OhR6N+fiBMnIF0657Lmxo1Qp47fvwv0qRIl4NtvYfNm6NcPdu3ipsGDnffgfP964SVfH6hfvz4bNmzg4Ycf5rnnnqNGjRr8/fffPt2GiL8pWUrYOXr0KPXr1+ftt9+ma9euzJkzh+zZs6dusPh4mDQJSpWCF18k97ffOu0lSzqXSENVyZLOd61btrB21Ci4M+FO+f/+FwoUgA4d4MsvwUdngfny5WPWrFmMHz+eb775hvvvv18JU0KKkqWElV27dnHPPfewYMEC3nvvPYYPH05UVCq/mv/2W6hSBVq3di5Zrl7N3+E224gxxJUs+e/ZcYsWzs1D06bB/ffD9dfD66/7aFOGxx57jLlz57JlyxZq1KjB/v37fTK2iL8pWUrYWLNmDXfddRfbtm3j008/pXPnzlc3YO/esGMHjB8Pq1dDxYq+CDO4VavmPPry118wdSrccQdceAa4b99Vb+LBBx9k3rx5/Prrr9SoUYN9PhhTxN+ULCUszJo1i/vuu48MGTLw1VdfUatWLe8HOX0ahgxxboIBJ0lu3gyPPfbvYx1pRebMzk1Bc+c6fxNwfi9WzHlO9OzZqxr+P//5D/Pnz+e3335TwpSQkMY+ASTcWGsZNGgQjz76KLfddhvffPON948nWAvz50OZMvDss84lSIAiRZznH9O685dob78d/vMf507aChVgzZqrGvb+++9n/vz5bNu2jerVq2v2Eglqqg0bJPUIg4nb++Tp9s+cOcPQoUP57LPPqFGjBr169SLDFW66SWrczL//TvGYGPJ8+y3HCxdma5cuHLzrrquO7WrXC9rasNaSd8UKbnznHdIfOsSO1q3Z+dhjnq2bjPNVfwoUKMDQoUPJnTv3VY13nmrD+n4M1YZNgy/Vhk2e2/vkyfb//vtve99991nA9uvXz8afr7vq7bjt2lmbPbu1Q4ZYe+qUT2LzxXqe9nftGD582NqnnrJ24kTnvQd//5RiyJw5sy1durTdu3fvVY114ZhuUW3Y8KsNq8uwEnK2bNlCpUqVWL16NZMmTeKVV17xrsbrqVPw55/O74MGwa+/Qo8eoHJsnsuRwyn116qV8/7dd+HRR//9vtdL1apVY8GCBezcuZPq1auzd+9eHwYrcvWULCWkfP3111SsWJFDhw7xxRdf0LJlS+8GOHwYatVyvns7fRquuQY0ybNvLFjgVDKKiUnV85n33XcfCxcuZNeuXUqYEnSULCVk/Prrr9SuXZvcuXPzzTffULVqVa/Wz7B/P9xzD6xaBX366EzSl7p2hQ0bnMdrunRxnk/dsMHrYe69914WLlzI7t27qVatGn/88YcfghXxnpKlhIS///6bhx9+GIDPPvuMG264wbsBNmzg9s6dYedO5wzokkLn4gPFi8OiRU7Fo23b4MCBVA1zzz338Nlnn/HHH39QrVo19qTy0q6ILylZStA7efIk9erVY9euXcydO5cSJUp4N4C10KULJj4eli93LsGKfxgDLVs6xRzOVzt66y1nlhYvVK1alc8++4y9e/dSvXp1JUxxnZKlBLX4+HjatGnDV199xaRJk6hcubK3Azgf4B99xPr33oPzU1WJf2XO7Pw8dcop7lCzpnM270XxgSpVqrBo0SL+/PNPqlWrxu7du/0Tq4gHlCwlqL3wwgtMmzaNN954g0aNGnm38tChzh2aZ89CwYKcuppZRyR1MmSAdeucGU5mzICbboJx45yzfQ9UrlyZRYsW8ddff1GtWjV27drl54BFkqZkKUHr/fffZ9CgQTzxxBM899xznq8YH+88CtKzJ0RGXnVpNrlKGTPCK6/ADz84VZKeegq2bvV49UqVKrF48WL2799PtWrV+P333/0YrEjSlCwlKC1atIhOnTpRq1Yt3nvvPc+fozx1Cpo1g2HDnDs0p01zPqzFfaVLw7JlTpm8G2902jw8w6xYsSJLlizhwIEDSpjiCiVLCTq//fYbjRo1okyZMkyfPt27KbZatoTp052bSt5+2zmzlOAREQFlyzq/jx/vFGs/dcqjVe+66y6WLFnCwYMHqVatGjt37vRfnCKXULKUoLJnzx769OlD9uzZmT9/PtmyZfNugGefhcmTnZ/eVPWRwPvnH+d7zHr14Ngxj1apUKECn3/+OYcOHaJatWrs2LHDvzGKJFAh9SAp3htM3Nqn48eP061bN/bs2cM777zj8SMiWbZtI9fatexu3PiK/dwsQh02hdR9rMCCBZQaMoQjN9/MTwMHctbDODZv3syzzz5Ljhw5SO7/YxVS9/0YKqSeBl8qpJ48N/bpzJkz9qGHHrKRkZF20KBBnq/45ZdOIfTrrrP24MErdnWzCHXYFVL3pRkzrE2Xztpy5aw9ftzj1VauXGmjoqJs/fr1kyykr0Lqvh9DhdRFXGSt5emnn2bhwoWMGDGCu+++27MVp051nt8rVAhWr4ZcufwbqPhHw4Ywbx40bgyZMnm8WpUqVXjzzTf55JNPGDx4sB8DFNF3lhIEhgwZQmxsLM8//zyPP/64ZysNG+bc9VqpEqxc6UzULKGrZk2nXi/A2rWwZYtHq3Xv3p2GDRvSp08fli9f7scAJa1TshRXzZw5k+eee47GjRvz+uuve75i7tzOmciiRTqjDCfnzjnTft1zD3z/fYrdjTGMHTuW4sWL06RJE/48P/WaiI8pWYprVq9eTatWrahcuTLjx48nIiKFw9FaZ+5JgMcecy7D6hnK8BIZCZ984swIU60afPVViqtkz56dmTNn8s8//9C0aVPOqgiF+IGSpbjit99+o27duhQqVIg5c+aQyZPvqgYPdirAfPed816PhoSnUqWcS+v58sEDD8CSJSmuUrZsWWJjY/nf//5H3759AxCkpDVKlhJwBw8e5OGHHyY+Pp4FCxaQN2/elFeaORN69YLoaLjtNv8HKe66/npYsQJKlID33/doldatW/P444/zxhtvMHfuXD8HKGmNkqUE1KlTp6hfvz47duxgzpw53Hi+7NmVfP218z1W5cpO1ZeULtdKeMifH/73P5g40Xl/+nSKqwwfPpw77riD1q1ba+Jo8Sl96kjAWGtp164dK1asYMKECVStWjXllf74A+rWhYIFYc4cfUeZ1uTM6fw3P3QI7roLhg+/YveMGTMyc+ZMjDG8/PLLnDx5MjBxSthTspSA6devH5MnT+b111+nadOmnq1UoAB06gQLFoAnl2slPGXODMWLQ/fu8OqrVyzAXqxYMSZNmsSvv/7K008/HbgYJawpWUpAjBs3jv79+9OhQwd69+6d8gqnT8OePc4l11degZIl/R+kBK8MGZwZZNq0gZdecqZfu0LCrF27Ns2bN2fMmDGMHz8+YGFK+PJiOgeR1Pnmm2944oknePDBBxkxYkTK021ZCx07OndB/vyzcylOJCoKxo6F7NmdohQ5cjiJMxnt2rVj7969dOrUidtvv53bdGOYXAWdWYpfHT9+nNatW3Pttdcybdo00qVLl+I610+c6NzU8eSTSpRysYgIZ+q1oUOhQ4crdo2MjGTKlCnkypWLhg0b8s8//wQmRglLSpbiV3369GHLli2MHz+enJ4kvg8/pNj48U7RgRdf9Hd4EoqMgWeecW76Onv2is9h5s+fn+nTp7N9+3batWuHTWOzLInvaIquIJkWJpj4ap/WrVvHs88+S3R0NF27dk2xf/aff6Zct24cvPlmNg4ejPXgLNQbmqLLN7EFk0LTp1Ni5Eh+eu01Dlxyd/WF+zV9+nRGjhxJp06daJzCVG6+oCm6NEVX2Lw0RVfyfLFPhw8ftoULF7YlS5a0x44d82ylI0es7dzZrpg796q3nxRN0XV12wpKJ05YW6GCtVmzWrthw0WLLtyv+Ph426BBAxsZGWlXrFjh97A0RZem6BLxSPfu3dmzZw8TJkwgc+bMV+584AAcOwbZssF773E2W7bABCmhL2NGmD0bsmaFevXg4MEkuxljGDduHMWKFaNx48b89ddfAQ5UQp2Spfjc3LlzGT9+PH369KFixYpX7nzyJNSpAw8/fMVHAUSSVbAgzJoFu3ZB27bJdsuRIwcff/wxhw4dolmzZiq4Ll5RshSf2r9/Px07dqRcuXL069fvyp3j450beVavhq5dVRhdUq9SJacUYgo3hd16662MHDmSpUuXpnx8ilxAz1mKz1hr6dSpE4cPH+bzzz8nffr0V16hb1+YPh3efBMefTQwQUr4atbs39937Ei2W5s2bVi1ahUDBw6kcuXK1K5d2/+xScjTmaX4zOTJk/n444959dVXKVu27JU7T5gAAwfCE0/As88GJkBJG0aNgtKlybZ5c7Jd3n33XW6//XZatWrF9u3bAxichColS/GJPXv20KVLFypXrsyzniS/e+5xar6+954uv4pvNWgA+fNzy4svwp9/JtnlfMF1gIYNG6rguqRIyVKumrWW9u3bc/r0aSZMmEBkZGTynffudW7kueEGGDHCKWEm4kvXXAOffEK6o0edy/unTiXZ7YYbbmDChAmsX7+ebt26BThICTVKlnLVRo0axaJFi3jrrbcoUaJE8h337oW774YePQIXnKRN5crxS69e8NVXcIWZR+rWrcvzzz/P6NGjmTx5cgADlFCjf9bLVdm6dSs9e/bkgQceoFOnTsl3PHbMeUTk4EFnImcRP9tfvTpERkIKE4z379+f5cuX061bN2rWrEmePHkCFKGEEp1ZSqqdO3eONm3akC5dOsaNG5f8bCLnzkHz5vDddzB1KtxxR2ADlbSrX79/75I9cSLJLlFRUYwaNYrDhw/Tq1evAAYnoUTJUlJtyJAhrFq1ivfee49ChQol37F/f5g715ktQrfpixs+/tiZE3XnziQXly1blp49ezJu3DiWL18e4OAkFChZSqr89NNPvPjiizRo0IAWLVpcuXOVKs4M9126BCQ2kcuUKQNHj0L9+nD8eJJd+vXrR9GiRXnyySc5ffp0YOOToKdkKV47ffo0rVu3JmfOnMTGxiZ/+fV8ObH//MeZrFePiIhbSpWCyZPhhx+gXbskSytmzpyZmJgYNm3axFtvveVCkBLMlCzFa6+99hrff/89o0eP5pprrkm605kzcP/9oA8dCRYPP+wUwpg2Dd54I5kuD9OwYUP69+/P1q1bAxygBDMlS/HKt99+y8CBA3nssceoV69e8h3/+19YvhwKFw5ccCIp6dULmjZ1LskmY/jw4aRLl46nnnpKk0VLIiVL8diJEydo3bo11113HcOHD0++4yefwODB8NRTzgeTSLAwBj76CAYMSLbLddddx+uvv86SJUuYOnVqAIOTYGbS2r+cjDF1gDoFCxbsGBsbGzSzcweT5Pbpvffe4+OPP2bw4MGUL18+yXUz/vEHdz7+OMcLFeK7d97BplRM3YvtXy03Z5n3dj1P+wfTDPPBJqX9yv7zzxR7/302vPYa5y7pd+7cOTp37sy+ffuYOHGi13+fcDyGvV3X130DdRxXr159nbX2zssWJDUjdFp4lS9fPqhm5w4mSe3Tl19+aQHbpUuXK688ZYq1+fJZu22bT7fvC27OMu/tep721zGcvBT363//szYqytpHHrH27NnLFq9bt85GRETYJ5980vfbTiU3j2Fv1/V130Adx8Bam0TO0GVYSdGRI0do27YtJUuW5I1kboxI1LQp/PYbFCsWmOBEUuvee+Gdd+DTT5OcB/OOO+6ga9eujBo1itWrV7sQoAQTJUtJ0TPPPMOuXbuYMGECmTNnTrrT9OnOd5UAYXhJT8LUk0/C4487d8l+/PFli1999VUKFizIE088wZkzZ1wIUIKFkqVc0bx58xg3bhy9e/emYsWKSXf6+Wdo29Z5ljI+PrABilwNY+Ddd6FChSSTZbZs2Xj33Xf56aefePvttwMfnwQNJUtJ1t9//03Hjh257bbbeOmll5LuFBcHDRs6Z5NTpkCEDikJMenTw2efOXfJJqF+/frUrVuXl19+mZ3JlMuT8KdPNklWly5dOHjwIBMnTiR9Une1WgtPPAGbNzvVUa67LvBBivhC7tzOWeauXbB27WWL3333XYwxdOnSRc9eplFKlpKkZcuWMW3aNPr27cutt96adKelS50k+corTrUekVBmrVM7tkmTy+rHFilShFdffZX58+cze/Zsd+ITVylZymXOnTtHjx49KFKkCM8991zyHatXh/nznWo9IqHOGKeYxrZtzj8AL9G1a1fKlSvH008/zZEjR1wIUNykZCmXWbJkCd999x2DBg0iU6ZMl3c4dMi5qccYeOQRfU8p4aN6dWjfHoYMceZfvcD5eS/37t3Li0k8aiLhTZ9ycpG4uDjGjBlDxYoVaZpUqTprnTtfq1SBw4cDHp+I3731FuTNCx07/jtzToK77rqLTp068d5777Fu3TqXAhQ3KFnKRd58800OHDjA0KFDk556a8gQmDMHXnoJcuYMeHwifpcrl/M4SblycOrUZYtff/118uXLxxNPPMG5c+cCH5+4QslSEu3atYvBgwdTo0YNKlWqdHmHlSuhd2949FHo1i3wAYoESqNGMGYMZMly2aIcOXIwfPhw1q1bR0xMjAvBiRuULCXRf//7X6y1dOzY8fKF+/c7dwkWKwZjx2oiZ0kb1q+HZ565bLLoRo0aUatWLV544QV2797tUnASSEqWAjjzVH744Yf06NGDAgUKXN4hRw4nWc6Y4fwukhZ89RW8/fZlBQuMMcTExHD27Fm66SpLmqBkKVhr6dGjB/nz56d3796Xdzh92qlyMnSo8z2OSFrRqRNUrOicXf7990WLbrjhBvr168esWbOYP3++SwFKoChZCjNnzmTVqlX079+fbNmyXbxw8WIoXRq2bHEnOBE3RUbC++/DP/9Ajx6XLe7Zsye33HILnTt35tixYy4EKIGiZJnGnTx5kl69enHrrbfStm3bixfu3g0tWjg3ORQq5E6AIm4rUwaefx4mTYJlyy5alD59ekaNGsXvv//OK0kUMpDwEeV2AOKu4cOHs2PHDj7//HMiIyP/XXDunDM35cmTzveUyU3NJZIWvPACXHMNVK582aIqVarQoUMHhg4dSosWLbjttttcCFD8TWeWadi+ffsYMGAAderU4f5La7sOGwarVsHIkVCqlDsBigSLjBmha1fnu/sk5rV84403yJ07N08++STxmqYuLClZpmH9+vXjxIkTvPXWWxcvsBY+/9wpKt2ihSuxiQSl9evhxhvhkuo9uXPnZsiQIXz99deMHj3apeDEn0xam27GGFMHqFOwYMGOsbGxZM2a9Yr94+LirtgnpeXBatu2bXTs2JH69evz9NNPX7QsLi6OrJkzE3nyJOdcuPzqr7+pL8ZN7Rjerudpf0/6hesxnBJ/7FdUXBwV2rThdO7crB85EnvBVxfWWnr27MmWLVsYOXIkhQsX9um2wd1j2Nt1fd03UMdx9erV11lr77xsgbU2Tb7Kly9vly5dalOSUh9Pxgg28fHx9oEHHrC5cuWyBw4cuHjhxx/bVTNnuhNYAn/9TX0xbmrH8HY9T/un1WPYE37br5kzrQVr33zzskW//PKLjYyMtA0aNPDLpt08hr1d19d9A3UcA2ttEjlDl2HToIULF7JkyRJeeuklcufO/e+CDRugWTOKjRnjXnAiwa5BA6hXz6mPvG3bRYtKlSpF27ZtmTdvHr///rtLAYo/KFmmMWfOnKFnz56ULFmSp5566sIF8NhjkCMH2554wr0ARYKdMRATA1FRTunHS/Tr1w+AV199NdCRiR8pWaYxo0aN4pdffuGtt94iXbp0/y4YNMi5eSE2ljOaTUTkygoWdG7y6d//skWFCxembt26jB8/ni0q5hE2lCzTkEOHDvHyyy9To0YN6tSp8++CH36AV1+F5s2dS0wikrIbb3TOMn//3Zlo4ALNmzcnY8aMvPTSSy4FJ76mZJmG9O/fn4MHD14+V2WhQs5Et++8415wIqEoLg5uvx26d7+oOXfu3HTv3p2pU6fyww8/uBOb+JSSZRrx66+/8u6779KuXbuLK4xYC3nywIgRzk8R8VzWrE6xgsmTYeHCixY9++yz5MyZk759+7oUnPiSkmUa8fzzz5MhQwb6X/gdy7p1UKUKbN/uXmAioa53b2eygU6dnDPNBDlz5qRXr17Mnz+f1atXuxig+IKSZRqwbNkyZs+eTZ8+ff6dq/LUKefu1507IVcudwMUCWUZMjgzk+zcCS++eNGirl27ki9fvsSJ1SV0KVmGuXPnztGjRw+KFCnCM8888++CV16BjRthzBjQ3a8iV6dKFXjqKThxwvlqI0GWLFno27cvy5Yt44svvnAxQLlaSpZhbuLEiXz33XcMGjSITJkyOY3ffgtvvAHt28NDD7kboEi4ePddiI117pC9wOOPP06RIkV0dhnilCzDWFxcHC+88AIVK1akadOm/y4YNMh5TmzIEPeCEwk3EQkfp99+S55VqxKbM2TIwEsvvcSaNWuYO3euS8HJ1VKyDGNvvvkme/fuZdiwYRc/KjJ5MixaBDlyuBecSLjq1YtSQ4ZcdLNP69atKVmyJH379uXcuXMuBieppWQZpnbt2sXgwYNp2rQpFStWdBo3b4ajR525+UqXdjdAkXA1aBDpDx1y5oRNEBUVxauvvsqGDRuYOnWqi8FJailZhqnz348MGjTIaTh+HOrUceaoFBH/qViR/ffcA2++eVFln0aNGnHbbbfx0ksvcSaJCaQluClZhqFvv/2WDz/8kB49enD99dc7jf/9L/z6K+gBaRG/296hg/MP1AEDEtsiIiIYMGAAv/32Gx988IGL0UlqKFmGGWstPXr0IH/+/PTu3dtp/N//YPhw6NIFqld3N0CRNOB4kSJOCbxrr72o/eGHH6ZSpUq8+uqrnDhxwp3gJFWULMPMzJkzWbVqFf379ydbtmzOTQbt2kHx4s5dsCISGEOGwPPPX9RkjOH1119nz549jBw50qXAJDWULMNIfHw8L7/8MjfffDNt27Z1Go8eheuvhw8+gCxZ3A1QJK2xFmbMgJ9+SmyqVq0aDzzwAAMHDuTo0aMuBifeULIMI3PmzOHnn3+mb9++REZGOo3XXgtffAH33ONucCJp0ZEj8MQT0KvXRc39+/fn77//5u2333YnLvGakmWYsNYyYMAASpQoQePGjZ3/STt2hL17L6soIiIBkiMH9OkDn30Gy5YlNt91113Ur1+fwYMHc/DgQffiE48pWYaJRYsWsW7dOnr37u2cVT73HIwbBzt2uB2aSNrWpYszZ+zzz19UN/a1117j6NGjvPnmmy4GJ55SsgwTAwYMoHDhwrRq1QoWL4bRo+HZZ6FSJbdDE0nbMmVyJi749luYNSuxuUyZMjRv3px33nmHvXv3uhigeELJMgwsX76clStX0qtXL9KfOOEUSC9d2vkfVETc17o1PPwwpE9/UfPLL7/MmTNnGHDB85gSnJQsw0D//v3Jnz8/7du3h5dfdr6nnDDBKWsnIu6LioJPP3WqaF2gRIkStG/fntGjR7NDX5kENSXLELdmzRqWLFlCjx49nCm4XngBpk6FChXcDk1ELnXyJLzzjlPdJ0Hfvn2JiIjgFV0JCmpKliFuwIAB5MqVi06dOjk3D+TNCw0buh2WiCRl/Xro1s2pqJWgUKFCdO7cmYkTJ7Jp0yYXg5MrUbIMYT/99BNz5syhW7duZNu2zbmZZ/Nmt8MSkeRUruxcin3jDThwILG5d+/eZM6cmX79+rkYnFyJSWszdxtj6gB1ChYs2DE2NpasWbNesX9cXNwV+6S03J9ee+01Vq9ezdSpU6kycCDZN2zgm8mTOZst21WN6+Y++XP7vhg3tWN4u56n/T3pF8zHsD+5uV9X2naW7du5s0MHdjdsyG+dOiW2jxs3jkmTJjFq1ChKlizp9bi+iM2X6/q6b6CO4+rVq6+z1t552QJrbZp8lS9f3i5dutSmJKU+nozhD5s3b7YRERH2+eeft3b5cmvB2kGDfDK2W/vk7+37YtzUjuHtep72D+Vj2N/c3K8Ut92mjbUZMli7c2di0+HDh22uXLnsQw89lPpxfRGbj9b1dd9AHcfAWptEztBl2BA1aNAg0qdPzzPduzsVQq69Fp5+2u2wRMQTr7zifG1yQW3YHDly0Lt3bxYuXMjKlStdDE6SomQZgnbu3MmkSZPo2LEj+b//HlatgpdegsyZ3Q5NRDxRpAgsXQq33HJRc5cuXShQoEDi5O0SPJQsQ9Bbb72FMYbnnnvOmZ9y9GhnGi4RCS1//QWxsYlvM2fOTN++fVmxYgWLFy92MTC5lJJliNm7dy9jxozhscceo3DhwpAhg1MwPV06t0MTEW+NHw+dOsEFl107duxI0aJFdXYZZJQsQ8zQoUM5c+YMzz/zDFSrBnPnuh2SiKTW00879xtcUGQ9ffr0vPzyy6xfv55ZF9SSFXcpWYaQAwcOMHLkSJo2bUqJpUvhf/+7rNakiISQzJmdEpVffQXz5iU2t2zZkptuuol+/foRHx/vXnySSMkyhAwfPpxjx47xQvfu8NprcN99ULOm22GJyNVo1w5KlnTuaj93DoDIyEj69u3Lzz//zMKFC10OUEDJMmT8888/vPvuu0RHR3PzkiXOjQEDB2piZ5FQFxXl/L98220XPUrSuHFjChUqxODBg10MTs5TsgwRI0aM4PDhw/R7+mmnVFa9epqrUiRcNGgAkydDzpyJTenSpaNr164sW7aM9evXuxebAEqWIeH48eMMGzaMWrVqUa5aNZg40fmXqIiEl59+gtmzE98+/vjjZMuWjSFDhrgYlICSZUh4//332b9/Py+88IJz2bVePWdyZxEJL336OJO3HzoEOFV9OnTowLRp09i1a5fLwaVtSpZB7tSpU7z55pvce++9VJ01C/r3dzskEfGXAQPg8GHnq5YE3bp1A+Cdd95xKSgBJcugN2HCBP744w9ef+wxZ9LY/fvdDklE/OW226BFC2e+yz17ALj++utp2LAho0eP5tixYy4HmHYpWQaxs2fPMmjQICpUqEDlzz6DjBnhhRfcDktE/OnVV51HSF5+ObGpZ8+eHDlyhE8//dS9uNI4JcsgNmXKFLZv385bTZpgZsyAnj0hXz63wxIRfypWzPl/vUCBxKo+FSpU4N5772XWrFmcPXvW5QDTJiXLIBUfH8/AgQMpW7Ys9y5aBHnyOP8DiUj4GzjQKTxywXPUPXv25K+//mLmzJkuBpZ2RbkdgCRt9uzZbNq0iSlTpmBKlYKdOyF7drfDEpFAsRY++wyKF4eSJalduzaFCxdm8ODBNGnSBKOCJAGVqjNLY0wWY0ykr4MRh7WWAQMGcOONN9KoUSO4/XaoX9/tsEQkkA4dgoYNE5+pjoiIoGHDhqxbt47ly5e7HFza41GyNMZEGGOaG2M+NcbsA34B9hpjNhpj3jLG3OjfMNOWhQsX8t133zHywQeJbNcO/vnH7ZBEJNBy53aeufzoI9i9G4AHH3yQvHnzqkiBCzw9s1wKFAf6AAWstYWttfmAe4CvgUHGmJZ+ijFNsdbSv39/ihUuTI3PP4dvv4UsWdwOS0Tc0KMHxMfD228DkDFjRp566inmzZvH5s2b3Y0tjfE0Wf7HWvuatfZHa23ifDHW2oPW2o+ttY8C0/wTYtqybNkyVq9ezZh778Vs3uw8pBylr5ZF0qSiRaFJExg1yilWAHTu3JkMGTIwbNgwV0NLazxNlsOMMZWv1MFae8YH8aR5AwYMoEi+fFRbtgwqVIDoaLdDEhE3Pfcc5MoFW7YAkC9fPlq1asWECRPYryIlAeNpsvwVGGKM2WGMecMYU86PMaVZX3/9NV988QXj776biD17YNAgTcElktaVKwfbt8NddyU29ejRg5MnTzJixAj34kpjPEqW1trh1tpKwH3AQeADY8wmY0w/Y0xJv0aYhgwYMIDcuXNz1+DBMGwY1KjhdkgiEgwiI+H0aTIl3OhTunRpHnnkEWJiYjhx4oTLwaUNXj06Yq3daa19w1p7O9AciAY2+SWyNOaHH35g/vz5dOvWjSwlS0L37m6HJCLBpEEDyv73v04pPJwiBfv37+fDDz90ObC0watkaYxJZ4ypY4z5CFgIbAEe9Utkaczrr79OsSxZ6LNyJWzc6HY4IhJsWrUi865dMHcuANWqVeP2229n6NChxMfHp7CyXC1Pn7N8wBgzDtgNPA4sAIpba5tYaz/xY3xpwpYtW5gxYwaTbrqJdF9+qbtfReRyjz7KiWuvdabvshZjDD179uSXX35h4cKFbkcX9jw9s/wvsBooba2tY639yFqruWJ8JCYmhuKRkVT+6Sdo2xZKlXI7JBEJNlFR7GrUCL75BlauBKBx48YUKlSIwYMHuxxc+PP0Bp/q1tr3gUPGmJbGmH4Axpgixpi7UlhdriAuLo7x48cztnBhp9bjSy+5HZKIBKk/H3oI8uaFGTMASJcuHd26dWPZsmWsX7/e5ejCm7e1YUcAlYBmCe+PAjE+jSiN+fDDD7n2yBHu2bEDnn4aChVyOyQRCVLxGTM6Vb2GD09s69ixI9myZVMJPD/zNlneba3tDJwEsNYeAtL7PKo0wlpLTEwM+cqWhTffdB4+FhG5kmLFnOevT50CIEeOHHTo0IFp06axa9cul4MLX94myzMJs41YAGPMNYBuw0ql5cuXs2HDBh7r1g3z7LOa2FlEPDNnDhQsmFhgvVu3bgC88847bkYV1rxNlu8As4F8xpgBwErgdZ9HlUbExMTQOksWWp496xRLFhHxxG23ObViEy7HXn/99TRq1IjRo0dz5MgRd2MLU54+OhIFYK39COgFDAT2AvWttTP8F1742rNnD7M+/pg3M2Qgw9ixEJGqqUVFJC1KosB6z549OXLkCGPGjHE1tHDl6Sf0t+d/sdb+Yq2Nsda+Z61V9Z5UGj16NPfFx5P/4EHo3NntcEQk1Dz3HBw9CrGxANx5553ce++9DB8+nLNnz7ocXPjxNFmqmrcPnT59mtGjR/Na/vyQJ4/zL0QREW+UKwcPPuhcij3jTPrUs2dPfv/9d2bOnOlubGHI01Ix1xhjeiS30Fo71EfxpAmzZs0i3Z9/UikiwvnXYcaMbockIqHozTfBWkiXDoDatWtTsmRJBg8eTJMmTZxnt8UnjLU25U7G7AVGkswZprX2FR/H5TfGmDpAnYIFC3aMjY0la9asV+wfFxd3xT4pLU9K165dKbh3L9Ny5+bnV17hZIECXq3vb6nZp1DYvi/GTe0Y3q7naX9P+vnjGA4Fbu6Xm8fw3LlzGTZsGMOGDaNcuXI+jc2bdX3dN1DHcfXq1ddZa++8bIG1NsUXsN6TfqH0Kl++vF26dKlNSUp9PBnjQt9//70F7ODBg71aL5C83adQ2b4vxk3tGN6u52l/N47hUOHmfgX0GD540NpWraydM8daa+3x48dt3rx5bZ06dXwemzfr+rpvoI5jYK1NImfoO8sAi4mJ4db06Wn3qCZrEREfyJYNVq2CgQPBWjJlysRTTz3FvHnz2Lx5s9vRhQ1Pk+X9fo0ijTh06BAfffQRs3LkIJeSpYj4QlQU9OwJX3+dWGC9c+fOZMiQgWHDhrkcXPjwtJD6QX8HkhaMHz+em44fp/j+/fDYY26HIyLhok0bp8D6m28CkC9fPlq1asWECRPYv3+/u7GFCT0JHyDx8fGMGDGCV/LlgyxZlCxFxHcyZ3YmYpg/P3Hy+B49enDy5ElGjBjhcnDhwatkaYxpZIzJlvB7X2PMLGPMHf4JLbwsXryYg1u38tChQ9CqFeTI4XZIIhJOOneGHj0gd24ASpcuzSOPPEJMTAwnTpxwObjQ5+2Z5YvW2qPGmKpATWACziMlkoKYmBiaZs9O5JkzqtgjIr6XJw8MGQLXXpvY1LNnT/bv38+HH37oYmDhwdtkeS7h5yPASGvtHDRFV4q2b9/Op59+Su6uXWHrVihTxu2QRCRcffkljBsHQLVq1bj99tsZOnQo8Zqs4ap4myz3GGNGAY2BBcaYDKkYI80ZOXIkkcbwxBNPQPHibocjIuFszBjo3h0OH8YYQ8+ePfnll19YuHCh25GFNG8TXWNgEVDLWnsYyA1oxuIrOHHiBGPHjuWr/PkpNFRVAUXEz84XWB81CoDGjRtTsGBB3n33XZcDC21eJUtr7XFr7Sxr7a8J7/daaxf7J7TwMHXqVHIfPEiFvXshZ063wxGRcHf77fDAA/D223DqFOnSpaNDhw4sXryYHTt2uB1dyNIlVD+y1vLee+/xYu7c2HTp4PHH3Q5JRNKC55+HP/+EhBt72rVrB8DYsWPdjCqkKVn60TfffMPm9etpcuIEpmFDCLKC6SISpmrUgEceSZxUvkiRItSqVYtx48Zx7ty5FFaWpChZ+lFMTAztM2Qgw4kT0KWL2+GISFphjFOgoG3bxKaOHTvyxx9/8M0337gYWOhKMVkaYx4wxrxvjCmX8F7XEj2wb98+pk+fTo4mTWDwYKhUye2QRCStOXsWFiwAa6lduzYFChTg008/dTuqkOTJmeVTOHe8tjTG1ADK+TWiMDFmzBhOnz5N8z59nCLHmoRVRAJt/HjncuyqVaRLl462bdvy9ddfs3v3brcjCzmeJMv91trD1tpngQeBCn6OKeSdPXuW2NhYRpUowU179rgdjoikVc2bO5V9Egqst2/fnvj4eD744AOXAws9niTLxHN2a21vYKL/wgkP8+bN4+yuXXTYvt25BCIi4obzBdbnzYOff6Z48eLccccdjB07VhV9vJRiskwoaXfhez3ZmoKYmBiey54dEx8PnTq5HY6IpGWdO0OmTDB8OAC1a9dm586dLFmyxOXAQou3s47caYyZbYxZb4z50RjzkzHmR38FF4o2bdrE8i++oGN8POahh6BECbdDEpG0LG9eaNgQvvsOrKVKlSrkzZuX0aNHux1ZSInysv9HODf7/AToHD4JI0aMoFFkJFnj4jS7iIgEhxEjnHl0jSF9+vQ89thjDB8+nL/++ov8+fO7HV1I8PY5y/3W2rnW2u3W2p3nX36JLAQdPXqUCRMmULViRahWDWrVcjskERHImtW5Iz+hIEGHDh04e/Ys48ePdzeuEOJtsnzJGDPGGNPMGNPg/MsvkYWgSZMmcfToUe4YMgSWLk2sniEi4roZM+C664j65x9uuukm7r33XsaMGaMbfTzk7ad5W5znLGsBdRJetX0cU0iy1hITE0PLm27irnLl3A5HRORiJUrAvn3k+9//AKeiz9atW1m2bJm7cYUIb5PlbdbaO621j1lr2ya82vklshCzbNky/vj5Z8Zt24Z5/nm3wxERuVi5cnDLLeRf7EwU9eijj5IzZ07ef/99d+MKEd4my6+NMTf7JZIQFxMTw1OZM5Pu9Gl47DG3wxERuZgx0LIlOTZuhG3byJQpE61atWLWrFn8/fffbkcX9LxNllWB740xm/XoyL92797NnNmzeSZ9eqhc2ZlPTkQk2LRo4fxMmLqrY8eOnD59mkmTJrkYVGjwNlnWAm7EKXt3/vvKOr4OKtSMGjWK/8THk/fwYc0uIiLBq3BhtnbuDHWcj+2yZctSsWJF3n//fay1LgcX3LxKlhc+LqJHRxynT59m9OjRPH399ZA/Pzz6qNshiYgka3fDhhdd/erYsSObNm1i1apVLkYV/Lyt4DPBGJPzgve5jDHjfB5VCFm+fDn79u0jMjYWvv4a0qd3OyQRkStbvRqmTwegSZMmZMuWTTf6pMDby7C3WmsPn39jrT0EpOkv6D755BNuLF6cBx58EIoWdTscEZGUDRvmfGV05gxZsmShRYsWzJgxg8OHD7sdWdDyNllGGGNynX9jjMmN9yXzwsZ3333Hbxs38u0//xCR8IW5iEjQa9kS9u+HhMdIOnbsyIkTJ/joo49cDix4eZsshwBfGWNeM8a8CnwFvOn7sEJDTEwMraKiyPn331CkiNvhiIh4plYtZ57LhLtg77jjDu644w5Gjx6tG32S4e0NPhOBR4G/gP1AA2ttmrzn+MSJE0yfNo1eWbLALbfAffe5HZKIiGfSp4cmTWDOHDhyBHDOLn/88UfWrFnjcnDByevipdban62171lr37XW/uyPoEJBpkyZ2PrRR5T45x9ndhFj3A5JRMRzLVs681xu3AhA8+bNyZw5s270SYYqfV+FfNOnczZLFmjVyu1QRES8U7Ei7N0LlSoBkD17dpo2bcqUKVM4evSoy8EFH4+SpTGmkjE6dbpM9+5s6d7dmf5GRCSUGAMZMoC1cPo04FyKPXbsGFOmTHE5uODj6ZnlY8A6Y8xUY0wbY0wBfwYVMu68k33/+Y/bUYiIpM7Ro1CyJAwfDsDdd99NmTJldCk2CR4lS2vtk9baO4CXgVzAeGPMamPM68aYe40xkf4MUkRE/CBbNsibN/GuWGMMHTt2ZO3atXz//ffuxhZkvL0b9hdr7TBrbS2gBrASaAR844/gRETEz1q1gp9+gh+dOTFatmxJxowZdXZ5iVTf4GOtPWGtXQCst9be6cOYREQkUBo3hqioxJlIcufOTcOGDfnwww85duyYy8EFD1/cDfuKD8YQERE35M0LDz0EkyfDuXOAc6PPkSNHmDFjhsvBBQ+PStVdYc5KA+T3XTgiIhJwzz8P+/Y5d8YC99xzD6VKleL999+nTZs27sYWJDyt65ofqAkcuqTd4JS8ExGRUFWlykVvz9/o8+yzz7Jx40ZuueUWlwILHp5ehp0PZE1iLssdwDK/RSciIoGxZw/07w/HjwPQunVr0qVLpxt9Enj66Eh7a+3KZJY1921IIiIScL/+Ci++CHPnAnDNNdcQHR3NpEmTOHnypMvBuU/l7kREBO69FwoXTnzmEpwbfQ4ePMisWbNcDCw4KFmKiAhERECLFrBokXOzD1CjRg1uuOEGXYpFyVJERM5r2dJ5fGTaNAAiIiLo0KEDy5YtY8uWLS4H5y4lSxERcdxyC9x1l3OzT4I2bdoQGRnJmDFjXAzMfZ4+OiIiImnBV19B5L/lvq+99lrq1KnD+PHjeeCBB1wMzF06sxQRkX+dT5RxcYlNjz/+OPv37+err9LuY/VKliIicrFnn4Xbbkus6PPggw9SpEgR5s+f73Jg7lGyFBGRi912G2zbBqtWARAZGUm7du1Yu3Yt27dvdzk4dyhZiojIxaKjIXPmxJlIANq1a0dERARjx451MTD3KFmKiMjFsmZ1Eub06XDqFACFCxemQoUKfPDBB5w9e9blAANPyVJERC7XsiUcOgQLFyY21a5dmz/++IPPPvvMxcDcoWQpIiKX+89/IDYWqlZNbKpYsSJ58uTho48+cjEwdyhZiojI5aKi4IknnMmhE5uiaNSoEXPmzCHugkdL0gIlSxERSdq5czBmDFzwyEiLFi04ceIEc+bMcTGwwFOyFBGRpEVEwJAh8OabiU2VK1emSJEiTJ482cXAAk/JUkREkmYMtGoFK1bAjh2AU1y9WbNmLFq0iP3797sbXwApWYqISPKaN3d+XnBTT4sWLTh37hwzZsxwKajAC8lkaYy5wRgz1hgz84K20saYWGPMTGNMJzfjExEJG0WLOhNDf/hhYvm7smXLUqZMmTR1KTbgydIYM84Ys88Ys+GS9lrGmM3GmK3GmN5XGsNau81a2/6Stk3W2ieBxsCdvo9cRCSNatkSsmcn6ujRxKbmzZuzatUqdiRcng13bpxZjgdqXdhgjIkEYoCHgJuBZsaYm40xZY0x8y955UtuYGNMXWAl8IX/whcRSWM6dIBvvuFs9uyJTc2aNQNgypQpbkUVUAFPltba5cDBS5rvArYmnDGeBqYC9ay1P1lra1/y2neFsedaaysDLfy3ByIiaYwxAETGxUFCqbuiRYtSpUqVNHMp1tiEa9AB3agxRYH51toyCe8bArWstR0S3rcC7rbWdklm/TzAAOABYIy1dqAxphrQAMgA/GitjUlivceBxwHy589ffsyYMWTNmvWKscbFxV2xT0rLQ5Hb++Sv7fti3NSO4e16nvb3pF9aPIbB3f0Kx2M42y+/UK5rVza+9hoH774bgE8++YThw4czduxYbrjhhlRvJ5iO4+rVq6+z1l7+VZ61NuAvoCiw4YL3jXCS3vn3rYB3/RlD+fLl7dKlS21KUurjyRihxu198tf2fTFuasfwdj1P++sYTp6b+xWOx7A9dcqezp7d2ubNE5v27dtnIyMj7fPPP39V2wmm4xhYa5PIGcFyN+xuoPAF7wsBf7gUi4iIXCp9evZVqwazZ0PCjT7XXHMNDz74IFOmTCE+Pt7d+PwsWJLlGuBGY0wxY0x6oCkw1+WYRETkAn/95z9w4gR88kliW4sWLfj999/56quv3AssANx4dGQKsBooZYzZbYxpb609C3QBFgGbgOnW2o2Bjk1ERJJ35JZboEgRmDo1sa1evXpkypQp7G/0iQr0Bq21zZJpXwAsCHA4IiLiqYgIeP99KPzvt2ZZs2alXr16TJ8+neHDh5MuXToXA/SfYLkMKyIioeDBB6F06YuamjdvzoEDB1i8eLFLQfmfkqWIiHhn2TIYMCDxbc2aNcmdO3dYX4pVshQREe8sXQovvgh79wKQPn16GjVqxCeffMKxY8dcDs4/lCxFRMQ7TZo4RdVnJs5lQfPmzTl+/HjYTgrtSgUfNxlj6gB1ChYs2DE2NjZoqkYEE7f3KRyrn6iCT+Cpgo/vx7hw3Tvbt+dcpkx89957AMTHx9OsWTNuuOEGBg4cqAo+4fJSBZ/kub1P4Vj9RBV8Ak8VfHw/xkXrvv66tWDtjh2JTb169bJRUVF2//79quAjIiJCkyZw/fWwfXtiU/PmzTl79iwzL7g8Gy6ULEVExHs33OAkymrVEptuvfVWbr75Zj766CP34vITJUsREUkdY5wpuxLugDXG0KJFC1auXMmff/7pcnC+pWQpIiKpc+yYU81nyJDEpvOTQi9dutStqPxCyVJERFInSxYoVQqmTHEeJQGKFStGpUqV+Pzzz10OzreULEVEJPWaNoVffoGffkpsat68Odu2beOnC9pCnZKliIik3qOPQmSkc3aZoHHjxkRERDDlgrZQp2QpIiKpd8018MADzrRdCZdi8+XLx5133snkyZOxYVL4RslSRESuzosvwtixickS4P7772fnzp1hMym0kqWIiFydypWhRg1nvssEVatWDatJoVUbNkjqEQYTt/cpHOtqqjZs4Kk2rO/HuNK6mX//nXyff86Oxx6DyEji4uIYOnQo69evZ+bMmURFRV1VTKoNq9qwQcftfQrHupqqDRt4qg3r+zGuuO7UqU6t2IQ+S5cutXPmzLGAXbBgwVXHpNqwIiIS+mrXdp67nDo1salWrVrkypUrLMrfKVmKiMjVy5IF6tZ15rg8cwYIr0mhlSxFRMQ3mjaFAwfgiy8Sm5o3b86xY8eYN2+ei4FdPSVLERHxjZo1oUgR2LUrsemee+6hYMGCIX9XrJKliIj4RoYMzrRdHTsmNkVERNCsWTMWLlzIgQMHXAzu6ihZioiI70REgLVEnDyZ2NSiRYuQnxRayVJERHzr/vu56Y03Et/edtttlC5dOqQvxSpZioiIb5UuTZ7VqyEuDnAmhW7evDnLly9n1wXfZ4YSJUsREfGtpk2JPHUK5s5NbDo/KXSozkSiZCkiIr5VpQonr7nmogIFxYsXp2LFiiF7KVbJUkREfCsigv3Vq8Nnn8GhQ4nNzZs354cffmDjxo0uBpc6KqQeJMV7g4nb+xSORahVSD3wVEjd92N4s67dtIn8e/bw9z33EJ8hAwAHDx6kUaNGNG/enPbt23s1rgqpq5B60HF7n8KxCLUKqQeeCqn7fgxv1k2u74MPPmiLFStm4+PjvRpXhdRFRCQ87d8Pb70Ff/2V2NSiRQu2b9/O119/7WJg3lOyFBER/9i3D3r1coqrJ6hfvz4ZM2YMuRt9lCxFRMQ/brkFypS56K7Y7NmzU6dOHaZNm8aZhNlJQoGSpYiI+E/TprBy5UXF1Vu0aMH+/fv54oLZSYKdkqWIiPhPkybOz+nTE5tq1apFzpw5Q+pSrJKliIj4T4kScPfd8PvviU0ZMmSgYcOGzJ49mxMnTrgYnOeULEVExL9WrIDhwy9qatiwIXFxcXz++ecuBeUdJUsREfGvdOmcn6dPJzZVr16dHDlyMHv2bJeC8o6SpYiI+F/PnlChQuLb9OnTU7t2bebOncu5c+dcDMwzSpYiIuJ/xYvDjz/Chg2JTQ0aNODAgQP8+OOPLgbmGSVLERHxv4YNISLiomcua9asScaMGVmxYoWLgXlGyVJERPwvXz64/34nWSZM4JElSxZq1arFypUrsUE+qYeSpYiIBEbTpvDbb7BuXWJTdHQ0+/fvZ+3atS4GljJN0RUk08IEE7f3KRynN9IUXYGnKbp8P4Y36ybVN+roUQosXMhfDz7ImZw5AThy5AgNGjSgSZMmdOzYMdXb1hRdmqIr4Nzep3Cc3khTdAWepujy/Ri+mKIrKeXLl7clS5a8aNoub8fTFF0iIhI+Tp6ESZOcO2MTVK1alS1btrBp0yYXA7syJUsREQmcs2fhiSdg1KjEpqpVqwIEdYECJUsREQmcrFmhTh2YMcNJnEDevHmpWLEis2bNcjm45ClZiohIYDVtCvv3w9KliU0NGjRg/fr17Ny508XAkqdkKSIigfXQQ5A9+0UFCqKjo4HgvRSrZCkiIoGVMSNER8OvvyYWKChRogRly5ZVshQREUk0ahQsXw7GJDZFR0ezYsUK9u3b52JgSVOyFBGRwMuQwfl5QWGcBg0aYK1l7ty5LgWVPCVLERFxx/DhcNNNEB8PwK233kqxYsWC8q5YJUsREXFHvnywZQvZf/4ZAGMMDRo04IsvvuDIkSMuB3cxJUsREXHHww9DunTkXbkysSk6OprTp0+zYMECFwO7nJKliIi4I0cOqFGDa1asSPzuslKlSuTPnz/oLsUqWYqIiHuio8n0xx+wYQMAERER1K9fnwULFnDy5EmXg/uXkqWIiLinXj1+b9bMOctMEB0dzbFjx1iyZImLgV1MyVJERNxToADbHn8cihRJbKpevTo5cuQIqgIFSpYiIuIqc/YsLFkCe/cCkD59emrXrs3cuXM5m1Bs3W3GXvBAaFpgjKkD1ClYsGDH2NhYzTKfBLf3KRxnmfd2PU/7e9IvLR7D4O5+heMx7O263vSN37qVGh07svWpp9jdqBEAy5cv56WXXmLo0KHcfvvtATuOq1evvs5ae+dlC5KaETotvMqXL69Z5pPh9j6F4yzz3q7naX8dw8lzc7/C8Rj2dl2v+956q7VVqya2xcXF2YwZM9ouXbp4NJ6v/ubAWptEztBlWBERcV90NKxaBX/9BUCWLFmoVasWs2fPJj6hwo+blCxFRMR90dHOs5YX1IWNjo5mz549rF271sXAHEqWIiLivltvhWLF4IsvEptq165NVFRUUNwVq2QpIiLuMwa+/BI++iixKXfu3FSrVo1Zs2ZhXb4ZVclSRESCQ9GiEBl5UVODBg3YsmULO3fudCemBEqWIiISPF57DXr0SHxbr149AFZeUGzdDUqWIiISPHbtgvffh4S6sNdddx0VK1Zk+fLlroalZCkiIsEjOhri4i660adBgwb8+uuvrl6KVbIUEZHgUaMGZMsGF9wBGx0dDeDqXbFKliIiEjwyZIBHHnGetzx3DoASJUpwww03KFmKiIgkatECHnoIjhxJbKpatSorVqxg3759roSkZCkiIsGldm2YMAFy5Upsuueee7DWMveCCj+BpGQpIiLBx1rYuNH5CRQvXpxixYoxa9YsV8JRshQRkeAzZQqUKQPffw+AMYYGDRrwxRdf8M8//wQ8HCVLEREJPg88ABERl90Ve/r0aRYsWBDwcJQsRUQk+FxzDVStelGyrFSpEgUKFHDlrlglSxERCU7R0bBhA2zdCkBERAT16tVjwYIFnEyo8BMoSpYiIhKc6td3fl5yKfbYsWMsWbIkoKEoWYqISHAqWhQWLoQnn0xsql69Ojly5Aj4pVjj9hxhgWaMqQPUKViwYMfY2FiyZs16xf5xcXFX7JPS8lDk9j75a/u+GDe1Y3i7nqf9PemXFo9hcHe/wvEY9nZdX/e9sM+AAQP49ttvmTVrFpEJU3r56m9evXr1ddbaOy9bYK1Nk6/y5cvbpUuX2pSk1MeTMUKN2/vkr+37YtzUjuHtep721zGcPDf3KxyPYW/X9Vnf06etHTjQ/vTaa4lNH3/8sQXsl19+martXQmw1iaRM3QZVkREgldUFIwdy3Vz5iQ21axZk4wZMwa0QIGSpYiIBC9jIDqanN99B4cPA5AlSxZq1arF7NmziY+PD0gYSpYiIhLcoqOJOHcOPv30gqZo9uzZw9q1awMSgpKliIgEt7vv5lSePBc9QlK7dm2ioqICdleskqWIiAS3iAj233efU1Q94QmO3LlzU61aNWbNmoUNwFMdSpYiIhL0tnbpAh9/7HyHmaBBgwZs2bKFTZs2+X37SpYiIhL8zifJY8cSm+rVqwcQkLtilSxFRCQ0vPaaU9Xn7FkArrvuOipVqhSQ7y2VLEVEJDTccgv8/TcsX57YFB0dzfr16/nzzz/9umklSxERCQ01a0LGjJcVVgdYuXKlXzetZCkiIqEhSxYnYX7ySeJdsSVKlKBs2bKsW7fOr5uO8uvoIiIivhQdDXPmwNq1UKECAPPmzWNrwpyX/qIzSxERCR116sDw4c6NPgmuv/76xNlH/EVnliIiEjpy54auXQO+WZ1ZiohIaDl6FMaPhx07ArZJJUsREQkt//wDbdvClCkB26SSpYiIhJZChZybewJURB2ULEVEJBRFR8OaNbB7d0A2p2QpIiKhJ6EYAZ98EpDNKVmKiEjouekmKF0afvwxIJvToyMiIhKavv4asmcPyKZ0ZikiIqEpQIkSlCxFRCSUde0KLVv6fTNKliIiErqshY8/JuLECb9uRslSRERCV3Q0nDxJ7jVr/LoZJUsREQld994LlStjEqbs8hdj/byBYGOMqQPUKViwYMfY2FiyZs16xf5xcXFX7JPS8lDk9j75a/u+GDe1Y3i7nqf9PemXFo9hcHe/wvEY9nZdX/cN1HFcvXr1ddbaOy9bYK1Nk6/y5cvbpUuX2pSk1MeTMUKN2/vkr+37YtzUjuHtep721zGcPDf3KxyPYW/X9XXfQB3HwFqbRM7QZVgREZEUKFmKiIikQMlSREQkBUqWIiIiKVCyFBERSYGSpYiISAqULEVERFKgZCkiIpICJUsREZEUKFmKiIikQMlSREQkBUqWIiIiKVCyFBERSYGSpYiISAqULEVERFKQ5iZ/Ps8Ysx84DPyTQtccKfTJC/zto7CCRUr7HKrb98W4qR3D2/U87e9Jv7R4DIO7x3E4HsPeruvrvoE6jq+31l5zWWtSk1ymlRcw+mr7kMxEoaH88uTvEorb98W4qR3D2/U87a9j2L//vYNt224ew96u6+u+bh/Haf0y7Dwf9Qk3bu+zv7bvi3FTO4a363naX8dw8tzc73A8hr1d19d9XT2O0+xlWF8xxqy11t7pdhwiqaVjWMKBv4/jtH5m6Quj3Q5A5CrpGJZw4NfjWGeWIiIiKdCZpYiISAqULEVERFKgZCkiIpICJUsREZEUKFn6kTHmBmPMWGPMTLdjEfGUMSaLMWaCMeZ9Y0wLt+MRSQ1ff/4qWSbDGDPOGLPPGLPhkvZaxpjNxpitxpjeVxrDWrvNWtvev5GKpMzL47kBMNNa2xGoG/BgRZLhzXHs689fJcvkjQdqXdhgjIkEYoCHgJuBZsaYm40xZY0x8y955Qt8yCLJGo+HxzNQCNiV0O1cAGMUScl4PD+OfSrK1wOGC2vtcmNM0Uua7wK2Wmu3ARhjpgL1rLUDgdoBDlHEY94cz8BunIT5PfoHtQQRL4/jn325bf2P4J2C/PsvbnA+VAom19kYk8cYEwvcbozp4+/gRLyU3PE8C3jUGDOStFtXVkJHksexrz9/dWbpHZNEW7IlkKy1B4An/ReOyFVJ8ni21h4D2gY6GJFUSu449unnr84svbMbKHzB+0LAHy7FInK1dDxLOAjIcaxk6Z01wI3GmGLGmPRAU2CuyzGJpJaOZwkHATmOlSyTYYyZAqwGShljdhtj2ltrzwJdgEXAJmC6tXajm3GKeELHs4QDN49jzToiIiKSAp1ZioiIpEDJUkREJAVKliIiIilQshQREUmBkqWIiEgKlCxFRERSoGQpEuKMMeeMMd9f8Lps6jhjzLKEKYzqGmNiEvr9bIw5ccF6DS9ZJ4sx5oAxJscl7Z8YYxobY5okTIk039/7KOI21YYVCX0nrLXlPOjXwlq7loTqJgmzN8xPbl1r7TFjzGKgPjAhYZ0cQFWgubX2uDHmL+DZq90BkWCnM0sROX8WOc4Ys8YY850xpl7Coik45cPOiwY+s9YeD3yUIu5RshQJfZkuuQzbJBVjvAB8aa2tAFQH3jLGZAE+A8obY/Ik9GuKk0BF0hRdhhUJfZ5ehr2SB4G6xpjzl1QzAkWstZuMMXOBhsaYj4FywOKr3JZIyFGyFBFw5gR81Fq7OYllU4C+CX3mWGvPBDQykSCgy7AiAs6MDU8bYwyAMeb2C5YtBW4EOqNLsJJGKVmKhL5Lv7MclIoxXgPSAT8aYzYkvAfAWhsPfAzkAZb7JGKREKPLsCIhzlobmcr1dgBlEn4/ATxxhb7dgG6p2Y5IONCZpUjacBAYb4yp66sBE+66HQEc8tWYIsFKkz+LiIikQGeWIiIiKVCyFBERSYGSpYiISAqULEVERFKgZCkiIpKC/wPgd6AfS2GMuwAAAABJRU5ErkJggg==\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "plt.figure(figsize=(7, 7))\n", - "plt.loglog(E, sed, 'k-')\n", - "plt.loglog(E, A*sed, 'r--')\n", - "plt.xlabel(\"E [TeV]\")\n", - "plt.ylabel(\"1/ cm$^2$ s TeV)\")\n", - "plt.grid(which='both')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Function inside agnpy" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This notebook demostrate how to use Absorption class to calcuare tau in BLR, which 26 shells following the Finke paper from 2016" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": {}, - "outputs": [], - "source": [ - "def get_absorbtion_tau(tau):\n", - " \"Function return absorption\"\n", - " return np.exp(-tau)" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": {}, - "outputs": [], - "source": [ - "# Set BLR parameters\n", - "L_disk=6e45*u.Unit(\"erg s-1\")\n", - "R_line = 2.45*1e17*u.cm\n", - "r_blob_bh = 1e18*u.cm\n", - "z = 1\n", - "xi_line = 0.1\n", - "\n", - "E = np.logspace(-1, 1, 20)*u.TeV\n", - "nu = E.to(\"Hz\", equivalencies=u.spectral())\n", - "\n", - "blr_one_line = SphericalShellBLR(L_disk, xi_line, \"Hbeta\", R_line)\n", - "abs_blr_one_line = Absorption(blr, r=r_blob_bh, z=z)\n" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/home/pgliwny/anaconda3/envs/agn/lib/python3.8/site-packages/astropy/units/quantity.py:611: RuntimeWarning: invalid value encountered in sqrt\n", - " result = super().__array_ufunc__(function, method, *arrays, **kwargs)\n" - ] - } - ], - "source": [ - "# Set BLR parameters\n", - "L_disk=6e45*u.Unit(\"erg s-1\")\n", - "R_line = 2.45*1e17*u.cm\n", - "r_blob_bh = 1e18*u.cm\n", - "z = 1\n", - "xi_line = 0.1\n", - "\n", - "E = np.logspace(-1, 1, 20)*u.TeV\n", - "nu = E.to(\"Hz\", equivalencies=u.spectral())\n", - "\n", - "blr = SphericalShellBLR(L_disk, xi_line, \"Hbeta\", R_line)\n", - "abs_blr = Absorption(blr, r=r_blob_bh, z=z)\n", - "tau_all_lines = abs_blr.tau_blr_all_lines_cubepy(nu)\n", - "A_blr = get_absorbtion_tau(np.array(tau_all_lines))" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 9, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAD8CAYAAAC7IukgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxiklEQVR4nO3dd3hVZbrw/++dTiCFBBIgCSQUITGEZgI6AziCM4DYG7ajgFiOjmXOmfPTOXNef+87nvGM75wZdcaxotgGVNSxDGIbFQu9hyYdEkoCgVBT9/3+sYMnbBOyd7KTtcv9ua59hb3Ws591b69l7jx1iapijDHG+FOE0wEYY4wJPZZcjDHG+J0lF2OMMX5nycUYY4zfWXIxxhjjd5ZcjDHG+F2U0wEEgm7duml2drbTYRhjTFBZvnz5AVXt3tQ5Sy5AdnY2y5YtczoMY4wJKiKys7lz1i1mjDHG7yy5GGOM8TtLLsYYY/zOkosxxhi/s+RijDHG7yy5GGOM8TubitwGB49Vs2znIeKiI4mNijjjz8gIcTpcY4zpMJZc2mDD3qPc/spyr8pGRcgPkk5sdCRJnaJ47Nph9EiKa+dojTGm41hyaYOhvZP54Oc/prqunupaF1UeP6vrXFTV1jf7s6q2nk83lPHe6lJuG9PP6a9jjDF+Y8mlDbrERpGfkdSmOiY+/hWfri+z5GKMCSk2oO+wC3PTWLazgkPHa5wOxRhj/MaSi8PG5abjUvh8U5nToRhjjN9YcnHY4Iwk0hJi+XTDfqdDMcYYv7Hk4rCICGFcbjpfbiqnuq7e6XCMMcYvLLkEgAvz0jheU8+ibRVOh2KMMX5hySUAnNevG52iI/l0vXWNGWNCgyWXABAXHcnoAd34dMN+VNXpcIwxps0suQSI8Xnp7K2sYt2eI06HYowxbWbJJUBcMCgNEWzWmDEmJIRschGRviIyU0TmOh2LN7p1iWV4766WXIwxIcGr5CIi94pIsYisE5H7mimTLCJzRWSjiGwQkXNbG5SIvCAiZSJS7HF8gohsEpEtIvLAmepQ1W2qOr21MThhfG46xaVH2Ft50ulQjDGmTVpMLiKSD8wAioAhwGQRGdBE0ceB+ao6qKHcBo960kQkweNY/2YuOwuY4FE2EngSmAjkAdeJSF7DucEi8oHHK62l7xZoLsxzh/zpBlutb4wJbt60XHKBRap6QlXrgC+ByxsXEJFEYAwwE0BVa1T1sEc9Y4F3RSSu4TMzgCeauqCqLgA8F30UAVsaWiQ1wBzg0obya1V1sscr6H5D9+vehezUeJuSbIwJet4kl2JgjIikikg8MAnI8ijTFygHXhSRlSLyvIh0blxAVd8E5gNzROQGYBpwjQ+xZgC7G70vaTjWpIZ4nwaGiciDzZS5WESerays9CGM9iMijM9NZ+HWgxyrrnM6HGOMabUWk4uqbgB+B3yCOzmsBjx/80UBw4GnVHUYcBz4wZiIqj4KVAFPAZeo6jEfYm3qUY7NLgpR1YOqeoeq9lPVR5op876q3paU1LZt8/1pfF46NfUuvvqu3OlQjDGm1bwa0FfVmao6XFXH4O6u2uxRpAQoUdXFDe/n4k42pxGR0UA+8A7wkI+xlnB6iykT2ONjHQHvnD5dSeoUzSc2a8wYE8S8nS2W1vCzN3AFMLvxeVXdB+wWkYENh8YB6z3qGAY8h3ucZCqQIiIP+xDrUmCAiOSISAwwBXjPh88HhajICC4YlMbnG8uod9lqfWNMcPJ2nctbIrIeeB+4S1UPAYjIPBHp1VDm58BrIrIGGAr81qOOeOBqVd2qqi7gZmBnUxcTkdnAQmCgiJSIyPSGyQR3Ax/hnon2hqqu8/aLBpPxuekcOlHLil2HnA7FGGNaxavHHKvq6GaOT2r071XAOWeo4xuP97W4WzJNlb2umePzgHktRxzcxpzVjehI4dP1+ynMTnE6HGOM8VnIrtAPZglx0Yzqm2rjLsaYoGXJJUCNz01nW/lxtpb7MqHOGGMCgyWXADUu171a/zNrvRhjgpAllwCV2TWe3J6JfLo+6DYaMMYYSy6B7MLcNJbtrKDieI3ToRhjjE8suQSw8XnpuBQ+32itF2NMcLHkEsDyeyWRnhhrz3gxxgQdSy4BLCJCGJebzpfflVNVW+90OMYY4zVLLgHuwtx0TtTUs2jbQadDMcYYr1lyCXDn9kulU3SkdY0ZY4KKJZcAFxcdyZizuvHp+jJUbSNLY0xwsOQSBMbnprPvSBXr9hxxOhRjjPGKJZcgcMGgNETgE3v8sTEmSFhyCQKpXWIZ0bsrn2205GKMCQ6WXILE+Lx0ikuPsLfypNOhGGNMiyy5BInxuekAfLrBVusbYwKfJZcg0a97Z3K6deZTG3cxxgQBSy5BQkQYn5vGwq0HOVZd53Q4xhhzRpZcgsj43HRq6l189V2506EYY8wZWXIJIiP6dCWpU7Q9/tgYE/AsuQSRqMgILhiUxucby6irdzkdjjHGNCtkk4uI9BWRmSIy1+lY/Gl8bjqHTtSyYtdhp0MxxphmeZVcROReESkWkXUict8ZykWKyEoR+aAtQYnICyJSJiLFHscniMgmEdkiIg+cqQ5V3aaq09sSRyAac1Y3oiPFNrI0xgS0FpOLiOQDM4AiYAgwWUQGNFP8XmBDM/WkiUiCx7H+zdQzC5jgUTYSeBKYCOQB14lIXsO5wSLygccrraXvFowS4qIZ1TfVpiQbYwKaNy2XXGCRqp5Q1TrgS+Byz0IikglcBDzfTD1jgXdFJK6h/AzgiaYKquoCoMLjcBGwpaFFUgPMAS5tKL9WVSd7vFpcbSgiF4vIs5WVlS0VDSgX5qWz7cBxtpYfczoUY4xpkjfJpRgYIyKpIhIPTAKymij3GPBvQJMjzar6JjAfmCMiNwDTgGt8iDUD2N3ofUnDsSY1xPs0MExEHmwmpvdV9bakpCQfwnDeuFOr9a31YowJUC0mF1XdAPwO+AR3clgNnLaKT0QmA2WquryFuh4FqoCngEtU1Zc/vaWpKs9wrYOqeoeq9lPVR3y4TsDLSO5EXs9EG3cxxgQsrwb0VXWmqg5X1TG4u6s2exT5EXCJiOzA3V11gYi86lmPiIwG8oF3gId8jLWE01tMmcAeH+sIGePz0lm+8xAVx2ucDsUYY37A29liaQ0/ewNXALMbn1fVB1U1U1WzgSnAP1T1Ro86hgHP4R4nmQqkiMjDPsS6FBggIjkiEtNwnfd8+HxIuTA3HZfC5xttI0tjTODxdp3LWyKyHngfuEtVDwGIyDwR6eVlHfHA1aq6VVVdwM3AzqYKishsYCEwUERKRGR6w2SCu4GPcM9Ie0NV13l57ZCTn5FIemKsdY0ZYwJSlDeFVHV0M8cnNXHsC+CLJo5/4/G+FndLpql6r2vm+DxgXosBhwH3RpbpvLOylKraeuKiI50OyRhjvheyK/TDwfjcdE7U1LNsxyGnQzHGmNNYcgliw/t0BWDVbksuxpjAYskliCV1iqZvt86sLgmuRaDGmNBnySXIFWQmsabksNNhGGPMaSy5BLmCzGT2H6lm/5Eqp0MxxpjvWXIJckOy3FvXrN592NlAjDGmEUsuQS6vZxKREcIaG3cxxgQQSy5BrlNMJGelJ7Daxl2MMQHEkksIGJKZxNrSSlSb3cfTGGM6lCWXEFCQmczhE7XsrjjpdCjGGANYcgkJBZkNg/rWNWaMCRCWXELAwB4JxERF2HoXY0zAsOQSAqIjI8jrmWgr9Y0xAcOSS4gYkplEcWkl9S4b1DfGOM+SS4goyEzmRE09W8t9eXK0Mca0D0suIcJW6htjAokllxDRt1sXusRG2Up9Y0xAsOQSIiIihPyMRJsxZowJCJZcQsiQzGQ27D1KTZ3L6VCMMWHOkksIKchMpqbexcZ9R5wOxRgT5kI6uYhIXxGZKSJznY6lI/zPSn0bdzHGOMvr5CIi94pIsYisE5H7mjifJSKfi8iGhjL3tjYoEXlBRMpEpLiJcxNEZJOIbBGRB85Uj6puU9XprY0j2GR27URK5xjW2IwxY4zDvEouIpIPzACKgCHAZBEZ4FGsDvgXVc0FRgF3iUieRz1pIpLgcax/E5ecBUxoIo5I4ElgIpAHXCcieSIyWEQ+8HilefPdQomINDz22FouxhhnedtyyQUWqeoJVa0DvgQub1xAVfeq6oqGfx8FNgAZHvWMBd4VkTgAEZkBPOF5MVVdAFQ0EUcRsKWhRVIDzAEuVdW1qjrZ41Xm5XcLKQWZyWwuO8qJmjqnQzHGhDFvk0sxMEZEUkUkHpgEZDVXWESygWHA4sbHVfVNYD4wR0RuAKYB1/gQbwawu9H7En6YwBrHkSoiTwPDROTBJs5fLCLPVlaGzl/6BRlJuBTW7bFBfWOMc7xKLqq6Afgd8Anu5LAadzfYD4hIF+At4D5V/cFvOFV9FKgCngIuUVVf9iuRpsI7Q9wHVfUOVe2nqo80cf59Vb0tKSnJhxACW4Gt1DfGBACvB/RVdaaqDlfVMbi7rDZ7lhGRaNyJ5TVVfbupekRkNJAPvAM85GO8JZzeYsoE9vhYR0hLS4ijZ1KcjbsYYxzly2yxtIafvYErgNke5wWYCWxQ1T80U8cw4DngUmAqkCIiD/sQ71JggIjkiEgMMAV4z4fPhwX3oP5hp8MwxoQxX9a5vCUi64H3gbtU9RCAiMwTkV7Aj4CbgAtEZFXDa5JHHfHA1aq6VVVdwM3ATs8LichsYCEwUERKRGQ6QMNkgruBj3BPGHhDVdf58oXDQUFmMjsOnqDyRK3ToRhjwlSUtwVVdXQzx08lkD00PSbSuOw3Hu9rcbdkPMtdd4Y65gHzWoo3nA3JTAZgTelhRg/o7mwwxpiwFNIr9MPV4IaV+jbuYoxxiiWXEJTUKZqcbp1txpgxxjGWXEKUrdQ3xjjJkkuIKshMZt+RKsqOVDkdijEmDFlyCVFDbIdkY4yDLLmEqLN7JREZIbbexRjjCEsuIapTTCQD0rpYy8UY4whLLiHs1Ep91Wa3XzPGmHZhySWEFWQmc/hELSWHTjodijEmzFhyCWGnVuqvtnEXY9pFvUutZ6AZllxC2MAeCcRERth6F2PawYpdhyj6z0957qttTocSkCy5hLCYqAhyeyXaSn1j/OybLQe48fnFHDxew0fr9jsdTkCy5BLihmQmUVxaSb3Lmu7G+MPH6/Yx9cWlZHWN58rhmawpOUxVbb3TYQUcSy4hriAzmeM19Wwr9+WBn8aYpvxtZSl3vraC3F6JvH77KC4q6EFtvbJy12GnQws4llxCnK3UN8Y/Xlm4g/vfWMXInBReu3UkyfExjOiTgggs2V7hdHgBx5JLiOvbvQudYyJtpb4xbfCXL7bwH++uY9ygNF64pZAuse5HYSV1imZQj0SW7DjocISBx5JLiIuMEPIzkqzlYkwrqCr/9eFGHp2/iUuH9uKpG0cQFx15WpmROSms2HmY2nqXQ1EGJksuYWBIVjIb9hyhps5ufmO85XIpv/5bMU9/uZUbRvbmj9cMJTryh78yi3JSOFlbT3Gp/QHXmCWXMFCQmURNvYtN+446HYoxQaG23sX9b6zitcW7uGNsPx6+LJ+IiKaf4l6YnQLYuIsnSy5hwFbqG+O9qtp67nx1Oe+u2sO/TRjIAxMHIdJ0YgHonhBL326dWbrDkktjllzCQGbXTnSNj7ZBfWNacKy6jqkvLuWzjWX85rJ8/vn8/l59rignhSXbK3DZerLvWXIJAyLC4Mxk2wbGmDM4fKKGG59fzJIdFfzhmiHcNKqP158tyknhSFUdm/Zb1/MpIZ1cRKSviMwUkblOx+K0IZlJfLf/KCdrbCWxMZ7KjlRx7TOLWL/3CE/fOILLh2X69Hkbd/khvyYXEblXRIpFZJ2I3NeGel4QkTIRKW7i3AQR2SQiW0TkgTPVo6rbVHV6a+MIJQWZybgU1u2x1osxje2uOMHVzyxk96ETvHhLIRfmpftcR2bXTvRKimOJjbt8z2/JRUTygRlAETAEmCwiAzzKpIlIgsexpjo1ZwETmrhGJPAkMBHIA64TkTwRGSwiH3i80vzyxUKErdQ35oe2lB3j6qcXcvhELa/dOpIf9e/WqnpE5PtxF9uC382fLZdcYJGqnlDVOuBL4HKPMmOBd0UkDkBEZgBPeFakqguApv4EKAK2NLRIaoA5wKWqulZVJ3u8yloKWEQuFpFnKytD/xduWmIcPRLjbFDfmAZrSg5zzTMLqXMpc24bxbDeXdtUX2FOCuVHq9lx8ISfIgxu/kwuxcAYEUkVkXhgEpDVuICqvgnMB+aIyA3ANOAaH66RAexu9L6k4ViTGmJ5GhgmIg96nlfV91X1tqSkJB9CCF7uxx6HfiI1piXfbDnAdc8uIj4mkjfvOJfcnoltrnNkjnvcZamNuwB+TC6qugH4HfAJ7gSyGqhrotyjQBXwFHCJqvqyXW9Tk82bbYOq6kFVvUNV+6nqIz5cJyQNyUpm+4HjVJ6sdToUYxwzb+1e95b5KfG8ded55HTr7Jd6+3XvQkrnGBZbcgH8PKCvqjNVdbiqjsHdrbXZs4yIjAbygXeAh3y8RAmnt4YygT2tDDfsFDSMu6y11osJU68u2sldf11BQWYSr992LumJcX6rW0QozO5qm1g28PdssbSGn72BK4DZHueHAc8BlwJTgRQRediHSywFBohIjojEAFOA9/wRezgoyEgGbKW+CT+qyhOfbebXfyvmgoFpvDJ9JEnx0X6/TlFOKrsrTrK38qTf6w42/l7n8paIrAfeB+5S1UMe5+OBq1V1q6q6gJuBnZ6ViMhsYCEwUERKRGQ6QMNEgbuBj4ANwBuqus7P3yFkJcVHk50ab4P6Jqy4XMr/fn89f/jkO64YnsHTN42gU0xkyx9shVPjLrbeBaL8WZmqjm7h/Dce72txt2Q8y113hjrmAfNaG2O4K8hMtj2QTNioqXPxr2+u5r3Ve5gxOocHJ+Y2uwGlP+T2TKRLbBRLtldw6dBm5xqFhZBeoW9+qCAzib2VVZQdrXI6FGPa1YmaOm59eRnvrd7DAxMH8e8X5bVrYgH385NG9Olqf8BhySXsFDTskLxmtw3qm9B16HgN1z+3mK83l/PolQXcMbZfh127KCeF7/Yfo+J4TYddMxBZcgkz+RmJRAg27mJC1p7DJ7n6mYXf7xN2TWFWyx/yo6JT613CvPViySXMxMdEMSAtwbaBMSFpS9kxrnrqW/ZXVvHytCJ+enaPDo+hIDOJmKiIsF9MacklDBVkJrG2tNL2QDIhZdXuw1z99LfU1Ctzbh/FqL6pjsQRGxXJsKzksN/E0pJLGCrISqbieA0lh2wuvgkNX20u5/rnFpEQF81bd57L2b2c3dKpKCeF4tJKjlX/YJOSsGHJJQyd2iHZ9hkzoeCDNXuYNmspvVPimXvHufRJ9c92Lm1RlJOCS2HFTs+lfuHDkksYGtQjkZjICBvUN0FNVXnp2x38fPZKhmV15fXbzyXNj9u5tMXw3l2JjJCwXkzp10WUJjjEREWQ2zPBtoExQavsSBW/eqeYTzfsZ3xuOn++fhhx0e2z6r41OsdGkZ+RFNbJxVouYaogM5ni0iO4XDaob4KHqvLOyhIu/OMCvtpczq8vyuWZm0YEVGI5pSi7K6tKDlNVG56PFrfkEqYKMpM4Vl3HtgO+PPHAGOeUHa1ixsvLuf/11fRP68K8e0dz6+i+RLbzqvvWKspJpabOFbZjm9YtFqaGZCUDsHp3Jf3TEs5c2BgHqSrvrtrDQ++to6q2nl9flMvUH+UEbFI5pTDb/WTLJdsPfr+wMpxYyyVM9evehfiYSBvUNwGt7GgVt72ynPteX0W/7p0DvrXSWHJ8DAPTE8L24WHWcglTkRFCfkaSrdQ3AUlVeW+1u7Vysqaef5+Uy7QfB35rxVNRTgpvryihrt5FVGR4/S0fXt/WnKYgI4n1e49QU+dyOhRjvld+tJrbX1nOvXNWkdPN3VqZMSY4WiueinJSOF5Tz/q9R5wOpcNZcgljBVnJ1NS52LTvqNOhGNMwtlLKhX/8ki++K+dXkwYx947z6Ne9i9OhtVpRGD88zJJLGBuZk0KEwEfr9jkdiglz5UeruePVRq2Ve0Zz25h+QdlaaSw9MY4+qfGWXEx4SU+MY+xZ3Xlz+W7q6q1rzHS8U2MrP/3jl3y+qZwHJ7pbK/3Tgre14qkoO4WlOyrCbk2ZJZcwd21hb/YfqWbB5nKnQzFhpqbOxd2zV3LP7JX0Se3MvHt+zO1jg7+14qkwJ4VDJ2rZUh5ea8osuYS5cblpdOsSw5wlu50OxYQRVeVX76zl72v28sufDWTuHeeG7HqrkWE67mLJJcxFR0Zw5fBMPttYRtnRKqfDMWHiyc+3MHd5CfeNH8BdP+kf0tN0e6fEk54Ya8nFhJ9rCrOodylvLS91OhQTBt5bvYfff/wdlw/L4N5xA5wOp92JCIXZKSzZXhFWD+gL2eQiIn1FZKaIzHU6lkDXr3sXCrO78say3WF185uOt3xnBf/65mqKslP4rysHIxJa4yvNGZmTwr4jVWH1gD6vkouI3C8i60SkWERmi8gPHprgTRkvr/WCiJSJSHET5yaIyCYR2SIiD5ypHlXdpqrTWxNDOLq2sDfbDxwPu6a76Tg7Dx5nxsvLyUjuxDM3jSA2KvB2Mm4vRTnuRy6H01YwLSYXEckA7gHOUdV8IBKY0ooyaSKS4HGsfxOXnAVMaCKOSOBJYCKQB1wnInkiMlhEPvB4pbX0vczpJg3uQUJsFK8vtYF943+HT9QwddZSXKq8cEshXTvHOB1ShxqQ1oXk+GiWbD/odCgdxttusSigk4hEAfHAnlaUGQu8e6pFIyIzgCc8K1HVBUBT6b0I2NLQIqkB5gCXqupaVZ3s8Srz8nuZBvExUVw8tBfzivdSebLW6XBMCKmpc3HHq8spqTjJszedQ0435x9D3NEiIoRz+qSwdEf4PPa4xeSiqqXA74FdwF6gUlU/bkWZN4H5wBwRuQGYBlzjQ6wZQOM/q0sajjVJRFJF5GlgmIg82EyZi0Xk2cpK27wRYEphFlW1Lt5b3dTfDsb4TlV58O21LNpWwaNXFYTl1vOnjMxJYfuB45QdCY9Zmd50i3UFLgVygF5AZxG50dcyAKr6KFAFPAVcoqq+rCpqauSv2dFnVT2oqneoaj9VfaSZMu+r6m1JSUk+hBG6BmckkdszkdeX7nI6FBMi/vyPLby1wj3l+LJhzf4tGBa+32dsR3iMu3jTLTYe2K6q5apaC7wNnNeKMojIaCAfeAd4yMdYS4CsRu8zabp7zrSSiDClMIvi0iMUl1przrTNu6tK+e9PwmfKcUvO7pVIfEwkS8NkUN+b5LILGCUi8eKeNzgO2OBrGREZBjyHu4UzFUgRkYd9iHUpMEBEckQkBveEgfd8+LzxwmVDM4iJiuCNZTawb1pv2Y4KfvnmGopywmvK8ZlERUYwok/XsJkx5s2Yy2JgLrACWNvwmWcBRGSeiPQ6U5lG4oGrVXWrqrqAm4GdntcTkdnAQmCgiJSIyPSGOOqAu4GPcCeuN1R1ne9f2ZxJUnw0E/N78M7KUqpq650OxwShHQeOM+PlZWR07cQzN4bXlOOWFGWnsGn/USpPhP6kGa+eRKmqD9FEN5aqTmqpTKPz33i8r8XdkvEsd90Z6pgHzPMmZtN61xZm8e6qPcwv3hf2/eTGN4dP1DBt1lIUwnLKcUsKc1JQhWU7KxiXm+50OO0qZFfom9YblZNK75R45tjAvvFBTZ2L219ZTsmh8J1y3JKhWcnEREaExWJlSy7mByIihGsLs1i0rYIdB447HU67UFVeWbiDYf/nYxZ8Z48baCtV5YG317B4u005PpO46EiGZCWFxbiLJRfTpKtGZBIhhOTAftmRKm55cSn/8e46Dp2otV0J/ODP/9jC2ytKuX/8WdaV2oLC7BSKSys5UVPndCjtypKLaVJ6Yhw/GZjGm8tLQuoplR+u3ctPH1vA4u0H+c2lZ3P9yN58vqnMJi+0wakpx1cMy+CecU3t6GQaK8pJoc6lrNx12OlQ2pUlF9OsawuzKD9azeebgr/b6EhVLb94YxV3vraC3inx/P2e0dx0bjaT8ntyoqaeL61rrFWWNppy/IhNOfbKiD5diZDQ38TSkotp1k8GpdE9ITbou40WbzvIxMe+4m8rS7nngv68ded59Ovufkb7yL4pJMdHM794n8NRBp91eyptynErJMRFk9crMeQ3sbTkYpoVHRnBVSMy+XxTGfuDcD+k6rp6Hpm3gSnPLSI6Uph753n84qcDiW701MPoyAguzE3n0w37qakLne6/9rZq92Gue3YRnWOimDXVphz7qig7lZW7Dof0PWfJxZzRNee4n1I5d3mJ06H4ZNO+o1z25Lc8s2AbUwp78/d7RjO8d9cmy04c3IOjVXV8s/VAB0cZnJbuqODG5xeTHB/D67ePok+qTTn2VVFOCtV1LtaWHnY6lHZjycWcUU63zozMSQmap1S6XMrzX23j4j99TfnRKmbefA6PXDGYzrHNrxf+Uf9uJMRGMX+tdY215NstB/inmUtIS4jljdvPJbNrvNMhBaXCbPcfOqE87mLJxbRoSlEWOw+eYNG2wP4fofTwSW54fjEP/30DYwd256P7xni1Cjo2KpILctP4eP2+kJoZ529fbCpj6qylZKV0Ys7to+iR1KqHzRogtUss/dO6hPQmlpZcTIsm5vckIS4qYLfiV1X+trKUCY8tYE3JYR69soBnbxpBapdYr+uYmN+DQydqw2LldGt8vG4fM15eRv+0Lsy57VzSEiyxtFVRTgrLdhyi3hX4PQKtYcnFtCguOpLLhmbwYfG+gNtw7/CJGu6evZL7Xl/FwPQEPrx3DNcUZvk8JXbsWWl0io7kQ5s19gMfrNnDP7+2grN7JfHXW0eRYoP3fjEyJ4Wj1XVs2HvE6VDahSUX45VrC7OornPx7upSp0P53pfflfOzxxbwUfE+fvmzgbx++7n0Tm3dGECnmEjOH9idj9btwxWif0m2xtsrSrhn9kqG9U7mlelFJMVHOx1SyCjMdm+Rs3BraE5JtuRivJKfkcTZvRKZs8T5NS/Hquv41TtrufmFJSTGRfO3u37EXT/pT2RE2xbwTcjvQdnRalbsCp/nnJ/J7CW7+Jc3V3Nuv1RemlZEQpwlFn/qldyJoVnJPP7Z5pBsvVhyMV6bUpjF+r3OPqVy8baDTHx8AbOX7OL2MX15/+c/Jj/DP4+pvmBQGjGREdY1Brz07Q4efHst55/VnZk3FxIf49XTOYyPnrpxOF1io5g2ayn7KoNvLdmZWHIxXrtkaAaxURGObMVfVVvPwx+sZ8pzi4gQ4c3bz+XBSbnERftvVXhCXDSjB3RjfvG+oJh23V6e+XIrD723jp/mpfP0TSP8+t/YnK5nUideuKWQIydrmTZrKceqQ2czS0suxmtJnaK5aHBP3l25h5M1HbfR4+rdh5n8p695/uvt3DiyDx/eO5pzsttnS/cJ+T0oPXyStQ62zpyiqjzx2WYe+XAjkwt68uQNw21Llw6Q1yuRJ28Yzqb9R7nrtRUhMx3ekovxyTWFWRytruPD4r3tfq2aOhd/+HgTVzz1Lcer63hlehG/uSy/XbtoLsxLJypCwq5rTFX5vx9t4g+ffMcVwzN4fMqw07bJMe3r/IFp/ObSfL78rpz/eHddSLSc7e4xPhmZk0J2ajxz2nkzy437jnD5X77hiX9s4bKhGcy/bwyjB3Rv12sCJMfHcG6/1LDqGlNVHv77Bv7yxVauK+rN768a0ubJEcZ314/szZ3n92P2kl08s2Cb0+G0mSUX4xMR4drC3izZXsG28mN+r7/epTz1xVYu+dM37D9SxTM3jeC/rxlCUqeOm6k0Ib8H2w8cZ9P+ox12Tae4XMp/vFvMzK+3c8t52fz28nwiLLE45pc/Hcjkgp7814cb+WDNHqfDaRNLLsZnV47IIDJCeN3PT6ncfuA4Vz/9Lb+bv5FxuWl8dN8YfnZ2D79ewxs/zeuBCHwY4nuN1bvcjyZ+ddEu7hjbj4cuzrPnsTgsIkL4/dVDOKdPV37xxmqW7QjeHSMsuRifpSXEccGgNN5aXkqtHwYfXS7lpW93MPHxBWwpO8bjU4bylxuG+7R9iz91T4ilMDslpJ/xUu9Sfvnmat5YVsK94wbw/00YaIklQMRFR/LcP51DRnInZry8jO0HjjsdUqtYcjGtMqUwiwPHqvnHxrI21VNy6AQ3zlzMQ++tY2ROKh/fP5ZLh2Y4/otuYn4PNu0/2i5df06rq3fxizdW8fbKUv7lwrO4/8KzHP/vbU7XtXMML95SiIgw9cUlVByvcTokn4XsyigR6Qv8O5Ckqlc5HU+oGXtWd9IT3U+p9Oy6UlUqT9Zy4Fg15UdrOHi8mgNHqzl4vOb0Y8eq2V9ZTXSk8MgVg5nSij3B2suE/B787/fX82HxPu76Seg8F76u3sV9r6/igzV7+bcJA/nn80Pnu4Wa7G6dee6fzuG65xYx4+VlvHbryKBac+RVchGR+4FbAQXWAlNVtcqjTDLwPJDfUG6aqi5sTVAi8gIwGShT1fxGxycAjwORwPOq+l/N1aGq24DpIjK3NTGYM4tqeErlU19s5Revr+LA8RoOHnMnjIPHaqhrYn+uCIGUzrF06xJDty6x9OkdT3piHDeO6kNWSmA9F6RnkntrjvkhlFxq613cO2cl89bu48GJg7h9bD+nQzItGNGnK49dO5R/fm0F//LGav503bCgmXDRYnIRkQzgHiBPVU+KyBvAFGCWR9HHgfmqepWIxADxHvWkASdV9WijY/1VdUsTl50F/Bl4uVHZSOBJ4EKgBFgqIu+p6noRGQw84lHHNFVtW5+NOaMphb15a3kpi7dXkNolhvTEOM7ulUhql1i6dfmfJHLq38nxMUE1xXVifg8e+XAjuytOBFzy81VNnYt7Zq9k/rp9/PqiXG4d3dfpkIyXJg3uya8mDeK38zaSmdKJByfmOh2SV7ztFosCOolILe6kcdocORFJBMYAtwCoag3g2Uk4FrhTRCapapWIzAAuByZ5XkxVF4hItsfhImBLQ4sEEZkDXAqsV9W1uFs6pgNlpcSz6FfjnA6j3UzM78kjH25kfvE+ZowJ3l/GNXUu7vrrCj5Zv5//NTmPaT/OcTok46MZo/uyq+IEz3y5jayu8dw4qo/TIbWoxQF9VS0Ffg/sAvYClar6sUexvkA58KKIrBSR50Wks0c9bwLzgTkicgMwDbjGh1gzgMZzX0sajjVJRFJF5GlgmIg82EyZi0Xk2crK8Nvqw7Ssd2o8eT0TO2Q3gvZSXVfPna8u55P1+/k/l55tiSVIiQj//8Vn85OB3flf7xbzeRsn0nSEFpOLiHTF3ULIAXoBnUXkRo9iUcBw4ClVHQYcBx7wrEtVHwWqgKeAS1TVl6k4TfWnNLuEWlUPquodqtpPVT27zE6VeV9Vb0tK8s+uuib0TMzvwYpdh4Nyx9qq2nrueGU5n20s4+HL8vmnc7OdDsm0QVRkBH++fji5PRO5668rHN2d3BveTEUeD2xX1XJVrQXeBs7zKFMClKjq4ob3c3Enm9OIyGjcA/7vAA/5GGsJkNXofSYe3XPG+NvEwe6ZcB+tC641L1W19dz2ynI+31TOI1cMDopuFNOyzrFRvHBLIUmdopn+0lL2HD7pdEjN8ia57AJGiUi8uOeJjgM2NC6gqvuA3SIysOHQOGB94zIiMgx4DncraCqQIiIP+xDrUmCAiOQ0TBiYArznw+eN8Vn/tAT6p3UJqq6xkzX1zHh5GV9tLufRKwu4rqi30yEZP0pPjOPFqYUcr65n2qylVJ4MrEePn+LNmMti3C2RFbinIUcAzwKIyDwR6dVQ9OfAayKyBhgK/NajqnjgalXdqqou4GZgZ1PXFJHZwEJgoIiUiMh0Va0D7gY+wp3c3lDVdb58WWNaY2J+D5Zsr+DgsWqnQ2nRyZp6pr+0lK+3HOD/XjWEawqzWv6QCTqDeiTy1I3D2VJ2jML//JSbZi7m2QVbWb/nSMBsuCqBEoiTzjnnHF22bJnTYZgAtW5PJRc98TWPXDE4oFsBJ2rqmDZrKUu2V/Df1wzh8mGZTodk2tnKXYf4YM1evtpcznf73UPY3brEMnpAN0YP6MaP+3cjLTGu3a4vIstV9ZymzoXsCn1j/CWvZyK9U+L5sHhfwCaX49V1TH1xKct2VvDHa4dy6dBmJ1KaEDKsd1eG9e4KwL7KKr7ecoCvNpez4Lty3llZCsCgHgkNyaY7RTkpHbbK35KLMS0QESbm92Dm19upPFFLUnzHbf/vjWPVddzywhJW7j7M41OGcfGQXi1/yIScHklxXDUik6tGZOJyKev3Hvk+2bz07U6e+2o7MVERFGWnfJ9sBvVIaLcV/9YthnWLmZat3HWIy//yLf999RCuHBE43U1Hqmq55YUlrCmp5InrhjFpcE+nQzIB6GRNPYu3H+TrzQf4avOB759V1K1LLFN/lN3qLY6sW8yYNhqSmUzPpDg+LN4XEMlFVfl6ywF+N38jG/ce5c/XD2dCfsc/+8YEh04xkZw/MI3zB6YBsP9IFV9tPsDXm8vpEts+acCSizFeiIgQfnZ2D/66ZBfHquva7X/IlhyvruPtFSXM+nYHW8uP061LDE/fOILxeemOxGOCU3ri/3ShtRdLLsZ4aWJ+D2Z9u4PPN5Z1+LjGjgPHeXnhTt5ctpuj1XUMyUzij9cOYdLgnsRGBc827CZ8WHIxxkvnZKfQrUsM84v3dUhyUVW+2nzAndA2lREpwkUFPbnlvOzvZwgZE6gsuRjjpcgI4adn9+BvK0upqq1vtymdxxq6vl76vusrlnsuGMANI3u365oFY/zJkosxPpiY34O/Lt7Fl9+V/+AJnG2148BxXlq4g7nLStxdX1nJPHbtUCYO7mFdXyboWHIxxgej+qaS1Cma+cX7/JJcXC7lqy0HmPXNdr74rpyoCOGiwT252bq+TJCz5GKMD6IjI7gwL52P1u2jps5FTJQ3e7/+UFVtPW+tKGHm19vZ1tD1de+4AVw/sjdpCdb1ZYKfJRdjfDRpcA/mLi/hm60H+EnDugFvHThWzSsLd/LKop1UHK+hIDOJx64dyqTBPVudqIwJRJZcjPHRj/p3IyE2ivlr93mdXLaUHWPm19t4a0UpNXUuxuemM2N0DkU5KbifZGFMaLHkYoyPYqMiuSA3jY/X7+M/6/OJimy6xaGqLNpWwfNfbeOzjWXERkVw1YhMpv84h37du3Rw1MZ0LEsuxrTCxPwevLtqD0u2V3Be/26nnautdzFv7V6e+2obxaVHSO0cw33jB3DTqD6kdol1KGJjOpYlF2NaYexZaXSKjuTD4n3fJ5ejVbW8vnQ3L3y9nT2VVfTt3pnfXj6YK4ZndNg258YECksuxrSCeyPA7ny0bh+3j+3LS9/uYM4S99YsI3NS+M1l+fxkYFq7bWduTKCz5GJMK03I78GHxfsY/ejnRIh7fcqto3MoyEx2OjRjHGfJxZhWGp+bznn9UsnrmcjUH+eQkdzJ6ZCMCRiWXIxppc6xUfx1xiinwzAmINmqLWOMMX5nycUYY4zfWXIxxhjjd5ZcjDHG+F3IJhcR6SsiM0VkrtOxGGNMuPEquYjI/SKyTkSKRWS2iDS5J7iIRIrIShH5oC1BicgLIlImIsUexyeIyCYR2SIiD5ypDlXdpqrT2xKHMcaY1mkxuYhIBnAPcI6q5gORwJRmit8LbGimnjQRSfA41r+ZemYBEzzKRgJPAhOBPOA6EclrODdYRD7wePm2F7oxxhi/8bZbLAroJCJRQDywx7OAiGQCFwHPN1PHWODdU60eEZkBPNFUQVVdAFR4HC4CtjS0SGqAOcClDeXXqupkj1dZS19KRC4WkWcrKytbKmqMMcYHLS6iVNVSEfk9sAs4CXysqh83UfQx4N+AhCbOoapvikgOMEdE3gSmARf6EGsGsLvR+xJgZHOFRSQV+E9gmIg8qKqPNBHT+8D7InK5iBwGzpRlkpo53w040HL4Aae57xPo12pLXb5+1tvy3pRrqcyZzgfjPWb3l//KB/L91afZM6p6xhfQFfgH0B2IBv4G3OhRZjLwl4Z/nw98cIb65gBHgO4tXDcbKG70/mrg+UbvbwL+1FL83r6AZ1tzHljmrxg68tXS9w3Ua7WlLl8/6215b8q19v5qOBd095jdX/4rH6z3lzfdYuOB7aparqq1wNvAeR5lfgRcIiI7GpLHBSLyqmdFIjIayAfeAR7y4tqNlQBZjd5n0kT3XBu838bzwaYjv48/r9WWunz9rLflvSln91dwXMvuLz+RhszVfAGRkcALQCHubrFZuDPdn5opfz7wr6o62eP4MGA27nGZ7cCrwDZV/XUz9WTjbgHlN7yPAr4DxgGlwFLgelVd1/LXbD8iskxVz3EyBhPa7B4z7am97q8WWy6quhiYC6wA1jZ85tmGoOaJSC8vrxUPXK2qW1XVBdwM7GyqoIjMBhYCA0WkRESmq2odcDfwEe4ZaW84nVgaPOt0ACbk2T1m2lO73F8ttlyMMcYYX4XsCn1jjDHOseRijDHG7yy5GGOM8TtLLu3INs80/iYinUXkJRF5TkRucDoeE3r89XvLkkszbPNM01F8vNeuAOaq6gzgkg4P1gQlX+4xf/3esuTSvFl4uXmmbZxp2mgW3m/Umsn/bINU34ExmuA2Cx82A/aHFvcWC1equqBhIWdj32+eCSAic4BL1b1v2WSMaQVf7jXcO1VkAquwPw6Nl3y8x9b745p2c/qmqc0zM5orLCKpIvI0DZtntndwJqQ0d6+9DVwpIk8RelvGmI7V5D3mr99b1nLxjTRxrNlVqKp6ELij/cIxIazJe01VjwNTOzoYE5Kau8f88nvLWi6+ae/NM405xe41097a9R6z5OKbpcAAEckRkRjcT+R8z+GYTGiye820t3a9xyy5NCPINs80QczuNdPenLjHbONKY4wxfmctF2OMMX5nycUYY4zfWXIxxhjjd5ZcjDHG+J0lF2OMMX5nycUYY4zfWXIxxhjjd5ZcjDHG+J0lF2OMMX73/wAEDhOu7uiWOAAAAABJRU5ErkJggg==\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "plt.loglog(E, A_blr)" - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": {}, - "outputs": [], - "source": [ - "abs_nu_agnpy = abs_blr.absorption(nu)" - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "[]" - ] - }, - "execution_count": 13, - "metadata": {}, - "output_type": "execute_result" - }, - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEFCAYAAAAPCDf9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfeUlEQVR4nO3de5QU5Z3/8fd3hiEGRLwwKnKZwSxG0ESFEURcRHdDwKCoMSs6hJWYIApodM2KwXhBMWajiSJE5Rg2URFDNApeMT9jTPDKjIKAaILcxQvq8YIYDfL9/fE069AMTM3Y3dVd9Xmd06enq57q/lbq+KFS9dTzmLsjIiLpUBZ3ASIiUjgKfRGRFFHoi4ikiEJfRCRFFPoiIinSKu4CGtOhQwevrq6OuwwRkZJRX1//trtXNtWuKEO/urqaurq6uMsQESkZZrY6Sjtd3hERSRGFvohIiij0RURSRKEvIpIiCn0RkRRpMvTNbIaZvWVmS3aw3sxsipktN7MXzaxXg3WDzeyVzLoJuSw828yZ0KEDmG37Ki8P761ahfcOHcKrrGzbv6urw3fMnBn+brhs6/c3trypddk1RmnX0va53j5f3yUiMXL3nb6AAUAvYMkO1h8HPAwYcATwbGZ5OfAqsD/QGlgE9Gzq99yd3r17e3PccYd7RYU7fLFXRYV769bbLmvTxv3ss8N79vI77givHa3LrjFKu5a2z/X2+fouEckPoM4j5Kt5hKGVzawaeMDdD25k3S3An919VubzK8BAoBq43N2/mVl+ceYfmZ829Xs1NTXenH761dWwOlIP1ZYxC1GXrX378P7++42vO/fczz9PmRKtXUvb53r7KN9VVQWrVjXvu0QkP8ys3t1rmmyXg9B/ALjG3ednPj8GXEQI/cHu/v3M8u8Cfd193A5+YzQwGqBr1669VzcjxcvKGg/luJl9/vfO6mvYrqXtc719lO8ygy1bmvddIpIfUUM/FzdyG4sQ38nyRrn7dHevcfeaysomnyTeRteuzWrebOXljS+vqgqvHa3bsuXzV9R2LW2f6+2jfFe+/3cXkdzLReivA7o0+NwZWL+T5Tk3eTJUVHzx76mogNatt13Wpg2MHh3es5dPnhxeO1qXXWOUdi1tn+vtm/ougOHDm/9dIhKzKBf+CZdqdnQj91tseyP3uczyVsAKoBuf38g9KMrvNfdGrnu4qbjXXtvfnC0rC+/l5eF9r73Cy2zbv6uqPr8xW1W17bKt39/Y8qbWZdcYpV1L2+d6+x19V5cu4bXrru719S3/ThHJHXJ1I9fMZhFuzHYA3gQuAyoy/2DcbGYGTAUGA5uAUe5el9n2OOB6Qk+eGe4e6TyzuTdypfDWr4d+/eCTT+Cpp2D//eOuSCTdcnojt9AU+qVh2TLo3z886/DUU+FdROJRyBu5klI9esD998PatTB0KHz0UdwViUhTFPryhfTvD7NmwYIF4cbu5s1xVyQiO6PQly/sxBNh2jR44AEYM6Y4n5kQkaAoZ86S0jNmDLz2Glx1FXTqBFdcEXdFItIYhb7kzKRJoVfPpEmw335w1llxVyQi2RT6kjNmcPPN8MYbcM450LEjnHBC3FWJSEO6pi85VVEBs2dD797hxu7TT8ddkYg0pNCXnGvbFh58MFzbHzoUXn457opEZCuFvuRFZSXMmxcmrxk8OFzrF5H4KfQlb/bfHx56CN55B447rvEx+UWksBT6kle9e8M998DSpXDyyWGsHhGJj0Jf8m7QIJgxA/70Jxg1ShOviMRJXTalIL773XBdf8KE0If/2mvjrkgknRT6UjD//d+wbh1cd13o2XP++XFXJJI+Cn0pGDO4/vrw8NYFF4SHtzT7lkhhKfSloMrL4fbb4a23YORI2HtvOPbYuKsSSQ/dyJWC22UXuO8+OOAAOOkkWLQo7opE0kOhL7HYYw945BHYbTcYMgRWr467IpF0UOhLbDp3DsH/8cfhqd133om7IpHkU+hLrA46CObMgZUrw0TrXbtCWRlUV8PMmXFXJ5I8Cn2J3YABYez9v/89zLfrHi73jB6t4BfJNYW+FIU5c7ZftmkTTJxY+FpEkkyhL0VhzZrmLReRllHoS1Ho2rV5y0WkZRT6UhQmT4Y2bbZdVlYW5tsVkdxR6EtRqK2F6dOhqioM19ChQxiN89VX465MJFkU+lI0amth1aoQ9hs2wBlnwFVXwfz5cVcmkhwKfSlaU6ZAt24wYgS8917c1Ygkg0Jfila7dqGf/rp1MHZs3NWIJEOk0DezwWb2ipktN7MJjazfw8zuNbMXzew5Mzu4wbpVZrbYzBaaWV0ui5fk69sXLr8c7rwT7rgj7mpESl+ToW9m5cA0YAjQEzjNzHpmNfsxsNDdvw6MBG7IWn+Mux/q7jU5qFlS5uKL4aij4JxzwnANItJyUc70+wDL3X2Fu38K3AUMy2rTE3gMwN1fBqrNbJ+cViqpVV4ezvLLysLN3s2b465IpHRFCf1OwNoGn9dlljW0CDgZwMz6AFVA58w6Bx41s3ozG72jHzGz0WZWZ2Z1GzZsiFq/pERVFdx8Mzz9dOjTLyItEyX0rZFlnvX5GmAPM1sIjAdeALaej/V3916Ey0NjzWxAYz/i7tPdvcbdayorKyMVL+kyfHiYYH3SJHjqqbirESlNUUJ/HdClwefOwPqGDdz9A3cf5e6HEq7pVwIrM+vWZ97fAu4lXC4SaZGpU8NZf20tfPBB3NWIlJ4oob8A6G5m3cysNTAcmNuwgZntnlkH8H3gL+7+gZm1NbN2mTZtgUHAktyVL2mz226hG+faterGKdISTYa+u28GxgHzgGXAbHdfamZjzGxMplkPYKmZvUy4jHNeZvk+wHwzWwQ8Bzzo7o/keickXfr1g0svDTd377wz7mpESou5Z1+ej19NTY3X1alLv+zY5s0wcCAsXhwmVq+ujrsikXiZWX2UbvF6IldKUqtWnz+sNWKEunGKRKXQl5JVXQ033QRPPgk//Wnc1YiUBoW+lLTTTw89ea64IvThF5GdU+hLyZs2Dbp0UTdOkSgU+lLy2rcP1/dXr4bx4+OuRqS4KfQlEfr3h5/8BG67De66K+5qRIqXQl8S45JLQh/+MWPCWb+IbE+hL4mxtRvnli1hjJ7PPou7IpHio9CXRNl//3Bj969/hWuuibsakeKj0JfEGTECTjsNLrsMnn027mpEiotCXxLHDH71K+jcOXTj/PDDuCsSKR4KfUmk3XcP1/dXroRzz427GpHiodCXxDrqKJg4EX7zG5g9O+5qRIqDQl8S7Sc/gb594ayzYM2auKsRiZ9CXxKtoiJMurJ5M4wcqW6cIgp9SbyvfCVMs/jEE/A//xN3NSLxUuhLKowcCaeeGmbcWrAg7mpE4qPQl1QwC2Pvd+wYhmPeuDHuikTiodCX1Nhjj9CN89VX4bzzmm4vkkQKfUmVAQPg4othxgy4++64qxEpPIW+pM7ll8Phh8Po0bB2bdzViBSWQl9Sp6IC7rwTPv1U3TglfRT6kkr/8i9w443w5z/DtdfGXY1I4Sj0JbXOOANOOSVMvlJfH3c1IoWh0JfUMoPp02HffUM3zo8+irsikfxT6Euq7bEH3H47/P3v8MMfxl2NSP4p9CX1Bg6Eiy6CW2+FykooK4Pq6jBmj0jStIq7AJFicOCBIezffjt8Xr06dOmEMBGLSFLoTF+EMLXili3bLtu0KYzHL5IkkULfzAab2StmttzMJjSyfg8zu9fMXjSz58zs4KjbihSDHY21rzH4JWmaDH0zKwemAUOAnsBpZtYzq9mPgYXu/nVgJHBDM7YViV3Xrs1bLlKqopzp9wGWu/sKd/8UuAsYltWmJ/AYgLu/DFSb2T4RtxWJ3eTJ0KbN9stHjix8LSL5FCX0OwENRyhZl1nW0CLgZAAz6wNUAZ0jbktmu9FmVmdmdRs2bIhWvUiO1NaGPvtVVaH/fufOof/+TTfBqlVxVyeSO1FC3xpZ5lmfrwH2MLOFwHjgBWBzxG3DQvfp7l7j7jWVlZURyhLJrdraEPBbtoSB2J54IkyzePzx8OGHcVcnkhtRQn8d0KXB587A+oYN3P0Ddx/l7ocSrulXAiujbCtSrA44AGbPhmXLwj8IGphNkiBK6C8AuptZNzNrDQwH5jZsYGa7Z9YBfB/4i7t/EGVbkWL2jW/A9dfD/fer+6YkQ5MPZ7n7ZjMbB8wDyoEZ7r7UzMZk1t8M9ABuM7PPgJeAM3e2bX52RSQ/xo6FpUvhZz+Dnj11c1dKm7k3eok9VjU1NV5XVxd3GSL/55//hG9+E558MgzH3K9f3BWJbMvM6t29pql2eiJXJIKKCvj976FLFzjxRD20JaVLoS8S0V57hWv7//gHnHACbNwYd0UizafQF2mGHj3gd7+DxYvDtf3s8XpEip1CX6SZBg+G666De++FSy+NuxqR5tHQyiItcN55oUfP5MmhR8/pp8ddkUg0OtMXaQEzmDYNBgyA730Pnnsu7opEolHoi7RQ69Zwzz3QsSMMGwbr1sVdkUjTFPoiX0CHDqFHz8aNoSvnpk1xVySycwp9kS/o4INh1ix4/nk44wz16JHiptAXyYGhQ8MwDb//PVx5ZdzViOyYeu+I5MiFF4YePZdfHnr0fOc7cVcksj2d6YvkiBnccgsceST8539CfX3cFYlsT6EvkkNf+lJ4aKuyMvToef31uCsS2ZZCXyTH9t4b5s6F994Lwf/xx3FXJPI5hb5IHhxyCNxxByxYAGeeCUU4grmklEJfJE9OPBGuvjp057z66rirEQnUe0ckjyZMCD16Lrkk9Og56aS4K5K005m+SB6Zwa23Qt++MGIELFwYd0WSdgp9kTzbZRe47z7Yc88w+cobb8RdkaSZQl+kAPbdF+bMgbffDpd4/vGPuCuStFLoixRIr15w++3wzDMwerR69Eg8FPoiBfTtb8OkSSH8f/7zuKuRNFLvHZECu+SS0KNnwgQ48MBwnV+kUHSmL1JgZvC//wu9e0NtbZhkXaRQFPoiMfjyl0OPnnbt4PjjYcOGuCuStFDoi8SkU6fQo+fNN+Hkk+GTT+KuSNJAoS8So8MPh9/8BubPh7PPVo8eyT/dyBWJ2amnhhu7V14Zpl684IK4K5Iki3Smb2aDzewVM1tuZhMaWd/ezO43s0VmttTMRjVYt8rMFpvZQjOry2XxIklx+eWhO+ePfgQPPRR3NZJkTYa+mZUD04AhQE/gNDPrmdVsLPCSux8CDASuM7PWDdYf4+6HuntNbsoWSZayMvjtb8OQzMOHw0svxV2RJFWUM/0+wHJ3X+HunwJ3AcOy2jjQzswM2BV4F9ic00pFEq5t23Bjt02b0KPn7bfjrkiSKErodwLWNvi8LrOsoalAD2A9sBg4z923ZNY58KiZ1ZvZ6C9Yr0iidekSgv+11+CUU+DTT+OuSJImSuhbI8uy+xh8E1gI7AccCkw1s90y6/q7ey/C5aGxZjag0R8xG21mdWZWt0GdliXF+vaFX/8anngCxo9Xjx7JrSihvw7o0uBzZ8IZfUOjgD94sBxYCRwI4O7rM+9vAfcSLhdtx92nu3uNu9dUVlY2by9EEqa2Fi6+GKZPhxtvjLsaSZIoob8A6G5m3TI3Z4cDc7ParAH+DcDM9gG+Cqwws7Zm1i6zvC0wCFiSq+JFkuyqq8LE6uefD48+Gnc1khRNhr67bwbGAfOAZcBsd19qZmPMbEym2ZXAkWa2GHgMuMjd3wb2Aeab2SLgOeBBd38kHzsikjRlZWFy9YMPhv/4D3j55bgrkiQwL8ILhjU1NV5Xpy79IgCrV4cnd9u3h2efDTNwiWQzs/oo3eI1DINIkauqgnvvhTVr4F//NXwuK4Pqapg5M+7qpNQo9EVKQP/+cMYZ4aGtNWtCj57Vq8MMXAp+aQ6FvkiJmDdv+2WbNsHEiYWvRUqXQl+kRKxZ07zlIo1R6IuUiK5dG1/euXNh65DSptAXKRGTJ4dxebLtsotm3pLoFPoiJaK2NjyhW1UV5tmtqgrDNKxdC0ccoX78Eo1CX6SE1NbCqlWwZUt4nzIFHn8cNm6Efv3C3yI7o9AXKXFHHBEe2tpvPxg0KEy/KLIjCn2RBKiuhiefhIEDYdSo0I1zy5amtpI0UuiLJMTuu4epFn/wA7j6ajj9dPj447irkmKjidFFEqSiAm65Bbp3h4suCk/tzpkDe+8dd2VSLHSmL5IwZmGC9bvvhkWLwjX/ZcvirkqKhUJfJKFOPjnMvrVpU+jZ89hjcVckxUChL5Jghx8eevZ06QKDB4dpGCXdFPoiCVdVBfPnw7HHwve/DxMmqGdPmin0RVKgfXt48EEYMwZ+9jM49VT17Ekr9d4RSYlWreBXvwo9ey68MIzOOXcu7LNP3JVJIelMXyRFzOCCC+APf4AlS6BvX1i6NO6qpJAU+iIpdOKJ8Je/wCefwJFHwh//GHdFUigKfZGU6t0bnnsuDOEwZEgYwVOST6EvkmJduoSePYMGwVlnhYe61LMn2RT6IinXrl24oTt2LFx7LZxySnigS5JJoS8itGoFN94I118P990HRx8Nr78ed1WSDwp9EQFCz57zzgsDtC1bFnr2LF4cd1WSawp9EdnG8cfDX/8Kn30G/fvDI4/EXZHkkkJfRLZz2GFhzJ6vfAWGDoWbboq7IskVhb6INKpz53DGP3gwnHNOeKjrs8/irkq+KIW+iOzQrruGa/znngu//GUYrvmjj+KuSr6ISKFvZoPN7BUzW25mExpZ397M7jezRWa21MxGRd1WRIpbeTnccEPo3fPAAzBgAKxfH3dV0lJNhr6ZlQPTgCFAT+A0M+uZ1Wws8JK7HwIMBK4zs9YRtxWREjBuXOjP/7e/hZ49ixbFXZG0RJQz/T7Acndf4e6fAncBw7LaONDOzAzYFXgX2BxxWxEpEd/6VniC1x2OOipMxC6lJUrodwLWNvi8LrOsoalAD2A9sBg4z923RNwWADMbbWZ1Zla3YcOGiOWLSKEdckgYs+eAA0L3zqlT465ImiNK6Fsjyzzr8zeBhcB+wKHAVDPbLeK2YaH7dHevcfeaysrKCGWJSFz22y+M0jl0KIwfHx7qUs+e0hAl9NcBXRp87kw4o29oFPAHD5YDK4EDI24rIiWobdswLv/558OUKWG45o0b465KmhIl9BcA3c2sm5m1BoYDc7ParAH+DcDM9gG+CqyIuK2IlKjycvjFL8KMXA89FHr2vPZa3FXJzjQZ+u6+GRgHzAOWAbPdfamZjTGzMZlmVwJHmtli4DHgInd/e0fb5mNHRCQ+Z58d5uBdvhz69IEXXoi7ItkRc2/0EnusampqvK6uLu4yRKSZFi8OPXzefRdmzQo3eqUwzKze3WuaaqcnckUkZ772tTBmT48e4Rr/lClxVyTZFPoiklMdO8Kf/wzDhoVePePHw+bNcVclWyn0RSTn2raFu++GCy8M/fiHDYMPP4y7KgGFvojkSVkZ/PzncMstMG9eeIJ37dqmt5P8UuiLSF6NHh26c65aFcbsqa+Pu6J0U+iLSN4NGgRPPgmtW4e+/HPmxF1Rein0RaQgDj4YnnkGDjoITjopjM9fhD3GE0+hLyIFs+++oWfPySeHmbi+8Q2oqgrX/6urYebMuCtMvlZxFyAi6dKmDcyeHXr0PPDA58tXrw7X/wFqa+OpLQ10pi8iBVdWFp7ezbZpE0ycWPh60kShLyKxWLOm8eWrV8M77xS2ljRR6ItILLp23fm688+HdesKV09aKPRFJBaTJ4fr+w21aQPXXAPf/naYiH3//eHMM8O8vJIbCn0RiUVtLUyfHnrvmIX36dPhoovgttvCMM1nnQV33gkHHgjf+Q48/3zcVZc+Da0sIkXtrbfghhtg2jR4//3woNfFF8PRR4d/LCTQ0Moikgh77x0uBa1eHS79LFoExxwDRx4ZnuzdsiXuCkuLQl9ESkL79uHSz8qVYXrGN98MY/Z//etw++3wz3/GXWFpUOiLSEn58pfD9Ix/+xvccUe4xDNyJHTvHi4Bffxx3BUWN4W+iJSkVq3CzeBFi+D++6FTJxg3Lgzn8NOfhuv/sj2FvoiUtLIyGDoU5s+HJ56AXr3gxz8Off0nTAiXgeRzCn0RSQSzMGzzww+Hrp2DB4dJXKqq4Jxzwr0AUeiLSAIddhj87nfw8svhev+vfx2u+Y8YAUuWxF1dvBT6IpJY3buHB75WrIAf/hDuuw++9jU4/nh46qm4q4uHQl9EEq9TJ7j22jDI2xVXwNNPQ//+4QGvRx5J12QuCn0RSY0994RLLw0Pev3yl+H/AQwZAr17hzH+P/ss7grzT6EvIqnTtm243PPqqzBjRhjH/9RToUcPuPVW+OSTuCvMH4W+iKRW69YwahQsXQp33w277QY/+EEY3fMXv4CNG+OuMPcU+iKSeuXlYTjnBQvg0Ufhq1+F//qv0Nf/ssuSNamLQl9EJMMsTNb+pz/BM8+Efv+TJiVrUpdIoW9mg83sFTNbbmYTGln/IzNbmHktMbPPzGzPzLpVZrY4s07jJYtISejbN3TxXLIETjll20ldXnkl7uparsnQN7NyYBowBOgJnGZmPRu2cfefu/uh7n4ocDHwhLu/26DJMZn1TY71LCJSTA46CH7723DTd+ukLj16hEld6uvjrq75opzp9wGWu/sKd/8UuAsYtpP2pwGzclGciEixqKoKZ/urV4dJXP74R6ipCZO6PP546fT1jxL6nYC1DT6vyyzbjpm1AQYD9zRY7MCjZlZvZqN39CNmNtrM6sysbsOGDRHKEhEpvK2TuqxZEyZ1efFFOPZY6NevNCZ1iRL6jU1ItqN/044Hnsy6tNPf3XsRLg+NNbMBjW3o7tPdvcbdayorKyOUJSISn912C5O6rFoFN90UpnUshUldooT+OqBLg8+dgfU7aDucrEs77r4+8/4WcC/hcpGISCLssguMGRMmdZk5Mwz1vHVSl6lTi29SlyihvwDobmbdzKw1IdjnZjcys/bA0cCcBsvamlm7rX8Dg4CUj3EnIknUqhWcfvq2k7qMHx/uBVx9Nbz3XtwVBk2GvrtvBsYB84BlwGx3X2pmY8xsTIOmJwGPuvtHDZbtA8w3s0XAc8CD7v5I7soXESkuZttO6tK7N0ycGMJ/wgR4442Y6/MivOVcU1PjdXXq0i8iyfDCC+Gm7913Q0UFfO978KMfQbduufsNM6uP0i1eT+SKiORZU5O6zJwZ5vYtKwvvM2fmrxad6YuIFNj69WFAt5tvho8+CmP/NBzWuU2bMPlLbW3079SZvohIkdpvv88ndWnffvtx/DdtCvcB8kGhLyISkz33hA8+aHzdmjX5+U2FvohIjLp2bd7yL0qhLyISo8mTwzX8htq0CcvzQaEvIhKj2tpw07aqKvTxr6pq/k3c5miVn68VEZGoamvzF/LZdKYvIpIiCn0RkRRR6IuIpIhCX0QkRRT6IiIpUpRj75jZBmB1MzbpALydp3KKXZr3HdK9/2ned9D+Z+9/lbs3Oe1gUYZ+c5lZXZSBhpIozfsO6d7/NO87aP9buv+6vCMikiIKfRGRFElK6E+Pu4AYpXnfId37n+Z9B+1/i/Y/Edf0RUQkmqSc6YuISAQKfRGRFCmZ0DezwWb2ipktN7MJjaw3M5uSWf+imfWKo858ibD/A83sfTNbmHldGked+WBmM8zsLTNbsoP1ST/2Te1/ko99FzN73MyWmdlSMzuvkTaJPP4R9735x97di/4FlAOvAvsDrYFFQM+sNscBDwMGHAE8G3fdBd7/gcADcdeap/0fAPQCluxgfWKPfcT9T/Kx7wj0yvzdDvhbWv7bj7jvzT72pXKm3wdY7u4r3P1T4C5gWFabYcBtHjwD7G5mHQtdaJ5E2f/Ecve/AO/upEmSj32U/U8sd3/d3Z/P/P0hsAzolNUskcc/4r43W6mEfidgbYPP69h+56O0KVVR962fmS0ys4fN7KDClFYUknzso0r8sTezauAw4NmsVYk//jvZd2jmsS+VmbOskWXZfU2jtClVUfbtecLYGxvN7DjgPqB7vgsrEkk+9lEk/tib2a7APcAP3f2D7NWNbJKY49/Evjf72JfKmf46oEuDz52B9S1oU6qa3Dd3/8DdN2b+fgioMLMOhSsxVkk+9k1K+rE3swpC6M109z800iSxx7+pfW/JsS+V0F8AdDezbmbWGhgOzM1qMxcYmbmTfwTwvru/XuhC86TJ/Tezfc3MMn/3IRzbdwpeaTySfOyblORjn9mvXwPL3P0XO2iWyOMfZd9bcuxL4vKOu282s3HAPEJPlhnuvtTMxmTW3ww8RLiLvxzYBIyKq95ci7j/pwBnm9lm4GNguGdu75c6M5tF6KXQwczWAZcBFZD8Yw+R9j+xxx7oD3wXWGxmCzPLfgx0hcQf/yj73uxjr2EYRERSpFQu74iISA4o9EVEUkShLyKSIgp9EZEUUeiLiBRAUwPnZbW9wMxeygwg95iZVWWWH9NgcLWFZvYPMzuxWXWo946ISP6Z2QBgI2GcoIObaHsMYeC4TWZ2NjDQ3U/NarMnoZtqZ3ffFLUOnemLiBRAYwPnmdlXzOwRM6s3s7+a2YGZto83CPJnCE8ZZzsFeLg5gQ8KfRGROE0Hxrt7b+BC4FeNtDmTMHR0tuHArOb+YEk8kSsikjSZgdSOBH6fGUkB4EtZbUYANcDRWcs7Al8jPKXfLAp9EZF4lAHvufuhja00s38HJgJHu/snWav/A7jX3f/Zkh8VEZECywyTvNLMvgP/N+3jIZm/DwNuAU5w97ca2fw0WnBpB9R7R0SkIBoOnAe8SRg470/ATYSpESuAu9x9kpn9P8Llm62jha5x9xMy31MNPAl0cfctza5DoS8ikh66vCMikiIKfRGRFFHoi4ikiEJfRCRFFPoiIimi0BcRSRGFvohIivx/75w9UPtCY+MAAAAASUVORK5CYII=\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "plt.plot(nu, abs_nu_agnpy, 'bo-')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": 22, - "metadata": {}, - "outputs": [], - "source": [ - "E = np.logspace(-1, 1, 20)*u.TeV\n", - "lp = log_parabola(\n", - " E,\n", - " amplitude=1e-12/((u.cm**2)*(u.s)*(u.TeV)),\n", - " reference=1 * u.TeV,\n", - " alpha=2.3,\n", - " beta=0.5\n", - " )\n", - "sed = E*E*lp" - ] - }, - { - "cell_type": "code", - "execution_count": 23, - "metadata": {}, - "outputs": [], - "source": [ - "# Set BLR parameters\n", - "L_disk=6e45*u.Unit(\"erg s-1\")\n", - "R_line = 2.45*1e17*u.cm\n", - "r_blob_bh = 1e18*u.cm\n", - "z = 1\n", - "xi_line = 0.1\n", - "\n", - "E = np.logspace(-1, 1, 20)*u.TeV\n", - "nu = E.to(\"Hz\", equivalencies=u.spectral())\n", - "\n", - "blr = SphericalShellBLR(L_disk, xi_line, \"Hbeta\", R_line)\n", - "abs_blr = Absorption(blr, r=r_blob_bh, z=z)\n", - "tau_all_lines = abs_blr.tau_blr_all_lines_cubepy(nu)\n", - "A_blr = get_absorbtion_tau(np.array(tau_all_lines))" - ] - }, - { - "cell_type": "code", - "execution_count": 24, - "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAGtCAYAAAB5r18AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABTXElEQVR4nO3deZxN9R/H8dd3Zuz7EsoSESlKSdkq9Csq28i+Za1EiBK/pI2oLKmGIWQpe2SJUPGzpLK0kUiWkCJLjJ35/v44Y7LMmHvHvffce+f9fDzuY+Z+z/d8z+dMp/txzj3n8zXWWkRERCR5EW4HICIiEuyULEVERFKgZCkiIpICJUsREZEUKFmKiIikIMrtANySN29ee80115AlS5Yr9jt27NgV+6S0PBS5vU/+2r4vxk3tGN6u52l/T/qlxWMY3N2vcDyGvV3X130DdRyvW7fub2vtNZctsNamyVf58uXt0qVLbUpS6uPJGKHG7X3y1/Z9MW5qx/B2PU/76xhOnpv7FY7HsLfr+rpvoI5jYK1NImfoMqyIiEgKlCxFRERSoGQpIiKSAiVLERGRFChZioiIpEDJUkREJAVKliIiIilQshQREUmBkqWIiEgKlCxFRERSoGQpIiKSAiVLERGRFChZioiIpEDJUkREJAVKliIiIilIs5M/iwSLP/74g1WrVrFy5UpWrlzJ5s2byZQpE9deey25c+cmT5485MmTJ8nfd+zYwZ9//knu3LlJnz6927siEraULEUCKD4+nl9++SUxMa5atYpt27YBkDlzZipWrEiHDh3YsmULGTJk4MCBA/zyyy8cPHiQAwcOcObMmWTHzpYt22UJtWzZspQoUSJQuycStkIyWRpjbgBeAHJYaxsmtNUHHgHyATHW2sXuRSjiOHXqFOvWrWPKlCkMHTqUVatWcfDgQQDy5ctH1apV6dKlC1WrVqVcuXKkS5cOgGXLllGtWrWLxrLWEhcXl5g4Dxw4wIoVKyhQoAAHDhy4qP3gwYNs376dadOmYYxh7NixtG3blvr165MpU6ZA/xlEQl7Ak6UxZhxQG9hnrS1zQXstYDgQCYyx1g5Kbgxr7TagvTFm5gVtnwCfGGNyAYMBJUsJuEOHDvHVV18lnjmuWbOGU6dOAVCqVCmio6OpWrUqVapUoUSJEhhjPB7bGEO2bNnIli0b119/PQDp0qW7LKleaPv27bzyyissW7aM5s2bkzNnTpo1a0bbtm258847vdq+SFrmxpnleOA9YOL5BmNMJBADPADsBtYYY+biJM6Bl6zfzlq77wrj900YSyQgdu7cydChQ/nyyy/ZsGEDAFFRUZQvXz7xrNFaS3R0dMBjK1asGG3atGHcuHEsXbqUDz74gA8++ICRI0dyyy230K5dO4oVKxbwuERCjbHWBn6jxhQF5p8/szTGVAJettbWTHjfB8Bae2mivHScmRdchjXAIGCJtfbzZPo/DjwOkD9//vJjxowha9asV4w1Li7uin1SWh6K3N4nf23fF+NeOMY///zD5MmTmT17NgDlypWjbNmylC1blptuuomMGTOmetue9vek36V94uLi+PLLL/nss8/YtGkTkZGRVKxYkYceeoi7776bqKiQ/HbmMm4ex6FyDPtzXV/3DdRncfXq1ddZa++8bIG1NuAvoCiw4YL3DXEuvZ5/3wp47wrr5wFigd+APgltXYF1Ce1PphRD+fLl7dKlS21KUurjyRihxu198tf2fTHu0qVL7fHjx+2gQYNsjhw5rDHGtmnTxv7+++8+3ban/a/2GN64caNt0qSJzZ8/vwVsvnz5bM+ePe2GDRs8CzSIuXkcB/sxHIh1fd03UJ/FwFqbRM4Ilucsk/riJNlTXmvtAWvtk9ba4jbh7NNa+461tnxCe6zfIpU069y5cyxYsIAbb7yR3r17U7VqVX744Qc++OADChcu7HZ4qXLzzTfz5JNPsmvXLubOnUvlypUZPnw4ZcqU4e677yY2NpbDhw+7HaaI64IlWe4GLvy0KQT84VIsIhex1jJv3jxuvfVW3nrrLQoVKsSyZcuYP38+ZcuWdTs8n0iXLh116tRh9uzZ7Nmzh6FDh3L8+HE6derEtddeS/Pmzdm0aZPbYYq4JliS5RrgRmNMMWNMeqApMNflmERYvXo19913H3Xr1uXMmTO8/PLLiW3hKl++fDzzzDP8+OOPrFmzhnbt2rFgwQJuv/12hg0bRnx8vNshigRcwJOlMWYKsBooZYzZbYxpb609C3QBFgGbgOnW2o2Bjk3kvM2bN/Poo49SuXJltmzZwsiRI9m4cSP33XdfmnncwhjDnXfeSUxMDJs3b6ZmzZr06NGDGjVqsGPHDrfDEwmogCdLa20za+211tp01tpC1tqxCe0LrLUlE76HHBDouEQA9u7dy5NPPsktt9zC4sWLefXVV9m6dStPPvlkYsGAtCh//vx88sknjBs3jvXr13Prrbcybty48zfciYS9YLkMK+KqI0eO8OKLL1KiRAnGjh3LU089xW+//caLL74Ydo8GpZYxhrZt2/Ljjz9yxx130L59e+rVq8dff/3ldmgifqdkKWna6dOneeeddyhevDj9+/enbt26/PLLL7zzzjvky5fP7fCCUtGiRfnyyy8ZNmwYixcvpkyZMnz88cduhyXiV0qWkmbNmTOHm266iW7dunHrrbeyZs0apkyZQvHixd0OLehFRETQvXt31q9fz/XXX0/Dhg1p3bq1HjORsOVKBR83GWPqAHUKFizYMTY2NmiqRgQTt/fJ39VPrLVMnDiR8ePHc8MNN/DEE09QoUIFj27cSW1swVTB52pju9TZs2f58MMPmTRpEnny5OH555+nfPnyqR7PV1TBx/djqIJPGnypgk/y3N4nf1Y/OXHihG3RooUFbOvWre3JkycDEtvSL7+0Nj7eefPHH9Z+8421n39u7axZ1o4fb+2IEf92njHD7mje3NoBA6wdPtzasWOtnTHj3+U7dli7ebO1e/bY5fPnW3v27FXF7Ku/97fffmtLlSplAdulSxd77Ngxn4ybWqrg4/sx0nIFn/AoAinigcOHD3P//ffz1VdfMWDAAPr06eObx0D274d334WePSFHDpgwAWJj4cgR53X0KPed/z1rVhg8GIYOvXycxx+HyEj48kuKTJ0KFz7PmDMnNGzo/P7cczBjBgD3nF9evDhs3er83qULrF/vbKtUKbKXKgX33Qd+fuSlQoUKfPfdd/Tp04fhw4ezePFiJk6cyN133+3X7YoEgpKlpAk///wzTz31FIcOHWLatGk0btz46gc9cwZiYuDll51E2K6dkyyjoiBbNihYELJnh2zZ2HnoEEXPJ6s2baB6dWfZha+IhFsIRozgf40aUa1yZYiLc14nT/673WeegXr1IC6OrT/8QIn8+eHCOSpz5XIS5dGjMGYMd5w8CZ9+CgsXXv0+pyBTpky8/fbb1K1bl7Zt21K5cmX++9//8uKLL5I+fXq/b1/EX5QsJewtWbKERo0aERERwbJly3xzprN4MXTvDps2Qc2aMGwYFC3qLGvRwnldYMeyZRTNksV5U7as87oSYyBDBueVJ8/FyypVcl7A7mXLKHHpfJavvfbv70ePsun11yldrpzz/sQJJ966daFZMyeh+0GNGjX48ccf6d69O/379+fTTz9l4sSJlClTJuWVRYKQ7oaVsBYbG8tDDz1EkSJFGDFihG8S5Zkz0LkznD4Nc+c6Z2ylS1/9uP6QLRt/1awJTZo47/fscc5Sn3sOCheGGjUo8Omnzpmxj+XIkYMPPviA2bNns3v3bsqXL8/gwYM5d+6cz7cl4m9KlhKWzp07xzPPPEOnTp2oWbMmK1eupECBAqkf8OhR6N+fiBMnIF0657Lmxo1Qp47fvwv0qRIl4NtvYfNm6NcPdu3ipsGDnffgfP964SVfH6hfvz4bNmzg4Ycf5rnnnqNGjRr8/fffPt2GiL8pWUrYOXr0KPXr1+ftt9+ma9euzJkzh+zZs6dusPh4mDQJSpWCF18k97ffOu0lSzqXSENVyZLOd61btrB21Ci4M+FO+f/+FwoUgA4d4MsvwUdngfny5WPWrFmMHz+eb775hvvvv18JU0KKkqWElV27dnHPPfewYMEC3nvvPYYPH05UVCq/mv/2W6hSBVq3di5Zrl7N3+E224gxxJUs+e/ZcYsWzs1D06bB/ffD9dfD66/7aFOGxx57jLlz57JlyxZq1KjB/v37fTK2iL8pWUrYWLNmDXfddRfbtm3j008/pXPnzlc3YO/esGMHjB8Pq1dDxYq+CDO4VavmPPry118wdSrccQdceAa4b99Vb+LBBx9k3rx5/Prrr9SoUYN9PhhTxN+ULCUszJo1i/vuu48MGTLw1VdfUatWLe8HOX0ahgxxboIBJ0lu3gyPPfbvYx1pRebMzk1Bc+c6fxNwfi9WzHlO9OzZqxr+P//5D/Pnz+e3335TwpSQkMY+ASTcWGsZNGgQjz76KLfddhvffPON948nWAvz50OZMvDss84lSIAiRZznH9O685dob78d/vMf507aChVgzZqrGvb+++9n/vz5bNu2jerVq2v2Eglqqg0bJPUIg4nb++Tp9s+cOcPQoUP57LPPqFGjBr169SLDFW66SWrczL//TvGYGPJ8+y3HCxdma5cuHLzrrquO7WrXC9rasNaSd8UKbnznHdIfOsSO1q3Z+dhjnq2bjPNVfwoUKMDQoUPJnTv3VY13nmrD+n4M1YZNgy/Vhk2e2/vkyfb//vtve99991nA9uvXz8afr7vq7bjt2lmbPbu1Q4ZYe+qUT2LzxXqe9nftGD582NqnnrJ24kTnvQd//5RiyJw5sy1durTdu3fvVY114ZhuUW3Y8KsNq8uwEnK2bNlCpUqVWL16NZMmTeKVV17xrsbrqVPw55/O74MGwa+/Qo8eoHJsnsuRwyn116qV8/7dd+HRR//9vtdL1apVY8GCBezcuZPq1auzd+9eHwYrcvWULCWkfP3111SsWJFDhw7xxRdf0LJlS+8GOHwYatVyvns7fRquuQY0ybNvLFjgVDKKiUnV85n33XcfCxcuZNeuXUqYEnSULCVk/Prrr9SuXZvcuXPzzTffULVqVa/Wz7B/P9xzD6xaBX366EzSl7p2hQ0bnMdrunRxnk/dsMHrYe69914WLlzI7t27qVatGn/88YcfghXxnpKlhIS///6bhx9+GIDPPvuMG264wbsBNmzg9s6dYedO5wzokkLn4gPFi8OiRU7Fo23b4MCBVA1zzz338Nlnn/HHH39QrVo19qTy0q6ILylZStA7efIk9erVY9euXcydO5cSJUp4N4C10KULJj4eli93LsGKfxgDLVs6xRzOVzt66y1nlhYvVK1alc8++4y9e/dSvXp1JUxxnZKlBLX4+HjatGnDV199xaRJk6hcubK3Azgf4B99xPr33oPzU1WJf2XO7Pw8dcop7lCzpnM270XxgSpVqrBo0SL+/PNPqlWrxu7du/0Tq4gHlCwlqL3wwgtMmzaNN954g0aNGnm38tChzh2aZ89CwYKcuppZRyR1MmSAdeucGU5mzICbboJx45yzfQ9UrlyZRYsW8ddff1GtWjV27drl54BFkqZkKUHr/fffZ9CgQTzxxBM899xznq8YH+88CtKzJ0RGXnVpNrlKGTPCK6/ADz84VZKeegq2bvV49UqVKrF48WL2799PtWrV+P333/0YrEjSlCwlKC1atIhOnTpRq1Yt3nvvPc+fozx1Cpo1g2HDnDs0p01zPqzFfaVLw7JlTpm8G2902jw8w6xYsSJLlizhwIEDSpjiCiVLCTq//fYbjRo1okyZMkyfPt27KbZatoTp052bSt5+2zmzlOAREQFlyzq/jx/vFGs/dcqjVe+66y6WLFnCwYMHqVatGjt37vRfnCKXULKUoLJnzx769OlD9uzZmT9/PtmyZfNugGefhcmTnZ/eVPWRwPvnH+d7zHr14Ngxj1apUKECn3/+OYcOHaJatWrs2LHDvzGKJFAh9SAp3htM3Nqn48eP061bN/bs2cM777zj8SMiWbZtI9fatexu3PiK/dwsQh02hdR9rMCCBZQaMoQjN9/MTwMHctbDODZv3syzzz5Ljhw5SO7/YxVS9/0YKqSeBl8qpJ48N/bpzJkz9qGHHrKRkZF20KBBnq/45ZdOIfTrrrP24MErdnWzCHXYFVL3pRkzrE2Xztpy5aw9ftzj1VauXGmjoqJs/fr1kyykr0Lqvh9DhdRFXGSt5emnn2bhwoWMGDGCu+++27MVp051nt8rVAhWr4ZcufwbqPhHw4Ywbx40bgyZMnm8WpUqVXjzzTf55JNPGDx4sB8DFNF3lhIEhgwZQmxsLM8//zyPP/64ZysNG+bc9VqpEqxc6UzULKGrZk2nXi/A2rWwZYtHq3Xv3p2GDRvSp08fli9f7scAJa1TshRXzZw5k+eee47GjRvz+uuve75i7tzOmciiRTqjDCfnzjnTft1zD3z/fYrdjTGMHTuW4sWL06RJE/48P/WaiI8pWYprVq9eTatWrahcuTLjx48nIiKFw9FaZ+5JgMcecy7D6hnK8BIZCZ984swIU60afPVViqtkz56dmTNn8s8//9C0aVPOqgiF+IGSpbjit99+o27duhQqVIg5c+aQyZPvqgYPdirAfPed816PhoSnUqWcS+v58sEDD8CSJSmuUrZsWWJjY/nf//5H3759AxCkpDVKlhJwBw8e5OGHHyY+Pp4FCxaQN2/elFeaORN69YLoaLjtNv8HKe66/npYsQJKlID33/doldatW/P444/zxhtvMHfuXD8HKGmNkqUE1KlTp6hfvz47duxgzpw53Hi+7NmVfP218z1W5cpO1ZeULtdKeMifH/73P5g40Xl/+nSKqwwfPpw77riD1q1ba+Jo8Sl96kjAWGtp164dK1asYMKECVStWjXllf74A+rWhYIFYc4cfUeZ1uTM6fw3P3QI7roLhg+/YveMGTMyc+ZMjDG8/PLLnDx5MjBxSthTspSA6devH5MnT+b111+nadOmnq1UoAB06gQLFoAnl2slPGXODMWLQ/fu8OqrVyzAXqxYMSZNmsSvv/7K008/HbgYJawpWUpAjBs3jv79+9OhQwd69+6d8gqnT8OePc4l11degZIl/R+kBK8MGZwZZNq0gZdecqZfu0LCrF27Ns2bN2fMmDGMHz8+YGFK+PJiOgeR1Pnmm2944oknePDBBxkxYkTK021ZCx07OndB/vyzcylOJCoKxo6F7NmdohQ5cjiJMxnt2rVj7969dOrUidtvv53bdGOYXAWdWYpfHT9+nNatW3Pttdcybdo00qVLl+I610+c6NzU8eSTSpRysYgIZ+q1oUOhQ4crdo2MjGTKlCnkypWLhg0b8s8//wQmRglLSpbiV3369GHLli2MHz+enJ4kvg8/pNj48U7RgRdf9Hd4EoqMgWeecW76Onv2is9h5s+fn+nTp7N9+3batWuHTWOzLInvaIquIJkWJpj4ap/WrVvHs88+S3R0NF27dk2xf/aff6Zct24cvPlmNg4ejPXgLNQbmqLLN7EFk0LTp1Ni5Eh+eu01Dlxyd/WF+zV9+nRGjhxJp06daJzCVG6+oCm6NEVX2Lw0RVfyfLFPhw8ftoULF7YlS5a0x44d82ylI0es7dzZrpg796q3nxRN0XV12wpKJ05YW6GCtVmzWrthw0WLLtyv+Ph426BBAxsZGWlXrFjh97A0RZem6BLxSPfu3dmzZw8TJkwgc+bMV+584AAcOwbZssF773E2W7bABCmhL2NGmD0bsmaFevXg4MEkuxljGDduHMWKFaNx48b89ddfAQ5UQp2Spfjc3LlzGT9+PH369KFixYpX7nzyJNSpAw8/fMVHAUSSVbAgzJoFu3ZB27bJdsuRIwcff/wxhw4dolmzZiq4Ll5RshSf2r9/Px07dqRcuXL069fvyp3j450beVavhq5dVRhdUq9SJacUYgo3hd16662MHDmSpUuXpnx8ilxAz1mKz1hr6dSpE4cPH+bzzz8nffr0V16hb1+YPh3efBMefTQwQUr4atbs39937Ei2W5s2bVi1ahUDBw6kcuXK1K5d2/+xScjTmaX4zOTJk/n444959dVXKVu27JU7T5gAAwfCE0/As88GJkBJG0aNgtKlybZ5c7Jd3n33XW6//XZatWrF9u3bAxichColS/GJPXv20KVLFypXrsyzniS/e+5xar6+954uv4pvNWgA+fNzy4svwp9/JtnlfMF1gIYNG6rguqRIyVKumrWW9u3bc/r0aSZMmEBkZGTynffudW7kueEGGDHCKWEm4kvXXAOffEK6o0edy/unTiXZ7YYbbmDChAmsX7+ebt26BThICTVKlnLVRo0axaJFi3jrrbcoUaJE8h337oW774YePQIXnKRN5crxS69e8NVXcIWZR+rWrcvzzz/P6NGjmTx5cgADlFCjf9bLVdm6dSs9e/bkgQceoFOnTsl3PHbMeUTk4EFnImcRP9tfvTpERkIKE4z379+f5cuX061bN2rWrEmePHkCFKGEEp1ZSqqdO3eONm3akC5dOsaNG5f8bCLnzkHz5vDddzB1KtxxR2ADlbSrX79/75I9cSLJLlFRUYwaNYrDhw/Tq1evAAYnoUTJUlJtyJAhrFq1ivfee49ChQol37F/f5g715ktQrfpixs+/tiZE3XnziQXly1blp49ezJu3DiWL18e4OAkFChZSqr89NNPvPjiizRo0IAWLVpcuXOVKs4M9126BCQ2kcuUKQNHj0L9+nD8eJJd+vXrR9GiRXnyySc5ffp0YOOToKdkKV47ffo0rVu3JmfOnMTGxiZ/+fV8ObH//MeZrFePiIhbSpWCyZPhhx+gXbskSytmzpyZmJgYNm3axFtvveVCkBLMlCzFa6+99hrff/89o0eP5pprrkm605kzcP/9oA8dCRYPP+wUwpg2Dd54I5kuD9OwYUP69+/P1q1bAxygBDMlS/HKt99+y8CBA3nssceoV69e8h3/+19YvhwKFw5ccCIp6dULmjZ1LskmY/jw4aRLl46nnnpKk0VLIiVL8diJEydo3bo11113HcOHD0++4yefwODB8NRTzgeTSLAwBj76CAYMSLbLddddx+uvv86SJUuYOnVqAIOTYGbS2r+cjDF1gDoFCxbsGBsbGzSzcweT5Pbpvffe4+OPP2bw4MGUL18+yXUz/vEHdz7+OMcLFeK7d97BplRM3YvtXy03Z5n3dj1P+wfTDPPBJqX9yv7zzxR7/302vPYa5y7pd+7cOTp37sy+ffuYOHGi13+fcDyGvV3X130DdRxXr159nbX2zssWJDUjdFp4lS9fPqhm5w4mSe3Tl19+aQHbpUuXK688ZYq1+fJZu22bT7fvC27OMu/tep721zGcvBT363//szYqytpHHrH27NnLFq9bt85GRETYJ5980vfbTiU3j2Fv1/V130Adx8Bam0TO0GVYSdGRI0do27YtJUuW5I1kboxI1LQp/PYbFCsWmOBEUuvee+Gdd+DTT5OcB/OOO+6ga9eujBo1itWrV7sQoAQTJUtJ0TPPPMOuXbuYMGECmTNnTrrT9OnOd5UAYXhJT8LUk0/C4487d8l+/PFli1999VUKFizIE088wZkzZ1wIUIKFkqVc0bx58xg3bhy9e/emYsWKSXf6+Wdo29Z5ljI+PrABilwNY+Ddd6FChSSTZbZs2Xj33Xf56aefePvttwMfnwQNJUtJ1t9//03Hjh257bbbeOmll5LuFBcHDRs6Z5NTpkCEDikJMenTw2efOXfJJqF+/frUrVuXl19+mZ3JlMuT8KdPNklWly5dOHjwIBMnTiR9Une1WgtPPAGbNzvVUa67LvBBivhC7tzOWeauXbB27WWL3333XYwxdOnSRc9eplFKlpKkZcuWMW3aNPr27cutt96adKelS50k+corTrUekVBmrVM7tkmTy+rHFilShFdffZX58+cze/Zsd+ITVylZymXOnTtHjx49KFKkCM8991zyHatXh/nznWo9IqHOGKeYxrZtzj8AL9G1a1fKlSvH008/zZEjR1wIUNykZCmXWbJkCd999x2DBg0iU6ZMl3c4dMi5qccYeOQRfU8p4aN6dWjfHoYMceZfvcD5eS/37t3Li0k8aiLhTZ9ycpG4uDjGjBlDxYoVaZpUqTprnTtfq1SBw4cDHp+I3731FuTNCx07/jtzToK77rqLTp068d5777Fu3TqXAhQ3KFnKRd58800OHDjA0KFDk556a8gQmDMHXnoJcuYMeHwifpcrl/M4SblycOrUZYtff/118uXLxxNPPMG5c+cCH5+4QslSEu3atYvBgwdTo0YNKlWqdHmHlSuhd2949FHo1i3wAYoESqNGMGYMZMly2aIcOXIwfPhw1q1bR0xMjAvBiRuULCXRf//7X6y1dOzY8fKF+/c7dwkWKwZjx2oiZ0kb1q+HZ565bLLoRo0aUatWLV544QV2797tUnASSEqWAjjzVH744Yf06NGDAgUKXN4hRw4nWc6Y4fwukhZ89RW8/fZlBQuMMcTExHD27Fm66SpLmqBkKVhr6dGjB/nz56d3796Xdzh92qlyMnSo8z2OSFrRqRNUrOicXf7990WLbrjhBvr168esWbOYP3++SwFKoChZCjNnzmTVqlX079+fbNmyXbxw8WIoXRq2bHEnOBE3RUbC++/DP/9Ajx6XLe7Zsye33HILnTt35tixYy4EKIGiZJnGnTx5kl69enHrrbfStm3bixfu3g0tWjg3ORQq5E6AIm4rUwaefx4mTYJlyy5alD59ekaNGsXvv//OK0kUMpDwEeV2AOKu4cOHs2PHDj7//HMiIyP/XXDunDM35cmTzveUyU3NJZIWvPACXHMNVK582aIqVarQoUMHhg4dSosWLbjttttcCFD8TWeWadi+ffsYMGAAderU4f5La7sOGwarVsHIkVCqlDsBigSLjBmha1fnu/sk5rV84403yJ07N08++STxmqYuLClZpmH9+vXjxIkTvPXWWxcvsBY+/9wpKt2ihSuxiQSl9evhxhvhkuo9uXPnZsiQIXz99deMHj3apeDEn0xam27GGFMHqFOwYMGOsbGxZM2a9Yr94+LirtgnpeXBatu2bXTs2JH69evz9NNPX7QsLi6OrJkzE3nyJOdcuPzqr7+pL8ZN7Rjerudpf0/6hesxnBJ/7FdUXBwV2rThdO7crB85EnvBVxfWWnr27MmWLVsYOXIkhQsX9um2wd1j2Nt1fd03UMdx9erV11lr77xsgbU2Tb7Kly9vly5dalOSUh9Pxgg28fHx9oEHHrC5cuWyBw4cuHjhxx/bVTNnuhNYAn/9TX0xbmrH8HY9T/un1WPYE37br5kzrQVr33zzskW//PKLjYyMtA0aNPDLpt08hr1d19d9A3UcA2ttEjlDl2HToIULF7JkyRJeeuklcufO/e+CDRugWTOKjRnjXnAiwa5BA6hXz6mPvG3bRYtKlSpF27ZtmTdvHr///rtLAYo/KFmmMWfOnKFnz56ULFmSp5566sIF8NhjkCMH2554wr0ARYKdMRATA1FRTunHS/Tr1w+AV199NdCRiR8pWaYxo0aN4pdffuGtt94iXbp0/y4YNMi5eSE2ljOaTUTkygoWdG7y6d//skWFCxembt26jB8/ni0q5hE2lCzTkEOHDvHyyy9To0YN6tSp8++CH36AV1+F5s2dS0wikrIbb3TOMn//3Zlo4ALNmzcnY8aMvPTSSy4FJ76mZJmG9O/fn4MHD14+V2WhQs5Et++8415wIqEoLg5uvx26d7+oOXfu3HTv3p2pU6fyww8/uBOb+JSSZRrx66+/8u6779KuXbuLK4xYC3nywIgRzk8R8VzWrE6xgsmTYeHCixY9++yz5MyZk759+7oUnPiSkmUa8fzzz5MhQwb6X/gdy7p1UKUKbN/uXmAioa53b2eygU6dnDPNBDlz5qRXr17Mnz+f1atXuxig+IKSZRqwbNkyZs+eTZ8+ff6dq/LUKefu1507IVcudwMUCWUZMjgzk+zcCS++eNGirl27ki9fvsSJ1SV0KVmGuXPnztGjRw+KFCnCM8888++CV16BjRthzBjQ3a8iV6dKFXjqKThxwvlqI0GWLFno27cvy5Yt44svvnAxQLlaSpZhbuLEiXz33XcMGjSITJkyOY3ffgtvvAHt28NDD7kboEi4ePddiI117pC9wOOPP06RIkV0dhnilCzDWFxcHC+88AIVK1akadOm/y4YNMh5TmzIEPeCEwk3EQkfp99+S55VqxKbM2TIwEsvvcSaNWuYO3euS8HJ1VKyDGNvvvkme/fuZdiwYRc/KjJ5MixaBDlyuBecSLjq1YtSQ4ZcdLNP69atKVmyJH379uXcuXMuBieppWQZpnbt2sXgwYNp2rQpFStWdBo3b4ajR525+UqXdjdAkXA1aBDpDx1y5oRNEBUVxauvvsqGDRuYOnWqi8FJailZhqnz348MGjTIaTh+HOrUceaoFBH/qViR/ffcA2++eVFln0aNGnHbbbfx0ksvcSaJCaQluClZhqFvv/2WDz/8kB49enD99dc7jf/9L/z6K+gBaRG/296hg/MP1AEDEtsiIiIYMGAAv/32Gx988IGL0UlqKFmGGWstPXr0IH/+/PTu3dtp/N//YPhw6NIFqld3N0CRNOB4kSJOCbxrr72o/eGHH6ZSpUq8+uqrnDhxwp3gJFWULMPMzJkzWbVqFf379ydbtmzOTQbt2kHx4s5dsCISGEOGwPPPX9RkjOH1119nz549jBw50qXAJDWULMNIfHw8L7/8MjfffDNt27Z1Go8eheuvhw8+gCxZ3A1QJK2xFmbMgJ9+SmyqVq0aDzzwAAMHDuTo0aMuBifeULIMI3PmzOHnn3+mb9++REZGOo3XXgtffAH33ONucCJp0ZEj8MQT0KvXRc39+/fn77//5u2333YnLvGakmWYsNYyYMAASpQoQePGjZ3/STt2hL17L6soIiIBkiMH9OkDn30Gy5YlNt91113Ur1+fwYMHc/DgQffiE48pWYaJRYsWsW7dOnr37u2cVT73HIwbBzt2uB2aSNrWpYszZ+zzz19UN/a1117j6NGjvPnmmy4GJ55SsgwTAwYMoHDhwrRq1QoWL4bRo+HZZ6FSJbdDE0nbMmVyJi749luYNSuxuUyZMjRv3px33nmHvXv3uhigeELJMgwsX76clStX0qtXL9KfOOEUSC9d2vkfVETc17o1PPwwpE9/UfPLL7/MmTNnGHDB85gSnJQsw0D//v3Jnz8/7du3h5dfdr6nnDDBKWsnIu6LioJPP3WqaF2gRIkStG/fntGjR7NDX5kENSXLELdmzRqWLFlCjx49nCm4XngBpk6FChXcDk1ELnXyJLzzjlPdJ0Hfvn2JiIjgFV0JCmpKliFuwIAB5MqVi06dOjk3D+TNCw0buh2WiCRl/Xro1s2pqJWgUKFCdO7cmYkTJ7Jp0yYXg5MrUbIMYT/99BNz5syhW7duZNu2zbmZZ/Nmt8MSkeRUruxcin3jDThwILG5d+/eZM6cmX79+rkYnFyJSWszdxtj6gB1ChYs2DE2NpasWbNesX9cXNwV+6S03J9ee+01Vq9ezdSpU6kycCDZN2zgm8mTOZst21WN6+Y++XP7vhg3tWN4u56n/T3pF8zHsD+5uV9X2naW7du5s0MHdjdsyG+dOiW2jxs3jkmTJjFq1ChKlizp9bi+iM2X6/q6b6CO4+rVq6+z1t552QJrbZp8lS9f3i5dutSmJKU+nozhD5s3b7YRERH2+eeft3b5cmvB2kGDfDK2W/vk7+37YtzUjuHtep72D+Vj2N/c3K8Ut92mjbUZMli7c2di0+HDh22uXLnsQw89lPpxfRGbj9b1dd9AHcfAWptEztBl2BA1aNAg0qdPzzPduzsVQq69Fp5+2u2wRMQTr7zifG1yQW3YHDly0Lt3bxYuXMjKlStdDE6SomQZgnbu3MmkSZPo2LEj+b//HlatgpdegsyZ3Q5NRDxRpAgsXQq33HJRc5cuXShQoEDi5O0SPJQsQ9Bbb72FMYbnnnvOmZ9y9GhnGi4RCS1//QWxsYlvM2fOTN++fVmxYgWLFy92MTC5lJJliNm7dy9jxozhscceo3DhwpAhg1MwPV06t0MTEW+NHw+dOsEFl107duxI0aJFdXYZZJQsQ8zQoUM5c+YMzz/zDFSrBnPnuh2SiKTW00879xtcUGQ9ffr0vPzyy6xfv55ZF9SSFXcpWYaQAwcOMHLkSJo2bUqJpUvhf/+7rNakiISQzJmdEpVffQXz5iU2t2zZkptuuol+/foRHx/vXnySSMkyhAwfPpxjx47xQvfu8NprcN99ULOm22GJyNVo1w5KlnTuaj93DoDIyEj69u3Lzz//zMKFC10OUEDJMmT8888/vPvuu0RHR3PzkiXOjQEDB2piZ5FQFxXl/L98220XPUrSuHFjChUqxODBg10MTs5TsgwRI0aM4PDhw/R7+mmnVFa9epqrUiRcNGgAkydDzpyJTenSpaNr164sW7aM9evXuxebAEqWIeH48eMMGzaMWrVqUa5aNZg40fmXqIiEl59+gtmzE98+/vjjZMuWjSFDhrgYlICSZUh4//332b9/Py+88IJz2bVePWdyZxEJL336OJO3HzoEOFV9OnTowLRp09i1a5fLwaVtSpZB7tSpU7z55pvce++9VJ01C/r3dzskEfGXAQPg8GHnq5YE3bp1A+Cdd95xKSgBJcugN2HCBP744w9ef+wxZ9LY/fvdDklE/OW226BFC2e+yz17ALj++utp2LAho0eP5tixYy4HmHYpWQaxs2fPMmjQICpUqEDlzz6DjBnhhRfcDktE/OnVV51HSF5+ObGpZ8+eHDlyhE8//dS9uNI4JcsgNmXKFLZv385bTZpgZsyAnj0hXz63wxIRfypWzPl/vUCBxKo+FSpU4N5772XWrFmcPXvW5QDTJiXLIBUfH8/AgQMpW7Ys9y5aBHnyOP8DiUj4GzjQKTxywXPUPXv25K+//mLmzJkuBpZ2RbkdgCRt9uzZbNq0iSlTpmBKlYKdOyF7drfDEpFAsRY++wyKF4eSJalduzaFCxdm8ODBNGnSBKOCJAGVqjNLY0wWY0ykr4MRh7WWAQMGcOONN9KoUSO4/XaoX9/tsEQkkA4dgoYNE5+pjoiIoGHDhqxbt47ly5e7HFza41GyNMZEGGOaG2M+NcbsA34B9hpjNhpj3jLG3OjfMNOWhQsX8t133zHywQeJbNcO/vnH7ZBEJNBy53aeufzoI9i9G4AHH3yQvHnzqkiBCzw9s1wKFAf6AAWstYWttfmAe4CvgUHGmJZ+ijFNsdbSv39/ihUuTI3PP4dvv4UsWdwOS0Tc0KMHxMfD228DkDFjRp566inmzZvH5s2b3Y0tjfE0Wf7HWvuatfZHa23ifDHW2oPW2o+ttY8C0/wTYtqybNkyVq9ezZh778Vs3uw8pBylr5ZF0qSiRaFJExg1yilWAHTu3JkMGTIwbNgwV0NLazxNlsOMMZWv1MFae8YH8aR5AwYMoEi+fFRbtgwqVIDoaLdDEhE3Pfcc5MoFW7YAkC9fPlq1asWECRPYryIlAeNpsvwVGGKM2WGMecMYU86PMaVZX3/9NV988QXj776biD17YNAgTcElktaVKwfbt8NddyU29ejRg5MnTzJixAj34kpjPEqW1trh1tpKwH3AQeADY8wmY0w/Y0xJv0aYhgwYMIDcuXNz1+DBMGwY1KjhdkgiEgwiI+H0aTIl3OhTunRpHnnkEWJiYjhx4oTLwaUNXj06Yq3daa19w1p7O9AciAY2+SWyNOaHH35g/vz5dOvWjSwlS0L37m6HJCLBpEEDyv73v04pPJwiBfv37+fDDz90ObC0watkaYxJZ4ypY4z5CFgIbAEe9Utkaczrr79OsSxZ6LNyJWzc6HY4IhJsWrUi865dMHcuANWqVeP2229n6NChxMfHp7CyXC1Pn7N8wBgzDtgNPA4sAIpba5tYaz/xY3xpwpYtW5gxYwaTbrqJdF9+qbtfReRyjz7KiWuvdabvshZjDD179uSXX35h4cKFbkcX9jw9s/wvsBooba2tY639yFqruWJ8JCYmhuKRkVT+6Sdo2xZKlXI7JBEJNlFR7GrUCL75BlauBKBx48YUKlSIwYMHuxxc+PP0Bp/q1tr3gUPGmJbGmH4Axpgixpi7UlhdriAuLo7x48cztnBhp9bjSy+5HZKIBKk/H3oI8uaFGTMASJcuHd26dWPZsmWsX7/e5ejCm7e1YUcAlYBmCe+PAjE+jSiN+fDDD7n2yBHu2bEDnn4aChVyOyQRCVLxGTM6Vb2GD09s69ixI9myZVMJPD/zNlneba3tDJwEsNYeAtL7PKo0wlpLTEwM+cqWhTffdB4+FhG5kmLFnOevT50CIEeOHHTo0IFp06axa9cul4MLX94myzMJs41YAGPMNYBuw0ql5cuXs2HDBh7r1g3z7LOa2FlEPDNnDhQsmFhgvVu3bgC88847bkYV1rxNlu8As4F8xpgBwErgdZ9HlUbExMTQOksWWp496xRLFhHxxG23ObViEy7HXn/99TRq1IjRo0dz5MgRd2MLU54+OhIFYK39COgFDAT2AvWttTP8F1742rNnD7M+/pg3M2Qgw9ixEJGqqUVFJC1KosB6z549OXLkCGPGjHE1tHDl6Sf0t+d/sdb+Yq2Nsda+Z61V9Z5UGj16NPfFx5P/4EHo3NntcEQk1Dz3HBw9CrGxANx5553ce++9DB8+nLNnz7ocXPjxNFmqmrcPnT59mtGjR/Na/vyQJ4/zL0QREW+UKwcPPuhcij3jTPrUs2dPfv/9d2bOnOlubGHI01Ix1xhjeiS30Fo71EfxpAmzZs0i3Z9/UikiwvnXYcaMbockIqHozTfBWkiXDoDatWtTsmRJBg8eTJMmTZxnt8UnjLU25U7G7AVGkswZprX2FR/H5TfGmDpAnYIFC3aMjY0la9asV+wfFxd3xT4pLU9K165dKbh3L9Ny5+bnV17hZIECXq3vb6nZp1DYvi/GTe0Y3q7naX9P+vnjGA4Fbu6Xm8fw3LlzGTZsGMOGDaNcuXI+jc2bdX3dN1DHcfXq1ddZa++8bIG1NsUXsN6TfqH0Kl++vF26dKlNSUp9PBnjQt9//70F7ODBg71aL5C83adQ2b4vxk3tGN6u52l/N47hUOHmfgX0GD540NpWraydM8daa+3x48dt3rx5bZ06dXwemzfr+rpvoI5jYK1NImfoO8sAi4mJ4db06Wn3qCZrEREfyJYNVq2CgQPBWjJlysRTTz3FvHnz2Lx5s9vRhQ1Pk+X9fo0ijTh06BAfffQRs3LkIJeSpYj4QlQU9OwJX3+dWGC9c+fOZMiQgWHDhrkcXPjwtJD6QX8HkhaMHz+em44fp/j+/fDYY26HIyLhok0bp8D6m28CkC9fPlq1asWECRPYv3+/u7GFCT0JHyDx8fGMGDGCV/LlgyxZlCxFxHcyZ3YmYpg/P3Hy+B49enDy5ElGjBjhcnDhwatkaYxpZIzJlvB7X2PMLGPMHf4JLbwsXryYg1u38tChQ9CqFeTI4XZIIhJOOneGHj0gd24ASpcuzSOPPEJMTAwnTpxwObjQ5+2Z5YvW2qPGmKpATWACziMlkoKYmBiaZs9O5JkzqtgjIr6XJw8MGQLXXpvY1LNnT/bv38+HH37oYmDhwdtkeS7h5yPASGvtHDRFV4q2b9/Op59+Su6uXWHrVihTxu2QRCRcffkljBsHQLVq1bj99tsZOnQo8Zqs4ap4myz3GGNGAY2BBcaYDKkYI80ZOXIkkcbwxBNPQPHibocjIuFszBjo3h0OH8YYQ8+ePfnll19YuHCh25GFNG8TXWNgEVDLWnsYyA1oxuIrOHHiBGPHjuWr/PkpNFRVAUXEz84XWB81CoDGjRtTsGBB3n33XZcDC21eJUtr7XFr7Sxr7a8J7/daaxf7J7TwMHXqVHIfPEiFvXshZ063wxGRcHf77fDAA/D223DqFOnSpaNDhw4sXryYHTt2uB1dyNIlVD+y1vLee+/xYu7c2HTp4PHH3Q5JRNKC55+HP/+EhBt72rVrB8DYsWPdjCqkKVn60TfffMPm9etpcuIEpmFDCLKC6SISpmrUgEceSZxUvkiRItSqVYtx48Zx7ty5FFaWpChZ+lFMTAztM2Qgw4kT0KWL2+GISFphjFOgoG3bxKaOHTvyxx9/8M0337gYWOhKMVkaYx4wxrxvjCmX8F7XEj2wb98+pk+fTo4mTWDwYKhUye2QRCStOXsWFiwAa6lduzYFChTg008/dTuqkOTJmeVTOHe8tjTG1ADK+TWiMDFmzBhOnz5N8z59nCLHmoRVRAJt/HjncuyqVaRLl462bdvy9ddfs3v3brcjCzmeJMv91trD1tpngQeBCn6OKeSdPXuW2NhYRpUowU179rgdjoikVc2bO5V9Egqst2/fnvj4eD744AOXAws9niTLxHN2a21vYKL/wgkP8+bN4+yuXXTYvt25BCIi4obzBdbnzYOff6Z48eLccccdjB07VhV9vJRiskwoaXfhez3ZmoKYmBiey54dEx8PnTq5HY6IpGWdO0OmTDB8OAC1a9dm586dLFmyxOXAQou3s47caYyZbYxZb4z50RjzkzHmR38FF4o2bdrE8i++oGN8POahh6BECbdDEpG0LG9eaNgQvvsOrKVKlSrkzZuX0aNHux1ZSInysv9HODf7/AToHD4JI0aMoFFkJFnj4jS7iIgEhxEjnHl0jSF9+vQ89thjDB8+nL/++ov8+fO7HV1I8PY5y/3W2rnW2u3W2p3nX36JLAQdPXqUCRMmULViRahWDWrVcjskERHImtW5Iz+hIEGHDh04e/Ys48ePdzeuEOJtsnzJGDPGGNPMGNPg/MsvkYWgSZMmcfToUe4YMgSWLk2sniEi4roZM+C664j65x9uuukm7r33XsaMGaMbfTzk7ad5W5znLGsBdRJetX0cU0iy1hITE0PLm27irnLl3A5HRORiJUrAvn3k+9//AKeiz9atW1m2bJm7cYUIb5PlbdbaO621j1lr2ya82vklshCzbNky/vj5Z8Zt24Z5/nm3wxERuVi5cnDLLeRf7EwU9eijj5IzZ07ef/99d+MKEd4my6+NMTf7JZIQFxMTw1OZM5Pu9Gl47DG3wxERuZgx0LIlOTZuhG3byJQpE61atWLWrFn8/fffbkcX9LxNllWB740xm/XoyL92797NnNmzeSZ9eqhc2ZlPTkQk2LRo4fxMmLqrY8eOnD59mkmTJrkYVGjwNlnWAm7EKXt3/vvKOr4OKtSMGjWK/8THk/fwYc0uIiLBq3BhtnbuDHWcj+2yZctSsWJF3n//fay1LgcX3LxKlhc+LqJHRxynT59m9OjRPH399ZA/Pzz6qNshiYgka3fDhhdd/erYsSObNm1i1apVLkYV/Lyt4DPBGJPzgve5jDHjfB5VCFm+fDn79u0jMjYWvv4a0qd3OyQRkStbvRqmTwegSZMmZMuWTTf6pMDby7C3WmsPn39jrT0EpOkv6D755BNuLF6cBx58EIoWdTscEZGUDRvmfGV05gxZsmShRYsWzJgxg8OHD7sdWdDyNllGGGNynX9jjMmN9yXzwsZ3333Hbxs38u0//xCR8IW5iEjQa9kS9u+HhMdIOnbsyIkTJ/joo49cDix4eZsshwBfGWNeM8a8CnwFvOn7sEJDTEwMraKiyPn331CkiNvhiIh4plYtZ57LhLtg77jjDu644w5Gjx6tG32S4e0NPhOBR4G/gP1AA2ttmrzn+MSJE0yfNo1eWbLALbfAffe5HZKIiGfSp4cmTWDOHDhyBHDOLn/88UfWrFnjcnDByevipdban62171lr37XW/uyPoEJBpkyZ2PrRR5T45x9ndhFj3A5JRMRzLVs681xu3AhA8+bNyZw5s270SYYqfV+FfNOnczZLFmjVyu1QRES8U7Ei7N0LlSoBkD17dpo2bcqUKVM4evSoy8EFH4+SpTGmkjE6dbpM9+5s6d7dmf5GRCSUGAMZMoC1cPo04FyKPXbsGFOmTHE5uODj6ZnlY8A6Y8xUY0wbY0wBfwYVMu68k33/+Y/bUYiIpM7Ro1CyJAwfDsDdd99NmTJldCk2CR4lS2vtk9baO4CXgVzAeGPMamPM68aYe40xkf4MUkRE/CBbNsibN/GuWGMMHTt2ZO3atXz//ffuxhZkvL0b9hdr7TBrbS2gBrASaAR844/gRETEz1q1gp9+gh+dOTFatmxJxowZdXZ5iVTf4GOtPWGtXQCst9be6cOYREQkUBo3hqioxJlIcufOTcOGDfnwww85duyYy8EFD1/cDfuKD8YQERE35M0LDz0EkyfDuXOAc6PPkSNHmDFjhsvBBQ+PStVdYc5KA+T3XTgiIhJwzz8P+/Y5d8YC99xzD6VKleL999+nTZs27sYWJDyt65ofqAkcuqTd4JS8ExGRUFWlykVvz9/o8+yzz7Jx40ZuueUWlwILHp5ehp0PZE1iLssdwDK/RSciIoGxZw/07w/HjwPQunVr0qVLpxt9Enj66Eh7a+3KZJY1921IIiIScL/+Ci++CHPnAnDNNdcQHR3NpEmTOHnypMvBuU/l7kREBO69FwoXTnzmEpwbfQ4ePMisWbNcDCw4KFmKiAhERECLFrBokXOzD1CjRg1uuOEGXYpFyVJERM5r2dJ5fGTaNAAiIiLo0KEDy5YtY8uWLS4H5y4lSxERcdxyC9x1l3OzT4I2bdoQGRnJmDFjXAzMfZ4+OiIiImnBV19B5L/lvq+99lrq1KnD+PHjeeCBB1wMzF06sxQRkX+dT5RxcYlNjz/+OPv37+err9LuY/VKliIicrFnn4Xbbkus6PPggw9SpEgR5s+f73Jg7lGyFBGRi912G2zbBqtWARAZGUm7du1Yu3Yt27dvdzk4dyhZiojIxaKjIXPmxJlIANq1a0dERARjx451MTD3KFmKiMjFsmZ1Eub06XDqFACFCxemQoUKfPDBB5w9e9blAANPyVJERC7XsiUcOgQLFyY21a5dmz/++IPPPvvMxcDcoWQpIiKX+89/IDYWqlZNbKpYsSJ58uTho48+cjEwdyhZiojI5aKi4IknnMmhE5uiaNSoEXPmzCHugkdL0gIlSxERSdq5czBmDFzwyEiLFi04ceIEc+bMcTGwwFOyFBGRpEVEwJAh8OabiU2VK1emSJEiTJ482cXAAk/JUkREkmYMtGoFK1bAjh2AU1y9WbNmLFq0iP3797sbXwApWYqISPKaN3d+XnBTT4sWLTh37hwzZsxwKajAC8lkaYy5wRgz1hgz84K20saYWGPMTGNMJzfjExEJG0WLOhNDf/hhYvm7smXLUqZMmTR1KTbgydIYM84Ys88Ys+GS9lrGmM3GmK3GmN5XGsNau81a2/6Stk3W2ieBxsCdvo9cRCSNatkSsmcn6ujRxKbmzZuzatUqdiRcng13bpxZjgdqXdhgjIkEYoCHgJuBZsaYm40xZY0x8y955UtuYGNMXWAl8IX/whcRSWM6dIBvvuFs9uyJTc2aNQNgypQpbkUVUAFPltba5cDBS5rvArYmnDGeBqYC9ay1P1lra1/y2neFsedaaysDLfy3ByIiaYwxAETGxUFCqbuiRYtSpUqVNHMp1tiEa9AB3agxRYH51toyCe8bArWstR0S3rcC7rbWdklm/TzAAOABYIy1dqAxphrQAMgA/GitjUlivceBxwHy589ffsyYMWTNmvWKscbFxV2xT0rLQ5Hb++Sv7fti3NSO4e16nvb3pF9aPIbB3f0Kx2M42y+/UK5rVza+9hoH774bgE8++YThw4czduxYbrjhhlRvJ5iO4+rVq6+z1l7+VZ61NuAvoCiw4YL3jXCS3vn3rYB3/RlD+fLl7dKlS21KUurjyRihxu198tf2fTFuasfwdj1P++sYTp6b+xWOx7A9dcqezp7d2ubNE5v27dtnIyMj7fPPP39V2wmm4xhYa5PIGcFyN+xuoPAF7wsBf7gUi4iIXCp9evZVqwazZ0PCjT7XXHMNDz74IFOmTCE+Pt7d+PwsWJLlGuBGY0wxY0x6oCkw1+WYRETkAn/95z9w4gR88kliW4sWLfj999/56quv3AssANx4dGQKsBooZYzZbYxpb609C3QBFgGbgOnW2o2Bjk1ERJJ35JZboEgRmDo1sa1evXpkypQp7G/0iQr0Bq21zZJpXwAsCHA4IiLiqYgIeP99KPzvt2ZZs2alXr16TJ8+neHDh5MuXToXA/SfYLkMKyIioeDBB6F06YuamjdvzoEDB1i8eLFLQfmfkqWIiHhn2TIYMCDxbc2aNcmdO3dYX4pVshQREe8sXQovvgh79wKQPn16GjVqxCeffMKxY8dcDs4/lCxFRMQ7TZo4RdVnJs5lQfPmzTl+/HjYTgrtSgUfNxlj6gB1ChYs2DE2NjZoqkYEE7f3KRyrn6iCT+Cpgo/vx7hw3Tvbt+dcpkx89957AMTHx9OsWTNuuOEGBg4cqAo+4fJSBZ/kub1P4Vj9RBV8Ak8VfHw/xkXrvv66tWDtjh2JTb169bJRUVF2//79quAjIiJCkyZw/fWwfXtiU/PmzTl79iwzL7g8Gy6ULEVExHs33OAkymrVEptuvfVWbr75Zj766CP34vITJUsREUkdY5wpuxLugDXG0KJFC1auXMmff/7pcnC+pWQpIiKpc+yYU81nyJDEpvOTQi9dutStqPxCyVJERFInSxYoVQqmTHEeJQGKFStGpUqV+Pzzz10OzreULEVEJPWaNoVffoGffkpsat68Odu2beOnC9pCnZKliIik3qOPQmSkc3aZoHHjxkRERDDlgrZQp2QpIiKpd8018MADzrRdCZdi8+XLx5133snkyZOxYVL4RslSRESuzosvwtixickS4P7772fnzp1hMym0kqWIiFydypWhRg1nvssEVatWDatJoVUbNkjqEQYTt/cpHOtqqjZs4Kk2rO/HuNK6mX//nXyff86Oxx6DyEji4uIYOnQo69evZ+bMmURFRV1VTKoNq9qwQcftfQrHupqqDRt4qg3r+zGuuO7UqU6t2IQ+S5cutXPmzLGAXbBgwVXHpNqwIiIS+mrXdp67nDo1salWrVrkypUrLMrfKVmKiMjVy5IF6tZ15rg8cwYIr0mhlSxFRMQ3mjaFAwfgiy8Sm5o3b86xY8eYN2+ei4FdPSVLERHxjZo1oUgR2LUrsemee+6hYMGCIX9XrJKliIj4RoYMzrRdHTsmNkVERNCsWTMWLlzIgQMHXAzu6ihZioiI70REgLVEnDyZ2NSiRYuQnxRayVJERHzr/vu56Y03Et/edtttlC5dOqQvxSpZioiIb5UuTZ7VqyEuDnAmhW7evDnLly9n1wXfZ4YSJUsREfGtpk2JPHUK5s5NbDo/KXSozkSiZCkiIr5VpQonr7nmogIFxYsXp2LFiiF7KVbJUkREfCsigv3Vq8Nnn8GhQ4nNzZs354cffmDjxo0uBpc6KqQeJMV7g4nb+xSORahVSD3wVEjd92N4s67dtIn8e/bw9z33EJ8hAwAHDx6kUaNGNG/enPbt23s1rgqpq5B60HF7n8KxCLUKqQeeCqn7fgxv1k2u74MPPmiLFStm4+PjvRpXhdRFRCQ87d8Pb70Ff/2V2NSiRQu2b9/O119/7WJg3lOyFBER/9i3D3r1coqrJ6hfvz4ZM2YMuRt9lCxFRMQ/brkFypS56K7Y7NmzU6dOHaZNm8aZhNlJQoGSpYiI+E/TprBy5UXF1Vu0aMH+/fv54oLZSYKdkqWIiPhPkybOz+nTE5tq1apFzpw5Q+pSrJKliIj4T4kScPfd8PvviU0ZMmSgYcOGzJ49mxMnTrgYnOeULEVExL9WrIDhwy9qatiwIXFxcXz++ecuBeUdJUsREfGvdOmcn6dPJzZVr16dHDlyMHv2bJeC8o6SpYiI+F/PnlChQuLb9OnTU7t2bebOncu5c+dcDMwzSpYiIuJ/xYvDjz/Chg2JTQ0aNODAgQP8+OOPLgbmGSVLERHxv4YNISLiomcua9asScaMGVmxYoWLgXlGyVJERPwvXz64/34nWSZM4JElSxZq1arFypUrsUE+qYeSpYiIBEbTpvDbb7BuXWJTdHQ0+/fvZ+3atS4GljJN0RUk08IEE7f3KRynN9IUXYGnKbp8P4Y36ybVN+roUQosXMhfDz7ImZw5AThy5AgNGjSgSZMmdOzYMdXb1hRdmqIr4Nzep3Cc3khTdAWepujy/Ri+mKIrKeXLl7clS5a8aNoub8fTFF0iIhI+Tp6ESZOcO2MTVK1alS1btrBp0yYXA7syJUsREQmcs2fhiSdg1KjEpqpVqwIEdYECJUsREQmcrFmhTh2YMcNJnEDevHmpWLEis2bNcjm45ClZiohIYDVtCvv3w9KliU0NGjRg/fr17Ny508XAkqdkKSIigfXQQ5A9+0UFCqKjo4HgvRSrZCkiIoGVMSNER8OvvyYWKChRogRly5ZVshQREUk0ahQsXw7GJDZFR0ezYsUK9u3b52JgSVOyFBGRwMuQwfl5QWGcBg0aYK1l7ty5LgWVPCVLERFxx/DhcNNNEB8PwK233kqxYsWC8q5YJUsREXFHvnywZQvZf/4ZAGMMDRo04IsvvuDIkSMuB3cxJUsREXHHww9DunTkXbkysSk6OprTp0+zYMECFwO7nJKliIi4I0cOqFGDa1asSPzuslKlSuTPnz/oLsUqWYqIiHuio8n0xx+wYQMAERER1K9fnwULFnDy5EmXg/uXkqWIiLinXj1+b9bMOctMEB0dzbFjx1iyZImLgV1MyVJERNxToADbHn8cihRJbKpevTo5cuQIqgIFSpYiIuIqc/YsLFkCe/cCkD59emrXrs3cuXM5m1Bs3W3GXvBAaFpgjKkD1ClYsGDH2NhYzTKfBLf3KRxnmfd2PU/7e9IvLR7D4O5+heMx7O263vSN37qVGh07svWpp9jdqBEAy5cv56WXXmLo0KHcfvvtATuOq1evvs5ae+dlC5KaETotvMqXL69Z5pPh9j6F4yzz3q7naX8dw8lzc7/C8Rj2dl2v+956q7VVqya2xcXF2YwZM9ouXbp4NJ6v/ubAWptEztBlWBERcV90NKxaBX/9BUCWLFmoVasWs2fPJj6hwo+blCxFRMR90dHOs5YX1IWNjo5mz549rF271sXAHEqWIiLivltvhWLF4IsvEptq165NVFRUUNwVq2QpIiLuMwa+/BI++iixKXfu3FSrVo1Zs2ZhXb4ZVclSRESCQ9GiEBl5UVODBg3YsmULO3fudCemBEqWIiISPF57DXr0SHxbr149AFZeUGzdDUqWIiISPHbtgvffh4S6sNdddx0VK1Zk+fLlroalZCkiIsEjOhri4i660adBgwb8+uuvrl6KVbIUEZHgUaMGZMsGF9wBGx0dDeDqXbFKliIiEjwyZIBHHnGetzx3DoASJUpwww03KFmKiIgkatECHnoIjhxJbKpatSorVqxg3759roSkZCkiIsGldm2YMAFy5Upsuueee7DWMveCCj+BpGQpIiLBx1rYuNH5CRQvXpxixYoxa9YsV8JRshQRkeAzZQqUKQPffw+AMYYGDRrwxRdf8M8//wQ8HCVLEREJPg88ABERl90Ve/r0aRYsWBDwcJQsRUQk+FxzDVStelGyrFSpEgUKFHDlrlglSxERCU7R0bBhA2zdCkBERAT16tVjwYIFnEyo8BMoSpYiIhKc6td3fl5yKfbYsWMsWbIkoKEoWYqISHAqWhQWLoQnn0xsql69Ojly5Aj4pVjj9hxhgWaMqQPUKViwYMfY2FiyZs16xf5xcXFX7JPS8lDk9j75a/u+GDe1Y3i7nqf9PemXFo9hcHe/wvEY9nZdX/e9sM+AAQP49ttvmTVrFpEJU3r56m9evXr1ddbaOy9bYK1Nk6/y5cvbpUuX2pSk1MeTMUKN2/vkr+37YtzUjuHtep721zGcPDf3KxyPYW/X9Vnf06etHTjQ/vTaa4lNH3/8sQXsl19+martXQmw1iaRM3QZVkREgldUFIwdy3Vz5iQ21axZk4wZMwa0QIGSpYiIBC9jIDqanN99B4cPA5AlSxZq1arF7NmziY+PD0gYSpYiIhLcoqOJOHcOPv30gqZo9uzZw9q1awMSgpKliIgEt7vv5lSePBc9QlK7dm2ioqICdleskqWIiAS3iAj233efU1Q94QmO3LlzU61aNWbNmoUNwFMdSpYiIhL0tnbpAh9/7HyHmaBBgwZs2bKFTZs2+X37SpYiIhL8zifJY8cSm+rVqwcQkLtilSxFRCQ0vPaaU9Xn7FkArrvuOipVqhSQ7y2VLEVEJDTccgv8/TcsX57YFB0dzfr16/nzzz/9umklSxERCQ01a0LGjJcVVgdYuXKlXzetZCkiIqEhSxYnYX7ySeJdsSVKlKBs2bKsW7fOr5uO8uvoIiIivhQdDXPmwNq1UKECAPPmzWNrwpyX/qIzSxERCR116sDw4c6NPgmuv/76xNlH/EVnliIiEjpy54auXQO+WZ1ZiohIaDl6FMaPhx07ArZJJUsREQkt//wDbdvClCkB26SSpYiIhJZChZybewJURB2ULEVEJBRFR8OaNbB7d0A2p2QpIiKhJ6EYAZ98EpDNKVmKiEjouekmKF0afvwxIJvToyMiIhKavv4asmcPyKZ0ZikiIqEpQIkSlCxFRCSUde0KLVv6fTNKliIiErqshY8/JuLECb9uRslSRERCV3Q0nDxJ7jVr/LoZJUsREQld994LlStjEqbs8hdj/byBYGOMqQPUKViwYMfY2FiyZs16xf5xcXFX7JPS8lDk9j75a/u+GDe1Y3i7nqf9PemXFo9hcHe/wvEY9nZdX/cN1HFcvXr1ddbaOy9bYK1Nk6/y5cvbpUuX2pSk1MeTMUKN2/vkr+37YtzUjuHtep721zGcPDf3KxyPYW/X9XXfQB3HwFqbRM7QZVgREZEUKFmKiIikQMlSREQkBUqWIiIiKVCyFBERSYGSpYiISAqULEVERFKgZCkiIpICJUsREZEUKFmKiIikQMlSREQkBUqWIiIiKVCyFBERSYGSpYiISAqULEVERFKQ5iZ/Ps8Ysx84DPyTQtccKfTJC/zto7CCRUr7HKrb98W4qR3D2/U87e9Jv7R4DIO7x3E4HsPeruvrvoE6jq+31l5zWWtSk1ymlRcw+mr7kMxEoaH88uTvEorb98W4qR3D2/U87a9j2L//vYNt224ew96u6+u+bh/Haf0y7Dwf9Qk3bu+zv7bvi3FTO4a363naX8dw8tzc73A8hr1d19d9XT2O0+xlWF8xxqy11t7pdhwiqaVjWMKBv4/jtH5m6Quj3Q5A5CrpGJZw4NfjWGeWIiIiKdCZpYiISAqULEVERFKgZCkiIpICJUsREZEUKFn6kTHmBmPMWGPMTLdjEfGUMSaLMWaCMeZ9Y0wLt+MRSQ1ff/4qWSbDGDPOGLPPGLPhkvZaxpjNxpitxpjeVxrDWrvNWtvev5GKpMzL47kBMNNa2xGoG/BgRZLhzXHs689fJcvkjQdqXdhgjIkEYoCHgJuBZsaYm40xZY0x8y955Qt8yCLJGo+HxzNQCNiV0O1cAGMUScl4PD+OfSrK1wOGC2vtcmNM0Uua7wK2Wmu3ARhjpgL1rLUDgdoBDlHEY94cz8BunIT5PfoHtQQRL4/jn325bf2P4J2C/PsvbnA+VAom19kYk8cYEwvcbozp4+/gRLyU3PE8C3jUGDOStFtXVkJHksexrz9/dWbpHZNEW7IlkKy1B4An/ReOyFVJ8ni21h4D2gY6GJFUSu449unnr84svbMbKHzB+0LAHy7FInK1dDxLOAjIcaxk6Z01wI3GmGLGmPRAU2CuyzGJpJaOZwkHATmOlSyTYYyZAqwGShljdhtj2ltrzwJdgEXAJmC6tXajm3GKeELHs4QDN49jzToiIiKSAp1ZioiIpEDJUkREJAVKliIiIilQshQREUmBkqWIiEgKlCxFRERSoGQpEuKMMeeMMd9f8Lps6jhjzLKEKYzqGmNiEvr9bIw5ccF6DS9ZJ4sx5oAxJscl7Z8YYxobY5okTIk039/7KOI21YYVCX0nrLXlPOjXwlq7loTqJgmzN8xPbl1r7TFjzGKgPjAhYZ0cQFWgubX2uDHmL+DZq90BkWCnM0sROX8WOc4Ys8YY850xpl7Coik45cPOiwY+s9YeD3yUIu5RshQJfZkuuQzbJBVjvAB8aa2tAFQH3jLGZAE+A8obY/Ik9GuKk0BF0hRdhhUJfZ5ehr2SB4G6xpjzl1QzAkWstZuMMXOBhsaYj4FywOKr3JZIyFGyFBFw5gR81Fq7OYllU4C+CX3mWGvPBDQykSCgy7AiAs6MDU8bYwyAMeb2C5YtBW4EOqNLsJJGKVmKhL5Lv7MclIoxXgPSAT8aYzYkvAfAWhsPfAzkAZb7JGKREKPLsCIhzlobmcr1dgBlEn4/ATxxhb7dgG6p2Y5IONCZpUjacBAYb4yp66sBE+66HQEc8tWYIsFKkz+LiIikQGeWIiIiKVCyFBERSYGSpYiISAqULEVERFKgZCkiIpKC/wPgd6AfS2GMuwAAAABJRU5ErkJggg==\n", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - }, - "output_type": "display_data" - } - ], - "source": [ - "plt.figure(figsize=(7, 7))\n", - "plt.loglog(E, sed, 'k-')\n", - "plt.loglog(E, A_blr*sed, 'r--')\n", - "plt.xlabel(\"E [TeV]\")\n", - "plt.ylabel(\"1/ cm$^2$ s TeV)\")\n", - "plt.grid(which='both')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "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.8.8" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -}