{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cell-based cluster model \n",
    "To illustrate an example from a cell model, we will use Model B from [Reynolds et al. 2020](https://ui.adsabs.harvard.edu/abs/2020ApJ...890...59R/abstract), _\"Astrophysical Limits on Very Light Axion-like Particles from Chandra Grating Spectroscopy of NGC 1275\"_. It is an approximate model for a turbulent field in the Perseus cluster. \n",
    "\n",
    "The total magnetic field strength at a distance $z$ along the line of sight is set according to\n",
    "\n",
    "$$\n",
    "B(z) = 7.5\\,\\mu{\\rm G} \\left[ \\frac{n_e(z)}{n_e(25\\,{\\rm kpc})} \\right]^{0.5}\n",
    "$$\n",
    "\n",
    "where the electron density $n_e$ is set by the analytic profile given by Churazov et al: \n",
    "\n",
    "$$\n",
    "n_e(z) = \\frac{3.9\\times10^{-2}}{\\left[1+ (z/80\\, {\\rm kpc})^2\\right]^{1.8}} + \\frac{4.05\\times10^{-3}}{\\left[1+ (z/280\\, {\\rm kpc})^2\\right]^{0.87}}\n",
    "~{\\rm cm}^{-3}.\n",
    "$$\n",
    "\n",
    "The cell sizes $\\Delta z$ are set according to a probability distribution function $p(\\Delta z)$, given by a power-law so that \n",
    "\n",
    "$$\n",
    "p(\\Delta z) \\propto \\Delta z^{-1.2}\n",
    "$$\n",
    "\n",
    "spanning $3.5-10$\\, kpc. In their model B, Reynolds et al. 2020 scale these minimum and maximum cell sizes linearly with radius as $z/50{\\rm kpc}$, because coherence lengths are expected to grow with distance from the cluster centre. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt \n",
    "import numpy as np \n",
    "import alpro \n",
    "alpro.util.set_default_plot_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we initialise the survival class and set up the ALP parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = alpro.Survival(\"1275b\")\n",
    "s.init_model()\n",
    "s.set_params(1e-12 * 1e-9, 1e-13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now initialise a random cell model with a given random number seed and plot $B_\\perp$ as a function of line of sight distance $z$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$z$ (kpc)')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYFUlEQVR4nO3de7hddXng8e/rGYMabnUEuShgCvEGNWCrIBiONjwwMj5atTKMgHRqOjaKjRcQ2hEZrJV5RjHxQu1EZ1CkgzNt9RG5aKOEFEREhgPhGjHgtFzkMphIuETiO3+sdcLOZu9z9jln7b32Oev7eZ71ZO912+/5cdjv+V3W7xeZiSRJz6o7AEnScDAhSJIAE4IkqWRCkCQBJgRJUsmEIEkCTAiSpNK/qvPDI2IEOAs4AdgTuA+4EDgrM5+a4LobgN2AOwcQpiTNFfsDD2bmwZ0O1poQgI8C7wPeDawDfgf4KvAk8IkJrtttl1122XvRokV79z9ESZobxsbG2LhxY9fjdSeE1wEXZ+bF5fu7I+LbwGsnue7ORYsW7b1mzZq+BidJc8no6ChXXnll15aVuvsQrgLeEBEvA4iIVwBvBC7tdHJErImINcCigUUoSQ1Rdw3hvwA7AbdGxNYynk9m5nn1hiVJzVN3DeE44CTg3wOHlK+XRcQfdzo5M0czcxQYG1iEktQQddcQ/ivw6cy8qHy/LiL2Bc4AvlJfWJLUPHXXEJ4HbG3bt5X645Kkxqm7hnAxcHpE3AXcAhwMfAj4Wq1RSVID1Z0QTqF43uA8YHeKB9NWAWfXGZQkNVGtCSEzfwUsL7eBWrV2AytWr2fzlq3MnzfC8iULWbp4waDDkKSh0di2+vFkALB5y1ZWrF5fc0SSVK/GJoTxZNDtvSQ1TWMTgiRpeyYESRJgQpAklUwIkiTAhCBJKpkQJEmACUGSVDIhtFi1dkPdIUhSbUwILXxaWVKTmRBa+LSypCYzIUiSABOCJKlkQpAkAfUvkDNw+51+Sd0hSNJQsoYgSQJMCJKkUuMSwt3nHFt3CJI0lBqXECRJnZkQJEmACUGSVDIhSJIAE4IkqWRCkCQBJgRJUsmEIEkCTAiSpFIjE8L8eSMdX0tSkzUyISxfspD580aYP2+E5UsW1h2OJA2Fxk1/DbB08QKWLl6w7f0nL72txmgkaTg0soYgSXomE4IkCTAhSJJKJgRJEmBCkCSVTAiSJMCEIEkqNfI5hIm88szL2bxl67aH1lqfV5CkucwaQpvNW7Zu+3fF6vU1RyNJg2NCmMB4cpCkJjAhSJKAIUgIEbFnRHw1Ih6MiCci4taIOLLuuCSpaWrtVI6IXYGrgauAY4EHgQXAA3XGJUlNVPcoo9OA+zLzpJZ9d9UVjCQ1Wd1NRm8Fro2Ib0TEAxExFhHvj4jodHJErImINcCigUYpSQ1Qd0JYACwDNgBHAyuBc4D31RmUJDVR3U1GzwJ+kplnlO9viIgDKBLCF9pPzsxRKGoKgB3PklShumsI9wG3tu27DdinhlgkqdHqTghXAy9t27cQ+HkNsUhSo9WdED4LHBoRfxER+0fEHwIfAL5Yc1yS1Di1JoTMvI5ipNE7gZuBTwIfA86rMy5JaqK6O5XJzEuAS+qOQ5Karu4mo6H3yjMvZ9XaDXWHIUl9Z0KYhNNgS2oKE0IX8+eNbHvtNNiSmsCE0MUtZx9TdwiSNFAmBEkSYEKQJJVMCJIkwIQgSSqZECRJgAlBklQyIUiSABOCJKlkQpAkASYESVLJhCBJAkwIkqSSCUGSBJgQJEklE4IkCZjimsoRcShwDHAosBfwXOAh4A7gSuBbmflI1UEOg/1Ov4T580ZYvmQhSxcvqDscSapcTzWEiHh3RKwDfgh8EHge8FPgWuAR4LXAl4F7IuL8iHhJn+IdqNZV08DlNCXNbZPWECLiJmA34GvAScBYZmaH83YB/i3wLuDWiDg5M79RcbwDtXzJQlasXr/dEpoupylpruqlyegrwN9k5hMTnZSZG4ELgQsj4lXAHhXEV6ulixdsax7a7/RLao5Gkvpr0oSQmSunetPMvBG4cVoRSZJq4SgjSRLQQ0KIiJ0i4hMRcWTb/uf2LyxJ0qD1UkN4D/AR4OHxHRExAmyMiLGIWBURfxIRB/UrSElS//XSqfwHwFcz8+YO124B3gr8MbA1Il6WmT+rOEZJ0gD0UkM4ELisw/4E3puZuwEHAOsohqVKkmahXhLCjrQ0F7WI8RdlreDrwBsqikuSNGC9NBk9DLy4dUdmbo2IFwG/bNl9O/DyCmOTJA1QLzWEa4B3tu/MzHsz87GWXY8DO1cVmCRpsHpJCF8E3hIRJ05y3kJg48xDkiTVYdKEkJnfB1YC/yMiPhURO7afU+77EHB19SFKkgahp+mvM/ODEfFr4FTgfRFxKXADRY1gX4oJ7V4InNCvQCVJ/dXzegiZeVpE/D3w58Bb2L5f4efAH2TmdRXHJ0kakCktkJOZ11L0J8wHFgAvAB4Ebuk0JbYkafaY6oppJwL/B7gtM9e1HdshM5+sMjhJ0uBMKSEAX6V4QvnJiLiFoh/hBmAM+N2IOC4zj6g4RknSAEw1ITwfOLjcDgEOB/4DTz+1vKm60CRJgzTVPoRfAleUGwARsRvwAYq1lo+vNDpJ0sDMeIGczHwwMz9G0Zw0o4QQEWdEREbEF2YalyRpaqpcMe17FFNhT0tEHAr8CXBTZRFJkno2pYQQER+LiDdFxJ4dDu8ObJ5OEBGxC3AhRX/EI9O5hyRpZqbaqXwqxXTYGREPUAxBvQHYSrGy2kenGcd/A/4uM6+IiI93Oyki1pQvF03zcyRJXUw1IexCMYndIeX2amAZsGt5/NyIOAG4HvhJZv7DZDeMiKXA/jjthSTVaqqjjBK4o9z+5/j+iFhAkRzGE8VSitrCyET3i4iXAn8FHJGZv+7h80fL69YAR04ldknSxKZaQ+goMzcAG4D/Pb4vIvbp4dLDKKa/uCVi2wJsI8DiiHgvMN+nnyVpMCbtVI6Ib0fEwb3eMCKeExEfAt7Uw+nfAg6i6BMY334CXFS+3tLr50qSZqaXGsLdwI8iYoxiJNBVwE2Z+dT4CRGxF/Aa4M3A24B7gT+a7Mblg26ty3ASEZuB/5eZN/f4M0iSKjBpQsjMD0TESmA5cBZFx3JGxCbgSYoO5XkU01f8uDzv65m5tV9BS5Kq1+sCOT8DTomID1O0+78W2At4DvAwcDuwNjN/PtOAxjuOJUmDNdVRRluAK8tNkjSHVDl1hSRpFjMhSJKAPiWEiHBqCUmaZfpVQzivT/eVJPVJvxJCTH6KJGmYTHvqivJhtN/tdIinJ7uTJM0SM5nL6HnAi7scmzeD+0qSajDthJCZdwJ3djoWEe+adkSSpFpUMttpB3O+D2HV2g2sWL2ezVuKGTrmzxth+ZKFLF28oObIJGl6+tWp/Ik+3XdotCYDgM1btrJi9foaI5KkmelLQsjMS/tx32HSmgwm2idJs4VPKkuSABOCJKnUy4pp+0XElRGxKSJ+HBGLy/3zI+JtEfGuiHhJ/0OVJPVTLzWEc4FDgTXl+ZdGxGHALRRrKF8ArI+IOd+RLElzWS/DTg8HTs3MzwFExOcoEsFm4B3AE8AfAn8eEddn5rf6FawkqX96SQgvAK5vef8Z4P3AOzPzm+W+yyLi2cAyoNEJYb/TLwF8LkHS7NNLk1EAv255f0/5791t5/0d8OoKYpp15s8becY+n0uQNNv0OsooO7z+Tds59wO7zDiiWWj5koVdk4IkzRa9Tl3x/YhYB9xI0ZmcwLM7nDfnp6zoZOniBds1DY03G0nSbNJLQlgKHAwsAk4Adiz3Xx0RPwNuKren+hKhJGkgJk0ImfmV1vcRcQBFclhEkSgOB94+fnrVAUqSBmPKs51m5k+Bn1IMPQUgInYHDgFeVV1okqRBqmT668x8ALi83NSitT/BoaiShplzGfVBpxFH4FBUScPNhNAH3YahgkNRJQ2vfq2Y1mjtw1DBoaiShp81BEkSYEKQJJVMCJIkwIQgSSqZECRJgAlBklQyIUiSABOCJKnkg2k1cH4jScPIGsKAOL+RpGFnQhgQ5zeSNOxsMhoQ5zeSNOysIUiSABOCJKlUa0KIiDMi4rqI2BQRD0bExRFxYJ0xSVJT1V1DGAXOA14HvBF4ClgdEc+vMyhJaqJaO5Uz8+jW9xFxIrAROBy4uJagJKmhhm2U0U4UtZZHOh2MiDXly0WDCkiSmmLYEsJKYAy4pu5AJtLv4aKr1m5gxer1z3g+waeaJfVT3X0I20TEucARwNszs+OTWpk5mpmjFEljoLo9VNZt/0x0SgbgU82S+msoEkJEfBY4HnhjZm4Y9Oe3f6l3+pLv9KTx+F/sVZvoyWWfapbUL7U3GUXESuA44A2ZeXsdMSxfsnDbX+XdvuQ7PWk8CHefcyywfTPV+GubkCRVqdaEEBFfBE4E3go8EhF7lIcezcxHBxVHXV/2UzF/3sgzagfjTUjDHruk2aHuJqNlFCOLvg/c17J9pM6ghlG3yfFsQpJUlbqfQ4g6P382aa/FODGepKrVXUOQJA0JE4IkCTAhSJJKJgRJEjAEzyFo5to7mH0+QdJ0WEOYpSaaMsMpLiRNhwlhlur2XMI4n0+QNFU2Gc1S3Z6u7uX5hG6zqU7EZihp7rOG0EBTTQZgM5TUBCaEBppuc5LNUNLcZpNRw43PpjqRTjOtToXNTdLsYA1Bk5rpIkA2N0mzgwlBk5psRFMvbG6Shp9NRprUTNaLmGlzUzubn6T+MSHMYcMwRXanhX1mwkWBpP6xyWiOmUrTzkybgXpRRXNTO5ufpP6whjAkqvprvnV96Il0Wzu6alUuTzoMNR5pLjMh1Giy5pTp/GU9G9aHrsJMk4N9EdIz2WRUo4maUwb1F/xsUmXTk0NhpWeyhlCjpvw1X5Vem8N6ZV+EtD0TgmaNqhKofRFSZzYZSZIAE4IkqWRCkCQBJgRJUsmEIEkCTAiSpJIJQZIE+ByCNDRWrd1Q6YN3VU/PUXV8E3FqkXpYQ5CGRNVftlVPzzGoZABOLVIXE4I0JPrxZVt1ghkkpxYZPJuMpGnqZxPK3eccO6Prq16prt1M45tIv2MfdnU2l1lDkKapX8mgilld+7n4Ub8XVhrEwk3DrM7mMmsI0jT1KxlUMe151TPDjhvEtOz9in022bxla0+1o6praiYEqQL9bEKZjtk8tfpsjn2mXnnm5bUmQpuMJGlI9GMN8qmwhiBJQ6Lu2pE1BEkSYEKQJJVsMlKjNXGcu9SNNQQ1TtWddk0fN6+5YygSQkQsi4i7IuKJiLg+Il5fd0yau6ocyTGIcfnSoNTeZBQRxwErgWXAVeW/l0XEKzLz/9YanOakukdySMNqGGoIHwLOz8xVmXlbZp4C3Af8ac1xSVKj1FpDiIh5wKuBT7cd+h7wug7nrylfLupvZJLUPHXXEF4AjAC/aNv/C2CPwYcjSc1Vd0KYkswczcxRYKzuWCRprqk7ITwEbAVe2Lb/hcD9gw9Hkpqr1oSQmVuA64Gj2g4dBfxw8BFJUnPVPuwUOBe4ICJ+DFwNvBfYC/jSBNfsPzY2xujo6ADCk6S5YWxsDGD/bscjMwcXTbcgIpYBpwF7AjcDH8zMtROcfwOwG3DnND5ufISS/RDVsUyrZ5lWzzItksGDmXlwp4NDkRAGaXzoatk5rQpYptWzTKtnmU6u7k5lSdKQMCFIkoAGNhlJkjqzhiBJAkwIkqSSCUGSBDQsIbgQT+8i4qyIyLbt/pbjUZ5zb0Q8HhFrIuKVbff4rYi4ICI2ltsFEbHr4H+aekTE4oj4dkTcU5bfyW3HKynDiDgoIq4s73FPRJwZETGAH3HgeijT8zv83v6o7ZwdIuLzEfFQRGwu7/eitnP2iYiLy+MPRcTnytmZ57TGJISWhXj+CjiYYmqMyyJin1oDG253UDwsOL4d1HLsNODDwCnA7wEPAP8YETu1nPO3wCHAMeV2CHBB/8MeGjtSPGj5Z8DjHY7PuAwjYmfgHylmCP698rNOpVhnZC6arEwBVrP97+2b2o6vAN4OHA+8HtgZ+E5EjACU/14C7FQePx54B/CZKn+QoZSZjdiAa4FVbft+Cnyq7tiGcQPOAm7uciwoFjH6i5Z9zwV+BfzH8v3LgQQObznniHLfS+v++Wooz0eBk6suQ4qFpDYBz2055z8B91COIpyrW3uZlvvOB74zwTW7AFuAd7XsezHwG+Do8v2/Kd+/uOWcE4AngJ3r/rn7uTWihtCyEM/32g51XIhH2ywomzPuioiLImJ83cmXUKxXsa08M/NxYC1Pl+dhFP/Dtk5SeDWwGcscqivDw4B/Kq8d912K+cD260fgs8AREfFARKyPiFURsXvLsVcDz2b7cv9n4Da2L9Pbyv3jvgvsUF4/ZzUiIeBCPNNxLXAyRTPFUopy+mFE/GueLrOJynMPijlTtj3oUr5+AMscqivDPbrco/UzmuRy4CTg9yma414D/CAidiiP70Ex5f5Dbde1l3t7mY5P1T+ny3QYZjvVEMrMy1rflx1zG4B3Az/qeJFUs8y8qOXtuoi4Hvg5cCzwD/VENXs0pYbgQjwzlJmPArcAB/B0mU1UnvcDu7WOdilf745lDtWV4f1d7tH6GY2VmfcC/0LxewtFmYxQtBq0ai/39jIdb2WY02XaiISQLsQzYxHxHOBlFB2hd1H8j3FU2/HX83R5XkMxIuSwltscBszHMofqyvAa4PXlteOOAu4F7u5H4LNJRLwA2Jvi9xaK74Ffs325v4iiA7+1TF/eNhT1KODJ8vq5q+5e7UFtwHEUowveQ/EffyVFh92+dcc2jBvwaeBIis7P1wLfoRjNsm95/KPARuBtwIHARRRfQju13OMyYB3Fl9hh5euL6/7ZBliGO1LMwb8IeAw4s3y9T1VlSDFq5v7y2gPLe20CPlz3zz/oMi2Pfbosp/2AUYov939pK9O/LvctoRiCfgXFGgkj5fGRspx/UB5fQjFq6/N1//x9L9+6AxjwL9Myir+axjP94rpjGtat5ctpS/k/w98Dr2g5HhRDU++jGI53JXBg2z1+C/h6+QW1qXy9a90/2wDLcJRiiGj7dn6VZUjxfMja8h73AR9njg45nahMKYbtfpei030LRd/B+bQMHy3vsQPweeDhMqlc3OGcfSj+CHqsPO9zwA51//z93pztVJIENKQPQZI0OROCJAkwIUiSSiYESRJgQpAklUwIkiTAhCBJKpkQJEmACUHqWbmM4nda3o8vMzqQWYMjYnlErIsI/79VX/iLJfUgIn4beC/FVBN1+RtgN4opyKXKmRCk3iwHbszMn9QVQBaron0N+EhdMWhuMyGoUSLi8LKZp9P2pS7X7ECxpu7f9nD/YyLi0Yj4wnjTTkvT0kERcUVEPBYR90XE2e3NPxHxqoj4ZkQ8HBGPR8QdEXFGyykXAa+ICJchVeVcMU1Nczvbry8AcDrFUqHf6HLNocCuwD9NdOOIOAn4MnB2Zv5lh1O+Bfx34FPA0cDHKBZzP6u8/jXAGuBO4IM8vbDL77TcYwz4VRmv60qoUiYENUpmPkwxnTEAEXE2xZfzWzPzii6XHUoxxfJN3e4bEacBnwT+NDO/3OW0VZl5Tvn6exGxM/DhiFiRmb+kmMv/YeDQzHysPO8HbfH/JiJuLGOSKmWTkRorIs6hWIj9zZl5+QSn7gVsymLlvU4+C/xn4B0TJAOA/9X2/iKKRV0OjIjnAYcDF7Ykg24eLGOSKmVCUCNFxLnA+4E3ZebqSU5/DsWiSt0cD9wMTHafX3R5vzfFQjjPomgmmszjFIvBSJUyIahxIuLzFEupHp2ZV/ZwycMUfQjd/D7FCluXRcSOE5zXvnD7+Pt7gEco+hP27iGe5wMP9XCeNCUmBDVGFL4EnAgclZlX93jp7cC8tkXXW91CsbTjAUycFN7Z9v7fUazrva5sJroKOCEiJvvr/yXAHb0ELk2FncpqkpUUNYM/o8gPrR2zt2bmpi7XrS3/fQ1dmnQy87aIGKVYsP27EXFMZv6q7bSl5TDT6yg6st8DnJWZG8vjH6FYV/maiPhM+VkLgEWZeQpF0LsCCyk6oKVKWUNQI0REUDzhOwJ8Abimbduj27WZeTfwY+DNE31GZt4BHAnsy9OjiFq9BTgK+DbFcw1/CXyi5frrKDqW/5liEfhLgVPZPgkdS7GA/DcnikWajsjMumOQhl5EnExRw9izh1FA7deeBXwceHZmPjXDOC4DHsrME2dyH6kTawhSb74O3AssqyuAiFgEvJFiiKtUOROC1IPyL/s/AqZUO6jYHsDJmXlnjTFoDrPJSJIEWEOQJJVMCJIkwIQgSSqZECRJgAlBklQyIUiSABOCJKlkQpAkAfD/AQtRuUeyDzk6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Lmax = 1800.0\n",
    "i = 0\n",
    "s.domain.create_box_array(Lmax, i, s.coherence_func, r0=0)\n",
    "plt.step(s.domain.rcen, s.domain.B * 1e6, where=\"mid\")\n",
    "plt.ylabel(\"$B_\\perp (\\mu G)$\")\n",
    "plt.xlabel(\"$z$ (kpc)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use ALPRO to compute the survival probability "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "energies = np.logspace(3,4,1000)\n",
    "P, Pradial = s.propagate(s.domain, energies, pol=\"both\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the survival probability first using the inbuilt method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEZCAYAAAAzL+qdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeZxcZZnvv09X73u6O+nO1lkhIQmQCAQIi0EviDgo6lxBRpZhxCsi4zIzMnqd0Rm38Tojcp1xVK4zgo6CO25sIhGGLbIEEggJ2dNJOkmn053el6r3/lHbe06d6q69Tqef7+eTT6rOeevUqe6u8zvPLsYYFEVRFMVPlBT7BBRFURTFjYqToiiK4jtUnBRFURTfoeKkKIqi+A4VJ0VRFMV3lBb7BKY6IvIiMBPYUexzURRFmUIsBY4aY9Z47VRxyp6ZDQ0Nc1evXj232CeiKIoyVdi0aRO9vb1J96s4Zc+O1atXz92wYUOxz0NRFGXKsH79ev7whz8k9ThpzElRFEXxHSpOiqIoiu9QcVIURVF8h+/ESUQuFpFfisgBETEicmMKrzldRP4gIkOR1/29iIhrzbtF5FURGYn8/07XfhGRz4rIwchxNojIyhx/PEVRFCUFfCdOQC2wBfgIMDTZYhGpBx4BDgPnRF73N8DHrTXnA/cB/wWsjvz/YxE51zrUJ4C/Am6LHOcI8IiI1GX/kRRFUZR08J04GWN+a4z5lDHmJ0AohZf8GVAN3GCM2RJ53ZeBj1vW00eBx4wxXzDGbDXGfAHYENlOZN1HgX8yxvzUGLMFuAGoA67N5edz81rniXweXlEUZUpyMqSSnw88YYyxrayHgM8BC4HdkTVfd73uIeDDkceLgDbg4ehOY8yQiDwOrAO+5X5TEdkQebg60xP/7eZDfPgHL/CBi5dw++XLsD2RXf0jbNh2lO6BER7Y0klAhKWzajm1tY6KshLe/YZ5VJYFMn1rRVEUX3MyiFMb0OHadtjatzvy/2GPNW3WOpKsyUtx7fN7u/nIvS8SMvDNP+zkvMVNrF82C4DB0XHe862n2XV0wPGa5/Yejz3+5aaD/PDm8ygpcYTWlAJhjOHRrUdorC7j7IVNxT4dRTnp8J1bb6pgjFlvjFkPbMrk9SvnNHDGvMbY899tjevig1s6E4TJzbO7u7n/pQOZvLWSA770wGu8/57n+NNvPs1PnnffGymKki0ngzh1Aq2uba3WvonWdLrWTbQmp1SWBfjIm0+JPX+5I97G4/5NBx1rG6vL+MTly5jTUOnYftfju9FJxoVnPBjiu0/uiT3/54e2YYxhT9cAL+3v0d+JouSAk8Gt9zTwZRGpNMYMR7ZdChwE9lhrLgW+Yr3uUuCpyOPdhEXoUuCPACJSCVxEOPMvL6ya2xB7/HJHL0/u6GLF7Hqe3NEV2/6VPz2DK06fTU1FKTdftJjO3mEuveMPDI+FePXQCbYf7mdZmyYUFpJ93YOMBuO5Op0nhvnSA69x1xO7MAbefuYc7rxmNa5qBkVR0sB3lpOI1IrIahFZTfj82iPP2yP7vyQij1ov+QEwCHxXRFaJyLuAvwW+auK3sHcCbxKRvxWR5SLySeAS4GsAkXVfA24XkXeJyCrgu0B/5Ph5oammnNmWNfRn/+9Zzvr8I4yHwqd95rwG/ufZ86mpCN9DlAVKmN9UzZuWz4q95onXj+br9JQkvH6kP2Hbtx8PCxPAL186yIbt+ntRlGzwnTgBZwMvRv5VAf8QefyPkf2zgSXRxcaYXsIWzxzgOeDfgH8BvmqteQq4BrgReBm4HrjaGPOs9b7/B7gj8vrnIu9zmTGmL9cf0GblnHrH85DlEXr3WfM8X3Ph0pmxx7aVpRSGHR7i5OY7T+z23H60b4TP/fpV/vfPN7O7a+K4oqJMZ3zn1jPGbACS+kOMMTd6bNsMXDzJcX8C/GSC/Qb4bORfwbjyzDn8buuRhO2LWmq45px2z9ecv6Q59vjFSIxDXUiFY2cK4vTfO7rYe2yABc01sW1jwRDXfedZXusM3+889MphfvuXFzKrvjLZYRRl2uJHy2la8fYz5/Cld53OGfMaHNv/9xWnUV7q/etZ2FxNY3UZAD2DY+w5Npj381Ti7DgaFyc7SWXxzBouOqUl9vxHz+13vO7ejftiwgThWrY7frc9j2eqKFMXFaciIyK8d2079996AZ+6YjlrFzbxj+9Yyf9Y4U4cdL5m9fx4GvqL+44nXavkhuGxIPdu3Mdj24443Ho//MB5fODixVxzznzu/vO1vO+8BbF9P3qug+GxIBCui/r+M/sSjvvT5w/Q1T+S/w+gKFMM37n1pisiwgcuXsIHLl4y+WJgzfwZbNgWDrq/sO8473qDd3xKyQ1/+9OX+YVHin97UzWfuuK02La2hkpm1VVwpG+Eo30j3P3UHv7XG5fwUkcv2w6HraaqsgBzZ1Sx40g/o8EQ9286yF9cuKign0dR/I5aTlOUsxbMiD1+asexIp7Jyc/IeDBBmABOmVWbEOsrC5Rw25uWxp7/6+93cKh3iPv+GLea3nbGbG66IC5GP9UiXkVJQMVpinL2whlUloV/fbu6BtjfrXGnfJEsO2/prFrP7desbWfxzHAiRN/IOB+7b5OjsPqac+bztjNmx2KKrx46wdZD2gBYUWxUnKYolWUBzlscz9r7g9bV5I2DPcOe209t9S5+LguU8E/vOiP2/Jld3QyOhmNPp8yq5awFM2ioKuMyK674wOZDOTxjRZn6qDhNYS4+JV7v9LiHOP3shQ6u+86z/HBjYiBeSZ1kCQt2UoqbtYuauM5Kjohyw7qFMVfgFafPjm1/+FV3z2FFmd6oOE1hLj41Lk5P7zpG0KrgPdY/wu0/fZknXu/iUz/fzIGeSec2Kkk42pcoTotaahyNe724/a3LWW61llq7qIn3ro3Xrl186kzKA+Gv4GudfXQPjObojBVl6qPZelOYJTNraK2v4PCJEfqGx3nlYG/sgrntcB9jwbBYGROOm8xtrCrm6U5Zjg/GRePMeQ2saZ/B+85rJzDJuJLailJ+cesFPLDlEGPjhivPnON4TW1FKavm1vPCvh4Ant97nEsnKCEoBmPBEGUBvYdVCo/+1U1hRITzrbjTbyJxi4GRcV6MXPCidPaq5ZQpgyPB2OOrz2nns29fydJZqTXbrSwL8M4183jPOfOpKk8cDmnPgnp+r7/q1f7+/i0s/7sHw3PHQtppXSksKk5TnLesbIs9vvupPdy/6QBnff4RvvLQNse6zl4t9MyUwbG4OFV7CEw22J1Bth/OaxvHtHj9cB/3PL2XYMhw/6aD/HaLJmwohUXFaYrzlpVtrJgdbh47PBbiI/duYngslLCud2is0Kd20jA4Mh577GX9ZIOdjr7z6OQ9+wrFM7u7Hc+f2K4NhpXCouI0xSkpET5z5YpJ150YVnHKlGgaOEBNeW7DtAuba4jW8e7vHoy1Oyo22zqddVc7fCScyvRAxekk4NzFzZxtdYzw4oRaThlju/VybTlVlgWYNyOcqBIysNcnTXwPuWq7dhzp1wm/SkFRcTpJ+NAlzp58DVVljufq1sucodG4Wy/XMSeABU3xsRoHfZLyf7jPKU69Q2Mc01R3pYBoKvlJwpuWt/L1967hyR1dvOec+ayZ38hLHb1c9W9PAnBieHySIyjJsN16+RCnNmvsxqFe724UhebIicQEmj1dA7TUVhThbJTpiIrTScSVZ87hyjPnxJ43VZfHHqtbL3OGRvPn1gPnTCg/pPwHQ8azK0bnCX8Ip5IdI+NB/vPJPRw4PsS73jCXNe0ThwSKhYrTSYzt2lNxypwBy62X64QIgLaGeHH0QR9YTsf6R/Aqa+r0wbkp2fPJn23mZy8cAOC+P+7n/g9fwGmRjF8/oTGnk5jayviFtG9k3NHeSEmNUMg4UvOrynJvOc12WE7FF4DDHi498I/LUcmcTft7YsIEMBoM8a+P7SjiGSVHxekkJlAi1NkCpenkaTNkZepVlpVQMknLokxorbfEyQeus2RlB344NyU7fvzc/oRtD27ppHfQf9cGFaeTnPpK27WnSRHp4kyGyI8XvKU2Hhv0Q/PXfqvouK4i/pn9YNUpmTMeDPHgls6E7cGQ4amd/iuyVnE6ybHjTppOnj6OZIg8uPQAZtTExen44GjR3a+DVoxtsdXBQsVparNxT3esHKC1voJbrfKTJ3aoOCkFpr4qfuerXSLSZ3AsvzVOEB5OGL2JMAZ6BotrPfVbjW6XtMRrsA6fGNYGsFOYZ3bFW1JduqKVdUtaYs9f7ujxeklRUXE6yXG69VSc0sXh1qvIX3Jrc41/XHt2L8Hm2nIaq8N/Q+MhQ9eANhCeqjy3Jy5Oaxc1s2pOvOnwts4+Rsb90ToriorTSY669bLDHpdRnSe3HkCTJU7F7sQwYIlTTUUpbfX+yiZU0mcsGHKM0Tln4Qwaqstob6qO7De8fthf/RNVnE5y6u1aJ3Xrpc1gnlsXRWnykeVku/Vqykt9l+qupM+rB0/EMk/nNlYxO1Jbd9rs+FyyHUdUnCZFRD4kIrtFZFhEnheRiyZZf6uIbBWRIRHZJiLXu/ZvEBHj8e8Va82NSdZUJr7j1OFksJyMMRwpUhrzUB6bvto0Wxl7xzy6MxQSW5BrKkp92V5JSY8/Wi69cxbGO0IsaoknvOzqGijoOU2G7zpEiMjVwJ3Ah4D/jvz/gIisMMbs81h/C/Bl4GbgWWAtcJeIHDfG/Cqy7F1AufWyCmAz8CPX4QYBRwdVY8yU/jbWW3VOXv3S/M7RvhGu/vbT7Do6wA3nL+Af3rGqoO+f7756Ufzk1ut3uPUCzG2Md7A44JPGtEp6PLcnPmXZnr682Ep42e0zcfKj5fRx4LvGmLuMMVuNMbcBh4Bbkqy/DrjLGPNDY8wuY8y9wLeB26MLjDHdxpjO6D/gQqAa+A/XsYy9LrJ2StNQHbecfvx8B3uPDdA/Mu6buUGT8Y0NO9h1NPylufvpvQUfyFeIOieAppp4Q9Viu/Xc86vmzrDE6biK01TDGMNze23LyRKnmbY4+cut5yvLSUTKgbOAf3btehhYl+RlFYDbuhkC1opImTHGy5d1M/CgMcZdLl0lInuBALAJ+DtjzItJznVD5OHqJOflC5prnF2kr/vORg70DNFcU84DH7mIZh93mR4ZD/KT5zoc257d1c2SmbVJXpF78j0uI0qzby2nUscNTsdxf8ybUlJnz7FBuvrDf1P1laWcYtWuLbItp6MDGGMQyX0XlEzwm+XUQlgYDru2HwbakrzmIeAmETlHwpwNvB8oixzPgYicCrwRuMu1axtwE/AO4L2EBe9JETklw8/iC85d3OR4vq97kGDIcKRvhN96VIv7ic0dvfSNOLtaFPrubqAIbr3ufj9l66lbb6pjx5vOXtjkaMHVVFMec/0PjAY50ucf17/fxCkTPgf8BngKGAPuB+6O7At5rL+ZsJvwN/ZGY8zTxpi7jTGbjDFPAFcDO4HbvN7UGLPeGLOesIXlWypKA/zi1gs89+30WXaOm43WlypKR4HdSs5xGfl06zm7RBQTh1uvopTW+kpKIxe0rv5Rx89E8T8bd9vi5ByPISIstKyn/d3+sYz9Jk5dQBBodW1vBTxv840xQ8aYmwjHkBYC7cAeoA84aq+NuA1vAP7TGDNhozljTBB4DpjSlhPgMONt+nw+gPCl/YlV64XOFpueqeTOESGBEmGOWk9TFttyOndRU8J+v1rGvhInY8wo8DxwqWvXpYQto4leO2aM6YiIyjXAr40xbsvpKsKuvu9Mdi4SdryeQdjKmtJUlweoKE38VRf7Dn0y9h5LvIsrdGufQmXrzah2Wk7GFK9N0KDLrQfOC5jGnaYOh08Mx75HFaUlnD63MWGNfeNxsMc/ycm+EqcIXwVuFJH3i8hpInInMAf4JoCI3CMi90QXi8ipInKdiJwiImtF5F5gFfApj2N/AHjUGLPLvUNEPiMibxGRxSKymrCAnRF936mMiDgC7lH8PELDGMM+DxfD8QK39i9E41cI11BFjz8WNA7rJRte2t+TVt+0UMi44mxhV+aC5urYtp1H/ZVyrCTHdumtaW+k3OMm1WkV++fGw1fZegDGmPtEpBn4NDAb2AJcYYzZG1nS7npJgHD6+TLCMafHgHXGmD32IhFZDLyJsFXlRSPhFPQ2oBd4EbjYGLMx28/kB5pqyxOmrPrZrdfVPxqzWqrLA7HHJ4bHCIYMgTzMVfLCHX/JJ0015TG3yvGBMeqsvoiZcO/GffztzzYD8JkrV/DnFyya9DWDY04xjv6cl7XFOwm8duhEVuelFI5ndx+LPV67MNGlB06r2E+Wk+/ECcAY8w3gG0n2rXc93wqsSeGYu5jAUjTGfAz4WFonOoVoqklMGfezONlW06KWGvZ3D3JieBxjwp0umjwswXxgx5zy2SECYEZNWUycugdHabeslUz45h92xh5/9eHtvO+8BZQFJnaWDLrSyKMsb4uP8X6tsy+r81IKgzGGx16Lh93PW9zsuc4Rc/JRHZsf3XpKHvBy6/m5157dw21OY1XCzKNCUaiYE7jiTlkmRfQMjrLHitn1jYyn1NjT3R0iynLLctp2uI/Rca9EWKVYDI0GHTdSAFsP9cVuduoqSznHIxkCcBRZH9SECKXQeFka/SPjvp3P02X1l5tZV0GjdeHuKWDcySFOZfl1NORybMbmA70J21KJPbm7Q0SZUVMe62A9Oh7iVXXt+YZ7nt7Dms89zKrPPMSHf/BC7KbzwVfiCc6XLJuV1GqeUV1GZVl4X9/IuG9uWlWcpgle4mQMDIz607Vni1NLbQUzrC4FhczYK1TjV0iciJsNrx5MFI+XOhIFy00yywngDe3xTK8X9h5HKT4bd3fzmV++wvBYiJCBX798iGvveoYjfcPc98d4K9K3rEzWwyCcMGWPRSlWk2U3Kk7TBC+3HmQXd9p+uI+ndnTlxfpyWE615TRa3dULmbFXqDongKbq3FlOXsXKr3VObu24O5LbvGFBvIDzhX0qTn7gXx/bgbvqYMuBE6z9wqMcjjR6nllXwaUr3KWjTmZZ4nTYJw2iVZymCTNyLE6/eukgl93xONf+v2e55+k9SdftOzbIX//4Jb70260E0xCxLquFT0tthWP0R6Em+gZDhuGxeGwln6nkkFvL6VBvojht7+ybtH7KMcvJLU7tcXGyB9cpxeFQ7xBPvB5OeBCBWy9ZgldbvFvXL/FMIbeZVRdPmDrSp5aTUkCSW06ZXei/+sj22OOfvXgg6bpbf/ACP3m+g289votHXnW3TEyObTk111ZkPDTxQM9QxgPyhlxp1SV5Tl/PZZeIAx4pwQOjwUnbPzn66rksxeVtdbHYxIGeIQ77xP0zXXlie1fMajpvUTN/85blfPndZ8RaTQG8dVUb152/cNJjtfrQcvJlKrmSe9qbqykRcBsv7saqqWJn9bycJJZhjHEE5l852Mvlq5L7vm2cMady6ittyym1c/7t5kP85Q9fJFAi/OeN57BuaQvjwRClk6RTRymkSw/c2XrZWYf272fJzJpY4ez2w33Mb0qeou4e0W5TGijhjHmNscLOF/Ye562nz87qPJXMeXJnV+zxG5fNBOA9Z8/nvEXNPLWzi3kzqlm3pDmlmyqH5eQTcVLLaZowq66SL7zzdC5b0cpCq34mU7deict/4BV3cotIOu64rj7LrVdXQX1V/EKZquX0fx99nfGQYWQ8xJ2Pvs4dj2xnxd8/xAe/93xKcTJn09f8i5PDcsrCrdc/Mh6belxeWsK6JfHm/JPVKA2MeGfrRVkzP54UsVXrnYqKHfdbtyRew9TeXM01a9u58JSWlK19h+Wkbj2l0Lx3bTvfvv5sRzFeJm69E8NjDpcXQL9H1p/bd51qO/6BkfHY8ctLS6irKHVZTpOf88DIuONC/Ozubu589HVGgyEefKWTP7x+dIJXhylkjROEi3CjZFPndMiymuY0VLJ8drxGafvhiQVlooQIwGF1HfJRTcx0o294jP3d4Z9/WUAcRdKZYFtOR9VyUopFnTW6PRPLyatQr9cjg86dVZfqex2zkiFm1lYgImnHnLZOUodj9xxLRqGm4EZxN3/NNAvS7iw9p7GKZa1WAe0k1o6dSl5bkSjIjlY3HkkXSmGwb7yWzKydNOFhMmbVx8VJLSelaNg92zKxnLxanPR6WDP9I85t7gr2ZBztj385WmrDF+x0Y05bPIpQbVLpDzdUYMupLFASu3EImcw7eNj90WY3VHGq1d1h59F+xoLJuzvYMScvQbabhB7yUR+26YZ987VidnZWEzhTyY+cGClqV/woKk7TkGwtJ6+ZL17i5D72YIpD6uyA7My68Jcm3ZjTFo8iVJu9kd59HccH2bi72zPNvdAJEZCbLhF2GvncxkrqK8uY0xD+OY4FzYTW08AkjW5nN8YvYgd6hnxxEZuO2OJ0Wg7Eqa6iNJaJOTQWzDhRKpeoOE1DnJZTBuKUouXkPrY7TpUMOzYVdTekG3OazHLqOD7Enq4BLrvjcd7zrae59b9eSFgzWKApuDYzciBO9s3D7Iilc7bVkfpXLx9M+lr3iHY39ZVl1EZEa2Q8VNBWUkocu0+iHVPMFBFxJEX4IWNPxWka4rSc0r+4dKRoOblnEk1mOf30+Q4+/YvNjiykaKDWcc6T9AQcHB3n9UlG0I+Oh/i/v389dk4PvtKZUA9lt3Zy1/zki5Zay/ef4QXCdrdF3XBvOyOe8v3dJ/fwfJL2Q5NZTgDNtf4ZKT9dsZv6LmyumWBl6jjTyYvvslVxmobUVdgusvxZTv1uy2kCcdrc0ctf/fglvv/MPu7fFL+znxn5wpQGSmICYYx3dmCU3209EnPTLWut4+JTZ3que2Bzp+O5uyWPbfnZ4phPou438O7ykAp2okL0eJcsm8Ups2qBsMXzwe8/7xkDHHAkRHh/ZruVVE+BunUocQZGxmN1gGUBccQBs2GWz9LJVZymIdm69bw6Lni5dxItp/GkMYrfbfXuHjGrLv6FqU+hhdHwWJAvP/Ba7PmbT5vF165ezY3rFvK5d6zkcqsBptvNuN81ede2KrMd/JcqbQ3xC00mnS1CIcMh63VRt155aQnfuu4sGiMNdI/2jfDgls6E1zsTIrytxQYrq9ArS1PJL3stq2n+jOqcDd70WyGuitM0xLYC0u1TFwoZR/eG2HE83IPubSETvmv3Ys8x79Hf9hcmlYy9B7d0xmIuTTXl/K+Ll9BUU85n376S685fOOEAv/3H3eJUeMtptsNySl+cjg2MxmYt1VeWOqyfxTNrufmixbHnf9ieWOuViuXU4LCc1K1XaPZ1x78rC7IcSGnjtxZGKk7TEEdNw4lhxidILXbTOzTGuEe8x8sCc7v1ILlrL9mQM7vo087Y6xka9cywe3pnfCz1n69bSEO10+KZPyO5CyRa1BjFKU6FsZxmZ+nWs1/j5e453+ok4G47ZYxxxJyS1XbZbj21nAqPHW9akKN4E/iv+auK0zSkurw09oc4HjKOupjJ8LKawDuxwu3WAxhMkrHX6RGArasodcxxaqiKu5OuvetZLv4/jyW4vnYejSdCnGm12okyb4K+ch0JlpPt1iuU5ZSdW293V/yuep6HENs1MXuPDThqnkbGQzHBLw+UJC3sbKzWmFMxsd167RP8PafLTLtLRIrdXPKJitM0ZWFL/I5rt8ultvfYAD94dp/nXfHRpOLkYTl5iNNQkkQGr5jV3BlViNXDr9Wy+CCcMv39Z/Y6ttkX58UzE+8q589I/mV2u9FOFMGt19pgV+qPpDVmBMKjuaOc2pqYYlxZFoj9HEPGabE64k0eaeRRHG49tZwKjl0qMFET33Sx47sqTkrRsJu/7rEu6CPjQd7970/zqZ9v5upvP51wcbTnLC22BM7TcvIQLK908vFgyFPc/vyChY7ntk88ysY98TZEvYNjHIvUBlWWlTCnIdFy8LIm7HOzBdU+p/oCufUqSgOxrhjBkEnbvZJKcaZ9t227Midr+hql0U6IUMup4Ng9DW03cLY43XoqTkqRsC0nOxlh77HBmOvutc6+hGLWLuuPdpFDnBLFxStN3Uuc7HWlJcI/vH0l37nhbN5z9nzHOq8Yyj7LxbHPyrZb0FTj2ZG5sizg+BK6se8Yi+HWA6drL1kszs3urgE+eu+LjiSH05IUZ9rWo50EYtd1JUuGAKflpOKUPcaYtDpt2BZ+rtLIIeyuLY+Mk+kfGU+53Vi+UHGapthWj+0Kc7u23MPp7JjTZOLk7q0H3gkRPVYh59wZVdywbiFvPq3V4dIDWDU30RLoPDHMcCSOZffka53gjnKmS5xWW7Epu/iwGAkRAPOb4hccd5KGF2PBEDf8x0Z+YdWHzaqrYHFLref6eQ7LyRKnFN16dueIAR+0uZnKvNzRw/p/3sC5X3yUDduOTLr+xPBYzLqvLCtxxGSzRUR8FXdScZqmtDfFhcW+QHW6MsTcGWP2H+wCS5z6R8YdLsCxYMgx4jzKPU/v4e/v38LBniGCoXCfN/uYdiaYm2WtdVywtDlhe9Riso8zsza5dfTONXNjj9+7tp02u21L5BihkHG4+AppOdmWzT5X7ZUXT7x+NGHdX1y4KOksH/vz2j8zO1NvIsvJdvml2i8xG17rPMEnf7aZHz23P+NO7X4kFDJ89N5N7D02yJG+ET5636ZJrZWDjnEoVQk3cNnS4iPXnk7Cnaa0NSRekCHRAnI3ebUtp9a6CmorSmMX8f6R8ZjLJ9kd9WPbwm6ngz1DzKgu58fPdzj22wWebkSE7910LnuODfCpn2/mmV3heNPeY4Oc2lrnFKcJXHc3XbCIEhEGRsa5+eLFfOE3W2P7osc4NhBPVW+sLqMsxem5uSCZZZOMTft6HM+/dvVq3rF6TtL1LVb7oWNW/75UCnDBZTnl2fXTOzTGe7/9DMcHx/jhxn30DI7ygYuX5PU9C8VLHT3ssrwWPYNj/PfrXVy2Mvm0aLs1ld2EN1f4qRDXt5aTiHxIRHaLyLCIPC8iF02y/lYR2SoiQyKyTUSud+2/UUSMx79K17q03neqMqO6jLJA+K6rb3g85m5zu93cMQ87IaKlriJpn77JOk/8bjbGczgAACAASURBVOuRBGGCiS0ngJISYfHMWhbPjLus9ntZThOIU0mJcNOFi7jtzackxKCiQn3Ycu+11uX+IjARdsJCKpbTZisu+LWrV3PVmrkT3lE3W1alfbPRP8GIdhu7/mlwJL+W029ePuSYC/afT+5JO4PRr/xxT+JMsejNWzLs1lSzPRJ+ssUxdLDItU6+FCcRuRq4E/gisAZ4CnhARNqTrL8F+DLwj8BK4DPAv4nIla6lg8Bs+58xZtg6TlrvO5UREUeT0ehFyl2H5K6Bsi9mM2vd4hS/uGUaKE/Vh+4V1LfT3CcSJzd2UfKRvmGMMY7Egln1qR8rF9iFwu6YnxevWONBVs1tmHS97fK0BzsOptAdApxuvXxbTu4GtYd6h3l217Ekq6cWmw8kjnXZuHviz+Zo6pvDTL0odjp5sd16vhQn4OPAd40xdxljthpjbgMOAbckWX8dcJcx5ofGmF3GmHuBbwO3u9YZY0yn/S/L953SNHhMl3VbTnbMyRjjuJi11FY4EgXsVkiZilPDJJZTFK+kgVRjTm7cQeBP/XwzX3loW2zbKbOyH0mQDuH6rvDjQ71DsXZEXowFQ7GLSIk4SwSSYXcVP9ofHyyXSncIgCrL5Tc4GszrTKdXDiaOPnlyZ1fe3q+QbO7oSdi2u2tgwgbJDssph5l6UWb6KObkO3ESkXLgLOBh166HgXVJXlYBuG3QIWCtiNhXuyoR2SsiHSLyaxFZk+n7isgGEdkArJ7kI/kWr1517i9GV/9oLBvuxNA4o5GOAjXlAarKA7m3nGqSx5xsbMsp2tnh8IkMLSfrbvHJHV38cON+x/417YmdJvJJRWkglrTgLpR1Y98sNNdWUJpCbKymopSqsrDAjI6HYu68gUlGtEcpLy2JpRwHQyZpv8RsGR4Leo4+eXZXojtsqjE8FowNvCwRmBsRmpCBbYeTD4N0xJzyYjmpOE1ECxAA3G2qDwPJIoUPATeJyDkS5mzg/UBZ5HgA24CbgHcA7yUsZk+KyClZvO+Uxu5VFxUTr/ZC0btX223WEpuzZHU4H8necpoxQUKEjV0Zv/fYIOPBkCNO1JbGF9cWMnc4Y+3CJi5fVfhff7JaJDeZWou29RSNI042ot3GTjXPV8Zex/HBWHzJnqf1UkdP0WtwsqXj+CBRg3N2Q5XjBui1Q8mnODunHOch5lSvqeS55nPAbwjHiMaA+4G7I/tCAMaYp40xdxtjNhljngCuBnYCt2XyhsaY9caY9cCmLM+9aNR7uvUSv/T/8d97gMR4EyQf+W63tbGzwyYjVcupqaY81oZnaCzIlx98LXYHX1MemDBm4qa5phyv/IHb3rSU77//3IJm6kWZZ7ktJ0qKsDtIpGMt2uPgowMD+0dSSyUHV9wpT7VOdrxz1dwGTm0NJ8GMBQ0v7E10iU0l9jqat1azxErwcbcTi2KM4aDHOJRc4mxhpAkRbrqAINDq2t4KJA6gAYwxQ8aYm4BqYCHQDuwB+gDP9BdjTBB4DohaTmm/71THa/S5113wQ690MjQadIhTyyTiZFtO6WQVpVNUeO6ieM3TXU/sjj2eqADXi9JACc01iRf2t66anbT5ab5J1mLITaoZim7sm4DjkXRy2xqZKFsPnKnm+bKc3B3W1y6Kj5rf6JHpNpVwi5Nd0L63y/tmxB6HUucah5IrmmvjN2rHBkbTmliQa3wnTsaYUeB54FLXrksJW0YTvXbMGNMREZ5rgF8bYzx/uhLOtT2DcMJDVu87VXFaTpGYk4dbbzxk2Hm039G6qKUufHFzCNywt1svHd94qm49gA++cUks9mGTSer3XFfNSGVZSexOvRhk5NZLQ5yarJ9z90DUckqtQwQ4xStfGXsHHePmKzlnYVyc/rh7aouT/Tttb6pJ2k7Mxpmpl3urCaAsUBL72zDGWTpSaHwnThG+CtwoIu8XkdNE5E5gDvBNABG5R0TuiS4WkVNF5DoROUVE1orIvcAq4FPWms+IyFtEZLGIrAa+Q1icvpnq+55s1HsMHbQTIuy7951H+501TpNYTicyFKfGNCynFXPq+eVtFyRsd3cvTwX7zhXgzHmNKSUX5Iv5KRbiHvVwtabCDA+33mCKHSLAWYibr1qnQ66aHttyenH/8QmzGP2OfVPR1lDhbMR8bMAzA9Lx88hDAW4Uv7Qw8qU4GWPuAz4KfJpwTOdC4ApjTHQ+QnvkX5QA4TTwl4BHgEpgnTFmj7WmkXB6+VbCGXhzgYuNMRvTeN+TCq+Yk32BskdOHO0bcfyhRsWpPsnId9tyakvjLi9dV8XytnrWLXG2NMpkxs2pbc508bMXzkj7GLmkPVVxytStZ90EdA+Ef1epdogAqCrLv+VkZ1/ObqhkdkNVrKv88FiILR5p5lMFt4u8sbo8dmM2PBbynERr973MRwFulFmOdl7Fizv5UpwAjDHfMMYsNMZUGGPOMsY8bu2LJiNEn281xqwxxlQbYxqMMVcZY7a5jvcxY8yCyPFmGWPeYox5Op33Pdlw1Dl5xJwWWlM2j/SNOLP1PC2nuCDZ47vnpHGXl0mvsEuWzXI8TzYqYiIuXNrieH6B63mhmVVXQWmkN97xwbFYOr+bXMScoo13U+1KDmG3Z5R8pZIftxoCN0XOd+1J4trz8kLY3ze7GXMUu8YpHwW4UWwLvJjp5L4VJyX/2MISjTnZF8EFlqvhyR1d/P61eNfk6B1sKkW4qd7leTV1TYU/OXN2rG6nqaacC09JX1jOmNfILeuXUF5awrXntnP+4szOJVeUlDg7eCRzr2TaFcMr5uSY5zSJOFWUxi2rZMKZLbY4RWOR51iuvWemcKcIr+Qi27W3rztRnOyYUzqlEunil3Rybfw6jXFn6xljHBlb9p2c3SIHoD3yRbJrpRxuPSuVfLKuBTOqy+gbHucvLlyU5icIM7uhip/ccj6/fOkgV62em/F4i9svX84n3rIs552eM6Wlrjw2vr6rfyRh6qkxxtGcc6I5VW68Yk6O3nqT1Dk5LKc8iVPPQPxvKCpOtgv3yZ3H6BseK+g4k1wwFgzFSi1E4i7Wduv75lU+4M5ezBfOQtziufVUnKYx7mSGkfFQrAi1PFCSNOi6ZGZNTNi8svVCIUOfdaFrmqB2aU5DJY/+1XoGR8cdDUnTZeWcBlbOmbyv3GT4RZgAV+/DxKypgdFgLLuysqwkrXid/TvpjqQoRxMMAiXiEB8vKsvillM+3HpjwVDsb6hE4n+rC5prWDG7nlcPnWB0PMSjW49w1Zq5bNh2hK8+sp3DJ4a5du0C/vLNS331u7Sxu3o0VZfHEm/aXYXlbg7muTtEFEd/vSJ2JldxmsY4ujsMjzky9SrLShx/pDafftsK6xiJllPf8His+r22opTSQAnnLW6KjbiwmVlXQVWkFZLiZDK3njvelM7F2E7ZPz445mpdVDrpsSqs+q98uPVsl15jdbljNtXbzpjNq5EuCj/YuI+R8SCf/Nnm2I3VHb/bTmt9Bdes9We/Zi+XHkycBBMMGUcHlHwmRKTaX++xbUf4+QsHaK2v4IKlLax3xX6zRWNO0xhbWPpHxh0B8eryUsdYjSifv2oVlyyP/xHWlJcSvW4MjgYZC4Yc9U7R9/j4pctYOaeeG9ctdByvokxFKRn2RcK+oEXJtHUROFP2ewZHHTHCVCww23LyGiqZLXaHEXd5wTtWz4kli2zc3c3tP92c0Hbqjt9tZ6yIBaQT4RCnuvhNgh3j3esSp67+EcYjH3JGdVleb+ZmpZhKvqWjl1++dJC7ntjteeOZLSpO05iyQEnMfRMyTndDdXkAEUmwntyZcCUl4rLAxl3jzcMXurWLmvjNX17EZ9++0vF6v15A/IDXSBObTDP1IPy7j/5u3M1lU5n668zWy4PlNJCYDBFl3oxqrj030SpabNWqHT4x4ttsPq9MPQjfYEQt0p7BMccNQ6HSyCExISJZ1/lOe+ZZHsbKqDhNc+yYkW3CR+/M3CMsvLom2EkRJ4bGHCnlXsFqe0rrh9YvzeCspwctjuasXuKUWV+9KHbcyQ7Ap2I5ObP1cn+DYQ8Y9Gpp9akrTuN/nBbvNHbZilbu//AFXH/+gti2R7a6ezj7g2RuvZISSeraO2SPZ89jAS6EvSbRv4HRYChpE2e7FqutPvfnpDGnaU5dZann9NdoavbSWbUx/354feKFIixw4S/PieExh+VU73EXfvvly6mtKGXFnHouXeFuZahEsV11XX2JCRHO7hDpXxxmVJfHAu8OcUrTcspHzKnHFXNKfP8Ad11/Fh3HhwiUSCx77U3LZ3HP0+Ga+Y1+tZysm8BmV1Pk9qbq2JiQfd2DseGRdsPXfKaRR5lZVxHL3jzaN+L5O7Az+WblQZyytpxE5CEROS0XJ6MUHltsjtjiFLGcrloTt3L+5i3LPI/hngtlj87wErM5jVV84Z2n82fnLkjYp8RxtJHJsVsPsrOc8p2tN5nlBOHMyvlN1Y606jXt8c4er3X25a0GKxuODXi79SBeogHOjD3bcsq3Ww9SS4rId/ZgLiynDYTnIv0A+IwxZupWxk1D7PiCbaZH29dcsmwWd11/NiPjQd52+mzPYzjcesNjscGF7uMr6eGIOXlcII5kKU52LMd2IaXyO8t3Ee5kllMyGqrKWDKzhp1HBwiGDFsO9HK21VXCD3iNnoliu/XsGwY75pRvtx5MXus0Mh6fUhAokbRq7FIla8vJGPMlYAXhcRXbRORjIqJXpCmCbfUctv4Io8PmRIRLV7TyJ2fMSZpe7O4SMVnMSUmNhqoyApGstL6R8QQRsC2nTC4OTTXx303aMSfbrZcHy6lnMLEAN1VOnxuvd9t+OHGSbrHx6lEZxc7Ys28Y7N/PvBnp945Ml8lqnTotsWytS20Cc7rk5IjGmM7IPKXLgKuAV0Xk7bk4tpJfkllOlWmkeLtT0r2y9ZT0KSkRx1DA7gFn3Clbt16jq9YpSm3F5DcUlXm2nGzXsG2Zp4I9uG/XUf+JkyNbry4x5hRlb6SFkTHGMUbD7tySLyZrYXSgJ//dKnIiTiJSIiJLgdnAr4By4Oci8pqI/EBEPpGL91Fyjy0edsxpsq7UjmNUOMXpxCQJEUrq2F0z7FT/YMg4YhfuwHoqJOvckW5CRD5iTv1p9Plzs9gSp50+E6dgyNA9YCVEuIZc2lbRwZ5hxoIhjg/Gk4xqygNpTZbOFNsSt1PG7XOLki9xyvrKISL/DtwEDAA7gNeB70b+PwGcDrwh2/dR8oPtdrMvdumIk30x6x8eV7deDkmWTt49MEowUpTZUFXmiAGlSjJ3WSrTiO33y0dvPXfHinRYMituWezy6O5dTI4PjsYKhhuqyhImLVeWBWirr6TzxDDBkOFgz5Dje7mguaYgbZnspAu7Bi6KI0EjTzGwXNzWvgc42xizOcn+X+XgPZQ8kcztlk4Fuu0Gcrv10nXJKE6SFeLaaf+Z1pgks5xmTNALMUq+U8kH0mhC62Zhcw0i4Umu+7sHGRkPZiTe+cC2fpNZu+3N1TFrZV/3oOP3vrAl//EmiE8dAKcLL4o9vmOuj916LRMIk+Jz3EW2UarTiDnZllPfyLirfZFaTtlgx5zsWIUjIJ1hGq+dEOHYnkICQr5TyfuzsJwqywKxi2vIeDdRLRbJCnBt7LjTnmOD7OlyjnQvBG0NlbG2ZEf6RhK6gDjTyP0rTg9qndPUJdkXJB3LyRFzStK+SMmMljo75hS/sHU6LKfM0niTufUm6iIfJd+NXx2WU0X6Vo+dNDDRJOFCM1EaeRQ7oWNb5wl2HInHzRYVyHIqC5TQGrHIjXHOkgKnqy9fqe25EKcNhOuc/lVEijuhTUmb1iQuoXQsntqEbD21nHKFbTnZsYdcuPWSWc2piFO+G7+mM/jQizkNE7ulioUzjdz757xiTrx/5asHT/CKNY4+kynPmWK76+yfoTHGIU6+detpndPUJlnDxsYkFy4vatVyyhupxJwydeuVBkoSBKq8tCSlZBinWy+YtDloJoyOhxiNNAQuLRGHlZYqc+2YyXH/iFOypq82p82uiz1+YV8PeyJuyUCJcGprnedr8kGyn2HP4BgDkfE61eWBpDc52aJ1TtOchqoyz+ys+gzFqWdolMHIH64I1KYZzFacJBs4eMDRziZzt4rbSmqqLk8pGyxQIrFxKiEDY8HciZPTpTf5bCkv7Lv5Dh9ZTscc4zK8xWlWXaXnTePSmbVp1R9mi+NneDzuGt1t1VzlM3tQ65ymOSLed2Pp3A0lK+StrSh1DIlT0qc5SSp5roLk7huTVFx6URwtjHI4NiObZIgo9l2/Vyp0sbB/h80T/KwvWNqSsO38JYWNmrgTM6LsPhoXJ3tMSa7JRePXfyfcknoj8BlgDeE6p+uBvwG2oHVOvubNpyVOsEznIpUsJlCv8aasscWpe2CUUMgwPBaMpfKWiPMiki7uzhJ2CvFkOApxcxh3sodeZpIMAa54iV/dehN09bj4lJkJ2960PLeTZidjySyr00ZXPCljt1U7tiiP4qR1TgpXnjmHL/72Nce2dILQ0aGF7sC4xpuyp6I0QF1lKX3D4wRDht6hMbr6R4iGeObOqEoo5EyHcDeF+NyjdIQuX81f3W69TIimQodMPBXaD7VOqWTrAVy2spWmmvJYy6q5jVWct7iwlpOdNbjzyAChkKGkRBzitLCYlpOI/IWIPCgiT4nIHSIyx7VE65ymOLMbqnjrqrbY8zeemnjXNhle/dhUnHKDbd0c6h12tOTJts/a8janS9fOFJuMijxNw7VbF2Xq1isLlDiyGO26sGJhjHEU4SZLiIBw4+VvXXcWy9vqOLW1ln+9dk1WNyGZ0FRTHnP7Do0FY+ULvrCcRORm4FvWpvOAa0TkAmPMLgCTyzQdpWh88q2nsbtrgPLSEj55xfK0X19XWZowrVXTyHPDvBnV7Ir4+Q/0DPFSRzy1eEWWqcXrT51FbUUp/SPjVJaVcHEaNyaVeZqGm013CJvWhsrYkL7DJ0ZYUICGqRNxYmg8loVYUx6YtJbwnIVNPPjRiwtxaklZMrOW5/YeB8J9ClvrK/0hTsCHgP2EXXcdwOXAV4F/Ad6Zt7NSCk57c3VWXwSvO1y1nHKDHQfqOD7IZkucTp/X4PWSlGmoLuN7f7GWx147wltWtU14N++mMm+WU/ZuPYBWa+zDYY/mpYWma2DyTD2/sXhmTUycdhzpZ3ZDFUMRF+6suoq0YtPpMpmduAT4ujHmWWPMAWPMd4DPAm8TkbxJpoh8SER2i8iwiDwvIhdNsv5WEdkqIkMisk1Ernftv1lEnhCR4yLSIyKPiciFrjWfFRHj+teZj893MqLilD9scdp1dIBN+3tiz8+c15j18de0z+Djly1j5Zz0hM52M+WyhZGz6WvmcSJ7nLkvxGmCOU5+ZVlb3DLfcsBZELwyDRdwJkwmTrXAQde2BwlbXN4zu7NERK4G7gS+SDjz7yngARFpT7L+FuDLwD8CKwlnDP6biFxpLVsP3Ae8CTgX2AY8JCKnuA63jXA6fPTf6bn5VCc/XmMWNFsvN9hjFH72QkfMspjfVJVWdl2ucXQmz5M4ZWM52TOJfCFOdtPXPFocuWT1/PgNy6b9x3m5wxan7Kz2ycjkN38k8n++ZPPjwHeNMXdFnt8mIpcDtwCf9Fh/HXCXMeaHkee7ROQc4HYimYLGmD+zXxARtKsIuylft3aNG2PUWsqAOk/LScUpF9i1JNHKfIBLls0qyPiEZNidG3KZSt6XF7de4sC8QtOVQgGu31g5p4FAiRAMGXYeHWBg5FBs3xsWZG+1T0Qq6R9vE5G3iog7Uprz1BERKQfOAh527XoYWJfkZRWA+7ZoCFgrIsmujuVAJXDctX2xiByMuBTvFZHFE5zrBhHZAKxOtmY64Wk56biMnLCsrc4R34ny1lWzi3A2cWy3XjTQnwuymeVkY/eN9BqYV2hS6UjuNyrLAqyeHxeh6M+xtERYuyi/qe2pCMx7gV8DnSKyB7gbMMAKEZmR4/NpAQLYhRdhDgNticsBeAi4SUTOkTBnA+8HyiLH8+LzQD/wS2vbs8CNhK2pmyPv95Q2s00Nr4tIY9XUcF34nbJASUJsaX5TFecuairSGYXJ18DBbJu+RmlriAvAEV+IU9ytN7MA02xzxWUrWhO2nb1wRlY3DqkwmTg1AJcAfw38EBgkfPEWwnGhLhHZLyK/EZEvicg1eT1bbz4H/IZwbGoMuJ+wgAIk3M6JyEeA/wW8yxhzIrrdGPOAMeZHxpiXjTG/A/6E8M/nBq83NcasN8asBzbl8LNMWbwsp8YUJqoqqXGp6wLxkTefWvTWUBV5GtXen6OEiFn1TrdesatepqLlBPAnZ86hPOCUivedtyDv7zuhOBlj+owxfzDG3GGMeZ8xZgVhwboI+BjwPaCHcMPX24H/yvJ8uoAg4JbqVsAzFmSMGYo0na0GFgLtwB6gDzhqrxWRjxK2mq4wxmyc6ESMMf3AK4A7aULxwCvmlK9uxdOR689fyLvfMI/G6jJuumAR71ozt9in5LhgjfowIaKuopSqSKPUobGgI5ZVDKZizAnC3Sn+5i3x/LfzFjdx+cpkjqzckfZv3hgzADwZ+QeAiFQBZ5JlDz1jzKiIPA9cCvzY2nUp8NNJXjtGuBaLiAX3a2NM7BsjIh8H/gF4mzHmvyc7FxGpBJYDj6X7OaYjXpaTilPuKC8t4V/ec2axT8NBviynXImTiNDWEC8aPdw7XNQM0qlqOQHcfPFizlnUxPHBUS5c2kJpIP/dKnLiNDTGDAHPRP5ly1eB74nIRsIC+EFgDvBNABG5J/Ke10een0o4PfwZYAbhbL9VWO44Efkb4AvA+4DtIhKV/SFjTG9kzT8Tzu7bB8wC/g6oIe4iVCbAq32RuvVObpyp5P7qSh5lVl1FXJxOjHBKAech2RhjUho06GfsxIhC4Lt0KmPMfZEkhE8TrjXaQtgNtzeyxF3vFCAsSMsIx5weA9YZY/ZYa24lnCBxn+u1dxNOggCYRziu1kLYHfgMcJ71vsoEuC8igRLJe8BUKS52Knlu3Xq5SYgAZ8ZeMWud+kbGYy2eqssD+t1IAV/+hIwx3wC+kWTfetfzrYSLdSc63sIU3rMYyRwnDe5uEI1VZUWtwVHyT0UhOkRkOazSHtpXzHRyO1twVl2FfjdSoLBtbpWTFvedoMabTn6c4pQbt54xJifznKLYllM+08mf3XWMezfu41Cv9+yoI1YRsJ1FqCTHl5aTMvVwi1FbFqPDlalBeR7cekNjQUKRjO+K0pKsA++FKMT92u+287XfhRvN1FWW8uMPns/yNmcDncN9TstJmRy1nJSc4E5+qCwr/mA3Jb/ko7deLpMhwB1zyn0Lo+f3dseECaBveJzP3P9KQk2V/d6tajmlhIqTkhNEhEuWxTtc/c+z5hXxbJRCkI/eerlMhgAcAwfz4da745HXE7Y9u7ubLQdOOLY53HpqOaWEuvWUnPHpP1lBU81OzpjXwOWr8l+kpxSX8jzEnHJV4xTF7kx+pG8kNmo8F/xxTzf/vaMLCGennj63ITbO5P5NBxyztmy3nlpOqaGWk5Izlsys5V/ecyY3rFuo2UjTANutl6vGr7lqXRSlsiwQi4eOhwzHBkYneUXq3PHI9tjjq1bP5SP/I95M5tHXjjjWHlXLKW1UnBRFyQhHh4icufVyazmB07WXq1qnp3Z28dTOY0DYarrtTUs5f3FzzNW5u2uAAz3xzL2DVhafZuulhoqToigZYffWy0dCRK7EKRdDB0Mhw5M7uvjhxn08uOUQf/eLLbF9f/qGeSxsqaGyLMBaq1P8kxGX33gwxKHe+PsWc0DkVEJjToqiZIRtOeUqldxOiMi2ADdKW5YZe2PBELd8/wV+t9U9yQeqygLc9ualsefrlrTwxOthUXpqRxfvOXs+h3qHCUby42fWVWgma4qo5aQoSkbko7dePtx62bYw+t7Tez2FCeDzV61i3ozq2PMLlsbHvz27uxtjDPu7B2Pb5qvVlDJqOSmKkhH5aF+U64QIcLYwSlecQiHDtx7fGXu+am49FaUB6ipLuXHdQtYvm+VYv2J2PXUVpfSNjHOod5h93YPsP26JU1M1SmqoOCmKkhH56BDhN8vpub3HY67A5ppyfnrLOofF6KY0UMI5i5r4fSRb75ldx9h5dCC2f4GKU8qoW09RlIzIh+Xk7KuXD3FKL+b06Gtxd95bVrVNKExRzrWSIp7Z1c22zr7Y82WutkZKctRyUhQlI/IRc+q3EyJ8YDm9sPd47PHFp8ycYGWc8xbH407P7DqG3cloWVtx5klNRVScFEXJiLJAvNB6LGhy0n0hH269ltpySgRCBo4NjDI6HnK4JJMxOh7ipY7e2POzFsxI6f1WzqmntqKU/kjcKUp5aQkLm9Wtlyrq1lMUJSNExDlwMAddIpx1TrlJiCgNlDjGoh/pS8162nGkPxZLm99UxcwUOzuUBko4f0lzwvY18xsLMt78ZEF/UoqiZEyum78O5LgreRR7hEuqcafXj1ixotb0YkVXnjknYZs7s0+ZGBUnRVEyptyOOwWzjzvlw60HMKsu/e7krx/ujz0+tbU2rfe7bEWro4deoES44nRthpwOKk6KomRMri2nfCREQGbj2rcfjltOp6QpTpVlAe68Zg1zGiopD5Tw2StXsKC5Jq1jTHc0IUJRlIxxNH/NQTp5viynTFoY7TgSt5xOmZV+lt35S5p54vY3EQyZlBIwFCcqToqiZIzd/DXbQtxgyDA0FrecqnPYg641zaGDo+Mh9hwLF8+KhMfBZEKgRAjkaH7UdEPlXFGUjKkoy12tk6MAtzyQs6GAAK1WQkQqbr2DPUNEerXSVl9JVbk2ay00Kk6KomRMLrtE9A/HxamusiyrY7lJt7+ePYtJR1wUBxUnRVEypiKH/fUcTV8rcxtxaK1LL+bUYTVrnduoPhju4QAAFN9JREFU4lQMVJwURcmYXFpOfcP5qXECaKwuiyUl9I+MO4TQi47jtuWkXR2KgW/FSUQ+JCK7RWRYRJ4XkYsmWX+riGwVkSER2SYi13usebeIvCoiI5H/3+naLyLyWRE5GDnOBhFZmevPpignC7nsr2cLRl2OLScRYbYVdzpkue28cIqTWk7FwJfiJCJXA3cCXwTWAE8BD4hIe5L1twBfBv4RWAl8Bvg3EbnSWnM+cB/wX8DqyP8/FpFzrUN9Avgr4DbgHOAI8IiIaLdGRfEgl2Mz+vNoOYHTPWeLjxcHrP1zVZyKgi/FCfg48F1jzF3GmK3GmNuAQ8AtSdZfB9xljPmhMWaXMeZe4NvA7daajwKPGWO+EDnmF4ANke2IiEQe/5Mx5qfGmC3ADUAdcG0ePqOiTHly69Ybiz3OhzjZFlDHpJZTPOakbr3i4DtxEpFy4CzgYdeuh4F1SV5WAbhTcIaAtSISTfs53+OYD1nHXAS02WuMMUPA417vG3H5bSBshSnKtMTZISKXbr3cZusBzG2Mi4wtPm5Gx0OOdPM5jZVJ1yr5w3fiBLQAAeCwa/thwuLhxUPATSJyTiRudDbwfqAscjwir53omG3WtlTfV1GmNeX5SojIccwJnJbTgQncep29w7Eap9b6ipQGDCq5x4/ilAmfA35DODY1BtwP3B3Zl5sRnS6MMeuNMeuBTfk4vqJMBewLdy5TyevyEXOakVrMqaNH08j9gB/FqQsIAq2u7a1Ap9cLjDFDxpibgGpgIdAO7AH6gKORZZ2THLPT2pbS+yrKdCdfRbh5t5wmiDlpGrk/8J04GWNGgeeBS127LiVsGU302jFjTIcxJghcA/zaGBP9xjw9yTF3Exah2BoRqQQumux9FWW64nTr5S7mlI+EiLb6ylifu6N9IwwniZFpGrk/8Gvj168C3xORjcCTwAeBOcA3AUTkHgBjzPWR56cC5wLPADMIZ/utIpxtF+VO4HER+VvgF8A7gUuACyPHMiLyNeBTIvIasB34NNAP/CCfH1ZRpiq57BDRl8cOERCeUNtWXxmzmg72DLHYo6GrppH7A1+KkzHmPhFpJiwOs4EtwBXGmL2RJe56pwBhQVpGOOb0GLDOGLPHOuZTInIN8HnC9VA7gauNMc9ax/k/QBXwb4RF7lngMmNMH4qiJOBs/JqtWy+eSp6PmBOELaGoOHUc9xYnTSP3B74UJwBjzDeAbyTZt971fCvhYt3JjvkT4CcT7DfAZyP/FEWZBHtkRtbilOdUcohYQrvDj5PFnWy3niZEFA/fxZwURZk62MMGs3br5TkhApyWkFet03jQWeOkMafioeKkKErG5LS3Xp7bF4FTbPZ1J1pOnSeGCUaKnFpqK6jM4cBDJT1UnBRFyZhcFeEGQ4b+0fyL06KWmtjjPV0DCfs1GcI/qDgpipIxuapz6hsew0S6MtRVluZttPnC5rg47e4awETfNIKmkfsHFSdFUTImV+LUMxjP1Guszk8yBEBLbXksE7B/ZJyj/c7Bg44JuJoMUVRUnBRFyZjyCRq/jgdDdA+MpnSc3iFLnKrKc3NyHogIi2barj1nUoQzjVzFqZioOCmKkjGO3nrBuOV0YniMd37jKd7wuUf419+/PulxeixxaqjKn+UEzrjT7q5+xz7bctKYU3FRcVIUJWOcIzPi4vSjP+5n84FeAO589HV6Bie2oOz9DXl064Ez7rTLlRRxQPvq+QYVJ0VRMsauc7JjTht3d8cejwUNL3X0Tngcp1svv+K02HLr7T4aF6dQyHCwJ17jpAW4xUXFSVGUjKkI2CMz4jGnHUed7rJtnScmPE6hEiLA6dbbaZ1n54nhmGtyRnUZNXlKZ1dSQ8VJUZSM8bKcQiFDh6vA1bZIvLDFKd8xpyUza5FIpvqeY4Ox7uS7LCvKq+eeUlhUnBRFyRh3bz1jDF0DI47kCIBDvcnnJwF0D8RTuhur85etB1BTUcqCpnA8KRgy7DgStp5sK2qxZV0pxUHFSVGUjCkpEcoC8YLZsaDh+MBYwrpDvRNbTl398YSImXUVuTvBJCxvq4893noo7HLcZYnTkllqORUbFSdFUbLC3V/vuEdm3mRuvaN9cctpZm0BxGl2Xezxa53hiTh25p5aTsVHxUlRlKxwd4nwShvv6h+ZsDGs3amh0JbTlkjK+84jajn5CRUnRVGyotw1Dff4YKJbD/B09wGMBUMxa0sEmmryG3MCWD2/MfZ40/4ejvaNcDDieiwLCO1NWuNUbFScFEXJCrfl5OXWAzg2MOK5vXtgNNb0tam6nLJA/i9LbQ2VMQEaGQ/x/Wf2xvadNru+IOegTIz+BhRFyQp3zKknTcvpyIm4aLUUIN4U5ZyFTbHH3358V+yxbVUpxUPFSVGUrEhw6yVp9tqdxKKym60Wsp/duYvi4jRkNa1d067i5AdUnBRFyYpEt17cQrItoWSitbc7Lk6FjPW8+bRZlLrmRpWWCJcsm1Wwc1CSo+KkKEpWOLpEjDmz9ew+dseSiNM+S5zmF1CcmmsruGxlq2Pb+mUz814ErKSGipOiKFnh7BLhrHNaYolTMstpvyVOCwqcJXf75ctpqQ2LUUNVGX/71uUFfX8lOdrZUFGUrHDMdBoPORIillg96rxiTsaYWBEswMKWworTguYaHv/EJby4r4cz5zdSq81efYP+JhRFyYpKy603OBp0DA5cPInl1HliONYdoqY8wKKWwhe/VpeXcsHSloK/rzIx6tZTFCUrqi1r43DfMMFQuGippjzArLrK2D6vke2b9vXEHp8+r4GAK0FBmb6oOCmKkhU15XG3nj1JtrG63NHtwUucHn71cOzxmvYZeTpDZSriS3ESkQ+JyG4RGRaR50XkoknWXysim0RkUEQ6ReT7ItJm7d8gIsbj3yvWmhuTrKn0fldFUSDsFotyoCcuTjNqyhzidHxwFBNtBQEMjQZ56JXO2PO3nT47z2eqTCV8J04icjVwJ/BFYA3wFPCAiLQnWX8B8D3gbmAlcBWwAvgva9m7gNnWv4VAH/Aj1+EGXetmG2MmbqesKNOcmgpvy2lGdTmVZQGqI5bVWNDQPzIe2//oa4cZHA0Xvy6eWcPKOfFmrIriO3ECPg581xhzlzFmqzHmNuAQcEuS9ecDHcaYO4wxu40xzwBfB86NLjDGdBtjOqP/gAuBauA/XMcy9rrIWkVRJiCZ5RStF5pR7e3au3/Twdjjt585BxGNNylxfCVOIlIOnAU87Nr1MLAuycueBGaLyJUSpgW4BvjtBG91M/CgMWa/a3uViOwVkQ4R+bWIrJngXDeIyAZg9QTvoygnPbblFLWEAGZUh8ete8WdhkaD/GH70dj2t585J9+nqUwxfCVOQAsQAA67th8G2hKXgzHmacJi9F/AKHAUEOAGr/UicirwRuAu165twE3AO4D3AsPAkyJySiYfRFGmC7blZBO1nGxxOhaZeLtxTzej4+FR7ktn1bJ4ps5PUpz4TZzSRkRWEHbjfY6w1XU5YSH7VpKX3EzYTfgbe6Mx5mljzN3GmE3GmCeAq4GdwG1eBzHGrDfGrAc25eJzKMpUpSaJOEUtJ3t4YFdkqOBTO7ti2y46RWuMlET8VoTbBQSBVtf2ViBZ/OeTwEZjzFciz18WkQHgCRH5lDGmI7ow4ja8AbjLGDPucawYxpigiDwHqOWkKBNQZaWS20RjTXbz16g4vXrwRGybPbpCUaL4ynIyxowCzwOXunZdSjhrz4tqwoJmE33u/nxXEXYdfmeyc5FwdPYMwlaWoihJsGNONo0ellO0G8TWQ3FxOm22ZukpifjNcgL4KvA9EdlIONnhg8Ac4JsAInIPgDHm+sj6XwF3icgtwEOEU8C/BrxgjNnnOvYHgEeNMbtc2xGRzwDPAK8D9cBfEhanZFmCiqIwkVsvajnFY05H+0c40jdMVyT2VF0eKHizV2Vq4DtxMsbcJyLNwKcJC80W4ApjTHSOcrtr/XdFpA74MPAvQC/we+B2e52ILAbeRDh5wotG4NuE41W9wIvAxcaYjbn4XIpyslI9iVvPEXPqG3W49Ja31VGiLYsUD3wnTgDGmG8A30iyb73Htq8TToqY6Ji7mMCNaYz5GPCxtE5UURRqK70vI00Ri2mmFXM62j/C1kPxLuTq0lOS4auYk6IoU4+K0gBVZU7rqbo8EBs/4bScRhzxphXaFUJJgoqToihZE01+iGILUkNVGWWBsOuub2ScF/cfj+1Ty0lJhoqToihZ01DlFKdZljiJiCOdfH/3UGR7OOakKF6oOCmKkjV2/zxwWk4ArfWJzf0XNdck7S6hKCpOiqJkTVuDU3zsIYPgnIgbRV16ykSoOCmKkjWz3eJU77Scls5K7J135vyGvJ6TMrVRcVIUJWsWNjstI3csaYlHY9ezFujkWyU5Kk6KomTN+UuaidbSVpcHeINr5PoKlwuvvLSElXPUclKSo+KkKErWzG+q5vNXnc4bT53Jt647KzYuw96/pr0x9vztZ86hssy7s4SigE87RCiKMvW49tx2rj23Pen+O69ew1ce3kZjVRl/fdmyAp6ZMhVRcVIUpSC0N1fz9fcmHS6tKA7UracoiqL4DhUnRVEUxXeoOCmKoii+Q8VJURRF8R0qToqiKIrvUHFSFEVRfIcYY4p9DlMaEeloaGiYu3r16mKfiqIoypRh06ZN9Pb2HjDGzPPar+KUJSLyIjAT2AE0A8fSPERU1Tbl8ryUpGTyO5oK+PVzFeO88v2e+Th+Lo6Z7TEKff1aChw1xngWv6k45RAR+bYx5gNpvmYDgDFmfT7OSXGSye9oKuDXz1WM88r3e+bj+Lk4ZrbH8Nv1S2NOueVXxT4BZVJO1t+RXz9XMc4r3++Zj+Pn4pjZHsNXf0NqORUZtZwURZmq5PP6peKkKIqi+A516ymKoii+Q8VJURRF8R0qToqiKIrvUHFSFEVRfIeKk08RkUYReU5ENonIFhG5udjnpCiKkg4iUi0ie0Xkn9N9rU7C9S99wMXGmEERqQG2iMjPjDF+7AKgKIrixf8GnsnkhWo5+RRjTNAYMxh5WgFI5J+iKIrvEZFTgOXAA5m8XsUpT4jIxSLySxE5ICJGRG70WPMhEdktIsMi8ryIXOTa3ygiLwEdwFeMMV0FOn1FUaYxubh+Af8MfDLTc1Bxyh+1wBbgI8CQe6eIXA3cCXwRWAM8BTwgIu3RNcaYHmPMmcAi4FoRaS3EiSuKMu3J6volIu8Athtjtmd6AtohogCISD/wYWPMd61tzwIvG2Nutra9DvzEGJNwtyEi3wB+b4z5SQFOWVEUBcjs+iUiXwLeBwQJC10Z8C/GmH9M9X3VcioCIlIOnAU87Nr1MLAusqZVROoijxuAi4FthTxPRVEUN6lcv4wxnzTGzDfGLAT+GrgrHWECFadi0QIEgMOu7YeBtsjjBcATkZjTE8DXjTGbC3eKiqIonqRy/coaTSX3KcaYjcQHeSmKokxJbHdgOqjlVBy6CPti3QkOrUBn4U9HURQlZQpy/VJxKgLGmFHgeeBS165LCWe9KIqi+JJCXb/UrZcnRKQWWBp5WgK0i8hqoNsYsw/4KvA9EdkIPAl8EJgDfLMY56soihLFD9cvTSXPEyKyHnjMY9fdxpgbI2s+BHwCmE24puBjxpjHC3WOiqIoXvjh+qXipCiKovgOjTkpiqIovkPFSVEURfEdKk6KoiiK71BxUhRFUXyHipOiKIriO1ScFEVRFN+h4qT8//buJ0SrKozj+PeXhmkiaFbWIlpUmBrVTEhGQhmkqxpKqMWoTEstIlsEbVoFRX9gKIIWY5nhphL6R5oxBSFSkgUpFFEqRiJZWQRNU/m0OOel4eXe1zvvTHMvM78PvNzhnj/P2T1zzj33XDOzxnFyMjOzxnFyMjOzxnFyMpsGJK2Q9Lek9sM4/694d0oalXTlVMSzmcfJyaxGkhZIOiMpOvxurtDVs8C+iNjbxRhey3FKvx+m5Iik05LmRsSbwJfAk+ONZ1aFTyU3q1cPIGAn8F5JnQOdOpC0ivS5gr4uxzAErAcGgAdL6twKXA68GBF/5HuDwHZJyyPicJexzQr54FezGknaCjwDrI2I97vsYwewDrg0Iv7qov05wFFgXu5jtCRGP7AyIg7ke/NJn+beFhEPdDN2szJe1jOrVy8QnGV2VEbSbNKM6YOixCRpjqRHJR2WNJKX5d6WdH2rTkScAV4GLgDuKOhjAXA3cKiVmHK734GPSbMus0nl5GRWrx7gGDBL0uL2X4X2vcB84NP2AknnAruBx4D9wEPAE8AyYJ+kG8ZUf4mUJAcKYtwLzCUt/7XbDyyRtLTCWM0q8zMns5rkZbGrSP8k/lhQ5QTp66KdLMvXbwvK7gduAdZFxJ4xcV8gfRzu6VxORByR9CGwVtIlEXFiTD8DwCjwakGMVtzlwFdnGatZZU5OZvW5jpSYBoF3Csp/qdDHhfn6c0FZPylhfFYwC9sLbMo771obHIaANcBG8i68PCO6EXg9Ik4VxPgpXy+qMFazypyczOrTm69vRcRwl320djSpoOxq0nJc0aysZTFwPP+9CzhNmim1tojfl6/bStq34npnlU0qJyez+vTk60S2YbcSz6KCMpHeRdpaoT0RMSJpJ7BZ0k3AJ8AG4HtgT0n7VtxOCdBs3JyczOrTC5yKiJMT6ONQvhad1PANadlvOO/Iq2II2EyaPS0ClgCPd2h/Rds4zCaFd+uZ1UDSPGApE5s1AXwO/EZ6LtTuFVJyKZw5Sbq4/V5EHAS+AO4BtpCW68qW9MhxT0bE1+MbtllnnjmZ1eNaYBaApP6SOu9GRMdNERHxj6RdQJ+kORHx55jiQdLJEU9JWgMMkxLZZcBtwAjp5Id2Q8BzpBd7P4qI74pi592Gq+mcvMy64hMizGogaQvwfIcqASyMiF8r9LWS9HxofUS80VY2m7RMt4H/tp3/QHovanvRqRSSFuY65wEbI2JHSdxNpJd3r4kIL+vZpHJyMpsGJO0Gzo+I1VMY8yBwNCLumqqYNnP4mZPZ9PAwsErS7VMRTFIfsAJ4ZCri2czjmZOZmTWOZ05mZtY4Tk5mZtY4Tk5mZtY4Tk5mZtY4Tk5mZtY4Tk5mZtY4Tk5mZtY4Tk5mZtY4/wLgzgc+i+jB9gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = s.default_plot(plot_kwargs = {\"lw\": 3}, mode=\"survival\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the radial profile at the highest energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEBCAYAAACaHMnBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxc1Xnw8d+j0WattqzFq7zbGNtg4wVsNuPgLCwhJG8CScpSSpoAIaG8bUOahAQCtGkaUpqEkNCkQNISmqQNAcKLTcAGA8ZgEMbYeMGWbXnR4kWWrF163j/uzOjOaEaakWaVnu/no4905557dUZjzzPnnOecI6qKMcaYkS0j2RUwxhiTfBYMjDHGWDAwxhhjwcAYYwwWDIwxxgCZya5AtETkbaAM2J3suhhjTBqZCdSr6qJQJ9MuGABlxcXFExcuXDgx2RUxxph0UVVVRWNjY9jz6RgMdi9cuHDiunXrkl0PY4xJGytXrmT9+vVhe1RszMAYY4wFA2OMMRYMjDHGEGEwEJELROSPInJQRFREro/gmgUisl5EWr3X3SkiElTmUyKyTUTavd+vHOTzMMYYMwSRtgwKgK3AV4HWgQqLSBGwFqgFlnqv+zvgdleZ5cATwH8CC73ffysiZ0dRf2OMMTEQUTaRqv4J+BOAiDwSwSWfB/KA61S1FdgqIqcBt4vI/eoslXob8KKq3uu95l4Rucj7+GejexrGGGOGIl6ppcuBl72BwOc54LvAVGCvt8yPgq57DvhyqBuKyDrvjwtjWdFEUlU6unto7+qhvbOHts5u2rt6v7d3ddPe6Xxv8373n+/sCSjrvqYt4Brn52xPBmWFOZQX5lJWmNP7VeB8Ly/KoTAnk6CeO2PMCBWvYDAOqAl6rNZ1bq/3e22IMuPiVCe/7h4NfMPt7KGtK5o33N437aiu6eohkdtHvH+kqd/zOZkZfYJEqABSWpBNTqYnQbU2xiRD2kw6U9WV4G8hXBjt9YcbW1l9/0u0d3XT2W0b+gC0d/VQc7yVmuMDDgNRPCrLHzTKiwKDR1lhDlNK8qkcm5eAWhtj4iFeweAIUBH0WIXrXH9ljhAHWZ4Mmtu74nHrKOsh5GZ6yMnKIMf1PTcrg5xM98+933MyM8jNcn33lvU9Fuo+rR091De3Ud/UTn1TO3Xe7/VN7dQ3t1N3sp3Wzu6I693Y2kljaye765rDljl/Vim3rprFsmklsfhTGWMSKF7B4DXgeyKSq6pt3sdWA4eAaleZ1cD3XdetBl6NR4Vys3q7OUTwvyH3vjEHveGGfJMNOu+6PneAsrlZHrIzM/BkJLKPvrjfs6fau4ICRRv1zX2DxtFTHXT3DNyaenlXAy/vauCc6SV8ZdUsls8Ya2MSxqSJiIKBiBTgrHgHTjpqpYgsBI6p6n4R+Udgmap+yFvmv4BvA4+IyD3AbOAO4C7t3XT5AeAlEbkD+ANwJXARcF4Mnlcf+dke3v3Oh8nJ9JDlEXuTAvJzMsnPyWRqaX6/5bp7lOMtHb1BwhsofIGktrGNN/cdwxcvNu45xsY9r7N4yhi+8qFZXDCr1P7exqS4SFsGS4AXXcd3eb8eBa4HxgMzfCdVtVFEVgM/Ad4EjgM/AO53lXlVRK4G7gHuBj4ArlLV1wf7ZPojIhTmZsXj1sOeJ0MoLcihtCCHueNDl6luOMWD63bzP28dpMsbFTbvO851v9zEmZOK+cqHZrHqtHILCsakKNFEprfEgIisu/DCCy+0VUtT04FjLfx0/Qf89s0DfQbq500o4tZVM/nw6ePISGh3mTHGu2rpel8yTjBbm8jE1OSSPO67cgHr/+4irls+hezM3n9i7x06yZd+/RYfe+BlnnrnUETjEANp6+xm++GTHG1uH/K9jBnJ0ia11KSXCaNHcdcV87nlopn8/KU9/Pr1fbR19gCwo7aJWx9/m399fidfXjWTy8+YQKZn4M8lqsrBE628vf8Eb+0/ztv7T7Dt0Ek6unvIEDh3ZikfP3MCH5k/jiLrEjQmKtZNZBKiobmdf395L4+9Vk1LR2BK69Sxedx80UyuXDSRLFdQaO3o5t2Djd43fufNv65p4BZAdmYGq+aUc8XCCVx0WnlAJpkxI9VA3UQWDExCHT/VwS9f2csjr1TTFDTvY9KYUXzu7EqONLbx9v4TbD980j8Y3Z8JxbkcamwLea4wJ5OPzB/HFQsnsHz62IhaIMYMRxYMTEpqbO3kkVeq+eUre2ls7Yz4uoKcTBZOHs2iytGcVTmGhZNHMyY/m0MnWnl6yyGerDrEe4dOhry2tCCHy84Yz8cXTmDR5NGW2WRGFAsGJqU1tXXy2Gv7+MWGvRw71dHn/MzyAs6qHM2iyjGcVTmGmeUFA07c213XxB+rDvHkO4fYd7QlZJl5E4p4+NolTBg9KibPw5hUZ8HApIWWji4e33SArQcbmTI2j0XeT/3FowY/EKyqbKlp5MmqQzy15RD1QeMN580s5Vd/tcxaCGZEGCgYWDaRSQl52Zn81XnTYnpPEeHMyaM5c/JovnHpXDbuOcof3j7I79+qoUdhw+4G/uetg3xq8aSY/l5j0pGNppkRwZMhnDuzlO9/+kz+8tzeoPPdZ7bRYHMUUsqp9i7u+P0WLr5/PS/uqEt2dUYMCwZmxLl99WwmescKTrR0cvdT25JcI+Nz8EQrn/rpq/zmjQPsrmu21yaBLBiYESc/J5N7r5zvP/7jO4d44f3gfZZMom3ed5wrfrwhYFOmvQ2nONw48H4bZugsGJgRaeWccq5cNNF//M3/3ZoS+12MVH94+yCffXgjDc19M8pe++BoEmo08lgwMCPWty47nZL8bAAONbbxL8/tSHKNRp6eHuVfntvBbU9U0dHlLFcyJi+Lj87r3f3WgkFiWDAwI1ZJfjZ3Xna6//jR16p5a//x5FVohGnp6OKW/3qLH7+42//YrPICnrzlPG48v3eQ/7U9FgwSwYKBGdGuWDiBC2aXAaAKd/x+i/8TqomfI41tfOZnr/Hs1t5dbi+cXcbvb15B5dg8zpg0mlHeNaVqjrdy4FjoyYMmdiwYmBFNRLj3E/P9bzw7a5t5aP0HSa7V8NbR1cPn/30jWw/2Lhvyl+dO5RfXLfGvNpudmcGSqWP85611EH8WDMyIN7kkj7/9yBz/8Y9f2M3uuqZ+rjBD8Wb1MT6oPwU48z/uvXI+3758Xp9FBJfPGOv/eaMFg7izYGAMcP2KqZw5qRiAju4evv4/79ITg813TF9vHzjh//nTiyfx+bOnhCy3fLorGHxwlHRbOifdWDAwBucT6j996gwyvYvgvVF9nP/atD/JtRqe3nEFg0WVo8OWmz+xmPxsp/vuUGMb+23cIK4sGBjjNXd8EV+6cIb/+J+efZ8jYfZJMIOjqlS5gsGZk8MHgyxPBkunlfiPLcU0viwYGOPy5VUzmV6aD0Bzexff/MNW656IocONbf7d6vKzPcwqL+y3vLuryAaR48uCgTEuuVke7vvkAv/x89trA9IfzdC4WwULJhUPuDeFexD5NRs3iCsLBsYEOWf6WD67rNJ/fOeT73GyLfLd2Ex47vGChZPH9FPSMW9CMYW5zkr7dU3t7G04Fbe6jXQWDIwJ4Y6PnUZ5YQ4ADc3t/PiF3QNcYSLxdkAwKB6wvCdDONs9bmBdRXFjwcCYEIpHZfGNS+f6j//jlb32qXSIurp7eLem0X8cScsAnJaajw0ix48FA2PC+PiZE1gyxXnD6uxW7n3G1tYfil11zbR2dgMwriiXccW5EV0XOPnsmI0bxIkFA2PCEBHuvLx3Ibvnt9fx0s76JNYovVUFdBGFTykNNndckX8v7IbmdnbXNce8bsaCgTH9OmPSaD7t2iP5u09vo7PbFrIbjKr9kc0vCJZh4wYJEXEwEJGbRWSviLSJyGYROX+A8reIyHYRaRWRHSJybYgyXxWR971lakTkJyJSMJgnYky8/N1H5/hnwu6qa+Y/N+5Lco3S0zs1g2sZQN8UUxN7EQUDEbkKeAC4D1gEvAo8KyKVYcrfBHwPuBuYB3wb+ImIXO4q8zngn4F7gbnAtcAl3t9jTMooL8zly6tm+Y9/+Pwujp/quyOXCe9Uexc7a53F/zIEzpg0cCaRmzsYvL73mK0bFQeRtgxuBx5R1YdVdbuq3gocBm4KU/4a4GFVfVxV96jqb4CfA19zlVkBbFTVX6lqtaq+ADwGnD24p2JM/Nxw3lSmjM0DoLG1kx8+vzPJNUovW2oa8b1/zyovJD8nM6rrZ5cX+nelO3aqg522qmzMDRgMRCQbWAysCTq1BucNPZQcIHhRl1ZgmYhkeY83AAtF5Bzv76kEPg78KUw91onIOmDhQHU2JtZyMj38wyW9qaa/3riPHUfsDSlSgx089snIEM6ZbusUxVMkLYNSwAPUBj1eC4zrWxyA54AbRGSpOJYANwJZ3vvhbS38A/CSiHQC+4B3CWw9GJMyPnx6BefOdLorehTufvo9S3OMUMDM435WKu3PcptvEFfxyib6LvAMzthCJ/Ak8Kj3XA+AiFwIfAu4GTgL+CSwErgr1A1VdaWqrgSq4lRnY/olItx52Tx8y+m8svsoa7cFf0YyoQy1ZQCBk89s3CD2IgkGDUA3UBH0eAUQcgUvVW1V1RuAPGAqUAlUA02AL1H7HuBxVf13VX1XVf8Xp6Xw9yISXYeiMQkyZ1xhwGYs9/5pO+1d3UmsUeo70tjGkZNOr/GoLA+zygeXMDizvIDSAmeJkMbWTrYdPjnAFSYaAwYDVe0ANgOrg06txvnk39+1napao6rdwNXA06rqS9LOwwkybt1A/8sYGpNkt6+e7Z8Ete9oC//xSnVyK5TiglcqDd7eMlIigeMGthVmbEX6qtwPXC8iN4rIXBF5AJgAPAQgIo+JyGO+wiIyW0SuEZFZIrJMRH4DzMf55O/zFPDXInK1iEwTkdU43UtPq2pXLJ6cMfEwJj+bv7m4N9X0xy/spq7JNsEJxx0MFg2yi8jH9kWOn4iCgao+AdwGfBOnz/484BJV9c2+qfR++Xhw0lHfAdYCucAKVa12lbkH+AFOANgG/BInQ+nGQT4XYxLm8+dMYaa3u6O5vYt/eW5HkmuUuqoOHPf/PNjxAp/lQeMG3TZuEDMRt9dU9UFVnaqqOaq6WFVfcp3zDe76jrer6iJVzVPVYlX9hKruCLpfl6repaqzVHWUqk5W1ZtV9TjGpLgsTwbfuqx33aLfbq4JWJHTOLp7NODvEs0yFKFMK82nosgZN2hq6+K9Q/Y3jxVbm8iYQbpwdhkfOq0cAFW46ylLNQ22u66ZUx3O0GB5YQ7jI1ypNBwRsRTTOLFgYMwQfOPSuWR5nJyHN/cd56kth5Nco9QS3EUkMvT8kIB1imzcIGYsGBgzBNPLCvjLc6f5j//pT9tp7RgZqaaNrZ1c9bPX+PiPN7CnPvSy0u7B46F2Efm45xu8WX2cLltFNiYsGBgzRF9eNZOx3nVzDjW2cffT20ZEd9EvNuzl9b3H2FLTyG1PVIUczK060NunP9RMIp/Kkjz/uEFzexfv27IgMWHBwJghKsrN4u8/Osd//Pim/SNi7sEL7/fOvt5S08ijr1YHnG/p6GLHEWdimIgzxyAWRISlU3vnG7xRfSwm9x3pLBgYEwOfWTKZj585wX98zzPbePH9uiTWKL7qTrax9WDgDOAfrNnBoROt/uOXdta7ViotoDA3i1hZNs2CQaxZMDAmBkSEf/4/Z7DIuwhbj8Ktj789bFc2fXFH30B3qqObO590Mqpe3FHHV3/Tu4zY4iklfcoPxRLX/TbtPT4iuuXizYKBMTGSm+Xh59csYeLoUYDTn33DI2/Q0Nye5JrF3p+39waDj83vXbz4+e21fOvJrXzh0Tdp73IGdssKc7h55YyY/v454wopzHWWMGtobmff0ZaY3n8ksmBgTAyVFebwi+uX+LfJPHiilb9+7E3aOodPhlF7Vzcbdjf4j//vh+fw2WWT/ce/3rifLm//0MTRo/jtF5czuSQvpnXwZAhLpozxH2+yrqIhs2BgTIydNq6IH31ukX+p67f2n+CO328ZNl0Zm/Yeo8WbPltZkseMsnzu+Ohc/4qiPjPK8vndTcuZWpofl3osdY8b7LVgMFQWDIyJg1WnVfCNS3uXq/hD1SF+/MLuJNYodl5wDYyvOq0cEaE4L4vvfLz3+c6fWMR/f3E544tHxa0ellEUW7ZvgDFxcsO5U9ld18zjm/YD8IO1O5leVsClZ4xPcs2G5sWgYOBz2RkTyMzIoOZ4C1cvq6Qgyn2Oo3XGpGKyMzPo6Oqh+mgLdU1tlBcObbmLkcxaBsbEiYhw9xXzWOFaPuH2/64K2AIy3eypb6baO1ibl+3h7OmBWUIfnT+OG8+fHvdAAM6+1Asn9U5ke7Pa1rgcCgsGxsRRlieDn35+MdO9/ebtXT3c+NibAfn46cTdRXTuzFJyMj1JrA0sneYaRLZxgyGxYGBMnBXnZfGL65f6d0erb2rnxkff5FR7+u3h5J5f4O4iSpYlNm4QMxYMjEmAaaX5PPQXi8n0phhtO3yS256oSqtN3ZvbuwI+fV80J/nBYPGUMfgWQt1++CRNbZ3JrVAas2BgTIIsnzGWez4x33+8dlst33vu/STWKDobdtXT2e0Er9PHFzFuiHsTxEJRbhZzxxUBzqzvt/an73hMslkwMCaBrl5WyRfO713y+mfr9/Dfbx5IYo0iF5xSmiqW2XyDmLBgYEyC3fGxuVw8t/fN9Bv/+y5bD6b29o09PcqLO+r9xxelUDBYMtVmIseCBQNjEsyTITxw9SJOG1cIQGe39ln+OdW8d+gk9U3OGksl+dlD3tg+lpa5BpGrDpygvWv4LP2RSBYMjEmC/JxM/vGTC/zHz713hI6u1N2xy91FdOHsMjwZQ9++MlbKi3KZMtZZ+6ijqyflW1mpyoKBMUmycPJo/wqnJ9u6eOWDhgGuSJ4XXCmlqdRF5ONemmLTXpt8NhgWDIxJEhHhkgW9yz8/s+VwEmsTXkNzO1tqnCwdT4Zw4ayyJNeor6WucQObbzA4FgyMSaJLz+jdHW1NinYVrdtRj2/B1cWVYyjOi92OZbHibhm8WX0sreZvpAoLBsYk0ZmTigO7inanXleRe6/jVOwiAmdSX2lBNuD8HXfWDc8d5uLJgoExSSQiAauYPvNuanUVdXb38PLO3gCVSvML3EQkcElrm28QNQsGxiTZJQt6g0GqdRW9WX2cJu8aShNHj2J2RUGSaxSee52iTbaCadQsGBiTZKncVbRuZ28W0co5ZYikTkppsGVBLYPhsrNcokQcDETkZhHZKyJtIrJZRM4foPwtIrJdRFpFZIeIXBuiTJGI/JuIHBKRdhHZLSKfGcwTMSZdpXJX0XrXrOOVKbAwXX/mji/07z195GQbNcfTc5nwZIkoGIjIVcADwH3AIuBV4FkRqQxT/ibge8DdwDzg28BPRORyV5ksYC0wC/gMMAe4Htg7yOdiTNpKxa6iI41tvH/EGYjN8kjAJj2pKNOTwVlTLMV0sCJtGdwOPKKqD6vqdlW9FTgM3BSm/DXAw6r6uKruUdXfAD8HvuYq85dAGXCFqm5Q1Wrv9zcG+VyMSVup2FW03tVFtHRqCfkJ2L1sqJbZ/gaDNmAwEJFsYDGwJujUGmBFmMtygLagx1qBZd4WAcAngFeAH4nIERHZJiLfcZ0Prsc6EVkHLByozsakm+CuoqdTYALauoAuotSbaBZKwCCyZRRFJZKWQSngAWqDHq8FxvUtDsBzwA0islQcS4AbgSzv/QCmA5/2PnYp8C3gS8A/RvUMjBkmLnV3FW1LbldRZ3cPG3b1tk5SfbzAZ1HlaLI8ziD3B/WnONrcnuQapY94ZRN9F3gGZ2yhE3gSeNR7zvcvPAOoA76gqptV9ffAncBNEiJlQVVXqupKoCpOdTYmqc6YVMykMU5XUVNbFxt21w9wRfy8ta83pXRCcS6zylM3pdQtN8vDgonF/uM391mKaaQiCQYNQDdQEfR4BXAk1AWq2qqqNwB5wFSgEqgGmgDfv/DDwE5Vda83u917TSnGjDAiEtA6eGZLyP9eCbF+Z28gunBOeUqnlAZbapvdDMqAwUBVO4DNwOqgU6txPvn3d22nqtZ43/CvBp5WVV/L4BVgpoi46zAbaMEJQMaMOJcEdRUla23+dBwv8Fk6xQaRByPSbqL7getF5EYRmSsiDwATgIcAROQxEXnMV1hEZovINSIyS0SWichvgPnAP7ju+VOgBHhAROaIyEeAu4AH1WaLmBEquKsoGVlFdSfb2Hb4JACZGcK5M9Oroe7e+WzroZO0dHQlsTbpI6JgoKpPALcB38Tpsz8PuERV93mLVHq/fDw46ajv4MwlyAVWqGq1654HgA/jZCpV4QSWXwLfGPzTMSa9BXcVJSOraJ2ri2jJ1DEUpEFKqdvovGzmVDi7yHX3KG/vP5HkGqWHiF9lVX0QeDDMuZVBx9txJqcNdM+NhE9PNWZEumTBeH720h4A1m6rpb2rm5xMT8J+fzrNOg5n6bQx7Kh1Jsxt2nss7Vo3yWBrExmTYvpkFe1KXFdRV3cPL+9K3/ECH/cKpq+m8A5yqcSCgTEppk9WUQLXKnr7wAlOtjl97OOKcv3dLenm3Jml+BKgNu87zrFTHcmtUBqwYGBMCnLPRvZ1FSXC+qAsonRKKXUrLcjhrEpnILlH4YX36wa4wlgwMCYFLZiYnK6i4CWr09nFc3unRj2/LXgBBRPMgoExKSgZXUV1TW1sPdibUroizQddV5/eO/j90q562jqTM2cjXVgwMCZFBXQVvRf/rqKXXNtbnjVlDEW5qbfxfTRmlBUwrTQfgJaObl7bczTJNUptFgyMSVEBXUXt8e8qWrdj+HQRgdO6unhub+vAuor6Z8HAmBQVvKz1fX/azo/+vIvN+47R2R3bFU2dlFLXKqWz03N+QbCAcYPttbYVZj/Sa2qhMSPMpQvG87P1zgS0D+pP8YO1O/nBWsjL9rBsWgnLp49lxYxSTp9QhCdj8Jk/79Q00tjaCUBFUQ5zx6dnSmmwxVPGMDovixMtndSebGfrwZMsmFQ88IUjkLUMjElhCyYW88mzJvZ5vKWjm3U76vnHZ9/n8h9vYNHda/jJi7sH/XvWu7qILpydvimlwTI9GaxyzaJeuy15K8GmOgsGxqQwEeH+zyxk49c/xA+vOpNPL57k3x7T7WRbF99/bgdVBwa3Do97PaILh0kXkc/Fp/d2Fa3dbvMNwrFuImPSwLjiXK5cNIkrF00C4MCxFl79oIFXPzjKhl0NHPXOsP3h2p08esOyqO7d0NzOlppGADwZwnmz0julNNgFs8vI9mTQ0d3D9sMnqTnewqQxecmuVsqxloExaWhySR5XLa3kgasX8cQXl+MbLli/s57NUe7u9ZKrVXBW5WiKR6V3SmmwgpxMzpkx1n/8Z2sdhGTBwJg0N7O8gI+fOcF//K/P74zq+nXDYJXSgaw+PTCryPRlwcCYYeArH5rlbx28vKsh4h2+uns0YJXSC2en//yCUNzzDTbuOcrJts4k1iY1WTAwZhiYXlbAJxb1Zh39cG1krYMtNSc43uK8MZYV5jBvQlFc6pds44tHMX+i89w6uzWga8w4LBgYM0x8ZdUs/1yDVz84ysYIll9wdxENp5TSUKJZuE5VR9wENQsGxgwTU0vz+ZRrTsL9a3f2+4a27+gpfvnKXv/xcO0i8nEHgxferws7i/ulnfWc8Z01rP7hSzS2jJzuJAsGxgwjt66aRaa3dbBp7zFe+yB066Cts5ubfv0WTd6NbCaOHhXwZjkczZtQxPjiXMCZlxFqXKWhuZ3bnqiiqb2L3XXNI2qw2YKBMcPI5JI8Pr1kkv/4h8+Hbh3c9dR7bDvsLFed5REe/PxZjMpO3D7LyeAsXOfuKgpMMVVV/uF/3g3YFe3IybaE1S/ZLBgYM8zcctFMsjxO6+CN6uNs2B242unvNtfw+KYD/uM7LzudMyePTmgdkyVwNvKRgED5h6qDrAkaS6i1YGCMSVeTxuTxmSWT/cc/dI0dbD98km/+4V3/uSsWTuAvzpmS8DomyznTSyjIcRZeOHCslV11zQAcbmzlziff61P+SKMFA2NMGrvloplke5z/3m/tP8H6nfWcbOvk5v98i7ZOZ+B0ZnkB9125YFhnEAXLyfQEDJSv3eYsa/2137/rHz/Jzep9W6xtak94HZPFgoExw9CE0aO4ellg6+Brv9vC3oZTgLME9kN/cRb5OSNvebKLXdthPr+9lsc3HfDPOxCBez6xwH++zrqJjDHp7uaVM8nOdP6Lv1PTyLNbe5dv/qdPncHM8uGxZ0G0LppT7p+PUXXgBPc8s81/7q/OncZlrg2F6pra6ekZGfMNLBgYM0yNK87lc8sq+zx+7fIpAWsZjTSj87JZMmUMAKrO3hDgdJv97UfmkJvlYXSes1hfd4/6V4RNlu6exEyAs2BgzDB288oZ5GT2/jc/c/JovnHp3CTWKDW4F64DZ+nuH3z6THKznPTaisJc/7lkZhRt3HOUhXetYck9z3Pfn7azp745br/LgoExw1h5US63r54NOBPLfvK5ReRkDu/5BJH4UNAEu1tWzghIry0vyvH/XNeUvGDw+Kb9NLV3cfRUBz9/aQ+rfrCeq372Gi/uiP0y3BEHAxG5WUT2ikibiGwWkfMHKH+LiGwXkVYR2SEi1/ZT9rMioiLydDSVN8YM7IsXzuDlv7+INX9zgW3q4jWtNJ+PzHMCwtnTSvjyqlkB5yuK3C2D5GUUHT7RNxC9vvcY+7yJALEUUSqBiFwFPADcDGzwfn9WRE5X1f0hyt8EfA/4AvA6sAx4WESOq+pTQWWnA98HXh7KEzHGhDe5xIJAsB999iyqj55iRlmBf0DZp8LVMkhmN5G7VbJsWgmb9x0nM0P8O97FUqR5ZbcDj6jqw97jW0Xko8BNwNdDlL8GeFhVH/ce7xGRpcDXAH8wEJEs4HHgG8BFwPDab88Yk7KyMzOYXRE6oyoVWgaqGvC7H752CW2d3bxz4ATFebHfjW7AbiIRyQYWA2uCTq0BVoS5LAcIDqetwDJvAPC5F6hW1UcjqMc6EVkHLHj5A1AAABHkSURBVByorDHGDEW5awA5WXMNmtu7aO10Mp1yMjMoys2koiiXD88bF5ffF8mYQSngAYKX76sFwtXqOeAGEVkqjiXAjUCW936IyIeBzwBfHEzFjTEmXsYV9waDZC1WV+ea/VxRlBv3meLxyib6LvAM8CrQCTwJ+D7994hIGfAIcJ2qnojkhqq6UlVXAlUxr60xxrgEjhkkp5uozvV7ywtz+ikZG5EEgwagGwhe7LwCONK3OKhqq6reAOQBU4FKoBpoAuqBecB44M8i0iUiXcC1wCXe4znRPxVjjImN0oIcfB/Ej55qD7sRTjy5B4/dqa7xMmAwUNUOYDOwOujUapxP/v1d26mqNaraDVwNPK2qPcAbwAKc/n/f1x9xMooWAnvD3NIYY+Iuy5PB2HznDVjV2fQm0QJbBrn9lIyNSLOJ7gd+JSKbgFeALwETgIcAROQxAFW91ns8Gzgb2AiMwclGmg9c5y13Ctjq/gUicgLIVNWAx40xJhkqinL8QaD2ZDvji0cl9PcnumUQUTBQ1SdEZCzwTZzuna3AJaq6z1skeAEUD04AmIMzZvAisEJVq2NRaWOMibeKolzeO+TsBpeMuQa1KdoyQFUfBB4Mc25l0PF2YFE0FVHV66Mpb4wx8eQeRE5Geqm7ZVCRCmMGxhgzEpUXJnfimTu1NBEtAwsGxhgTQuAs5CS0DFIwtdQYY0acgLkGCd7+sqWji+Z2ZxvObE+Gf3+FeLJgYIwxIbhbBokeM3C3CsoKcxKyT7UFA2OMCcGdzpnoJSnc3VKJSCsFCwbGGBNSaX6Of2nrEy2dtHkXjUuEwMFjCwbGGJM0GRkS8EZcn8Bxg+BF6hLBgoExxoRRnqSMIvcYhbUMjDEmySoKk7N6aaLnGIAFA2OMCStZcw0SvS4RWDAwxpiwAucaJC4YJHpdIrBgYIwxYZUHzDVIYDeRpZYaY0zqGGw3UWNrJ7/euI+tBxuj/p1tnd2cbHNmH2dmCCV52VHfYzAiXrXUGGNGmsDtLyMPBt/7f+/zX6/vRwQevmYJF58evFFkeO4U1rLCHDIy4j/7GKxlYIwxYVUUDq6b6I29xwBnl7Sv/OZttnn3RYhEbRLSSsGCgTHGhDU6L4tsj/M22dTexSnv4nEDOXaqw/9zS0c3Nz76RkCGUH8C0koTNOEMLBgYY0xYIhIwgBtJV1F3j3KspSPgsUONbfz1Y5sjWtIiGRPOwIKBMcb0K3AQeeCuohMtHag6P2d5BF+Xf9WBE/ztb99BfSfDqE3ChDOwYGCMMf0a504vjaCr56iri2hySR53Xna6//jpLYf51+d39Xt9wKY2CUorBQsGxhjTr2i7iY429waD0vwcrlsxlWvOmeJ/7Ecv7OJwY2vY6xO997GPBQNjjOlHtN1ER0/1linJz0ZE+PblpzNvQhEAPQpbasLPP6hLwuxjsGBgjDH9inaugTuTaGyBM2Es05PBsmkl/sd31TaFvT5gXSIbQDbGmNQQ7VyDBlc30dj83tnDsysK/T/vqmsOeW17VzfHWzoByBAYW2DBwBhjUkLAngYRDCAfc3UTud/MZ5UX+H/eWRs6GLhnH5cW9O60lggWDIwxph/B3UQDpYa6B5BLXC2DWeW9LYMP6pvp7ul7n8AJZ4lrFYAFA2OM6VdBTiZ52R4A2jp7/IvIheMOBr4xA4DivCz/GEBHVw/7j7X0udbdDVWRwMFjsGBgjDH9EpGoVi91ZxONzQ/8dB8wbhBiEDkZm9r4WDAwxpgBlBdGnlF0NEQ2kc9M17hBqEFkd8ugLFVbBiJys4jsFZE2EdksIucPUP4WEdkuIq0iskNErg06/wUReVlEjovICRF5UUTOG+wTMcaYeIl0rkFXdw8nvNlAIjAmaC+CaFoGiZxwBhEGAxG5CngAuA9YBLwKPCsilWHK3wR8D7gbmAd8G/iJiFzuKrYSeAJYBZwN7ACeE5FZg3omxhgTJ5HONXAvUDcmL7tPNtCsiv4zipKx3aVPpJvb3A48oqoPe49vFZGPAjcBXw9R/hrgYVV93Hu8R0SWAl8DngJQ1c+7L/AGkE8AHwX6X7zDGGMSqCJg+8t+gsGp0JlEPu70Ul9GkTtgBGQTJXDCGUTQMhCRbGAxsCbo1BpgRZjLcoDgv1grsExEssJckw3kAsfD1GOdiKwDFg5UZ2OMiaVIu4mOhplw5jM6L5sy75t8e1cPB4IyiupTfAC5FPAAtUGP1wLjwlzzHHCDiCwVxxLgRiDLe79Q7gGagT9GUCdjjEmYiggnnvU3eOwzuyL0IHJnd4//ehFn0lkixSub6LvAMzhjC53Ak8Cj3nM9wYVF5KvAF4FPqmrI/eFUdaWqrgSq4lFhY4wJxz1m0N+SFEebw6eV+rgnn+10DSI3NLf790EYm59NliexyZ6R/LYGoBsI3tG5AjgS6gJVbVXVG4A8YCpQCVQDTUC9u6yI3IbTKrhEVTdFUXdjjEkI92BuXVMbPSFmD8PAYwYQOIi829UySGZaKUQQDFS1A9gMrA46tRrnk39/13aqao2qdgNXA0+rqr9lICK347QiLlXVDdFW3hhjEmFUtoeiXCffprNbOR60raWPe5G60jDdROFaBu7B40SnlULk2UT3A78SkU3AK8CXgAnAQwAi8hiAql7rPZ6Nky66ERiDk400H7jOd0MR+TvgXuAvgJ0i4ht/aFXV8It9G2NMElQU5XKyzfkkX3uyPeSKoscC9jII100U2DLwZRTVJmnvY5+IOqVU9QngNuCbOH325+F06+zzFqn0fvl4cALAO8BanCyhFapa7SpzC86A8hPAYdfXA4N8LsYYEzeRLEkRbl0itzH52f7B4fauHmqOOxlFdUna+9gn0pYBqvog8GCYcyuDjrfjTE7r735TI/3dxhiTbJFsfxmwsU2YMQNwWgcN3sHmXbXNTBmbH5BWmoxuIlubyBhjIhDJXIOG5tB7GQRzp5furGvqc8+UHEA2xhgDFe7F6kLMNejo6l3eOkNg9Khw82thpmuNot3eZSmSuWIpWDAwxpiIuFsGRxr7BgN3hlFJfjYZ/exSNru8b8sgYC+DImsZGGNMSppckuf/eU9930Xmwu1wFsosd8ugrpnO7p6ALqayBM8+BgsGxhgTkZnlBYj3w/7+Yy20dXYHnO9vU5tgJfnZ/nkIbZ09bKk5gW8e25i8LLIzE//WbMHAGGMikJvlodLbOuhRZ9VRt4DZx2HSSt3cG91s2HXU/3My0krBgoExxkTMPXt4V9B+BAGzjwfoJoLAjW5e2d3g/zkZg8dgwcAYYyIWkBIatFNZJLOP3dwzkd/a37tyv7UMjDEmxQVsWxm0h3Eks4/d3IPIXa6F76xlYIwxKS5gQ/uglsHRCGcf+7hbBm4VSViXCCwYGGNMxGaWF+CbPrAvKKPoaISzj91lQgWN8iTMMQALBsYYEzF3RpFq4H4ERyPYyyDYzBCtg2SsWAoWDIwxJiqzAsYNeruKjkWwl0Ew9xiETzJmH4MFA2OMiUrAHsbe9NL2rm6a2p11iTwZQlFu+HWJ3Ny7nvmUWcvAGGNSX+BOZU4wCN7usr91icLdC6AoN5PcLE8Mahk9CwbGGBMF96d5XzdRQFpphOMFwfeC5HURgQUDY4yJyoyy3oyi/cdaaO3oDkwrjXC8AKC0ICdgsDlZcwzAgoExxkQlN8vDlLH5gJNR9EF9c0BaaSSzj93cGUXJmn0MFgyMMSZqAZPP6poi3u4yFPeAdLLSSsGCgTHGRC1wjaLmgEXqog0GlywY7/951WnlQ6/cIGUm7TcbY0yaClijqLYpoN8/ktnHbitmlLLmby5ACJzDkGgWDIwxJkoBS1nXNTOzrLelEOnsY7dQk88SzYKBMcZEaXpZPhnibHKz/1gLo1xzAyKdfZxqbMzAGGOiFJxRtMO1gulgWgapwIKBMcYMgnsJau3djiDqMYNUYcHAGGMGIVQ/f5ZHKMpNz953CwbGGDMIoRaZK8nPRiSydYlSjQUDY4wZhOBF5iD62cepJOJgICI3i8heEWkTkc0icv4A5W8Rke0i0ioiO0Tk2hBlPiUi20Sk3fv9ysE8CWOMSTRfRpFbumYSQYTBQESuAh4A7gMWAa8Cz4pIZZjyNwHfA+4G5gHfBn4iIpe7yiwHngD+E1jo/f5bETl70M/GGGMSJDfLw1RvRpFPumYSQeQtg9uBR1T1YVXdrqq3AoeBm8KUvwZ4WFUfV9U9qvob4OfA11xlbgNeVNV7vfe8F1jnfdwYY1Je8LjB2OHcTSQi2cBiYE3QqTXAijCX5QBtQY+1AstExLcF0PIQ93wu3D1FZJ2IrMNpRRhjTNIFjxtEs3x1qomkZVAKeIDaoMdrgXFhrnkOuEFElopjCXAjkOW9H95ro7mnMcaklL4tg+EdDAbju8AzOGMLncCTwKPecz2DuaGqrlTVlUBVLCpojDFDFTzXIF0nnEFkwaAB6AYqgh6vAI6EukBVW1X1BiAPmApUAtVAE1DvLXYkmnsaY0yqCc4oGtYDyKraAWwGVgedWo3zyb+/aztVtUZVu4GrgadV1dcyeG0w9zTGmFSRk+lh7vgiADwZwuSSUUmu0eBFOm/6fuBXIrIJeAX4EjABeAhARB4DUNVrvcezgbOBjcAYnGyk+cB1rns+ALwkIncAfwCuBC4CzhvaUzLGmMS5+4p5/Nufd3Px3PKkbls5VBEFA1V9QkTGAt8ExgNbgUtUdZ+3SPB8Aw9OAJiDM2bwIrBCVatd93xVRK4G7sGZj/ABcJWqvj74p2OMMYm1eEoJj96wLNnVGLKIV1RS1QeBB8OcWxl0vB1nctpA9/wd8LtI62CMMSY+bG0iY4wxFgyMMcZYMDDGGIMFA2OMMVgwMMYYA4i6N+9MAyJSU1xcPHHhQluvzhhjIlVVVUVjY+NBVZ0U6nw6BoO3gTJg9yBvMRY4GsfrIinXX5lozwU/5ouSyV7DabB/51jeL5pr7HVz2OsW3ePp9LrNBOpVNXTav6qOqC/g5/G8LpJy/ZWJ9lzwYzh7QqxL179zLO8XzTX2utnrNpjHh9PrNhLHDJ6K83WRlOuvTLTnBvt84i3W9RrM/aK5xl43h71ug3s82YZcr7TrJjL9824AhAbNCjepzV639DScXjcLBsYYY0ZkN5ExxpggFgyMMcZYMDDGGGPBwBhjDBYMRhwR+V8ROS4ito9EmhCRySKyTkS2icgWEfl0sutk+icio0XkTRGpEpGtIvKFZNdpIJZNNMKIyEqgELhOVf9PkqtjIiAi44EKVa0SkXE4e5LPVtVTSa6aCUNEPECOqraISD7O7pBLVDWWs7tjyloGI4yqrgOakl0PEzlVPayqVd6fjwANQElya2X6o6rdqtriPcwBxPuVsiwYpBERuUBE/igiB0VEReT6EGVuFpG9ItImIptF5PwkVNW4xPJ1E5HFgEdVD8S73iNZLF4zb1fRO0AN8H1VbUhQ9QfFgkF6KcBpbn4VaA0+KSJXAQ8A9+HsQf0q8KyIVCaykqaPmLxuIlICPAb8dbwrbIb+mqnqCVU9E5gGfE5EKhJR8cGyMYM0JSLNwJdV9RHXY68DW1T1C67HdgG/U9Wvux5b6b3WxgwSbLCvm4jkAGuBh1X1V4mt9cg2lP9rrnMPAi+oasombljLYJgQkWxgMbAm6NQaYEXia2QiEcnrJiICPILzZmKBIMkifM0qRKTQ+3MxcAGwI5H1jJYFg+GjFPAAtUGP1wLjfAci8jzwW+ASEakRkeWJq6IJIZLX7VzgKuAT3lTFKhFZkMA6mkCRvGZTgJe9YwYvAz9S1XcTV8XoZSa7AiaxVPXiZNfBREdVN2Af3NKKqm6id+ObtGD/wIaPBqAbCB6kqgCOJL46JkL2uqWfYfmaWTAYJlS1A2cy0uqgU6txMh1MCrLXLf0M19fMuonSiIgU4OxjCk4grxSRhcAxVd0P3A/8SkQ2Aa8AXwImAA8lo77GYa9b+hmRr1my9+60r6j2OV0JaIivR1xlbgaqgXacTy8XJLveI/3LXrf0+xqJr5nNMzDGGGNjBsYYYywYGGOMwYKBMcYYLBgYY4zBgoExxhgsGBhjjMGCgTHGGCwYGGOMwYKBMcYYLBgYY4wB/j895Z7xl+TipQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s.domain.rcen, 1.0 - Pradial[:,-1])\n",
    "plt.semilogx()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the radial profile in every energy bin using a colormesh plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$z$ (kpc)')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAEXCAYAAAAN0FvQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9fbxsV1nn+f2tvXdVnXPPvUlIIIHQARzpFtHpKDoaGyROG7VRp1u7p2FsFWxRASdo47tNK9jazcdBBB3zgUHtCDaaVrQbbd50NCqQgGZE5R0hQGPeIbn3npeq2nuvZ/5Ya1fVqXte6uSem9Q99/l+PvtTp/Z69lprr11Vz9lrP+v3yMxwHMdxnGUgPNwdcBzHcZwOd0qO4zjO0uBOyXEcx1ka3Ck5juM4S4M7JcdxHGdpKB/uDiwzkv4SeCTwtw93XxzHWXo+F7jXzL7owVYg6ZXA1Qc45L1m9v0Ptr1lxJ3S3jyyoLryhC65EgAJcgj9NJS+exVgCEERQCHZS11Rek0HQ68CIFaB0BixFFaIWIjYAwvJPowhNIaFVE8sQLlJtalMMdm2PdFWYBWYjFCLUCebUBvtIGABQgOxAiIoghXTuizfO4cG1MRJW7EKxFKTU7UA1p0PEFogGrLUP4N0Hx47+zwO0dJwBBGLXF13PhGKYUSjGpo2j3eENuZGNDfk2j7GXZnF7TZBEAIUBbEXsDKNswWwMo0pRT7ANBmXNCbT/mFGaEFtOkELM+cF6fqVaTwtn3tRQzE2wrCZnlMRsCJgRep7LKHti3Y1sjIYc2lvk1ItW7HHVlsxbgskY1A09EJDpUhENDHQWEFrIhJooxi3Jc24IIyFmtz/7hpo++u2sezI5VZ0m0FhhMIoQqQMkUKRkAclmogmmljQRGFRmOnMesn75veHmQ9AN/4thEaoTddAcftnyorcx7D9PCbXrAXF/L0wJt/ZWWNZd4xBtGQTu53pM2NFgCBimTYru2vO5HsCTNoZ3fl3xPFw/sQPytXA059+zcq+hn9yy9bZtrWUuFPam789zsVXfmn/a9IPXVFANKxtsbZNFpZ/MBXAIioKwsoK9CpUVVCWyUm1MX3YAdqW+JjLABhfMqB6YMT4EX3q4yXjNbH5KNEO0pdv5V4YfDbSrKQvRn0MQg0Iqg2j/0Ck3IxYKdYfXbL1KBhdZlhh9D8TGNwLg/sjg880nL6qR7MC/ZPG8BEijKHaNMZroto0+ictOSug90BL71SNskMYXTZgdFHyIrGAZiX/qBeAQf9UcnyhgWIUsQBtL1DUhgma1UAsoRilH4tmRYyPp/EIdaqjHBnHPzGk96nPYqdOoarCRmNscxOLRhj005C36cddVQVFQGWZrk0eWxvX0LYQclmvgpUBtrbC6FHHGF1SMj4umoEYn4DxRRBX8nVsRRiLagPKdeidTk4FktPurUeq0w1EI/YLmpXsrCWaFTG8WNTHoVmBYgSrd8OxuxtWP3EK3X8ayhJbG9CeWKFZLYmVGF9UcPqxgfUnj3ni4+7iGZe/n0eU63xk6wr+duNR3Ll5nH7R8Li1+3l0/ySPKDcYWsn99TEeqFfZaHsM24r1usedp05w8p41BndU9B5I17oYpc9LLKGtRKymP6qy6Y9q56TbCpo1qNeMdi3CSkN/bcyJlSEnBkPWqjGDoiEoMmwrNpuK+7dWOb3VZzwuaccFxO2OyaJQK2g1+UeFAFZFVEZUZEdfF2i9oDoZqDagOp0+o2pT/9q+aFah7UM7mP6TpgbKLdIx60Y5NMphRLUlxwPT758ZoTXCOBJGDRq1hFENdZOcUwjYoKJdG9AeKxldVDK6ODA+IZoBtCupfQvJuYWxKEbwqf/0S2x++mNnPavyldcM+MPfecy+dv/4m/+OP73lrJ3g0uFOyXEcZ4kwoO3+2d3H7ijiTslxHGepMOJCLudouiV3So7jOEuEARG/U3Icx3GWAMOoF5q+O5puyZ2S4zjOEmFAu4DDOZouyZ2S4zjO0rHYM6WjiTslx3GcJcIM2gVSCh3VrEMuM+Q4jrNkxAW2gyDpBZJulzSUdJukp+1h+2hJb5D0IUmtpBt3sHmOJNthGzzYdjvcKTmO4ywREWNs+2+LTvFJeibwKuA/AF8EvAt4i6SrdjmkD9wHvAx49x5VbwKPnt3MbLKa90G0C7hTchzHWToO+U7pRcCNZvZaM/ugmV0P3Ak8fydjM/uEmb3QzG4EPrtHvWZmd81uZ9Nuhzslx3GcJcIQ7QJbUpHkakk3d9t8XZJ6wFOAt88VvR34irPs6oqkT0r6tKTflzQRoj2bdj3QYR8kTbXu6gYFgQKhKrG2RUVvoomnKumsaTDoDk5PI0fjJMbZCYr2emg4Rq3Ri5Fw/wZ9O065VdE7VdBbT5dF0Vi5a0hxeoj1qyT42S8J46S7F7ZqtDGE4QiKwMrfnaB+xAqjiyvaHvRPNVQnx5SnRugzJ+nffTFxtaI4PaQ9PgBLdcRjPdRGwvow6X5VBWFzDKc3oK6hjazeucLK8dUkbGpGXO1hRUhCrb2CYtgQ6vT/m4YNMsOqAmJErRFXKprVKrXTGG2/YHhpRX0saejFnth6hGirFY4PHsnKp3pYGSAEwvoW9sBJVBSp/brGhqPUVlXBoJ+uQdtCE9JXtU1aeGYRtoYwGqHT6wzuvZ9BVUK/T1ztp/F6RMXooqSF14luKiY9NQtJMy50Uof5EhajlvLkiH6MWAhQiHa1YnC8ol4rqI8lw2rdKIb5/9q6xk6egvsDZb9H2etBWbKy0mfl3uMUwz4fO/lY3vQFga+54oNc1f8Ml1Qb3LN6giDj7w/u4vP6d/AFVcNaGDCymk+3Iz5eX8zHxpfzd6NL+ED5aD5Ql9QbBWoCxUgTAdZmFeq1tLUDIw4i1ouoNBSS/lxZtRRFZLVfs9YfcVFvmIRgiwZIAqzjWHKq7nNqNOCBzRWGWz3akz3KBwqqTVhdT3qGynp6sYK2B7Gf9ff6RuwbcSVSrjYcP77FpasbrFVjyhBZr3vct3mM0xsD1k/2CZsFoc6fk0GL1mp6qzUnVob0ypZhU7Kx1ef0Z1fo31MyuFf0H4BiLEInoFwFmkGgPibaXtb/UxZvbZLuYqiTxh4i6SIeF/UajC41mkfWnLh0gyuPn+Zxa/dzWW+dSi2jWHL36ASfPH0Jn/ydBj599r85RtaGXcBuAS4DCuDuuf13A199oI5t58PAvwb+CjgOfB/wTkn/0Mw+ejbtulNyHMdZMtqJovm+vNfMrj2HXdkRM7sFuKV7L+ldwHuB64EXnk3d7pQcx3GWiLR4dn+ntOCd0n1AC1w+t/9yYP4Z0IPGzFpJfwE88Wzb9WdKjuM4S4Qhagv7braI4zIbA7cB180VXUeKhjsUJAn4n0mBDGfVrt8pOY7jLBHpTmn/+4UDrJ19BfB6Se8B3gk8D3gM8GoASa8DMLNv7w6Q1GW/PQHE/H5sZh/I5T8J3Ap8NNu8kOSUZiPr9mx3N9wpOY7jLBEpCe4id0GL1mc3SboUeDFpPdH7gGeY2SezyU7rhv5y7v03Ap8EHp/fXwz8P8AVwMls/5Vm9p4DtLsj7pQcx3GWCi0Y6LBwMARmdgNwwy5l1+6wb8/KzezfAP/mbNrdDXdKjuM4S0TKPHuo03fnFe6UHMdxlghD1BQL2R1F3Ck5juMsEX6n5DiO4ywNhogLrVPyOyXHcRznIWCRkPCjijslx3GcJcLQgtN3fqfkOI7jnGMMiIe7ePa8wp2S4zjOEmEmxrZA9N0CC2zPR9wpOY7jLBEpdYXfKTmO4zhLgRYMdPA7JcdxHOcck9YpHVrqivOOhz3uUNJXSnqTpL+TZJKeM1d+Y94/u906Z9OX9IuS7pO0ket77JzNVZJ+L5ffJ+kXcspex3GcpSGtUwr7bkc1+u5hd0rAGkk99vuArV1s/pCkMtttz5grfyXwz4H/A3gaSUr99yUVAPn1v5PS9j4t2/0L4OcO80Qcx3EOg9bCvttR5WGfvjOzNwNvhnRXtIvZyMx2zFYo6SLgO4HvMLM/yPu+jSSz/tXA24CvAZ4MPM7M/ke2+WHglyX9WzM7dXhn5DiO8+CJJuoFou8WSW9xPnK+uNunSrpH0kckvVbSo2bKngJUwNu7HdnxfBD4irzrGuCDnUPKvA3o5+O3IelmSTcDV8+XOY7jnGtawr7bUeVhv1NagLcCvwPcTkow9dPAH0l6ipmNSEmmWlJO+FnuzmXk17vnyrsc8lfgOI6zJBhaLMnfEX2mtPROycx+c+bt30i6jTQ19/UkZ3Uu2rwW0h0T8PRz0YbjOM5OnIN06OcV5909oJndAXwaeGLedRdQAJfNmV6eyzqby+fKL8vH7fisynEc5+FBRAv7bkd1ndJ555QkXQZcCdyZd90G1MB1MzaPBZ4EvCvvugV40lyY+HXAKB/vOI6zFBhQW7HvdlTvlB726TtJa8Dn5rcBuErS1cBn8/YS4I0kJ/R44D8C9wC/C2BmJyX9CvCzku4BPgO8AvhrUig5pCCI9wOvk/QDwKXA/wW81iPvHMdZJgzRXsD5lJbhTulLgL/M2wrw0vz3T5ECEb4Q+G/AR4BfAz4MXGNmp2fq+H6Sk7oJeCewDnyjmbUA+fXrgc1cfhPJ0f3gOT43x3GcA2HGQtN3dkRvlR72OyUzu5m9J0e/doE6RsD1edvN5lPANxy0f7MoCBRQkdcQKH8qJMBAQgr5UxXTaxuxtoWmgXycQkAbI2gawriC0+uEskCjiqJXUmxVKILaSHHvKVjfQGUJEkVVQV2ndusa2xpi4xqCCFtD+ifX6K2tYFWBRjXaGsPWkHj/A2hri3IwwIZDyn4/9Xs0JvR7IGFbW0gB9XvQttjmFtY00LawuYlOnoYiQN0QiiKNQxGg10vnagYhpHNvmnRRzbBxTSgC/ZWVZD8aUxYF/RNrtBetEHsFzUrB1iMrmoEYnyjpnRiguqU53ies9Sn6PRiNoVehNqLNrTSm/R524hhxtYdJKEbURCwE4qDAgrAiEHuBWAkZlOsNxVadrlcbqU43lMOABbAgMCOM07W1SjSDQCwh5GFvjpU0x0pUR8qtlrDVoO4XwqAYRarTRu+BMeU9J7EHThJPradr36tQCGlcmwYVBYot5ek+a3dWhDZwxwOP5ZfXriRWBgIrDVZbVo4PufjYFo9dO8kVg1McK0eshjER8ZnxGncNT/Dx+x/B6L4VVu4LDO6DY3e19B+oCeP0eYxVIPYLYi/Q9EU7KGn60PZFfQzqEzBaM7YuqRleXFJd1HKiN+TxK5/hyv79XFxsMtCY2koeaFe5s76Yu0YnuH39Uu48dYKNjQHD0xVhvaBaF2EM5Sb0T0IxMkIDoQEQFgraqiD2Vriv9wjuDekrpQbKofGIIZSbaYwVDRO0g8D4eJ96dcBw7Tjrx6DtAQEGLRTDdBnanmgGqf1iaBSnRwzGLTQxXasYoTXUtOnvuoa6mX5viwJ6FawMaNcG1Jf0GV18EfevXczdF11FfQzaPrQDw4rUb9s6rJ/TxfIpHdVnSg+7U3Icx3GmpHxKrn3nOI7jLAFJkHV/maGDOCVJL5B0u6ShpNskPW0P20dLeoOkD0lqd1LakfRdkv5M0v2SHpD0x5KeOmfzkh10S/eNdnan5DiOs0QYWjD6brHpO0nPBF4F/Afgi0hRyW+RdNUuh/RJ4gIvA969i821pGfz/yvwZaRn/W+T9MQ5uw+zXbf0C/frr0/fOY7jLBO2oK7d4rdKLwJuNLPX5vfXS/o64PnAj51RrdkngBcCSPoXOzZt9q9m30t6PvDPgK8DPjpT1OymW7obfqfkOI6zRBwwdcXVnVZnVqDZRk7P8xRmtEEzb2eqDXoY9IABcP/c/s+RdEeeOvxNSZ+zX0XulBzHcZaILsnfftuCN0qdcs289uesNuhh8NOkpThvmtn3buA5pLun78rtvUvSpXtV5NN3juM4S8YB0lK8t9PqfLiQ9H3A9wBfPStGYGZvmbO7Ffg48GySwMGOuFNyHMdZIixr3y1itwBdNoR57c9ZbdAHjaTvB/498E/M7D172ZrZuqT3M9Ut3RGfvnMcx1kizKC2sO+2iKKDmY1J+p7XzRVdx1Qb9EEh6UUkh/T1ZvaOBewHwOcx1S3dEb9TchzHWSoWu1M6gKLDK4DXS3oPSWbtecBjgFcDSHodgJl9+6TmpD8KcAKI+f3YzD6Qy38I+BngW4GPSOqeT22Z2cls83Lg94BPAY8C/h1wjCQXtyvulBzHcZaIw1Z0MLObcnDBi0lrhd4HPMPMPplNdlqv9Jdz77+RlMfu8fn995Iyft80Z/drpOAGgMcCv0EKtrgXuBX48pl2d8SdkuM4zhLRRd8tYrdwnWY3ADfsUnbtDvv27ICZPX6BNp+1YPe24U7JcRxnqTj06bvzCndKjuM4S4QZNItE3x1RRVZ3So7jOEtECgm/cJP8uVNyHMdZMhabvjuauFNyHMdZIozFFB2O6OydOyXHcZzlQguFhHugg+M4jnPOsQVTV3igg+M4jnPOMaCJi2jfHU3cKTmO4ywVi0Xf+fTdBYpZBIsAqCigKEDCxjUKwqKhNtlY06R76vEYa9vuPjxXFKEsUQjYaAx1DW1MdW1uEoqANkskEYqQjjPDNjax0QgkMENlibUtCuk/qdROJA4brG7Qxibh2Crq9dIxAFI6bjiCNmLjMba+AQpgEY16UBRYXUMIKCgdGwIqS+j3Uz1BEC21Oa4n/6mpKidtqdeDIKxuoGmwtk1/A2FcT8YOQJtbFKdXKYpA1etRnTpOfVGPYtRSnBrCuKZqDSsD6sayaWcvznScpO3f0QAWROwVxEo0qwX1qmj7IjQFoe5jAdq+iFU+jxZCYxQj6K1HwigSxkahSKhFURvFZpuudxCxCrT9AhPIoB0UNKuBti+0AooVxXofbfVR2EiflapCx1ZhZYD1KqwqqC8esP73Bpz8n8TwipbBozY53q8ZjirqcUlsRNFrWenX9IuGcSy4Y3gRTQyMY8FWXXF63Gd9q8/oMyv07y4Z3Acrn4n0TjeEUQsGVgorArFQ6nOEYmSohXIE5RCqDdEMRHN/j/Hxio+trfGxlZZbVx7PYFCz2htTFpEgo24LhnXJqC4ZbfRgs0RjUYxEMU4XI/agLqBZhdAIYmoXgRUQq7yVaZ9aKEbQboi2BxYCFtJ1QalvbT8dYwVgEOruXKDcgmrDKIdpS30IxF4PC3kM8ucktBDGMW2jljCqUZ2/tyEQ+xVxpWR8UcXo4oLRRaJehfo4tKtGrAwrDEVBrUO7dTlsmaHzDXdKjuM4y8Thp0M/r3Cn5DiOs0R4SLjjOI6zRPgzJcdxHGdJMIN2kei7I3qr5E7JcRxnyVhs8ezRxJ2S4zjOEuGCrI7jOM5SsU+OvSONOyXHcZwlwqPvHMdxnOXBtFCgA0f0bsqdkuM4zhJhLBZZ53dKjuM4zkOCR985juM4S4MHOjiO4zhLwYWeT+nCTQTvOI6zpJjtvx0ESS+QdLukoaTbJD1tD9tHS3qDpA9JaiXduIvdP5f0AUmj/PpNc+WS9BJJd0jaknSzpCfv11e/UzoobQsKKX0DBZDTRwCSYbQTO4uWUlrkV7UtFMVkttjaNqVyaFtsNE77u5QThbKEvqZ1QU7VELGcdmLSdhAqipReIxp0aTTalFLDRiPieJzK0wGEQYnFmFJL5PQS1kYsbk4/9RIyw2JMdeVUGV1/FIQ19cz7zVRP7ucEBeLW1tSuKFDTYFtbSClVR7HZR2Zo2MDpDahrlG0ZjWBrOKnO2jYd17ZofUi5NU7n3OTUA2VB0e9jvQLrlcSVHs2xkvpYgRXT/0KL2rAgQm2E2ii3IsVWQ7lRo+EYjWdSZbQtjMYQU8oRcpoO2hbaSNWr6K+t0K72iIMCtYZVBVpdQaNjsDVM6UE2t9J5mWH0CHVMqRY2RLEeGK72aJpAMyxhs6TYCNCKB3qrfLYfsX5EwbBWaBwIo0AxFGEMa6eh/wCs3tfSe6Cmun8LbY4mY5LGosopPYqc1kHEQoRahDqngBiJclM0pwvafkGzWrEx6LPej6g0VMT0eTdhUVgdUEwpPBSV0oDUoAaKMYRm+j60llJphJQ2pO2T0lSU03QUVqS0F81AYIFinI8pUl9jD5pjUB832tWIlQYy1ATCSIRhoBhCMSwITUqJgTGRi5Plvo2gHBq905HeyZJysyasD9HWmGJzRLFeErYGhHEfxSpFxq1AXRm21lKtjVkZpHQe9xyrD/BDsjuGiAsl+Vtsik/SM4FXAS8A3pFf3yLp883sUzsc0gfuA14GfPcudV4D3AT8JPA7wDcDvyXpH5nZu7PZDwM/ADwH+DDwE8AfSPoHZnZ6t/76nZLjOM4ykafv9tsOEH73IuBGM3utmX3QzK4H7gSev2PzZp8wsxea2Y3AZ3ep8/uBPzazn8l1/gxwc96PJOW/X2ZmbzSz9wHPBo4D37JXZ90pOY7jLBkHmL67Ok+L3Szp5vl6JPWApwBvnyt6O/AVZ9HFa3ao820zdT4BuGLWxsy2gD/dr113So7jOEtEWqek/bfFqruM9Jzh7rn9d5OcxoPlin3qvGJm34Ha9WdKjuM4S4UWDAkXwHvN7Npz25+HFndKjuM4y8ThpkO/D2iBy+f2Xw7cdcCezXLXPnXeNbPvU7vY7MiBpu8kfXkO8XurpL+W9FFJt0i6UdJ3SLrkIPU5juM4O2ALbItUYzYGbgOumyu6DnjXWfTwln3qvJ3kfCY2kgbA0/Zrd6E7JUnPBn4QeDJwGvgr4KPAFvAI4MuAbwN+SdJ/AV5qZrcvUrfjOI4zpXumtIjdgrwCeL2k9wDvBJ4HPAZ4NYCk1wGY2bd3B0i6Ov95Aoj5/djMPpD3vwr4U0k/CvxX4JuArwKemusySa8EflzSh4CPAC8G1oE37NXZfZ2SpL8GHgm8Dvh20hzmGeMh6SLgG4B/BXxA0nPM7Kb96nccx3G2c5hqDWZ2k6RLSU7h0cD7gGeY2SezyVU7HPaXc++/Efgk8Phc57skPQv4aeCngI8Bz5xZowTws8AK8EvAJcC7ga/Za40SLHan9CvAa8xsuJeRmZ0E/jPwnyX9Q84ussNxHOfCxBbUvjuA4zKzG4Abdim7dod9+3bAzH4b+O09yg14Sd4WZl+nZGavOkiF+Zi/Ik3xOY7jOAdCC+ZKOpqirQcNdPgRSb+4S9kvSPqhw+mW4zjOhYvF/bejykEXz34H8Ne7lL03lzuO4zgPErMFF88eUZXwg65TuooUdbcTHwced3bdcRzHcY5sWtkFOKhT2gSu3KXsscDooB2Q9JWkcPOnkMIUvyMLAXblIinRfjfTCI7vNbP3z9hcAvwC8L/lXW8CrjezB2ZsvhD4v4H/hSQy+Brg3+8USbi9fwEUUIgzqt9NVuhm2310p9id3uT9CttfARuPt9la22LDUVKfbpppu0FJ8burO6uCKygpIhcF6hVJyTuEpPjdttjWVvp3S5qogqvfJ+R6un1WNxCySvmk3qS8jcV03kUBVUVYGSSbuoGQzyXG9HfbYuPxzKkboRJ0iuVBqFM/71TNB304dgyqEoKIqz2aEwMwo9rMit/K59+OYWMLy2NDPk96PRSL1NeWpBDeqaNLUDcoRtRE1Boat5TrASwpdAPEKmBBKFqyqSNh3KBxA21M45PPVwD9XupDNwZ1DeM6j1m6JoUZYRTQuEHrW0ndfFxProeKIvWxaVAIhHFLOYz0TwaQqLf6WNFnMIZiCCEPbeyJ2Euq3bHMHZJBVuXuVLohP5IQWFZ5p65hOEISxemSoghJNbwMUKTX2CuJg4JmUFCvFYxPaKL4HUsRe0KFUfRa+oMxvaqlCOlz3sZA3RS0baBpAuO6wJqANUKjpGRebiUV8jAWalJfLSRFcCtSn2VTRfFOwbvcioTGkEEbRVmkcTIBUYRRgZWGBdK1bUVo0/m3g8lXaFI/uY1iCKFNfUBASMrlVpbdlx+rCuJKmdTU85iHBoqtQAuM2z71RgVAMz48LYILOcnfQafv/gz4IUn92Z35/Q/k8oOyRgpR/D7Suqd5Ovnz64EvBe4hyZ8fn7F5A/DFwNfl7YuB18/07wTwByTdpS/Nbf0QST3XcRxnuTikxbPnIwd17S8hrcb9iKRfB/6OdOf0rcClpLwZB8LM3gy8GWA+mdS8/Hne92ySY/oW4DWSnkRyRE81s1uyzfcAf5bzdnyYtHZqFXh2Vqp9n6TPA14k6RX73S05juM8tPid0kLkUO+vIi2i+hHSdNiPkCQlrs3lh8ki8ufXkFYJz0pXvBPYmLP5s3xsx9tI04WPn290Rgb+6vkyx3Gcc4oBcYHtiP4rfeBJUDN7D/CVklZIz3jun/uxP0z2kj+/csbm3tm7nSxxcQ/bZdQ/vUMdXZlLIjmOsyRc2OuUHtSTufyM5gtIjuHTkt63n3TE+US3wjnfLT39Ye2M4zgXHBfyA4UDJ/mT9BPA/yAFNdxEmir7tKQXH3LfYLv8+SzzEumPzM+fuj4KeNSczU51zLbhOI7z8LNIkMMRDnY4qKLDS0nBDjeRJMm/EPhq4L8AL5X0kkPu3yLy57eQIviumTnuGuDYnM3T8rEd1wF3AJ845D47juOcHab9tyPKQafvvgv4OTOblRN6P/BHkk6S1hK95CAVSloDPje/DcBVWSb9s2b2qf3kz83sg5LeSorE++5cz2uA38+Rd2TbnwRulPTTwN8HfpSUYuOI/r/hOM75io6wjNB+HHT67iJS1NpOvDWXH5QvIcmk/yVJ5vyl+e+fyuU/C/w8Sf78L0jS6/Py599CEoB9W97+ipTfCZgomF9Hirb7i1zXz5HyjDiO4ywPxmJ3Skf03+mD3im9m7T49A93KPvSXH4gzOxm9ggjWUT+3MzuJ62V2qudvwG+8qD9cxzHecg5og5nEQ7qlF4I/K6kBvgtUlj15cC/BP418E+lqZ6O2VHWsnUcxzlHuFNamE4h/GV5m0XA38y8twdRv+M4juNOaWF+igt6uBzHcc4x3TOlReyOIAdySmb2knPUj/OTRWYns9r2NhRIstbztjZR1TaLKIYdV9FZNBRi+uC2qT5rc1/CVBkcae70+DgAACAASURBVKJmbXU9adsmssn1GfVCUhFXWU3Pr22xMUkV3Cy9b9vUV0vq6TavLK4wlWeOllTEu/dFgHGN2Ej9KwpC01I1qT1tDFN/25gUrpWVzLs6gtK5mWGjcVLAbmMqb9vUdhFQVaXXskSDPkW/T1ztY4MS6+rIm4XtEU9WCBGS0ngzmiqQx5jOxyLUDdY0SW09n3sY9NHKAK0MoNeDQT+pma+uorpOfc1jSHedANWRcsuIlYgltCtJObtdSSrWoQa1STG8+72yArCkvF1upa13yug/EOl/dkx5ckR44DR2eh0bjrBxut6qyjQmvV5Sa18ZwKBCVUEsA7En2oFo+6LtZaXtvmH9SNFrKXsNVdlShZZudWAILUUvYiaaNlBXBXVd0IxKaIWGuf+j3Nf8Xq0ldfOYFLvD2CjGkXIzUoxais0ajWqUPxtWBqxXEvslzbGKZjXQrATa3rS/sUpjE0uwcqpE3imEK5JUynO7KKmgt4MCkyjKgOoqX5g03sUoMrgfqnXRfyCkdsqQ+pPH4L51HTxNwi549N0hkcO7HcdxHOdBcdDFs7+wR9kau4eLO47jOAvQpcnad3u4O3qOOHA6dEk/Nr9T0ippndLfO5ReOY7jXKgstEbp6Ko6HDTQ4X8H/puku8zsP8E2h/QEfB2Q4zjO2XNEgxgW4aD5lN5Kkhp6taRvyOkr3kKSCbrWzD52DvroOI5zYXHIYqySXiDpdklDSbdJeto+9k/PdkNJH5f0vLnyT0iyHbb/PmPzkh3K9xXAfjD5lF4n6QqSCOvfAI8jOaSPHrQux3EcZw5bMPpuQcck6ZnAq4AXAO/Ir2+R9Plm9qkd7J9Aygb+qySlnKcCN0i6t8sATlLwKWYOezRwG8kvzPJh4NqZ9zuEHW9nX6c0q9Aww8uBxwLPAv4xKT16AFdxcBzHOWsOd/ruRcCNZvba/P56SV8HPB84I0YAeB5wh5ldn99/UNKXAT8IvBHAzO6dPUDSdwKnONMpNWZ2oPRAi9wpNew+RALeO/PeVRwcx3HOEi3ulK7OyUiBaYLSST1SD3gK6UZilrcDX7FLndfk8lneBjxbUmW2fYFjzl/3ncCv75CF/HMk3QGMSNqoP25mH9/rhBZxIK7i4DiO81ByeJF1l5Gm2e6e2383KRfeTlzBmaLbd5P8xWXAnXNl15EC3V47t//dwHOAD5GSrr4YeJekJ5vZZ3br8L5OyVUcHMdxHkIWDWRINu+dvzt6GPgu4M/N7K9md5rZW2bfS7oV+DjwbPZIG+RTbY7jOEvGIcoM3UcKLrh8bv/lpKzeO3HXLvZNrm+CpEcB/xT43v06Ymbrkt4PPHEvu31DwiW9SdIX7Wc3Yz+Q9KL5EELHcRxnfw5T0cHMxqSouOvmiq4D3rXLYbfsYv8X88+TSNNzI+A39uuLpAHweZw5/beNRdYpfQK4VdK7Jb1Q0hdL2naHJekxkv6ZpF/JDX4n8P8tULfjOI4zyyJrlA62VukVwHMkPVfSkyS9ipSF+9UAkl4n6XUz9q8GrpT0ymz/XJLz2RYskQMcngv8ppmtzzcq6eV5vdMTcvTebwPHgF/bq7OLPFN6YT6J7ydlf70IMEmnSB7yYqBHctzvyXa/bmb7xqOfF8yqdB802n0mml5h5v+arKy9vZ2ZunOZhZhUqeftOqXvHKKjgqScDVmhOyloJ9XxTtl7h74rTNpVUFYf1w59jQiwri5ISuEzKuG7jk2nkj5rY0rH11nRPEYIAWVFc8Y11E1SAW8apJBUtXObMiWF7XwcrU3L2xZokRUYoDaNVXdGoQhEgKpI/5JZgKz0bCEPSRnS+YakUC4pFXTjGiPqzr9p8mkaKvIYtS1qWihnxiQoqaNbOR3DTqmcdClDa4RaFMP03pSVteukoA2dMnguL3KfWyBy5o9UAIoijV8ISX0+ZpX37iPQvcZIkCiDkJVYIWIRiKVyu8KKQCwL6k78vQwURSTkz2HMD+cloyySYjg0NFG0jVAbJi1aSOekVlO17ghFBVaE9NFWUhEPZlhhyAwrA7FfEnsFsVLqV1et5XHKlwlLw0KR4wYCE9tYpTbbCoiiqKCtlBTN2zD5b92CiL1AOwg0g0AzEM1AtP10KWOR+xk5XHnrQwwtM7ObJF1KCjR4NPA+4Blm9slsctWc/e2SngH8PCls/A7ghTNrlDquJU3F7Zb1+7GkO6jLgHuBW4Evn2l3RxZ6ppSVGq6X9AOkcMEvI3naAfAZUnTFn+7XmOM4jrM/BwgJXwgzuwG4YZeya3fY9yfAF+9T5x+zxyyimT3rYL1MHDSf0hj4k7w5juM4zqHi0XeO4zjLxCHLDJ1vuFNyHMdZNo6ow1kEd0qO4zjLhjslx3EcZ1k47ECH84mDpkO/QZIn8nMcxzlXHP46pfOKg0bWPw94q6T51b5IWpP0uYfTLcdxnAuTw1R0OB95MMu93gb8V0lfNbf/yaSETo7jOM7ZEBfYjigP5pnSfyRlnH2TpH9iZu845D45juNcuNiCz5SO6PTdgwp0MLOfyMmj3izpa8zs1kPul+M4zoXLEXU4i/Cgo+/M7EezY3qrpN2SRZ33GCQ9tU7DrXvdhaluXHGm9p2S/hhtO9Gay4VJyy3roE200eLM3/u2lw+NEYWQNNgg6clJqCqxtk06eR1ti7VMNO+wyESxcPZftXGdviOddt2MZt60D8X0fCwSx1PpQxWp0W4cVRRQVhAjZgbRUDVEoxGEgNU1NhoR64YwLrGimOgBTvrYtth4nOoqClQUqNfbrpHXNKnfdZM06sY1Go4InW6eWRr3qoReD4qAhQCFUpuFoCqwQfqaWFFgZcDKrNfWtITNMWFzlPT6zGDQx1Z6NCdWaI6VtP0wnfy3pOVGELGAZiXQrIjRxWJ0CdQnjNhrIYBqUWwEyi0ITda/q5M+XqfxptYoh0a50VKdrik2xujUJmxsEDc2oWlou2tbFGhlhdCr0rmWReprv6Jd7RH7Bc1qQX2soM36bs0gacQBSY+vCdjJ3uRc6gbaOvVPTdady5epsLQppn4XtVGMjVAboY7TxaHKGn/ZVnWkGLaUGzUa1WhrDKNxvlYB9fuo6aFYTcZVMWDBUj0RFC2NT2OEJqLW0tZE1Mb0Wrdo3Ew/G0363ABJJ7AqYWWQxudYD6zKeffSdzq0SRevazM0Rhjv9A19kLhTenCY2YskVaTUuf/2cLrkOI5zYXMhh4Qf1Ck9Ffjo7A4zuz47pl84tF45juNcqLjM0OKY2Y5JoczseZLGwLcfSq8cx3EuZI6ow1mEQ8sAYmYvNLOLD6s+x3GcC5YLdOEsuMyQ4zjOUiEWWxh7VBfPulNyHMdZJha9Ezqid0vulBzHcZYMj75zHMdxlocjLCO0H+6UHMdxlgmXGXIcx3GWiiPqcBbBnZLjOM6S4c+UHMdxnOXhAnZKh7Z41nEcxzl7Fkrwt+hzp65O6QWSbpc0lHSbpKftY//0bDeU9HFJz5srf4kkm9vumrNRtrtD0pakmyU9eb+++p3SQdlDIXyiYD2xafMhyurb7aTM2qmKNgpZsbvdXn8RJirhSVV8h7YVUFVOFca39cem7URL/Yuz9eW/i2K7AvhsH/ZRRZ8cY/EMxfJ5ZfXpOOQyU+p3V320pJCtgHXjFA3lb9/kXHJdXbnVo7nu5LK2s41ZIT0rsfcqKPOYSViZFcHzmNEaogULWAFqIjSREMfJPttqWMPmFrY1TArTRYFiRHVNuT5MXy6LU6V3KbUbQlIk75XEXklzUY/RxRXDSwLNSkEskvJ277RRbRjVRkuoDQtZmboSls9R0Yi9wPgRfbioh644RjG+lDBqCaOGMKpTP9s4vcZZWd1CwCTURIKgika51WJFqj8Wwsp0vrEUbY9cli6dBZLKdxCxTIri3atlNfqJSvhYFCNRjI1iFJJieGOoNorGCOOWMI65vw0ajZN6d92kse3GcDhCp0QRAoVErxtXKSnwKymiT8a4DEl5PoT0Oct/22qBHR8Qq5C2fqAZBNp+UkdvBqLtp/OxKp/rrOr5KG3VplGtR6r1lmI4850+Ww4x+k7SM4FXAS8A3pFf3yLp883sUzvYPwF4M/CrwLeSNE9vkHSvmb1xxvTDwLUz7+cH4IeBHwCek21/AvgDSf/AzE7v1l+/U3Icx1kyDvMuCXgRcKOZvdbMPmhm1wN3As/fxf55wB1mdn22fy3wa8APztk1ZnbXzHbvpP+SgO8HXmZmbzSz9wHPBo4D37JXZ90pOY7jLBOL6N5NVR+uztNiN0u6eb6qnPPuKaT0QrO8HfiKXXpwzQ72bwO+JGeE6PicPDV3u6TflPQ5M2VPAK6YrcfMtoA/3aNdwJ2S4zjO0iGzfbcFuQwogLvn9t9Ncho7ccUu9mWuD+DdpGm5rwO+Kx/zLkmXztTRHbdou4A/U3Icx1k+Fp+ee6+ZXXvuOrIzZvaW2feSbgU+Tpqie8XZ1O13So7jOMtElxp+n21Bx3UfKQDh8rn9lwN3nWkOef9O9k2u78wum60D7weeOFNHd9yi7QLulBzHcZYKsWBI+AJ1mdkYuA24bq7oOmDHpK3ALbvY/4WZ1Tv2WRoAn0cKoAC4neR8rpuzedoe7QLngVPaLx5+kVh4SZdIer2kk3l7vSRPSOg4znJyuEn+XgE8R9JzJT1J0quAxwCvBpD0Okmvm7F/NXClpFdm++eSnh+9vDOQ9PK8lukJkr4M+G3gGClKDzMz4JXAj0j6ZklfANwIrANv2Kuz58szpb3i4ReJhX8DcBXpoRzALwOvB77x3HXZcRznwXGYMkNmdlMOQHgx8GjgfcAzzOyT2eSqOfvbJT0D+HlS2PgdwAvn1ig9FvgNUuDDvcCtwJfP1Anws8AK8EvAJaTgiK/Za40SnD9OqTGzM+Yh52Ph875nA/eQYuFfI+lJJGf0VDO7Jdt8D/Bn2XF9+KE6CcdxnH05B0n+zOwG4IZdyq7dYd+fAF+8R33PWqBNA16St4VZ+um7zG7x8IvEwl9DumWcncd8J7DBLvHyMzH/Vx/qWTiO4yzAIS+ePa84H+6Uunj4DwGPIt2Cvis/N9orFv7K/PcVwL3ZawPJg0u6h33i5R3HcR56DMULN6HS0julfeLhbz1HbV6b27oZePq0YEYHrtN0m3RsTi9utr521maq3zZbX6d7pzCjTwdZzyuAbFo2oyk3qbvTncsaeApzN8EF2HicNOqKYtqXGY25Mwci7vz3Dv3ffpjNvsmafdPHgNv0+Dq6fVZtPx+LqfNBqCiSFl5IGnIqCtTVN6vtpjDVQSNr+xVF0giEpKfWdNdBU82+aEljLca0MTeOMenYqdMZNMPGdRrXNgJ1rnIFjg+Igx6USb9QdZs06MbjSduqKop+lfXyItV6SewlnTnLTYQ2jYUFI9RpXscaEau0Nf1AzEOmCMXYsCoQykDoFYR+iVb6qY0mQtui1qCNqG5Q3Ux15cy2X5OgpCE3ow9oRQGB9Jo18lDqM0ET3TwL03OQJY0+NenHVm1MfTAmi0BNAkFc6cFKD2wlHxPRuIEmoiZr4GVNO+tXtMd61MdK2pVAvRqyXl2qiwhFbYQayqFRjCNhFFM4NUnHL/ZE2wvUq6JZSXp3bS/p96UPF4QRFGMoN6HaMnqnW6qTNcVmQ3F6Cza2YDxGm9s1GB8052D67nzifJm+mzAXD79ILPxdwCOlqWJp/vtR7BMv7ziO81AjFluntEhI+PnIeeeU5uLhF4mFvwVYIz1b6riGFL64Z7y84zjOw8LhhoSfVyz99J2klwO/B3yKdHfz78jx8PnZ0CuBH5f0IeAjpGdOk1h4M/ugpLeSIvG+O1f7GuD3PfLOcZylY9FAhiPqmJbeKbF/PPwisfDfAvwiSekW4E3A/3nuu+44jnNwFgt0OJosvVPaLx5+kVh4M7uflKzKcRxn+blwfdLyOyXHcZwLCp++cxzHcZaKxfMlHTncKTmO4ywRnUr4InZHEXdK+2J7LyLdb/+8jQIK2rbAVDMLQycLS0mLbmWa7p/rx2RRarcgNlWWFhfOLv6Upv95KUwWl4piUqfNStyeJfPnN7vYt+tDWgQ7s9i1o91hHLsFtApgdRoX0uyFQpgsjk22M+1quvhTVQlVBWWR9rdxuki2ZfK3tW1aQBqE8gJja1toGizG6XiGAG2L1c3k2qgo8gJlQ3VNmFlkS3cNJ4uELfW5CITTA3rHVolrfaxXEnvpXFRHQt2iUY1GzXThaFVg/YrYK9JWhcmC0bR4FdqVQDsIqC1Qa4TWUn1NWrw6O04WQv4lTMdCl6+nW1Sb980sslVrMGwJeUFxYZbGRUoLW6sC6xa55roxSz+2EpYXCCcbEQtNfmUVITRpQXEB0ARUGBYKCIHYr4grJfWJitGJwPh4oB1AMwArwZTrGEM5FMUoL+BthQqlzK3dYt5WhFGk3Jz2KZ0rqDZCYxTDhmKzJmyOYTiG0QjbGkLTELvvWa83Wax9KFy4N0rulBzHcZYKIzn9BeyOIu6UHMdxlo0j6nAWwZ2S4zjOMuHRd47jOM7yYAtG3x1Nr+ROyXEcZ8k4yvmS9sOdkuM4zrLhTslxHMdZBrRg9N1RvZtyp/RQ0625ySgn8Ztd1zNJ8he6tUjd2p/tmUa6Y1WkJHjd+pppRTZZr5T+7taZRJipa6c1ULP9gHaaCG+mfPL3Don+9mNbfZOq8nqrmfU8FtNTXymt89nWJ3KyOW0/50miutnEh3ltCm1e0zSXzA9Ix0VLbVvACnIiwSKtiWqaZFcUqZ+9CsoGza5h6lWo34eVwXTtVIwwGqPxGNvK655y+6KCpoGtYVrzUxaEKq3HMU37b2WBYoSc+I7QokJAQSxF7AcsMEmuN/nB6tYGdWm0o217kJ7WEDFZnwNM7dsZ21yPokEEtWnNE9G2J+oLMwn/uiSL0aYJ+6JBm46ZnN/Mx0dNTAkRx01Kxjgap/GJMV3HqqIY9AlbFWHUp1yv6K8U0+SIQZN+htoI45he65ZQ50SHs9e8W1slTfsRSWPdWnrtEgzGmD5T/X66xkEoBKwqsV4JH+t1eR7PGrmig+M4jrMUXOCZZ90pOY7jLBt+p+Q4juMsC0f1edEinHfp0B3HcY40WVtwv+0gd1OSXiDpdklDSbdJeto+9k/PdkNJH5f0vLnyH5P055JOSbpX0u9J+oI5mxsl2dx263599Tulh5q9ggIsnhHw0ImW2m6ZKC2mp9SdoGoMMw9s8+R0FwTR1b/tSXjuk7Lw7Iz46UQwdaY/84ENZwQs5DIVTJ+eZzHWbaKvWZh29tyUBUqR0kPuXP9sG5PgjqpM9jmIYSKkOhuEEcM02CFGzCKKmgqvzou/xu1CpUTDYpOEV7vgBCVBT7M4HXNmhG8hBUoMR1Mh3LbF6hrqBmuaqSirAkad6pwMcpUCOEpLGqmtTR/4d8eFgAZ9Ql0lsdW6JFYBC+lhvUmT4AWbHbsuYCGyPTih043thFPhzB88Ces+S0p9iGEmqGTGLrWdxokcdKAssKq2TecU45kP87sAg6aFup68zl4rZSFd1TVFE9GopRjm8y8DVijV0/2w5+CGMG6gbnMAw1ygQ5GPK4o0HoWSSG0gicB2usEzgrTbxiXb29xX4aw4xOk7Sc8EXgW8AHhHfn2LpM83s0/tYP8E4M3Ar5KSoz4VuEHSvWb2xmx2LXAD8OekT8RPAX+Y6/zsTHV/CHzbzPvxfv11p+Q4jrNsHO703YuAG83stfn99ZK+Dng+8GM72D8PuMPMrs/vPyjpy4AfBN4IYGZfO3uApG8DTgL/CPi9maKRmd11kM769J3jOM4SIct3evtuAFwt6eZuO6MuqQc8BXj7XNHbga/YpQvX7GD/NuBLJFW7HHOc5E/un9v/VEn3SPqIpNdKetRu593hd0oPNRa3rTfalldojm6qKv0dd53C27b+SDadVuqm2toz7RV2X1+047ThLjmlrGXuqWx75pTe7LRkN23XjcPslOJoBOM61ztdp2Sx2d6uLK3bmm1jJufUdJ3XzPRhaFJOpNnpqW3rnubOreunZvuYy5pmulZprxw63VRkr4f6PWhjWutkNpkKpCig34dBD+uVaRpJSmtkxs10Kqvt1smkfFlpGAzqlmCGxnnadr4/ZpPpLGKcTJ2dkXeq23bA5var63+3bSvUdI3V7Lofm7Y7eR4yu65uto0QoCwn+yRBMWc/sU05mWJV0A4CsReIXWqxvNYqtKT1SuNIGLfbpxEnecbSNBwhTeVZN6Un7bCGy6b5pjoHEuPhBicc3vTdZaSFkXfP7b8b+OpdjrmCNO02b1/m+u7c4ZhXAe8FbpnZ91bgd4DbgccDPw38kaSnmNlotw67U3Icx1k2FsgZmnmvmV177jqyP5JeQXru9FSz6ZNjM/vNGbO/kXQb8Eng60nOakfcKTmO4ywT+c5rEbsFuI80V3L53P7Lgd2e9dy1i32T65sg6eeBZwFfZWYf37u7doekTwNP3MvOndISsWeE3Rw7yf7sOBW4LQX7NH26xXDmFF2XsryLlJud3pv9e74/c+/POLabRjwj4q7dFs2Xpuxm5IV2amvWdqfpt2396KLcLKVQn4l625WubBIt2O5cHg3TTHr0jrrB6mYyzZYkoEKapuum62DbtNHk+hVKacRDyJI8EZqcun1cp+jBrp6ynEoVNSDitjon02fbUrKT07/bdApu27lNz2P2eM1M7VlX96xEUCfVM5numo34S7JCqV2hGDEiaoFC02k8y/1vI8R2mq6+i3qENG1ZlimlfVWlCDmy5FGd2rJxnMobWZfyPCaJoXGTJIyaHOHXtKmtmMZC3RgEQciRoPk6psi+IkXmlWGa6r0QsXutDvER/SFN35nZON+hXAf81kzRdeSghR24BfimuX3XAX9hZhMhJUmvAp5Jckgf2q8vki4DrmTn6b8JHujgOI6zTBj5H4h9tsX91iuA50h6rqQnZWfyGODVAJJeJ+l1M/avBq6U9Mps/1zgOcDLOwNJvwR8B/AtwP2SrsjbWi5fk/RySddIeryka0lRefcAv7tXZ/1OyXEcZ8k4TEFWM7tJ0qXAi4FHA+8DnmFmn8wmV83Z3y7pGcDPk8LG7wBeOLNGCdJaJ4D/d665lwIvIU0xfCHw7cDFpLujPwb+pZmd3qu/7pQcx3GWikXVGhZ3XGZ2A2mx605l1+6w70+AL96jvj2XCpvZFvC1e9nshjulh4sdnhMdTrW2S0j2AjO1h9mnnZ47zaXomNhMnnEtFgKfDpmWzStD7Nj+1Hix/u/1DE1hkrbgjDDq/Jxj2/Okbl+XKmP2eU5X3rRorPS8Q4K2RcM6KTmMx0kJIkasjVP1h6ZBdZmO3xZSrenzEZh5bmPbn9PslOKjO7fJe02fs5DCs7eFiHfDpLk2OzWLrv3ZUPB2ZgxmQ+07+8kzzvwMbvahf/dcME7X61hrqCmIvbRZJSxADAEqkBU5lUWZntU1RuieL9Uzz69mr2ER0vO9MhD7VQo17+c0Ib1AW4lYkd6XYEG0f3WIkg4uyOo4juMsBUYKDlnE7gjiTslxHGfJ8CR/juM4znKwU7j+bnZHEHdKjuM4y8Yez1OPOu6UloDdHtBv05XbaYFqR7foddtxM/tm0lAoKJXPBhccNjMLUM+ujp0Xrh6oz/OBDbv0aa+givlxtBZkmi7OnW8vCBUFpma6WBamgQZdqg1IAQJVhbaGqCy2L64d11jTYONxWpCb+5FSfBST19lAhF2ZDXTo0lbMPNxXF4wxWRQ7/15opyCIebv5v7eN40zgQ4xn7uv6CSk1yGy/Z8nnq1FaRKuyQFVFqMqkH9gFKFRpkWuXSsMkKANW8P+3d/ehktV1HMffn5m5u2qKooYPlFmpSFKuiqaS5j9GYWFQZNRWG2SUBGVEJRWEfyiVBSpFDwhuhj0gJRpYZiqKj6lsulJa9IRZ5pqLT5u7c+fbH+fM3DNnzjzeuTPnnvm84DL3nqffb2b3nu/9fc/vgViqpfMAprMntLpTZtFZviIZFNtaqiWdGpaUvpIOnE06O0QDYpqjPivaChqFg5KZWZlkg/aw4yrIQcnMrGycvrOZyi8DkUu/9T12lOv1KaOdqsunCjvjl4ZdY9TxPf2uNa4B1xgp1Vbwc8+cfEPKKbp2e/7ATvq0Z6mLQNQJlmmvUpsc1uoal5NdCiN27+msuts19ikiOa+d0uq8NmFPMzeuqztlVpQS7juOa9BKwv32t7cVLT3RVZFMmcNutPl6FtS7s/JwTUQ7hdlorHx+jUbyGdZEPTPuqktRPWqZMVmN+spcd40atJJxTlquUWuK1jLJaz2oLYmoQ6sOag5+eyMLRvt/WdG45aBkZlYq05/RYT1xUDIzKxun72yu1mjKoanKLDExkUHT9qTb2qmxrhRStidfJlXWt47DrPZ9FClMCdbTb3vrFNkUHtDuZdh535kUniKgttLLr726bvI9dFbhbfV+vgPTnIN6MhZ9Pp1rrkwDhWorq8G2V9KF7imM8qv2jlJ+T3VWUnaohpYaSe/D9jIW7VVqG43uVF0rkqUpmpFMT9RspkuetNJlTzK9D9s9GpcaSY++DRvQxiViqU4NaAlUF+kCF8lqtF1LdMDwdzKidHXbkY6rIAclM7NSGbH3XUWjkoOSmVmZBCN2CV/zmsyFg9Ki6DdTeDYlNiz9tZo04wSptZ7egsOM8h5WO7C3z3nZnnDAymzWXadGd2ot3zMwWkQmTZXvNaaaUK3Rk/6L5eWVQdZdx2ZWku3zOfZsjxY9g5a7LCezvQ86ZCTpCsOZ6/SkbdU9c3ynhxwkM6Y3m9BsgnZn6h+d/Z3ei+3t+TRi9rrp1D6xvJystAvQWkZ7GujlJmrUqS/Vk0G19aR3HrVkIC21ZNXd2u5pzrJf0YgzAgclM7NSce87WxSDHvKXobPFqHVbTWeFcd7ngAf+vZuj8/C/GO5r4QAACG1JREFUsFWSbRXl94/yXmqCVkFLMLMmVfZaY7UwxzWsw8kE8uOfulqLyqxd1Xltr++Uvt/O56K080UNWi1UqyWtJSlZ+ymbLci2vLLjrCK61m2ivoxajWTNq3qNqDc79ajVk2tErUatOaXPPBit9101Y5KDkplZ2UQZ/kicEwclM7Myaa/OO8pxFeSgtGjK/BfYoCXTi46dUX0GlpdJYXVSeD2X6b15JOmqes92opU+/B+xJ8E807HjplFHnLZK+bFGtXZnDvV0IBn1ttyZBb2RGwuWmf28K5XXPr6zfPuezj5l6i9Y6ZCyZ1rzDDFil/BqclAyMysVd3QwK481eIi+KuN0k+/XyaHwtEhPWYP1rGbZGp5S1/pO13BFT3fzns8o939j4P50X6el1Skms3ZS0WSytQn+/02rpRTFM4EUHVdFDkpWTmVOM/YzSp1zN9SeGbBHKmZw7751rSDF1xPAB6UNB/0xkw1Ew276RQs4DjPNZzwVfV40CgclM7MyyXZJH3ZcBTkoma21IanIcVJ4lW4l5Q0Ye5XtUJKdqaSro0l2wtjOtvwsGQP+bfqtETWt4wdY03FmJeegZLbWhgSNwil7Rn2mVpWAtJrAnQle3Tfz5a6XwlMHFTruZztJyq+w3IJZ1fsdV0EOSoMd9Tw7eSBun3c9bNFU837T37D3uw5i7/PsBDhqtdd5lqf5XfPWocft5OnVFlVKDkqD7b9Mc/dOdtwz74qscwcBz8y7EhMqS91nVY+1KGda11zNdSY9d5zzTgP2n6CMrG0wVsDZtsrySkdR0SbgNEi6HSAizppvTdY3Sd+PiI/Pux6TKEvdZ1WPtShnWtdczXUmPXec83y/mI6SDAaxirtx3hVYhbLUfVb1WItypnXN1Vxn0nPL8u+/MNxSGsB/+ZjZqHy/mA4HJTMzKw2n78zMrDQclMzMrDQclMzMrDQclMzMrDQclCYk6QBJD0jaJmm7pPPnXSczKzdJ+0j6u6TL5l2XsvKMDpN7HjgzIl6S9Apgu6SfR0QZRv+bWTl9Cbh33pUoM7eUJhQRyxHxUvrjRpKVkddgxTYzqwJJRwPHAjfNuy5ltrBBSdKZkm6Q9E9JIWlLwTEXSPqrpP9JelDSGbn9B0j6PfAE8I2I2DGj6pvZDE3jfgFcBlw0kwqvYwsblIB9ge3Ap4Fd+Z2SzgMuBy4BTgDuBm6SdET7mIjYGRHHA68FPiDpkFlU3MxmblX3C0nnAo9HxOMzq/E65RkdAEkvAJ+KiKsz2+4DHo6I8zPb/gRcFxE9f+1I+g5wa0RcN4Mqm9mcTHK/kHQpsJlkdad9gSXgmxFx8Uwrvw4sckupL0kbgJOAm3O7bgZOT485RNJ+6ff7A2cCj82ynmY2f6PcLyLiooh4dUQcCXwO+IEDUjEHpWIHA3Xgqdz2p4BD0+9fA9yZPlO6E7gyIh6ZXRXNrCRGuV/YiNwlfEIRcT+wad71MLP1JZv2s15uKRXbQZL7zXdcOAT49+yrY2Yl5vvFFDkoFYiI3cCDwNm5XWeT9KoxMwN8v5i2hU3fSdoXOCr9sQYcIWkT8N+I+AfwLeAaSfcDdwGfAA4HvjuP+prZ/Ph+MTsL2yVc0lnAbQW7tkbElvSYC4DPA4eRjFG4MCLumFUdzawcfL+YnYUNSmZmVj5+pmRmZqXhoGRmZqXhoGRmZqXhoGRmZqXhoGRmZqXhoGRmZqXhoGRmZqXhoGRmZqXhoGRmZqXhoGSVJ+ljkqLP14uS6kPOv0LSL8cs83pJz0ra2Gf/fmnZV6c/f0bSI5L8O2kLzb8AtghOAF4ETiv4Oj0ilvudKOn1JJNrfnXMMrcCBwDv7LP/vcA+6XEA3wNeCXxkzHLMKsVz31nlSboLaETEmyc490rg1Ig4eczzNgBPAndFxLkF+28DXgccGekvoaSvA+dExHHj1tOsKtxSskqTJOBNwMMTnLsR2AxcW7DveEk3pCm6XZLuknRGe3+6xs6PgXdIOih37hHAW4Frovuvwp8Ab5B0+rh1NasKByWruqOBfYFHJTVyXwOfJQGnkqTg7sxulHQiyeJtBwLnA+8BngFukXRS5tCtwBLw/tx1NwMCfpjbvg14Hnj7qG/OrGqcvrNKk/Q+4Kd9dm+LiBMGnPsF4FJgr7Tl097+W5IF3I5vb08D3HbgsYh4d+bYR4EXsqlDSX8AdkbEaQVl3gnsioi3jfE2zSpjYVeetYWxCQjgLcDu3L5nh5x7OPBcLiDtTZJ6uwRoScr+Dt0CfDB3ja3A1yQdExGPSzoFOBb4ZJ8ynwaOGVIvs8pyULKqOwH4S0TcPcG5ewEv57YdCNSBr6RfPSTVIqKV/vgjktbWh4Evp68v07/1tgvYe4K6mlWCg5JV3Sbg3gnPfYbkmVLWTqAFfJveZ0IAZAISEfGkpN8AmyVdDJwH3BgR/VppBwI7Jqyv2brnoGSVJekQ4FAm6HmX+iOwQdKrIuIJgIh4MX3uczzwUDYADbCVpAffpcDBrIxNKvJa4P4J62u27jkoWZW1OzE0JZ1asP/hiHhpwPl3pK+nAE9ktn823fdrSVcB/yIJNicC9Yj4Yu461wPPARcC/wF+VVSYpANIniddNqBOZpXmLuFWZZvS14uBe3JfdwMbBp0cEX8jabW8K7f9IeBkkvTeFcDNwOXAG1kJZNnjdwE/I+kGfm1ENPsUeQ5JZ4xfDH1nZhXlLuFmA0jaQhJwDhvSqppGWTcBOyLiQ2tZjlmZOSiZDZB2+X4EuCoi1iytJmkTcB9wXET8ea3KMSs7p+/MBkhTbR8F1rSVRNIhY4sDki06t5TMzKw03FIyM7PScFAyM7PScFAyM7PScFAyM7PScFAyM7PScFAyM7PScFAyM7PS+D8eO46rK4YY8gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolormesh(energies,s.domain.rcen, Pradial)\n",
    "plt.colorbar()\n",
    "plt.semilogx()\n",
    "plt.xlabel(\"$E$ (eV)\")\n",
    "plt.ylabel(\"$z$ (kpc)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting the coherence length / cell size distribution\n",
    "Model B uses a scaling of the coherence length with radius. This can be modified or turned off via the command `set_coherence_r0`. Setting it to `None` turns it off so the minimum and maximum scale lengths are constant "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$z$ (kpc)')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAETCAYAAAA23nEoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAfhklEQVR4nO3de7QddX338ffHYIIGCLUgCIoxlSgG9ABVQDScalik5OnSVquyBKXV9FG8xQuWi0YaG8wfisFbbaNPsWqrrVWrEgKNEqKIEXkIlxiJGMOj3KEa4Agcjd/nj5lJ5gx7n32bvWdfPq+19sreM7Nnf8+ck/nu310RgZmZ2eOqDsDMzPqDE4KZmQFOCGZmlnJCMDMzAPaqOoB2SLoeOBC4tepYzMwGyDOBeyPi6Fo7BzIhAAfOmTPn0LGxsUOrDsTMbFBs3ryZnTt31t0/qAnh1rGxsUM3bNhQdRxmZgNjfHycq666qm7NitsQzMwMcEIwM7OUE4KZmQFOCGZmlnJCMDMzwAnBzMxSTghmZgYM7jiE0q3ZuJ3V67cBsGzRfJYunFdxRGZmveUSQmr1+m1MTO5iYnLX7sRgZjZKnBBSE5O7aj43MxsVTghmZgY4IZiZWcoJwczMACcEMzNLOSGYmRnghGBmZqnKE4Kkp0j6nKR7JT0i6ceSTqo6LjOzUVPpSGVJ+wNXA98DlgD3AvOAe6qMy8xsFFU9dcV7gTsj4nW5bT+vd7CkDenTsW4GZWY2iqquMno5sEnSlyXdI2mzpLdKUsVxmZmNnKoTwjzgLGA7cApwMbAKeEutgyNiPCLGgc29CtDMbFRUXWX0OOBHEXFu+vp6SYeTJIRP9CqINRu39+qjzMz6VtUlhDuBHxe2bQUO62UQnt3UzKz6hHA18KzCtvnAbb0MwrObmplVnxA+Chwv6XxJz5T0l8DbgU9WHJeZ2cipNCFExLUkPY1eBdwMrATeD3yqyrjMzEZR1Y3KRMSlwKVVx1G0YPk6wMtpmtnoqLrKqG95OU0zGzVOCA24wdnMRoUTgpmZAU4IZmaWckIwMzPACcHMzFJOCGZmBvTBOIQqrNm4fXd30mWL5lccjZlZfxjJhLB6/bbd3Uk9zsDMLDGSVUb5sQUeZ2BmlhjJhGBmZo/lhGBmZoATgpmZpZwQzMwMcEIwM7OUE4KZmQFOCGZmlnJCMDMzwAnBzMxSTghmZgZUnBAkXSApCo+7qozJzGxU9cPkdrcA47nXnlzIzKwC/ZAQfhcRLhWYmVWsHxLCPEl3AI8Cm4DzImJ7rQMlbUifjvUoNjOzkVF1o/Im4ExgMbAUOBj4vqQ/rDIoM7NRVGkJISIuy7+W9ANgO/B64KIax4+nx20ATup+hGZmo6PqEsIUEfEQsAU4vOpYzMxGTV8lBEl7A88G7qw6FjOzUVP1OIQPSzpJ0jMkHQd8BZgNfK7KuMzMRlHVvYyeCvwbcABwL/AD4PiIuK1bHzj3nEu7dWozs4FWdaPya6r8fDMz26Ov2hDMzKw6I5cQdqxawo5VS6oOw8ys74xcQjAzs9qcEJqwZmPNmTTMzIaKE0ITVq/fVnUIZmZd54TQhIlJz8htZsPPCcHMzAAnBDMzSzkhmJkZ4IRgZmYpJwQzMwOcEMzMLOWEYGZmgBOCmZmlnBDMzAxwQjAzs5QTgpmZAU4IZmaWckIwMzOg4jWVB8mC5euYmNzF7JkzWLZoPksXzqs6JDOzUvVVCUHSuZJC0ieqjqUomwJ7YnKX10cws6HUNwlB0vHA3wA3Vh1LI14fwcyGUV8kBElzgC8Cfw38quJwzMxGUkttCOm3+MXA8cAhwBOA+4BbgKuAr0dEOzf0fwK+EhFXSvrANJ+/IX061sZnmJnZNJoqIUh6vaSbgO8D7wSeCPwU2ETyjf444DPA7ZIukfSMZgOQtBR4JvC+FmM3M7MSNSwhSLoROBD4F+B1wOaIiBrHzQH+F/Ba4MeSzoyILzc497OAC4EXRcRvG8USEePp+zYAJzU63szMmtdMldFngX+MiEemOygidpK0A3xR0vOAg5s49wnAAcAWSdm2GcBCSW8CZkfEo02cx8zMOtQwIUTExa2eNCJuAG5o4tCvAz8qbPtnkuqoC4HJVj/bzMzaU+nAtIj4NfDr/DZJE8D/RMTN1URlZjaaGjYqS9pX0gclnVTY/oTuhWVmZr3WTAnhjcB7gN0NxJJmADsl/Ri4Nn1cExE3dRpQ1nBsZma91UxC+HPgczWqcPYiqeN/OfAGYJekZ0fEz0qO0czMeqCZcQhHApfV2B7AmyLiQOBw4CaSbqlmZjaAmkkI+wD319i+u59oWir4AvAnJcVlZmY91kyV0f3A0/IbImKXpKcytYfQT4AjSozNzMx6qJkSwjXAq4obI+KOiPhNbtPDwH5lBWZmZr3VTEL4JPAySWc0OG4+sLPzkMzMrAoNE0JEfBu4GPhnSR+StE/xmHTbu4Cryw/RzMx6oamRyhHxTkm/Bc4G3iJpLXA9SYng6SQT2h0EnN6tQM3MrLuanroiIt4r6T+B84CXMbVd4TbgzyPi2pLjMzOzHmlpxbSI2BQRLyOZofR5wEuB5wLzImJtF+LrWwuWr2PNxu1Vh2FmVpqWEoKkMyQtAB6OiJsi4sqIuDkiQtKsLsXYM7Nnzpj2dd7E5C5Wr9/W7ZDMzHqm1TWVPwfcCDwk6VpJ/yTpzZJOAP5G0vfKD7F3li2avzsJzJ45g2WL5k97/MTkrl6EZWbWE61Of/0k4Oj0cQxwIvDX7Bm1/EB5ofXe0oXzWLpw3pRtK9durSgaM7PeaikhpOsXXJk+AJB0IPB2krWWTys1ugGwYPk6JiZ37S5RFBOKmdmgaLXK6DEi4t6IeD9JddLIJYSs2shtCmY26DpOCDlXkEyFPbLcpmBmg6zVXkbvl3SqpKfU2P1kYKKcsMzMrNdabVQ+m2Q67JB0D/B/SUYs7yJZWe1vyw3PzMx6pdWEMIdkErtj0sexwFnA/un+iySdDlwH/CgivlpWoGZm1l2t9jIK4Jb08W/ZdknzSJJDliiWkpQW6o/sGlJzz7nUPY7MbCCV0qgcEdsj4j8i4tyIOCVdVvMZZZx7ELnHkZkNooYJQdI3JB3d7Akl7S3pXcCpTRz7Fkk3SnogfVwjaUmzn9XP3OPIzAZNMyWEHcAPJG2S9HZJx0iaUtUk6RBJL5f0WeBO4A0kDc6N/JKkaukY4I+B7wBfl/TcVn4IMzPrXDML5LwdeA7wQ+AC4FrgEUn/I+lOSQ8DvwC+CiwAlgHPjYgfNnHu/4qIyyLi1ojYFhHnAw8CJ7T9E/URz4ZqZoOk2QVyfga8TdK7SW7WxwGHAHsD9wM/ATZGxG3tBiJpBvCXJN1av1/nmA3p07F2P6eXVq/f5oZlMxsYrfYymgSuSh+lkHQUcA1JcnmIZKGdm8o6f5XcjmBmg6TMqSvadQvJN/7jgH8APifpyFoHRsR4RIwDm8sOYrq1D8zMRkGrA9NKl5Y6bk1fXifp+SQzp76hl3E0WvvAzGzYdaWEIKmTOv7HAT1ffa2sun6XNMxsUHWryuhTzRwkaZWkF0uaK+koSR8CxoEvdimurpk9c0ZTq6yZmfWrblUZqfEhABwMfCH9dyfJ8px/GhGXdymurtmyYvHu515lzcwGUdsJQdIhJIPJHrOLPZPdTSsizmz3883MrFydlBCeCDytzr6ZHZzXzMwq0HZCiIhb2dM7aApJr207IjMzq0S3GpWbbUMwM7M+0a2E8MEundfMzLqkKwkhItZ247xmZtY9/TB1hZmZ9QEnhC5bsHydp8E2s4HQzIppcyVdla5o9kNJC9PtsyX9haTXShq45TLzU0x0c7oJL6dpZoOimRLCRcDxwIb0+LWSTgC2AP8BfB7YJmmgGpKXLZrfs+kmJiZ37S4lrNm4nQXL17nkYGZ9p5lxCCcCZ0fExwAkfYwkEUwArwQeIVnY5jxJ10XE17sVbJmWLpzX08VrssVyVq/ftnudBC+gY2b9pJmEcABwXe71R4C3Aq+KiK+l2y6T9HjgLGAgEkKvZUkgv2iOF9Axs37STJWRgN/mXt+e/rujcNxXgGNLiGloLVi+ruoQzMzqaraXUdR4/vvCMXcBczqOaIi5RGBm/azZuYy+Lekm4AaSxuQAHl/jOE9ZYWY2oJpJCEuBo0nWPT4d2CfdfrWkn5GsYXAj8LuuRDjk1mzc7oZlM+sLDRNCRHw2/1rS4STJYYwkUZwIvCI7vOwAh93KtVt3j1NYtmi+k4OZVabl6a8j4qfAT0m6ngIg6cnAMcDzygttdLgbqpn1g1KW0IyIe4B16cPa5EZnM6uS5zJqoBvTWsyeOYOD9ptV+nnNzDrhhNBAfoqLsmxZsZhN5y0q7XxmZmWoNCFIOlfStenEefdK+qakI6uMqWjpwnlsWbGYLSsW1z2m3WTRzUn1zMxaVXUJYRz4FPBC4CUkXVfXS3pSlUG1opPJ8bLSh5lZPyilUbldEXFK/rWkM4CdJF1Zv1lJUC3YsWpJR+/PJtibe86lJUVkZta+ShNCDfuSlFp+VWunpA3p07FeBdSO2TNntN1jaM3G7R6XYGaVqLrKqOhiYDNwTdWBdKKThuhseuz8wjrZGgpzz7nU6yiYWdf0TUKQdBHwIuAVEVHz63VEjEfEOEnS6FvNNETXU2t67PwaCl6Bzcy6pS+qjCR9FHgN8CcR4a+/BcXqp4nJXbvbHbJGbVctmVmnKi8hSLoYOA14SUT8pOp4Bo1LDGZWlkpLCJI+CZwBvBz4laSD010PRcRD1UU2WCYmd+1efMelBTNrV9UlhLNIehZ9G7gz93hPlUGVqdNxBs02IBcbos3MWlVpQogI1XlcUGVcZWp30Fpm5dqtLR2flRbcE8nMWlV1CWHoVVF945KCmbWjL3oZjbpmBrJlvYnyXVCnMzG5a/cgt4nJXe6NZGYNuYTQB+rNaXT+qUfsHuCW3cxbGduwcu3WKeMXVq7d6qokM6vLJYQ+UG9Oo2x7mbwqm5nV4xLCiPGqbGZWj0sIIygbszB71l5MPPo7wOMXzMwJYSTl2xUyWfdWJwWz0eUqI9tt5dqtHsNgNsKcEAZQN1dZ8xgGs9HlhDCA2ll6s5Xjs9HOLi2YjRYnhAHU6ngEaH0KDc+NZDZ6nBBGRLuNxe6majY63MvIGioOmMtGTgNe/9lsiDghWMvyVUlZCWLl2q2sXLvVcyaZDTBXGVlbsjaGWttXrt3KcReuryAqM+uEE4I9xo5VSzo+x90PPFpCJGbWS64ysq7LT8OdOWi/WWw6b1GFUZlZkUsIPZAfA9DNQWXd1EnctdZwcAnCrP84IfRANpAs3ztn0LQzGC7jrqtmg8FVRj3QjXUNei37GdZs3N7yOs/TqbWqG7g7q1kVXEKwlnRjwZ787Kur12/bvS3fvXXNxu2eTsOsyyovIUhaCLwHOBY4BPiriLik0qCGVDNrN/fSmo3bHxNPvdf5xLFy7VavFW3WBf1QQtgHuBl4B/BwxbEMtU7aAbqhlXmS6iWKbNyDSw9mnas8IUTE2og4LyK+Avx+umMlbZC0ARjrSXBDpp1J8Tox95xLHzPtRV6ZpRVPxmfWucqrjGyPfJVOVd/km/ncg/ab1bfdRqdLMlkDNuxZPjR/vV31ZKOu8hJCKyJiPCLGgc1Vx9INnXZP7SSJtPK5gzqgLN9Yfc+Dj05JHi5dmLmE0Fc66Z6a3czb6RI6e+aMnlYldduC5etqfttvVEVVb39xpLVLEzasnBBa0A9VOvVkN/R2EkIZg+Ua9WDqZQ+nrKG5k7EMtabbyJ8/O7fHS9gwGagqo6oNw4jjWmrdyFqdbqNRCaPMEkj2O2iknYbmbLzDyrVbp01g2Xmzz1i5dqt7ONnAqzwhSNpH0piksTSew9LXh1UdW1HWS2fLisVD/22wn5Nf9jtoVislk3qlAoDzTz1i2vPmE48H0tkgqjwhAH8MXJ8+ngD8Xfp8RZVBjbpRSn550yWPpQvnTVsyyb+31mhrs35XeRtCRGwAVHUc1r/a+Ya9Y9WSacdAtCpLBMsWzX9MA3O9dobi83y312KbQ719jeZ6yrrP1jqnWasqTwjWvmw8wEH7zao6lGl12gCf/4bdq8b8rHoof5OG2j3BiomnXgKrNf1Gdu5a+4rbay1dmk88+fc1k2zMipwQBli/jwfIf6vuRP6m1+65Wi1lZDfNdm6e9aqI6k2/UatnWNZQ3egctfavXr9tStz5pFLsHZUvYdQqbTiZjBYnhCHTTxPYdWNsQys3pPy16FU9fq0J+7LtvTIxuWvKjbxYfZVPNLWqtiBJUsWEVNz25H1nubpqyPRDo7KVqJkJ7IZhBbdm5EsT3UyS+WtYL/F0kpDa+R1l3Wa7+XNno73zEwzmE98LVq5n7jmXctyF64Hpe165V1Z/cAlhyGR13AuWr6s7iC6rm86eD6ulC+eVuphPPfkR4vVuwJ3cmLesWDztQLl+UayuuufBZL6rux949DHtLPmqq1rdd+s1qmfHek3u7nBCGFLT3fSrWsGtleqsfqr6aqQXiSf7nZXZc6ob8j2qGh1X75pNTO7iBSvXT5l8MH9u2JNkGq2yV6sNpFbPLbeXJJwQhlQVN/1GU3vkk1Sjm32xh40lepUoi912s99n/rN3rFoypSSa12mbTVa6aKRWYsl6W82etdeU8+QXViq+f813t085tl4JppiAau0b5ETihGBtqXXzb1QVlU9Sjb7pdvKNeNjbRZoZB1G2rIPAcReun9LVuV7irhXT7Jkz2GfvvXoydXq99pN616qYgKbr5VXvS02tHl6ZRuNMoD/GlDghWFtq3fy7USppdULBfpxqoxPFn7mZcRDdVKy3byVxZ0mlXqliUDSa46qWZseZZDqdnLFdTgjWll5VSbXaAF5GV9d2Shjd+KY+qMmtUXLI3wg7vW7duu5ln7NeaaKZ5JIvqXQ7SbjbqVWi2ZtuFXMqtXMT7saNu5s/8/mnHlFZ1Vr2O92xaknHCbyT615v1txOY8p3nc2609bS6qSLvZgfywnBKpGfTbU4i2jV2rkJD1pDYv6m3IvrP13y6SQxdbqgVDcSeXG9jDJKG8U2kW5VuTkhWCXK+Obfb43HxXg6XdK0W2q1S9SSNRx3OldWo6qv/JeDds9f6/l0sr+7biXyWu0Cg8AJwQZWv9Wv50eJd7IudqvvbeVG2sq5N523iB2rljQ1AGy6GBol/fyXg3aUuRZ5v33J6DU3KltfaGd50nYGhHWzm2atb5zNxJfvttvON9Z8w/s+e+/FQ48kXRdr/ZzdWju7lTEmZSte91b/Jmp1XCj+neR/R82ev9u9v+aecyk7Vi0p9ZxOCNYX2p1Oo9VE0sp/6F7p9CZdr+qjl91R8zFMN21KJ7qVzGtdv+LfYycJp55+HI3vhGB9od363FYTSa/mN+oH7U4V0ulNvFtzZW1ZsbipJFfGjbaVv8dswN1Dj/yupc/NvgjkB6flz1NFwnBCsIHW6D9umTe6dhQ/v5f/wVspDZV5E++ksbaMa9TrUmC+hDddwqr3tzjd9SpOH9LtQWpuVLahVqvBsdiI2M1GxU4bPDvRyo2jX9bQbjR9ezO/qyrjrxdT9vvv5G+hF78blxBsqDVTPwyPXSqzW59fa0I026PRfFe9amtq91z1Gvizb/at3tB7vUxuXyQESWcBZwNPAbYAyyLiu9VGZcOq1n/MXo6CLk693E1VV5l1olbsvWpravdcZY9t6PWaD5UnBEmvBi4GzgK+l/57maTnRMT/qzQ4s5L1elryQV4MqV/aNbp5rn5TeUIA3gVcEhFr0tdvk7QYeDNwbnVhmQ2+Qb55ldU4PWgloypV2qgsaSZwLHBFYdcVwAtrHL9B0gZgrPvRmdmgqrIxf5BVXUI4AJgB3F3YfjfgBVPNrC2DXDKq0kB1O42I8YgYBzZXHYuZ2bCpOiHcB+wCDipsPwi4q/fhmJmNrkoTQkRMAtcBJxd2nQx8v/cRmZmNrqrbEAAuAj4v6YfA1cCbgEOAT1calZnZiKk8IUTElyX9IfA+koFpNwOnRsRt1UZmZjZaKk8IABHxKeBTVcdhZjbKFBFVx9AySb+cM2fOoWNjHo5gZtaszZs3s3Pnztsj4qm19g9qQrgeOBC4tY23Z1nEXVfL42taLl/P8vmaJp4J3BsRR9faOZAJoRPpSGfS8QxWAl/Tcvl6ls/XtDlVj0MwM7M+4YRgZmbACFYZmZlZbS4hmJkZ4IRgZmYpJwQzMwOcEMzMLDVSCUHSWZJ+LukRSddJenHVMfUjSRdIisLjrtx+pcfcIenhdCW7BYVz/IGkz0vamT4+L2n/3v801ZC0UNI3JN2eXr8zC/tLuYaSjpJ0VXqO2yUtl6Qe/Ig918Q1vaTG3+0PCsfMkvRxSfdJmkjP99TCMYdJ+ma6/z5JH0tXdxx6I5MQJL0auBi4EDiaZHrtyyQdVmlg/esWkskGs8dRuX3vBd4NvA14PnAP8N+S9s0d86/AMcDi9HEM8Pnuh9039iGZqPEdwMM19nd8DSXtB/w3yQqDz08/62ySdcqHUaNrCrCeqX+3pxb2rwZeAZwGvBjYD/iWpBkA6b+XAvum+08DXgl8pMwfpG9FxEg8gE3AmsK2nwIfqjq2fnsAFwA319kn4E7g/Ny2JwAPAv87fX0EEMCJuWNelG57VtU/XwXX8yHgzLKvIfBm4AHgCblj3gfcTtqlfFgfxWuabrsE+NY075kDTAKvzW17GvB74JT09Z+mr5+WO+Z04BFgv6p/7m4/RqKEkBb3jgWuKOy6Anhh7yMaCPPS6oyfS/qSpGyB2mcAB5O7lhHxMLCRPdfyBJL/sPlFjq4GJvD1hvKu4QnAd9P3Zi4nWU9kbjcCHwAvknSPpG2S1kh6cm7fscDjmXrdfwFsZeo13Zpuz1wOzErfP9RGIiEABwAzSIrWeXeT/Me0qTYBZ5JUUywluUbfT9etyK7XdNfyYJIJtHaPekyf34OvN5R3DQ+uc478Z4ySdcDrgJeSVMe9APiOpFnp/oNJluy9r/C+4nUvXtNsqd+hv6Z9sR6C9ZeIuCz/Om2Y2w68HvhBzTeZVSwivpR7eZOk64DbgCXAV6uJarCMSgkhy/AHFbYfBNz12MMtLyIeArYAh7Pnek13Le8CDsz3dkmfPxlfbyjvGt5V5xz5zxhZEXEH8EuSv1tIrskMkhqDvOJ1L17TrIZh6K/pSCSEiJgErgNOLuw6mal1tFaDpL2BZ5M0hP6c5D/GyYX9L2bPtbyGpEfICbnTnADMxtcbyruG1wAvTt+bORm4A9jRjcAHiaQDgENJ/m4huQf8lqnX/akkDfj5a3pEoSvqycCj6fuHW9Wt2r16AK8m6WHwRpI/gItJGu2eXnVs/fYAPgycRNL4eRzwLZLeLE9P9/8tsBP4C+BI4EskN6F9c+e4DLiJ5CZ2Qvr8m1X/bD28hvuQLMoyBvwGWJ4+P6ysa0jSa+au9L1Hpud6AHh31T9/r69puu/D6XWaC4yT3Nx/Wbim/5BuW0TS/fxKkkVzZqT7Z6TX+Tvp/kUkvbY+XvXP35NrXHUAPf6DOovkm1OW7RdWHVM/PnI3p8n0P8N/As/J7RdJ19Q7SbrjXQUcWTjHHwBfSG9QD6TP96/6Z+vhNRwn6SJafFxS5jUkGR+yMT3HncAHGNIup9NdU5Juu5eTNLpPkrQdXEKu+2h6jlnAx4H706TyzRrHHEbyJeg36XEfA2ZV/fP34uHpr83MDBiRNgQzM2vMCcHMzAAnBDMzSzkhmJkZ4IRgZmYpJwQzMwOcEMzMLOWEYGZmgBOCWdPSpRS/lXudLTXak1mDJS2TdJMk/7+1rvAfllkTJP0R8CaS6Saq8o/AgSTTkJuVzgnBrDnLgBsi4kdVBRDJymj/ArynqhhsuDkh2EiRdGJazVPr8ek675lFsq7uvzZx/sWSHpL0iaxqJ1e1dJSkKyX9RtKdklYUq38kPU/S1yTdL+lhSbdIOjd3yJeA50jyUqRWOq+YZqPmJ0xdYwDgHJLlQr9c5z3HA/sD353uxJJeB3wGWBERf1/jkK8D/wf4EHAK8H6SBd0vSN//AmADcCvwTvYs7vLc3Dk2Aw+m8XptCSuVE4KNlIi4n2RKYwAkrSC5Ob88Iq6s87bjSaZZvrHeeSW9F1gJvDkiPlPnsDURsSp9foWk/YB3S1odEb8mmc//fuD4iPhNetx3CvH/XtINaUxmpXKVkY0sSatIFmP/s4hYN82hhwAPRLLyXi0fBf4OeOU0yQDg3wuvv0SysMuRkp4InAh8MZcM6rk3jcmsVE4INpIkXQS8FTg1ItY3OHxvkkWV6jkNuBlodJ6767w+lGQxnMeRVBM18jDJgjBmpXJCsJEj6eMkS6meEhFXNfGW+0naEOp5KckqW5dJ2mea44qLt2evbwd+RdKecGgT8TwJuK+J48xa4oRgI0OJTwNnACdHxNVNvvUnwMzCwut5W0iWdzyc6ZPCqwqvX0OyrvdNaTXR94DTJTX69v8M4JZmAjdrhRuVbZRcTFIyeAdJfsg3zP44Ih6o876N6b8voE6VTkRslTROsmj75ZIWR8SDhcOWpt1MryVpyH4jcEFE7Ez3v4dkbeVrJH0k/ax5wFhEvI0k6P2B+SQN0GalcgnBRoIkkYzwnQF8Arim8Di43nsjYgfwQ+DPpvuMiLgFOAl4Ont6EeW9DDgZ+AbJuIa/Bz6Ye/+1JA3LvyBZCH4tcDZTk9ASkkXkvzZdLGbtUERUHYNZ35N0JkkJ4ylN9AIqvvcC4APA4yPidx3GcRlwX0Sc0cl5zGpxCcGsOV8A7gDOqioASWPAS0i6uJqVzgnBrAnpN/u/AloqHZTsYODMiLi1whhsiLnKyMzMAJcQzMws5YRgZmaAE4KZmaWcEMzMDHBCMDOzlBOCmZkBTghmZpb6/yEdtJbIAeRGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "s.set_coherence_r0(None)\n",
    "s.domain.create_box_array(Lmax, i, s.coherence_func, r0=0)\n",
    "plt.step(s.domain.rcen, s.domain.B * 1e6, where=\"mid\")\n",
    "plt.ylabel(\"$B_\\perp (\\mu G)$\")\n",
    "plt.xlabel(\"$z$ (kpc)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The power-law function can also be customised using the `set_coherence_pl` function. For example, if we want a PDF of the form $p(\\Delta z) \\propto \\Delta z^{-3}$ chosen between 50 and 100 kpc, we can do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$z$ (kpc)')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAETCAYAAAA23nEoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAW7UlEQVR4nO3de7RkZXnn8e/PVkAbbGJAO6CoREjMYGzQUQgInQSWRiajmTheRlCSkYwhUfGuyZigxsgfaiCaiRk0gagJmTHR8QJqUBoygkQJLZcgiIhLuSgw2AgirfjMH3ufdHVxLlV1qs6uc+r7WWuvrtq3euo9p+s576XeN1WFJEkP6DoASdJ0MCFIkgATgiSpZUKQJAHwwK4DGEWSy4C9geu6jkWSVpHHAbdW1cHzHVyVCQHYe8OGDftu2rRp364DkaTVYuvWrWzbtm3B46s1IVy3adOmfbds2dJ1HJK0amzevJkLLrhgwZYV+xAkSYAJQZLUMiFIkgATgiSpZUKQJAEdJ4QkpySpvu2WLmOSpFk1DcNOrwE29zy/r6M4JGmmTUNC+FFVrVit4DFv+OSCx9bvso6Tjz6QE4/cf6XCkaSpMQ19CPsnuSnJ15OcnWTBT+MkW5JsATZNIpC7t9/HaeddO4lbS9LU6zohXAKcADwDOBHYCFyU5Ce7Cuju7bZYSZpNnTYZVdW5vc+TfAG4Hngx8K55zt/cnrcFOGqU17zh1GPn3b9YU5IkzYKuawg7qaq7gKuAA7qORZJmzVQlhCS7AT8L3Nx1LJI0a7r+HsI7khyV5LFJngp8GFgPnNVlXJI0i7oedvpI4G+BvYBbgS8Ah1bVNzqNSpJmUNedys/v8vUlSTtMVR+CJKk7JgRJEmBCkCS1TAiSJMCEIElqmRAkSYAJQZLUMiFIkgATgiSpZUKQJAEmBElSy4QgSQJMCJKklglBkgSYECRJLROCJAkwIUiSWiYESRJgQpAktUwIkiTAhCBJapkQJEmACUGS1DIhSJIAE4IkqWVCkCQBJgRJUsuEIEkCTAiSpJYJQZIEmBAkSa2pSghJ3pikkryn61gkadZMTUJIcijwW8DlXcciSbNoKhJCkg3Ah4DfBO7oOBxJmklTkRCA/wl8uKrOX+ykJFuSbAE2rUhUkjRDHth1AElOBB4HHNd1LJI0yzpNCEl+Bvhj4Iiq+uFS51fV5va6LcBREw1OkmZM1zWEw4C9gKuSzO1bBxyZ5KXA+qq6t6vgJGmWdJ0QPgp8qW/fXwFfpak5bF/xiCRpRnWaEKrqu8B3e/cluRv4f1V1ZTdRSdJsmpZRRpKkjnXdZHQ/cx3HkqSVZQ1BkgSYECRJLROCJAkwIUiSWiYESRJgQpAktUwIkiTAhCBJapkQJEmACUGS1DIhSJIAE4IkqWVCkCQBJgRJUsuEIEkCTAiSpJYJQZIEmBAkSS0TgiQJMCFIklomBEkSYEKQJLVMCJIkAB44zMlJDgWeARwK7AM8GLgNuAa4APhoVd0x7iAlSZM3UA0hyYuTXAFcBLwSeAjwVeAS4A7gqcD7gBuTnJnksROKV5I0IUvWEJJcDuwN/DXwImBrVdU8520A/gPwQuBfk5xQVX835nglSRMySJPR+4G/qKofLHZSVW0DPgR8KMkTgY1jiE+StEKWTAhVdfqwN62qLwNfHikiSVInHGUkSQIGSAhJ9kjy1iRH9e1/8OTCkiSttEFqCC8BXgPcPrcjyTpgW5KtSc5I8ltJnjCpICVJkzdIQvg14KyqurJv/wOB7cCzgfcC/5Lkp4d58SS/k+TyJHe228VJjh3mHpKk8RgkIRwEnDvP/gJeWlV7AwcAV9AMSx3Gt4DXA4cATwY+B3w0yc8PeR9J0jINkhB2p6e5qEfmHlTV14APAr84zItX1f+pqnOr6rqquraqfh/4HnDYMPeRJC3fIN9DuB14VO+OqrovySOB7/bs/grw+FEDafsl/jNNArpogXO2tA83jfo6kqT5DVJDuBh4bv/Oqrqpqr7fs+se4KHDBpDkCUnuAu6l6Yv4taq6Ytj7SJKWZ5CE8GfAs5Icv8R5BwLbRojhGpq/+J8K/DlwVpKD5juxqjZX1WZg6wivI0laxJIJoao+C5wO/FWStyfZvf+cdt+rgM8PG0BVbW/7EC6tqjfSfNi/ctj7SJKWZ6Dpr6vqlUl+CLwW+J0k5wCX0dQIHk0zod0jgOPGENMDgF3HcB9J0hAGXg+hql6X5O+B3wOexc79Ct+gafv/4jAvnuRU4JPAN4E9gP8CbAb8LoIkrbChFsipqkto+hPWA/sDewG3AlfNNyX2ADbSDFfdSFPbuBz4lar69Aj3kiQtw7Arph0P/Atwdf9IoCS7VtW9w9yvqk4Y5nxJ0uQMlRCAs2i+oXxvkqto+hEuo+kIfnKS51XVEWOOUZK0AoZNCA8DDm63Q4DDgd9kx7eW7xxfaJKklTRsH8J3gfPbDYAkewMvpxkq+oKxRidJWjHD1hDup6puBd6U5GE0CeFTy45KkrTixrli2mdopsKWJK1Cw44yehNwKXBZVd3cd/jhwN3jCmzanXHh9Zx23rXcvf2+ga9Zv8s6Tj76QE48cv8JRiZJoxm2yei1NLORVpLv0AxBvQy4j2ZltdePN7zpNWwyALh7+32cdt61JgRJU2nYhLCBZhK7Q9rtScBJwJ7t8XclOY6mFvGlqvqHcQU6bYZNBsu9TpImbdhRRkUzO+k1wN/O7U+yP01ymEsUJ9LUFtaNLdIpdsOpS8+08Zg3fHIFIpGk0S17lBFAVV0PXA/877l9SfYbx70lSStjyVFGST6W5OBBb5hktySvAp65rMgkSStqkGGnNwBfSHJJkpcnOSTJTjWLJPskeXaS9wM3A/+VpsNZkrRKLNlkVFUvT3I6cDJwCk3HciW5k2bZyz2BXWimr/jn9rwPVpW9p5K0igy6QM7XgJcleTVwGM1yl/sAuwG3A18BLqyqb0wqUEnSZA07ymg7cEG7SZLWkHFOXSFJWsVMCJIkYEIJIcmmSdxXkjQ5k6oh/I8J3VeSNCGTSghZ+hRJ0jQZeeqKJPsAT57vEDsmu5MkrRLLmcvoIcCjFji2yzLuK0nqwMgJoaquA66b71iSF44ckSSpE/YhSJKAySWEt07ovpKkCZlIQqiqcyZxX0nS5PhNZUkSYEKQJLWWHGWU5DHAWcDBNNNcv6aqLkyyHng68GDgoqr6+gTjXFMGXV95/S7rOPnoAznxyP0nHJEkDVZDeBdwKLClPf+cJIcBV9GsofwB4NokdiQvYv0u64a+5u7t93HaeddOIBpJur9BEsLhwGur6j9W1ZOBv6RJBPcCzwGOpUkKv5fk2ROLdJU7+egDR04KkrQSBvli2l7ApT3P3wn8LvDcqvpIu+/cJA8CTgI+Ot4Q14YTj9x/qKafQZuVJGlcBqkhBPhhz/Mb239v6Dvvw8CThnnxJG9M8sUkdya5NcnHkxw0zD0kSeMx6Cijmufxj/vOuQXYMOTrb6aZKvsXgF8CfgScl+RhQ95HkrRMg85l9NkkVwBfpulMLuBB85w31JQVVfX0nS5Ojge20fRbfPx+N0+2tA9dgEeSxmyQhHAizZDTTcBxwO7t/s8n+Rpwebv9aAzx7EFTa7ljDPeSJA1hyYRQVe/vfZ7kAJrksIkmURwO/Prc6cuM53RgK3DxArFsbmPYAhy1zNeSJPUYevrrqvoq8FWaoacAJHk4cAjwxFEDSfIu4AjgiKpyrKUkrbDlLJDzb6rqO8Cn2m1oSf4EeD7wi1V1/ThikiQNZywJYTmSnA48jyYZfKXreCRpVnWaEJL8GXA88GzgjiQb20N3VdVd3UUmSbOn6xrCSe2/n+3b/2bglJUNZQe/JSxpFnWaEKpqapbaXL/LupHmDRplfiJJmkauh9AaZfK5uempJWkt6LrJaGoMO/mcJK011hAkSYAJQZLUMiFIkgD7EFaFQYbBuv6ypOWyhjClhh3x5PrLkpbLhDClRhkG6/rLkpbDJqMpNcwwWL9ZLWkcrCFIkgATgiSpZUKQJAH2IcycMy68ntPOu3aoDmiHtEqzwRrCjBk2GYBDWqVZYUKYMaMOTXVIq7T22WQ0w2449dglz3FIqzQ7TAhrjB/gkkZlk9EaMMqqba70JqmfCWENGHaaC1d6kzQfm4zWgJVa7W3Y5iiHq0qrizUELWo5TUsOV5VWFxOCFjXKrKu9HK4qrR42GWlRozZHOdpJWn2sIUiSABOCJKllQpAkASYESVLLhCBJAkwIkqSWCUGSBJgQJEmtzr+YluRI4DXAk4B9gN+oqjM7DUqdGmWZT3DuJGm5pqGGsDtwJfAK4J6OY9EUGCUZgHMnScvVeQ2hqs4BzgFIcuZi5ybZ0j7cNNmo1KXlzH/k3EnS6DpPCNJiBlnmE3aeO2nQeZRsYpJ2Ng1NRgOrqs1VtRnY2nUsmi6jzMhqE5O0s1WVEKSFjDpNt01M0g42GWlNGHaabpuYpPuzhqCZZBOTdH+dJ4QkuyfZlGRTG89+7fP9uo5Na5dNTNL9TUOT0ZOB83uev7ndzgJO6CIgrX3LaWKS1qrOE0JVbQHSdRySNOs6bzKSJE2HzmsI0lrkfExajUwI0gQsZz6mt51zNW875+qxxGGC0TBsMpImYFpGIzlUVsOwhiBN2KDzMY3azLSUaUlOmn4mBGlKDDsUdikOldWwTAjSDFgqOdjXILAPQVqzhvkmtn0NAmsI0pp18tEHDtUnMUxfw3L6O6yNTC8TgrRGDdonMUpfw3I6v0cdWmsimTybjCQNrYuRSzZrTZ41BEnLMuiwWlj+0FqH0E6WCUHSihl1aK1DaFeGCUHSqjJqchilD2JSXxYcl2FqZ4OwD0HS1BtlMaN+o/RBTHMymAQTgqSpN+oKd/2G/XCfpWQANhlJQ7M9e+Utd1qPcfzMxt08M42sIUgDGPWv03H8VSutFGsI0gCG/dYv7OjE1HSxhrcwE4I0gHHPRKqVtX6XdcvqD5iVmp4JQRPnX2Srx1r9WY1Sw5szSzU9E4ImYrl/kc3dQ5O3nJ/VavkZWcMbjJ3KmojlDhOcpb/Kujbqz8qf0dpjDUET4V9kq4c/K82xhiBJAkwIkqSWCUGSBJgQJEktE4IkCTAhSJJaqaquYxhakm9t2LBh302bNnUdiiStGlu3bmXbtm03VtUj5zu+WhPCZcDewHVLnDqXMbZONqKZZNlOhuU6OZYtPA64taoOnu/gqkwIg0qyBaCqNncbydpj2U6G5To5lu3S7EOQJAEmBElSa003GUmSBmcNQZIEmBAkSS0TgiQJMCFIklprNiEkOSnJ15P8IMmlSZ7WdUzTLMkpSapvu6XneNpzbkpyT5ItSf5d3z1+IskHkmxrtw8k2XPl3023khyZ5GNJbmzL8YS+42MpyyRPSHJBe48bk/xBkqzAW+zMAGV75jy/x1/oO2fXJO9OcluSu9v7PbLvnP2SfLw9fluSP02yywq8xU6tyYSQ5HnA6cAfAwcDFwHnJtmv08Cm3zXAT/VsT+g59jrg1cDLgH8PfAf4xyR79JzzN8AhwDPa7RDgA5MPe+rsDlwJvAK4Z57jyy7LJA8F/hH4dnuPVwCvBV415vcybZYqW4Dz2Pn3+Jl9x08Dfh14AfA04KHAJ5KsA2j//SSwR3v8BcBzgHeO841MpapacxtwCXBG376vAm/vOrZp3YBTgCsXOBbgZuD3e/Y9GPge8N/a548HCji855wj2n0/0/X767Bc7wJOGHdZAr8N3Ak8uOec/w7cSDucfK1v/WXb7jsT+MQi12wAtgMv7Nn3KODHwNPb57/SPn9UzznHAT8AHtr1+57ktuZqCG217knAZ/oOfQb4hZWPaFXZv23G+HqSs5PMLbT7WGAjPWVaVfcAF7KjTA+j+Q96Uc/9Pg/cjeXea1xleRjwT+21cz4N7AM8ZhKBryJHJPlOkmuTnJHk4T3HngQ8iJ3L/5vA1exctle3++d8Gti1vX7NWnMJAdgLWEdTle71bZr/iJrfJcAJNM0TJ9KU1UVJfpId5bZYmW6kmTTr377p2D7+DpZ7r3GV5cYF7tH7GrPoU8CLgF+maZZ7CvC5JLu2xzcC9wG39V3XX/79ZXtbe92aLtsHdh2ApkNVndv7vO2Iux54MfCFeS+SpkxVnd3z9IoklwLfAI4F/qGbqFaPtVhDmMvkj+jb/wjglvufrvlU1V3AVcAB7Ci3xcr0FmDv3lEu7eOHY7n3GldZ3rLAPXpfY+ZV1U3At2h+j6Epm3U0LQm9+su/v2znWh7WdNmuuYRQVduBS4Fj+g4dw85tslpEkt2An6XpAP06zX+EY/qOP40dZXoxzQiQw3pucxiwHsu917jK8mLgae21c44BbgJumETgq1GSvYB9aX6Pofls+CE7l/8jaTrye8v28X1DUY8B7m2vX7u67tWexAY8j2YkwUtoftCn03TSPbrr2KZ1A94BHEXT6flU4BM0o1ge3R5/PbAN+E/AQcDZNB8+e/Tc41zgCpoPr8Paxx/v+r11UJa70yzGsgn4PvAH7eP9xlWWNKNlbmmvPai9153Aq7t+/12VbXvsHW15PQbYTPPh/q2+sv3zdt/RNMPSz6dZNGdde3xdW96fa48fTTN6691dv/+Jl2/XAUzwF+ckmr+U5rL6kV3HNM1bz4fS9vaX/++Bn+s5HpqhqTfTDL+7ADio7x4/AXyw/WC6s328Z9fvrYOy3EwzRLR/O3OcZUnzPZEL23vcDPwha3zI6WJlSzN899M0ne/bafoOzqRn+Gh7j12BdwO3t0nl4/Ocsx/NH0Xfb8/7U2DXrt//pDenv5YkAWuwD0GSNBoTgiQJMCFIklomBEkSYEKQJLVMCJIkwIQgSWqZECRJgAlBGli7jOInep7PLTu6IrMGJzk5yRVJ/H+rifAXSxpAkp8GXkoz5URX/gLYm2ZKcmnsTAjSYE4GvlxVX+oqgGpWR/tr4DVdxaC1zYSgmZLk8LaZZ77tvQtcsyvNmrp/M8D9n5HkriTvmWva6WlaekKS85N8P8nNSd7S3/yT5IlJPpLk9iT3JLkmyRt7Tjkb+LkkLkuqsXPFNM2ar7DzOgMAb6BZOvTvFrjmUGBP4J8Wu3GSFwHvA95SVX80zykfBf4SeDvwdOBNNIu5n9Je/xRgC3Ad8Ep2LOzy8z332Ap8r43XdSY0ViYEzZSqup1mOmMAkryF5sP52VV1/gKXHUozxfLlC903yeuAtwG/XVXvW+C0M6rq1PbxZ5I8FHh1ktOq6rs0c/nfDhxaVd9vz/tcX/w/TvLlNiZprGwy0sxKcirNQuy/WlWfWuTUfYA7q1mNbz5/ArwZeM4iyQDgf/U9P5tmUZeDkjwEOBz4UE8yWMitbUzSWJkQNJOSvAv4XeCZVXXeEqfvRrPQ0kJeAFwJLHWfby/wfF+aBXEeQNNMtJR7aBaDkcbKhKCZk+TdNMurPr2qLhjgkttp+hAW8ss0K2ydm2T3Rc7rX7h97vmNwB00/Qn7DhDPw4DbBjhPGooJQTMjjfcCxwPHVNXnB7z0K8AufYuu97qKZmnHA1g8KTy37/nzadb6vqJtJvq/wHFJlvrr/7HANYMELg3DTmXNktNpagavoMkPvR2z/1pVdy5w3YXtv09hgSadqro6yWaaBds/neQZVfW9vtNObIeZfpGmI/slwClVta09/hqa9ZUvTvLO9rX2BzZV1ctogt4TOJCmA1oaK2sImglJQvMN33XAe4CL+7aNC11bVTcA/wz86mKvUVXXAEcBj2bHKKJezwKOAT5G872GPwLe2nP9F2k6lr9Jswj8OcBr2TkJHUuzgPxHFotFGkWqqusYpKmX5ASaGsZPDTAKqP/aU4A/BB5UVT9aZhznArdV1fHLuY80H2sI0mA+CNwEnNRVAEk2Ab9EM8RVGjsTgjSA9i/73wCGqh2M2UbghKq6rsMYtIbZZCRJAqwhSJJaJgRJEmBCkCS1TAiSJMCEIElqmRAkSYAJQZLU+v8VJ1xtGcNFKAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "s.set_coherence_pl(n=3, xmin=50.0, xmax=100.0, r0=None)\n",
    "s.domain.create_box_array(Lmax, i, s.coherence_func, r0=0)\n",
    "plt.step(s.domain.rcen, s.domain.B * 1e6, where=\"mid\")\n",
    "plt.ylabel(\"$B_\\perp (\\mu G)$\")\n",
    "plt.xlabel(\"$z$ (kpc)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more control over the distribution of cell sizes, any function that chooses random numbers can be passed to the coherence_func object - a convenient way to do this is using a [scipy stats function](https://docs.scipy.org/doc/scipy/reference/stats.html) like rayleigh using the inbuilt `rvs` method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$z$ (kpc)')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAETCAYAAAA23nEoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAW/UlEQVR4nO3dfZQldX3n8ffHUVCHp7iiCIo4C2NUjAMmCgGHiRkOLGyOJBoNRzAkcbKGqBkfMJJEZTVG/lBkfFqz6AZFE8wm0SPyoEGZISIisgwC8ijCJjzIw+KMDOCY8bd/VDVz56Zvd93ue2/d7vt+nXMPfavqVn+76KlP1+/3q/qllIIkSY9ruwBJ0ngwECRJgIEgSaoZCJIkAB7fdgFzkeRqYE/g1rZrkaQFZH/gvlLKQdOtXJCBAOy5++6777NixYp92i5EkhaKjRs3smnTpp7rF2og3LpixYp91q9f33YdkrRgrFq1ig0bNvRsWbEPQZIEGAiSpJqBIEkCDARJUs1AkCQBBoIkqWYgSJKAhXsfwpzt987zp12+dKclrF29nDUrl424IkkaD14h1LZs3caZF9/cdhmS1BoDocOWrdvaLkGSWjNxTUa3n37sf1jWqxlJkiaJVwiSJMBAkCTVDARJEmAgSJJqBoIkCTAQJEk1A0GSBBgIkqSagSBJAgwESVLNQJAkAQaCJKlmIEiSAANBklQzECRJgIEgSaoZCJIkwECQJNUMBEkSYCBIkmoGgiQJMBAkSTUDQZIEGAiSpJqBIEkCDARJUs1AkCQBBoIkqWYgSJKAMQiEJM9I8pkk9yV5NMn3kxzRdl2SNGke3+Y3T7IHcBnwTeBY4D5gGXBvm3VJ0iRqNRCAdwB3l1Je17Hsh20VI0mTrO0mo+OAK5J8Icm9STYmeWOSTLdxkvVJ1gMrRlqlJE2AtgNhGXAycBtwFLAOOB344zaLkqRJ1HaT0eOA75ZSTq3fX53kAKpA+Fj3xqWUVVBdKQB2PEvSALV9hXA38P2uZTcA+7ZQiyRNtLYD4TLguV3LlgN3tFCLJE20tgPhw8AhSf48yf5Jfht4M/DxluuSpInTaiCUUq6kGmn0auA64P3Au4BPtFmXJE2itjuVKaWcD5zfdh2SNOnabjKSJI0JA0GSBBgIkqSagSBJAgwESVLNQJAkAQaCJKlmIEiSAANBklQzECRJgIEgSaoZCJIkwECQJNUMBEkSYCBIkmoGgiQJMBAkSTUDQZIEGAiSpJqBIEkCDARJUs1AkCQBBoIkqWYgSJIAA0GSVDMQJEkAPL7tAsbNfu88v+e6pTstYe3q5axZuWyEFUnSaHiFQHWib2LL1m2cefHNQ65GktphIABrVy/vKxQkaTGyyQhYs3LZrM1AMzUlSdJi0FcgJDkEOBo4BNgbeBJwP3ATsAH4UinlwUEXKUkavkZNRkl+N8m1wLeAtwBPBm4BrgAeBF4KfAq4M8nZSZ4zpHolSUMy6xVCku8BewKfBV4HbCyllGm22x34r8Brge8nOamU8oUB1ytJGpImTUafBv66lPLoTBuVUjYBnwc+n+RFwF4DqE+SNCKzBkIpZV2/Oy2lXANcM6eKJEmtcNipJAloEAhJdk3yviRHdC1/0vDKkiSNWpMrhNcDbwcemFqQZAmwKcnGJGcl+cMkLxxWkZKk4WvSqfybwGdKKddN89mtwHHAHwDbkvxiKeUHA65RkjQCTa4QDgQunGZ5Ad5QStkTOAC4lmpYqiRpAWoSCLvQ0VzUIVNf1FcFnwN+bUB1SZJGrEmT0QPAszoXlFK2JXkm8OOOxTcCzxtgbZKkEWpyhXA58OruhaWUu0opD3csegTYbVCFSZJGq0kgfBx4RZITZ9luObBp/iVJktowayCUUr4OrAP+JskHkuzSvU297K3AZYMvUZI0Co0ef11KeUuSnwGnAH+c5ALgaqorgmdTPdDu6cAJ8ykmyanAXwEfL6W8cT77kiT1p/F8CKWUdyT5R+DPgFewY7/CHcBvllKunGsh9VwLfwh8b677kCTNXV8T5JRSrqDqT1gKLAOeCtwHXD/dI7Gbqh+d/Xng94H3zLDd+vrLFXP9XpKk6fX1cLskJyZ5AfBIKeXaUsolpZTrSiklyc7zqON/Av9QSrlkHvuQJM1Dv3Mqf4bqDuWfJrmeqh/hamAj8MtJXlNKObyfHSZZA+xPg/6HUsqq+jPrgSNm3FiS1Jd+A+EpwEH162DgMKpmnqm7ljf3s7Mkz6XqRD68lPKzPmuRJA1Qv30IPwYuqV8AJNkTeDPVXMvH9/n9D6Xqh7g+eexJGEuAlUneACwtpfy0z31Kkuag3yuE/6CUch/wriRPoQqEi/r4+JeA73Yt+xvgFqorh63zrU+S1My8A6HD14DP9vOB+oqj83lIJNkC/L9pHrctSRqifkcZvSvJMUmeMc3qpwFbBlOWJGnU+r1COIXqcdglyb3A/6EaZbSNama1P51vQVMjiSRJo9VvIOxO9RC7g+vXi4GTgT3q9WckOQG4CvhuKeWfBlWoJGm4+h1lVICb6tffTS1PsowqHKaCYg3V1cKSgVUqSRqqgXQql1JuA24D/vfUsiT7DmLfkqTRmLVTOcmXkxzUdIdJnpjkrcAx86pMkjRSTUYZ3Q58O8kVSd6c5OAkO1xZJNk7yXFJPg3cDfwBVYezJGmBmLXJqJTy5iTrgLXAaVQdyyXJZuCnVB3KO1E9vuI79XafK6VsG1bRkqTBazpBzg+ANyV5G9XjJl4K7A08EXgAuBG4tJRyx7AKlSQNV7+jjLYCG+qXJGkR6etOZUnS4mUgSJKAIQVCEqe4lKQFZlhXCJ8Y0n4lSUMyrEDI7JtIksbJnB9dkWRv4JenW8X2h91JkhaI+TzL6MnAs3qs22ke+5UktWDOgVBKuRW4dbp1SV4754okSa2wD0GSBAwvEN43pP1KkoZkKIFQSrlgGPuVJA2PdypLkgADQZJUazJj2n5JNiTZnOQ7SVbWy5cm+a0kr03ynOGXKkkapiZXCGcAhwDr6+0vSHIocD3VHMrnADcnsSNZkhawJvchHAacUkr5CECSj1AFwRbgVcCjwG8Df5bkqlLKl4ZVrCRpeJoEwlOBqzrefwh4I/DqUsoX62UXJnkCcDJgIEjSAtSkySjAzzre31n/9/au7f4BePEAapIktaDpKKMyzdc/79rmHmD3eVckSWpF02cZfT3JtcA1VJ3JBXjCNNv5yApJWqCaBMIa4CBgBXACsEu9/LIkPwC+V7/+fSgVSpJGYtZAKKV8uvN9kgOowmEFVVAcBrxyavNBFyhJGo2+H39dSrkFuIVq6CkASZ4GHAy8aHClSZJGaT4T5DymlHIvcFH9kiQtQD7LSJIEGAiSpNpAmowEZ116G2defDNbtm6bcbulOy1h7erlrFm5bESVSVIzXiEMSJMwANiydRtnXnzzCCqSpP4YCAPSJAzmsq0kjYpNRkNw++nHTrt8v3eeP+JKJKk5rxAkSYCBIEmqGQiSJMBAkCTVDARJEuAoo9bMNOLIm9cktaHVK4Qkpya5MsnmJPclOS/JgW3WNExLd1rSaDtvXpPUhrabjFYBnwB+FXg51SQ7Fyd5SptFDcva1cv7CgVJGqVWm4xKKUd1vk9yIrCJatKd81opaojWrFw2azOQN69Jasu49SHsSnXV8uB0K5Osr79cMaqCJGlSjFsgrAM2Ape3Xcg4sONZ0ii13YfwmCRnAIcDryylTNuAXkpZVUpZRRUai5Idz5LaMhaBkOTDwPHAy0spt7VdT5vseJbUltabjJKsA14D/Fop5ca262mbHc+S2tJqICT5OHAicBzwYJK96lUPlVIeaq8ySZo8bTcZnUw1sujrwN0dr7e3WZQkTaK270NIm99fkrRd21cIkqQxYSBIkoAxGGWk+fHmNUmD4hXCAuTNa5KGwSuEOWj7PoC1q5dz5sU3N7oxzZvXJDVlIDS0dKcljU6uTf96nw9vXpM0DDYZNdTkkRJTbfaStBB5hdBQk7/KJWkhMxAm1FmX3taoH8KRStLksMloQvXTKe1IJWkyGAgTqp/RR45UkiaDTUbi9tOPnXZ550glb4CTFj+vENSTN8BJk8VAUE/O3iZNFpuM1JM3wEmTxSsESRJgIEiSagaCJAmwD0EDNJ+hqd45LbXPQNC8NH0K7NTQ1F4n8n7unH7/BTfw/gtu6LvWKYaKND2bjDQvgxqaOsphq943IU3PKwTNyzCGpva6c7pps1ITW7Zu8+5rqYuBoAVjEI8gf8G7LxpIE5e0GNlkpIni3ddSb14haKJ497XUm4GgkfJkK40vA2ECtH0Sbjo0dWpbSe0wEBappifhUZyA165e3tdNZ9rOG/Y0SgbCItXkJDyqE/AgRgctRE1O5rOdyPud6tQ7wTUfBsIiNakn4XHS5AQ824m836lO59s86HDbyWYgSEPS9GTedLteN+w1vbeiKYfbTi4DQRqB6U7mTeesnk2/fTS9/vof1RzaNl+NLwNBmoNBPEZjUKOvBtU8OKgHFc5mUP0iTRg+/fFOZWkO+gmDXifzpndNj6rzf1R3cffbLzIf/YbPpPMKQZqDfsKg18l83Dr++72LexDNSr36RZo0oQ3yYYcw3n0nvX7WXsdvrgwEaZ4G/Y9ynI2qWamJfq/Srn/v0dOuaxJybTcpDTL4ZmIgSGqsaQc2DGYYLMz/TvvZmtyahFzb93mM6urFQJDUWJNmpabDYGfqr+i3w73XX/9NNA25mdaPsqN8mFekBoKkgRrEXfKjfNzJbCE3qOHBU8a5r8JAkDRQg+gsH6cO90FdrYzqPo/5cNipJM1gUMOD+xnS29YQWK8QJGkGg7pa6bdDvg0GgjSDtueS0OKxEGbrMxCkLk7oo3HQRjiMRR9CkpOT/DDJo0muSvKytmvS5Bq3R0pocsz2ezfsP0Bav0JI8hpgHXAy8M36vxcmeX4p5f+2Wpwm0jiNcNFkmamfYRR/gLQeCMBbgbNLKWfV79+U5Gjgj4BT2ytLkkar7T9GWg2EJDsBLwY+2LXqa8CvTrP9+vrLFcOtTJImT9t9CE8FlgA/6lr+I2Cv0ZcjSZOr7UDoSyllVSllFbCx7VokabFpOxDuB7YBT+9a/nTgntGXI0mTq9VAKKVsBa4CjuxadSTwrdFXJEmTaxxGGZ0BnJPkO8BlwBuAvYFPtlqVJE2Y1gOhlPKFJP8J+AvgGcB1wDGllDtm+Nj+GzduZNWqVaMoUZIWhY0bNwLs32t9Simjq2ZAklwN7Anc2mDzqSGqdkQPlsd1ODyug+cx3W5/4L5SykHTrVyQgdCPqXsX6tFJGhCP63B4XAfPY9pc26OMJEljwkCQJAET0GQkSWrGKwRJEmAgSJJqBoIkCTAQJEm1RR0ITs3ZnySnJSldr3s61qfe5q4kjyRZn+QFXfv4hSTnJNlUv85Jssfof5r2JFmZ5MtJ7qyP4Uld6wdyHJO8MMmGeh93Jnl3kozgRxy5Bsf07Gl+d7/dtc3OST6a5P4kW+r9PbNrm32TnFevvz/JR+p5WybCog2Ejqk5/wo4iOpheRcm2bfVwsbfTVSPEJl6vbBj3TuAtwFvAn4FuBf45yS7dmzzt8DBwNH162DgnOGXPVZ2oXoEy58Aj0yzft7HMcluwD9TzR3yK/X3OoVqBsLFaLZjCnAxO/7uHtO1/kzglcDxwMuA3YCvJFkCUP/3fGDXev3xwKuADw3yBxlrpZRF+QKuAM7qWnYL8IG2axvXF3AacF2PdQHuBv68Y9mTgJ8A/61+/zygAId1bHN4vey5bf98LR3Th4CTBn0cqaaY3Qw8qWObvwDupB5Ovlhf3ce0XnY28JUZPrM7sBV4bceyZwE/B46q3/+X+v2zOrY5AXgU2K3tn3sUr0V5hdAxNefXulZNOzWndrCsbsr4YZJzk0xN8PocqlnsHjumpZRHgEvZfkwPpfrH2vno8suALXjcpwzqOB4K/Ev92SlfpXpS8H7DKHwBODzJvUluTnJWkqd1rHsx8AR2PO7/CtzAjsf0hnr5lK8CO9efX/QWZSDg1JxzdQVwElUTxRqqY/Wt+mm0U8dtpmO6F9WDsx6727H++l487lMGdRz36rGPzu8xSS4CXgf8OlVz3EuAbyTZuV6/F9VkXPd3fa77uHcf06lJvCbimLb++GuNj1LKhZ3v606524DfBb497YekMVBKObfj7bVJrgLuAI4F/qmdqhaexXqF4NScA1BKeQi4HjiA7cdtpmN6D7Bn50iX+uun4XGfMqjjeE+PfXR+j4lVSrkL+Deq312ojskSqtaDTt3HvfuYTrU2TMQxXZSBUJyacyCSPBH4RapO0B9S/aM4smv9y9h+TC+nGg1yaMduDgWW4nGfMqjjeDnwsvqzU44E7gJuH0bhC0mSpwL7UP3uQnU++Bk7HvdnUnXgdx7T53UNRT0S+Gn9+cWv7V7tYb2A11CNKng91f/0dVQddc9uu7ZxfQEfBI6g6vh8KfAVqpEsz67X/ymwCfgt4EDgXKoT0K4d+7gQuJbqBHZo/fV5bf9sIz6Ou1BNyrICeBh4d/31voM6jlSjZu6pP3tgva/NwNva/vlHfUzrdR+sj9N+wCqqk/u/dR3T/1EvW001FP0SqklzltTrl9TH+Rv1+tVUo7Y+2vbPP7Lj3HYBQ/4lOpnqr6WphF/Zdk3j/Oo4MW2t/yH8I/D8jvWhGpp6N9VQvA3AgV37+AXgc/XJaXP99R5t/2wjPo6rqIaIdr/OHuRxpLpH5NJ6H3cD72GRDjmd6ZhSDdv9KlWn+1aqvoOz6Rg+Wu9jZ+CjwAN1qJw3zTb7Uv0h9HC93UeAndv++Uf18vHXkiRgkfYhSJL6ZyBIkgADQZJUMxAkSYCBIEmqGQiSJMBAkCTVDARJEmAgSI3V0yl+peP91JSjI3lqcJK1Sa5N4r9bDYW/WFIDSf4z8AaqR0605a+BPakeRy4NnIEgNbMWuKaU8t22CijV7GifBd7eVg1a3AwETZQkh9XNPNO9PtnjMztTza37tw32f3SSh5J8bKppp6Np6YVJLknycJK7k7y3u/knyYuSfDHJA0keSXJTklM7NjkXeH4SpyTVwDljmibNjew4zwDAO6mmDf1Cj88cAuwB/MtMO07yOuBTwHtLKX85zSZfAv4X8AHgKOBdVJO6n1Z//iXAeuBW4C1sn+Dllzr2sRH4SV2vc0xooAwETZRSygNUjzUGIMl7qU7Ox5VSLunxsUOoHrX8vV77TfIO4P3AH5VSPtVjs7NKKafXX38tyW7A25KcWUr5MdUz/R8ADimlPFxv942u+n+e5Jq6JmmgbDLSxEpyOtWE7L9RSrlohk33BjaXaia+6XwY+O/Aq2YIA4C/73p/LtXkLgcmeTJwGPD5jjDo5b66JmmgDARNpCRnAG8EjimlXDzL5k+kmmSpl+OB64DZ9vOjHu/3oZoQ53FUzUSzeYRqUhhpoAwETZwkH6WaWvWoUsqGBh95gKoPoZdfp5pp68Iku8ywXfcE7lPv7wQepOpP2KdBPU8B7m+wndQXA0ETI5VPAicCR5ZSLmv40RuBnbomX+90PdUUjwcwcyi8uuv971DN831t3Uz0TeCEJLP99f8c4KYmhUv9sFNZk2Qd1ZXBn1DlQ2fH7PdLKZt7fO7S+r8voUeTTinlhiSrqCZu/2qSo0spP+nabE09zPRKqo7s1wOnlVI21evfTjW/8uVJPlR/r2XAilLKm6iK3gNYTtUBLQ2UVwiaCElCdYfvEuBjwOVdr716fbaUcjvwHeA3ZvoepZSbgCOAZ7N9FFGnVwBHAl+muq/hL4H3dXz+SqqO5X+lmgz+AuAUdgyhY6kmkv/iTLVIc5FSSts1SGMvyUlUVxjPaDAKqPuzpwHvAZ5QSvn3edZxIXB/KeXE+exHmo5XCFIznwPuAk5uq4AkK4CXUw1xlQbOQJAaqP+y/z2gr6uDAdsLOKmUcmuLNWgRs8lIkgR4hSBJqhkIkiTAQJAk1QwESRJgIEiSagaCJAkwECRJtf8PNTgwFSl/O3MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import rayleigh\n",
    "s.coherence_func = rayleigh(50).rvs\n",
    "s.domain.create_box_array(Lmax, i, s.coherence_func, r0=0)\n",
    "plt.step(s.domain.rcen, s.domain.B * 1e6, where=\"mid\")\n",
    "plt.ylabel(\"$B_\\perp (\\mu G)$\")\n",
    "plt.xlabel(\"$z$ (kpc)\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
