docs/examples/example_crb_xr.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "88e92827",
   "metadata": {},
   "source": [
    "# The Metagalactic X-ray Background"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7174c853",
   "metadata": {},
   "source": [
    "In this example, we'll compute the Meta-Galactic X-ray background over a series of redshifts ($10 \\leq z \\leq 40$). To start, the usual imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5eb78166",
   "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 c, ev_per_hz, erg_per_ev"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "83dae27c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize radiation background\n",
    "pars = \\\n",
    "{\n",
    " # Source properties\n",
    " 'pop_sfr_model': 'sfrd-func',\n",
    " 'pop_sfrd': lambda z: 0.1,\n",
    " \n",
    " 'pop_sed': 'pl',\n",
    " 'pop_alpha': -1.5,\n",
    " 'pop_Emin': 2e2,\n",
    " 'pop_Emax': 3e4,\n",
    " 'pop_EminNorm': 5e2,\n",
    " 'pop_EmaxNorm': 8e3,\n",
    " 'pop_rad_yield': 2.6e39,\n",
    " 'pop_rad_yield_units': 'erg/s/sfr',\n",
    " \n",
    " # Solution method\n",
    " 'pop_solve_rte': True,\n",
    " 'tau_redshift_bins': 400,\n",
    "\n",
    "\n",
    " 'initial_redshift': 60.,\n",
    " 'final_redshift': 5.,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f89d7d2a",
   "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 power-law spectrum with index $\\alpha=-1.5$, given by ``pop_sed`` and ``pop_alpha``, extending from 0.2 keV to 30 keV.\n",
    "* A yield of $2.6 \\times 10^{39} \\ \\mathrm{erg} \\ \\mathrm{s}^{-1} \\ (M_{\\odot} \\ \\mathrm{yr})^{-1}$ in the $0.5 \\leq h\\nu / \\mathrm{keV} \\leq  8$ band, set by ``pop_EminNorm``, ``pop_EmaxNorm``, ``pop_yield``, and ``pop_yield_units``. This is the $L_X-\\mathrm{SFR}$ relation found by [Mineo et al. (2012)](http://adsabs.harvard.edu/abs/2012MNRAS.419.2095M).\n",
    "\n",
    "See the complete listing of parameters relevant to :class:`ares.populations.GalaxyPopulation` objects [here](../params_populations.html).\n",
    "    \n",
    "Now, to initialize a calculation and run it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "66b26a14",
   "metadata": {},
   "outputs": [],
   "source": [
    "mgb = ares.simulations.MetaGalacticBackground(**pars)\n",
    "mgb.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dc95ce5",
   "metadata": {},
   "source": [
    "We'll pull out the evolution of the background just as we did in the [UV background](example_crb_uv) example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "775cca21",
   "metadata": {},
   "outputs": [],
   "source": [
    "z, E, flux = mgb.get_history(flatten=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08e9da0c",
   "metadata": {},
   "source": [
    "and plot up the result (at the final redshift):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f45f6a27",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$J_{\\\\nu} \\\\ (\\\\mathrm{erg} \\\\ \\\\mathrm{s}^{-1} \\\\ \\\\mathrm{cm}^{-2} \\\\ \\\\mathrm{Hz}^{-1} \\\\ \\\\mathrm{sr}^{-1})$')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEMCAYAAADHxQ0LAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAv1klEQVR4nO3deXQUVfr/8fdDIAjCgMhyRJYAYR+UfScy7CCL8EVB/aqgguI283P5jjPCoDM44jbOODqijIIbqywCgqjsmyNhEQyb7AQFIiKDCwrk+f1Rldi0SUg61V3d6ed1Th86t7uqPulzePrm1q1boqoYY4wp+or5HcAYY0xkWME3xpg4YQXfGGPihBV8Y4yJE1bwjTEmThT3O0A0qlixoiYlJfkdwxhj8m3Dhg1fqWqlvN5jBT8HSUlJpKam+h3DGGPyTUQOXOg9NqQTQET6icgrJ0+e9DuKMcZ4zgp+AFWdr6ojy5Ur53cUY4zxnBV8Y4yJE1bwjTEmTljBN8aYOGEF3xhj4oQVfGOMiRNW8AOEOi3z+++/57HHHmP16tVhSmaMMYVnBT9AqNMyExISePTRR1m+fHl4ghljjAes4HugZMmSlCtXjoyMDL+jGGNMrqzge6RSpUocO3bM7xjGGJMrK/geqVy5shV8Y0xUC3nxNBGpD9QGSgEZwCZV/darYLGmcuXK7N692+8YxhiTqwIVfBFJAu4C/heoAkjAy2dFZDUwAZipcXZ39MqVK7Nu3Tq/YxhjTK7yPaQjIk8DnwH1gT8CvwbKASWBy4A+wFpgPLBZRJp7njaKVapUiYyMDDIzM/2OYowxOSpID78sUE9Vv8jhtaPu4yNgtIhcCzQENhY+YsGIyDXA1UBl4EVV/cBtvxhYCYxV1QVeH7dy5cpkZmby9ddfU7FiRa93b4wxhZbvHr6q3plLsc/pvTNV9e2ChhGR10TkmIh8FtTeS0R2ishuEXn4Aseeq6ojgGHAkICXfg/MKGim/KpcuTIAR48eDdchjDGmUKJtls5koFdgg4gkAC8CvYFGwPUi0khEmojIgqBH5YBNR7vbISLdgG04f4WERbVq1QBIT08P1yGMMaZQCnWLQxEZoaoTvQqjqivdE8OBWgO7VXWve8xpwABVfQLom0MmwTmPsEhVs4aUfgNcjPOF8YOILFRVTwfba9SoAcChQ4e83K0xxnimsPe0HQN4VvBzcTkQWEXTgTZ5vP9eoBtQTkSSVXWCqj4CICLDgK9yKvYiMhIYCT8X74KoWrUqxYoV4+DBgwXe1hhjIuGCBV9EtuT2Es7UzHCTHNpynfKpqs8Dz+fy2uQ8tntFRL4E+iUmJrYoaMjixYtTtWpV6+EbY6JWfnr4VYCewImgdsGZhhlu6UD1gJ+rAfk6eVxQqjofmN+yZcsRoWxfvXp16+EbY6JWfgr+IqCsqm4OfsG90Crc1gN1RaQWcBgYCtwQgeMWWI0aNdi4MeIzUY0xJl8uOEtHVYep6qpcXrvOyzAiMhVYB9QXkXQRuU1VzwL3AIuB7cAMVU3z8rgBxw9pPfwsWT18u/jKGBONCjwtU0QGhCMIgKper6qXqWoJVa2mqq+67QtVtZ6q1lHVx8N4/JDWw8+SnJzMjz/+yBdfhGXEyRhjCiWUefhPeJ4iShS2h1+vXj0Adu7c6WUsY4zxRCgFP6dZM0VCYXv4WQV/165dXsYyxhhPhFLwi+wqmIXt4VetWpXSpUtbwTfGRKVoW1rBV4Xt4YsI9erVs4JvjIlKVvA9Vr9+fRvDN8ZEpVAK/inPUxQhjRo1Yu/evXz33Xd+RzHGmPMUuOCral7r2MS0wo7hAzRt2hRVZevWrR4mM8aYwrMhnQCFHcMHp+ADbN682ZtQxhjjkVAuvCohIkdFpHE4AsW66tWrc8kll7Bp0ya/oxhjzHlCGdI5A5xzH0WKF0M6IkLTpk2th2+MiTqhDun8G7jPyyDRwIshHXCGdbZs2cKZM2c8SmaMMYUXasGvCtwoIjtE5G0ReSXw4WXAWNS2bVtOnz5tvXxjTFQJ9Y5XdYCsdYCrBr1WZK/Eza8OHToAsGbNGlq1auVzGmOMcYRU8FX1N14HKUouv/xykpKSWLNmDb/73e/8jmOMMUCIQzoi0l1EOgT8PEJE1ovIZBEp6128yPLipG2WDh06sHr1alTj/g8eY0yUCHUM/ymgIoCI1ANeBFKBlsDT3kSLPK9O2gJ07NiRI0eOsGfPHg+SGWNM4YVa8OsAn7nPBwIfqeooYATQ14tgsa5Lly4AfPjhhz4nMcYYR2GutM0aq7gK+MB9fhi4tFCJioi6detSq1Yt3n//fb+jGGMMEHrB3wKMEpEUoAs/F/zqQIYXwWKdiNCzZ0+WLFnCTz/95HccY4wJueA/DNwKLAPeUNVtbns/YL0XwYqCXr168d1337FmzRq/oxhjTGgFX1VXA5WAS1V1ZMBLE4G7vQhWFHTp0oXExETmz5/vdxRjjAl9DF9VM1X1m6C2Pap6pNCpfOLltEyAsmXL0qNHD9555x0yMzM92acxxoTKlkcO4OW0zCzXXnsthw4d4pNPPvFsn8YYEwor+GHWv39/EhMTmTlzpt9RjDFxzgp+mJUvX54ePXowY8YMzp0rcitKG2NiiBX8CLj55ptJT0+3i7CMMb4qVMEXkRFeBSnKBgwYQMWKFZk4caLfUYwxcaywPfwxnqQo4hITE7nllluYN28eR48e9TuOMSZOXbDgi8iWXB5bgSoRyFgk3H777Zw9e5bXXnvN7yjGmDglF1q+V0SOAj2BE8EvAWtVNfgGKDGvZcuWmpqa6vl+u3fvTlpaGvv27aNkyZKe798YE79EZIOqtszrPfkZ0lkElFXVA0GP/cBqL4J6SUSuEZGJIvKuiPRw2zqLyCoRmSAinf3K9tBDD/Hll18yZcoUvyIYY+LYBQu+qg5T1VW5vHadl2FE5DUROSYinwW19xKRnSKyW0QevkDeuao6AhgGDMlqBr4FLgLSvcxcEN27d+fKK6/kmWeesRujGGMiLtqmZU4GegU2iEgCzg1WegONgOtFpJGINBGRBUGPygGbjna3A1ilqr2B3wOPhf23yIWI8NBDD7Ft2zbmzJnjVwxjTJy64Bj+LzYQGaCq74YpDyKSBCxQ1V+7P7cDHlXVnu7PfwBQ1Sdy2V6A8cCHqvpR0GuJwBRVHZzDdiOBkQA1atRoceDAAc9+p0Dnzp2jSZMmAGzdupWEhISwHMcYE1+8GsMPlmOhDaPLgUMBP6e7bbm5F+gGDBaROwFEZJCIvAy8CbyQ00aq+oqqtlTVlpUqVfImeQ4SEhL4y1/+wvbt23nrrbfCdhxjjAlWPIRtxPMUBT9ern+WqOrzwPNBbbOB2Rc8kEg/oF9ycnJBMxbIoEGDaNGiBWPHjmXIkCFcdNFFYT2eMcZAaD38SJ9tTMe5k1aWasAX4ThQOFbLzImI8NRTT3HgwAGefjpm7/lujIkx0XbSNifrgboiUssdgx8KzAvHgbxeDz8vXbp04brrruOvf/0r+/btC/vxjDEmqgq+iEwF1gH1RSRdRG5T1bPAPcBiYDswQ1XTwnH8SPXwszz77LMkJCTw29/+1qZpGmPCLpSCf8rzFC5VvV5VL1PVEqpaTVVfddsXqmo9Va2jqo+H6/iR7OEDVKtWjbFjxzJ//nymT58ekWMaY+JXgadlxoNwLa2Qk7Nnz9KpUyd27txJWloal112WUSOa4wpWsI1LdN4qHjx4rz++uucPn2a22+/3YZ2jDFhE1LBF5G9IvJvESkR1F5RRPZ6Ey3yIj2kk6VevXqMHz+ehQsX8vzzz194A2OMCUFIQzoikokzNXI3MFBVT7jtVYAvVDWmLx+N5JBOFlVl4MCBvPfee6xcuZJ27dpF9PjGmNgWziEdBbq7//5HRMJ7pVIcEBEmT55MjRo1uPbaa8nIyPA7kjGmiAm14AtwHKforwQ+FpFOnqXyiV9DOlnKly/PO++8w/Hjxxk4cCCnT5/2JYcxpmgqTA8fVT2rqrcDT+HMk7/Fq2B+iPQ8/Jw0a9aMN954gzVr1jB8+HAyMzN9y2KMKVoK08PPpqpPATdg97j1xLXXXsv48eOZNm0af/rTn/yOY4wpIkJZPA3gN8DXgQ2qOldE2gCtCp3K8H//93/s3buXxx9/nAoVKnD//ff7HckYE+MKVPBF5IaAH69zlp7/hTOFSuSjSK2Wmc8svPjii5w4cYIHHniAiy66iLvuusvvWMaYGFagaZnudMxAyi+XL1ablumdM2fOMHjwYObNm8fEiRO5/fbb/Y5kjIlCnk/LVNVigQ/geyA5qD2mi320KVGiBDNmzKBXr16MGDGC5557zu9IxpgYZUsrxICSJUsyd+5cBg8ezP3338+YMWNsCQZjTIFZwY8RJUuWZNq0adx2222MGzeOUaNGcfbsWb9jGWNiSKizdIqkaDppm5OEhAQmTpxIxYoVefLJJ9mzZw8zZszgkksu8TuaMSYGeNHDLzJjC9Fw4dWFiAjjx4/n1VdfZcWKFbRr147PP//c71jGmBhQ0GmZHwQ1XQS8ISI/BDaqao/CBjN5u/XWW0lOTmbQoEG0bt2a119/nf79+/sdyxgTxQrawz8c9HgLZ8XM4HYTASkpKaxfv57atWszYMAAHnzwQc6cidnLIIwxYVagHr6qDg9XEBOaWrVqsXbtWh544AGeffZZ1qxZw7Rp06hZs6bf0YwxUcZm6RQBJUuW5IUXXmD69OmkpaVxxRVXMHnyZJu6aYw5jxX8IuS6667j008/pVmzZgwfPpwBAwZw5MgRv2MZY6KEFfwAfq+H74VatWqxdOlSnnvuOT788EMaN27MW2+9Zb19Y4wV/ECxMC0zP4oVK8bvfvc7Nm3aRN26dbnpppvo2rUrO3bs8DuaMcZHVvCLsAYNGrBmzRpeeuklNm3axBVXXMGYMWP44YcfLryxMabIsYJfxCUkJHDnnXeyY8cOhgwZwrhx42jQoAFTp061YR5j4ky+Cr6I1BeRriJSNqi9b3hiGa9VqVKFN998k2XLllGhQgVuuOEG2rZty5o1a/yOZoyJkAsWfBG5G3gXuA9IE5GBAS//OVzBTHh07tyZ1NRUJk2aRHp6Oh07dmTw4ME2vm9MHMhPD/8OoIWqDgCuAkaLSNb99nK85ZWJbgkJCQwbNoxdu3bx5z//mffff5/GjRtz00032bo8xhRh+Sn4xVX1OwBV3Qd0BnqKyN+wgh/TLr74YsaMGcO+fft44IEHmDVrFg0bNmT48OHs3bvX73jGGI/lp+AfEZGmWT+o6ingaqAi0CRMuUImIteIyEQReVdEerhtxUTkcRH5p4jc4nfGaFOpUiWeeuop9u3bx3333ce0adOoV68eN910E1u3bvU7njHGI/kp+DcD512uqapnVfVmIMXLMCLymogcE5HPgtp7ichOEdktIg/ntQ9VnauqI4BhwBC3eQBwOc4N1tO9zFyUVKlShb/97W/s2bOHe++9lzlz5nDFFVfQp08fVqxYYbN6jIlxFyz4qpquqtkFX0QGBLzm9RSPyUCvwAYRSQBeBHoDjYDrRaSRiDQRkQVBj8oBm452twOoD6xT1fuBUR5nLnKqVq3Kc889x8GDBxk3bhypqal07tyZtm3bMnPmTFuR05gYFco8/Cc8T+FS1ZXA10HNrYHdqrpXVX8CpgEDVHWrqvYNehwTx5PAIlXd6O4jHTjhPj+X07FFZKSIpIpIakZGRhh+u9hToUIFHnnkEQ4cOMCECRP4+uuvue6666hVqxbjxo3j6NGjfkc0xhRAKAU/0idqLwcOBfyc7rbl5l6gGzBYRO5022bjnGj+J7Ayp41U9RVVbamqLStVquRB7KKjVKlS3HHHHezYsYN58+bRuHFjxowZQ/Xq1bnxxhtZt26dDfcYEwNCKfiR/p+d0xdMrhlU9XlVbaGqd6rqBLfte1W9TVXvVdUXc9u2KCyeFk4JCQn069ePxYsXs3PnTu666y4WLFhA+/btadGiBf/61784ceLEhXdkjPFFLCytkA5UD/i5GvBFOA5UVBZPi4R69erx97//ncOHD/PSSy+RmZnJ3XffTdWqVbnxxhtZsmQJmZmZfsc0xgSIhYK/HqgrIrVEJBEYCswLx4Gsh19wZcqU4c4772Tz5s1s2LCB2267jYULF9KtWzfq1KnDn//8Zw4cOOB3TGMMoRX8U56ncInIVGAdUF9E0kXkNlU9C9wDLAa2AzNUNS0cx7cefuE0b96cF154gS+//JIpU6aQnJzM2LFjSUpKIiUlhQkTJvDVV1/5HdOYuCV2su1nItIP6JecnDzClhjwxv79+3n77bd5++232b59O8WLF6dnz57ccMMN9O/fnzJlyvgd0ZgiQUQ2qGrLPN9jBf+XWrZsqampqX7HKFJUlS1btjBlyhSmTp3KoUOHKF26NAMGDODaa6+lZ8+elC5d2u+YxsQsK/ghsoIfXpmZmaxZs4YpU6Ywc+ZMjh8/TunSpenduzf/8z//w9VXX82vfvUrv2MaE1PCWvBFpCfQFahM0LkAd9mFmGNDOpF35swZVqxYwezZs5kzZw5HjhwhMTGRHj16MGjQIPr378+ll17qd0xjol7YCr6IjAP+CGzBWWfnvJ2oau8C7zSKWA/fH+fOnWPdunXMmjWL2bNnc/DgQRISErjqqqvo27cv/fr1Izk52e+YxkSlcBb8Y8CDqvpGqOGimRV8/6kqGzZsYNasWcyfP5+0NGdiVv369enbty99+/alQ4cOlChRwuekxkSHcBb8I0BHVd0darhoZEM60Wvfvn0sWLCABQsWsHz5cn766SfKly9Pr1696Nu3L71796ZChQp+xzTGN+Es+H8CSqtqnksVxyrr4Ue3U6dO8dFHHzF//nzee+89jh07RrFixWjTpg09e/akR48etGrViuLFi/sd1ZiICWfBF2ABzjIHW3DWmc+mqrcWeKdRxAp+7MjMzCQ1NZUFCxawePFi1q9fj6pSvnx5unXrRo8ePejZsyc1atTwO6oxYRXOgv8X4BFgGzmftO1e4J1GESv4sev48eMsWbKExYsXs3jxYg4fPgxAgwYNsov/VVddxcUXX+xzUmO8Fc6CfwL4f6o6OcRsUcnG8IsWVWX79u3ZxX/FihWcPn2axMRE2rVrR5cuXejSpQutW7cmMTHR77jGFEo4C/5RoENRO2mbxXr4RdPp06dZtWoVH3zwAcuWLWPjxo2oKqVLl6ZTp07ZXwDNmjUjISHB77jGFEi4h3SKq+ofQg0Xzazgx4cTJ06wYsUKli5dytKlS7OnfpYrV47OnTtnfwE0btwY57SVMdErnAX/38D/APuAT/nlSduRBd5pFLGCH5+OHDnC8uXLs78A9uzZA0ClSpXo3LkzV111FSkpKTRu3JhixWJhZXETT8JZ8Jfl9bqq/qbAO40iVvANwIEDB1i2bBlLlixhxYoVHDrk3GmzQoUKdOrUiZSUFFJSUmjatKlNATW+s8XTCshO2prcqCoHDhxg5cqVrFixgpUrV7J7t3MKq2zZsnTo0IGUlBSuuuoqWrZsaSeBTcSFs4c/A9isqn8Nan8YaKqqQwu80yhiPXyTH1988QUrV67MfmSdAyhVqhRt27bNHgJq06aNLf1swi7cs3R6qurmoPYrgfdV9bIC7zSKWME3ocjIyGD16tXZfwVs3rwZVaVEiRK0bt2alJQUOnXqRPv27bG7qhmvhbPgnwZ+HTwtU0TqAltV9aIC7zSKWME3Xjh58iRr1qzJ/gJITU3l7NmzFCtWjCuvvDL7HEDHjh2pXLmy33FNjAtnwU8DXlDVl4La7wLuVdWGBd5pFLGCb8Lhu+++4+OPP2bVqlWsXLmSjz/+mB9++AFwrgTOOhHcqVMnatas6XNaE2vCWfDvBp4AxgIf4iyt0BN4FHhEVf9Z4J1GESv4JhJ++uknNmzYkP0FsHr1ak6ePAlAjRo1zpsJVL9+fbsWwOQp3He8egx4CCjpNv0IPKuqY0LaYRSxgm/8cO7cOT777LPsL4BVq1Zx5MgRwLkWoFOnTtlfAldeeaVdDWzOE/ZpmSJSGmjs/rhNVb8LeWdRwKZlmmiiquzevTu7+K9cuZJ9+/YBP08FzfoCaNWqFSVLlrzAHk1RZvPwQ2Q9fBOt0tPTWbVqVfYXQNZU0JIlS9KmTZvscwDt2rWjbNmyPqc1kWQFP0RW8E2sOH78OKtXr87+Ati4cSPnzp0jISGBZs2aZX8BdOzYkYoVK/od14SRFfwQWcE3serbb79l3bp1580E+vHHHwFo1KhR9hdASkoK1apV8zmt8ZIV/BBZwTdFxY8//khqamr2eYDVq1dz6tQpAJKSks77Aqhbt67NBIphVvBDZAXfFFXnzp1jy5Yt550IzsjIAKBKlSrnXQvQpEkTmwkUQ6zgh8gKvokXqsquXbvO+wI4cOAA4NwXoFOnTnTt2pWuXbvy61//2v4CiGKRmJY5QlUnhryDKGUF38SzgwcPZhf/ZcuWkTVFuUqVKnTp0oVu3brRtWtXuxo4ykSi4B9U1Roh7yAMROQa4GqgMvCiqn4gIp2AG4HiQCNVbZ/XPqzgG/OzgwcPsmTJEpYsWcJHH33E0aNHAUhOTs7u/f/mN7+xWUA+86Tgi8iW3F4C6qmqZ1d7iMhrQF/gmKr+OqC9F/APIAH4t6qOz8e+LgGeUdXbAtquAaqo6st5bWsF35icqSrbtm3LLv7Lly/n1KlTiAhNmzale/fu9OrViw4dOtg9ASLMq4J/FGednBPBLwFrVbVqoVKef6wU4FvgjayCLyIJwC6gO5AOrAeuxyn+TwTt4lZVPeZu9yzwtqpuDNj/DOB2Vf1vXjms4BuTP2fPnmX9+vXZXwBr167lzJkzlC1blq5du9K7d2969+5N9erV/Y5a5OWn4OfnvmyLgLLBa9+7B1gdYrYcqepKEUkKam4N7FbVve4xpwEDVPUJnL8GgjMJMB5YFFTsawAnL1TsjTH5V7x4cdq1a0e7du0YPXo0p06dYunSpSxatIhFixYxd+5cABo3bkzv3r3p06eP9f59FHWzdNyCvyCghz8Y6KWqt7s/3wS0UdV7ctn+PuAWnL8ENqvqBLf9MWCxqq7NZbuRwEiAGjVqtMiaqWCMCY2qsn37dhYuXMiiRYtYtWoVZ86coUyZMnTr1o3+/fvTt29fKlWq5HfUIiEsJ21FZICqvluoZHnvP4nzC/61OHfXCiz4rVX13nBlsCEdY7wX2PtfuHAhhw4dolixYnTs2JFrrrmGAQMGULt2bb9jxqz8FPxiIew3eNw83NKBwAHAasAX4TiQiPQTkVey1iQ3xninbNmyDBgwgAkTJnDgwAE2btzI6NGj+eabb7j//vupU6cOV155JWPHjmXTpk1E2+hDURBKD397OO9olUMPvzjOSduuwGGcoZobVDUtXBmsh29MZO3du5d3332XuXPnsnr1ajIzM6lZsyZDhgxh6NChNG3a1C76uoBwDelsU9VGhUqW+76nAp2BisBRYKyqvioifYC/48zMeU1VHw/T8W09fGN8lpGRwfz583nnnXf48MMPOXv2LPXq1WPo0KEMHTqUhg1j+g6qYRNzBT9aWA/fmOjw1VdfMXv2bKZNm8by5ctRVa644gqGDh3KDTfcYFf7BrCCX0DWwzcmen355Ze88847TJs2jbVr1yIidO3aleHDhzNw4EBKlSrld0Rfhavg/0dV2xQqWZSzHr4x0W3//v28/vrrTJ48mf3791OuXDmuv/56hg8fTqtWreJyvN9Wyywg6+EbE1syMzNZvnw5kyZNYtasWfzwww80btyYu+66i5tuuimubvNoBT9E1sM3JvacPHmS6dOn8/LLL7Nx40bKli3LsGHDuPvuu6lfv77f8cIuXPPwEZHuItIh4OcRIrJeRCaLSPx8pRpjoka5cuUYOXIkqamprFu3jv79+zNhwgQaNGhAjx49eP/99+N+bn9IBR94CmfqJCJSD3gRSAVaAk97E80YYwpORGjbti1vvfUWhw4dYty4caSlpdG7d2+aN2/O9OnTOXfunN8xfRFqwa8DfOY+Hwh8pKqjgBHksKBZrLArbY0pWqpUqcIjjzzCvn37eO211/jhhx8YOnQo9evX5+WXX+b06dN+R4yoUAs+QNbfRlcBH7jPDwOXFiqRj1R1vqqOLFeunN9RjDEeSkxMZPjw4aSlpTFr1iwuueQS7rzzTurUqcOECRM4c+aM3xEjItSCvwUY5a5f34WfC351IMOLYMYY47WEhAQGDRrEJ598wkcffUStWrUYNWoUDRs2ZMqUKWRmZvodMaxCLfgPA7cCy3BuVrLNbe+Hs9ZNTLIhHWPiQ9ZFW6tWrWLBggWUKVOGG2+8kWbNmrFs2TK/44VNSAVfVVcDlYBLVXVkwEsTgbu9COYHG9IxJr6ICFdffTUbN25k6tSp/Pe//6VLly4MGTKEQ4cO+R3PcyGP4atqpqp+E9S2R1WPFDqVMcZEULFixRg6dCjbtm3jscceY968eTRo0IDHH3+cH3/80e94ninMSVtjjClSSpUqxZ/+9Ce2b99Or169GD16NM2bN2f9+pgdqT6PFXxjjAmSlJTErFmzeO+99zh58iRt27bl97//fcxP47SCH8BO2hpjAvXp04e0tDRuvfVWnnrqKdq1a8euXbv8jhUyK/gB7KStMSZYuXLlmDhxIvPnz+fQoUM0b96cN9980+9YIbGCb4wx+dC3b182b95MixYtuPnmm7nnnnti7oKt4qFsJCKf8/OVtoEUOI1zD9qJqvpBDu8xxpiYVK1aNZYuXcrDDz/MM888w44dO5gxYwYVKlTwO1q+hNrDnwlUBo4DC9zHV27bUqACsMhdX94YY4qMhIQEnn76aSZNmsSqVato06YNe/bs8TtWvoRa8MsB/1LVdqp6v/toj7NqZglV7Qo8A4z2KqgxxkSTYcOGsXTpUr7++mtSUlLYsWOH35EuKNSCPxSYlEP768AN7vM3gQYh7t8YY6Jehw4dWL58OefOnSMlJYVPP/3U70h5CrXgJwD1cmivF7DPH4GYWonIpmUaYwqqSZMmrFy5kpIlS9K5c+eoLvqhFvwZwL9FZJiINBSRBiIyHHgFmOa+px0Q/X/jBLBpmcaYUNSrV49Vq1Zx8cUX06dPn6hdhyfUgn8fzonbCTg3QtkGvAS8A/zWfc9G4PbCBjTGmFiQlJTEokWL+Pbbb+nduzfffPON35F+IdTVMk+r6n04s3GaAU2BCqr6W1X90X3PZ6qa5llSY4yJck2aNGHOnDns2rWLgQMHRt3CawUu+CJSQkSOikhjVf1eVbe4j+/DEdAYY2JJly5dmDRpEsuXL+fBBx/0O855ClzwVfUMcM59GGOMCXLjjTdy//3388ILLzB9+nS/42QLdQz/3zjj+MYYY3Iwfvx42rZtyx133EF6errfcYDQC35V4EYR2SEib4vIK4EPLwMaY0wsKlGiBG+++SZnzpxhxIgRqOa0Gk1khVrw6+DMwvkSp/jXDXgkexPNGGNiW3JyMk8++STvv/9+VKywKdHwreMlEbkGuBpnXZ8XVfUDEakBvICz3s8uVR2f1z5atmypqampYc9qjCn6MjMzad++PQcOHGDnzp386le/CstxRGSDqrbM6z2FXh5ZRC4VESnsftx9vSYix0Tks6D2XiKyU0R2i8jDee1DVeeq6ghgGDDEba4HvKeqtwKNvMhqjDH5UaxYMf75z39y5MgR/vKXv/ibJZSNRCRBRB4TkRPAUaCW2z5eRO4oRJ7JQK/gY+EsytYbp1hfLyKNRKSJiCwIelQO2HS0ux3AJmCoiCwFlhUinzHGFFirVq0YPnw4zz//vK8ncEPt4f8euAVnps5PAe2bcHrWIVHVlcDXQc2tgd2quldVf8JZumGAqm5V1b5Bj2PieBJYpKob3X0MB8aqahec4Z5fEJGRIpIqIqkZGRmh/grGGJOjsWPHkpmZyfjxeY4oh1WoBf8W4E5VfZPz5+NvJedF1QrjciBwYYp0ty039wLdgMEicqfb9j5wn4hMAPbntJGqvqKqLVW1ZaVKlQqf2hhjAtSsWZNbb72ViRMn+tbLD7Xg1wC259B+FigVepwc5XR+INczzar6vKq2UNU7VXWC2/aZqg5223K99M1WyzTGhNMf//hHzp07xz/+8Q9fjh9qwd8PXJlDe3e8XyEzHage8HM14AuPjwHYapnGmPCqWbMm11xzDZMmTeL06dMRP36oBf9fwD9EpIf7c10RuQt4HPinJ8l+tt7dfy0RScS5+co8j48BWA/fGBN+o0aN4vjx48ycOTPixw55Hr6IPAo8xM9DOKeBv6rquJDDiEwFOgMVcWb/jFXVV0WkD/B3nBuvvKaqj4d6jPywefjGmHBRVRo0aECFChVYt26dZ/vNzzz8Ql14JSKlgMY4fymkqep3Ie8sCrg3Xe+XnJw84vPPP/c7jjGmiHrmmWd46KGH+Pzzz0lO9mZxAk8vvBKRFsFtqvqDqqaq6ieBxV5ELhKRhgWL6z8bwzfGRMLgwYMBmDNnTkSPW5Ax/HdFZI6I9BSRHLcTkcvdK2E/Bzp4ktAYY4qYpKQkmjdvzuzZsyN63IIU/Po4tzN8CzgpIitFZKqITBKRd0VkN3AQZ6bO9ar67zDkDSs7aWuMiZRBgwbx8ccfc/jw4YgdM98FX1W/U9UxONMibwJSgYuAy4D/4ixj0FhVu6rq6nCEDTcb0jHGRMqgQYMAmDt3bsSOWeRWy/SCzdIxxkRCnTp1aNKkiSdFPyKrZRYlNqRjjImkjh07snbt2ojdHMUKfgAb0jHGRFKHDh3IyMhg9+7dETmeFXxjjPFJ+/btAVi7dm1EjmcF3xhjfNKoUSPKly/PmjVrInI8K/gBbAzfGBNJxYoVo127drFX8EWkontRVk6raMYEG8M3xkRa+/bt2bZtGydOnAj7sbzs4c8D+gN/E5EPgm43aIwxJgctWjir1qSlpYX9WJ4VfFVtr6p3q2pXnOWTF4tI9QttZ4wx8axOnToA7N27N+zH8nQMX0RKiEhNnKWNFwFLvdy/McYUNUlJSYgIe/bsCfuxinu1IxHJcPf3pfv4Apjl1f4jIWB5ZL+jGGPiRGJiItWrV49ID9+zgg/UVNXvPdxfxKnqfGB+y5YtR/idxRgTP+rUqRORHn6BhnREZJGIjBORgcHj87Fe7I0xxi+1a9eOyh7+GWAY8EdAReQ4sAFn5cwNwAZVPeRpQmOMKeLq1KnD0aNH+fbbbylTpkzYjlOgHr6q9lfVajhLIvfDuWH5jzhfArOB/SLypbtGfhOvwxpjTFFUu3ZtAPbt2xfW44Q0hq+qR4GF7gMAd959S6AV0BtYLyLdVXWVF0GNMaaoypqauWfPHpo0CV9f2ct5+MdUdaGqPqaqbYG/Ak94tX9jjCmqsnr44R7HD+daOm8AMbXMgq2lY4zxQ4UKFShfvnxMF/yDOOP8McPW0jHG+KV27dphn5rp5Tz886hqJrA8XPs3xpii5JJLLuHUqVNhPYYtj2yMMXHCCr4xxsQJK/jGGBMnrOAbY0ycsIJvjDFxosgVfBG5RkQmisi7ItLDbWskIjNE5CURGex3RmOM8UNUFXwReU1EjonIZ0HtvURkp4jsFpGH89qHqs5V1RE46/sMcZt7A/9U1VHAzeHIbowx0S5s8/BDNBl4AecqXQBEJAF4EegOpOOs0TMPSOCXSzfcqqrH3Oej3e0A3gTGikh/4NKwpTfGmCgWVQVfVVeKSFJQc2tgt6ruBRCRacAAVX0C6Bu8DxERYDywSFU3uvs9BtztfnnMzunYIjISGAlQo0YNb34hY4yJIlFV8HNxORC4xn460CaP998LdAPKiUiyqk5wv0T+CFwMPJ3TRqr6CvAKOLdrFJEDBchYEfiqAO+PFpY7six35MVcdqfPGlLumhd6QywUfMmhTXN7s6o+Dzwf1LYft/eeH6paKb/vBRCRVFVtWZBtooHljizLHXmxmj1cuaPqpG0u0oHA2ylWw7lBujHGmAKIhYK/HqgrIrVEJBEYCszzOZMxxsScqCr4IjIVWAfUF5F0EblNVc8C9wCLge3ADFVN8zNnDl7xO0CILHdkWe7Ii9XsYcktqrkOhxtjjClCoqqHb4wxJnys4BtjTJywgl8IBVnyIVJEZL+IbBWRzSKS6rZVEJEPReRz999LAt7/Bzf/ThHpGdDewt3PbhF53r2gzcucv1hGw8ucIlJSRKa77f/J4YI+L3M/KiKH3c98s4j0icLc1UVkmYhsF5E0Efmt2x7Vn3keuWPhM79IRD4RkU/d7I+57f595qpqjxAeOEs77AFqA4nAp0CjKMi1H6gY1PYU8LD7/GHgSfd5Izd3SaCW+/skuK99ArTDuQ5iEdDb45wpQHPgs3DkBO4CJrjPhwLTw5j7UeDBHN4bTbkvA5q7z8sCu9x8Uf2Z55E7Fj5zAcq4z0sA/wHa+vmZ+1qcYvnhfviLA37+A/CHKMi1n18W/J3AZe7zy4CdOWXGmQnVzn3PjoD264GXw5A1ifMLp2c5s97jPi+Oc9WihCl3bsUnqnIHZXsXZ32qmPjMc8gdU585UBrYiLNKgG+fuQ3phC6nJR8u9ylLIAU+EJEN4qwPBFBFVb8EcP+t7Lbn9jtc7j4Pbg83L3Nmb6PO1N6ThHfhvHtEZIs75JP1J3pU5nb/7G+G0+OMmc88KDfEwGcuIgkishk4Bnyoqr5+5lbwQ1egJR8iqIOqNsdZEvpuEUnJ4725/Q7R9ruFkjOSv8NLQB2gKfAl8OwFMviWW0TKALOA36nqf/N6ay45fMmeQ+6Y+MxV9ZyqNsVZIaC1iPw6j7eHPbsV/NBF5ZIPqvqF++8xYA7OaqNHReQyAPffrCWkc/sd0t3nwe3h5mXO7G1EpDhQDvg6HKFV9aj7HzsTmIjzmUddbhEpgVM031bVrFVjo/4zzyl3rHzmWVT1G2A50AsfP3Mr+KGLuiUfRORiESmb9RzoAXzm5rrFfdstOOOguO1D3TP9tYC6wCfun5mnRKStOxvg5oBtwsnLnIH7GgwsVXeg02tZ/3ldA3E+86jK7R7nVWC7qv4t4KWo/sxzyx0jn3klESnvPi+Fs4rvDvz8zL08MRFvD6APzqyBPcAjUZCnNs5Z/k+BtKxMOGN6S4DP3X8rBGzziJt/JwEzcYCWOP+J9uDclMbrk1hTcf4UP4PTS7nNy5zARcBMYDfODIfaYcz9JrAV2OL+B7wsCnN3xPlTfwuw2X30ifbPPI/csfCZXwFscjN+BvzJ6/+PBc1uSysYY0ycsCEdY4yJE1bwjTEmTljBN8aYOGEF3xhj4oQVfGOMiRNW8I0xJk5YwTfGmDhhBd+YICLykIhs9XifZdz121t5ud+A/b8sIs+EY9+m6LCCb8wvNce5otNLvwdSVXV9fjcQkXdF5JNcXksUka9EZJzb9GdglIjU9iCrKaKs4BvzS54WfBG5CBgFvFzATV8GWonIlTm8Ngi4BPg3gKoexrlM/65CRDVFnBV8YwK4y/AmA1+LyAwR+a+IHBWROwLeM1NEFgf8XElEvheRFrnsthdQCvggh+PdKyI7ROS0e8u7R9xVDwHeBw4AI3LY5wic9dX3B7TNAf63AL+uiTNW8I05XzOc/xf3AVNw1lt/A3jBXYEUgpaxVdUMnJtyXJPLPq8CNqlzg4psIvIo8CDOnY4aAr8F7gDGuvvNxOnB/6+72mLWdnWA3wCvBB3nP0AVEWlYgN/XxBEr+MacrxnOSphDVHWuqu7FKfjFcYZQ4JfrlgN8y893LgpWCzgc2CAipYH/A+5Q1Tmquk9VFwKjgXsD3voqcDFwbUDb7cBRfrkcd9ZdkWwc3+TICr4x52sOLFHVXQFtdYHvOf+mE2VEpBxk33ugI7Ail32WAk4HtTV222eJyLdZD5xx+3IiUgmyb4G3AHdYxx3uGQZMCv6LIeAYpTAmB8Uv/BZj4kpzYEYObVvcIRb4uSddHeceok8AR4DZ5CwDqBDUltXZuhbnngrBAu9a9DKwyB2qaQhUwT1ZGyTrGBm55DBxzgq+MS53Nk1DYGPQS82D2rIKfjV3Xv1NQEdVDe7FZ9kI3BPUlobTI6/tDuXk5QNgP04vvyHwkTvUFKwJcA7nphvG/IIVfGN+dgXO/4kNQe3Nce6pmuUwkIlzkrU10ENV0/LY7yLgWRGprqqHAFT1WxH5K/BX5651fOgeuwnQTFV/n7WxqmaKyEScMf+ywJBcjtMZWK1535zcxDEbwzfmZ82AL1T1aFaDiFyOM4SS3cN3x86PAlcCXS90MZWqbse5gfVNQe1/Af4fzknYT4HV7s/7c9jNazgnb78ih/sLu/c6vYGCz/U3ccRucWhMBIhIJ2AaUFdVvw/D/q8DxgBNVfWc1/s3RYP18I2JAFVdBTyGM0UzHEoCw63Ym7xYD98YY+KE9fCNMSZOWME3xpg4YQXfGGPihBV8Y4yJE1bwjTEmTljBN8aYOPH/AYBi6IhAW3jTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.semilogy(E, flux[0] * E * erg_per_ev, color='k')\n",
    "pl.xlabel(ares.util.labels['E'])\n",
    "pl.ylabel(ares.util.labels['flux_E'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6e4020a",
   "metadata": {},
   "source": [
    "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)\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 = -1.5$, $\\beta = 0$, $z=5$, and $z_i=60$,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "952fbee5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1873fdcd0>]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAAD4CAYAAAA5FIfVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoU0lEQVR4nO3deXRUZbb38e/OwIyICAmEQAggCDZCCPOgDDJzodEWaK4KMrSt2NrX9jaoV+2X7kXra0srerVB0EZlElEBRWRSRplnZAhzICQICBiGQLLvH6eCISYxkKo6VZX9WatWqp6cYddZ8vPJc55zjqgqxhhj3BXmdgHGGGMsjI0xJiBYGBtjTACwMDbGmABgYWyMMQEgwu0Crsett96qcXFxbpdhjDHXZcOGDd+rauWClgmqMI6Li2P9+vVul2GMMddFRA790jJBMUwhIr1FZMKZM2fcLsUYY3wiKMJYVeeq6ogKFSq4XYoxxvhEUISxMcaEOgtjY4wJABbGxhgTACyMjTEmAFgYG2NMAAiKMC7K1LZZs2Yxbtw4H1RljDHeExRhXJSpba+9dpLnn7/Z+0UZY4wXBdUVeDfi4sXb+fHHBFQVEXG7HGOMyVNQ9IyLokoVgHKkpv7odinGGJOvkA/jqlXDAdi9+5TLlRhjTP5CPoyrVy8JQFLSWZcrMcaY/IV8GMfHlwW+JyXFwtgYE7hCPow7dCgHVCYqaqfbpRhjTL78FsYi0ldEJorIZyLSJUd7WRHZICK9fLHfypWd+zmnpaX5YvPGGOMVhQpjEZksImkisj1XezcR2S0iSSIyqqBtqOqnqjocGAz0z/GrPwMzr7PuQitVqhSRkVNZuPA2X+3CGGOKrLDzjN8D3gCmZDeISDjwJnAPkAysE5E5QDgwNtf6D6tqdtf0Oc96iEhnYCdQ6gbrLxSRliQlHfPlLowxpkgKFcaqukxE4nI1NweSVHU/gIhMB/qo6ljgZ0MO4lxx8Xdgvqpu9DR3AMoCDYALIvKFqmbd0DcpQLlyP/DDD+W9vVljjPGaolyBFwMcyfE5GWhRwPKPA52BCiJSR1XfVtVnAURkMPB9XkEsIiOAEQA1atS4oUJvvfU8SUm1bmhdY4zxh6KEcV7XFmt+C6vq68Dr+fzuvQLWmyAiKUDvEiVKNL3eIgGqVctkz55ozp27SPnyPh0RMcaYG1KU2RTJQGyOz9UBnwzMFvUZeLfdpsAWdu60cWNjTGAqShivA+qKSC0RKQEMAOZ4pyzv6t8/E0jg/PlffFq2Mca4orBT26YBq4F6IpIsIkNV9QowElgAfAfMVNUdviiyKPczBoiNdTrwhw5ZGBtjApOo5jvMG3ASExN1/fr1173ehQsZlCmzlnvuOctXX/XwQWXGGJM/EdmgqokFLRMUl0MXtWdcunQJIiLi2L3bprcZYwJTUIRxUU/gAdx8cxonTlT0YlXGGOM9QRHGRe0ZA8TEpHPhQnUyM71+TYkxxhRZUISxN3rGt90mwM1s23bce4UZY4yXBEUYe0PbtqWBOezYccDtUowx5meKTRjff38M0Ifvv7/+2RjGGONrQRHG3hgzjo6OJioqio0bt3mxMmOM8Y6gCGNvjBkDiExj1qw/eKkqY4zxnqAIY2+pXj2S8+frcOFChtulGGPMNYIijL0xTAGQmBgBlGHBgv3eKcwYY7wkKMLYW8MUXbtGATB/fqo3yjLGGK8JijD2lh49agI/smLFFbdLMcaYaxSrMC5RIoz69Wdw9uxst0sxxphrFKswBnjwwTSSk/+XkydPul2KMcZcFRRh7K0TeACtW7cBbmPu3I2/uKwxxvhLUISxt07gAdSv3wzYyYQJRa/LGGO8JSjC2JuiokpTvvwetm6t7HYpxhhzVbELY4BmzU6Rnv4rvvvO7uBmjAkMxTKMH3igMhDOG2/scrsUY4wBimkYDxpUF5EzzJ8fPM//M8aEtmIZxpGRQvfub3Dy5MNkZNh9Kowx7guKMPbm1LZsv//9nZw9e5DFixd7bZvGGHOjgiKMvTm1Lds999xDqVJ/5MUXT3ltm8YYc6Mi3C7ALSVLlqRKld+ybl0MFy9mUKpUCbdLMsYUY0HRM/aV/v3DUK3KuHH2KCZjjLuKdRiPHn0HIud4662LbpdijCnminUYV6xYgkaNdnDkSEt277Z7HBtj3FOswxjg+eerAmuZMOFTt0sxxhRjxT6M+/WrSfv2L/Dppy+TmZnpdjnGmGKq2IcxwMiRI9m/P5233lrkdinGmGLKb2EsIn1FZKKIfCYiXTxtd4vIchF5W0Tu9lctufXt24+IiPU899xNqNol0sYY/ytUGIvIZBFJE5Htudq7ichuEUkSkVEFbUNVP1XV4cBgoH92M/AjUApIvu7qvSQyMpxu3ZI5c6YVU6ZscqsMY0wxJoXpCYpIe5zQnKKqd3jawoE9wD04QboOGAiEA2NzbeJhVU3zrPcP4ENV3SgiYaqaJSJRwKuqOqigOhITE3X9et/MCU5OPk9srBIdvZaUlA4+2YcxpngSkQ2qmljQMoXqGavqMiD3dcPNgSRV3a+qGcB0oI+qblPVXrleaeJ4CZivqhs9283ybOs0UDKfLzFCRNaLyPoTJ04UptwbUr16GTp23M7x4x2YNMkuAjHG+FdRxoxjgCM5Pid72vLzONAZuE9EHgEQkX4i8i/gfeCNvFZS1QmqmqiqiZUr+/bpHO+/3wiR0/z1ryts7NgY41dFCWPJoy3fBFPV11W1qao+oqpve9pmq+rvVLW/qn6d7458cNe2vFSrVppXX/2Mgwf/yNy5c326L2OMyakoYZwMxOb4XB04VrRy8uaLu7bl57HHBlG/fn0ee2wi6el2mbQxxj+KEsbrgLoiUktESgADgDneKeta/uoZA0RGRjJy5BSSkz9j4EC717Exxj8KO7VtGrAaqCciySIyVFWvACOBBcB3wExV3eGLIv3ZMwZ49NFmVKmyk7lzW7FhwyG/7NMYU7wVamqb20SkN9C7Tp06w/fu3euXfS5dmkrHjrcQFbWMY8c6EBZmFysaY26M16a2uc3fPWOADh2i6NNnC6mpnRg2bIHf9muMKZ6CIozd8tFHTalYcRsffPA5e/bscbscY0wIC4ow9ucJvJwiI4WtWytRrtxUBgwYwMWLNrvCGOMbQRHGbgxTZKtevRpTpkxh06ba3HXXPL/v3xhTPARFGLutV69e3HHHaNau7cfTT3/ldjnGmBAUFGHs1jBFTitWNKJs2UO88kpTpk/f6FodxpjQFBRh7OYwRbYKFSJYurQCYWHCoEE3s27dAddqMcaEnqAI40DRrNktTJ2aTlZWNXr1msLJkyfdLskYEyIsjK9T//6xTJ++nTNnxtK9e3fcHDoxxoSOoAjjQBgzzql//0Q++ugjNm7Mol69pfzww49ul2SMCXJBEcaBMGacW+/evXnooYmkpvalfv1VnD2b7nZJxpggFhRhHKgmTWrCb36zmdTULsTHf0tq6mm3SzLGBCkL4yKaObMxDz64lZMnO1Gnzk727TvudknGmCBkYewF//53I556aicXLpzjnnvutvtYGGOuW1CEcaCdwMvLK680YMWKipw7d5JmzXozefK3bpdkjAkiQRHGgXgCLy8tW7Zg7dq1qP4vQ4fezuOPf+F2ScaYIBEUYRxMatWqxbffNuemm77njTe60rz5bNLTL7hdljEmwFkY+0CDBuU5ciSO22/fzrp1/ahadQNr1iS5XZYxJoBZGPvITTeFs2PHnYwcuZ3z52+hU6e7mDZtmttlGWMClIWxD4nA+PF3sHdvee68M47f/vZhmjSZSkrKKbdLM8YEGAtjP6hVK5ZvvvmG/v0/YPPm31KjRiqvvbbc7bKMMQEkKMI4GKa2/ZKIiAimT7+X8eP3AZV48smWNG48l6NH7c5vxpggCeNgmdpWGCNH1ubw4Qo0bLidLVt6U7v2Uv7973+jqm6XZoxxUVCEcaipWrUk27c3YcKE/dSv/wmDBw+mdeu+fP21XblnTHFlYeyi4cPj2bjxfSZOnMjGjQPo0KEaLVt+waFD37tdmjHGzyyMXRYWFsawYcNYvrwL8fF7WLOmB7VqZdC//5ecP3/J7fKMMX5iYRwgmjevxL59Cbz//gEqVDjLzJndqFbtbaZOnUpmZqbb5RljfMzCOMD853/W4tSp+jzzzDaqVZvLoEGDqFu3P888s5zMzCy3yzPG+IiFcQASgb/97Vds3/4VM2bM4OTJBxk7th3ly2/j+edXkpVlMy+MCTV+C2MR6SsiE0XkMxHp4mkLE5G/ich4EXnIX7UEi7CwMO6//36OH+/J0KFruHz5VsaMaUO5cjv4r/9ayuXLl90u0RjjJYUKYxGZLCJpIrI9V3s3EdktIkkiMqqgbajqp6o6HBgM9Pc09wFigMtA8nVXX0yULh3OO++04MyZKIYM+ZbMzLKMG7eQOnXq8Oqrr5Oaag9ENSbYFbZn/B7QLWeDiIQDbwLdgQbAQBFpICK/EpF5uV5Vcqz6nGc9gHrAalX9L+D3RfkixUGZMhFMntyS9PSazJ7dlpo1a/LUU6upWjWDtm2/Yf36VLdLNMbcoEKFsaouA3Lf3aY5kKSq+1U1A5gO9FHVbaraK9crTRwvAfNVdaNnG8lA9lM885wyICIjRGS9iKw/ceLEdX/BUBQREcavf92DZcuW8e67o4mO3svKle1o1uxWYmPX8s9/brNxZWOCTFHGjGOAIzk+J3va8vM40Bm4T0Qe8bTNBrqKyHhgWV4rqeoEVU1U1cTKlSsXodzQNHhwI44da8GSJYdp1uwbjh6N549/DKNJkya88847/PDDebdLNMYUQlHCWPJoy7c7pqqvq2pTVX1EVd/2tJ1X1aGq+riqvpnfuqFwoyBf69AhjrVrO3LiRGnGjNkBKMOHP8Ett6RTr94G3nprHzZd2ZjAVZQwTgZic3yuDhwrWjl5C6UbBflapUplee65+9m8eTNz5y6mTp2N7NkTx6OP1qZ06eN067aBHTt+cLtMY0wuRQnjdUBdEaklIiWAAcAc75R1LesZXz8RoVevluzZ05Vjx4SHH55P6dJ7WbCgMQkJ3Rg0aBBTp67g+HG7kMSYQFDYqW3TgNVAPRFJFpGhqnoFGAksAL4DZqrqDl8UaT3joqla9RYmTerOmTPtWLhwJ8OGJfLFF18waNAeqlZV4uJ288ILhzh1yk76GeMWCYb76IpIb6B3nTp1hu/du9ftckLCxYsXeeONpUyYcI69exOAOkAGCQnbmT69PHXr1nW7RGNChohsUNXEApcJhjDOlpiYqOvXr3e7jJBz8uQpXn31G6ZMySA5eRvwNxITm3P58hTuu+9WhgypRExB82SMMQWyMDbXLTk5mRkzZvDuu4vYseM14DYAYmOP8pvflODJJysTG1vwNowx1wqZMLZhCnckJe3j7be/YcaMSyQnNwWaU6PG73jooSiaNRtIeHh9OnQQSpd2u1JjAlvIhHE26xm758iRI7z77kIWLZrKypVLycp6AXie8PAMmjY9y4ABN9O7dwR16rhdqTGBx8LY+ERaWhoffTSP9947wKZNUWRmdgFuIzz8IhMmzKJ3766kpVUmJgZuvtntao1xX8iEsQ1TBK709HQWL17Mhx+uYeHCZE6fnoKIULr0bi5cqE3Dhhfp2bM0nToJbdpAmTJuV2yM/4VMGGeznnFgy8rKYtOmTcybN49p05LZvTsG6AS0BCJp3fow8+aVp2LFiqxdC3feCSVLuly0MX5gYWxcdezYMT7//HPmzfuaRYsucf78UcLC1nLnnf/Bpk2fUKJEJi1ahNG+vdCuHbRuDeXLu121Md5nYWwCxpUrV1i7di0LFizgiy+WsWFDJVRbEx5+N1lZd6Iazj//eYInnqjM4cOwdi20awdRUW5XbkzRhUwY25hx6Dl16hSLFy/mq6++Yv785Rw9Wh3YRJ06txAV9f9YuXIgAHXqQMuWzuuBB+Cmm9yt25gbETJhnM16xqFJVdm1axcLFixgyZIlfP31Ks6dqw20o3z5HmRmNuXChZs4cOAMNWvezHvvwdatTkC3aAE1ajgPcTUmUFkYm6B05coVNm7cyJIlS1i6dCnLli3n4sWbCQtLJSEhAZF/sGVLGzIywgFnKKNTJ/jwQ2f9CxewC1FMQLEwNiHh0qVLrF27liVLlrBkyRJWr17N5ctKREQCNWveT5kyHYmOrsLs2RUoV64cTZvC999DQsJPr8REG3827rEwNiEpPT2dVatWXQ3nDRs2kJmZSXh4OE2bNqVcudFcudKMlJRokpLCUYXf/AZmznTWf+kluP12J6RjYmyIw/heyISxncAzBfnxxx9ZvXo133zzDcuWLWPNmjVkZGQgIjRs2Ip69e6nWbM7GDLkV4SFVSEqCrI899S/5RZo1Aj+8Af49a8hMxMuXbKLU4x3hUwYZ7OesSmMixcvsmbNGpYtW8ayZctYtWoV5887D2atX78+rVvfQ2xsLyIimnHkSEW2bIEnnoCBA2HzZqfHXLeuE9LZr3btnOA25kZYGBsDZGRksHHjxqvhvHz5cs6ePQtAfHw87du3p3379rRr146IiNr8+9/C1q3OjI19+0AVFiyALl1g9Wp4911o0OCnlw11mF9iYWxMHjIzM9m6dSvLli27OrRx8uRJAKpWrUq7du2uBnTNmg357rswGjRwrg784AN48knwLA447Vu2QK1aTs/62DEnpGvUgLCiPGXShAwLY2MKISsri127drF8+fKrvefk5GQAKlasSNu2ba8GdEJCApGRkZw4ATt3/vT6xz+gRAln7Hn8eGe7ZcpA/fpOML/7LkREOLM8ype3e3IUNxbGxtwAVeXQoUPXhPOePXsAKFOmDK1atboazi1atKBMjrN9p09fG9Lffee0rVnj/L5fP/jsM4iLg9tuc15NmsDgwf7/nsZ/LIyN8ZLU1FSWL19+NaC3bNmCqhIZGUliYuLVcG7Tpg03F3AT5y++cIJ5927Ys8d5NWz4U1jffbcT3tlBfdtt0Lixc4c7E7xCJoxtapsJND/88AOrVq26Gs7r1q3j8uXLiAiNGjW6ekKwXbt2REdH57sdVTh7FipUcD6/+CJs3OiE9f79cOUK3HsvzJrl/L5bN+eG/bVr//S6/XaoUsXnX9kUQciEcTbrGZtAdf78edauXXt1tkbO6XR169a9Gs7t27cnLi4OKcT0i8uX4eBBZ050vXpOMPfsCUlJcOiQMycanBOK48Y586Pvuw/i451XdljXqgWlSvnuu5tfZmFsjEsuX77Mpk2brobz8uXLOX36NAAxMTHXhPPtt99O2HVOu7hyBQ4fdqbeVa0Kd9zhzOLo0cNp+/HHn5Z95RV46ik4ehT+53+ccI6Lg5o1nZ8xMRAe7r3vbn7OwtiYAJGVlcXOnTuvmet87NgxACpVqkTbtm2vBnSTJk2IiIi44X2pwokTTijv3w9NmzqzOtatgz59ICXl2uVnznQuF9+0CV5//dqgjouD6tWdmSDmxlkYGxOgVJX9+/dfHXNevnw5SUlJAJQrV45WrVpdDefmzZtT2ou3obt40elVHzzoDHd07erMiZ43D373Oyesc8bCt986tyr96ivnznjZQZ0zsK1nXTALY2OCSEpKyjXhvG3bNlSVEiVK0Lx586vDGq1bt+YmH95l/9IlOHLECeuDB51ec4UKzlzpF15whjuy7+0Bzudq1WDyZJgzB2Jjr321amVhbWFsTBA7ffo0K1euvBrO69ev58qVK4SFhdG4cWM6duxI586dadu2LWXLlvVbXRkZkJzsBPXhw84TWMLDnYtdJkxwgvzMGWfZiAinJx4eDqNHw5Il1wZ1rVrODZrA6Y2H6mXlFsbGhJD09HS+/fZbli9fztdff83q1avJyMggMjKSVq1a0blzZzp16kSzZs2IjIx0tdZz55xQTktz5k4DvPYafP65037kCKSnO0MdBw86v+/Txxm3jo11etoxMc4c7OHDnd8nJzs99GB8aG1AhbGI9AV6AlWAN1X1KxFpBwwCIoAGqtq6oG1YGBvzk/Pnz7NixQoWL17MokWL2LRpE6pK+fLlueuuu+jUqROdO3emYcOGhZpK50+q8MMPcOqUM/0O4M03nQfRJic7M0OOHnXumLdihfP7Ro1g2zYnjKtVc14dO8Jzzzm///JL5xmJ1ao5M0wC6ZJzr4WxiEwGegFpqnpHjvZuwGtAOPCOqv69ENuqCLyiqkNztPUFolT1XwWta2FsTP5OnjzJ0qVLWbx4MYsXLyb7AqmoqCg6depE165d6dq1K1FB9MiTS5d+CtWPP3Zmhxw96oT1sWPOTJHXXnN+X7GiE/DZKlWCoUOdhwkAjBkDlStDdLQT1tHRzssfoe3NMG4P/AhMyQ5jEQkH9gD3AMnAOmAgTjCPzbWJh1U1zbPeP4APVXVjju3PBIap6tmC6rAwNqbwDh8+fLXXvGjRItLS0gBo2rQp3bt3p3v37rRo0YLwEDm7tnXrTyGd3bNu0cK578f5806POueJR4D//m8nrM+dg759fx7ULVs6TyjPjskb/QPDq8MUIhIHzMsRxq2AF1W1q+fzaABVzR3E2esL8HdgoaouytFeA/gfVR3+SzVYGBtzY7Kysti8eTPz589n/vz5rF69mqysLCpWrEiXLl3o3r073bp1C6pe8/W6cgVSU+H4ceeVkgK/+pUT2CkpztWLKSnO6+JFZ53XX4fHH3du+tSkiTOsciPnSgsTxkWZyh0DHMnxORloUcDyjwOdgQoiUkdV3/a0DwXezW8lERkBjACoUaNGEco1pvgKCwsjISGBhIQEnn32WU6fPs3ChQuZP38+X375JTNmzAAgISGBnj170rdvX5o0aRJwY81FERHhnBSMifn576pWhZUrnfeqTk85JeWnp7uULw9/+tONBXFhFaVn/Bugq6oO83x+AGiuqo/7qFbrGRvjA1lZWWzZsoX58+fzxRdfXO0116hRgz59+tC3b1/atWvn+gyNYFaYnnFRnkOQDMTm+FwdOFaE7eVLRHqLyIQz2ZMXjTFeExYWRpMmTXjmmWdYsWIFx48fZ/LkyTRu3JiJEyfSqVMnoqKieOihh/jkk09IT093u+SQVJSecQTOCbxOwFGcE3i/VdUdvinVesbG+Ft6ejpfffUVn376KXPnzuX06dOUKlWKnj17MmDAAHr27OnVS7VDlTdnU0wD7gZuBVKBF1R1koj0AP6JM4Nisqr+rahF57N/u5+xMS67fPkyK1asYPbs2Xz00UekpqZSrlw5+vTpw4ABA+jSpQslSpRwu8yAFFAXfXiD9YyNCQyZmZl88803TJ8+nVmzZnH69GkqVqxIv379GDhwIB06dLju24KGspAJY+sZGxO4MjIyWLRoEdOnT+fTTz/l3Llz1KxZk4ceeojBgwdTq1Ytt0t0XciEcTbrGRsT2C5cuMCcOXOYPHkyCxcuRFXp0KEDQ4YM4d57773m4a3Fia9nU/iNzaYwJjiULl2a/v37s2DBAg4ePMiYMWM4dOgQDz74INHR0Tz22GPs3LnT7TIDkvWMjTE+lZWVxfLly5k0aRIzZ87k0qVLdOzYkZEjR9K7d+8iPdUkWIRMz9gYE7zCwsK46667mDJlCkeOHGHs2LHs3buXfv36ER8fz0svvYT91WthbIzxo8qVKzNq1Cj279/PJ598Qt26dRk1ahQ1atRg9OjRHD9+3O0SXRMUYWxjxsaEloiICPr27cvixYvZsGED3bp146WXXiIuLo5HH32U/fv3u12i3wVFGKvqXFUdUaFCBbdLMcZ4WUJCAjNmzGD37t08+OCDTJo0ibp16/Lwww9z+PBht8vzm6AIY2NM6Ktbty4TJkzgwIED/OEPf2Dq1KnUrVuXJ5988uq9mENZUISxDVMYU3xUq1aNcePGsXfvXh544AHGjx9PfHw8Y8aM4WL2jYZDUFCEsQ1TGFP8xMbG8s4777Bz5066devG888/T4MGDZgzZw7BNCW3sIIijI0xxVe9evWYNWsWixcvpnTp0vTp04cePXqwZ88et0vzKgtjY0xQ6NixI5s3b+bVV19l1apVNGrUiJdffpkrV664XZpXWBgbY4JGZGQkf/zjH9m1axfdu3fnz3/+M61btw6JS6yDIoztBJ4xJqeqVasye/Zspk2bxv79+2natCkTJ04M6rHkoAhjO4FnjMlNRBgwYADbt2+nXbt2jBgxggEDBgTtpdVBEcbGGJOf6OhovvzyS8aOHcvHH39M8+bNg/LknoWxMSbohYWFMWrUKJYuXcqpU6do0aIFCxcudLus62JhbIwJGe3atWPdunXExsbSvXt3/vWvf7ldUqFZGBtjQkpcXBwrV66kW7duPPLII7zyyitul1QoFsbGmJBTvnx5PvnkE/r378/TTz/NX/7yl4CfaREUt9jP8UBSt0sxxgSJyMhIPvzwQ0qXLs2LL77IlStXGDNmjNtl5SsowlhV5wJzExMTh7tdizEmeISHhzNp0iQiIiL461//SkxMDI888ojbZeUpKMLYGGNuVFhYGG+99RYpKSk89thjxMTE0Lt3b7fL+hkbMzbGhLyIiAhmzJhBQkIC/fv3Z82aNW6X9DMWxsaYYqFs2bJ8/vnnREdH069fP06cOOF2SdewMDbGFBtVqlRh9uzZnDx5kkGDBpGZmel2SVdZGBtjipXGjRszfvx4Fi5cyGuvveZ2OVdZGBtjip1hw4bRu3dvnn32WXbt2uV2OYCFsTGmGBIRJkyYQJkyZRg2bFhAXBDitzAWkb4iMlFEPhORLp62GiIyR0Qmi8gof9VijDHR0dG8/PLLrFy5kqlTp7pdTuHC2BOWaSKyPVd7NxHZLSJJvxSmqvqpqg4HBgP9Pc23AZ+r6sNAg+sv3xhjbtyQIUNITEzk6aef5ty5c67WUtie8XtAt5wNIhIOvAl0xwnSgSLSQER+JSLzcr2q5Fj1Oc96AJuAASKyBFhalC9ijDHXKywsjPHjx5OSksL48eNdrUUKO1YiInHAPFW9w/O5FfCiqnb1fB4NoKpj81lfgL8DC1V1kaftT8BaVV0mIrNU9b481hsBjACoUaNG00OHDl3fNzTGmF/Qu3dvVq1axYEDB7jpppu8vn0R2aCqiQUtU5Qx4xjgSI7PyZ62/DwOdAbuE5Hsi8O/BP4gIm8DB/NaSVUnqGqiqiZWrly5COUaY0zeXnjhBU6dOsUbb7zhWg1FuTeF5NGWbzdbVV8HXs/Vth34WW/4Zzuyu7YZY3woMTGRHj16MG7cOJ566ilKlizp9xqK0jNOBmJzfK4OHCtaOXmzB5IaY3ztiSee4Pvvv+fjjz92Zf9FCeN1QF0RqSUiJYABwBzvlHUtEektIhOC9amvxpjA17lzZ2rXrs1bb73lyv4LO7VtGrAaqCciySIyVFWvACOBBcB3wExV3eGLIq1nbIzxtbCwMB555BFWrFjBtm3b/L//wiykqgNVtaqqRqpqdVWd5Gn/QlVvU9Xaqvo3XxVpPWNjjD8MGTKEiIgIPvjgA7/vOyguh7aesTHGHypVqkTHjh2ZPXu23y+RDoowNsYYf7n33ntJSkpi+/btv7ywFwVFGNswhTHGX/r06YOIMHv2bL/uNyjC2IYpjDH+EhUVRdu2bS2MjTHGbb169WLr1q2kpqb6bZ9BEcY2TGGM8ae2bdsCsHr1ar/tMyjC2IYpjDH+1LRpU0qUKMHKlSv9ts+gCGNjjPGnkiVLkpiYyKpVq/y2TwtjY4zJQ5s2bVi/fj0XL170y/6CIoxtzNgY429t2rQhIyODDRs2+GV/QRHGNmZsjPG3Vq1aAfht3DgowtgYY/ytSpUqxMbG+u1KPAtjY4zJR+3atdm/f79f9mVhbIwx+YiPj2ffvn1+2VdQhLGdwDPGuKF27docP36c8+fP+3xfQRHGdgLPGOOG+Ph4AL8MVQRFGBtjjBtq164NWBgbY4yrsnvG/hg3tjA2xph83HLLLVSoUMF6xsYY4yYR8duMCgtjY4wpQHx8vPWMs9nUNmOMW2rXrs2BAwfIysry6X6CIoxtapsxxi3x8fFkZGRw9OhRn+4nKMLYGGPcUrFiRQDOnTvn0/1YGBtjTACwMDbGmABgYWyMMQHAwtgYYwKAhbExxgQAv4WxiPQVkYki8pmIdPG0NRCRmSLylojc569ajDEm0BQqjEVksoikicj2XO3dRGS3iCSJyKiCtqGqn6rqcGAw0N/T3B0Yr6q/Bx68/vKNMSY0RBRyufeAN4Ap2Q0iEg68CdwDJAPrRGQOEA6MzbX+w6qa5nn/nGc9gPeBF0TkP4BKN/IFjDEmFBQqjFV1mYjE5WpuDiSp6n4AEZkO9FHVsUCv3NsQEQH+DsxX1Y2e7aYBj3mCfXZe+xaREcAIgBo1ahSmXGOMCTqF7RnnJQY4kuNzMtCigOUfBzoDFUSkjqq+7Qn4Z4CywP/PayVVnQBMABCREyJy6DpqvBX4/jqWDxTBWjcEb+1Wt38FXd0NGzaEG6+75i8tUJQwljzaNL+FVfV14PVcbQfx9HoLQ1UrF3ZZABFZr6qJ17NOIAjWuiF4a7e6/cvq/rmizKZIBmJzfK4OHCtaOcYYUzwVJYzXAXVFpJaIlAAGAHO8U5YxxhQvhZ3aNg1YDdQTkWQRGaqqV4CRwALgO2Cmqu7wXak3ZILbBdygYK0bgrd2q9u/rO5cRDXfYV5jjDF+YpdDG2NMALAwNsaYABCyYXw9l2r7qZ6DIrJNRDaLyHpP2y0islBE9np+Vsyx/GhP7btFpGuO9qae7SSJyOuei2m8XevPLn/3Zq0iUlJEZnja1+RxQZE3635RRI56jvtmEekRgHXHishSEflORHaIyBOe9oA+5gXUHdDHXERKichaEdniqfsvnnZ3j7eqhtwL55LsfUA8UALYAjRwuaaDwK252l4GRnnejwJe8rxv4Km5JFDL813CPb9bC7TCmec9H+jug1rbAwnAdl/UCjwKvO15PwCY4cO6XwT+lMeygVR3VSDB8748sMdTX0Af8wLqDuhj7tlHOc/7SGAN0NLt4+1aOPny5Tk4C3J8Hg2Mdrmmg/w8jHcDVXP8h707r3pxZqy08iyzK0f7QOBfPqo3jmtDzWu1Zi/jeR+Bc0WT+Kju/IIhoOrOVdtnOPd8CYpjnkfdQXPMgTLARpyrh1093qE6TJHXpdoxLtWSTYGvRGSDOPfbAIhS1RQAz88qnvb86o/xvM/d7g/erPXqOupMkTyDb28UNVJEtnqGMbL/9AzIuj1/zjbB6a0FzTHPVTcE+DEXkXAR2QykAQtV1fXjHaphfF2XavtJG1VNwLlt6GMi0r6AZfOrPxC/143U6s/v8RZQG2gMpAD/+IUaXKtbRMoBHwNPqurZghbNpw5Xas+j7oA/5qqaqaqNca4cbi4idxSwuF/qDtUwDrhLtVX1mOdnGvAJzl3vUkWkKoDnZ/ZtRvOrP9nzPne7P3iz1qvriEgEUAE45YuiVTXV8w8vC5iIc9wDrm4RicQJtA9VNfsOhgF/zPOqO1iOuafWH4CvgW64fLxDNYwD6lJtESkrIuWz3wNdgO2emh7yLPYQzpgbnvYBnjOytYC6wFrPn07nRKSl56ztgznW8TVv1ppzW/cBS9QzuOZt2f+4PH6Nc9wDqm7PfiYB36nqqzl+FdDHPL+6A/2Yi0hlEbnZ8740zt0kd+H28fbmAH4gvYAeOGd39wHPulxLPM7Z2C3Ajux6cMaQFgN7PT9vybHOs57ad5NjxgSQiPMf9z6cG/774gTSNJw/Ly/j/B9+qDdrBUoBHwFJOGej431Y9/vANmCr5x9I1QCsuy3On7Bbgc2eV49AP+YF1B3QxxxoBGzy1LcdeN7b/x5vpG67HNoYYwJAqA5TGGNMULEwNsaYAGBhbIwxAcDC2BhjAoCFsTHGBAALY2OMCQAWxsYYEwD+D3LkECylAXAaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Grab the GalaxyPopulation instance\n",
    "pop = mgb.pops[0] \n",
    "\n",
    "# Plot the numerical solution again\n",
    "pl.semilogy(E, flux[0] * E * erg_per_ev, color='k')\n",
    "\n",
    "# Compute cosmologically-limited solution\n",
    "e_nu = np.array([pop.Emissivity(10., energy) for energy in E])\n",
    "e_nu *= c / 4. / np.pi / pop.cosm.HubbleParameter(5.) \n",
    "e_nu *= (1. + 5.)**6. / -3.\n",
    "e_nu *= ((1. + 60.)**-3. - (1. + 5.)**-3.)\n",
    "e_nu *= ev_per_hz\n",
    "\n",
    "\n",
    "# Plot it\n",
    "pl.semilogy(E, e_nu, color='b', ls='--')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4229dc12",
   "metadata": {},
   "source": [
    "## Neutral Absorption by the Diffuse IGM\n",
    "\n",
    "The calculation above is basically identical to the optically-thin UV background calculations performed in the previous example, at least in the cases where we neglected any sawtooth effects. While there is no modification to the X-ray background due to resonant absorption in the Lyman series (of Hydrogen or Helium II), bound-free absorption by intergalactic hydrogen and helium atoms acts to harden the spectrum. By default, ARES will *not* include these effects.\n",
    "\n",
    "To \"turn on\" bound-free absorption in the IGM, modify the dictionary of parameters you've got already:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9e698887",
   "metadata": {},
   "outputs": [],
   "source": [
    "pars['tau_approx'] = 'neutral'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "477c4ba9",
   "metadata": {},
   "source": [
    "Now, initialize and run a new calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "e2a93922",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Loaded /Users/jordanmirocha/Dropbox/work/mods/ares/input/optical_depth/optical_depth_planck_TTTEEE_lowl_lowE_best_H_400x862_z_5-60_logE_2.3-4.5.hdf5.\n"
     ]
    }
   ],
   "source": [
    "mgb2 = ares.simulations.MetaGalacticBackground(**pars)\n",
    "mgb2.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e37ceac1",
   "metadata": {},
   "source": [
    "and plot the result on the previous results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "db3da79e",
   "metadata": {},
   "outputs": [],
   "source": [
    "z2, E2, flux2 = mgb2.get_history(flatten=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bdfbd472",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1876ce590>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD8CAYAAABJsn7AAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAoPUlEQVR4nO3deVxVdf7H8deHC6KhuYxbo4EoSDrgpCJpy4xbaZm5ZEWOu6lomWlWmi2mo+aY5o9+pqmDNv1Mh7RMHXNLKxtNxaU0ETRZQnOp3HMDvr8/LiYiKHqXw4HP8/G4j+79nnu+93O703sO3/M93yPGGJRSSnmfj9UFKKVUSaUBrJRSFtEAVkopi2gAK6WURTSAlVLKIhrASillEV+rC7gRlStXNrVq1bK6DKWUuiFbt2792RhTJW+7LQJYRNoD7UNCQkhISLC6HKWUuiEikpZfuy2GIIwxS40x/cuXL291KUop5Ta2CGCllCqONICVUsoiGsBKKWURDWCllLKIBrBSSlnEFgEsIu1FZOaJEydueN+tW7eyZcsWD1SllFKusUUAuzINbfTo0URFRdGqVStWr16Nrn+slCoqbBHArpg3bx6TJk0iMTGRBx54gMjISD766COysrKsLk0pVcIV+wC+9dZbGT58OCkpKcyePZvTp0/z+OOPc8cddzBz5kzOnTtndYlKqRKq2AfwJf7+/vTt25fdu3ezcOFCKlSowIABAwgODuYf//gHNzO+rJRSrigxAXyJw+Hg0UcfZfPmzaxZs4aIiAheeuklAgMDGTlyJIcOHbK6RKVUCVHiAvgSEaFVq1asWrWKhIQE2rRpw8SJE6lVqxYDBw7khx9+sLpEpVQxV2IDOLfGjRsTHx9PUlISPXr0IC4ujrp16xIdHc327dutLk8pVUx5LYBFpKOIzBKRT0XkgVztASKyVUQe9lYtBQkNDWXmzJmkpqYyfPhwli9fTqNGjWjbti1ffPGFTmFTSrlVoQJYROJE5IiI7MrT3lZEkkRkn4iMuFYfxpjFxph+QC/giVybXgLib7Buj7rtttuYOHEi6enpjB8/nu3bt9OiRQuaNWvGJ598QnZ2ttUlKqWKgcIeAc8F2uZuEBEHMA14EKgPPCki9UUkQkSW5XlUzbXrKzn7ISKtgd3AYRe/h0dUqFCBkSNHkpqayvTp0zl69CidO3fmT3/6E3PmzOHChQtWl6iUsrFCBbAx5ivg1zzNUcA+Y8x+Y8wFYAHQwRiz0xjzcJ7HEXGaCHxmjNmW00cLoCnQFegnIkVyTLpMmTLExMSQlJTE/PnzKV26NH369KF27dpMmTKFU6dOWV2iUsqGXAm8GsCPuV5n5LQVZDDQGugiIjEAxphRxpjngA+BWcaYq/62F5H+IpIgIglHjx51oVzX+fr6Eh0dzbZt21ixYgWhoaE8//zzBAUF8dprr2F1fUope3ElgCWftgLPUhljYo0xjY0xMcaYGXm2zTXGLCtgv5nGmEhjTGSVKlfd084SIkKbNm1Yt24dGzdu5K9//Stjx44lKCiIwYMHk5qaanWJSikbcCWAM4Dbc72uCRx0rZz8ubIamqc1bdqUTz75hN27dxMdHc2MGTMICQmhe/fu7Ny50+rylFJFmCsBvAUIFZFgESkFRANL3FPWlexwU8569eoRFxfH/v37efbZZ/nkk09o0KAB7du35+uvv7a6PKVUEVTYaWjzgY1AmIhkiEhfY0wm8AywEkgE4o0x33uiyKJ8BJzX7bffzpQpU0hPT2fMmDFs3LiR++67j3vvvZdly5bpFDal1O/EThcXREZGmoSEBKvLuCFnzpwhLi6Ot956i/T0dMLDw3nxxReJjo7Gz8/P6vKUUl4gIluNMZF524vktK+87HQEnFdAQACDBw9m3759fPDBBxhj6NGjByEhIbzzzjv89ttvVpeolLKILQLYDmPA1+Pn50e3bt347rvvWLp0KbfffjvPPvssQUFBjBkzhl9/zTvNWilV3NkigIsTHx8fHn74Yb7++mvWr19P06ZNef311wkMDGTYsGFkZGRYXaJSyktsEcB2HoK4lnvvvZelS5fy3Xff0alTJ2JjY6lduzZ9+vQhMTHR6vKUUh5miwAuDkMQ1xIREcEHH3zAvn37iImJYcGCBdSvX59OnTqxadMmq8tTSnmILQK4pKhVqxaxsbGkpaXx6quv8uWXX9K0aVNatGjBihUrdDlMpYoZWwRwcR2CKEiVKlUYM2YM6enpTJkyhb179/Lggw/SsGFDFixYQGZmptUlKqXcwBYBXNyHIApStmxZhg4dyv79+4mLi+PcuXM8+eSThIWFMWPGDM6ePWt1iUopF9gigEu6UqVK0bt3b3bv3s3HH39M5cqVGThwILVq1WLChAkcP37c6hKVUjdBA9hGfHx86NSpE9988w3r1q2jYcOGvPzyywQGBvLSSy/x008/WV2iUuoG2CKAS9oY8PWICM2bN2fFihVs27aNdu3a8dZbb1GrVi369+/P3r17rS5RKVUItgjgkjoGXBgNGzZk/vz5JCcn06dPH/71r38RFhbG448/ztatW60uTyl1DbYIYHV9derUYfr06aSmpjJixAhWrlxJZGQk999/P59//rlOYVOqCNIALmaqV6/O+PHjSU9PZ+LEiezatYvWrVsTFRXFokWLyMrKsrpEpVQOWwSwjgHfuPLly/Piiy+SkpLCe++9x/Hjx+nSpQv16tVj9uzZnD9/3uoSlSrxbBHAOgZ880qXLk3//v3Zs2cP8fHxlCtXjn79+hEcHMxbb73FyZMnrS5RqRLLFgGsXOdwOHjsscdISEhg1apV1K9fnxdeeIGgoCBGjRrFkSNHrC5RqRJHA7iEERHuv/9+1qxZw+bNm2nVqhUTJkwgKCiIp59+mpSUFKtLVKrE0AAuwZo0acLChQtJTEykW7duzJo1i9DQULp27cq3335rdXlKFXsawIqwsDBmzZpFSkoKQ4cOZenSpdx555089NBDfPXVVzqFTSkP0QBWv6tRowaTJk0iPT2dv//97yQkJPDXv/6Vu+++m08//VTv6KyUm9kigHUamnucOAFHj8Jvv8G1DmorVqzIqFGjSEtLY9q0aRw6dIiOHTsSHh7O+++/z4ULF7xXtFLFmC0CWKehFc78+fDcc/DQQ5CQcPX2Ll2galUICICPPrp6+8cfw7x58O23cP48lClThkGDBrF3717mzZuHn58fvXr1IiQkhKlTp3L69GmPfyelijNbBLC60rFjzqPYvN5/H/7nf+CzzyA19ertvr6Xn5crd/X2adOgWze480749NPc+/nStWtXduzYwfLlywkODmbo0KEEBQUxevRofv75Z1e/klIlkgawjcyeDffeC5UrXxmQl4SFXX6e37TeSpXgD38Af//8A/jHHy8/Dw+/enuHDkJc3IP06PElK1d+w7333ssbb7xBUFAQQ4YMIT09/ca/lFIlmTHGNo/GjRubkmzkSGOco7fGdO9+9fY1a4wZN86YhQuNOXDg2n1lZ1/dNmaMMR06GBMWZsyFC1duO3fOmFKlLn9+aqqz/fvvvzc9e/Y0vr6+xtfX1/To0cPs2rXrpr6fUsUVkGDyyTQxNppiFBkZaRLyG9wsRrKyYMUKOH4c/va3K7ft2gURESAC7drB0qXeq+ubb6BZM+fzsDDYs+fK7Xv3pvP44/tJTJzI+fOrad/+IUaMGMHdd9/tvSKVKqJEZKsxJjJvu29+b1bW2LsXWreG9HSoVg0efxz8/C5vDw93nii77z7nMIQ33XUX7NwJK1c6T+LllZwcyI4dgUBzatTI4L///TP33HMP9913HyNGjODBBx9ERLxbtFJFnI4BFyG1asGlGx4fPgyLF1/9nk6dvB++4DzqDg+H55+HmJirty9adPl51641SU9PZ+rUqaSmptKuXTsaNLiTDz/8UO/orFQuGsAWMAZWrYKMjCvb/fzgqaecJ8peeAEir/qDpejq2xeefhqqV4dHH4WAgACGDBnCvn37mDt3LgcPPsff/laT6tVfZsqU9/gtv2kcSpUwXhsDFpGOQDugKjDNGLNKRJoDY4HvgQXGmC+u1UdxGANetQqGD3f+OT9sGEyefOX2U6ecQVy6tDX1uSo723m0nHu0ITMTatY0HD58qbEDVapsZMiQIQwaNIiKFStaUqtS3lLQGHChjoBFJE5EjojIrjztbUUkSUT2iciIa/VhjFlsjOkH9AKeuNQMnAZKAxkF7FqsXLzoDF+AWbOcV6flVq6cfcMXwMfnyvAF2LEDfv7Z2XjbbYa1a4fTpEkTXnnlFQIDAxk+fDj79h3wfrFKWS2/qRF5H8BfgEbArlxtDuAHoDZQCvgWqA9EAMvyPKrm2m8y0CjnuU/OP6sB865Xh92moZ0/f3VbVpYx9eoZc8stxgwebMzRo96vywqHDhkzaZIx06ZdbtuxY4fp2rWrEfmTgRPmjjvWmSVL9ltXpFIeQgHT0Ao9BxeolSeAmwErc70eCYy8xv4CTARa57OtFLCwgP36AwlAQmBgoKf/PbnFiRPGjB9vTNWqxuQ3JXb7dmN++cXrZRVZvXsf/31+MXxqOnfubDZv3mx1WUq5TUEB7MpJuBpArmunyMhpK8hgoDXQRURiAESks4i8B3wA/G9+OxljZhpjIo0xkVWqVHGhXO+JiYGXX3ZejTZu3NXb77zTeVWackZucvLlNT6efPIX1q5dS1RUFK1atWL16tVkZdlnrrpSN8KVAM5vUmeB/6UYY2KNMY2NMTHGmBk5bR8bYwYYY54w1zgBZ7fV0IYOvfx806b8121QTiKwfj18+SU8+yz83//1Ji0tjUmTJrFnzx4eeOABKlTYSvPmKWzZond0VsWLKwGcAdye63VN4KBr5eTPFNHV0A4fhtjYq5d2bNIEBgyAf/7TecXYLbdYU59diMBf/uJcSMjHB2699VaGDx/O/v37GTVqMadPR/Lll8FERWXx9tvvc+7cOatLVsotXAngLUCoiASLSCkgGljinrKuVBSPgEeMgOBgGDIE1q69evuMGdCnz5VXsqkb4+/vj8PR4ffXf/jDaoYN60VwcDATJ06kKP3vQambUdhpaPOBjUCYiGSISF9jTCbwDLASSATijTHfe6LIongEfPw4nD3rfD52rKWlFGtvvOFch+Jvf4Plyx/i888/JyIighEjRhAYGEj79v9mwIAzpKVZXalSN84Wi/GISHugfUhISL+9e/d6/fMvXrz6SDYtDUJCoEEDePVV6NDh6vmvynO2bt3KxImT+Oij14F6iGQxffphBgz4o9WlKXUVly7EsJpVR8Cpqc6x3D//+fIaDZcEBcG2bc47T3TsqOHrbY0bN6Z//wVAPQCMOcPAgfWJjo5m+/bt1hanVCHZIoCtcPEiNG0KM2dCYiJ8+OHV77m0NKSyRsuW8J//QKtWMGCADy+8MIDly5fTqFEj2rRpw6efrmfqVMOpU1ZXqlT+bBHAVpyE8/ODZ565/HrNGq99tCokHx/n/e/WrIFp08oyceJE0tPTmTBhAt9++y0dOy5n6FChevULvPOO3tFZFT22GAO+xFOL8WzfDsnJ8MQTV7afOOE8+TN0qPNoS4927ePkyXMEBmZx4oRz8eLq1UcwblxdunXrRqlSpSyuTpU0th4D9pSjR6FtW2jUyDnWe/LkldvLl4dly5x/4mr42ou/f2kmTgygbl1DpUq/UbXq5/Tt25fatWszZcoUTp06RVLS1XO4lfImWwSwp4YgKlWClBTn8xMn4L333Nq9spC/v/P/VBMThc2bb2HHjs2sWLGC0NBQnn/+eWrWvId69bJo0uQiH3+sQaysYYsAdscsiGPHrl4A3eGAF190jiVGR0ObNi4WqoocHx+oUwdEhDZt2rBu3Tq++eYbqlQZjzEOtm7149lnk0lLS7W6VFUC2SKAXZGaCgMHQs2azqvX8urWzTn+O3++c06vKv7uuusuHn74Yfz8nCfmDh0aRkhICN27d2dnzmLNekSsvKHYn4TbuvXyrX38/Jw3vKxe3QPFKdv56Sfn9MLHH89g6tS3ee+99zhz5gzt2rXjxIn3adToDwwd6rxXn1KuKOgknC0C2NUr4e65BzZscB7hzpnjPOmmVF6//vor06ZNY8qULzh+/HMA/PyyOHBAqFKl2P+xqDzI1rMgXB0DnjjRudzhjh0avqpglSpV4tVXX6VXrxW/t128GE+LFg344IMPuHjxooXVqeLIFgHsqnvvdS53qFPJVGFMmeLHihXQsmU2Y8eWR0To0aMHISEhxMbGsnjxOebPv/rydKVulC2GIC4pDndFVvZjjGH58uW8+eabfP311zgcO8nKCicoKIuPPnLQpInVFaqiztZDEEpZSURo164d69evZ9KkXWRlhQOQlnaOuLhX+PHHH6/Tg1L5s0UAF8UF2VXJ1KvXn3j9dShfPpM77ljPrFlvUrt2bXr37k1iYiJnzsBBj9wXRhVHOgSh1E04cwYuXICTJ9OYPHkys2fP5uzZs4SHx5GU1JOePX148UUIDbW6UlUU6BCEUm4UEAAVK0JQUBCxsbGkpaUxatQbfP99Gy5e9GH2bHjnne+x0wGO8j4NYKXcoEqVKvTv/xqRkdUA8PE5wjvvNKZhw4bMnz+fTJ0yofKhAayUmwQGwqZNDr78Ej78sCJz5szg/PnzdO3albp16zJhwvs0bZrFRx9BVpbV1aqiQMeAlfKg7OxslixZwoQJE9i8+VHgRQA6d77AokW6LnFJoWPASlnAx8eHjh078tVX31Cx4nO/t3/2WXdefPFFDuqUiRLNFgGs09CU3fn7C3v2lGLUKLjzztM88oiDyZMnExwcTP/+/dm7dy/vvw+HD1tdqfImHYJQysuMcV4W/8MPPzB58mTi4uI4fz4U2EmpUtn06+dDbKxzLWNVPOgQhFJFxKU1SerUqcO7775LWloa4eHvA3Dhgg/x8V+zdu0ancJWAmgAK2WxatWqMXZsIxo1ck5Vy8oaz/3330+TJk1YuHAhWVlZHDumi8QXRxrAShUBHTtCQoIvW7bAgQMfM3PmTE6cOMFjjz3GHXfUo1GjIzRtms3ixZCdbXW1yl00gJUqIkScd28pXbo0/fr1Y8+ePcTHxyPSktTUqmze7MNjj2WSknLK6lKVm2gAK1VEORwOHnvsMQYPnv77/esyM2fSuPHtjBo1isM6ZcL2dBaEUjbw008wdSo0a/Yd8+aNZdGiRfj7+9OnTx9CQ1/j4sVqxMRAuXJWV6ryY/k94USkI9AOqApMM8asEhEfYCxwK5BgjHn/Wn1oACvllJyczKRJk5g7dx6Zmd8Dwdx6ayYrVvjSrJnV1am8XJqGJiJxInJERHblaW8rIkkisk9E8rnp+2XGmMXGmH5AL+CJnOYOQA3gIpBRmFqUUlC3bl1mzZrF1KkZQDAAJ08eZ/TozmzYsMHa4lShFXYMeC7QNneDiDiAacCDQH3gSRGpLyIRIrIsz6Nqrl1fydkPIAzYaIwZBgx05YsoVRI99VQlZs2C4OAsWrXaxdatX3HPPffQokUL1qxZwy+/GPbvt7pKVZBCBbAx5ivg1zzNUcA+Y8x+Y8wFYAHQwRiz0xjzcJ7HEXGaCHxmjNmW00cGcCznua4PpdQN8veHp56C5GQHy5Y1Jy0tjSlTppCcnMz9999PgwZzCA01dOtm2LPH6mpVXq7MgqgB5L4ZVkZOW0EGA62BLiISk9P2MdBGRN4BvspvJxHpLyIJIpJw9OhRF8pVqvjy9YXSpSEgIIChQ4eyf/9+3n77nxw69CjZ2cK8ecL06RvI0nUwixRXAji/m7wXeEbPGBNrjGlsjIkxxszIafvNGNPXGDPYGDOtgP1mAm8A20qV0uX7lCoMf39/OnToQ4sWzmkRfn5pxMbeR/369ZkzZw4XL160uEIFrgVwBnB7rtc1AY+srWeMWWqM6V++fHlPdK9UsRQcDGvW+LBpEyxdWpOPPvo3t9xyS87UtVBGj57H3Xdns3SpXuZsFVcCeAsQKiLBIlIKiAaWuKesK+lylErdvKgoaNPGQZcuXdi2bRv/+c9/+OMf/8gbb5xl40YfHnkEBgy4YHWZJVJhp6HNBzYCYSKSISJ9jTGZwDPASiARiDfGfO+JIvUIWCn3EBEeeughPvvsv5Qp0+v39gULOjFmzBiOHTtW8M7K7fRKOKVKqIMHYcoU+O9/j1OlSg+WLl1KuXLlGDRoEMOGDWPLlqo88AD4+Vldqf1ZfiWcK0SkPdA+JCSk3969e60uR6li5dIC8d999x3jx48nPj4eP7/mXLiwlsDATEaP9qV3b6urtDdbB/AlegSslOclJSXRsuUFDh6MAOCOOzaxZk1NatS41ixTdS22viOGnoRTyntCQ8N4+ukIKlTIQiSbvXt7UadOHZ599lkOHDhgdXnFii0CWE/CKeU9Pj7w8suQnu5g8WIf9u5dTvfu3Zk+fTp16tRh8OAhNG9+jnfegXPnrK7W3nQIQilVKCkpKYwfP564uNNkZ88HoE6diyQl+eFwWFxcEadDEEoplwQHBzNr1iyaNYv7vS01dTJDh+rQxM2yRQDrEIRSRcfq1WWIjYXatS/w5JOHcg1NDObAgQMkJ4MuOVE4OgShlLopl6avXRqamDt3Lg5HOXx9U6lR4xbGjfOlc2fnmHJJZ+shCKVU0SM5y3FdGppISkriT396lzNnbiU52Zc+fU7yyy+nrS2yiLNFAOsYsFJFX+3atenePZqAAOf4w6lTLxAeXofY2FjOnz9vcXVFky0CWMeAlbKH556DtDQHf/87fP11X8LDwxkyZAh169Zlzpw5fPhhFomJVldZdOgYsFLKo9asWcPLL7/Mli0ZiOxDpDRPPSVMniyULWt1dd6hY8BKKUu0bt2aTZs20br1Boy5hexsH+bN+4HExG3X37mY0wBWSnmciPCPf9SiZctsAByOkdx1VyS9evUq0XOIbRHAehJOKftr2NB5h46NGyE9fTYvvPAC8+fPp27duowePZo33zzPpk1WV+ldOgaslLJMSkoKI0aMID4+GdgK+NC7t2H6dMHf3+rq3EfHgJVSRU5wcDD//ve/iYxcx6U4Wrx4M8nJu6wtzEs0gJVSlvv3vyvQvr3B4cgiO3swjRo15IUXXuD06eJ9IYcGsFLKcrVrw5IlQmKig337ltOzZ0/eeust6tWrx6JFi5g713DmjNVVup8GsFKqyAgNhcqVKzN79mw2bNhA5cqV6dLlfXr3FurVy2L1aqsrdC8NYKVUkdSsWTPWrt3CrbfOA+DHHx2MHfsDdpo4cD22CGCdhqZUyVShgi+xseUoXz4LP79fWb++CR06dODgwYNWl+YWtghgXQtCqZJJBHr2hKQkB2vWlGfy5FdYvXo1ERERLF68GGMgM9PqKm+eLQJYKVWyVasGf/mLg2HDhvHtt99Sq1YtOnXqRMuW82jWLJu0NKsrvDkawEopW6lbty4bNmygT58JfPFFRxISfPjznzOx4zVaGsBKKdvx9/enVasROBxlADh5cj979iyxuKobpwGslLKlrl3hq698qF37IuHhr9O9ewdee+01srOzrS6t0DSAlVK2dffdkJzsx5Ytc+nTpw9jx46lY8eOnD59GjvMVtMAVkrZmsPhHJKYPXs206ZNY/ny5dx1V18aNLhIUpLV1V2br7c+SEQ6Au2AqsA0Y8wqEbkP+FtOHfWNMXd7qx6lVPEiIgwaNIiAgLr07n0HxvjRrFkWa9c6uPNOq6vLX6GOgEUkTkSOiMiuPO1tRSRJRPaJyIhr9WGMWWyM6Qf0Ap7IaVtvjIkBlgHv39Q3UEqpXGrWbI2//20AHD9+jp0791tcUcEKOwQxF2ibu0FEHMA04EGgPvCkiNQXkQgRWZbnUTXXrq/k7JdbV2D+TX0DpZTKpVUr+PJLB9WqXaRixe4MH96M3bt3W11WvgoVwMaYr4Bf8zRHAfuMMfuNMReABUAHY8xOY8zDeR5HxGki8Jkx5vebQYlIIHDCGHPSTd9JKVXCRUVBaqofGzZMwOFw0LJlSxKL4O2YXTkJVwP4MdfrjJy2ggwGWgNdRCQmV3tfYE5BO4lIfxFJEJGEo0ePulCuUqokKV0awsLCWLduHSJCmzZt2L//AIcPW13ZZa4EsOTTVuDED2NMrDGmsTEmxhgzI1f768aYDdfYb6YxJtIYE1mlShUXylVKlURhYWF89tln/PprFg0aHKBlyyxOnbK6KidXAjgDuD3X65qAR5Yo0tXQlFKuiIi4k5o1d3PmTBS7dzvo1csUiXnCrgTwFiBURIJFpBQQDXjkWkBdDU0p5QqHA0aOvJwfP/30NUXhgrnCTkObD2wEwkQkQ0T6GmMygWeAlUAiEG+M+d4TReoRsFLKVT17wksvGVq0mMbGjX9h5crlVpekt6VXSpUsZ8+epVmzZvz444/s2rWL2267zeOfaevb0usRsFLKXcqUKUN8fDy//fYbMTExlt7iyBYBrGPASil3qlu3LuPGjWPJkmVER+/g88+tqcMWAayUUu7WseMQypbdQXx8QwYNyuLiRe/XYIsA1iEIpZS7lSnjIDu7PgDJyQ7++U/v12CLANYhCKWUu912G4we7cDhuICPz2iaNvX+2pW2CGCllPKEIUNg8+aTBARMYdy4V7z++bYIYB2CUEp5QqlS0KhRZQYPHsyiRYu8vmCPLQJYhyCUUp40dOhQypQpw4QJE7z6ubYIYKWU8qTKlSszcOBA5s37hClTfiYz0zufqwGslFKAw/EK2dlpPP98ZT791DufaYsA1jFgpZSn+ftXACoBMG2ad1bqsUUA6xiwUsrT+vcHEQPso2rVZK8sV2mLAFZKKU+rWRO+/TaLypXvBkYj+d1yws00gJVSKkdEhC+PPtqZZcuWcfbsWY9/ngawUkrl0qVLF86cOcOKFSs8/lm2CGA9CaeU8pbmzZtTuXJlFi1a5PHPskUA60k4pZS3ZGb6Eh7+Eh9/3IKBAz17Js4WAayUUt6SkQFffDGcs2f7MnduNhcueO6zNICVUiqXkBAIDHReCnfunINt2zz3WRrASimVR0yML5UqvUdU1BjCwz33ORrASimVx8iR0L79RtLS3iUgwHPjwLYIYJ0FoZTytqioKA4fPkx6errHPsMWAayzIJRS3hYVFQXA5s2bPfYZtghgpZTytgYNGuDv7+/RAPb1WM9KKWVjn39eijJlljB9ehg1asBzz7n/MzSAlVIqHz/8AMePPwDA9u2e+QwdglBKqXxERFx+vmNHlkc+QwNYKaXy0bAhPP30duAupk71zNUYGsBKKZWPW2+FgQNLAZs5dGifRz5DA1gppQpQq1YtAFJSUjzSv9cCWEQ6isgsEflURB7IaQsUkSUiEiciI7xVi1JKFUZAQABVq1a1NoBzAvKIiOzK095WRJJEZN/1AtQYs9gY0w/oBTyR01wX+I8xpg9Q/8bLV0opzwoODiYlJYUsD5yHK+wR8Fygbe4GEXEA04AHcYbnkyJSX0QiRGRZnkfVXLu+krMfwHYgWkTWAutc+SJKKeVuq1ZBcvI81q2LZ8AA9/dfqHnAxpivRKRWnuYoYJ8xZj+AiCwAOhhjJgAP5+1DRAR4E/jMGHPplGJv4PWc/hcCc/LZrz/QHyAwMLBQX0oppdzh7Fk4dqwOAIcOub9/V8aAawA/5nqdkdNWkMFAa6CLiMTktK0AnhWRGUBqfjsZY2YaYyKNMZFVqlRxoVyllLox1apdfn7woPvHIFwJ4Pxu2lzgum3GmFhjTGNjTIwxZkZO2y5jTJectuEFfpCuhqaUskBEBLz22lKgJvPmpbq9f1cCOAO4PdfrmsBB18rJn66GppSyQkAANGvmBxzg2LHDbu/flQDeAoSKSLCIlAKigSXuKetKegSslLJK9erVATjkgUHgwk5Dmw9sBMJEJENE+hpjMoFngJVAIhBvjPne7RWiR8BKKetUyxkIPnzY/UfAhZ0F8WQB7cuB5W6tKB8i0h5oHxIS4umPUkqpK1w6+X/o0BGyssDhcF/ftrgUWY+AlVJWGTTIF/iVsWNfZYmbB1l1PWCllLqG8+cBKmIMHDvm3r5tcQSsJ+GUUlapWPHy81On3Nu3GOO5Wy67W2RkpElISLC6DKVUCfLrr9C+fQfgBP/97xc31YeIbDXGROZtt8URsFJKWaVSJahSRTh58he3922LANYhCKWUlcqVK8eZM2fc3q8tAlhnQSilrFS2bNmSG8BKKWWlgIAATp06z+nT7u1XA1gppa5h/Xp4++2JnD17nIcfdu+kBVsEsI4BK6WsUqYMZGc7L387eTLbrX3bIoB1DFgpZZWyZS89yyQ7uwQGsFJKWaVuXZg9+/8APxYtSndr3xrASil1DT4+ULasHwDnndclu69vt/amlFLFkL+/P1BCA1hPwimlrHQpgC9cuODWfm0RwHoSTillJWPKABU4eDAbdy6fo8tRKqXUdTzyyF+AY3Tu7FyeslQp9/RriyNgpZSykp/f5cNedw4DawArpdR1lCuXDRynfPlzZGa6r18NYKWUuo716/cDFXn33Y+vWKDdVRrASil1HToNTaehKaUsUqIDWKehKaWsVKIDWCmlrPTzz37A7Rw9WgZ3XouhAayUUtfxxBO3AOmMG9efnTvd168GsFJKWUQDWCmlrqNqVQOkU6HCSfz83NevBrBSSl3HwoVngSBeeWUWDRq4r18NYKWUsojXAlhEOorILBH5VEQeyGmrLyLxIjJdRLp4qxallCoKChXAIhInIkdEZFee9rYikiQi+0RkxLX6MMYsNsb0A3oBT+Q0Pwi8Y4wZCPS48fKVUsq+Crsc5Vzgf4F/XWoQEQcwDbgfyAC2iMgSwAFMyLN/H2PMkZznr+TsB/AB8LqIPAL84Wa+gFJKeVpKigD1OHSoEmfPOu+U7A6FCmBjzFciUitPcxSwzxizH0BEFgAdjDETgIfz9iEiArwJfGaM2ZbT7xHg6Zww//imv4VSSnlQz55lgN289RZER0Pjxu7p15UF2WsAP+Z6nQHcdY33DwZaA+VFJMQYMyMn1F8GAoBJ+e0kIv2B/jkvT4tI0k3UWh7wxEISrvZ7o/sX9v2Fed+13nOj2yoDPxeiLm/z1O/uat83s6+7fvub3V5Qe0n77ctHRt5Uv0H5thpjCvUAagG7cr1+DJid63V3nOO5he7TWw9gZlHs90b3L+z7C/O+a73nRrcBCVb/xt783V3t+2b2dddvf7Pbr9Feon57d/fryiyIDOD2XK9rAgdd6M+TlhbRfm90/8K+vzDvu9Z7bnZbUePJWl3p+2b2dddvf7Pb7fS7Q9H9b/4KkpPq13+jc7hgmTEmPOe1L5AMtAIOAFuArsaY791ZoCr6RCTBGBNpdR3K+/S3d01hp6HNBzYCYSKSISJ9jTGZwDPASiARiNfwLbFmWl2Asoz+9i4o9BGwUkop99JLkZVSyiIawEopZRENYKWUsogGsHI7EaknIjNEZKGIDLS6HuU9IhIgIltF5KqrYdXVNIBVodzIgkzGmERjTAzwOKBTlGzsJhbiegmI926V9qUBrAprLtA2d0OuBZkeBOoDT4pI/ZxtjwBfA597t0zlZnMp5O8uIq2B3cBhbxdpV66sBaFKEHMDCzIBu40xS4AlIvIf4EOvFqvc5gZ/97I413WpD5wVkeXGmGxv1ms3GsDKFfkuyCQizYHOgD+w3PtlKQ/L93c3xjwDICK9gJ81fK9PA1i5QvJpM8aYL4AvvFuK8qJ8f/ffnxgz13ul2JuOAStX2GlBJuU++ru7iQawcsUWIFREgkWkFBANLLG4JuV5+ru7iQawKhRdkKlk0t/ds3QxHqWUsogeASullEU0gJVSyiIawEopZRENYKWUsogGsFJKWUQDWCmlLKIBrJRSFtEAVkopi2gAK6WURf4f0ZqYnP0aYDgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the numerical solution again\n",
    "pl.semilogy(E, flux[0] * E * erg_per_ev, color='k')\n",
    "pl.loglog(E2, flux2[0] * E2 * erg_per_ev, color='b', ls=':', lw=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7b667ee",
   "metadata": {},
   "source": [
    "The behavior at low photon energies ($h\\nu \\lesssim 0.3 \\ \\mathrm{keV}$) is an artifact that arises due to poor redshift resolution. This is a trade made for speed in solving the cosmological radiative transfer equation, discussed in detail in Section 3 of [Mirocha (2014)](http://adsabs.harvard.edu/abs/2014arXiv1406.4120M>). For more accurate calculations, you must enhance the redshift sampling using the ``tau_redshift_bins`` parameter, e.g.,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e5e62ac0",
   "metadata": {},
   "outputs": [],
   "source": [
    "pars['tau_redshift_bins'] = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b77898d9",
   "metadata": {},
   "source": [
    "The default optical depth lookup tables that ship with *ARES* use ``tau_redshift_bins=400``, but tables with ``tau_redshift_bins=1000`` are also provided. You can run ``$ARES/input/optical_depth/generate_optical_depth_tables.py`` if you'd like to generate tables with a different resolution."
   ]
  }
 ],
 "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
}