{ "cells": [ { "cell_type": "markdown", "id": "1f4c5d87", "metadata": {}, "source": [ "# The Metagalactic UV Background" ] }, { "cell_type": "markdown", "id": "f080572d", "metadata": {}, "source": [ "One of the main motivations for ARES was to be able to easily generate models for the metagalactic background. In this example, we'll focus on the ultraviolet background , which is noteworthy given the ''sawtooth'' modulation (e.g., [Haiman et al. (1997)](http://adsabs.harvard.edu/abs/1997ApJ...476..458H)) caused by intergalactic hydrogen atoms.\n", "\n", "In order to model this background, we need to decide on a few main ingredients:\n", "\n", "* The spectrum of sources, which can be one of several pre-defined options (like a power-law, ``pl`` or blackbody, ``bb``), or a Python function supplied by the user.\n", "* How the background will evolve with redshift, which could be based on the rate of collapse onto dark matter haloes as a function of time, a parameterized form for the star-formation rate history, or more detailed models of star-formation (see [this example](../uth_pop_sfrd.html)).\n", "* What (if any) approximations we'll make in order to speed-up the calculation, aside from the assumption of a spatially uniform radiation background, which we make implicitly in ARES throughout.\n", "\n", "First things first:" ] }, { "cell_type": "code", "execution_count": 1, "id": "2af090a5", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline\n", "import ares\n", "import numpy as np\n", "import matplotlib.pyplot as pl\n", "from ares.physics.Constants import erg_per_ev, c, ev_per_hz" ] }, { "cell_type": "markdown", "id": "43e782d5", "metadata": {}, "source": [ "Now, let's set some parameters that define the properties of the source population:" ] }, { "cell_type": "code", "execution_count": 2, "id": "14938fb3", "metadata": {}, "outputs": [], "source": [ "alpha = 0. # flat SED\n", "beta = 0. # flat SFRD\n", " \n", "pars = \\\n", "{\n", " 'pop_sfr_model': 'sfrd-func',\n", " 'pop_sfrd': lambda z: 0.1 * (1. + z)**beta,\n", "\n", "\n", " 'pop_sed': 'pl',\n", " 'pop_alpha': alpha, \n", " 'pop_Emin': 1,\n", " 'pop_Emax': 1e2,\n", " 'pop_EminNorm': 13.6,\n", " 'pop_EmaxNorm': 1e2,\n", " 'pop_rad_yield': 1e57,\n", " 'pop_rad_yield_units': 'photons/msun',\n", "\n", "\n", " # Solution method\n", " 'pop_solve_rte': True,\n", " 'lya_nmax': 8,\n", " 'tau_redshift_bins': 400,\n", "\n", "\n", " 'initial_redshift': 40,\n", " 'final_redshift': 10,\n", "}" ] }, { "cell_type": "markdown", "id": "48413474", "metadata": {}, "source": [ "To summarize these inputs, we've got:\n", "\n", "* A constant SFRD of $0.1 \\ M_{\\odot} \\ \\mathrm{yr}^{-1} \\ \\mathrm{cMpc}^{-3}$, given by the ``pop_sfrd`` parameter.\n", "* A flat spectrum (power-law with index $\\alpha=0$), given by ``pop_sed`` and ``pop_alpha``.\n", "* A yield of $10^{57} \\ \\mathrm{photons} \\ M_{\\odot}^{-1}$ of star-formation in the $13.6 \\leq h\\nu / \\mathrm{eV} \\leq 100$ band, set by ``pop_EminNorm``, ``pop_EmaxNorm``, ``pop_yield``, and ``pop_yield_units``.\n", "\n", "See [this page](../params_populations.html) for a complete listing of parameters relevant to :class:`ares.populations.GalaxyPopulation` objects.\n", "\n", "Next, let's initialize an :class:`ares.simulations.MetaGalacticBackground` object (which will automatically create an :class:`ares.populations.GalaxyPopulation` instance):" ] }, { "cell_type": "code", "execution_count": 3, "id": "966cdd51", "metadata": {}, "outputs": [], "source": [ "mgb = ares.simulations.MetaGalacticBackground(**pars)" ] }, { "cell_type": "markdown", "id": "28c8c7a3", "metadata": {}, "source": [ "To run the thing:" ] }, { "cell_type": "code", "execution_count": 4, "id": "cc9af89d", "metadata": {}, "outputs": [], "source": [ "mgb.run()" ] }, { "cell_type": "markdown", "id": "32f81828", "metadata": {}, "source": [ "The results of the calculation, as in any ``ares.simulations`` class, are stored in an attribute called ``history``. Here, we'll use a convenience routine to extract the redshifts, photon energies, and corresponding fluxes (a 2-D array):" ] }, { "cell_type": "code", "execution_count": 5, "id": "f046a937", "metadata": {}, "outputs": [], "source": [ "z, E, flux = mgb.get_history(flatten=True)" ] }, { "cell_type": "markdown", "id": "29d35c32", "metadata": {}, "source": [ "Internally, fluxes are computed in units of $\\mathrm{s}^{-1} \\ \\mathrm{cm}^{-2} \\ \\mathrm{Hz}^{-1} \\ \\mathrm{sr}^{-1}$, but often it can be useful to look at the background flux in terms of its energy, i.e., in units of $\\mathrm{erg} \\ \\mathrm{s}^{-1} \\ \\mathrm{cm}^{-2} \\ \\mathrm{Hz}^{-1} \\ \\mathrm{sr}^{-1}$:\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "0d7d5834", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD4CAYAAADSIzzWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABKo0lEQVR4nO3dd1hUx9cH8O8A0osCGhVURBRRESkWLEhRrGDBgjHWaDRq7D9rjBhbYo0tJnaNGGPHThQRsUGsCCogiCIgIFUpC8vO+we6ryuLAnt37y7M53n2MczenTk34PFy7twZQikFwzAMo3hqfAfAMAxTU7EEzDAMwxOWgBmGYXjCEjDDMAxPWAJmGIbhiQbfAVSGqakptbCw4DsMhmGYSrl79+4bSmndT9tVKgFbWFjgzp07fIfBMAxTKYSQF9LaWQmCYRiGJywBMwzD8IQlYIZhGJ6wBMwwDMMTloAZhmF4whIwwzAMT1gCZhiG4UmNSMBbtmzB+fPn+Q6DYRhGQo1IwKtWrUJAQID4a7YGMsMwykAlnoQjhHgB8LKysqrS51+9eoWSkhIAwNq1a/HgwQMcPHgQhBAOo2QYhqkclbgCppSeoZR+Z2RkVKXPq6urQ1NTEwAgFApBKUVRURGXITIMw1SaSlwBc2nBggUAwK5+GYbhnUpcAXOJEAJCCF6+fIlRo0YhOztbYWNTSlFYWKiw8RiGUW41LgF/kJaWhrNnz+LevXsKG3Pr1q3Q0dFBQUGBwsZkGEZ51bgSxAdOTk54+fIlDAwMFDamo6MjvLy8oKOjAwCYPHky7O3tMWnSJIXFwDCM8qixV8AAxMn32rVrKC4ulvt4nTt3xunTpwEAxcXFiI2NRXJyMoDS8sSTJ0/kHgPDMMqjRidgAAgLC0P37t2xd+9euY+Vm5srTri1atVCUFAQli5dCgAIDQ1Fq1atJOYrMwxTvdX4BNyhQwccPHgQo0ePlvtYq1evRtOmTSXa1NRKvwW2trZYt24devbsCQC4cuUKDh8+DJFIJPe4GIbhB1Glp8KcnJyoKm9JFB4ejoiICEyYMOGLx44YMQJ37tzBkydPoKFRY0v1DFMtEELuUkqdPm2v8VfAH9y+fRs9evRATk6O3Mbo0KFDhZIvAPj7++Py5cvQ0NCAUCiEt7c3Ll68KLfYGIZRPJaA39PS0kJCQgKeP38utzGysrKQkJBQoWPV1NTQpEkTAEBSUhLi4uIgEAgAsLUsGKa6YAn4PXt7e8TExKBdu3ZyG2PNmjWwtrau9OeaNGmCiIgIeHt7AwC2b9+OIUOGIC8vj+sQGYZRIFZc/IiamhpKSkrw+vVrmJmZcd7/kCFDYGNjU6XPqquri/9bKBSiqKgIurq6AEqviNmj1QyjethNuE/07t0b2dnZuH37tlzHkdWHpJuVlYXevXtj7dq1cHFx4TsshmGkKO8mHLsC/sTkyZNRUlIil6vK9PR0ZGdno3nz5jL39SG2169fo7CwEIaGhjL3yVRMUVERSkpKxE80MkxVsStgBVq4cCE2bNggvpnGlY//sVi9ejXMzMwUMq+5pnJzc0NGRgYiIiIAAImJiTA3N2dlIKZcbBpaJeTk5GDXrl2crxk8fPhw7N+/n9M+gf+/Gi4pKcGlS5cQEhLC+RjM/5s8eTJ8fX0BAAKBADY2NuJlTgHIdSojU72wK2Apzp8/j379+uHChQvo3bu33MfjUklJCYqKiqCjo4Pnz5/j4cOHGDhwIN9hVVsFBQXw9/eHnZ0d2rdvj+fPn8PKygr+/v7w9fWFSCQSL4HK1FzsCrgSPD098d9//6FXr16c9vv69WtERkZy2uen1NXVxbXJdevWYcyYMcjIyJDrmDVNTEwMXr58CQDQ0dHBhAkT0L59ewCl88kXL16MDh06AAAuXboECwsLREVF8RYvo7xYApZCQ0MDTk5OnF+1bNq0CY6Ojpz2+TkbNmzAlStXYGJiAqD0QRBGdgMHDsTcuXOlvtewYUP8/PPPsLS0BFC64l779u3FX+/Zswe+vr5sYX4GAEvA5crNzcWCBQtw+fJlzvocOXIkDh8+zFl/X6KlpSVO+MePH0ezZs3w4MEDhY1fXf3222+YPXt2hY7t3Lkzjh07Jv6tJDc3F6mpqdDW1gYA7Ny5E3/99ZfcYmWUG6sBl0MoFKJRo0aYOXMm5s+fr5Ax5enZs2dYt24dNm3aBC0tLb7DYd5zcXGBsbExTp06BaC0ZOHk5IQ6derwGxjDqfJqwCwBf4ZAIOA0WSUlJSE1NRUODg6c9VkV+fn52LRpE+bOnYtatWrxGosqioqKgo6OjrisIAtKKXJzc2FkZITs7GzUq1cPM2bMwNq1a0Epxdu3b9kc72qA3YSrAq6vFLdu3QpnZ2dO+6yKs2fP4scff8StW7f4DkUlDR06VGLamSwIITAyMgIAGBkZ4fr165gyZQoAICIiAl999RVbBa8aU9iTcISQgQD6AagHYBul9F9CiA2AGQBMAQRRSrcrKp6KEIlEGDZsGBwdHbFw4UKZ+xs9ejS6du3KQWSyGTZsGNq2bYuWLVsCKD3PDwvDM1+2bds2uewlSAgRz54AAENDQ0yYMAFOTqUXToGBgQgMDISfnx+7Kq4mKvS3jhCyhxCSRgiJ/KS9NyEkmhDyjBDy2UsCSukpSulEAGMBDH/f9oRSOhnAMABlLs/5pqamBk1NTYmFcGRhY2ODfv36cdKXrD4k3/DwcDg4OCA+Pp7niFSHm5ubOCnKU9OmTbFlyxaYmpoCAB4+fIjjx4+LF2F68OABm2Ko4ip62bMPgMQTCYQQdQDbAPQB0ArACEJIK0KILSHk7Ceveh999Mf3n/vQjzeA6wCCZDgPuTl06BDmzZvHSV8vXrxQyl/7tbS02LoGlXD//n1ER0crfNx58+YhNjYWGhoaoJRi9OjR8PLyUngcDIcopRV6AbAAEPnR184AAj/6eiGAhZ/5PAHwK4Ae5bx/rpz27wDcAXCncePGlC9CoVDmPhYsWEA1NTU5iIZbIpFI/GdKSgrP0Sg/GxsbOnToUL7DoBEREfT69euUUkoFAgHt1q0bPXPmDM9RMdIAuEOl5DdZCn9mABI/+vrV+7by/ACgB4AhhJDJAEAIcSWEbCaE/AngvLQPUUp3UEqdKKVOdevWlSHcqqGUwsnJCdOnT5e5r/Hjx+Ps2bMcRMWtDw+crFixAu3atRPv3MxIt3PnTixevJjvMGBra4suXboAAJKTk1FcXCye1ZKVlQVVXriqppDlJpy0x8TKndNGKd0MYPMnbVcBXJUhBrkjhKB///6cTDlq3rw5J0tRysuQIUNQVFSEBg0a8B2KUvuQ9JSJhYUFbt26Jd6uaufOnZg/fz5iY2NhZWXFc3RMeWRJwK8ANProa3MA1fLSyc/Pj5N+4uPjkZiYiO7du3PSH9dsbGywfPlyAKVrF4tEInz11Vc8R6V8wsPDoaenh9atW/MdShkffpuZNGkSLCwsxMn3l19+gb6+PqZNm8ZneMwnZClB/AegOSGkKSFEE4AvgNPchCWJEOJFCNnB5zJ/hYWFMi8zuHPnTnh6enIUkfyIRCL06tULPj4+bANQKcaOHYtly5bxHcZnGRkZYdiwYQBKy2ihoaEICwsTv//27Vu+QmM+Jq0w/OkLwN8AUgAUo/TK99v37X0BxACIA7C4In3J8nJ0dOS+Ol4BBQUFVFdXly5ZskSmfuLi4ujVq1c5ikq+zp8/T0NDQ/kOQymFhYXRyMhIvsOotMLCQkoppc+fP6d6enr0yJEjPEdUc6Ccm3AVKkFQSkeU034e5dw8q060tbWxYsUKiUnyVWFpaclJLVkR+vTpI/7vmJgYtGjRgsdolIusPwd8+fBkp6amJkaNGiV+KjMuLg5FRUVV3jCWqTq2FoQCxcbGIiEhAT179uQ7lAoLDQ2Fm5sb/P39MXz4cL7DUQo3btyAvr4+7Ozs+A6FE6NGjcKZM2eQnJwsfsiD4RZbC0JGJSUlePr0KXJzc6vcx549e9C/f38Oo5I/Z2dn+Pn5oW/fvnyHojQmTpyIlStX8h0GZzZs2IB//vlHnHw3bdqEFy9e8BxVzaASV8CEEC8AXlZWVhNjY2N5iSE8PBwdO3bEiRMnMGjQoCr18eLFCyQnJyvFgjxVUVJSgry8vBq/DsH9+/ehq6sLa2trvkPh3MuXL9GiRQssW7asWizDqizYcpQyKigowOHDh+Hp6Qkzs889b1I9UUrRr18/UEpx/vx5tsdZNZaYmAgTExPo6uoiJCQEoaGhmDt3rngReabyWAlCRjo6Ohg3bpxMyffJkyc4d+4ch1EpDiEEPj4+rA4MIDg4uFo/ZdaoUSNxOeLixYvYvXs3m44oJ+wKuBLS09Nx7969Km/WuXDhQmzYsAECgYDjyBhFatWqFdq0aYMjR47wHYpCZGVloU6dOigpKcH333+P77//Hvb29nyHpVLYFTAH9u3bh969e1d5c8tp06Yp5WpolXX+/Hm4u7vX2I0ljx49il9++YXvMBTmw/ZIcXFxCAgI4GUluOpKYQuyy+Kjm3C8xjFs2DA4OztDT0+vSp83MzOrFvVjNTU15OTkID09HY0aNfryB6oZZXwEWRFatGiBZ8+eQV9fHwBw7NgxFBYWYuTIkeyeQBWxEoQCRUZGIjY2tsqzKJRJTd5F499//4WhoSE6derEdyi88vb2Rnp6Om7cuFFjfxYqipUgOBISEoIbN25U6bP+/v7w9fXlOCJ+qKmpoaCgAOvXr4dQKOQ7HIWaOXMmNmzYwHcYvDt16hROnz4NNTU15OXlYfv27SgpKeE7LJXCroArqV27djAzM6vSbIbXr1/jzZs3aNOmjRwiU7wzZ87A29sbFy5cQO/evb/8gWoiJiYG2traaNy4Md+hKI0dO3Zg8uTJCAsLQ/v27fkOR+mwecAciYqKgomJCerXr89rHMqAUoqHDx+iXbt2fIfC8IxSirt374r3youIiICtrS2rDb/HShAcad26dZWT74MHD3D48GGOI+IPIUScfLOzs3mNRZHOnj2L69ev8x2GUiGEiJNvfHw8OnToUK0e15YXlUjAyrAe8AcpKSnYuXMnUlNTK/3Zf/75B2PGjJFDVPwKCgqCubk57t69y3coCjFv3jxs3rz5ywfWUE2bNsXmzZsxceJEAKU3bBnpWAmikm7cuIGuXbvi/PnzEks2VkR6ejqys7OVeluiqsjNzcWcOXOwePFiWFhY8B2O3CUkJEBTUxMNGzbkOxSlRynF0KFDYWdnhyVLlvAdDm/KK0GoxDxgZeLo6IgXL15UaT5v3bp1wcfGovJmaGiInTt38h2GwtSEf2S4UlxcDENDQ/HcYUaSSpQglMmHu9/q6uqV/uzdu3dx4MABOUSlHBITE7F8+fJqv27AyZMncfXqVb7DUAmamprYs2cPZs6cCQAICwvDvXv3+A1KibAEXAUnT57EoUOHKv25Y8eOieti1dG///6LFStW4NGjR3yHIleLFy/G77//zncYKoUQAkopZsyYgTFjxrC68HusBlwF/fv3R3JycqX/Jc/KykJOTk61/RW2uLgYaWlp1eJx689JSkpCrVq1UK9ePb5DUTlpaWnIysqCtbW1eF+0mvAUHZsHzKGsrCwYGBhAQ+P/S+ju7u7o06cP/ve///EYmfIQCoUS/38Y5lM//fQTYmNjceDAAdSqVYvvcOSKzQPmUJ06dSSSC6UUwcHB+NJuHeHh4di1a5e8w+PdDz/8AC8vL5n72bp1q1LOJT1y5AguX77MdxgqT19fH/r6+jX6H2qVSMDKNA8YAJ49e4alS5ciKSlJov1Lv3qfPHkSU6dOlWdoSqFly5awt7eXeV2A1atXw9/fn6OouOPn54cdO3bwHYbKmzdvHnbs2AFCCFJSUsr8faoJWAmiCq5duwZXV1dcuXIFrq6uoJQiJSUFBgYGMDAwKPdzubm5ePfuHZs/WkFv3rwBAJiamvIciaS0tDRoaGjA2NiY71CqBUopXF1dkZ6ejoiIiGp5RcxKEBzq3LkzBAIBXF1dAZTe4d26dSsGDx782c8ZGhrWqOT733//VXnxeqA08QoEAgwZMgT379/nMDLZ1KtXjyVfDhFCsHHjRmzZsqVaJt/PYQm4CjQ0NCRuGlBKceLEiS+uh3Dz5k1s27ZNztEph5iYGHTo0AF79uypch8nT57ExYsXce/ePTx+/JjD6GRz8OBBXLhwge8wqhUHBwd4eHgAKF1rIyoqiueIFOTDVBBVeDk6OlJl8euvv9KDBw9SSikViUQUAPXz8/vsZxYsWEA1NTUVEZ5SOHLkCM3Jyany5z08PGiXLl2oUCjkMCrZ2djY0KFDh/IdRrVUWFhImzRpQvv06cN3KJwCcIdKyWmsBlxF7dq1g729Pfbu3QtKKTIzM6GjoyPeTVaavLw8FBYWwsTERIGRqq7c3FxQSmFkZASgdA0GZZhDnZ2dDXV19c/W+5mqi4+Ph6mpKQwNDfkOhTOsBsyxe/fuYe/evQAgvotrbW2NixcvlvsZPT29Gpd8z549i7Vr11bps4aGhuLku2nTJrRo0QKvXr3iMrwqqV27Nku+cmRpaQlDQ0MIhULMnz8fKSkpfIckNywBV9HHT+9QSvHXX39BR0fnswk2NDQUv/32mwKiUx6BgYHYsWNHlbYtOnr0KP7++28AwIABA/Drr7+idu3aHEdYeXv37sWZM2f4DqPai4mJwe+//169/19Lq0so2wuAF4AdVlZWHFdmqi4wMJBOmDCBlpSUsBrwZ+Tm5la5hvuhBqxsWA1YcV69esV3CJwAqwFza+vWrVi9ejWePHkCAwMDvHv3DpqamtDU1Cx3GxaBQICioiL262sF5efnA4C4ri4SiXDu3Dloa2ujZ8+evMWVl5cHNTU16Ojo8BZDTfPkyRP4+/tj+fLlKrnNEasBc2zatGlISkqCoaEhCCEwMDDAiBEjPpsYtLS0amTyvXXrFlq3bo2XL19W6nO6uroSNzUJIZg/fz42bdrEdYiVoqenx5Kvgp06dQo7d+5EcnIy36FwiiVgDlBKsWrVKmhra6N///7lHhccHIxff/1VgZEpB1NTU5iamiIjI6NSnzt06JDE+smEEJw9exYnT57kOsRK2bFjB06cOMFrDDXNggULEBERUf1W2pNWl1DWlzLNA05KSqJjx46lt27dYjVgOWE1YOZTIpGIbt++ncbExPAdSqWgnBow70m1Mi9lSsCvXr2i5ubm9NixY1QkElGBQECFQqH4ppw0xcXFVCAQKDhS5VFYWEjz8/MrfHxRUREtKioq037x4kXatWvXSvXFJYFAQIuLi3kZu6ZLS0ujJiYmdObMmXyHUinlJWBWgqgiMzMzJCYmwsfHB4QQaGpq4uDBg9DU1Cy3TqWhoQFNTU0FR6ocEhMTUbduXRw8eLDCn6lVq5bUdWJr1aoFgUDA2+pZmpqaNW7NAmVRt25d3Lp1C+vXr+c7FE6wBMwBSin8/PyQlZWFBQsWlJtkL1++jOXLlys4OuVgbm6OyZMnw87OrsKfOXDggNS1JNzd3REeHg4rKysuQ6ywbdu24ejRo7yMzQDNmzeHmpoasrKy8O+///IdjmykXRYr60uZShCUUjpnzhy6Zs0aVgOWky/VgAUCAS9lCFYDVg5jxoyhtWvXlmm9EUUBmwfMvUGDBsHMzAxbt24VtwmFQpSUlEBLS4vHyJRXQkIC1NTU0LhxY5n6SUlJga2tLVauXIlJkyZxFB2jSl69eoXXr1/DyanM9Fqlw+YBy8HJkyclkm9eXh40NTV5n6eqrAoLC2FjY4MNGzbI3Ff9+vUxatQotGnThoPIGFVkbm4uTr7p6ek8R1M1LAFzgFKKRYsW4fbt21i6dCk6d+4s9bjAwEAsWbJEwdEpD21tbfj7+1d4W6bdu3fjzz//lPreh0W8u3TpwmWIFbJx40bxGhUM/44ePYomTZogMjKS71AqT1pdQtleUMK1ICildMeOHXT48OFUJBLRWrVq0Z9//vmzxy9evJjq6ekpKDrV16tXL9q9e/fPHpOcnEyjo6MVE9B7bdu2pb6+vgodkylfeno6/eGHH2hmZibfoZQLrAbMvfXr1+P06dMICQkRtxUVFaGgoEC8jCIjqaSkBJcuXUK9evXg4OAgU1+UUjRu3BgdOnTA8ePHOYqQYbjHasByMGfOHInkCwD9+/dH7969eYpI+RFCMGLECGzfvp2Tvv7880/8/PPPHETGqLrnz5/Dx8dHpdYPZrPJOUApxdy5c9GrVy9MnjwZAoFA6nHnzp1DSEgI1qxZo+AIlYeamhqCg4PRvHnzLx77xx9/QCgUYtq0aeUe07dvXy7Dq5A1a9agfv36GD16tMLHZspXUlKC69evIyIiAg0aNOA7nAphV8AyCA0NhYeHBxISErBjxw7cvXsXgwcPxogRI6QeHx4ejt27dys4SuXTrl076OnpffG4c+fOISAg4IvHhYaGIjAwkIvQKuTIkSO4fPmywsZjKsbKygovXrxAr169+A6lwlgNWAahoaFYuHAhdu/eDWtrawCla/6+efMGDRs2VMl1SxWhoKAAO3bsgKOjI7p27Spzf927d0deXh6U6WeD4VdISAi6desmsXMNn1gNWA66deuG69evi5MvULp3mbm5uXgxcaasWrVqYdGiRZxt7b5r1y5cuXKFk74Y1RcUFARXV1eVuDHLasAcoJRi+vTp6Nu3Lzw9PVG7dm2p//KePn0awcHB2LhxIw9RKg8NDQ0kJibC2Nj4s8dt3boVxcXFmDVr1mePq0g9mUsrV65EgwYNMH78eIWOy1SMm5sb9uzZAy8vL75D+SJ2BSyDjIwMuLi44OTJkzh06BAePXqEdu3a4bvvvpO6Y8LDhw9x+PBhHiJVPl9KvkDplcylS5cq1N+BAwewc+dOWcOqkLNnz+L69esKGYupPDU1NYwbNw7a2tp8h/JFLAHLQFNTE+rq6lBTU0NGRgbmzZuH4uJivHjxAu/evStz/JIlS1Rqiow8PX/+HFOmTMGTJ0/KPebkyZM4f/58hfo7evRopZa6lMWtW7ekrtLGKJfr16+jd+/eSl0OZAlYBgYGBggODsbAgQPFbZGRkbCwsGB3yb+AUoqDBw8iNjaWk/4OHz6Mq1evctIXUz2IRCLEx8fj+fPnfIdSLpaAOUApxXfffYczZ86gadOm2L17N+zt7cscd+LECUyZMoWHCJVP06ZNkZ2dDW9v73KP2bBhQ4XnTOvp6Sls1snSpUvLXaOCUR4uLi548uQJWrduzXco5WI34WTUt29fdO7cGWfPnoW1tTW8vLzKvTnz9OnTCv9KXd0RQr6YMMPCwlBYWFjhPufPnw8zMzNMnz5d1vA+6+rVq2jZsqVcx2C4oa6uDqFQiFevXsHCwoLvcMpgCVhGpqamMDQ0lNiG6NmzZ9DX10f9+vUljl20aBEWLVqk6BCV1rlz57Bt2zYEBARI3Xron3/+qVR/ERERCqn3ffr4OaPchgwZgpiYGERGRirNvOAPWAKW0cfbpn/Qtm1bTJ06FWvXruUhItWRl5eH1NRUpKeno2HDhjL3d/78efbwC1PGlClTlPZGHHsSjgOUUowbNw4DBw7EwIEDceTIEbRo0QLt2rWTOO7o0aMIDAzErl27+AlUxaxZswbFxcVYvHgx36FIWLRoEczMzCq8rjHDqPSTcIQQL0LIjpycHL5DKWPKlCn4+uuvcfXqVSQkJAAAhg0bVib5AkB8fDyuXbum2ABV2MOHD3H//v0KH5+fnw9PT0+5TxELDw/H48eP5ToGw613795h27ZtnM264Qq7ApbRihUr8PbtW/z666/itri4OIhEIoU/oaWKZsyYAaFQiG3btsncF6UUPXv2xNdff82eUmMkvH79Go0aNcKqVavwv//9T+Hjl3cFzGrAMvrxxx/LtI0aNQp6enoVfoqrJqtVqxZndVtCCJt/zUhVv359REdHw9LSku9QJLArYA5QSjFy5Ej4+PjAx8cHISEh0NDQKLNf2eHDh3H+/HmpN+6YslatWoWioiL4+fnxHYqE//3vfzA3N8eMGTP4DoVRESpdA1Zm69atg62tLe7cuSOeita9e3epm0UmJSXh7t27ig5RZUVHRyM6OrpSn/nvv//QtGlT3L59W05RlT7tGB8fL7f+GfnZsmULBgwYwHcYYiwBy8jMzAyOjo6IiYnBDz/8AAB4+fKl1LVp58yZg6ioKEWHqNQSExPRpk0bqUsH7t+/v9K7D3/11Vfo2LEjtLS0uAqxjAsXLmDTpk1y65+RH6FQCDU1NRQVFfEdCgBWgpCLadOm4fDhw3jy5AkopahXrx7fISmtgoIC+Pr6Ytq0aejZs6fcx8vMzERhYSEn844ZpqJYCUKOKKXw8fHBkSNHAADff/89/vnnH3h7e6Nfv37i4w4ePIjhw4fzFaZS0tHRQUBAgNTk6+fnV+U5wCUlJVLbx44dK/M+cjNmzMD69etl6oPhV1ZWFpTh4pMlYBkdPHgQjRo1ws2bN5GRkQEAaN26NTw8PNC6dWu0atVKfOybN2+Ubh6ispD2lyEpKQmJiYmV7mvx4sXlPvdvaWkp853whIQEiUfPGdVy8eJFfPXVV3jw4AHfoZT+4KvKy9HRkSqbK1eu0HHjxtHMzExxW0pKCr1y5QoVCAQ8RqY6li9fThs3bsxZfwEBAXTx4sW0uLhY6vslJSW0qKiIs/EY1ZKRkUFnz55Nnz9/rrAxAdyhUnIaqwHLwa5duzBx4kQMGDAAcXFxePToEd8hKbWTJ08iJCQE69atg4aGfKemR0dHw97eHvv27cOwYcPkOhbDfMBqwHJEKYWXl5f4jn3fvn1x5coVhIWFITIyUnzc/v37MWjQIL7CVFqDBg3Cb7/9Vib5/vjjj5g3b16V+iwpKZG6AIu/vz8KCgrQrFmzKvULlD5+/vGTj4zqoZTiwYMHePnyJa9xsAQso+DgYHz11Vc4e/Ys3r59CwBo2LAh3NzcMGDAAIlHYnNzc9mWROWglJa5cZaRkSGuq1e2LxMTE6kPcOTn52PYsGFwdHSsaqhIS0tDZmZmlT/P8C87OxtOTk74448/+A1EWl1CWV/KWAN+/PgxnTRpEo2NjRW3ZWRk0MDAQIm6MFO+7OxsamBgQDdu3MhZn2vWrKH//vuv1PdEIhFNTU3lbCxGNZ05c0ZhPwdgNWDFuXLlCjw8PODl5YVLly4hOztbrg8GqDpKKWbPng1vb2+4ubnJfbylS5di1apVKCgokHvNmWEAVgOWK0opPD098ddffwEAHBwccP36dURHR0tsqbN7926JecFMKUIINm7cWCb5zp8/H7NmzapSnyKRCElJSRCJRBLty5YtQ0BAADZs2AChUFilvidMmIAVK1ZU6bOM8hCJRPjnn38QFBTEWwwsAcvo2bNnqFOnDi5duiR+vLF27dro0qULPD09MXnyZPHVb1FRkdTt6plSeXl5El8XFBSgoKCgSn3t3LkT5ubmZebrFhYWwtbWFj/88AO0tbWr1Hd+fn6l9qpjlJOamhoWL16MHTt28BYDK0HI6M2bN/j5558xatQotG/fHkBpIgkODoadnR3Mzc0BgG2V8wUzZszAoUOHkJ6ezkl/MTExCAoKgq+vL+rUqSPxHqUUGRkZoJSibt26nIzHqKYXL17A3Nwc6urqch2nvBIES8By8Pz5c1haWsLb2xtnzpzBw4cPYWtry3dYSu3SpUuIiIjA7Nmz5f6PlUgkgq6uLqZPn17hbe8ZRhasBixHlFK4urpi3759AEqnoYWHhyMlJUU8JQoAduzYoZAFZ1RRz549MWfOHInkO2fOHPEKc1WRkpKCpKQkibYlS5ZgwoQJ+OOPPzB06NAq9Tt69GgsW7asynExymXlypXiv7uKxhKwjPLz86GnpyexVbmWlhbat2+Prl27Ytq0aWzlrQqglCInJ6fKNV9pOnToUO5iPmPHjhWXjJia7dy5c7h+/TovY7MShIyEQiEWLFiAfv36ie/ii0QinDp1Ci1btkSLFi1ACJF7jUnVRUREwM7ODseOHYOPjw8nfR4/fhwNGzaEs7Nzmffevn2LmJgYODg4sPp8DScUCuU+HZHVgBVIJBJBXV0dPj4+OH78OKdJpbrKzc3Fzp070b9/f1hbW8t9vE2bNmHmzJlITU1l6zUzcsdqwHJEKUXnzp2xe/duAKXTWx48eIDs7GwAQMuWLQEA27dvR/fu3fkKU6kZGhpizpw5Esl3+vTpmDx5cpX7zMrKQlhYmMRSlwsXLsSoUaPg5eWFkydPQk9Pr9L9jhgxAkuWLKlyXIxyoZTim2++wZYtWxQ+tsISMCFkICFkJyEkgBDi+VG7HiHkLiGkv6Ji4ZqxsTFu3boFTU1NcZudnR0cHBwwa9YstG7dGgCgqakJfX19vsJUejk5OXj9+rX4ax0dHejo6FS5v7/++gudOnXCmzdvxG3a2trQ1dWFpaUlBg4cWKUErKurW+U5xIzyIYQgKytLvJaLQseuSAmCELIHQH8AaZTSNh+19wawCYA6gF2U0l8q0FcdAOsopd++//pnAHkAoiilZz/3WWUtQfz8889wcnKS2GkhICAAdevWhb29PdTV1SWSMyOdk5MT6tatiwsXLnDSX3x8PB4/fgx3d3fo6upKvCcSiXD//n2YmJiUu3g7w3BFphowIcQFwDsABz4kYEKIOoAYAD0BvALwH4ARKE3Gqz/pYjylNO3959YD8KeU3iOE9ABgCkAbwBtVTcDSWFhYwNbWFmfPnsXmzZtlmk5VU5w6dQra2tro3bu33McqKSmBlpYW5s+fj5UrV8p9PKZmk6kGTCm9BuDT9fc6AHhGKY2nlBYBOAxgAKX0EaW0/yevNFLqVwAXKKX33vfhBqATgK8BTCSEqGRNWiQSwc7OTuKRxqCgIPE6BJ07dwZQuiV2p06deIlRFQwcOFAi+X7//ff49ttvq9yfSCRCWFgYnj17Jm6bO3cuhg8fDnV1dZw5c6ZK/fv4+GDBggVVjotRPhkZGbC3t8eePXsUOq4scy/MAHy8YdcrAB0/c/wPAHoAMCKEWFFK/6CULgYAQshYlF4Biz79ECHkOwDfAUDjxo1lCFd+mjdvjvj4eBgYGIjbmjVrhnbt2qFNmzbitWcNDQ3RoEEDvsJUerm5uUhKSoK1tTXU1NRgYmIi0/bhhBC4uLhg5syZ4gXUjY2NIRAIAAB9+vSpUr/16tUr83gzo9qMjY1haWmJ2rVrK3TcCk9DI4RYADj7UQliKIBelNIJ778eBaADpVRuv2srawli69atMDMzk9jt4ty5c9DQ0ECnTp1Qq1atMjVIpqyNGzdi9uzZyMzM5CzBXb58udyNOJ89e4YXL17Aw8ODk7EYpjzllSBkuQJ+BaDRR1+bA6iRW8VOmzatTNuKFStgYGCAq1evYu7cuVi1ahUPkamWvn37okGDBpzOMOjRo0e5723evBn79u1Dbm4uZ+Mxqu1D2VBNTTHVUFlG+Q9Ac0JIU0KIJgBfAKe5CUu1CAQCWFlZYfv27eK248ePQ1tbG8XFxeI1gH/77Tc4ODjwFabSs7a2hq+vr3jq2cSJEzF69GiZ+oyKikJwcLD465kzZ2Lw4MEASucZX7lyBRX9LfADb29vzJ07V6a4GOVz+fJlGBkZKXQT3QpdARNC/gbgCsCUEPIKwFJK6W5CyDQAgSid+bCHUholjyAJIV4AvKysrOTRvcxcXV0RFxcnXnQHKF2Qx8HBAW3btkWXLl0AAKampmjevDlfYSo9gUCAmJgYmJmZwdjYGGZmZiguLpapz19//RXXrl1DQkICgNLvy4erm6r+PFlYWLBafjVkZWWFsWPHKrRcyB5F5sDBgwehoaEBX19fcdvFixeRl5cHV1dXqKurK7y4r4piYmJgbW2NgwcPYuTIkZz0+WFXEjs7uzLv5ebm4sqVK3B0dESjRo2kfJphuMHWglAwLy8vJCUl4c2bN/Dw8MDevXv5Dknp5efn4/z58+jYsaNCEmJ0dDRatmyJv/76C998843cx2NUQ05ODoyMjDjtk60FIUfv3r1D/fr1JZ4l37dvHxo2bIjExESMGjUKALB+/XrxY8lMWbq6uhgyZIg4+Y4ZMwYjRoyQqc/Xr1/j+PHj4nU5pk6dCi8vLwClpYS7d++Kv66o3r17Y+bMmTLFxSin6dOnK7RMqBJbwip7DXj48OFITU2VWPfXxMQEzs7OaN++Pdzd3QEAZmZm4jnBjHSPHj2ClpYWWrRoAWtra5nmAQPAvXv3MGTIENy6dQudOnWClZWVeD0OLS2tKt0UbdOmjdLOSWdk4+3tDUtLS5SUlChkCVlWguBAQEAA8vLy8PXXX4vbLl26hFevXqFfv36glOKrr77iMULV0axZM3Tq1An+/v6c9JeTk4OEhAS0aNFC6sI+Fy5cgI6ODlxdXTkZj2GkYTVgBRs7diyCg4NRv359GBsbc7bATHV39epVmJiYKGwPPTs7O1hYWCAgIEAh4zHKLzs7G0VFRZyuE81qwHKUk5MDPT09/Pbbb+K2LVu2oE2bNnj48CFmzZoFoHRKVIsWLXiKUjW4urqKk+/IkSMxZMgQmfqjlOLIkSO4e/cuAGDSpEkS602cOHECO3furFSfHh4eUh++YVQfpRSNGjVS2AJNKlEDVnaTJk1Cfn6+xLKGBgYGcHFxQZcuXeDpWbr8saWlJVxcXHiKUjXExcUhKSkJLi4usLOzk3keMCEEEyZMwLfffgtHR0e0adNGYiv6Zs2aVbrPDh06sGlr1RQhBFu2bFHIriwAK0Fw4tKlS0hPT5eoAV+9ehUPHjzA8OHDUVRUhCZNmvAYoeqYNm0aDh06hMzMTxffq7ro6OhyF9C5d+8e7t69i4kTJ3I2HsN8qrwSBCilSv8C4AVgh5WVFVUVc+fOpbq6urR///4UAF28eDHfIamE6OhoeuPGDYWN5+fnRwFQoVBIHz58SNu0aUOvXbumsPEZ5fPu3Tt6//59WlRUxFmfAO5QKblNJWrAlNIzlNLvuJ4czZXMzEwQQrB+/Xpx2/Lly+Hm5ibe7nr9+vVYtWoV233hC1q0aCFeP3n48OEYMGCAzH0GBQXh2LFjAIBvv/1WYvWz6dOnIyUlBWpqarh58yYiIyMl1g+WxsXFRaa96hjldurUKdjb2yMuLk7uY7EaMAfmz58PABI32LS1teHu7g4XFxekp6ejbdu20NPTk9i2iCkrJSUF9+/fh6urKzp27AihUChzn3/88QciIyMxZMgQODk5Sczh/bgs0alTJ0ycOPGLy1O6urrC3Nxc5rgY5dS9e3ccPXpUIet9sBowB65fv46XL19K1IBv3ryJoKAgjBkzBgUFBQor6qs6f39/fPPNN4iOjuZsxkhaWhrU1NRgampa5r2kpCScOHECgwYNYkmVkRs2DU2OunbtKpF8AeDatWv46aefMHv2bLRs2RJTpkzhKTrV4unpibCwME5nGdSrV09q8gWAFy9eYPr06YiMjAQAODo64ueff+ZsbEY1RUVFiX8m5IklYA68efMGhBCsWbNG3DZ37lz0798fISEhAIDt27dj+fLlbBnDL6hbty46dOgAHR0dDBo0iJOSTUREBDZu3IiioiKMHj0a3bp1E7/n5OSE9PR09OzZE+fPn8e9e/e+uEB7p06dZNqrjlF+I0aMwOLFi+U+jkrUgJV9LYjVq0s3gf746S0NDQ14enrCw8MD6enpaNOmDfT09CSWrGTKevfuHYKCgmBvbw8PDw+Z5wEDwI0bNzB79myMGDECLi4usLGxEb+nqakpvjo2NzfHtGnTMGPGjM/2179/f4l1P5jq5/fff1fIErKsBsyBu3fvIiYmRmLlrjt37iAgIAAjRoxAcXGx1PVombLi4uJgZWWF/fv3y7wbxgd5eXkoKiqCkZFRma1mKKXYuHEjHBwc2HoQjNywGrAcOTo6llk28f79+1i5ciUWLVqEdu3aybysYk3RqFEj3Lt3r9JLRH6Onp4e6tSpI3WfL0IIlixZgjNnzgAApkyZwp5WZJCamoqzZ88iLy9PvgNJmxysrC9HR0fOJkZzKS0tjQKgq1atEreJRCLar18/Wr9+fQqAAqA//fQTNTY25jFS1dK/f3/ao0cPmftJSUmha9asobGxsXTEiBG0Y8eOEu/n5uZSoVBIz549SwFQb2/vz/bn4OBAR48eLXNcjPI6duwYBUAfPHjASX8o50EMlagBK7tt27YBAOzt7cVthBD069cPvXv3RkZGBlq1agVdXV1286YCjh07BisrK/Tr14+TecCpqamYN28erKys4OnpiTdv3ki8b2BgAABo0qQJZs+ejdmzZ3+2v2HDhrGbqdWcq6srwsLC5L44O6sBcyAqKgpRUVEYMmSI+Nfc+/fv4/Dhw/D19YVAIECnTp14jlJ16OrqYurUqVi7di0n/QmFQhQWFkJXV1dqGWL//v0ASnfgYBh5YDVgOWrdujWGDRsm8Zc7OjoamzZtwurVq+Hs7CxeEY35srt374qfLuSChoYG9PX1pSZfoDQB79mzB8D/3/3Oz8/nbHxGNZ07dw5yv+CTVpdQtheUfDGe1NRUCoD6+flJtHt6elIzMzNxDXjx4sVUT0+PpyhVT69evWj37t1l7kcoFNKVK1fSq1ev0qFDh1J7e3uJ9wUCARWJRDQgIIACoC4uLvTdu3fl9mdra0tHjBghc1yMcqtXrx797rvvOOkLqlwDppSeAXDGyclJKdcMPHjwIACU2cVh2LBhEAqFyMzMRMuWLaGrqwtCCB8hqpRz585BS0sLQ4cO5aQGrKamhp9++gkLFiyAt7d3maUuNTU1AQBWVlZYuHAhpk2bBj09vXL7GzduHOrXry9zXIxyCwoKkvtWYqwGzIG4uDhERkaiT58+4r/Mjx49wp49ezB8+HAUFhaiW7duCtnkrzpo37496tati/Pnz3PWp0AggJaWltT3Lly4gLCwMPj5+XE2HsN8jNWA5ahZs2YYMGCAOPkCwMuXL7F7925s2bIFbm5uaNWqFY8RqpaTJ0+Kb4xxpbzkCwAhISHYuHEjgNKn5vT09BAcHMzp+IzqefjwofjegLywBMyBtLQ0EEKwaNEicVu/fv3Qvn17hIaGAgBiYmKwcOHCzyYCppS5uTnq1q2LHj16oGvXrpz0uX37dhw4cACDBw9G27ZtJd5bvXo1cnJyEBAQgK5du6JTp06ffdTYxsYGw4cP5yQuRnmdOnUK3377LSePw5dHJWrAyu706dMAILHGAACMHj0aJSUlyMrKEm+Lrqury0eIKuXq1atITk7G6NGjOakBA6UzHerVq4eRI0eWqQF/qMtbW1tj6dKlmDRp0mfn+U6dOlXutUGGf1OnTsXEiROhoSG/NMlqwBxITExEZGQkXF1doaOjAwB48uQJtmzZguHDh0MgEKB79+7s6reCxo0bh6CgILx8+ZKzPiml5d4AvX//Pg4cOIBFixaJN+z83PEMU1msBixH5ubm8PDwkKgBp6am4u+//8auXbvQq1cv6OvrQygUoqioiMdIVcP69etx584dFBcXc/br34dkKq3P+Ph47N69G+np6cjMzIS6ujrWrVtXbl9FRUWcXZkzyiszMxO///47YmNj5TYGS8AcSEtLg5aWlsT6oa6urnBwcMDt27cBlD6NtWTJEvFjr0z5jI2NUa9ePfTp0wdubm6c9Pn3339j1apVGD58OBwdHSXe8/HxQW5uLmJjY2FiYoJevXp99slFOzs7jBw5kpO4GOWVmZmJqVOn4tatW3IbQyVqwMq+HvClS5cAoMwWOuPHj4dIJEJWVhasrKygo6OjkDVGVV1YWBjCw8Mxfvx4zq40g4KCcPv2bSxbtqzcLe9tbGywcuVKjB8//rPzfGfOnIl69epxEhejvCwsLJCSkgITExO5jcFqwBx4/fo1oqKi4OzsLL7J9uzZM/zyyy8YOnQoiouL4e7uzm7AVZCfnx+WLVsGkUikkDrsixcvsGbNGkyePBm2trYQiUQQCoUSJSWGkQWrAcvRV199hQ4dOkg8aJGdnY3jx49j9+7d8PLygp6eHgQCAd6+fctjpKph7ty5yMzMREFBAedrMhQUFJRZ4zU3NxeHDx/G8+fP8fbtW5iZmX12D7+8vDwUFBRwGhejnHbv3o1z587JrX+WgDmQmpoKQ0ND/Pjjj+I2JycnODg44NGjR+K2pUuXiu+yM+XT19dHnTp14O3tzdkiRhcvXsS0adMwcuRIODs7S7xna2uLjIwMUEphaGiIXr16YcCAAeX25eTkhLFjx3ISF6Pc1q5diwMHDsitf5WoASu7a9euAUCZrecnT56MkpISZGdnw9LSEtra2mwNgQp49OgRLl68iJEjR0JbW5uzPv/55x/88ccfePfundRjWrdujXXr1uGbb7757DzfefPmsX9Ia4iwsLDPrgsiK1YD5kB6ejqePn0KBwcH8TcrISEBS5YswZAhQ1BSUoKePXuyGRAVtHfvXowfPx4JCQlo0qSJ3McrKCjAjBkz4OPjg169eqGoqAjFxcVy/YvH1CysBixHpqamaNWqlcQNo/z8fFy4cAH79++Hj48PDA0NkZeXh4yMDB4jVQ0jR47E27dvUbt2beTk5HDa99u3b5GdnS3RRgjBmTNn8PTpU2RkZKBnz57o3bt3uX1kZ2ezWn4N8e+//2LNmjVy658lYA68fv0apqamWLJkibitVatWsLOzQ0xMjLht+fLlMDMz4yNElaKpqQl9fX34+PigX79+nPQZFhaGiRMnok+fPmU23dTW1kZKSgosLCxgamoKZ2dnTJ48udy+nJ2dMWHCBE7iYpRbYGAgVq5cKbf+WQ2YAx8etvh0HvC0adNQUlKC3NxcWFhYQFtbG40bN+YjRJUSHx+PQ4cOYeDAgZzVzFNSUnDu3DnMnz+/3Pqtra0tNm/eDF9f38/WeBcvXgxTU1NO4mKU26pVqzjbGksaVgPmQFZWFmJjY9G6dWtx3TA5ORmzZs3CwIEDxTVgtoBLxQQFBaFHjx4ICQlR2Bbx06dPh52dHb799lvx9Dd5TsBnahZWA5aj2rVrw9zcHCKRSNxWVFSEK1euwN/fH6NGjUL9+vWRk5OD5ORkHiNVDa6urigqKkKrVq3K7GAsq8zMTKSlpZVp/++//xAZGYnk5GRMnz69zMp2H0tLSyv3aTqmeomJicGCBQuQmJgol/5ZAubA69evYWZmJlEDtrCwQNu2bfHq1Stx26pVq2BpaclHiCpFXV0dtWrVgq+vLwYOHMhJn3FxcRgzZgzat2+PHj16lHn/1q1bcHV1hZmZGVq1avXZul/37t3x/fffcxIXo9xevXqFDRs2cLoy38dUogas7GtB3Lt3D0DZecCzZs2CSCRCbm4uGjduDC0trTJ1Yqas9PR0bN26FR4eHpztJJKXl4eQkBD4+vrCyanMb4IAShfZ+fPPPzF48ODP1niXLVsGY2NjTuJilJurqysEAoHcHolnNWAO5ObmIj4+Hs2bNxfXgDMyMjBhwgR4eXmhpKQEPXr0gIWFBVtjtgKePXuG5s2b4+DBgwpbdWzZsmWglMLPzw/5+fnIzs5G/fr1y93KnmEqg9WA5cjAwAC1a9eWWGdWKBTi5s2bOHr0KL777jtYWlriwIEDePHiBY+RqoZmzZpBJBLBw8MDKSkpnPadlpaGpKSkMu3Pnj1DVFQUEhISsGnTJpiZmWH27NnYt29fmWOTkpKk1pGZ6qewsBBz5sxBUFCQfAaQtle9sr4cHR2pMkpOTqYA6NSpUyXa3d3dqb29PQUgfmlpafEUperx8PCgXbp04aSvN2/eUF9fX2pkZERtbW2lHnPq1CkKgO7bt4/++eef4u/Zp1q2bEmHDRvGSVyMchMIBFRPT4+uX79epn4A3KFScppK1ICVXWRkJACgZcuWEu3z5s0T7wl34cIFuLq6sm2JKiAvLw/Lly9Hly5d0LFjR076LCkpwd27dzFixIhyn3Kzt7fH3r170b9/f5iYmKBWrVpSb5quWrUKderU4SQuRrlpamqWu3YIF1gNmAN5eXl4+fIlmjRpIl7zNysrCyNHjsSAAQMgEAjg4eEBGxsbVlOsgMzMTDRs2BAbNmz47LKQXFqzZg0SExOxZcsW5OXlIS0tDebm5qhVq5ZCxmeqN1YDliNdXV2oq6tLrBFLKcX9+/cREBCAGTNmoE2bNvjtt9/kur9UdWFsbIzCwkJ4e3tzPv0nOTkZCQkJZdrT09MRHR2NmJgYHD58GJaWlhg9ejTWr19f5tiEhAQ2n7sG+emnn6TeC+CEtLqEsr6UvQY8efJkiXZ3d3fq5OTEasBVxGUNWCAQ0EGDBlEAX6wBHz16lO7bt4/VgBlKKaXt27en33//vUx9gNWA5Sc6OhoAyjw9tWjRIohEImRnZyMkJATOzs4Su2Yw0olEIsydOxeOjo5SH5qoCjU1NcTGxmLUqFHw9fWVeoyjoyP8/f3h7u4OY2Nj6OrqolGjRmWOW7t2LdvbrwYJDw+XW9+sBsyBgoICJCcno2HDhtDR0QFQOjd48ODBGDRoEN69ewc3Nzc4OjqyBFwBlFIYGRlh0aJFWLBggULG/P333xESEoJ//vkH+fn5SExMRKNGjdg+fgwnyqsBsytgDmhrayM/Px+5ubniBAwAjx8/hqamJi5cuAAAWLhwIb7++mu0adOGr1BVAiFE/HBLXFwcmjVrxlnfL168QFFREZo3by7R/uFG6qNHj5CQkABvb2+4ubmhY8eOWL16tcSxMTExbGW7GmTLli3IyMiAn58f951Lq0so60tZa8ApKSkUAJ0wYYJEu7u7O+3QoQOrAVcRlzVgSikdOnRohWrA58+fp/7+/qwGzFBKKR0/fjzt1auXTH2A1YDlJz4+HkDpnmIfW7p0KUQiEXJychAeHg4HBwc2Da2C5s6di1atWmHw4MGc9ZmWloYhQ4bgu+++k/p++/btcezYMTg7O6N27dowMjJCgwYNyhy3adMmGBkZcRYXo9x2794tt75ZDZgDAoEA6enpMDU1FW8iWVBQAE9PTwwaNAhv376Fi4sLunXrBg0N9m9eRTRv3hxDhw7FqlWrFDLegQMHsHfvXly5cgUCgQCxsbFo1KgRu9nGcILNA5YjTU1NpKamSqwRSwjB8+fPERQUBD8/P7i7u2PKlCnildOYz4uNjcXYsWPx9OlTTvuNi4vD48ePy7SLRCJkZGTgv//+w8OHD9G2bVu4uLjg22+/LXNsVFSU+Lcepvo7cuQIRo8eLZ/OpdUllPWlrDXg169fUwB07NixEu3u7u60U6dOrAZcRVzXgEePHl2hGnBwcDA9evQoqwEzlFJK169fT62tralIJKpyH1DlGrCyrwf8YbX8T2c3rFixQrwe8N27d2FnZ8dqwBW0cOFCWFhYYNy4cZz1WVhYiP79+2POnDlS3+/YsSNOnz4Ne3t7GBkZITAwUOo2Ur///jsMDQ05i4tRbrNnz8bs2bPl0jerAXOguLgYWVlZMDIyEi+2U1RUhG7dusHHxwdZWVlwdnZG7969oampyXO0qqFz585wcnLC5s2bFTLeqVOnsHbtWpw5cwYGBgaIiIiAubk528eP4QSrAcuRhoYGnj17JrFGrJqaGpKTk3Ht2jX88ssvGDBgAHx9fXHr1i0eI1UdN2/exKRJk/Do0SNO+3369CkePHhQpl1DQwOFhYW4fv06EhMT4eTkhHbt2qFfv35ljr1//7746Uem+gsNDcXgwYM5X5saAKsBcyE1NZUCoN98841Eu7u7O+3cuTOrAVcR1zXgyZMnV6gGfPv2bXr69Gnx90woFEocx2rANcv58+dpmzZt6LNnz6rcB1S5BqzsPvzL2LZtW4n2X3/9FSKRCO/evcPDhw/RunVrVgOuoKVLl6Ju3bqYNWsWZ32qq6ujZ8+e5T7R5OzsjMDAQNjY2KBjx464evUqTE1Nyzw+vmvXLujr63MWF6Pc+vTpgz59+silb1YD5kBJSQnevn0LPT098fqxIpEIDg4OGDp0KDIyMtC+fXsMGjRIPE+Y+byBAwfCxMRErpPgPxYcHIx58+bh0KFDsLKywq1bt2BmZoYmTZooZHymemM1YDlSU1PDgwcPJPYaI4QgKysL4eHh2LhxI77++mt4enri2rVrPEaqOk6dOoXp06fj/v37nPYbGRmJ//77r0y7lpYW1NTUEBoainfv3qFLly5o0aIF3NzckJubK3FseHg4oqKiOI2LUV7x8fHw8vLCzZs3ue9cWl1CWV/KWgNOT0+nAOjw4cMl2t3d3WmXLl1YDbiKuK4Bz507t0I14Hv37tHAwEDx9ywlJUXiOFYDrlni4uKovb09vXTpUpX7AKsBy096ejoAoF27dhLt69evB6UU7969w+PHj2Ftbc1qwBW0atUq6OjoYO3atZz1qaenBxcXF6xbt07q+126dEFwcDCsrKxgb2+PGzduwMTEBPXr15c4bv/+/dDT0+MsLka5WVpayu0JVlYD5oBIJIJAIICmpqbEDZu2bdti0KBBSE9PR7t27TBy5Ej2F7eCJkyYgLy8PPz9998KGe/evXsYO3YsduzYgU6dOuHmzZswNjYus9Eqw1QFqwHLESEEoaGhZfYay8/PR0REBLZv345JkybByckJly9f5idIFbNr1y7MmTOH890IHjx4ILWWp6urCwMDA1y7dk28mP6H2RAfnnT84MaNG3j48CGncTHKSygUolevXvjrr7+471xaXUJZX8paA37z5g0FQAcPHizR7u7uTrt27cpqwFXEdQ142bJlFa4B37hxQ/w9i4yMlDiO1YBrFpFIRJ2dnemuXbuq3AdYDVh+srOzAQAODg4S7Zs3b4ZIJEJ+fj5iYmLQrFkzEEJ4iFD1rF+/HkKhENu3b+esTxMTEzg7O+P333+X+n63bt1w48YNNG/eHPr6+ggPD4exsXGZHTkOHTrEtiqqQQgh8pkBAbYlEScsLS0hFArLJNdhw4ZhwIABSE1NRWRkJIRCIQIDA3mKUrWkpqbCyMgItra2nPU5depUTJ06Vep7z549Q//+/bFx40bo6+sjNDQUhoaGaNy4MXbt2oXo6GgMGjQInTt3hr29PWcxMTWctMtiZX0pawlCJBLR06dP06dPn0q0W1tbi7dC//CaPXs2T1Gqnps3b9Lr169z2md4eDgNDg4u056YmEi7d+9OFy5cSLOzs2mLFi0oANqoUSPx987NzY1SSumVK1fonTt3OI2LUW5eXl5048aNVf48yilBsJtwHMjMzIS3tzfmz58v0W5mZiaeovbBhg0bFBmaSluyZEmZ/6eyWL58OTp06IDp06eXec/c3ByzZs3C6tWrER8fj8OHDwOAxA04c3NzAMCUKVOwZs0azuJilB8hRC7lQ1aC4MCHJ6WcnZ0l2rdv3w5KKfLz83H69GmkpaVhxIgRfISocn755RcIhULs2rWLsz7r168PV1dXbNu2Ter73bt3x507d2BtbQ1dXV08fPgQ+vr6yMrKwps3b2BjYwMAOHbsmMTu10z1FxAQIJd+WQLmQNOmTUGlzKfu168fBg8ejKSkJNy5cwdCoRBLly7lIULV8/btW9StW5fTebgTJ07ExIkTpb4XGxuLXr16YcuWLXB0dERQUBBq164Nc3NzXLhwAU+ePIGOjg4aN25cZvNVhqkyaXUJZX0pcw342LFjNCoqSqLdxsamTA146tSpPEWpeq5duya1XiuLmzdv0n///bdM+6tXr6iHhwf93//+R7OyssQ1YGNjY/H3rlu3bpRSSgMDA+nt27c5jYtRbp6ennTt2rVV/jxYDVh+MjMzMWTIkDL1ygYNGpSpAZf36y9T1rJly/Djjz9y1t/SpUvRuXNnqVsSmZmZ4YcffsDatWvx/PlznDhxAgAkNlr9sCXWjBkzWC2/hjEyMpJL2YmVIDjw9u1bAICbm5tE+86dO0EpRWFhIS5duoSMjAx4e3vzEaLKWbFiBYRCIfbv389ZnxYWFujbty9+++03qe+7ubkhIiICVlZW0NHRwePHj6Gnp4ecnBxkZmaiefPmAErrgWxZ0ZrlyJEjcumXJWAOWFhYSK0Bu7u7w8fHBy9evMDNmzdRUlKCadOm8RCh6ikuLkbDhg3LPAQhi3HjxpW7yWdMTAzc3d2xfft22Nra4uLFi6hTpw6aNWuG48ePIzIyEt988w0aNmyIFi1acBYTU8NJq0so60uZa8D+/v40IiJCot3GxoYOGDBAogY8fvx4nqJUPZcvX5Zar5VFSEgIPXfuXJn2pKQk2qtXLzpz5kyakZEhrgFraWmJv3fOzs6UUkrPnDlDQ0NDOY2LUW5ubm509erVVf48yqkBq9RqaISQdAAvKni4KYA3cgyHb+z8VFd1PjeAnZ80TSildT9tVKkEXBmEkDtUyvJv1QU7P9VVnc8NYOdXGWwWBMMwDE9YAmYYhuFJdU7AO/gOQM7Y+amu6nxuADu/Cqu2NWCGYRhlV52vgBmGYZQaS8AMwzA8qZYJmBDSmxASTQh5RghZwHc8siCENCKEBBNCnhBCogghM963GxNCLhFCYt//WYfvWGVBCFEnhNwnhJx9/3W1OT9CSG1CyDFCyNP330fn6nJ+hJBZ738uIwkhfxNCtFX53AghewghaYSQyI/ayj0fQsjC93kmmhDSq7LjVbsETAhRB7ANQB8ArQCMIIS04jcqmQgBzKGU2gDoBGDq+/NZACCIUtocQND7r1XZDABPPvq6Op3fJgAXKaUtAdih9DxV/vwIIWYApgNwopS2AaAOwBeqfW77APT+pE3q+bz/e+gLoPX7z/z+Pv9UnLTH41T5BcAZQOBHXy8EsJDvuDg8vwAAPQFEA2jwvq0BgGi+Y5PhnMzf/2C7Azj7vq1anB8AQwDP8f6G90ftKn9+AMwAJAIwRum6MmcBeKr6uQGwABD5pe/Vp7kFQCAA58qMVe2ugPH/PxQfvHrfpvIIIRYA7AGEAfiKUpoCAO//rMdjaLL6DcA8AKKP2qrL+VkCSAew932JZRchRA/V4PwopUkA1gF4CSAFQA6l9F9Ug3P7RHnnI3OuqY4JWNrGTSo/144Qog/gOICZlNJcvuPhCiGkP4A0SuldvmOREw0ADgC2U0rtAeRBtX4lL9f7WugAAE0BNASgRwj5ht+oFErmXFMdE/ArAI0++tocQDJPsXCCEFILpcnXn1J64n1zKiGkwfv3GwBI4ys+GXUB4E0ISQBwGIA7IeQgqs/5vQLwilIa9v7rYyhNyNXh/HoAeE4pTaeUFgM4AaAzqse5fay885E511THBPwfgOaEkKaEEE2UFslP8xxTlZHSrVh3A3hCKf14G4bTAMa8/+8xKK0NqxxK6UJKqTml1AKl36srlNJvUH3O7zWAREKI9fsmDwCPUT3O7yWAToQQ3fc/px4ovcFYHc7tY+Wdz2kAvoQQLUJIUwDNAYRXqme+C95yKqL3BRADIA7AYr7jkfFcuqL015oIAA/ev/oCMEHpjavY938a8x0rB+fqiv+/CVdtzg9AOwB33n8PTwGoU13OD8AyAE8BRAL4C4CWKp8bgL9RWs8uRukV7refOx8Ai9/nmWgAfSo7HnsUmWEYhifVsQTBMAyjElgCZhiG4QlLwAzDMDxhCZhhGIYnLAEzDMPwhCVghmEYnrAEzDAMw5P/A4EDIZXn7h5oAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "pl.semilogy(E, flux[0] * E * erg_per_ev, color='k', ls=':')" ] }, { "cell_type": "markdown", "id": "ffc1c049", "metadata": {}, "source": [ "You should see the characteristic sawtooth modulation of an intrinsically flat spectrum.\n", " \n", "Compare to the analytic solution, given by Equation A1 in [Mirocha (2014)](http://adsabs.harvard.edu/abs/2014arXiv1406.4120M>) (the *cosmologically-limited* solution to the radiative transfer equation), which does *not* take into account the sawtooth modulation:\n", "\n", "$J_{\\nu}(z) = \\frac{c}{4\\pi} \\frac{\\epsilon_{\\nu}(z)}{H(z)} \\frac{(1 + z)^{9/2-(\\alpha + \\beta)}}{\\alpha+\\beta-3/2} \\times \\left[(1 + z_i)^{\\alpha+\\beta-3/2} - (1 + z)^{\\alpha+\\beta-3/2}\\right]$\n", "\n", "with $\\alpha = \\beta = 0$ (i.e., constant SFRD, flat spectrum), $z=10$, and $`z_i=40$," ] }, { "cell_type": "code", "execution_count": 7, "id": "8dfb2f65", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEMCAYAAADK231MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABalUlEQVR4nO2deVxU5ffH3w8iioILoOIOCu67YOa+a+5marmUWdpqm6apbZaZlfmrzFSyNDW3tNzN3dz3XdxwAREVBGRT9uf3BzhfRgaUYWbuzPC8X6/7cua5957ncwXO3Dn3POcIKSUKhUKhsF8ctBagUCgUCvOiHL1CoVDYOcrRKxQKhZ2jHL1CoVDYOcrRKxQKhZ3jqLUAa8TDw0N6eXlpLUOhUCiemGPHjt2VUpYxtE85egN4eXlx9OhRrWUoFArFEyOECM5pnwrdKBQKhZ2jHL1CoVDYOcrRKxQKhZ2jHL1CoVDYOcrRKxQKhZ2jHL1CoVDYOcrRKxQKhZ2jHL2JCAoKYurUqdy8eVNrKQqFQqGHcvQmIjAwkEmTJnH79m0AVJ1/hUJhLShHnwUhRC8hREBMTEyez+3VqxdJSUk0atSIqKgo2rRpw5YtW8ygUqFQKPKGUHee2fHz85N5KYGQnJzMpEmTdO+TkpJYtWoVfn5+1KhRwxwSFQqFHTNo0CD8/PzydI4Q4piU0uBJqtaNCUhNTeWXX37JNr5t2za2bdumgSKFQmHLNG7cOM+OPjeMdvRCiJpANcAZiABOSCnjTSXMlihWrBgJCQkG961cuZILFy7w8ccfW0xPSkoKQggcHdXnuEKhyGOMXgjhJYT4VggRBgQCG4CVwH9AlBBihxBioBBCmEGrTbJjxw7WrVtHcnKyxeZs3LgxY8aMsdh8CoXCunliRy+E+A44C9QEJgL1gJJAEaA80B3YD0wDTgohmphcrQ3y/fffs2/fPpycnCw2p5+fH7Vq1QLgv//+o3v37oSEhFhsfoVCYV3k5bu9K1BDShlmYN+dzG0b8LEQYgBQGzief4m2jbOzMwD379/n5MmTtGjRwuxzLliwQPf67t27hISE4OHhAcDFixfx9PSkZMmSZtehUCisA5V1Y4C8Zt08Ca+88gorV64kJCTE4k5WSsnDaNrTTz9NcnIyx44ds6gGhUJhXnLLulF59BZiwoQJrFu3ziJOvly5crz//vu691kfmfz8889MmzYNgLS0NMaOHculS5fMrkmhUGhHvtIyhBAjpZS/mkqMPePj44OPj49F5nr77bdp1KiRwX1NmzbVvT579iyzZ8+mefPmKt9fobBj8hW6EUKESCmrmFCPVWCO0M1Dpk6dSmxsrO6uWmsiIiJwd3fHwcGBX3/9lV27djF37lxcXFy0lqZQKPJAvhZMCSFO57QLKJcfYQWR0NBQoqOj9eLmWlKmzP+axsfGxhIZGUnx4sUBSE9Px8FBRfcUClvnsXf0Qog7QFcg+tFdwH4pZQUzadMMc97RW8J5urm5MXToUH766ac8n/vwAyguLg5/f3++/PJLBgwYYAaVCoXClOS3BMImwFVKedKA4b351FbgeOjkIyIicHV1pWjRoiafY/z48TRo0MCocx9+y4iJiaFatWpUrVoVUHf3CoUto9IrDWDOO3rIyGVv2LAhv/zyCyNGjDDbPKZk0qRJXL58mSVLlqjSCgqFFWLS9EohRJ/8SyrY1KhRg4kTJ9KqVSuz2E9OTiY1NdWkNkuWLImbm5ty8hYkNjZW9TVQmARjvot/bXIVBQwhBJ9++qnZUho9PT354IMPTGpz3LhxzJkzB4Br167RrVs3rl27ZtI5FP/j0qVLlCxZUtfT4O7du0RFRWmsSmGrGOPotU8VsRPOnDnDunXrTG73448/pmfPnia3+5CgoCAuXLhAoUKFzDZHQcfZ2Zk+ffro1l7MnDkTT09P4uMzCsRGRkaSlpampUSFDZHnGL0QIlBKWcdMeqwCc8foH9K3b1+OHz9OcHCwVaRa5oWUlBQKFy4MwA8//MDAgQOpUMHuErCshlOnTnHo0CFGjRoFwPPPP8+5c+c4c+YMoP/zUBRMVOMRK+X777+nZMmSJnfycXFxFC5c2CwZPQ956FRu3LjBpEmTiI2N5dNPPzXbfAWN1NRUHjx4gLOzM46OjjRs2JCGDRvq9g8ePJiIiAjd+xYtWtC8eXNmzpyphVyFlaPy5TSkevXquqqSpqRq1aqMGzfO5HYNUblyZU6dOsXEiRMBuHXrlgopmICDBw9SokQJdu3aZXB/7969eeWVV4CM1Nfu3bvz1FNPARl3982bN2fVqlWWkquwcoxx9HEmV1GA2b9/P6+88opJneOXX35J3759TWbvcfj4+ODo6EhiYiIdOnRgyJAhFpvbXvH29mb69On4+vo+9lgHBwcmT57M0KFDAQgPD8fFxYUiRYoAEBYWxptvvsnVq1fNqllhvag8egNYKkYPsGLFCkaPHs3+/fupXr26ReY0J/Pnz6datWq0bdtWaymKTDZt2sRzzz3H0aNHqV27NpcuXeL27du0atVKLYKzI3KL0StHbwBLOvqUlBQcHBxMmsESFRVFkSJFdDVrtGLBggV4enrSrVs3TXXYIikpKcTGxuLq6mqS7mT379+nWLFiALz33nvMnTuX8PBwXF1diYyMpHTp0srp2zimXjBVWAhxRwhRN//SFIULFzZ5mqKPjw8TJkwwqc28kpaWxi+//MLPP/+sFv0YwaFDh/Dw8GD37t0msffQyQNMmTKFLVu24OrqCsCwYcNo3769SeZRWCd5dvRSyhQgLXOzOoQQfYUQvwoh1gghumSO1RZCzBFCrBRCvKG1xkfZt28fTZs25fr16yaxN23aNPr3728SW8ZSqFAhduzYwZ9//okQgqSkJOXw80C1atX46aefqFmzpsltu7i40Lp1a937F198UVeKQ0rJ4MGDzbK+Q6Edxn5Xmwe8Y0ohAEKI34UQ4UKIs4+MdxNCXBRCBAkhPsrNhpRytZRyJDAcGJQ5dl5K+TowEDD41UZL3NzcKF68ONHRjxYINY5Ro0ZZRYzcxcWFkiVLkpKSQr9+/XjnHZP/ytgtFSpUYPTo0VSuXNnscz3//PO89NJLQMYK3NOnT+tSNxMTE/n3339VJpWNY2wefQVggBCiA3AMSMi6U0o5yki7C4CfgYUPB4QQhYBZQGcgFDgihFgLFCJ7OYYRUsrwzNcfZ5730E5v4KNM+1ZF7dq1TfYVHTJSHIsVK2Y1DcAdHR1p0KCBXTxsthRJSUlERkbi7u6uy56xBGXKlOHMmTOkp6cD8M8//zB48GB27Nihwju2jJQyzxuwM5dthzE2s9j2As5mef80sDnL+wnAhFzOF8A3QKcc9m/IYXwUcBQ4WqVKFakFaWlpMi0tLd92SpcuLUePHm0CRebh8uXLMiUlRWsZVs2ePXskILdu3aqpjqSkJLl69Wrd7+V3330nn332WZmUlKSpLkV2gKMyB79oVOhGStk+l62DMTZzoSJwI8v70MyxnBgNdAKeE0K8DiCEaCeE+EkIMRfYaOgkKWWAlNJPSumXteuSpdixYwdubm66Je354f/+7/8YNGiQCVSZnvDwcJo1a2axBV22io+PD3PnzqVWrVqa6nBycqJPnz66jJxChQrh6OioywTatm2byUKOCjOS0ydAbhsZYZSWWd6PBI6QEXpxNcZmFlte6N/RDwDmZXk/DJiZnzketzVt2tTYD1WjCQkJkaNGjZLnz5+3+NyW5pdffpHXrl3TWoYin8THx0sXFxf58ssvay1FIc1wRw98C3gACCFqkBELP0rGg87vjLSZE6FA1idSlYAwE8+hOZUrVzbZHVxISAiRkZEmUGUe3njjDby8vAByXOJf0ElMTOT69es8ePBAayk5Urx4cfbt26crfxEcHEyvXr24cOGCxsoUj2Kso68OPMyM6Qdsk1K+Qcadvanr4x4BfIUQ3kIIJ+B5YK2J5wBACNFLCBEQExNjDvNPxJ07d/Jto1GjRkyePNkEaszL6tWrad++vUrlM8DRo0fx9vZm3759WkvJlQYNGuhKKV+6dIkTJ07oFupFRUWZvAGOwjjysxTuYVJ0W2BL5uubgLuxBoUQS4EDQE0hRKgQ4hUpZSrwNrAZOA+skFKeM152zkgp10kpR2mVrTJv3jw8PT0JDQ3Nl52ff/7ZJurN9OrVi/nz59OjRw+tpVgdvr6+zJ8/nzp1bKcieOfOnQkODtalhI4ePZqmTZvqMngU2mFUCYTMpuAHgHVkOPkmUspAIURLYKmUsoppZVoWS5ZAyMrly5dZt24dL774olmqWloz0dHRhIWFUbeuWnBtL6xfv54bN27wxhsZaxRXrVpFt27dNC/NYa+YvNaNEKIVsAYoBfwmM/PmhRDTAF8ppbbLMvOJVo7eVAQFBVGiRAnKli2rtZQnpkuXLgQFBXHx4kXVQAN48OABN27coGLFinbhGC9cuEDt2rX59ttv+fDDD7WWY5eYtNYNgJRyL1AGcJf6i6N+Bd4yxqYig/j4eM6dy19kqlmzZkyZMsVEiizDd999x8KFC5WTz+TYsWPUrFmTAwcOaC3FJNSqVYs9e/boOmTt27ePOXPmkJKSorGygoGqXpkFIUQvoJePj8/Iy5cva6Lh1VdfZc2aNYSHhxvdeWrlypV4eXnh52d11R6eiOvXr+uycgoq4eHhbNu2jQ4dOuDp6am1HJMzevRo/vnnHy5fvoyzs7PWcuwCVaY4j2gZujly5Ai3b9+mR48eBbJs7M6dO+nSpQurVq2id+/eWstRmAkpJWFhYVSsWJH09HReeeUVRowYoVdsTZE3TB66UZgPf39/evXqlS8nHxgYSFiYbS41aNmyJePHj6dNmzZaS9GUhIQEzpw5Q3x8vNZSzIIQgooVMxa4h4aGsnPnTpNVb1VkRzl6K+TMmTOcOnXK6PNbtWrFtGnTTKjIcjg5OTFlyhRKlSpFenp6ga2aeOLECRo0aMDBgwe1lmJ2qlSpwvnz53UpwcuWLWPcuHEkJiZqrMx+UI7eCunXrx9fffWV0ecvWLBAV1/cVnnw4AFdu3bl+++/11qKJtSqVYuVK1dSv359raVYBGdnZ9232JMnT/Lff/+ZpLOWIgNjyxQDIIQYKaX81VRitCbLw1hNdSxYsIBy5coZfb49xLaLFi1KxYoVcXc3ev2dTePh4aF58xitmDZtGklJSTg4OJCQkMCbb77JJ598gtZ/l7ZMvh7GCiFCbH1xlCFsPY/+5MmTuLm5UaWK3f1oCgzx8fFcvHgRX19fSpQoobUczThw4ADPPPMMa9euLfDPbR5Hvh7GCiFO57CdAYy/7VTkSFRUFCtXrtR1+ckrHTp0YPr06SZWpR0bNmxg4cKFjz/Qjjh58iR+fn4cPnxYayma8vTTTxMSEqJz8gsXLsSWb8K04rF39EKIO0BX4NGi0wLYL6WsYCZtmqH1Hf2hQ4do3rw5q1evpk+fPnk+/99//6VSpUrUq1fPDOosT9euXYmLi2Pfvn1Gry2wNSIjI9m/fz/NmzdHi/4I1khycjJ16tShUaNGrFy5Ums5Vke+8uiFEAvIKHOwx8C+FVLKgSZRaUVo7egTExMJDAykdu3aajEJEBERQalSpdSqWQXR0dGkpaXh4eFBeHg4d+/etanCb+YkX6EbKeVwQ04+c5/dOXlroGjRojRp0sRoJ3/48GGuXr1qYlXaUaZMGQoXLkxycnK+K3vaCrGxsezfv5979+5pLcWqKF26tK7g3/jx42nRogValhW3FVR6ZRasoR79Q/bs2cPy5cuNOrdbt2788MMPphVkBTzzzDP079+fgrCa+/Tp07Rs2VLFo3Phm2++YeHChTwsK56cnKyxIuslz1k3Qog+Uso1ZtJjFWgdugF46aWX2LFjBzdu3Hj8wY+wY8cOKlSooHm/UVOzbt06HBwc6N69u93H6qOjozly5AhNmzYtsCmmeWHnzp2MGDGCDRs2FNhQjklr3QghAqWUdv0/aQ2O/vbt2zg5OeHm5qapDoXCFjh+/DiffPIJy5cvx8XFRWs5mmDqWjf2fStlJXh6ehrt5Pfu3YtW1TfNTVpaGj/88APLli3TWopZiYmJYdeuXURHP5rspjBEkyZN2LBhAy4uLqSmpvLNN99Ydb9dS2OMo7f/AKkVEBUVxQ8//EBgYGCez+3duzczZ840gyrtcXBwYNmyZWzYsEFrKWblzJkztG/fnmPHjmktxebYuXMnEyZMYPPmzVpLsRryVQJBYT4SExN5//33mTNnTp5jjmvXrs1XCQVrRgjB1q1bcXV11VqKWalfvz47d+6kYcOGWkuxOTp37syZM2d0bSkTExMpWrSoxqq0RWXdWCnly5fn7t27uo48AKtXr8bNze2xYZlWrVrh6+trboma8dDJx8bG2m11y5IlS9KuXTtKly6ttRSb5KGTv3jxIj4+PmzatEljRdpijKOPM7kKRTaEELi7u+tllxQqVIgyZcpQqFChXM/dsWMHFy5cMLdETTl79ixVqlRh9erV+bKTmJhInz59CAoKMo0wExEdHc2WLVuIjIzUWopN4+7uTtOmTalRo4bWUjQlz45eSvmUOYRYA9aURw8ZLQFnzZqlex8WFsalS5ceu5Dqueee45dffjG3PE2pXbs2Q4cOzfcf8M2bN1m7di379+83kTLTcPbsWbp27crJkye1lmLTeHh4sGbNGqpXrw7A5s2bC8Q6jEdRrQQNYA3plQBDhgzh+PHjnD9/HshIubx69SpNmjTJNeZ4+PBhPDw8qFatmqWk2iyJiYkcP34cHx8fypQpYzX5+bGxsZw9e5Y6depQqlQpreXYBdu3b6dTp07Mnz+f4cOHay3H5OSWXomUMs8bcBWYBxR+ZNwDuGqMTWvamjZtKq2BBw8eyPT0dN37qKgoWaJECfnjjz9qqMq6uHHjhvzrr7/ybWfmzJmyU6dOev/fCvsiPT1dLl++XKampmotxSwAR2UOPs3Yh7FeQDdgqxAi69OiQkBVI20qHqFo0aJ6d5gxMTGUKVOGmjVr5nrev//+y9mzZ80tzyr45ptvGDZsGLGxsUadn5CQwLp164iKiqJUqVJW06M1MjKSdevWGV2qWpEdIQQDBw6kUKFCREdH88knn5CSkqK1LItgrKOXQOfMfw8JIVTrFzNw+fJlxowZo2uavHnzZq5cuUKDBg1yPW/w4MEEBARYQKH2jB8/nnPnzhndnOPOnTv07t0bLy8v/vrrL6tJ2wwMDKR3796cPn1aayl2yfr16/n2228LzDoFo2L0Qoh0wBOIAuYAfYF+wCUgTEqZe1qIlWMtMfr9+/fTqVMnNm7cSLt27YiIiCAkJAQfHx9dISdDqA5TT05SUhJnz57Fy8sLd3d3bt26xf3793UP77QiLi6OS5cuFfgOU+YkODiYqlXtJwBhjhh9GlA2y/txwP3Mf9OMsWlNm7XE6NPT07PFjF977TVZsWJFjRRZJ3fu3JEjRoyQ//33X77spKSkyPLly8u+ffuaSJnCFti8ebOcMGGCzT+fwQwxer3UBCnlt8Bg4BMj7SkMIITQi9FfvXqVlJQU3nzzzVzPW7t2bYFKy3N1dWXLli1G5cLHx8ezatUqgoODcXR0JCAgwCraMN69e5dVq1YRHh6utRS7Z8uWLWzcuJGEhAStpZiPnD4BctuAtoCjgfE6wEvG2LSGDegFBPj4+OTjc9W0fPLJJ3L+/PlSSinnzJkjARkWFpbrOaVLl5ajR4+2gDrrwdhMiitXrkhA/vHHHyZWlD92794tAblt2zatpdg9aWlpMjY2VmsZ+YZc7ujzFKMXQgx+wg+PJXn5sLE2rCVGD9CsWTP8/f2ZNWsWkZGRhIWFUbFiRYoWLUqxYsUMnhMYGEipUqWoUMHu2vk+luTkZJycnPJ0/MWLF6lUqZKu3MDly5f58ccf+eabbyhevLi5pOZKQkICV69exdvbu8CW3bU0KSkpvPfeewwcOJC2bdtqLSfP5Bajz2tRs8WPvJdkL1ssAZt29NbE4cOHda/d3d0JDw/H3d2dJUuW8MILLxg8p6A2XnjxxReJjIzMU2VLJycn6tevrzd2584dfv/9dwYNGkTr1q1NLfOJKF68eDZdCvOSkJDAjh07qFixok06+tzIk6OXUurF9IUQcUBDKaX9NCi1Yi5fvsyePXv44osvaNSoUY7HrVy5Ei8vL/z8DD+At1f8/f3znAcfHx/P2rVrad68uW4lccuWLbl161aumU3mJjw8nG3bttGhQwc8PT0101GQKFWqFEeOHLHPb1A5xXSeZCOjwFm1/Niwxs1asm6klHLZsmVy5MiRUkoVozcHj4vRa5WJoWL02hIYGCjHjBljU5k4mCHrRmEhrl69yu7du5FSMmjQIC5evIiTkxN37tzJ8ZzDhw/z8ccfW1Cl9ZCenp6nyp2VKlXi4sWL9O3bV29cSkmfPn14//33TazwyWjatCkXL17k6aef1mT+gs7WrVtZuHAhwcHBWksxCfkqamavoRtrehhriCZNmlC+fHm777JkDB9++CG//PILEREROT6sflI++OADKlSowNixY02kTmErSCmJjIzEw8NDaylPjCkfxhpClb+0EBcvXuTAgQOMGTMm11/AJUuWUL16dZ56ym4rSufI0KFDadq0KQ4OT/ZlNS4ujlWrVtGqVSt8fPQrecyYMcMcEp+IO3fusGnTJrp06VIgs6e0RgiBh4cHUkrmz59Pz549KVu2rNayjCenmI6hDdjyyJYC7Hl0PC82rXGzphj9+fPn5XPPPSdPnTqlYvRm4Eli9JcvX7awKhWjtxZCQkJkkSJF5Keffqq1lMdCLjH6vN7R33zk/aPplgoTk5aWxtmzZ4mKimLw4MF07dqVYsWKERgYSO3atQ3WTz958qRm+d/WQEREBDt27GDgwIGPrS9fuXJlrl27luM3pE8++YQZM2Zw9+7dfIeC8oK/vz/Xrl2z296/tkLlypXZv39/rllutoBqPGIAa4/Rz5gxgzFjxhAdHa2aUhggICCA1157TfdhmB9OnjzJyZMnGTRo0GM7eynsm9jYWOLi4qhYsaLWUgySW4xeZd3YEIGBgQQEBNCuXTuWLl2a4wrQP/74g3379llYnfXQv39/jh079ti6/ZDxxxsQEMClS5cM7m/UqBHDhw+3uJO/desWAQEB3Lz56JdohRakp6fTsmVLXnrpJa2lGEdOMZ2CuGGFtW7S09Pls88+K+fPn69i9GbgSWrdxMTEyKVLl8rk5GSL6VIxeutj5cqV8sCBA1rLyBFMVeumoGBtoZsWLVowaNAgXn31VWJiYihdujRBQUFUqFABd3f3bMffunWLYsWKabqyU2vOnz/PihUrmDRpEo6OOT+KSk1NJTw8nFKlSuUYg1+7di19+vRh69atdOrUyVyS9UhKSiIyMhJ3d3eKFClikTkVT46U0mr6Cz8kt9CNcvQGsDZH/yjXrl2jWrVqdtvk2BQsX76cwYMHc+rUKerVq5cvWw8ePOD48eM0b96cQoVsuqeOwgRMnTqVK1eu8Ntvv2ktRQ8Vo7cTzp49y8yZMylRogR//fUX7du3N3hcQEAA//33n4XVWRe9e/cmOjr6sU4+JiaGmTNncv78+RyPcXZ2pmXLlhZ18mFhYcycOZMbN25YbE7Fk5GYmEhiYiKpqalaS3lycorpFOTNmvLopZRywoQJ8u2331YxejPwpPXob968Kb/44gsZHBxsEV0qRm+9WGv9G/Jb60YIUVMI0VEI4frIeE/Tf/QoHiUpKYnExEReeukl7t69S7ly5Thz5kyOmSJBQUF8/fXXFlZpfezatYuXXnqJtLS0HI+pWrUqd+/eZdCgQbnaiouL47PPPuPgwYOmlmmQ5s2bc/fuXdq0aWOR+RRPzsPY/I0bN9i2bZvGap6Mx8bohRBvAaOBi0Bj4F0p5T+Z+45LKZuYXaWFsfYYPYCPjw9PPfUUf/75p9ZSrJaFCxcyadIkDhw4QKVKlfJlS0pJVFSUwYffioJJt27duHDhAleuXLGKZzf5ehgrhDgNPC2lTBBCeAMrgT+llDOEECeklI1NL1lbrNXRnzp1im3btvHaa69x5MgR3N3dadCgQbbjZs2aRe3atenQoYMGKq2H9PT0x9a8uXfvHr/99hvdunWjbt26FlL2eEJDQ1m+fDnPPfccVatW1VqOwgAXLlygWLFiVKlSRWspQP4fxjpKKRMApJTXgHZAVyHEDLJ3l1KYgT/++IMWLVqwf/9+xo4dS1xcHO3btzfo5CFj2f7q1astK9IKeZLCZlFRUYwdO5Zjx4499tirV68ybNgwTp06ZQp5uXLt2jXGjh1rVMNzhWWoVauWzsk/7oZZa57E0d8WQjR6+EZKGQf0ADwA1evMAjg6OuLi4sKwYcOIjY2lXLlyBAUFsWfPHoPHBwcH8+2331pYpXXy008/5bqa0cvLi9jY2BzbMmbF2dmZbdu2ERISYkqJBnn66aeJjY21u5Z29kZ6ejrDhg1j4sSJWkvJlScpavYioJdHJKVMBV4UQsw1iyqFHkOGDGHIkCF6Y9OnT+fvv/8mPDw82/Gurq7ZxgoqsbGxhIeH57jAxcHB4Yn/v8qXL09YWJhFFso4Ojqqn6MN4ODggIuLC0WLFtVaSq7kecGUEKKPlHKNmfRYBdYaoz9+/Dj//vsv77zzDqGhoURHRxvsQDRjxgzq1atHly5dNFBpW0RHRzN79mx69eplVc24Q0JCWLx4MYMHD8bLy0trOQobwNQLplTenoXZs2cP/v7+LF68mEmTJhEXF0etWrVybDM3ZcoU1q9fb2GVtkl0dDSTJk3ixIkTT3T8wYMHadasWY6praYiODiYSZMmceXKFbPOozAdx44dy3NzekthjKNXD2AtTNGiRSlXrhwjR44kKSkJT09PwsLC2LRpEw8ePMh2/O3btzXtjmRt9O/fnwkTJhjc5+3tTVJSUrbQWE6ULl0aZ2dn4uLiTCkxGy1btiQpKSnH1c8K6+LcuXP4+flZXVmEhxjTStC6Hy/bIf7+/tnu0Ldv386LL77I5cuXs7XAy6l8cUGlbNmylC5d2uA+IUSe/r9q1qxpkfISDg4O6udoQ9StW5dFixbRu3dvraUYRNW6sSGOHj3K559/TlxcHF27duXgwYMGmyB88803bNq0SQOF1sns2bMZN26cwX3R0dF8/vnneU6ZNHc6XXBwMJ9//jlXr1416zwK0zF06FBKlCihtQyDKEefBSFELyFEQExMjNZS9Lh58yYNGzbk448/ZvLkycTHx1O2bFmeeuopgw0xlKM3jCHnHB0dzeTJk/Pk6P/44w/KlStHQkKCKeXpERISwuTJk7l27ZrZ5lCYnu3bt+cYJtQSYxy9eYOTGiKlXCelHGVtddydnZ3x9vZmzJgxSCkpX748MTExrFmzxmAHoqioKH766ScNlFon586do3Llyvz777/Z9lWrVg0pJS+++OIT26tWrRp9+/Y1q6Nv3bo1Uko6duxotjkUpufo0aMsWrSIe/fuaS1FD1WP3gDWml6ZlbNnz1K/fn1WrFjBgAEDtJZj1dy7d4+3336bt99+m+bNm2stR2HHPHjwgMKFC+fa7MZcqHr0dsLhw4eZOHEisbGx+Pj4cOzYMTp37pztuC+//JJ169ZpoNA6KVWqFIsXLzbo5KOiopg4cSLHjx/Ps9379++bQp5Brl+/zsSJE1UJBBvD2dkZR0dHpJRWVa9eOXobwd/fn4EDBzJ9+nQSEhIoWrQoTZo0oVSpUtmO/fnnn9mxY4flRVo5hv7w7t27x/Tp0wkMDMyTreHDh5v120FoaCjTp0+3SLkFhWkJDw+nXr16LF68WGspOoz+fiGE6Ap0BMryyAeGlPLJA56KJ6Jhw4a0adNGF0uWUrJq1Spq1KiRrbjZnTt3tJBo1Xz++ef8+OOPREVF6ZUwqFatGsnJyXm21717dxo2bGhKiXq0atXKKF0K7SlTpgz16tXDw8NDayk6jHL0QogpwETgNHAblVtvdubNm5dt7IUXXuDDDz/EycmJGjVqPFG1xoLK008/TXp6OikpKSbJTx84cGCu+y9evEj16tU1idUqtEUIwfLly7WWoYexnmEUMFxK2UhK2U1K+UzWzZQCFf/jwIEDjBkzhtjYWIQQnDp1isGDB1OnTh02btyoO+7jjz/mn3/+0VCp9dG1a1e++OKLbE4+MjKSMWPGYMzD9/v37xvMroiLi6NZs2b8/fffxsrl6tWrjBkzhsuXLxttQ6EtiYmJVrMOwlhHnw7sN6UQRe706NGDFi1aEBAQoHsIWKdOHZydnSlevLhejY2FCxeyf7/68TxKampqtloksbGxBAQE5Ll2TVJSEiVLluT//u//su1LTk4mPT2d6Ohoo7Xevn2bgIAAg+mzCtugW7duPP/881rLAIyP0f8CvAp8ZEItilzw8/OjS5cuvPvuu7qxjRs34uTkxK1bt/RamakHeNlJS0ujRIkSvPfee0ydOlU37u3tbVTdmiJFivD999/j7++fbZ+7uztRUVGEhYUZrbdFixZmr6ejMC8fffQRDg4OOZbItiRG5dGLDNXrgUpkxOlTsu6XUo4wiTqNsIU8esjIxHF1dWXnzp1MnTrVKlfkWRPTp0+nadOmFikU9tZbb7FkyZJsD38VCnNhjjz6L4BngEJAeaDyI5vCDOzdu5fRo0fzsETDP//8o1sBm7VJxfjx41m5cqUmGq2ZsWPHZnPyd+/eZfTo0Rw+fDjP9pKSkggMDMyWthkfH8/Bgwfp378/6enpRmkNCgpi9OjRXLx40ajzFdZBSEgICxYs0FqG0Y7+bWCElLKelLKTlLJz1s2UAhUZvPrqq7Ru3ZolS5boShNXqlSJYsWK4ebmRtmyZXXHrly50qiHi/ZOamoqoaGhemPx8fEsWbLEqIdmS5YsoW7dugQHB+uNJyUlcf36dZo3b64XUssLERERLFmyhNu3bxt1vsI6+Ouvv3j55Ze5ceOGpjqMDd3cAVpKKe1y2Z41hm7mzp3LnTt3+PTTT3Vj27dvJzIyUtdJytDiKcX/mDp1KpMmTSI+Pp7ixYvn215wcDB79+6lR48e2f7vpZRcv34dIYTqEFWAuXv3LjExMVSvXt3sc+UWujHW0X8JOEop7TIobI2O3hCDBg3i1KlTPHjwgPbt21vFV0Rr5uTJkxw+fJihQ4dSrFgxs84lpcTDw4MBAwYwZ84cs86lUIB5YvTlgdeFEMeFEPOFEAFZN+OlKnJj9+7djBo1Spe7PWvWLFatWsWdO3fw9fXVHffee++xbNkyjVRaL40aNWLUqFF6Tj4iIoJRo0Zx4MABo2xeuXIlW/mEuLg4XnvtNV588UXeeOMNo+xevnyZUaNGcf78eaPOV1gP586d4+2339a0zaCxjr46cBKIAbwA30c2hYn59NNPadu2LevXrycpKQkADw8PXYy+Vq1aumO3bNnC2bNntZJqtUgpuXHjhl6JiPv377N+/fpssfsnZcCAAYwdO1ZvLCUlhfXr11O/fn2jyyRERUWxfv167t69a9T5Cuvhzp07LFiwgHPnzmmmQZUpNoA1hm7WrFnDnj17mD59um5s//79nDp1ioEDB5KSkoKnp6eGCq2f1NRUihYtyvjx4/nqq69MYvO///7D1dWVJk2aZNsXGRnJyZMnadOmDYULFzbJfArbIzU1lbS0NIoUKWLWecwRo18BnJRSTn1k/COgkZTSOpaDGYk1OnpDTJgwgRkzZtCmTRsSEhLUatgnYPHixfm6084L8+fPZ8SIEQQFBVnkYZyiYGOOGH1bYKOB8U2Z+xRmYOfOnQwfPlwXo580aRKnTp3i9OnTtG37v//2t956y6pKpFoTQ4cO1XPy4eHhDB8+nH379hllLyoqii1btuitYo2Li2P48OE4OjqyY8cOypcvn2e7Fy9eZPjw4Xkun6ywTo4cOUKXLl00S7M01tGXBAw9WbgPlDZejiInZs2aRYcOHVizZo2ufK2LiwtFixbF2dmZZs2a6Y7du3evaliRA3fu3OHkyZO694mJiezatcvofPWDBw/StWtXzpw5oxtLSUlh165dpKen0759e6MyfGJiYti1a1e+6uUorAcnJydu376t2boIY0M354CfpZSzHxl/ExgtpaxtIn2aYI2hm71797J27Vo+++wzXQ74yZMn2bp1Ky+88AJJSUkqPPAEvPXWWyxdupSoqCiT2IuOjubs2bM0atRIb3XyQ/bu3YuTk5PeB7FCYQ7MEaN/C/ga+AzYSkY9+q7A58AkKeVMo9VaAdbo6A0xd+5cXn/9dfr168fu3btVhsYTcOrUKcLCwnjmGctU065Tpw41a9ZUZaMVZsfkMXop5Szg/4CvgFNkFDabAvxo607emtm+fTuDBw/WfZ0fPnw4QUFBnD59msGDB+uOGzlyJH/88YdWMq2ahg0b6jn5O3fuMHjwYPbs2WO0zV27dnHw4EHd+7i4OAYPHszmzZtZtmyZwVLGj+PChQsMHjxYpcnaEZs2bcLLy0uT0tNGtySSUn4GeADNM7cyUspPTCXMWIQQfYUQvwoh1gghumQZLy6EOCaE6KmlPmNZtmwZnTp1YunSpaSkZBQLLVKkCI6Ojjg4ONCxY0fdsSdPnlSlinMgISGBffv2ERkZCWTUjj969Gi+vg29/fbbfPPNN7r3KSkpHD16lIiICBo0aGBUCYT4+HiOHj1KbGys0boU1kX58uXx9/c3a1P5nLCqPHohxO9ATyBcSlkvy3g34EcyqmXOk1JOewJbpYHpUspXMt9/ASQA56SU63M71xpDN6dPn2bFihW88847ugJmQUFBLF++nN69e+Pg4ECdOnVUSdzHcOzYMfz8/Pjnn3/o27evSWyePXuW0qVLU7FixWz7rl69yt69e3nhhRdULr3CrJg8Rm8uhBBtyMjmWfjQ0QshCgGXgM5AKHAEeIEMp//1IyZGSCnDM8/7HvhTSnlcCNGJjG8fRYG7tujoDbF582a6detGjx492LBhA4mJiWZflGHrxMfHs3fvXpo2bUqZMmXMPt+8efMYOXIk169fp2rVqmafT2H9JCUlmeXv1Bx59GZBSrkbeDQdohkQJKW8KqVMBpYBfaSUZ6SUPR/ZwkUG3wCbpJTHM220JyO8NBgYKYSwqut+UjZt2kS/fv10MfpOnToRHBzMlStXeOutt3QlcV966SV+++03LaVaLS4uLnTr1k3n5O/cuUP//v3ZtWuX0TavXLnC77//ritNERsbS//+/dm0aRPPPvssly5dMni3nxuBgYH079+f06dPG61LYX3MmDEDd3d33e+KpbAFh1cRyLrKIDRzLCdGA52A54QQrwNIKSdJKd8DlgC/SimzdYMQQowSQhwVQhyNiIgwmXhTsW3bNrp3787q1at1jS4eOnZHR0e6deuGo2NGZ8igoCDCw8M102rtHDx4UJdLn5yczKVLl/IVC9+9ezevvPKK7iFbWloaly5dIjo6Gjc3N3x9fXU/myflwYMHXLp0SZN4rsJ8NGvWjHfffVfXU8JSWFXoBkAI4QWszxK6GQB0lVK+mvl+GNBMSjnaXBqsMXRz/fp1Fi1axLBhw3QP927fvs3s2bPp0qULzs7O1KtXDycnJ22F2gC1a9emfv36rFixwiT27t27R3R0NJUrV87m0NPS0li8eDF16tQx2F9WoTAVuYVujG0ObklC0W9PWAkwvuuyjeLl5cUnn+gnNUVERPDFF19w/Phx1q9fz7Vr11STiydg0aJFlCxZ0mT2SpUqlWPTFwcHB1577TVGjx6tHL0CyPjwj4yM1OsKZ27yFboRQow0lZBcOAL4CiG8hRBOwPPAWgvMa1Wkp6fz999/88wzz+hWddarV4+bN28SEhLCe++9p4s7v/DCC8ydO1dLuVaNn5+frn7/7du36dWrFzt27DDaXmJiIvPmzdOFg2JjY+nVqxfr169HCMGFCxf44osv8mTz3Llz9OrVi1OnThmtS2Gd9OjRg169ell0zvzG6E2aNy+EWAocAGoKIUKFEK9IKVPJ6FG7GTgPrJBSmqWwsxCilxAi4GHzbWvi4cO5f//9l7S0NACEEKSlpeHo6Ejnzp11pRFu376t8q9zITAwkPXrMxKvUlNTCQsLy1csXErJyJEj2bRpE5BxxxYWFkZCQgKQ8W3M2dk5TzaTk5MJCwuz+EM7hfl5/fXXee+99yw652Nj9EKInB77C6CGlNLu8vmsMUYfGRnJvHnz6NOnj67JSExMDNOmTaNDhw64urpSt25dg/VWFPqMGTOGOXPm6ByxKQgJCaFcuXIG0+Z2797NmTNneOutt0w2n0LxKPnKo89sBN4VeLSMngD2SykrmESlFWGNjt4Qt27domrVqvTp04eVK1dy4MABmjdvrrUsqyckJISYmBjq1atnkQVmY8aMYfbs2SqDRgFkfAMMDQ3F0dHRqBLWOZHfPPpNgKuUMviR7Tqw12QqFbmSnp7On3/+Sdu2bXXL98uXL8+1a9cICgpizJgx1KxZE4D+/fsza9YsLeVaNVWqVKF+/foIIbh16xZdunRh69at+bK5Zs0aXQ+AmJgYunTpwpo1awD4/PPP81xi4cyZM3Tp0oUTJ07kS5fC+khJScHb29uif6OPzbqRUg7PZd9Ak6rRGCFEL6CXj4+P1lKycevWLYYOHQpk3BE8REqJs7MznTp1onTpjFYA8fHxupr1iuzcvn2bXbt20alTJ9LT04mPj9etTTCW+fPnc+3aNYYOHYqUkvj4eF1NImPCaWlpacTHx+uexyjsBycnJxYuXEiDBg0sNmee8+iFEH2klGvMpMcqsMbQTUJCAgEBAXTq1In69esDGXcGH330EW3atKF06dLUqVMHDw8PjZVaP9u3b6dTp07s2rVLrzNXfrh37x7Ozs4GY/QhISEsWrSIoUOHqjIICrNh6hIIj9aXUViA4sWL8/777+ucPGTczQcEBLBq1Sratm2br1K7BYmnnnqKwMBAkzYDKVWqVI71S27dusXHH3+s2gIqdDzsIGapb97GOHpVHlED0tPTmTdvHv7+/roYvZOTExcuXCAwMJDx48fTqlUrAHr16sVPP/2kpVyrxsXFhdq1a+Ps7ExYWBjt2rVj8+bN+bJ5+vRpJk+eTFxcHDExMbRr107XbKRp06YkJibmqdnJ6dOnadeuHceOHcuXLoV1smHDBtq3b8+lS5csMp8xjt66aiYUEBISEhg5ciSGQkouLi507NjRItUY7YH09HQWL15sUid67tw5Pv/8c8LCsi/adnR0VFVFFXp07NiRLVu2WGwluzEx+kApZR0z6bEKrDFGn5KSwq+//kqLFi1o1KiRbvyDDz7A398fDw8PateuTaVKlbQTaSNIKSlatCjvv/8+06Y9trXBE5GcnIwQIsea89OmTaNOnTr07t3bJPMpFI9iM2WKtcaaV8YWLlyYN998U8/JA/z999/s3LmTLl26ULlyZQ4dOqSNQBviYVmCiRMnmsymk5NTro1FZs6cyZYtW4CMBuUPF70pCi5Hjhxh3759FpnLGEcfZ3IVVoKUcp2UcpQpC16ZCiklM2fOpE6dOno52Xv37uX48eNUqJCxbi0qKoquXbsa1ae0IOHt7U2JEiW4efMmLVq00JUvMJbY2Fg+/fRTDh06xL1792jRogUrV67U7Q8JCeHnn38GMnLuL168mKu9kydP0qJFC4OhOoV9MGbMGCZMmGCRufJcvVJK+ZQ5hCgezzvvvAOgt5pTCIGbmxv9+vXjv//+o3r16ri4uKhyxY9hzZo1CCHw9/fHxcUlz/XiHyU1NZUpU6bg6elJrVq1cHFx0bvDf9g7ADIawzyuFlGhQoVwcXHRO09hX8yaNUtXn8rcWF09emvAGmP0AL/99huNGzemSZMmurGPPvqI6tWrU6FCBWrWrIk1LvayRtq2bYsQIl+dpbIipURKiYOD4S/JS5Ys4fr167pwUXp6OkII1eNXYTJMHqMXQnQWQrTM8n6kEOKIEGKBEEJV1TITr7zyip6TB9i1axdHjx6lZ8+e+Pr65jtNsKCwcuVKXYkCUyCEyNHJA+zYsYPly5cDsGLFCgoXLkxQUJDJ5lfYHhERESxevJg7d+6YfS5jH8Z+S0azbYQQNYBZwFHAD/jONNIUjzJ9+nSqVKlC1laHf/31F4cPH6ZGjRoA3L9/n/bt2/Pdd+rHkBtlypShZMmShIaG0rRpU13Z4vwwffp0/vjjD6Kjo2natKnOsUNGk/CHteXnzp1Lenp6rqURTpw4QdOmTTl8+HC+dSmsk6CgIIYNG2aR5zDGBiarA2czX/cDtkkp3xBCPA38ZRJlGmDNtW4APvzwQ0A/3luoUCEqVapE69at2bNnD/Xq1cPT05MSJUpoJdMm2LNnD6dOneK5556jQoUKFCtWLN82//rrL3x9fenXrx8VKlTIMf7aunVrGjZsiKenZ462nJycqFChgsq/t2MaNWrE+fPn8fb2NvtcRsXohRCxQCMp5VUhxEZgi5TyByFEFeCilDJvXRasDGuN0S9btgxfX1+aNm2qG/vss88oUaIE1atXx9fXlzp16qi47xMwbtw4Zs6cabEmzQcOHGDBggVMmzaN0qVLk5qaSnp6unporjAZ5sijPw28IYRoA3QAtmSOVwYicjxLkS+ef/55PScPGWl4586do1+/ftSrV08vXKDImU8//VQvBGZubty4werVq7l37x4REREULlyYgIAAi82vsE42bNigK5VhTox19B8BI4CdwEIp5cNqTb3I6PGqMANTp07F3d1dz0HNnDmTI0eO6H0AtGzZkq+/VrXncsPFxQUXFxdu3LhB/fr1Wbs2/22IlyxZwoQJE4iOjqZ+/fosWbJEt2/gwIHcuXMHb29vXe70U0/lnKl87Ngx6tevz8GDB/OtS2G9/Pjjj3zzzTdmn8eoGL2Ucq8QogxQQkp5L8uuXwHT9WdT6DFp0iQAvZzvwoULU6NGDVq2bEnVqlXx8/PDx8fHoh3mbZEzZ86wfv16+vfvT40aNUzyTOPYsWNs3ryZCRMmUKNGDV1/gEfx9vZm3Lhx+Pv752jL2dmZGjVqmOTZgcJ6WbRokUVy6VUevQGsNUa/du1aypcvr+cgPv/8cxwcHHQPYZs1a5bvxT8FgcWLFzNs2DAuX75skbUHoaGhfPbZZ7z++uv4+/uTmJhIamoqLi4uZp9bUTBQtW7shN69e2e7CwwKCuLKlSv079+fli1b8ttvv2mkzrYYOHAgDx48oHr16haZLzk5mc2bN3Pr1i0A6tatyxtvvGGRuRXWy8WLF5k2bRr37t0z6zzK0WfBmouaQcbde+HChQkPD9eNffXVVxw4cEDXKcnZ2Rl/f3++/PJLrWTaBE5OThQtWpTQ0FBq1Khhkgdi//33Hy+//DIhISHUqFFD10MWoFq1aoSGhtK7d29GjhzJ1atXef7553O0dfToUWrUqMH+/fvzrUthvZw/f54JEyZw7do1s86jHH0WrLmoGcDkyZNJTU3Vq6FSpEgR/Pz86N+/Py+88AKtWrWiUaNGVKlSRUOl1s/NmzeZPHkywcHB+Pn5maQF482bN9m+fTsPHjzAz88vx/4APj4+fPzxx/To0SNHWy4uLvj5+an1EHbOM888Q3x8PI0bNzbrPCpGbwBrjdFv376dEiVK6IVvvvzyS+Li4mjevDnOzs506NBBLbJ5Ak6cOEGTJk1Ys2aNxWrEjxgxgg4dOjB06FASEhK4f/++ahajMBkqRm8ndOzYMVuM/tatW4SGhjJ69Gi6d+/ODz/8oI04G6Nhw4akpqZatBHI8ePHuXHjBgAjR46kZcuWjzlDURD47rvvTFKCIzeMLWp2WQhxycB2UQhxSgjxlxCii6nFFnQmTZqEEEIvRj9u3Dj279+v++rn7u5Ow4YN+fzzzzVSaRs4ODhQqFAhQkJC8PLyYtWqVfm2eeXKFYYOHcqOHTvw8vLijz/+0Nt/8uRJJkyYwPDhw1m9ejWffvppjrYOHz6Ml5eXxRpTKLTjhx9+YOPGjWadw9g8vL+AN4HzwIHMseZAbWAB0ADYJIToK6Vcl1+RigymTp0KoLds3tnZmXbt2vHUU0/h4eFBx44dOXHihCpX/BiSkpL44osvaNSoEe3atcu17kxebB48eJDnnnuOdu3aUbFiRYPHNWjQAG9vb4YOHZqjrZIlS9KuXbscc/EV9sPVq1fNHm41ttbNLCBGSjnxkfEpQCkp5dtCiG+AdrbYqMRaY/T79++ncOHCeuGbqVOnEhoaSseOHUlOTqZ3794Wa2ZgyyQnJ1O8eHG++uorxo0bZ5E5P/vsM9LS0pgyZQrx8fFERkZSuXLlXMsbKxRPijli9M8D8w2M/wEMzny9CFCNMU1IixYtssXoY2JiiIqKYsqUKQwePFiVPnhCnJycSElJsZiTh4x2gjdv3gQyShV7eXkRHx9vsfkV1snatWv55JNPzDqHsY6+EFDDwHiNLDaTgHQj7SsM8OGHHyKE0GtU8Prrr7N7926qVasGQKVKlahTpw4ff/yxVjJtiuDgYCpUqMBff+W/uvaDBw8YOHAg8+fPp0KFCvz+++96++fPn8/8+fMZOnQoP/74I7///nuO1SsPHTpEhQoV2LNnT751Kaybffv28euvv5p1DmNj9CuAeUKICcAhQAJPA1OAZZnHPA1cyLdCC2Lt9einT58OoBfPK168OD179sTPzw93d3e6devGhQsXqFevnlYybYbPP/8cT09PevbsSaVKlfJtz8HBgTNnztCyZUt69uyZY53xZs2aUadOHV5++eUcbbm5udGzZ0+T5PcrrJtp06aZvbCZsTH6omR0mRoFFAYEkAzMBcZJKZOEEPUAKaU8Z0K9FsFaY/THjh1DCKHXTvCbb77h7Nmz9O7dm+joaAYNGoS1LviyNqpXr06vXr0slpL6888/c+jQIRYtWsT9+/e5efMmlStXpmjRohaZX2HfmDxGL6VMlFK+A7gBjYFGgJuU8l0pZVLmMWdt0clbM02bNs3WMzY5OZkHDx7w/fff89prr6nSB3ngypUrFl13EBsbqysxvXPnTmrUqMGZM2csNr/COjl//jxvvfWWWcsg5NnRCyEKCyHuCCHqSinvSylPZ273zSFQ8T/effddhBDcvn1bNzZkyBB27typCz34+PhQvXp1PvroI61k2hTXr1/H3d2dZcuWPf7gJ2D48OFMnjwZd3d35s2bp7dv4sSJ/PvvvwwaNIgPP/yQRYsW4eXlZdDOgQMHcHd357///jOJLoX1cvfuXZYtW6YreGcO8hyjl1KmCCHSgDQz6FHkwty5c4GM3PmHuLq6MnjwYJo0aUL58uXp3r07165dw8/P4Dc4RRa+++47EhMTGTx4sO5hdn65ffs2np6eDB48GF9fX4PHtGvXDn9//1zz6MuUKcPgwYNNkt+vsG5at25NZGSkWecwNkb/BeAhpXzT9JK0x1pj9KdPn0ZKScOGDXVjM2bMYNeuXQwdOpQbN24wfPhw3N3dNVRpO3Tu3Nmkd/OP46+//mL27Nls3LgRBwcHgoKCKF++vFoUpTAJucXojc26qQAMEEJ0AI7xSFcpKeUoI+0qcqFBgwbZxhwcHHB0dGTWrFns3r2b69evM3PmTA3U2R5bt2616HxpaWmkpKSQkpJCeHg4devWZdGiRbne2Svsn/T0dN544w26du3Ks88+a5Y5jM2jrw4cB26R4fR9s2zWmZtoB7zxxhvZYvS9evVi69atutaBderUoUqVKnz44YdaybQprl27hqurq15/1/wwZswYRowYgaurqy7U9pDnn3+ePXv2MHz4cF544QWWL19O69atDdrZv38/rq6u7Nq1yyS6FNaLg4MDW7du5dKlS2abw9iese1NLUTxeJYuXQqg10e0ZMmSjBo1SleDvmfPnty8eVOvWbjCMHPmzOHs2bOMGjWKGjUMrf/LOw8ePCA9PZ1Ro0ZRt25dg8d069aNhIQEBg4cmKMdT09PRo0alWO9HIV9cfXqVbPaz3dzUSGEOxAlVWF7s3Po0CFSU1P1mlEsWbKE48eP4+/vj4uLC4ULF2bKlCkaqrQdjh8/zqVLl0x61/zLL7/kuG/v3r189NFH/P777/j6+nLy5Ek8PDwMLtaqVq0a33//vcl0KQo2xpYpLiSEmCyEiAbuAN6Z49OEEK+ZUqDif9SsWTPbXaKzszMlS5bk999/54svvsi19K1Cn4CAAIuGRhwdHXFyciI9PaMyiJ+fH7Nnz7bY/ArrZcaMGYwfP95s9o2N0Y8HXgLeIWNF7ENOAMPzqUmRA6+++ipCCL182/bt27Nx40ZKlSoFQOPGjSlXrhzvv/++Ripti6tXr+Lk5KTX3zU/fPvtt3To0AEnJ6dsd/fNmzdnx44djBs3jpYtW/L3338zbNgwg3b27t2Lk5MTO3bsMIkuhXVz5coVzp8/bzb7xoZuXgJel1L+K4TI+tt8BsPFzmwCa691s2HDBgC9MsSlS5dm7NixNGjQAB8fH3r16sXdu3dp1KiRRipth6VLl7Jq1SrGjh1LnTp1TGJTCIGzszNjx47NsQ/os88+y/3793PtblWpUiXGjh2rev8WEGbNmmXeCaSUed6AB0DVzNdxQLXM1zWA+8bYtKatadOm0hq5cuWKvHDhgt7YL7/8Ihs1aiSXLFkiP/jgAxkSEqKROtvjq6++knXr1rXYfFevXpXNmjWTmzdvllJKeeLEiWw/T4XCWICjMgefZmzo5jrQ0MB4Z2ysYqUtUa1aNWrWrKk3VrJkSapWrcqSJUuYMWMGVapU4f59VY3iSZg4cSJnz5612HyOjo64ubnpShMPGjSIQYMG8eabbxIUFGQxHQrrY926dfTu3ZsHDx6Yxb6xjv4X4McsfWF9hRBvAl8BarWOmXjxxRezxeibNWvGmjVr9MoiFC9enHfeeUcLiTbH1atXEUKwcOFCk9gLCAigRo0aCCH4+eef9fZVrlyZTZs28e2339KsWTN+/fVXPD09mT17Nv/++6/esXv27EEIwfbt202iS2HdxMbGEhISQlJSklnsG5tHPzMzrfIfwBnYBCQCU6WUhjpPKUzAwyYULi4uujE3Nzc+++wz6tWrR0xMDFWrVqVSpUrZOlEpsvPPP/8we/ZsJk2apFdWIj8ULVoUT09PXnjhBZo1a2bwmCFDhnD//n3atGlDuXLlWLp0abYVkVWqVOGzzz7Lsaa9wr4YMmQIQ4YMMd8EOcV0nmQjw8n7Ac2A4vmxZU2btcbob9y4Ia9cuaI3Nnv2bFm7dm25aNEiOXLkyGz7FTkzZ84cWatWLZmYmGiR+aKjo2XDhg3ln3/+KaXMiNEfP37cInMr7B9MEaMXQmRbaimlfCClPCqlPCylTMhybFEhRO38fwwpslK+fHkqVar08EMWAA8PD3x9fVm1ahW//vor1atX5+7du6Smpmqo1DZ47bXXCAwMRAhBWprpirFKKUlOTs5ms3DhwlStWpVixYqRkpLCO++8w3PPPcewYcM4cuSI3rHp6ekkJyfrcu4V9s3x48d55plnzJZimZcY/RohxD9CiK5CCIPnCSEqCiE+Ai4DLU2iUKFj+PDhFClSRK/WTaNGjVi7dq1ee8EyZcrwwQcfaCHR5rh27RpFihThzz//NIm9pUuXUrlyZYoUKZJtMVTx4sVZs2YNAQEBtGzZkh9++IGWLVuyePFiNm7cqHfsvn37KFKkCDt37jSJLoV1k5aWRmRkJImJiWaxn5cYfU3gI2AxUFQIcQK4SUZs3g2oS8YK2V3AC1LKvaaVqnhYOtnV1VU35ubmxldffUXt2rV58OABVapUoUKFCqrWzROwefNmvv76az788MMcc97ziqurK7Vr12bIkCE0b97c4DEvvfQS9+/fp0mTJnzyySfUqlUrW//YqlWr8tVXX1G9enWT6FJYN/7+/hw+fNh8E+QU08lpA4oAfYEZZDyM/RdYBLwP1MqrPWvcrDVGf/v2bRkcHKw3NnfuXFm9enX5+++/y+eff16eP39eI3W2x4oVK2SDBg1keHi4xeasV6+e/OGHH6SUUp48eVLu2bPHYnMr7BtMmUcvpUySUq6WUn4gpewnpewmpRwmpfw/KaXKoTcjpUuXpnTp0npx23LlylGnTh02bNjAsmXLqF27NsHBwWb7CmhPDBgwgBMnTlC0aFFSUlJMZjc9PZ24uDiSk5Oz7atXrx6urq4kJCQwZcoU+vfvT9++fdm8ebPecampqcTFxalnLQWEiIgIOnXqxLp168xi39g8eoUGvPzyy5QoUYI7d+7oxurWrcu6detwdPxfFM7Ly4tx48ZpIdHmuH79OiVKlNCVgM4vmzZtonz58pQoUYKAgIBs+5cuXcrKlStp3749X375JUOGDGHNmjVs2rRJ77gDBw5QokQJ1TO2gODo6EhiYqJJkwL07JvFqsIsPHwinzVG7+7uzvTp06lZsyZCCCpVqkT58uVVrZsnYPfu3Xz00Ue88847JnumUaJECfz9/fHx8aFlS8P5CK+++ioPHjygVq1avPXWW1SsWDFblylvb2+mT5+OtdZdUpiW0qVLs3evGR9r5hTTKcibtcboIyMjZVhYmN5YQECArFSpkpw7d67s3r27PHnypEbqbI+tW7dKf39/GRQUZLE5mzdvLidMmCCllPLMmTNy06ZNFptbYd9ghlo3Cg0oVqyYXj1zyKhy2KRJE3bu3MnGjRtp1KgR58+fJyEhIRdLCoBOnTpx4MABSpUqZdKl5+np6URGRhqsW+Lv70/ZsmWJiYlhzpw5DBgwgE6dOmVL70xJSSEyMtKkzw4U1k337t2ZMWOGWWwrR29DvPzyy3h4eOjF6H19fVm7di0ZH+gZ1KlThwkTJmgh0eYIDg7Gw8OD5cuXm8Te/v37cXNzw8PDg99++y3b/p9++oktW7bQuXNnPvzwQ8aPH8/27dvZsmWL3nEHDx7Ew8OD3bt3m0SXwvopWrQohQsXNottk8XohRAeQFPgtpTylKnsWhJrr0d//fp1AL1Wgh4eHvz000/4+vpSrFgxKlSogKenJ/Xr19dIpe1w7NgxXn/9dUaNGmWy2kAlS5akQ4cOeHp60qZNG4PHvP766zx48ICqVasyZMgQSpYsyaBBg/SOqV69Oj/99JPJetkqrJ+///7bfMZziunkdQP2A7OA7cAWoKypbFt6s9YYfUxMjIyIiNAbmzdvnixbtqz8+eefZevWreXhw4c1Umd7HDlyRLZs2dKizzW6du0qX375ZSmllIGBgXL58uUyPT3dYvMr7BcsEaOXUraQUr4lpexIRhnjzUKIyqayr4BChQplq6FStWpVnnrqKQ4cOMCePXto1qwZR44cISYmRkOltoGfnx+7du2iTJkyJq3hn5aWRlhYmMHnJC1atKBq1apERESwcuVKBg0aRKdOnfjhhx/0jktKSiIsLMxsZWsV1sdrr73Ga6+Zp+W2SWP0QojCQoiqZDQM3wSohpcm5JVXXqFixYqEh4frxry8vFi3bp3ewppmzZrxySefaCHR5ggJCaFixYqsXLnSJPbOnz+Pi4sLFStWZP787BW7P/30Uw4dOkSPHj0YNWoUX3/9NTt27MhW0+bw4cNUrFjRvCl3CquidOnSuLm5mcW2KWP0EZn2bmVuYcAqU9lXoHPwJUuW1I2VKVOGuXPnUr16dTw8PChfvjxly5Y1WQ9Ue+bixYsMGjSIF198Mce6NHnF1dWVfv364eLiQrt27QweM3r0aBITEylXrhwvvPACbm5u9OvXT+8YHx8f5s6dS61atUyiS2H9TJs2zXzGc4rp5HUDipnKltabtcboExIS5L179/TGfv/9d1myZEn5f//3f7Jx48aqdkoeuHTpkmzfvr3cu3evxeYcOHCg7Nq1q5RSysuXL8tFixbJ+/fvW2x+hf2CqWL0QohNQogpQoh+j8bfpZSqUamZSUtLIzo6Wi9MU61aNdq2bcuJEyc4ceIErVu3ZteuXURGRmqo1Dbw9fVl8+bNVKxYkfj4eJPZTUtL4/r168TGxmbb17p1a+rXr09YWBg7d+5k2LBhPP3000ycOFHvuMTERK5fv262HqIK6+Orr76idevWZrGd1xh9CjCcjJDMdSFEeKbz/1II0Vc9fDUvr7zyCt7e3kREROjGKlasyNq1a/Ue2rVv357JkydrIdHmuHHjBt7e3iZLbbt9+zZFihTB29vbYB/at99+m3PnztG3b18GDBjA999/z6lTpzh48KDecUeOHMHb25v9+/ebRJfC+ilXrpzZSl7kKUYvpewNIIQoR0bO/MNtODAJkEKIcDJKF8+QUp4xqdoCzsO7zqwx+rJlyzJ//ny8vb2pXLky5cqVo0yZMiq2+wSEhYXRs2dP+vfvT4sWLUxis1ixYgwdOhQHBwc6duxo8Jj333+fpKQkSpUqxYABA3Bzc6Nnz556x9SoUYP58+dTu7Zq1FZQePXVV3n11VfNYzynmE5eN6As0B34DDhIRkOS1qayb8nNWmP0SUlJ2eK5CxYskM7OzvLbb7+V1atXl1u3btVIne1x+/Zt2bVrV7llyxaLzTlq1ChZv359KaWUwcHB8vfff5d379612PwK+wUL5dGHSyk3SiknSymbA1OBr01lXwH379/nxo0bejF6Hx8fOnfuzIULF7hy5QqdO3dmw4YNeimYCsOUK1eOtWvXUrVqVYPxdGNJS0vj0qVLREdHZ9vXpk0bWrZsSUhICCdPnmTEiBE0adKEV155Re+4+/fvc+nSJZPm9yusmz///JNq1aoRFRVlctvmrHWzEGhoRvsFjpEjR1KzZk29GP1DZ5V1cU7Pnj2ZMmWKFhJtjtDQUGrWrMnq1atNYi81NRVHR0dq1qxpsA/tkCFDCA4O5rnnnqNjx47MnDmTkJAQTp8+rXfc0aNHqVmzZrbYvcJ+KVeunMlCiI9iznr0IUAvM9ovcDxcEVuqVCndWLly5fjzzz+pUqUKvr6+lC1bFg8PD3x9fTVSaTvExcXRtWtXunXrRqtWrUxis1ChQowaNYoHDx7QpUsXg8eMGzeOxMREihcvTt++fXFzc6Nbt256x9SqVYs///xTrYcoQHTq1IlOnTqZx3hOMZ2CvFlrjD4tLU2mpqbqjS1cuFAWKlRIfvXVV7JcuXISkP369dNIoW2RkJAge/fuLVevXm2xOceMGSPLli0rpZTy5s2bcs6cOfLmzZty7969sk+fPnLMmDHy9u3bFtOjsB9Q9ejtg+joaAIDA/VqlNeoUYMePXpw9epVXfnif/75h7CwMK1k2gzFihXjr7/+olq1agbj6caSmprKmTNnDK5laNu2Lc888wxBQUFcunSJ119/ndq1a9OqVSvWrFnD999/z40bN0hISODMmTMmze9XWDebNm2iatWquk5ypkQ5ehti5MiRNGjQgLt37+rG3N3dWbt2bTaHYNbl1HZEaGgoDRo0MGlT5sKFC9OgQQODfWh79epFeHg4gwcP5umnn2bu3LnZHgQXLVqUY8eO0aBBAw4fPmwyXQrrpmzZsnTo0AEXFxeT21Y9Y20IJycnIHuMfsWKFVSpUoVGjRoxbdo0GjduzIgRIzRSaTvEx8fTrl07evbsSdu2bU1m99133yUhISFb3P0hEydOJDk5mSJFitCrVy/c3d2pXLkyBw4coFy5crr1ECtXrqRevXom06Wwbpo2bWqwEJ5JyCmmU5A3a43RG+KPP/6QgJw8ebJ0dXWVgOzWrZvWsmyChIQEOWDAALlu3TqLzfnee+/JUqVKSSmlvHHjhvzxxx/ljRs35H///Se7dOki3377bXnr1i2L6VHYD6gYvX0QERHBsWPH9GL0tWrVonfv3oSEhBAXFwfAv//+S0hIiFYybYZixYqxaNEiypcvb9LaQKmpqRw7dszgWobOnTvTt29fzp8/T1BQEO+++y5eXl60bduWLVu28PPPP3Pz5k3i4uI4duyYSfP7FdbNhg0b8PT0JDAw0OS2laO3IV599VX8/Pz0YvRubm4GY/TTp0+3tDyb5ObNm/j5+bFhwwaT2SxcuDB+fn6sWLEi277u3btz584dXnrpJVq2bMmCBQv0GskAFClShBMnTuDn58fRo0dNpkth3VSoUIG+ffvqtQo1FSpGb0M8rHFTunRp3VjZsmVZsWIF3t7eNG/enB9//JH69eszatQorWTaDHFxcbRs2ZK+ffvSvn17k9kdP3489+/fp0ePHtn2paWlMWnSJNLT0ylcuDBdu3Zl7dq1VKxYkSNHjlCmTBmqVKlC+fLlWbt2LQ0aNDCZLoV107hxY+bMmWMe4znFdAryZosx+k8//VQWKlRIArJdu3Zay7IJ7t+/L4cMGSI3bNhgsTnfffddWbJkSSmllCEhIfK7776TISEhcseOHbJNmzZy1KhRMiwszGJ6FPYDKkZvH9y+fZv9+/eTnJysG6tduza9evUiNDRUFwLYtWsXV69e1UqmzeDs7My8efMoVaqUSWsDJSUlsXfvXm7dupVtX/fu3Xn22Wc5ffo0V65c4cMPP6RKlSp06NCB3bt3ExAQwM2bN4mJiWH//v2q928BYt26dbi5uXHu3DmT2xYZHwSKrGS2RQx+wsM9gLuPPcp2sefrs+drA3V9towx11ZVSlnG0A7l6POJEOKolNJPax3mwp6vz56vDdT12TKmvjYVulEoFAo7Rzl6hUKhsHOUo88/AVoLMDP2fH32fG2grs+WMem1qRi9QqFQ2Dnqjl6hUCjsHOXoFQqFws5Rjt5IhBDdhBAXhRBBQoiPtNaTX4QQlYUQO4UQ54UQ54QQ72aOuwkhtgohLmf+W/pxtqwVIUQhIcQJIcT6zPf2dG2lhBArhRAXMn+GT9vZ9b2f+Xt5VgixVAhR1JavTwjxuxAiXAhxNstYjtcjhJiQ6WsuCiG65nU+5eiNQAhRCJgFPAPUAV4QQth6c89UYIyUsjbQHHgr85o+ArZLKX2B7ZnvbZV3gazte+zp2n4E/pVS1gIaknGddnF9QoiKwDuAn5SyHlAIeB7bvr4FwKMNCwxeT+bf4fNA3cxzfsn0QU+McvTG0QwIklJelVImA8uAPhpryhdSyltSyuOZr+PIcBQVybiuPzIP+wPoq4nAfCKEqAT0AOZlGbaXaysBtAF+A5BSJksp72En15eJI+AshHAEigFh2PD1SSl3A1GPDOd0PX2AZVLKJCnlNSCIDB/0xChHbxwVgRtZ3odmjtkFQggvoDFwCCgnpbwFGR8GQFkNpeWHH4BxQHqWMXu5tmpABDA/MzQ1TwhRHDu5PinlTWA6EALcAmKklFuwk+vLQk7Xk29/oxy9cQgDY3aRpyqEcAFWAe9JKe2i64UQoicQLqU8prUWM+EINAFmSykbAwnYVhgjVzJj1X0Ab6ACUFwIMVRbVRYl3/5GOXrjCAUqZ3lfiYyvkjaNEKIwGU7+Tynl35nDd4QQ5TP3lwdMV+bRcrQEegshrpMRZusghFiMfVwbZPw+hkopD2W+X0mG47eX6+sEXJNSRkgpU4C/gRbYz/U9JKfrybe/UY7eOI4AvkIIbyGEExkPStZqrClfCCEEGTHe81LKGVl2rQVeynz9ErDG0tryi5RygpSykpTSi4yf1Q4p5VDs4NoApJS3gRtCiJqZQx2BQOzk+sgI2TQXQhTL/D3tSMYzJHu5vofkdD1rgeeFEEWEEN6AL3A4T5ZzKlSvttw3oDtwCbgCTNJajwmupxUZXwdPAyczt+6AOxkZAJcz/3XTWms+r7MdsD7ztd1cG9AIOJr581sNlLaz65sMXADOAouAIrZ8fcBSMp43pJBxx/5KbtcDTMr0NReBZ/I6nyqBoFAoFHaOCt0oFAqFnaMcvUKhUNg5ytErFAqFnaMcvUKhUNg5ytErFAqFnaMcvUKhUNg5ytErFAqFnaMcvULxCEKID4UQZ0xs00UIcVMI4W9Ku1nszxVCTDeHbYXtoxy9QpGdJmSsDDYl44GjUsojT3qCEGKNEMLgUnchhJMQ4q4QYkrm0BfAG0KIaibQqrAzlKNXKLJjUkcvhCgKvAHMzeOpcwF/IURDA/ueJaPMwTzQlfLdDryZD6kKO0U5eoUiC5llmn2AKCHECiFErBDijhDitSzH/CWE2JzlfRkhxH0hRNMczHYDnIEtBuYbndn+LzGzhdykzOYaAP8CwcBIAzZHAlullNezjP0DFKTyvYonRDl6hUKfxmT8XbwDLCGjWNhC4OfMZh7wSNlYKWUEGU1a+uZgsy1wQkqZmnVQCPE5MBaYANQmo9Xha8BnmXbTybhjHyqEcM5yXnWgPRDwyDyHgHJCiNp5uF5FAUA5eoVCn8ZkVBQcJKVcLaW8SoajdyQjVALZ64MDxJNzhyNv4GbWASFEMTI6Xr0mpfxHSnlNSrkR+BgYneXQ34DiwIAsY68Cd8heGjs0818Vp1fooRy9QqFPEzIaNF/KMuYL3Od/zR5CARchREmAzDv9VsB/Odh0BhIfGaubOb5KCBH/cCMjLl9SCFEGdC3l1pMZvskM6wwH5j/6DSHLHM4oFFlwfPwhCkWBogmwwsDY6cxQCvzvzrkyEAN8Ddwmo/ORISIAt0fGHt5kDSCjr8GjZG0cPRfYlBmSqQ2UQ7/J+UMezhGRgw5FAUU5eoUik8zsmNrA8Ud2NXlk7KGjr5SZFz8MaCWlfPSu/SHHgbcfGTtHxh14tcyQTW5sAa6TcVdfG9iWGVJ6lPpAGnDiMfYUBQzl6BWK/9GAjL+JR5uINyGjl+5DbgLpZDw8bQZ0kVKey8XuJuB7IURlKeUNACllvBBiKjA1ozseWzPnrg80llKOf3iylDJdCPErGTF9V2BQDvO0A/ZKO2nqrjAdKkavUPyPxkCYlPLOwwEhREUyQiW6O/rM2PgdoCHQ8XGLoKSU54FdZNz5Zx3/EnifjIerp4C9me+vGzDzOxkPZe9ioDdqZi/VweQ9V19RAFCtBBUKCyCEaA0sA3yllPfNYH8g8AnQSEqZZmr7CttG3dErFBZASrmHjAbX3maaogjwsnLyCkOoO3qFQqGwc9QdvUKhUNg5ytErFAqFnaMcvUKhUNg5ytErFAqFnaMcvUKhUNg5ytErFAqFnfP/SmXpq1tWAPkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Grab the GalaxyPopulation instance\n", "pop = mgb.pops[0] \n", " \n", "\n", "# Compute cosmologically-limited solution\n", "zi, zf = 40., 10.\n", "e_nu = np.array([pop.Emissivity(zf, energy) for energy in E])\n", "e_nu *= (1. + zf)**(4.5 - (alpha + beta)) / 4. / np.pi \\\n", " / pop.cosm.HubbleParameter(zf) / (alpha + beta - 1.5)\n", "e_nu *= ((1. + zi)**(alpha + beta - 1.5) - (1. + zf)**(alpha + beta - 1.5))\n", "e_nu *= c * ev_per_hz\n", " \n", "# Plot numerical solution again\n", "pl.semilogy(E, flux[0] * E * erg_per_ev, color='k', ls=':')\n", " \n", "# Plot analytic solution\n", "pl.semilogy(E, e_nu, color='k', ls='-')\n", "pl.xlabel(ares.util.labels['E'])\n", "pl.ylabel(ares.util.labels['flux_E'])\n", " \n", "\n", "pl.savefig('ares_crte_uv.png')" ] }, { "cell_type": "markdown", "id": "957779ef", "metadata": {}, "source": [ "**NOTE:** In reality, the ionizing background before reionization should be heavily damped. This example is unphysical in some sense because while it treats the opacity of HI and HeI Lyman lines (which produce the sawtooth modulation) it ignores the continuum opacity at energies above 13.6 eV. This will be treated more carefully by setting ``pop_approx_tau='neutral'`` in [Example: X-ray background](example_crb_xr)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.7.6" } }, "nbformat": 4, "nbformat_minor": 5 }