{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os.path import join\n",
    "\n",
    "import csiborgtools\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "import scienceplots  # noqa\n",
    "from cache_to_disk import cache_to_disk, delete_disk_caches_for_function  # noqa\n",
    "\n",
    "import plt_utils\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Field evaluated at radial shells"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: This is a little dodgy\n",
    "\n",
    "def plot_field_shells(field, MAS, grid, to_save=True):\n",
    "    folder = \"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells\"\n",
    "\n",
    "    # with plt.style.context(\"notebook\"):\n",
    "    if True:\n",
    "        cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "        lw = plt.rcParams['lines.linewidth']\n",
    "        plt.figure()\n",
    "\n",
    "        # CSiBORG2 main\n",
    "        fname = join(folder, f\"csiborg2_main_{field}_{MAS}_{grid}.npz\")\n",
    "        file = numpy.load(fname)\n",
    "        dist, mean = file[\"distances\"], file[\"mean\"]\n",
    "        mean /= dist\n",
    "        mean /= 70\n",
    "        for i in range(len(mean)):\n",
    "            plt.plot(dist, mean[i], c=cols[0], label=\"CSiBORG\" if i == 0 else None)\n",
    "        \n",
    "        # # BORG2\n",
    "        # fname = join(folder, f\"borg2_{field}_{MAS}_{grid}.npz\")\n",
    "        # file = numpy.load(fname)\n",
    "        # dist, mean = file[\"distances\"], file[\"mean\"]\n",
    "        # for i in range(len(mean)):\n",
    "        #     plt.plot(dist, mean[i], c=cols[2], label=\"BORG\" if i == 0 else None)\n",
    "\n",
    "        # # CSiBORG2 random\n",
    "        # fname = join(folder, f\"csiborg2_random_{field}_{MAS}_{grid}.npz\")\n",
    "        # file = numpy.load(fname)\n",
    "        # dist, mean = file[\"distances\"], file[\"mean\"]\n",
    "\n",
    "        # mu = numpy.mean(mean, axis=0)\n",
    "        # std = numpy.std(mean, axis=0)\n",
    "\n",
    "        # plt.fill_between(dist, mu - std, mu + std, alpha=1/3, color=cols[1])\n",
    "\n",
    "        # for i in range(len(mean)):\n",
    "            # plt.plot(dist, mean[i], c=cols[1], label=\"Random\" if i == 0 else None, zorder=0, lw=lw/2)\n",
    "\n",
    "        # Plot settings\n",
    "        plt.legend(loc=\"lower right\")\n",
    "        plt.xlabel(r\"$r ~ [\\mathrm{Mpc} / h]$\")\n",
    "\n",
    "        if field == \"radvel\":\n",
    "            plt.ylabel(r\"$\\langle v_r \\rangle ~ [\\mathrm{km} / s]$\")\n",
    "            plt.axhline(0, c=\"k\", ls=\"--\",)\n",
    "            plt.ylim(-0.1, 0.1)\n",
    "            plt.xscale(\"log\")\n",
    "        elif field == \"overdensity\":\n",
    "            plt.ylim(-0.5, 0.5)\n",
    "            plt.axhline(0, c=\"k\", ls=\"--\",)\n",
    "            # plt.xlim(0, 200)\n",
    "            plt.ylabel(r\"$\\langle \\delta_r \\rangle$\")\n",
    "            # plt.axhline(-0.1, c=\"k\", ls=\"--\")\n",
    "        elif field == \"density\":\n",
    "            plt.axhline(277.5 * 0.304, c=\"k\", ls=\"--\",)\n",
    "            plt.ylim(50, 100)\n",
    "\n",
    "        plt.xlim(0, dist.max())\n",
    "        # plt.xlim(0, 100)\n",
    "\n",
    "        if to_save:\n",
    "            fout = join(plt_utils.fout, f\"field_shells_{field}_{MAS}_{grid}.png\")\n",
    "            print(f\"Saving to `{fout}`.\")\n",
    "            plt.savefig(fout, dpi=plt_utils.dpi, bbox_inches=\"tight\")\n",
    "\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_16251/4241837727.py:58: UserWarning: Attempt to set non-positive xlim on a log-scaled axis will be ignored.\n",
      "  plt.xlim(0, dist.max())\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAG+CAYAAAC+mCFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABnc0lEQVR4nO3deVyU5fo/8M/MwAw7yDqgICgo7igq4r5guJRZeo6Z5ZJHW9RS6pR6SltO6clKMzVP38zsdzLNtEWPx1IUcSFREFcE95VVhGEdhpn5/UFMTcwMDMwKn/frNa/iea7nmWvUmbl47vu5boFarVaDiIiIiExCaO0EiIiIiFoSFldEREREJsTiioiIiMiEWFwRERERmRCLKyIiIiITYnFFREREZEIsroiIiIhMiMUVERERkQmxuCIiIiIyIRZXRERERCZk18XV+vXrERoaCicnJ8TExCA1NVVv7IULFzBp0iSEhoZCIBBgzZo1TTpnVVUV5s2bBx8fH7i5uWHSpEnIy8sz5csiIiIiO2a3xdX27duRkJCA5cuXIz09Hb169UJ8fDzy8/N1xldUVKBDhw5YuXIlpFJpk8+5aNEi7N69Gzt27MDhw4dx7949PP7442Z5jURERGR/BPa6cHNMTAz69euHdevWAQBUKhWCg4OxYMECLF682OCxoaGhWLhwIRYuXGjUOUtKSuDn54etW7di8uTJAIBLly6hS5cuSElJwYABA0z/QomIiMiuOFg7gaaorq5GWloalixZotkmFAoRFxeHlJQUs50zLS0NCoUCcXFxmpjIyEiEhIToLa7kcjnkcrnmZ5VKhaKiIvj4+EAgEDQpVyIiIrIstVqN0tJSBAUFQSg0PPBnl8VVYWEhlEolAgICtLYHBATg0qVLZjtnbm4uxGIxvLy86sXk5ubqPO+KFSvw1ltvNSknIiIisi23b99Gu3btDMbYZXFlT5YsWYKEhATNzyUlJQgJCcHt27fh4eFhxcyIiIiosWQyGYKDg+Hu7t5grF0WV76+vhCJRPXu0svLy9M7Wd0U55RKpaiurkZxcbHW1StDzyuRSCCRSOpt9/DwYHFFRERkZxozpccu7xYUi8WIjo5GYmKiZptKpUJiYiJiY2PNds7o6Gg4OjpqxWRlZeHWrVtNfl4iIiJqWezyyhUAJCQkYMaMGejbty/69++PNWvWoLy8HLNmzQIATJ8+HW3btsWKFSsA1E5Yv3jxoub/7969i4yMDLi5uSE8PLxR5/T09MTs2bORkJAAb29veHh4YMGCBYiNjeWdgkRERATAjourKVOmoKCgAMuWLUNubi6ioqKwb98+zYT0W7duac3mv3fvHnr37q35+YMPPsAHH3yAYcOGISkpqVHnBIDVq1dDKBRi0qRJkMvliI+Px4YNGyzzoomIiMjm2W2fK3slk8ng6emJkpISzrkiIiKyE8Z8f9vlnCsiIiIiW8XiioiIiMiEWFwRERERmRCLKyIiIiITYnFFREREZEIsroiIiIhMiMUVERERkQmxuCIiIiIyIRZXRERERCbE4oqIiIjIhFhcEREREZkQiysiIiIiE2JxRURERGRCLK6IiIiITIjFFREREZEJsbgiIiIiMiEWV0REREQmxOKKiIiIyIRYXBERERGZEIsrIiIiIhNicUVERERkQiyuiIiIiEyIxRURERGRCbG4IiIiIjIhFldEREREJsTiioiIiMiEWFwRERERmRCLKyIiIiITYnFFREREZEIsroiIiIhMiMUVERERkQmxuCIiIiIyIRZXRERERCbE4oqIiIjIhFhcEREREZmQXRdX69evR2hoKJycnBATE4PU1FSD8Tt27EBkZCScnJzQo0cP7N27V2u/QCDQ+Vi1apUmJjQ0tN7+lStXmuX1ERERkf2x2+Jq+/btSEhIwPLly5Geno5evXohPj4e+fn5OuOPHz+OqVOnYvbs2Th9+jQmTpyIiRMn4vz585qYnJwcrccXX3wBgUCASZMmaZ3r7bff1opbsGCBWV8rERER2Q+BWq1WWzuJpoiJiUG/fv2wbt06AIBKpUJwcDAWLFiAxYsX14ufMmUKysvLsWfPHs22AQMGICoqChs3btT5HBMnTkRpaSkSExM120JDQ7Fw4UIsXLiwSXnLZDJ4enqipKQEHh4eTToHERERWZYx3992eeWquroaaWlpiIuL02wTCoWIi4tDSkqKzmNSUlK04gEgPj5eb3xeXh7++9//Yvbs2fX2rVy5Ej4+PujduzdWrVqFmpoavbnK5XLIZDKtBxEREbVcDtZOoCkKCwuhVCoREBCgtT0gIACXLl3SeUxubq7O+NzcXJ3xW7Zsgbu7Ox5//HGt7S+++CL69OkDb29vHD9+HEuWLEFOTg4++ugjnedZsWIF3nrrrca+NCIiIrJzdllcWcIXX3yBadOmwcnJSWt7QkKC5v979uwJsViMZ599FitWrIBEIql3niVLlmgdI5PJEBwcbL7EiYiIyKrssrjy9fWFSCRCXl6e1va8vDxIpVKdx0il0kbHHzlyBFlZWdi+fXuDucTExKCmpgY3btxA586d6+2XSCQ6iy4iIiJqmexyzpVYLEZ0dLTWRHOVSoXExETExsbqPCY2NlYrHgD279+vM37Tpk2Ijo5Gr169GswlIyMDQqEQ/v7+Rr4KIiIiaons8soVUDs8N2PGDPTt2xf9+/fHmjVrUF5ejlmzZgEApk+fjrZt22LFihUAgJdeegnDhg3Dhx9+iPHjx2Pbtm04deoUPvvsM63zymQy7NixAx9++GG950xJScGJEycwYsQIuLu7IyUlBYsWLcJTTz2FNm3amP9FExERkc2z2+JqypQpKCgowLJly5Cbm4uoqCjs27dPM2n91q1bEAp/vzA3cOBAbN26Fa+//jqWLl2KiIgI/PDDD+jevbvWebdt2wa1Wo2pU6fWe06JRIJt27bhzTffhFwuR1hYGBYtWqQ1p4qIiIhaN7vtc2Wv2OeKiIjI/rT4PldEREREtorFFREREZEJsbgiIiIiMiEWV0REREQmxOKKiIiIyIRYXBERERGZEIsrIiIiIhNicUVERERkQiyuiIiIiEyIxRURERGRCbG4IiIiIjIhFldEREREJsTiioiIiMiEWFwRERERmRCLKyIiIiITYnFFREREZEIsroiIiIhMiMUVERERkQmxuCIiIiIyIRZXRERERCbE4oqIiIjIhFhcEREREZkQiysiIiIiE2JxRURERGRCLK6IiIiITIjFFREREZEJsbgiIiIiMiEWV0REREQmxOKKiIiIyIRYXBERERGZEIsrIiIiIhNicUVERERkQiyuiIiIiEyIxRURERGRCbG4IiIiIjIhuy6u1q9fj9DQUDg5OSEmJgapqakG43fs2IHIyEg4OTmhR48e2Lt3r9b+mTNnQiAQaD3GjBmjFVNUVIRp06bBw8MDXl5emD17NsrKykz+2oiIiMg+2W1xtX37diQkJGD58uVIT09Hr169EB8fj/z8fJ3xx48fx9SpUzF79mycPn0aEydOxMSJE3H+/HmtuDFjxiAnJ0fz+Oabb7T2T5s2DRcuXMD+/fuxZ88eJCcnY+7cuWZ7nURERGRfBGq1Wm3tJJoiJiYG/fr1w7p16wAAKpUKwcHBWLBgARYvXlwvfsqUKSgvL8eePXs02wYMGICoqChs3LgRQO2Vq+LiYvzwww86nzMzMxNdu3bFyZMn0bdvXwDAvn37MG7cONy5cwdBQUEN5i2TyeDp6YmSkhJ4eHgY+7KJiIjICoz5/rbLK1fV1dVIS0tDXFycZptQKERcXBxSUlJ0HpOSkqIVDwDx8fH14pOSkuDv74/OnTvj+eefx/3797XO4eXlpSmsACAuLg5CoRAnTpzQ+bxyuRwymUzrQURERC2XXRZXhYWFUCqVCAgI0NoeEBCA3Nxcncfk5uY2GD9mzBh89dVXSExMxL/+9S8cPnwYY8eOhVKp1JzD399f6xwODg7w9vbW+7wrVqyAp6en5hEcHGz06yUiIiL74WDtBGzJE088ofn/Hj16oGfPnujYsSOSkpIwatSoJp1zyZIlSEhI0Pwsk8lYYBEREbVgdnnlytfXFyKRCHl5eVrb8/LyIJVKdR4jlUqNigeADh06wNfXF1euXNGc488T5mtqalBUVKT3PBKJBB4eHloPIiIiarnssrgSi8WIjo5GYmKiZptKpUJiYiJiY2N1HhMbG6sVDwD79+/XGw8Ad+7cwf379xEYGKg5R3FxMdLS0jQxBw8ehEqlQkxMTHNeEhEREbUQdllcAUBCQgL+7//+D1u2bEFmZiaef/55lJeXY9asWQCA6dOnY8mSJZr4l156Cfv27cOHH36IS5cu4c0338SpU6cwf/58AEBZWRn+/ve/49dff8WNGzeQmJiIRx99FOHh4YiPjwcAdOnSBWPGjMGcOXOQmpqKY8eOYf78+XjiiScadacgERERtXx2O+dqypQpKCgowLJly5Cbm4uoqCjs27dPM2n91q1bEAp/rx0HDhyIrVu34vXXX8fSpUsRERGBH374Ad27dwcAiEQinD17Flu2bEFxcTGCgoLw0EMP4Z133oFEItGc5+uvv8b8+fMxatQoCIVCTJo0CWvXrrXsiyciIiKbZbd9ruwV+1wRERHZnxbf54qIiIjIVrG4IiIiIjIhFldEREREJsTiioiIiMiEWFwRERERmRCLKyIiIiITYnFFREREZEIsroiIiIhMiMUVERERkQmxuCIiIiIyIRZXRERERCbE4oqIiIjIhFhcEREREZkQiysiIiIiE2JxRURERGRCLK6IiIiITIjFFREREZEJsbgiIiIiMiEWV0REREQmxOKKiIiIyIRYXBERERGZEIsrIiIiIhNicUVERERkQiyuiIiIiEyIxRURERGRCbG4IiIiIjIhFldEREREJsTiioiIiMiEWFwRERERmRCLKyIiIiITYnFFREREZEIsroiIiIhMyMHaCRARUfOk3SzC05tS4eQoQkyYNwaG+yI6pA06S90hEgoAAEqVGgA0PxOR+bC4IiKyY1UKJWZ/eQoV1UpUVCvxv/O5+N/5XACAk6MQ3YM8UVFdg+y8Mni7OmLHcwPR3sfVylkTtWwcFiQismOvfncWxZUKnfuqFCqcuvkAF3NKUaNSI7+0GiM/SMKpG0UWzpKodWFxRURkp45cLsBPZ+4ZdYxSDUzemILvTt02U1ZEZNfF1fr16xEaGgonJyfExMQgNTXVYPyOHTsQGRkJJycn9OjRA3v37tXsUygUeO2119CjRw+4uroiKCgI06dPx7172h9coaGhEAgEWo+VK1ea5fUREelTJq/BvK/Tm3z8K9+dxeRPj0GlUpkwKyIC7Li42r59OxISErB8+XKkp6ejV69eiI+PR35+vs7448ePY+rUqZg9ezZOnz6NiRMnYuLEiTh//jwAoKKiAunp6XjjjTeQnp6OXbt2ISsrCxMmTKh3rrfffhs5OTmax4IFC8z6WomI/uzvO85AVlXTrHOculmMDkv/h3/uuYASPUOLRGQ8gVqtVls7iaaIiYlBv379sG7dOgCASqVCcHAwFixYgMWLF9eLnzJlCsrLy7Fnzx7NtgEDBiAqKgobN27U+RwnT55E//79cfPmTYSEhACovXK1cOFCLFy4sFF5yuVyyOVyzc8ymQzBwcEoKSmBh4dHY18uEZFG4sU8zP7qVKPjJQ5CSD0kuFVUCUMf+C+ODMdLcZ14RyGRDjKZDJ6eno36/rbLK1fV1dVIS0tDXFycZptQKERcXBxSUlJ0HpOSkqIVDwDx8fF64wGgpKQEAoEAXl5eWttXrlwJHx8f9O7dG6tWrUJNjf7fHlesWAFPT0/NIzg4uBGvkIhIt+KKary0PcOoY+Q1Ktz8rbAyVDetPXgFj68/ijxZVbNyJGrt7LK4KiwshFKpREBAgNb2gIAA5Obm6jwmNzfXqPiqqiq89tprmDp1qlaF+uKLL2Lbtm04dOgQnn32Wbz33nt49dVX9ea6ZMkSlJSUaB63b3MSKRE13cvfnkGZvOnDgaoGxirO3JVhxKpDSMzMa/JzELV27HOlg0KhwF//+leo1Wp8+umnWvsSEhI0/9+zZ0+IxWI8++yzWLFiBSQSSb1zSSQSnduJiIy158w9JF7SPa8UAGYObI/o9t7wd5fA38MJ/u4SuEocUKNUoaRSgY/2Z2P7yduoaaDCqlCoMHvLKcwc2B5LxnWBxEFk6pdC1KLZ5ZUrX19fiEQi5OVp/2aVl5cHqVSq8xipVNqo+LrC6ubNm9i/f3+D46oxMTGoqanBjRs3jH8hRESNVFAqxys7zujdL3V3xJsTuuORXkGI6eCDMF9XuEpqf392EAnh4ybBu4/1wKV3xmD+iI6Nmlf15fGbmPDJMVwtKDPZ6yBqDeyyuBKLxYiOjkZiYqJmm0qlQmJiImJjY3UeExsbqxUPAPv379eKryusLl++jAMHDsDHx6fBXDIyMiAUCuHv79/EV0NEZJharcai7RmoqtHfNmH3i0MbdS4HkRCvxEfi+OKReKRnYIPxWXmlGLMmGbvS7zQ6X6LWzm6HBRMSEjBjxgz07dsX/fv3x5o1a1BeXo5Zs2YBAKZPn462bdtixYoVAICXXnoJw4YNw4cffojx48dj27ZtOHXqFD777DMAtYXV5MmTkZ6ejj179kCpVGrmY3l7e0MsFiMlJQUnTpzAiBEj4O7ujpSUFCxatAhPPfUU2rRpY50/CCJq8Xam3cHRK4V698+IDYGfu5NR5wzwcMInT/bB1P6FeOPH87haUK43VqFUI+HbM/jlQg4+ntqHw4REDbDbVgwAsG7dOqxatQq5ubmIiorC2rVrERMTAwAYPnw4QkND8eWXX2rid+zYgddffx03btxAREQE3n//fYwbNw4AcOPGDYSFhel8nkOHDmH48OFIT0/HCy+8gEuXLkEulyMsLAxPP/00EhISGj2vyphbOYmICkrlGLQyEdVK3R/VjkIg+91xEAia3j6hukaFz49ew0e/ZMHAxTEAgIujEKv+0gvjewY1+fmI7JEx3992XVzZIxZXRGSMNQeysebAZb37f3whFr1CvE3yXLeLKvD3HWfw6/WG1x70cxNj6fgumBjVtlmFHZG9aPF9roiIWotvUm/p3dct0N1khRUABHu74Ju5A/DeY90bjC0oq8ai7WfQffnP+PCXLCiVXEaHqA6LKyIiG3XrfgXyZHK9+3+YN8jkzykQCPBkTHscSBgKF8eGr0iVVyvxycEr6LJsH9786QLkCqXJcyKyNyyuiIhs1NbUm3r3LR0XCUczTiwP93fHwVdGIsBd3Kj4aqUaXx6/gcg39mHRttMoreJahdR6cc6VhXHOFRE1VvQ7v+B+ef0ixdVRiAvvjLVIDvmyKiz/6QL+d173ahaGCAW1D4FAAKVKrekO38bFEY/2CsKk6GB0CXSHg4i/55Pt44R2G8biioga42pBGUZ9eFjnvkOvDEOYr5tF88nOK8UniZex+2yOyc/t7eKAoZ398cKwjugk5eci2SazFVc//fST0cmMHj0azs7ORh/XUrG4IqLGeGfPRWw6er3e9tiOPvhmzgArZFTrRmE5lv14HsmX9ffdag4hAD8PMboGemJCryA81qedWZ6HyFhmK66EQuMu3QoEAly+fBkdOnQw6riWjMUVETVErVaj99u/oLiy/gLNl98dC0cbGEa7eK8Ez/8nHTeLKsz+XF0DPfDUgBA80jMI7s6OZn8+Il2M+f42ukN7bm5uo5d6cXd3N/b0REStXlZeqc7C6p8Tu9lEYQUAXYM8cfjVEcgtqcT2k7fxXdod3H5QaZbnupgjw9Lvz2Pp9+fRxsURb03ojkd6BbK/Ftkso4qrGTNmGDXE99RTT/HqDBGRkbae0N3b6qkBoZZNpBGkns54Ka4TXorrhOKKany0PxvbT96GvEYFAYAgL2d0DfRAgKcE5fIaXC8sw+2iKjyoqNZMcDfGgwoFXtx2Gh8nZuNfk3qib6jhPl+yKgWOZBeijasjYsJ8GrVgNVFzcUK7hZl7WLBKoURBqRzB3i4mPzcRmZ9arUbPN39GqVy7X1R8V3/8e3o/K2VlnMpqJS7myBAR4AYPJ/3DeNU1KuSUVOCTg1fw69X7yJdVodrIXqSDw30wOToYHs4OcJM4wk3iAFexCFcKyrDr9F38ciEXit+WDvJ1E+PhnkGYEBWE3sFevPJFRrHI3YLPPPMMhg4dipkzZwIAbt68iYsXL2LgwIHw9PRsyilbBXMWVwcv5eH5/6SjZztP7HhuoEnPTUSWce5OCR5Zd7Te9lP/iIOve+PWMLV3VwtK8c2JW9h09AbM9dt/uzbOeGlUBCb2bmszQ61k2yyy/M3evXsRGRkJACguLkZ0dDQmTpyIrl27Iisrq6mnpWaI8HeHvEaF9FvFkLGBH5Fd0tc4tLUUVgDQ0c8drz/cDWfffAh/6W2eBaLvPKjE3787iz5v78enh6+iip3lyYSaXFyVlJSgbdu2AICdO3dCKpVCJpNhypQpWLJkickSpMYL9nZBB19XKFVqHL9y39rpEJGR1Go1fjx9t972h3sEWCEb63N3csSqKb3x88KhCGnjZJbnKJXX4F//u4Ruy3/Gwm3pKCzTv9wQUWM1ubgKDg7G9eu1PVh27NiBmTNnQiKR4LnnnsOxY8dMliAZZ2gnPwBA8uUCK2dCRMZKv1WMCkX9SUdvPNLwQsotWWepOw6/OhIrH++BQE/zXMFTqtT4ISMHff95AANXHMA3J/QvPaSPWq2GpaYxq5pyNwBZjNGtGOrMnDkTL774Ih555BEkJiZi3bp1AACVSoWysjKTJUjGGRLhiy+P30BydgHUajUnbBLZEX1f6AEe5rlqY08EAgGe6B+CJ/qH4MytYiTsyMDVgnKtGD9XR4zrGQQXsQhlciVkldUoqqjGveJKXCuoaPT8rXslciz5/jz+8f15DO3sh39N6oEAj9/vlH9QLkdSVgESL+Xj4r0SVFarUFWjrF1PUQ0EeDrBz90JwW2c8WhUEHzcxLicV4aB4b5o18b4m40u3pPh6xM3sefMPZTKa7TusnQUAWO6BeL9yb3gLDbfWpNknCYXV0uWLIFarcYvv/yClStXIjw8HABw8uRJhISEmCxBMs6ADj5wFAlw50ElbtyvQJivq7VTIqJGUKrU2KNjaZmHe0itkI1t6xXihQMJw/Bjxl28vTsTXi6OWPF4D8R08NF7jEqlRlJWAT45mI3Tt0sa9TwqAElZBYh576BR+d0trsLd4ipk3C6ut1yQr5sYk6Lb4bHebRFpYKkftVqNHzPuYcXeTOSV6h+qVCiB3WdzsPtsDkJ9XPDRlCj0CWljVL5keiZvxbBq1SpUVVXhjTfeMOVpWwxLdGif+tmvSLl2H29N6IYZA0PN8hxEZFq/XruPJz77td72Y4tHoq0XlxDTpylX6MvkNfjh9F1sOX4dl/PLGz7ATLoGumPmoDDEdvBBWy9nCIUCXM0vxedHr2P3mXsokzdtkr2XswNmD+6Aib3bsi2PCZm1FcOyZcvw6KOPIjo6ullJtlaWKK4+TbqKf+27hFGR/tg00z764hC1di99k44fz9S/cnVj5XgrZNN6lFYpsDPtDr5Lu42LOaVNamxqKmKRANVK0yYQ4CHB3wZ3QN/QNoiUenDosBnMuvzNnTt3MHbsWIjFYjzyyCOYMGECRo0aBbFY3OSEybSGdvLFv/YBKdfuo7pGBbEDe7gQ2bIapQr7zufW2z6uld4laEnuTo6YOSgMMweFAQDyZVVYvOssUq7eR6WOmwv+rJ2nBL3be0EAIQpK5SiqqEZZlQJl1TUoqTTuylNjCysBgAh/N5RV16BAVgVDaebJ5Hh3b6bmZ6EAaO/jgvhuUoztHohewV5G5UiN06RhQZVKhWPHjmH37t348ccfkZOTg9GjR+PRRx/Fww8/DG9vw8sRtGaWuHKlUqnR/71EFJbJsXVODAZ29DXL8xCRaRy5XICnN6XW25709+EI9eG8SWuprlHhYGYedqbfhqODCE4OIjwor4azWITXxnZGex+3Bs/xoEKOr47fxOZj13WuF2mMSKk7FoyMwJBOvprO9yqVGkeuFODDX7Jw9o7M6HM+0isIr8Z35vBhI1ikQ/sfZWZmagqtU6dOISYmBhMmTMDUqVM1vbColiWKKwBYtD0D35++i+eHd8RrYyLN9jxE1Hwv/OcU9p7Pq7edQ4Ity4MyOV7adhrJRvYhFAFY9FAnPDuso8Fu8mXyGqzadwlbT9w0eDXrzxyFAswYGIp5I8LRxtXyo1BqtRrZeWU4cf0+ztwuQU5JJRRKFfzdnRDu74YQbxc4iASQK1RQKFWoVipRVK5ArqwKihoV3CQO8HB2hLxGCaFAgEBPZ7hJHOAgEkAgAFQqwMlRhK5BHgj1cWnyXfQWL67+qKCgQFNoDRkyBK+88oopT2/3LFVcfX/6DhZtP4NuQR7474tDzPY8RNQ81TUqdFv+P/y5QfjoLn74vxn9rZMUmVWVQonvT9/FxsNXcfN+Rb39IgEgcRBiULgvYjv6Iq5LAEJ8Gn9lSa1W4+iVQvxzz0Vk5TW+NVIbF0esnhKF4Z39G31Mc9y6X4Gtqbfww+m7yJVVWeQ5pR5OiOngjdgOPhjZxR/+7o1vc2Kx4qqqqgpnz55Ffn4+VCrtMnnChAlNPW2LZqniqrBMjr7/PAAAOPmPOPi1oqUziOzJwUt5eObLU/W27180FBEB7lbIiCxFrVbjxv0KXCsow/XCchRXKDCyi79JF5UurqjG6gPZ+PH0PRRXNm5ZtH6hbfDC8HBEh7YxuPB2U5XLa7D24GV8cfS6ZlFtc2vj4ohyuRLVyt9rFZFQgPhuAVgwMgJdAhv+PrZIcbVv3z5Mnz4dhYWF9U8qEECp5DpNuliquAKA8WuP4MI9GVZP6YXHercz63MRUdPM3nISiZn59bZzSJBMTVGjxIf7s7Hx8LVGxYsEwJBOfng0KggP9wwyyQLXaTeL8NK2DNx5UNnsczVFoIcEAzr64Gp+Gc7e/X2O2sSoICyM64RQA70hLVJcRURE4KGHHsKyZcsQEMA7WhrLksXVv/ZdwqdJV/FY77ZYPSXKrM9FRMarUijRbdk+/PmX92GdfbFlVox1kqIWL7+0CjM2/YrM3Mb3+GrXxgkvjuqEyX3aQSg0/qraqRtFWP7TBVy4Z/yke3MQCQUY2MEbSjVw/GrtHDihABjR2R8Te7fF6K4BcHLUblthkeLKw8MDp0+fRseOHZtyeKtlyeLq+NVCPPl/J+DrJkbq0rgmvSGIyHz2nc/Fc/9Jq7d9z4LB6N7W0woZUWuhVKnx9k/nseXXW0YdFxPWBu9P7oX2jbyLNe1mEV7ZcQbXC+vPLbMlf+4x1sbFEdNjQzFnaAe4SWq7Vpm1z1WdyZMnIykpicWVDevb3hsuYhEKy6qRmStDtyB+WBPZkm9Sda8lyMKKzE0kFOCtiT0wPNIfr353FgVl1Y067sT1Bxj5QRL+NqQDXoqLgItYdxmRnSvD05tSDS7d0xwh3s5YOq4LooLbwN9dgkpFDd797yXsTLuFpjS2/3OPsQcVCnyceBmfJV9FfDcpBoX7wsuh8a00mnzlqqKiAn/5y1/g5+eHHj16wNFRe9Lbiy++2JTTtniWvHIFALO/PInES/l4bUwknh/OQpjIVlRW1w4J/vmO+QFh3tj2bKxVcqLWa9+5HLz730zcLm78XCgBapfwWTAyAoMifPGgXIEjlwvw/r6LKKkyoheEETycHPDeYz0wvmegzkn/arUaaTcf4PvTd3G7qAK3H1Tg7oMqrYnsTaWSV+D2mr+a98rVN998g19++QVOTk5ISkrSepECgYDFlY0Y2skPiZfykZxdwOKKyIYcysqvV1gBwN/jO1s8F6IxPQIxpkcgAODAxVws2p6B0gYuAakBXMgpxXNfp5s9vwh/Nyx7uCsGR/gavJNSIBCgb6g3+oZqNzNPu1mE1b9k4+hV43qMNVWTi6t//OMfeOutt7B48WIIhVxexVYN7eQHADh1swjl8hq4Spr8V05EJvTtyds6t/dp38bCmRBpi+sqxbm3xuBGYTmW/3Qeh7PrdwWwBHcnERbFdcL02FA4NPNOxej23vjPnAG4W1yJjw9kY1f6XdSYcSHJJmdbXV2NKVOmsLCycaE+Lgj2doZCqcav1yxTsRORYVUKJZIvF9TbHtXO02T9jYiaK9TXFVueicGBhKHoKm14qR9TaevphI1P9cG5N8fgmcEdml1YaZ3byxnvT+6F7H+OwdonohDuZ57lpZqc8YwZM7B9+3ZT5kJmIBAIMCSi9upVcnb9D3MisrykrHzo+qU54aFOlk+GqAHh/u7Yu3AYts0dAE8n041+CACIHQRwk4jQJdAdfxsciqOvjcCxJaMwpnugyZ5HF6FQiAlRbXHg5eG4/O5YvPlIV3i5mK5hapP/lJRKJd5//338/PPP6NmzZ70J7R999FGzkyPTGBrhh60nbuHIZetc2iUibT+cvqtz++BwPwtnQtR4Azr44Myb8UjOyseSXedwt8T4JWukHhLMGBiKmQPD4CwWNXyABTiKhJg5KAwzB4Xh6OUCvPB1OmRVzVtku8nF1blz59C7d28AwPnz55uVBJnXwHAfiIQCXCssx+2iCq5+TmRFVQolDmTWX6Q5MsCNvejILgzt7I9jS0bhyOUCLNl1Fnce6C+yfN3E6BTgjqGd/DA43Nfm24wMjvDDqddHY/Ox6/g06Wqjlwz6syYXV4cOHWrqoWRhHk6O6BPihZM3HiD5cgGmxbS3dkpErVZydgFqdNwm+GJchOWTIWqGIRF+OJAwHBuSrmLTkWsor669uzDE2wWP92mLpwe0h4+b/a1rK3YQ4tlhHTF3aAfceVCJm/fLce5uCQ6cuQHdt6HU1+Q5V998843efX//+9+belqjrF+/HqGhoXByckJMTAxSU1MNxu/YsQORkZFwcnJCjx49sHfvXq39arUay5YtQ2BgIJydnREXF4fLly9rxRQVFWHatGnw8PCAl5cXZs+ejbKyxq86bi1DOe+KyCbsPnNP5/aHukotnAlR8zk5ipAwuhPOvRmP44tHYv+ioTj89+FYGNfJLgurPxIIBAj2dsHgCD88Pzwcm2f1b/SxTb5y9fzzz8PLywtjx47V2r5o0SJs27YNq1atauqpG2X79u1ISEjAxo0bERMTgzVr1iA+Ph5ZWVnw9/evF3/8+HFMnToVK1aswMMPP4ytW7di4sSJSE9PR/fu3QEA77//PtauXYstW7YgLCwMb7zxBuLj43Hx4kU4OTkBAKZNm4acnBzs378fCoUCs2bNwty5c7F161aj8i8vL4dIVH+8WSQSaZ6rLk4foVAIZ2fnRsX2C6m90+P4lftQKFVQyKugr3+sQCCAi8vvQ4cVFRWNjq2srIRKpb9Zm6ura5Niq6qqDC4Gbkysi4uL5o4suVyOmhr9Y+vGxDo7O2vunq2uroZCof9ysjGxTk5Omn8rxsQqFApUV+vvuiyRSODg4GB0bE1NDeRy/V2XxWKxZg6mMbFKpRJVVfqHFxwdHSEWi42OValUqKzU3xjRmFgHBwdIJLVfGGq1GhUV+pf00BVbpVDiv+k3oPrTbPZQXzetO6IMvZfN9Rnx51hj3vf8jOBnhEQiQZBX7b+flvwZ0WjqJtqzZ4/a09NTfeTIEc22+fPnq4OCgtSZmZlNPW2j9e/fXz1v3jzNz0qlUh0UFKResWKFzvi//vWv6vHjx2tti4mJUT/77LNqtVqtVqlUaqlUql61apVmf3FxsVoikai/+eYbtVqtVl+8eFENQH3y5ElNzP/+9z+1QCBQ3717V+fzVlVVqUtKSjSP27dvq1Hbe03nY9y4cVrHu7i46I0dNmyYVqyvr6/e2L59+6p7vfWzuv1re9Qnr99Xt2/fXm9s165dtc7btWtXvbHt27fXiu3bt6/eWF9fX63YYcOG6Y11cXHRih03bpzBP7c/mjx5ssHYsrIyTeyMGTMMxubn52tiX3jhBYOx169f18S+8sorBmPPnz+viV2+fLnB2NTUVE3s+++/bzD20KFDmth169YZjN2zZ48mdvPmzQZjv/32W03st99+azB28+bNmtg9e/YYjF23bp0m9tChQwZj33//fU1samqqwdjly5drYs+fP28w9pVXXtHEXr9+3WDsCy+8oInNz883GDtjxgxNbFlZmcHYmJHa73tDseb8jPgjfkbU4mdELX5GqNUlJSVqAOqSkhJ1Q5o8LDh+/Hhs2LABEyZMQFpaGl544QXs2rULhw4dQmRkZFNP2yjV1dVIS0tDXFycZptQKERcXBxSUlJ0HpOSkqIVDwDx8fGa+OvXryM3N1crxtPTEzExMZqYlJQUeHl5oW/fvpqYuLg4CIVCnDhxQufzrlixAp6enppHcHBw0160CQwO9wXAoUEiW1P3Gz8RtQxNXluwzoYNG5CQkAA/Pz8cOnQI4eHhpspNr3v37qFt27Y4fvw4YmN/X4Pr1VdfxeHDh3UWOmKxGFu2bMHUqVO1cn/rrbeQl5eH48ePY9CgQbh37x4CA3/vr/HXv/4VAoEA27dvx3vvvYctW7YgKytL69z+/v5466238Pzzz9d7XrlcrnXJUyaTITg4GPfu3dO5NpE5L/nvvlCIV787i17BXvhmVm9e8gcv+XNY0HLDgg9kpej/zwP1FogN83XBzwkjGv2+57Cg7lh+RvAzwtyfEcasDWzUnKuEhASd2/38/NCnTx9s2LBBs419rmpJJBLNB+wfubq6ar3Z9WlMTGNj6ya1n71TDLnaAW1cGzeG/McPxob88cPZlLF//DIxZay+v5/mxorF4kaP0Zsr1tHRsV7/OVPEOjg4aD5ETRkrEoka/e/dmFihUGiWWIFAYFRsRk4lakQSCP801fLVR6Lq/Zs15fu+qbHGvO/5GWF8LD8jjI+1hc+IxjKquDp9+rTO7eHh4ZDJZJr95l6+wdfXFyKRCHl52r1i8vLyIJXqvuNGKpUajK/7b15entaVq7y8PERFRWli8vPztc5RU1ODoqIivc9rS6SeTugc4I6svFIcvVKIR3oFWTslolZj95kcndvH9zBvJ2oisjyjiitb6W0lFosRHR2NxMRETJw4EUDtZb3ExETMnz9f5zGxsbFITEzEwoULNdv279+vGVYMCwuDVCpFYmKippiSyWQ4ceKEZrgvNjYWxcXFSEtLQ3R0NADg4MGDUKlUiImJMc+LNbEhEb7IyitFcnYBiysiC6muUWHvufrFVfcgw0MLRGSf7HbV5YSEBPzf//0ftmzZgszMTDz//PMoLy/HrFmzAADTp0/HkiVLNPEvvfQS9u3bhw8//BCXLl3Cm2++iVOnTmmKMYFAgIULF+Kf//wnfvrpJ5w7dw7Tp09HUFCQpoDr0qULxowZgzlz5iA1NRXHjh3D/Pnz8cQTTyAoyD4KlaGdaocGj1wu1DtHgohM69jVQsh1dA5dPNa8N/8QkXUYdeXq7Nmz6N69u2aSXUMuXLiAzp07N3o81RhTpkxBQUEBli1bhtzcXERFRWHfvn0ICAgAANy6dUsrz4EDB2Lr1q14/fXXsXTpUkREROCHH37Q9LgCaifEl5eXY+7cuSguLsbgwYOxb98+rfH5r7/+GvPnz8eoUaMgFAoxadIkrF271uSvz1z6h3lD4iBErqwKl/PL0CnA3dopEbV4e/Q0Dh302x28RNSyGHW3oEgkQm5uLvz8Gre4qIeHBzIyMtChQ4cmJ9jSGHO3gblM/yIVydkFeH18F/xtCP9uiMxJoVShx5s/o0qhfeUqOqQNdr4w0EpZEZGxzHa3oFqtxhtvvNHoO0MM3bZJ1jM0whfJ2QU4nF3A4orIzI5fvV+vsAKAF0eZv20NEVmHUcXV0KFD6/V4MiQ2NtaoW2nJMoZ28gP+m4nU60WoUijh5Fh/GR4iMo3/nuWQIFFrY1RxlZSUZKY0yJIi/N0g9XBCrqwKqdeLNJPcici0FEqVzvlW/UPbaK0lSEQtC9/drZBAIMDQTlwKh8jcfr12HxU6hgRnDgqzQjZEZCksrlqpuqtVyZdZXBGZy249dwmOjPS3cCZEZEksrlqpweG+EAqA7Lwy5JToX1OJiJqmRqnCHh3zrXoHe3KeI1ELx+KqlfJyEaNnOy8AwJHsQusmQ9QCnbhehIrq+kOCU/u3t0I2RGRJJimuioqKDK5cTrZpaETtvKvDHBokMrndZ+7q3P5QtwALZ0JEltbk4urixYtYuXIlBg4cCD8/P/j7+2P69OnYuXMnysvLTZkjmUndvKtjVwqhVHEpHCJTUarU2KNjoeYuUnd4uYitkBERWZJRxVVWVhZefvllREREYMCAATh58iSee+455OXlYe/evWjfvj3efvtt+Pr6YuzYsfj000/NlTeZQFSwF9ydHFBcocC5uyXWToeoxUi9XoSyamW97ZP7BlshGyKyNKOKq+PHj6O8vBxr165FYWEhdu7cienTp8PX1xf9+/fHO++8gzNnziAzMxNjxozBrl27zJU3mYCDSIhBHdmSgcjU9A0Jju0utXAmRGQNRq0tSM1nC2sL/tHWE7ew9Ptz6Nu+Db57nuucETWXUqVGrzd/rnflqr23Cw6/OsJKWRFRcxnz/W2SCe0qlQpffPGFKU5FFjbkt0ntp28XQ1alsHI2RPbv1A3dQ4KP92lrhWyIyBpMUlwJhUL8+9//NsWpyMKCvV3QwdcVSpUax6+wJQNRc+3RMyQ4rkeghTMhImsxWZ+rvn37Yt26daY6HVnQ793aWVwRNUeVQoldp+s3DvVzEyPc380KGRGRNZisuLpz5w4++ugjhIaG4sknn8SKFSuwZ88eU52ezOiP6wxyCh5R0/33bA7KdQwJTuzdFgKBwAoZEZE1NLm4unnzptbPP/74I65du4bz58/jpZdegp+fHw4cONDsBMn8BnTwgVgkxJ0HlbheyB5lRE216eg1ndvHdOeQIFFr0uTiKjIyEgkJCSgqKtLa7ubmhpiYGPztb3/DmjVrmpsfWYCL2AF9Q9sAYEsGoqY6fesBLuaU1tvuKhGhd7CX5RMiIqtpcnGVnJyMM2fOoEOHDnjvvfdQWcnFf+3ZkAjOuyJqjk+Trujc/nCPQAiFHBIkak2aXFz169cPiYmJ2L59O3bu3Inw8HB89tlnXGPQTtXNu0q5eh/ymvpzRohIv4JSOfZfzNe5byzvEiRqdZo9oT0+Ph5paWn44IMPsGrVKnTt2pWd2e1QF6kHfN0kqFQokXbzgbXTIbIrHx/Ihq5bQcQiAWI7+lg8HyKyLpPdLfjoo4/iq6++gre3N/7yl7+Y6rRkIUKhAEMj6u4a5NAgUWNV1yix/dRtnftGRgZA4iCycEZEZG0OTT3wiy++wMWLFzWPO3fuAABCQkLw8MMPmyxBspyhnfyw6/RdJGcXYPHYSGunQ2QX1hy4DIVSdwuT8T05JEjUGjW5uFqyZAl69OiB7t27Y9KkSejevTu6d+8OV1dXU+ZHFjT4tytXF3NkKCiVw89dYuWMiGybUqXGF8eu69wnFADDO/tZOCMisgVNLq7y8vJMmQfZAF83CboFeeDCPRmOXC7A433aWTslIpu27uAVVCl038QzoIMP3J0cLZwREdkCo+Zc5ebmQi6XNzr+2jXdDfXIdtUthXOELRmIDKqsVmLj4at69z/cM8iC2RCRLTGquPruu+/g7e2Nxx57DJs3b0ZBQf2GkydOnMDSpUvRrVs39OrVy2SJkmUMjagrrgqgUnEpHCJ9NiRdQaVCd9sSAYDRXQMsmxAR2Qyjiqv58+fjzJkzGDJkCL788ku0a9cOgwcPxnvvvYc5c+YgMDAQEydORH5+PlauXKmz+CLbFt2+DVzFIhSWVeNijsza6RDZpPtlcvz7sP4r84PDfThnkagVM3rOVXh4OBISEpCQkID79+9jz5492Lt3L0JDQ7Fz507ExsZygVI7JnYQIrajDw5k5iP5cgG6t/W0dkpENufjxMuoVupvmDwpOtiC2RCRrWnyhHYA8PHxwYwZMzBjxgxT5UM2YEiEX21xlV2AF4aHWzsdIptyo7AcX/96S+9+R5GAQ4JErZzJmohSy1E3qT3t5gOUy2usnA2RbXn/50tQqvXPRxzbPRCukmb93kpEdo7FFdUT6uOCYG9nKJRq/HrtvrXTIbIZ6bceYO+5XIMxj/Vpa6FsiMhWsbiiegQCgeauweRs3pRABABqtRor9mYajHF3csCQcF8LZUREtorFFelUNzSYzH5XRACAXy7m4eQNw4uaT4wKgoOIH6tErR0/BUin2I4+EAkFuF5YjttFFdZOh8iqFEoV/vW/Sw3G8S5BIgLstLgqKirCtGnT4OHhAS8vL8yePRtlZWUGj6mqqsK8efPg4+MDNzc3TJo0SWsJnzNnzmDq1KkIDg6Gs7MzunTpgo8//ljrHElJSRAIBPUeubmG52DYIw8nR/QJ8QIAHObQILVy20/exrXCchjqMuPvLkGvdmxdQkR2WlxNmzYNFy5cwP79+7Fnzx4kJydj7ty5Bo9ZtGgRdu/ejR07duDw4cO4d+8eHn/8cc3+tLQ0+Pv74z//+Q8uXLiAf/zjH1iyZAnWrVtX71xZWVnIycnRPPz9/U3+Gm3BH7u1E7VWZfIarDmQDQAwcJMg/to3mD3+iAgAIFCrDX1c2J7MzEx07doVJ0+eRN++fQEA+/btw7hx43Dnzh0EBdVfz6ukpAR+fn7YunUrJk+eDAC4dOkSunTpgpSUFAwYMEDnc82bNw+ZmZk4ePAggNorVyNGjMCDBw/g5eXVpPxlMhk8PT1RUlICDw+PJp3DUs7cLsaj64/BXeKA9GWj4ci5JNQKfbQ/G2sTL0PiIIS8Rn/j0KRXhiPU19WCmRGRJRnz/W1335YpKSnw8vLSFFYAEBcXB6FQiBMnTug8Ji0tDQqFAnFxcZptkZGRCAkJQUpKit7nKikpgbe3d73tUVFRCAwMxOjRo3Hs2DGD+crlcshkMq2Hveje1hNtXBxRKq9Bxu1ia6dDZHEKpQqbj10HAIOFVaiPCwsrItKwu+IqNze33jCcg4MDvL299c59ys3NhVgsrne1KSAgQO8xx48fx/bt27WGGwMDA7Fx40bs3LkTO3fuRHBwMIYPH4709HS9+a5YsQKenp6aR3Cw/Ux4FQkFGMyWDNSKZdwuRmlVDcQiw8N902LaWygjIrIHNlNcLV68WOdk8T8+Ll1q+G4dUzh//jweffRRLF++HA899JBme+fOnfHss88iOjoaAwcOxBdffIGBAwdi9erVes+1ZMkSlJSUaB63b9+2xEswmSERtT17WFxRa3Tkt3/3CpXh2RNsHEpEf2QzazS8/PLLmDlzpsGYDh06QCqVIj8/X2t7TU0NioqKIJVKdR4nlUpRXV2N4uJiratXeXl59Y65ePEiRo0ahblz5+L1119vMO/+/fvj6NGjevdLJBJIJJIGz2Or6ia1n71bgqLyani7iq2cEZHlHP6tz5uhmamRUnf4utnve5yITM9miis/Pz/4+fk1GBcbG4vi4mKkpaUhOjoaAHDw4EGoVCrExMToPCY6OhqOjo5ITEzEpEmTANTe8Xfr1i3ExsZq4i5cuICRI0dixowZePfddxuVd0ZGBgIDAxsVa4+knk7oHOCOrLxSHLtSiEd61b9hgKglelBejbN3ihuMmzUo1Oy5EJF9sZniqrG6dOmCMWPGYM6cOdi4cSMUCgXmz5+PJ554QnOn4N27dzFq1Ch89dVX6N+/Pzw9PTF79mwkJCTA29sbHh4eWLBgAWJjYzV3Cp4/fx4jR45EfHw8EhISNHOxRCKRpuhbs2YNwsLC0K1bN1RVVeHzzz/HwYMH8csvv1jnD8NChnbyRVZeKZKzC1hcUatx9EqhwStWdfieIKI/s7viCgC+/vprzJ8/H6NGjYJQKMSkSZOwdu1azX6FQoGsrCxUVPzeWXz16tWaWLlcjvj4eGzYsEGz/7vvvkNBQQH+85//4D//+Y9me/v27XHjxg0AQHV1NV5++WXcvXsXLi4u6NmzJw4cOIARI0aY/0Vb0dBOfvi/I9eRfLkAarWavXyoVWjMPMOebT3hIrbLj1EiMiO763Nl7+ypz1WdKoUSvd76BfIaFX5eOBSdpe7WTonIrNRqNWJXHESurMpg3KfT+mBsj5Y7LYCIftei+1yR5Tk5ihDTwQcA7xqk1uFyflmDhZUAwOiuAZZJiIjsCosrapShdS0ZuBQOtQKNGhJs5wkHrlpARDrwk4EaZVin2kn9qdeLUKVQWjkbIvNqzGLlf4/vbIFMiMgesbiiRgn3d0OgpxPkNSqcuF5k7XSIzKZKoURqA//GhQAGhftaJiEisjssrqhRBAKBpqEo511RS5Z6vcjgOoIA0KOdJ++aJSK9WFxRow3pxKVwqOVrzL/vf4yPtEAmRGSvWFxRow0O94VQUHsn1b3iSmunQ2QWDc23EgLoH8YhQSLSj8UVNZqXixg923kBAI7+tuYaUUuSU1KJy/llBmO6tfW0UDZEZK9YXJFRhv521+BhtmSgFuhIdsO/NCwew7sEicgwFldklGG/zbs6erkQShWb+1PL0tCQoADAoIiGF5gnotaNxRUZpVc7L7g7OaCkUoGzd4qtnQ6RyShVaiRl5xuM4dJPRNQYLK7IKA4iIQZ1rLtrkPOuqOU4d7cE5XLDDXIXjAy3UDZEZM9YXJHR6uZdcSkcakmSsgxftQKAMd25SDMRNYzFFRlt6G/zrjJuF6OkUmHlbIhM43/ncg3uD/NxgUjIxqFE1DAWV2S0dm1c0MHPFUqVGilXOTRI9k9WpUB2XqnBmNlDOlgoGyKydyyuqEnqlsI5zHlX1AIcvVyIhu59ndi7rUVyISL7x+KKmmToH5bCUavZkoHs2670Owb3dwt0h5vEwULZEJG9Y3FFTTKggw/EIiHuFlfiWmG5tdMhajK1Wo2Uq/cNxrwwgncJElHjsbiiJnERO6BvaBsAwBEu5Ex27GpBGcqr9bdgEAIYy7sEicgILK6oyX5vycB5V2S/vv71lsH9fUPbQMi7BInICCyuqMnqJrWnXL0PeY3h5otEturni4ZbMCwe28VCmRBRS8HiiposUuoOXzcJKhVKpN14YO10iIxWWa3EveIqvfsdRQL0ad/GghkRUUvA4oqaTCgUYGhE7V2Dh9mtnezQ1tSbBvf3C/W2UCZE1JKwuKJm0cy7Yr8rskM7ThluwfDmI10tlAkRtSQsrqhZBv925SozR4b8Uv3DK0S2RqlSI8tAV3axSIBOUg8LZkRELQWLK2oWXzcJuret/QI6yrsGyY4cuJgLQ/1v+3JIkIiaiMUVNVvdXYPJ7HdFduTTw1cN7l88trOFMiGilobFFTXbkN+KqyOXC6FScSkcsn1KlRpn75To3e8gFKBnO94lSERNw+KKmi26fRu4ikW4X16Nizkya6dD1KBjVwpg6PeA3iFeFsuFiFoeFlfUbGIHIWI7+gAAktmSgezAxweuGNz/tyFhFsqEiFoiFldkEr+3ZGBxRbZNqVIj43ax3v1CATAqMsByCRFRi8PiikyiblJ72s0HKJfXWDkbIv32X8yF0sBtgpGBHnAQ8aORiJqOnyBkEu19XBDs7QyFUo2Uq/etnQ6RXp8mGb5L8K/R7SyUCRG1VCyuyCQEAsHvLRk474psVI1ShXMG7hIEgId7BVkoGyJqqVhckclw3hXZuq2pt6AysL+tlxN83SQWy4eIWia7LK6Kioowbdo0eHh4wMvLC7Nnz0ZZWZnBY6qqqjBv3jz4+PjAzc0NkyZNQl5enlaMQCCo99i2bZtWTFJSEvr06QOJRILw8HB8+eWXpn55dmtgRx84CAW4cb8Ct+5XWDsdonr+X4rhhZrH9+RVKyJqPrssrqZNm4YLFy5g//792LNnD5KTkzF37lyDxyxatAi7d+/Gjh07cPjwYdy7dw+PP/54vbjNmzcjJydH85g4caJm3/Xr1zF+/HiMGDECGRkZWLhwIf72t7/h559/NvVLtEvuTo7oE1LbeJFDg2RryqpqcDnf8C9h43sEWigbImrJHKydgLEyMzOxb98+nDx5En379gUAfPLJJxg3bhw++OADBAXV/82zpKQEmzZtwtatWzFy5EgAtUVUly5d8Ouvv2LAgAGaWC8vL0ilUp3PvXHjRoSFheHDDz8EAHTp0gVHjx7F6tWrER8fb+qXapeGdvJF6o0iJGcX4KkB7a2dDpHGvxtY7sZVLEKPtp4WyoaIWjK7u3KVkpICLy8vTWEFAHFxcRAKhThx4oTOY9LS0qBQKBAXF6fZFhkZiZCQEKSkpGjFzps3D76+vujfvz+++OILqP9wy3ZKSorWOQAgPj6+3jn+SC6XQyaTaT1asrqlcI5fvQ+F0tDsFiLL+v70XYP7R0T6QygUWCgbImrJ7K64ys3Nhb+/v9Y2BwcHeHt7Izc3V+8xYrEYXl5eWtsDAgK0jnn77bfx7bffYv/+/Zg0aRJeeOEFfPLJJ1rnCQgIqHcOmUyGyspKnc+9YsUKeHp6ah7BwcHGvFy7072tJ9q4OKJMXoPTt4qtnQ4RAKBKocSdYt3v0Tpjuuu+Yk1EZCybKa4WL16sc0L5Hx+XLl0yaw5vvPEGBg0ahN69e+O1117Dq6++ilWrVjXrnEuWLEFJSYnmcfv2bRNla5tEQgEGR/CuQbIt+87r/sWrjkDw+1VXIqLmspk5Vy+//DJmzpxpMKZDhw6QSqXIz8/X2l5TU4OioiK9c6WkUimqq6tRXFysdfUqLy9P7zEAEBMTg3feeQdyuRwSiQRSqbTeHYZ5eXnw8PCAs7OzznNIJBJIJK3r1u6hEb7YfeYejlwuwCvxna2dDhG+STV8l2DvYC94OjtaKBsiaulsprjy8/ODn1/DvznGxsaiuLgYaWlpiI6OBgAcPHgQKpUKMTExOo+Jjo6Go6MjEhMTMWnSJABAVlYWbt26hdjYWL3PlZGRgTZt2miKo9jYWOzdu1crZv/+/QbP0RrV9bs6e7cEReXV8HYVWzkjas3UajXSbhYbjBndlUOCRGQ6NjMs2FhdunTBmDFjMGfOHKSmpuLYsWOYP38+nnjiCc2dgnfv3kVkZCRSU1MBAJ6enpg9ezYSEhJw6NAhpKWlYdasWYiNjdXcKbh79258/vnnOH/+PK5cuYJPP/0U7733HhYsWKB57ueeew7Xrl3Dq6++ikuXLmHDhg349ttvsWjRIsv/QdiwAA8ndA5wh1oNHL1SaO10qJXLyitFjUr/WoIAMCKSQ4JEZDo2c+XKGF9//TXmz5+PUaNGQSgUYtKkSVi7dq1mv0KhQFZWFioqfm9kuXr1ak2sXC5HfHw8NmzYoNnv6OiI9evXY9GiRVCr1QgPD8dHH32EOXPmaGLCwsLw3//+F4sWLcLHH3+Mdu3a4fPPP2cbBh2GdvJFVl4pkrMLMIHLiZAVbTl+w+B+f3cJOge4WyYZImoVBGq1geXhyeRkMhk8PT1RUlICDw8Pa6djNkcuF+DpTakI8JDg1yWjIBDwFneyjl5v/YKSSoXe/U/GhOC9x3pYMCMiskfGfH/b3bAg2Yd+od5wchQiTyZHdp7hrthE5nK/TG6wsAKAEZ39De4nIjIWiysyCydHEWLCfACwJQNZz9cnbhjc7ygSYGBHH8skQ0StBosrMpu6uwa5ziBZy3dphruyD+jgA1eJXU49JSIbxuKKzGZohC8A4MT1IlRWK62cDbU2VQolbhUZ7srOIUEiMgcWV2Q24f5uCPR0QnWNCieu37d2OtTKfH3CcONQoHY9QSIiU2NxRWYjEAgwVLMUDvtdkWVtP2l4qalQHxeE+bpaKBsiak1YXJFZ1c27OsJ5V2RBVQplg3epDueQIBGZCYsrMqtB4T4QCoDL+WW4V2x4/guRqew41fAC6SM5JEhEZsLiiszKy0WMnu28APDqFVnO1hO3DO53dhShf5i3hbIhotaGxRWZnaYlA+ddkQWUVCqQmVtqMGZQuA+cHEUWyoiIWhsWV2R2wzrVtmQ4eqUQygYW0CVqrh8zDPe2AjjfiojMi8UVmV2vdl5wd3JASaUCZ+8UWzsdauG+bmBIEGALBiIyLxZXZHYOIiEGh9deveLQIJlTbkkVshoYEuwc4I62Xs4WyoiIWiMWV2QRXAqHLOGnRg0J+lkgEyJqzVhckUUM+W0pnIzbxSipVFg5G2qptjWiBcND3aQWyISIWjMWV2QR7dq4oIOfK5QqNY5f4dAgmd6V/FJcKyg3GBPgIUHvYC/LJERErRaLK7IYzVI4HBokM/jh9L0GY8Z0k0IoFFggGyJqzVhckcUM+0O/K7WaLRnIdNRqNXadvtNg3NgegRbIhohaOxZXZDExHbwhFglxt7gS1woND98QGSP9VjHuFVcZjPFxE6NfKLuyE5H5sbgii3ERO6BvaBsAQHI2hwbJdBrTODS+mxQiDgkSkQWwuCKL+n0pHBZXZBoKpQq7zzQ832psd94lSESWweKKLKpuUvuv14ogr1FaORtqCY5eLsSDCsPtPTydHTGgg4+FMiKi1o7FFVlUl0B3+LlLUKlQ4tSNB9ZOh1qAHxoxJPhQ1wA4ivhxR0SWwU8bsiiBQKBpKMqWDNRc5fIa/HIhr8G4sT04JEhElsPiiizujy0ZiJpj/8U8VCoMDy+7SkQY9NvalkRElsDiiiyu7osuM0eG/FLDt88TGdKYuwRHdwmAxEFkgWyIiGqxuCKL83WToHtbDwDAEV69oia6XyZH8uWG//2M6c7GoURkWSyuyCq4FA4113/P5UCpMtzpXywSYHhnPwtlRERUi8UVWUVdv6ujlwuhauALkkiXH043PCQY1zUATo4cEiQiy2JxRVbRJ6QNXMUi3C+vxsUcmbXTITtz634F0m8VNxg3jmsJEpEVsLgiqxA7CBHbsbap42F2aycj1U1kN7SajUgAjOjsb6GMiIh+x+KKrIZL4VBTqNVqTeNQQyPKwzr7w1XiYKGsiIh+x+KKrKZuUnvazQcok9dYORuyFxfuyXC1oNzgVSsAeKQXhwSJyDpYXJHVhPq6IsTbBTUqNVKu3rd2OmQn6oYEDU1UFwAYGRlgoYyIiLSxuCKrGtqptqHoEbZkoEZQqtT46cw9AEBFtf7O7NGhbeDp7GiptIiItNhlcVVUVIRp06bBw8MDXl5emD17NsrKygweU1VVhXnz5sHHxwdubm6YNGkS8vJ+X5Psyy+/hEAg0PnIz88HACQlJencn5uba9bX25Jp+l1x3hU1wolr95Enk8PJ0fBH11+jgy2UERFRfXZZXE2bNg0XLlzA/v37sWfPHiQnJ2Pu3LkGj1m0aBF2796NHTt24PDhw7h37x4ef/xxzf4pU6YgJydH6xEfH49hw4bB31/7jqOsrCytuD/vp8aL7egDB6EAN+5X4Nb9CmunQzaubiK7u5Phieqju3JIkIisx+5upcnMzMS+fftw8uRJ9O3bFwDwySefYNy4cfjggw8QFBRU75iSkhJs2rQJW7duxciRIwEAmzdvRpcuXfDrr79iwIABcHZ2hrOzs+aYgoICHDx4EJs2bap3Pn9/f3h5eZnnBbYy7k6O6BPSBqk3inD4cgGe9mlv7ZTIRlUplPjfudqrxIWl1Xrjwv1c0cZVbKm0iIjqsbsrVykpKfDy8tIUVgAQFxcHoVCIEydO6DwmLS0NCoUCcXFxmm2RkZEICQlBSkqKzmO++uoruLi4YPLkyfX2RUVFITAwEKNHj8axY8cM5iuXyyGTybQepK1u3hWHBsmQpKx8lMpr4O0ihqGe/k8NYIFORNZld8VVbm5uvWE4BwcHeHt76537lJubC7FYXO9qU0BAgN5jNm3ahCeffFLralZgYCA2btyInTt3YufOnQgODsbw4cORnp6uN98VK1bA09NT8wgO5lyQP6vrd5Vy9T4USpWVsyFb9cuF2jmS3m6Gr0qN68kWDERkXTZTXC1evFjvhPK6x6VLlyySS0pKCjIzMzF79myt7Z07d8azzz6L6OhoDBw4EF988QUGDhyI1atX6z3XkiVLUFJSonncvn3b3Onbne5Bnmjj4ogyeQ1ON2JJE2p9lCo1DmXV3liSU1ypN87PTQJ/dydLpUVEpJPNzLl6+eWXMXPmTIMxHTp0gFQq1dy9V6empgZFRUWQSqU6j5NKpaiurkZxcbHW1au8vDydx3z++eeIiopCdHR0g3n3798fR48e1btfIpFAIpE0eJ7WTCgUYHCEH3afuYfk7AL0D/O2dkpkY9JvPcCDCgXcJA4GG87+pW87C2ZFRKSbzVy58vPzQ2RkpMGHWCxGbGwsiouLkZaWpjn24MGDUKlUiImJ0Xnu6OhoODo6IjExUbMtKysLt27dQmxsrFZsWVkZvv3223pXrfTJyMhAYCCHIZpraMRv867Y74p0OJBZOyQY5GX4qtTTsZxvRUTWZzNXrhqrS5cuGDNmDObMmYONGzdCoVBg/vz5eOKJJzR3Ct69exejRo3CV199hf79+8PT0xOzZ89GQkICvL294eHhgQULFiA2NhYDBgzQOv/27dtRU1ODp556qt5zr1mzBmFhYejWrRuqqqrw+eef4+DBg/jll18s8tpbsrp5V+fulqCovBrevNuL/iAxs/ZqdXGFQm+Ms6MQgZ7OevcTEVmK3RVXAPD1119j/vz5GDVqFIRCISZNmoS1a9dq9isUCmRlZaGi4ve+SatXr9bEyuVyxMfHY8OGDfXOvWnTJjz++OM6Wy1UV1fj5Zdfxt27d+Hi4oKePXviwIEDGDFihFleZ2sS4OGESKk7LuWW4sjlAjwa1dbaKZGNuHm/HFfyyyASCpBfKtcbN7wz+80RkW0QqNVqQ3c1k4nJZDJ4enqipKQEHh4e1k7Hpry3NxOfJV/D5Oh2+OAvvaydDtmITUev4509F9HR1xVXC8v1xiUmDEVHf3cLZkZErYkx3982M+eKqG4pnCOXC8Can+ok/jbfCgL9MQKAhRUR2QwWV2Qz+oa2gZOjEHkyObLySq2dDtkAWZUCqdeLAAB3H+hvwRAZyMKKiGwHiyuyGU6OIsSE+QBgt3aqdTirADUqNYK9nVFVo7/B7CujO1swKyIiw1hckU2pu2swObvQypmQLagbEgz0MHwX4Cgu1ExENoTFFdmUYb+tM5h6owiV1UorZ0PWVKNU4VBW7RXMBxX67xJ0dxJZKiUiokZhcUU2paOfG4I8nVBdo8KJ6/etnQ5ZUdrNByipVMDL2RHX8vXfJTi1X4gFsyIiaphd9rmilksgEGBIhB+2n7qN5OxC9i5qxRIv1TYO7RnsaXCYeNHoTpZKiajRlEolFAr9TW/J9jg6OkIkMs2VcBZXZHOGdvqtuOJSOK1a3ZI3Egf9H3YCAM5ifoyR7VCr1cjNzUVxcbG1U6Em8PLyglQqhUBgoPdLI/BTiWzO4HBfCAXAlfwy3CuuRJAXlzRpba4XluNaQTkchAJczS/TG9fR39WCWRE1rK6w8vf3h4uLS7O/pMky1Go1KioqkJ9fe8W8uWsGs7gim+Pp4ohewV44fasYydkFeKI/59S0NnV3CfYNbYNfrxXpjVs6JtJSKRE1SKlUagorHx8fa6dDRnJ2rv1FPj8/H/7+/s0aIuSEdrJJv3drZ0uG1qhuSDDM1/CVqZFdpZZIh6hR6uZYubi4WDkTaqq6v7vmzpdjcUU2qa7f1dErhVCquBROa1JSocDJGw8AAOXyGr1xYhE/vsg2cSjQfpnq746fTmSTerXzhLuTA0oqFThzp9ja6ZAFJWXnQ6lSI8LfDad+K7J0GdDR24JZERE1HosrskkOIiEGh9c2FOVSOK1LYmbthNKB4T64V1KlN+65oR0slRIRkVFYXJHN+n0pHBZXrYVCqUJSVm1x1cZFbDC2bygnDBPZkqSkJAgEArahAIsrsmF1xVXG7WKUVLAZX2tw6sYDyKpq4O0qRnZuqd44V7EIYgd+fBGZUm5uLhYsWIAOHTpAIpEgODgYjzzyCBITEwEAZ86cwYQJE+Dv7w8nJyeEhoZiypQpmvYFAwcORE5ODjw9PQH8XmzVPZydndGtWzd89tln9Z77+PHjGDduHNq0aQMnJyf06NEDH330EZRK7WXQ/ng+Dw8P9OvXDz/++GO981VXV2PVqlXo06cPXF1d4enpiV69euH111/HvXv3TP1HVw8/nchmtfVyRkc/V6jUwPGrvGuwNahrwTC8kx+OX9W//NHwzn6WSomoVbhx4waio6Nx8OBBrFq1CufOncO+ffswYsQIzJs3DwUFBRg1ahS8vb3x888/IzMzE5s3b0ZQUBDKy2uXpxKLxTobcGZlZSEnJwcXL17Es88+i+eff15TsAHA999/j2HDhqFdu3Y4dOgQLl26hJdeegn//Oc/8cQTT0Ct1r6pafPmzcjJycGpU6cwaNAgTJ48GefOndPsl8vlGD16NN577z3MnDkTycnJOHfuHNauXYvCwkJ88sknZvyTrMU+V2TThkT44WpBOZIvF2Bsj+Y1dSPbp1nypp0ndp2+qzduQq8gS6VE1CxqtRqVCussQu/sKGr03W8vvPACBAIBUlNT4er6ewuUbt264ZlnnkFSUhJKSkrw+eefw8GhtnQICwvDiBEjNLFJSUkYMWIEHjx4AC8vL812f39/zc8vvvgi1q5di/T0dIwaNQrl5eWYM2cOJkyYoHVF629/+xsCAgIwYcIEfPvtt5gyZYpmX10XdalUinfeeQcff/wxDh06hB49egAAVq9ejaNHj+LUqVPo3bu35riQkBAMGzasXrFmDiyuyKYN6+SHL4/fQHJ2IdRqNW9xbsGuFpThemE5HEUCKJQqg7F92vNOQbIPlQolui772SrPffHteLg0YnmooqIi7Nu3D++++65WYVWnrpipqanB999/j8mTJzfps1itVuPnn3/GrVu3EBMTAwD45ZdfcP/+fbzyyiv14h955BF06tQJ33zzjVZxVaempgabNm0CUHvVrM4333yD0aNHaxVWf2SJ7xEOC5JNi+ngDbFIiLvFlbhaUG7tdMiM6oYEB3TwwZEr+oeBnRyF8HOXWCotohbvypUrUKvViIzUv+LBgAEDsHTpUjz55JPw9fXF2LFjsWrVKuTl5TV4/nbt2sHNzQ1isRjjx4/H8uXLMXToUABAdnY2AKBLly46j42MjNTE1Jk6dSrc3NwgkUiwaNEihIaG4q9//atmf3Z2Njp37qx1zGOPPQY3Nze4ublh4MCBDebcXLxyRTbNReyAfmFtcOzKfSRnFyDc383aKZGZHPitBcOIzv5Y9fMlvXF9QrwslBFR8zk7inDx7XirPXdjNHaY7N1330VCQgIOHjyIEydOYOPGjXjvvfeQnJysGZLT5ciRI3B3d4dcLkdqairmz58Pb29vPP/880bnANQO+8XFxeHatWtYtGgR1q5dC29vw1ezN2zYgPLycqxduxbJycmNfq6m4pUrsnl1S+EkX2ZLhpaquKIaaTdrG4aGeLugUqF/WHBEZICl0iJqNoFAABexg1UejR3+ioiIgEAgwKVL+n+pqePj44O//OUv+OCDD5CZmYmgoCB88MEHBo8JCwtDeHg4unXrhlmzZuHpp5/Gu+++CwDo1KkTACAzM1PnsZmZmZqYOlKpFOHh4XjooYewefNmrTsW615PVlaW1jGBgYEIDw9vsAgzFRZXZPPqWjL8eu0+5DXWmRhK5pWUVQClSo3OAe7IztPfggEAotu3sVBWRK2Dt7c34uPjsX79es2df3+kr2+VWCxGx44ddR5jiEgkQmVlJQDgoYcegre3Nz788MN6cT/99BMuX76MqVOn6j1X//79ER0drSnWgNphw/379+P06dNG5WVKLK7I5kVK3eHnLkGVQmVwORSyX3ULNY/q4o9fLhqew9E10MMSKRG1KuvXr4dSqUT//v2xc+dOXL58GZmZmVi7di1iY2OxZ88ePPXUU9izZw+ys7ORlZWFDz74AHv37sWjjz5q8Nz5+fnIzc3FzZs3sWPHDvy///f/NMe4urri3//+N3788UfMnTsXZ8+exY0bN7Bp0ybMnDkTkydP1ppPpcvChQvx73//G3fv1t5hvGjRIsTGxmLUqFH4+OOPkZ6ejuvXr+Pnn3/G//73P4hEjRsubQ7OuSKbJxAIMCTCF7vS7yI5uwCDflsWh1oGhVKFw7914R8S4Yt/J1/TG+vrJoZTI+eREFHjdejQAenp6Xj33Xfx8ssvIycnB35+foiOjsann36KkJAQuLi44OWXX8bt27chkUgQERGBzz//HE8//bTBc9dNLndwcEBwcDCeffZZvPnmm5r9kydPxqFDh/Duu+9iyJAhqKqqQkREBP7xj39g4cKFDQ5vjhkzBmFhYXj33XexYcMGODk5ITExEWvWrMHmzZuxZMkSqFQqhIWFYezYsVi0aFGz/7waIlBbouEDachkMnh6eqKkpAQeHvwNvLF+zLiLl7ZlIFLqjn0Lh1o7HTKh41cK8eTnJ+DjKsbHU6Lw1BepemMf6hqAz6b3tWB2RI1XVVWF69evIywsDE5OTtZOh5rA0N+hMd/fHBYkuzA43BcCAXAptxT5Mv2L+ZL90dwlGOlvsAVDbQw7sxOR7WNxRXbBx02C7kG161UducylcFoKtVqNxEu1c6ziuvhrel3pE83moURkB1hckd0YElE714otGVqOqwVluHm/AmKREL3aeeGKgUaxQgHQ0Y99zojI9rG4IrtR15LhyOVCqFScKtgS1A0JDujogzN3SgzGSj2cIBJy+SOyfZzKbL9M9XfH4orsRp+QNnAVi1BUXo0L92TWTodMoG4YMK6LP5KzDV+R7B/GIUGybY6OjgCAiooKK2dCTVX3d1f3d9lUbMVAdkPsIERsR18cyMxD8uUC9Gjnae2UqBkelP/elX1kpD8+PXTVYPyoLuzMTrZNJBLBy8tL0y3cxcWFi83bCbVajYqKCuTn58PLy6vZvbBYXJFdGdaptrg6nF2AeSPCrZ0ONcOhrHyo1LVNYgUCAXIauAuUndnJHkilUgDQWo6F7IeXl5fm77A5WFyRXambd5V+8wHK5DVwk/CfsL1K/G2+VVyXABxroAWDg1CAQE/2DSLbJxAIEBgYCH9/fygUCmunQ0ZwdHQ0Wfd2fjORXWnv44oQbxfcKqpAytX7GN2VQ0X2qLrm967so7r44/Mj+ruyA0C7Ns4cXiG7IhKJLLLMCtkmTmgnuzO0028tGRqYAE22K/V6EcrkNfB1k6BnW88Ge5dxySMisid2WVwVFRVh2rRp8PDwgJeXF2bPno2ysjKDx3z22WcYPnw4PDw8IBAIdK7y3Zjznj17FkOGDIGTkxOCg4Px/vvvm/KlUSMMjagdGmS/K/tVt1DzyEg/XCkoh6yqxmA8r1ASkT2xy+Jq2rRpuHDhAvbv3489e/YgOTkZc+fONXhMRUUFxowZg6VLlzb5vDKZDA899BDat2+PtLQ0rFq1Cm+++SY+++wzk702alhsRx84CAW4eb8CN+/rbzpJtumPXdlHdQnA0QbmWwFAVLCXmbMiIjIdu5tzlZmZiX379uHkyZPo27d2AddPPvkE48aNwwcffICgoCCdxy1cuBAAkJSU1OTzfv3116iursYXX3wBsViMbt26ISMjAx999JHe4k4ul0Mul2t+LimpbZQok7FPU3P08Bcj7eYD/Hz6Gp7o397a6ZARruSX4mbOfTg6CNErQIwvDl2ESq6/L5DYQQhhTRVkXFOSiKyo7nu7UY1G1XZm06ZNai8vL61tCoVCLRKJ1Lt27Wrw+EOHDqkBqB88eGD0eZ9++mn1o48+qhVz8OBBNQB1UVGRzudbvny5GgAffPDBBx988NECHrdv326w1rC7K1e5ubnw9/fX2ubg4ABvb2/k5uaa9by5ubkICwvTigkICNDsa9Omfh+eJUuWICEhQfOzSqVCUVERfHx8GnX3U79+/XDy5EmjX09zmfN5TXnu5pyrqccae1xj42UyGYKDg3H79m14eHgYnVdLZ633grGskWdreL829Xi+X62jJb5f1Wo1SktL9Y6Q/ZHNFFeLFy/Gv/71L4MxmZmZFsrGdCQSCSQSidY2Ly+vRh8vEoms8sY15/Oa8tzNOVdTjzX2OGPjPTw8+GGtg7XeC8ayRp6t4f3a1OP5frWOlvp+9fT0bFSczRRXL7/8MmbOnGkwpkOHDpBKpfU639bU1KCoqKhZXVUbc16pVIq8vDytmLqfTdHRVZd58+aZ5bzWfF5Tnrs552rqscYeZ62/w5bGXv4crZFna3i/NvV4vl+tw17+HM2Vp0Cttq/luzMzM9G1a1ecOnUK0dHRAIBffvkFY8aMwZ07dxq8XJeUlIQRI0bgwYMHWleQGnPeTz/9FP/4xz+Ql5enWdRx6dKl2LVrFy5dumSeF0ythkwmg6enJ0pKSuziNz6i1ozvVzLE7loxdOnSBWPGjMGcOXOQmpqKY8eOYf78+XjiiSc0hdXdu3cRGRmJ1NRUzXG5ubnIyMjAlStXAADnzp1DRkYGioqKGn3eJ598EmKxGLNnz8aFCxewfft2fPzxx1pzqoiaSiKRYPny5fWGkYnI9vD9SgY1OOXdBt2/f189depUtZubm9rDw0M9a9YsdWlpqWb/9evX1QDUhw4d0mzTd9fe5s2bG31etVqtPnPmjHrw4MFqiUSibtu2rXrlypXmfrlERERkR+xuWJCIiIjIltndsCARERGRLWNxRURERGRCLK6IiIiITIjFFREREZEJsbgisgOPPfYY2rRpg8mTJ1s7FSJqwO3btzF8+HB07doVPXv2xI4dO6ydElkY7xYksgNJSUkoLS3Fli1b8N1331k7HSIyICcnB3l5eYiKikJubi6io6ORnZ0NV1dXa6dGFsIrV0R2YPjw4XB3d7d2GkTUCIGBgYiKigJQuzSar6+vpmE1tQ4srojMLDk5GY888giCgoIgEAjwww8/1ItZv349QkND4eTkhJiYGK3VBYjIskz5nk1LS4NSqURwcLCZsyZbwuKKyMzKy8vRq1cvrF+/Xuf+7du3IyEhAcuXL0d6ejp69eqF+Pj4eguJE5FlmOo9W1RUhOnTp+Ozzz6zRNpkQzjnisiCBAIBvv/+e0ycOFGzLSYmBv369cO6desAACqVCsHBwViwYAEWL16siUtKSsK6des454rIgpr6npXL5Rg9ejTmzJmDp59+2hqpkxXxyhWRFVVXVyMtLQ1xcXGabUKhEHFxcUhJSbFiZkSkS2Pes2q1GjNnzsTIkSNZWLVSLK6IrKiwsBBKpRIBAQFa2wMCApCbm6v5OS4uDn/5y1+wd+9etGvXjoUXkZU05j177NgxbN++HT/88AOioqIQFRWFc+fOWSNdshIHaydARA07cOCAtVMgokYaPHgwVCqVtdMgK+KVKyIr8vX1hUgkQl5entb2vLw8SKVSK2VFRPrwPUuNweKKyIrEYjGio6ORmJio2aZSqZCYmIjY2FgrZkZEuvA9S43BYUEiMysrK8OVK1c0P1+/fh0ZGRnw9vZGSEgIEhISMGPGDPTt2xf9+/fHmjVrUF5ejlmzZlkxa6LWi+9Zai62YiAys6SkJIwYMaLe9hkzZuDLL78EAKxbtw6rVq1Cbm4uoqKisHbtWsTExFg4UyIC+J6l5mNxRURERGRCnHNFREREZEIsroiIiIhMiMUVERERkQmxuCIiIiIyIRZXRERERCbE4oqIiIjIhFhcEREREZkQiysiIiIiE2JxRURERGRCLK6IiIiITIjFFREREZEJsbgiolZt+PDhEAgEEAgEyMjIsHY6Wl555RVMnDhR7/6ZM2dqcv/hhx8slhcRGcbiiohavTlz5iAnJwfdu3cH8HvR8txzz9WLnTdvHgQCAWbOnGmWXGbNmoXXX38dAJCRkYGePXvqjf3444+Rk5NjljyIqOlYXBFRq1FTU6Nzu4uLC6RSKRwcHDTbgoODsW3bNlRWVmq2VVVVYevWrQgJCTFLfkqlEnv27MGECRMAAGfOnEGvXr30xnt6ekIqlZolFyJqOhZXRNQi3bhxAwKBAN9++y2GDBkCiUSCn376qdHH9+nTB8HBwdi1a5dm265duxASEoLevXtrxQ4fPhzz58/H/Pnz4enpCV9fX7zxxhtQq9WaGJVKhffffx/h4eGQSCQICQnBu+++q3We48ePw9HREf369cOdO3dQWFgIABg9ejRcXFzQuXNnnDhxoil/HERkQSyuiKhFOnPmDABg1apVWLZsGS5cuIBRo0YZdY5nnnkGmzdv1vz8xRdfYNasWTpjt2zZAgcHB6SmpuLjjz/GRx99hM8//1yzf8mSJVi5ciXeeOMNXLx4EVu3bkVAQIDWOX766Sc88sgjWvO/1q9fj6VLl+LMmTMICQnB4sWLjXoNRGR5Dg2HEBHZn4yMDLi6umLHjh0IDQ1t0jmeeuopLFmyBDdv3gQAHDt2DNu2bUNSUlK92ODgYKxevRoCgQCdO3fGuXPnsHr1asyZMwelpaX4+OOPsW7dOsyYMQMA0LFjRwwePFjrHD/++CNWr16tyd/b2xvffvstfH19AQATJkzAv//97ya9FiKyHBZXRNQinTlzBhMmTGhyYQUAfn5+GD9+PL788kuo1WqMHz9eU+j82YABAyAQCDQ/x8bG4sMPP4RSqURmZibkcrnBK2eZmZm4d++eJiYjIwOPPvqo1vNdv34d4eHhTX49RGQZHBYkohYpIyMDw4cPb/Z5nnnmGXz55ZfYsmULnnnmmSadw9nZucGYn376CaNHj4aTkxOA2vwHDBigFZORkYGoqKgm5UBElsPiiohaHJlMhhs3btSbeN4UY8aMQXV1NRQKBeLj4/XG/Xmi+a+//oqIiAiIRCJERETA2dkZiYmJeo//8ccf8eijjwIASktLce3atXr5s7gisg8cFiSiFufMmTMQiUTo0aNHs88lEomQmZmp+X99bt26hYSEBDz77LNIT0/HJ598gg8//BAA4OTkhNdeew2vvvoqxGIxBg0ahIKCAly4cAGzZ89Gfn4+Tp06pbmbUVf+N2/exIMHD1hcEdkBFldE1OKcOXMGnTt31gyxNZeHh0eDMdOnT0dlZSX69+8PkUiEl156CXPnztXsf+ONN+Dg4IBly5bh3r17CAwM1DQp3b17N/r376+ZX5WRkVEv/9OnT8PLy6tZc8iIyDIE6j82YiEiamWGDx+OqKgorFmzxmrnmDBhAgYPHoxXX321SccLBAJ8//33BpfKISLL4ZwrImr1NmzYADc3N5w7d84qzz948GBMnTrV6OOee+45uLm5mSEjImoOXrkiolbt7t27miVuQkJCIBaLjT6HKa5+NUV+fj5kMhkAIDAwEK6urhZ9fiLSjcUVERERkQlxWJCIiIjIhFhcEREREZkQiysiIiIiE2JxRURERGRCLK6IiIiITIjFFREREZEJsbgiIiIiMiEWV0REREQmxOKKiIiIyIRYXBERERGZ0P8H8/AT2EPYfHsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_field_shells(\"radvel\", \"SPH\", 1024, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enclosed mass "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'a' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[2], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43ma\u001b[49m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'a' is not defined"
     ]
    }
   ],
   "source": [
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_enclosed_overdensity(to_save=True):\n",
    "    with plt.style.context(\"science\"):\n",
    "    # if True:\n",
    "        cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "        fig, axs = plt.subplots(2, 1, sharex=True, figsize=(2 * 3.5, 2 * 1.5 * 2.625), gridspec_kw={\"height_ratios\": [1, 0.8]})\n",
    "        fig.subplots_adjust(wspace=0, hspace=0)\n",
    "\n",
    "        # CSiBORG2 main\n",
    "        d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_csiborg2_main.npz\")\n",
    "        V = 4 / 3 * numpy.pi * d[\"distances\"]**3\n",
    "        V35 = 4 / 3 * numpy.pi * 135**3\n",
    "        rho_mean = 0.3111 * 277.53662724583074\n",
    "        boxsize = csiborgtools.simname2boxsize(\"csiborg2_main\")\n",
    "\n",
    "        dist = d[\"distances\"]\n",
    "        density = d[\"enclosed_mass\"] / V * 1e-9 / rho_mean - 1\n",
    "\n",
    "        density135 = d[\"mass135\"] / V35 * 1e-9 / rho_mean - 1\n",
    "        densitytot = d[\"masstot\"] / boxsize**3 * 1e-9 / rho_mean - 1\n",
    "\n",
    "        print(f\"CSiBORG2_main overdensity within 135 Mpc / h: {numpy.mean(density135)} +- {numpy.std(density135)}\")\n",
    "        print(f\"CSiBORG2_main density of the entire box:      {numpy.mean(densitytot)} +- {numpy.std(densitytot)}\")\n",
    "        mu = numpy.mean(density, axis=0)\n",
    "        y = numpy.copy(density)\n",
    "        for i in range(len(density)):\n",
    "            axs[0].plot(dist, density[i], c=cols[0], alpha=0.25, ls=\"dashed\")\n",
    "            y[i] /= mu\n",
    "            y[i] -= 1\n",
    "        axs[0].plot(dist, mu, c=cols[0], label=\"CB2_main\")\n",
    "        mu2, std2 = numpy.mean(y, axis=0), numpy.std(y, axis=0)\n",
    "        axs[1].fill_between(dist, mu2 - std2, mu2 + std2, alpha=0.5, color=cols[0])\n",
    "\n",
    "        # CSiBORG2 varysmall\n",
    "        d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_csiborg2_varysmall.npz\")\n",
    "        V = 4 / 3 * numpy.pi * d[\"distances\"]**3\n",
    "        V35 = 4 / 3 * numpy.pi * 135**3\n",
    "        rho_mean = 0.3111 * 277.53662724583074\n",
    "        boxsize = csiborgtools.simname2boxsize(\"csiborg2_varysmall\")\n",
    "\n",
    "        dist = d[\"distances\"]\n",
    "        density = d[\"enclosed_mass\"] / V * 1e-9 / rho_mean - 1\n",
    "\n",
    "        density135 = d[\"mass135\"] / V35 * 1e-9 / rho_mean - 1\n",
    "        densitytot = d[\"masstot\"] / boxsize**3 * 1e-9 / rho_mean - 1\n",
    "\n",
    "        print(f\"CSiBORG2_varysmall overdensity within 135 Mpc / h: {numpy.mean(density135)} +- {numpy.std(density135)}\")\n",
    "        print(f\"CSiBORG2_varysmall density of the entire box:      {numpy.mean(densitytot)} +- {numpy.std(densitytot)}\")\n",
    "        mu = numpy.mean(density, axis=0)\n",
    "        y = numpy.copy(density)\n",
    "        for i in range(len(density)):\n",
    "            axs[0].plot(dist, density[i], c=cols[2], alpha=0.25, ls=\"dashed\")\n",
    "            y[i] /= mu\n",
    "            y[i] -= 1\n",
    "        axs[0].plot(dist, mu, c=cols[2], label=\"CB2_varysmall\")\n",
    "        mu2, std2 = numpy.mean(y, axis=0), numpy.std(y, axis=0)\n",
    "        axs[1].fill_between(dist, mu2 - std2, mu2 + std2, alpha=0.5, color=cols[2])\n",
    "\n",
    "        # CSiBORG2 random\n",
    "        d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_csiborg2_random.npz\")\n",
    "        V = 4 / 3 * numpy.pi * d[\"distances\"]**3\n",
    "        V35 = 4 / 3 * numpy.pi * 135**3\n",
    "        rho_mean = 0.3111 * 277.53662724583074\n",
    "        boxsize = csiborgtools.simname2boxsize(\"csiborg2_random\")\n",
    "\n",
    "        dist = d[\"distances\"]\n",
    "        density = d[\"enclosed_mass\"] / V * 1e-9 / rho_mean - 1\n",
    "\n",
    "        density135 = d[\"mass135\"] / V35 * 1e-9 / rho_mean - 1\n",
    "        densitytot = d[\"masstot\"] / boxsize**3 * 1e-9 / rho_mean - 1\n",
    "\n",
    "        print(f\"CSiBORG2_random overdensity within 135 Mpc / h: {numpy.mean(density135)} +- {numpy.std(density135)}\")\n",
    "        print(f\"CSiBORG2_random density of the entire box:      {numpy.mean(densitytot)} +- {numpy.std(densitytot)}\")\n",
    "        for i in range(len(density)):\n",
    "            axs[0].plot(dist, density[i], c=cols[1], alpha=0.5, ls=\"dashed\", zorder=0)\n",
    "        mu = numpy.mean(density, axis=0)\n",
    "        std = numpy.std(density, axis=0)\n",
    "        axs[0].plot(dist, mu, c=cols[1], label=\"CB2_random\", zorder=0)\n",
    "        axs[0].fill_between(dist, mu - std, mu + std, alpha=1/3, color=cols[1], zorder=0)\n",
    "\n",
    "        # CSiBORG1\n",
    "        d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_csiborg1.npz\")\n",
    "        V = 4 / 3 * numpy.pi * d[\"distances\"]**3\n",
    "        V35 = 4 / 3 * numpy.pi * 135**3\n",
    "        rho_mean = 0.307 * 277.53662724583074\n",
    "        boxsize = csiborgtools.simname2boxsize(\"csiborg1\")\n",
    "\n",
    "        dist = d[\"distances\"]\n",
    "        density = d[\"enclosed_mass\"] / V * 1e-9 / rho_mean - 1\n",
    "\n",
    "        density135 = d[\"mass135\"] / V35 * 1e-9 / rho_mean - 1\n",
    "        densitytot = d[\"masstot\"] / boxsize**3 * 1e-9 / rho_mean - 1\n",
    "\n",
    "        print(f\"CSiBORG1 overdensity within 135 Mpc / h:    {numpy.mean(density135)} +- {numpy.std(density135)}\")\n",
    "        print(f\"CSiBORG1 density of the entire box:         {numpy.mean(densitytot)} +- {numpy.std(densitytot)}\")\n",
    "        mu = numpy.mean(density, axis=0)\n",
    "        y = numpy.copy(density)\n",
    "        for i in range(len(density)):\n",
    "            axs[0].plot(dist, density[i], c=cols[3], alpha=0.25, ls=\"dashed\", zorder=0.5)\n",
    "            y[i] /= mu\n",
    "            y[i] -= 1\n",
    "        axs[0].plot(dist, mu, c=cols[3], label=\"CB1\", zorder=0.5)\n",
    "        mu2, std2 = numpy.mean(y, axis=0), numpy.std(y, axis=0)\n",
    "        axs[1].fill_between(dist, mu2 - std2, mu2 + std2, alpha=0.5, color=cols[3], zorder=0.5)\n",
    "\n",
    "        \n",
    "        # BORG2\n",
    "        d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_borg2.npz\")\n",
    "        V = d[\"enclosed_volume\"]\n",
    "        dist = d[\"distances\"]\n",
    "        rho_mean = 0.3111 * 277.53662724583074\n",
    "        density = d[\"enclosed_mass\"] / V / rho_mean - 1\n",
    "\n",
    "        mu = numpy.mean(density, axis=0)\n",
    "        y = numpy.copy(density)\n",
    "        for i in range(len(density)):\n",
    "            axs[0].plot(dist, density[i], c=cols[4], alpha=0.25, ls=\"dashed\", zorder=0.5)\n",
    "            y[i] /= mu\n",
    "            y[i] -= 1\n",
    "        axs[0].plot(dist, mu, c=cols[4], label=\"B2\", zorder=0.5)\n",
    "        mu2, std2 = numpy.mean(y, axis=0), numpy.std(y, axis=0)\n",
    "        axs[1].fill_between(dist, mu2 - std2, mu2 + std2, alpha=0.5, color=cols[4], zorder=0.5)\n",
    "\n",
    "        # BORG1\n",
    "        d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_borg1.npz\")\n",
    "        V = d[\"enclosed_volume\"]\n",
    "        dist = d[\"distances\"]\n",
    "        rho_mean = 0.307 * 277.53662724583074\n",
    "        density = d[\"enclosed_mass\"] / V / rho_mean - 1\n",
    "\n",
    "        mu = numpy.mean(density, axis=0)\n",
    "        y = numpy.copy(density)\n",
    "        for i in range(len(density)):\n",
    "            axs[0].plot(dist, density[i], c=cols[4], alpha=0.25, ls=\"dashed\", zorder=0.5)\n",
    "            y[i] /= mu\n",
    "            y[i] -= 1\n",
    "        axs[0].plot(dist, mu, c=cols[4], label=\"B2\", zorder=0.5)\n",
    "        mu2, std2 = numpy.mean(y, axis=0), numpy.std(y, axis=0)\n",
    "        axs[1].fill_between(dist, mu2 - std2, mu2 + std2, alpha=0.5, color=cols[4], zorder=0.5)\n",
    "\n",
    "        # Plot settings\n",
    "        axs[0].set_ylim(-0.15, 0.15)\n",
    "        axs[1].set_ylim(-0.4, 0.4)\n",
    "        axs[1].set_xlim(-1, dist.max())\n",
    "        axs[1].set_xlabel(r\"$r ~ [\\mathrm{Mpc} / h]$\")\n",
    "        axs[0].set_ylabel(r\"$\\delta_r$\")\n",
    "        axs[1].set_ylabel(r\"$\\delta_r / \\langle \\delta_r \\rangle - 1$\")\n",
    "\n",
    "        axs[0].legend(fontsize=\"small\")\n",
    "        for i in range(2):\n",
    "            axs[i].axvline(135, c=\"k\", ls=\"--\")\n",
    "        fig.tight_layout(h_pad=0)\n",
    "\n",
    "        if to_save:\n",
    "            fout = join(plt_utils.fout, f\"enclosed_overdensity.png\")\n",
    "            print(f\"Saving to `{fout}`.\")\n",
    "            fig.savefig(fout, dpi=plt_utils.dpi, bbox_inches=\"tight\")\n",
    "\n",
    "        fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CSiBORG2_main overdensity within 135 Mpc / h: -0.04871875932547955 +- 0.00201984793460838\n",
      "CSiBORG2_main density of the entire box:      -5.9114960471373655e-05 +- 2.4011781815055744e-07\n",
      "CSiBORG2_varysmall overdensity within 135 Mpc / h: -0.050013774679444456 +- 7.933447782885623e-05\n",
      "CSiBORG2_varysmall density of the entire box:      -5.8827166581804094e-05 +- 1.3328610811057487e-07\n",
      "CSiBORG2_random overdensity within 135 Mpc / h: 0.006494700117062541 +- 0.03838291159994349\n",
      "CSiBORG2_random density of the entire box:      -5.874905109742312e-05 +- 4.066696396103268e-07\n",
      "CSiBORG1 overdensity within 135 Mpc / h:    -0.05300083843550532 +- 0.008948856840208618\n",
      "CSiBORG1 density of the entire box:         0.0008964738425259703 +- 8.270101657242762e-07\n",
      "Saving to `../plots/enclosed_overdensity.png`.\n"
     ]
    }
   ],
   "source": [
    "plot_enclosed_overdensity(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CSiBORG1\n",
    "d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_csiborg1.npz\")\n",
    "V = 4 / 3 * numpy.pi * d[\"distances\"]**3\n",
    "V35 = 4 / 3 * numpy.pi * 135**3\n",
    "rho_mean = 0.307 * 277.53662724583074\n",
    "boxsize = csiborgtools.simname2boxsize(\"csiborg1\")\n",
    "\n",
    "dist = d[\"distances\"]\n",
    "density = d[\"enclosed_mass\"] / V * 1e-9 / rho_mean - 1\n",
    "\n",
    "density135 = d[\"mass135\"] / V35 * 1e-9 / rho_mean - 1\n",
    "densitytot = d[\"masstot\"] / boxsize**3 * 1e-9 / rho_mean - 1\n",
    "\n",
    "\n",
    "density135_csiborg1 = density135\n",
    "\n",
    "\n",
    "# CSiBORG2 main\n",
    "d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_csiborg2_main.npz\")\n",
    "V = 4 / 3 * numpy.pi * d[\"distances\"]**3\n",
    "V35 = 4 / 3 * numpy.pi * 135**3\n",
    "rho_mean = 0.3111 * 277.53662724583074\n",
    "boxsize = csiborgtools.simname2boxsize(\"csiborg2_main\")\n",
    "\n",
    "dist = d[\"distances\"]\n",
    "density = d[\"enclosed_mass\"] / V * 1e-9 / rho_mean - 1\n",
    "\n",
    "density135 = d[\"mass135\"] / V35 * 1e-9 / rho_mean - 1\n",
    "densitytot = d[\"masstot\"] / boxsize**3 * 1e-9 / rho_mean - 1\n",
    "\n",
    "density135_csiborg2 = density135"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.008948856840208618, 0.00201984793460838)"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numpy.std(density135_csiborg1), numpy.std(density135_csiborg2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([15517, 15617, 15717, 15817, 15917, 16017, 16117, 16217, 16317,\n",
       "       16417, 16517, 16617, 16717, 16817, 16917, 17017, 17117, 17217,\n",
       "       17317, 17417])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)\n",
    "paths.get_ics(\"csiborg2_main\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAGdCAYAAADjWSL8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAs10lEQVR4nO3de1xU9b7/8TcgVxFGEBhQtIuWlzQ8pkRWupMTmnW6+DiPMrt79LTDXW486rFTlpabvcvSXWGeUxntvbPsftG9vZvuCrXsWKZm4tY0c7AkQEUB4fv7w8P8nARlYIbhC6/n4zGPnLW+a63vZxYz827Nd60VZIwxAgAAsExwoDsAAADQGIQYAABgJUIMAACwEiEGAABYiRADAACsRIgBAABWIsQAAAArEWIAAICV2gW6A41RU1OjH374QR06dFBQUFCguwMAABrAGKPDhw8rJSVFwcFNP45iZYj54YcflJqaGuhuAACARti3b5+6dOnS5PVYGWI6dOgg6eSLEBMTE+DeAACAhigrK1Nqaqr7e7yprAwxtT8hxcTEEGIAALCMr4aCMLAXAABYiRADAACsRIgBAABWsnJMDAAAzcUYoxMnTqi6ujrQXWnxQkJC1K5du2a7/AkhBgCAelRWVurAgQMqLy8PdFesERUVpeTkZIWFhfl9W4QYAADqUFNTo927dyskJEQpKSkKCwvjAqtnYIxRZWWlfvzxR+3evVs9evTwyQXtzoQQAwBAHSorK1VTU6PU1FRFRUUFujtWiIyMVGhoqL777jtVVlYqIiLCr9tjYC8AAGfg76MJrU1zvl7sGQAAYCV+TgIAwAv7S47p56OVzba9ju3D1NkR2WzbswkhBgCABtpfckyZT63VsarmO906MjREKycNIcjUgRADAEAD/Xy0UseqqjX35jR1T4z2+/YKDx7RxEWb9fPRSq9DjMvl0qxZs7RkyRLt379fiYmJSktL08SJEzVs2DCdc845+u677ySdHMeSlJSkESNGaPbs2erYsaMk6fjx47r33nu1adMmbd++Xddee63ee+89X5fZaIQYAAC81D0xWhd1jg10N+q1Z88eDR48WA6HQ08++aT69u2rqqoqLVu2TNnZ2frmm28kSTNnztS4ceNUXV2tb7/9VuPHj9f999+vP//5z5Kk6upqRUZG6v7779fbb78dyJLqRIgBAF8q2SeVH2r88lHxkiPVd/1Bm3TfffcpKChIGzduVPv27d3T+/Tpo3vuucf9vEOHDnI6nZKkzp07684779Rrr73mnt++fXs9//zzkqRPPvlEJSUlzVNAAxFiAMBXSvZJeYOkqiZc3TU0SsreSJBBoxUXF2vp0qWaNWuWR4Cp5XA46lxu//79+vDDD5Wenu7nHvoOIQYAfKX80MkAc9MLUqcLvF/+p2+ld8adXA8hBo1UWFgoY4x69ux51rZTp07VQw89pOrqah0/flzp6el6+umnm6GXvkGIAQBf63SBlJIW6F6gjTLGNLjt5MmTddddd8kYo3379unBBx/UyJEjtW7dOoWEhPixl77Bxe4AAGhFevTooaCgIPfg3TPp1KmTunfvrh49euiqq67S3Llz9emnn2rNmjXN0NOmI8QAANCKxMXFKSsrS3l5eTp69Ohp8880OLf26MuxY8f81T2f4uckAAC8VHjwSIveTl5engYPHqxBgwZp5syZ6tevn06cOKEVK1bo+eef1/bt2yVJhw8flsvlcv+cNGXKFCUkJOiyyy5zr2vbtm2qrKxUcXGxDh8+rM2bN0uS0tLSmlpekxFiAABooI7twxQZGqKJizY32zYjQ0PUsX2YV8ucd955+uKLLzRr1ixNmjRJBw4cUEJCggYMGOA+ZVqSpk+frunTp0uSEhISNHDgQC1fvlzx8fHuNtdcc437oniS1L9/f0nejb3xF0IMAAAN1NkRqZWThlhx76Tk5GQ999xzeu655+qcv2fPngatp6HtAoEQAwCAFzo7IrmPUQvBwF4AAGAlQgwAALASIQYAAFiJEAMAAKzkVYjJzc3VwIED1aFDByUmJuqGG27Qjh07PNoMHTpUQUFBHo97773Xo83evXs1cuRIRUVFKTExUZMnT9aJEyeaXg0AAGgzvDo7ae3atcrOztbAgQN14sQJPfjgg7r66qu1bds2jztljhs3TjNnznQ/j4qKcv+7urpaI0eOlNPp1KeffqoDBw7ojjvuUGhoqH73u9/5oCQAANAWeBVili5d6vE8Pz9fiYmJ2rRpk6688kr39KioKDmdzjrXsXz5cm3btk0rV65UUlKS0tLS9Nhjj2nq1Kl69NFHFRbm3QV9AABA29Sk68SUlpZKOnmfhlO9+uqr+stf/iKn06nrrrtODz/8sPtoTEFBgfr27aukpCR3+6ysLP3617/W1q1b3VcCPFVFRYUqKircz8vKyprSbQAAGq9kn1R+qPm2FxUvOVKbb3sWaXSIqamp0cSJEzV48GBddNFF7um33nqrunXrppSUFH311VeaOnWqduzYoXfeeUeS5HK5PAKMJPdzl8tV57Zyc3M1Y8aMxnYVAADfKNkn5Q2Sqsqbb5uhUVL2Rq+DjMvl0qxZs7RkyRLt379fiYmJSktL08SJEzVs2DCdc8457tsJBAcHKykpSSNGjNDs2bPVsWNHSdJHH32kOXPmaOPGjSorK1OPHj00efJkjRkzxudlNkajQ0x2dra+/vprffzxxx7Tx48f7/533759lZycrGHDhmnXrl06//zzG7WtadOmKScnx/28rKxMqamkUgBAMys/dDLA3PSC1OkC/2/vp2+ld8ad3K4XIWbPnj0aPHiwHA6HnnzySfXt21dVVVVatmyZsrOz9c0330iSZs6cqXHjxqm6ulrffvutxo8fr/vvv19//vOfJUmffvqp+vXrp6lTpyopKUmLFy/WHXfcodjYWF177bV+KdkbjQoxEyZM0OLFi7Vu3Tp16dLljG3T09MlSYWFhTr//PPldDq1ceNGjzZFRUWSVO84mvDwcIWHhzemqwAA+F6nC6SUtED3ol733XefgoKCtHHjRo8Tb/r06aN77rnH/bxDhw7u797OnTvrzjvv1Guvveae/+CDD3qs94EHHtDy5cv1zjvvtIgQ49Up1sYYTZgwQe+++65Wr16tc88996zL1N6yOzk5WZKUkZGhLVu26ODBg+42K1asUExMjHr37u1NdwAAwC8UFxdr6dKlys7O9ggwtRwOR53L7d+/Xx9++KH74EN9SktLTxsLGyhehZjs7Gz95S9/0cKFC9WhQwe5XC65XC4dO3ZMkrRr1y499thj2rRpk/bs2aMPPvhAd9xxh6688kr169dPknT11Verd+/euv322/Xll19q2bJleuihh5Sdnc3RFgAAmqiwsFDGGPXs2fOsbadOnaro6GhFRkaqS5cuCgoK0tNPP11v+zfeeEOfffaZ7r77bl92udG8CjHPP/+8SktLNXToUCUnJ7sfixYtkiSFhYVp5cqVuvrqq9WzZ09NmjRJo0aN0ocffuheR0hIiBYvXqyQkBBlZGTotttu0x133OFxXRkAANA4xpgGt508ebI2b96sr776SqtWrZIkjRw5UtXV1ae1XbNmje6++2698MIL6tOnj8/62xRejYk52wuTmpqqtWvXnnU93bp101//+ldvNg0AABqgR48eCgoKcg/ePZNOnTqpe/fu7uXmzp2rjIwMrVmzRpmZme52a9eu1XXXXac5c+bojjvu8FvfvcW9kwAAaEXi4uKUlZWlvLw8HT169LT5JSUl9S4bEhIiSe5hItLJ06xHjhypP/zhDx5nILcEhBgAAFqZvLw8VVdXa9CgQXr77be1c+dObd++Xc8884wyMjLc7Q4fPiyXy6UDBw5o48aNmjx5shISEnTZZZdJOvkT0siRI3X//fdr1KhR7rGwxcXFgSrNQ5Ou2AsAQJv007ctejvnnXeevvjiC82aNUuTJk3SgQMHlJCQoAEDBuj55593t5s+fbqmT58uSUpISNDAgQO1fPlyxcfHS5JeeeUVlZeXKzc3V7m5ue7lhgwZoo8++qjxdfkIIQYAgIaKij95Bd13xjXfNkOjTm7XS8nJyXruuef03HPP1Tl/z549Z11Hfn6+8vPzvd52cyHEAADQUI7Uk7cA4N5JLQIhBgAAbzhSCRUtBAN7AQCAlQgxAADASoQYAABgJUIMAABn4M1l/NG8rxchBgCAOoSGhkqSysvLA9wTu9S+XrWvnz9xdhIAAHUICQmRw+HQwYMHJUlRUVEKCgoKcK9aLmOMysvLdfDgQTkcDvctDPyJEAMAQD2cTqckuYMMzs7hcLhfN38jxAAAUI+goCAlJycrMTFRVVVVge5OixcaGtosR2BqEWIAADiLkJCQZv1yRsMwsBcAAFiJEAMAAKxEiAEAAFYixAAAACsRYgAAgJUIMQAAwEqEGAAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAlQgxAADASoQYAABgJUIMAACwEiEGAABYiRADAACsRIgBAABWIsQAAAArEWIAAICVCDEAAMBKhBgAAGAlQgwAALASIQYAAFiJEAMAAKxEiAEAAFYixAAAACsRYgAAgJUIMQAAwEqEGAAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAlQgxAADASoQYAABgJUIMAACwEiEGAABYyasQk5ubq4EDB6pDhw5KTEzUDTfcoB07dni0OX78uLKzsxUfH6/o6GiNGjVKRUVFHm327t2rkSNHKioqSomJiZo8ebJOnDjR9GoAAECb4VWIWbt2rbKzs7V+/XqtWLFCVVVVuvrqq3X06FF3m9/+9rf68MMP9eabb2rt2rX64YcfdNNNN7nnV1dXa+TIkaqsrNSnn36qV155Rfn5+Zo+fbrvqgIAAK1eO28aL1261ON5fn6+EhMTtWnTJl155ZUqLS3VSy+9pIULF+qqq66SJL388svq1auX1q9fr0svvVTLly/Xtm3btHLlSiUlJSktLU2PPfaYpk6dqkcffVRhYWG+qw4AALRaTRoTU1paKkmKi4uTJG3atElVVVXKzMx0t+nZs6e6du2qgoICSVJBQYH69u2rpKQkd5usrCyVlZVp69atdW6noqJCZWVlHg8AANC2NTrE1NTUaOLEiRo8eLAuuugiSZLL5VJYWJgcDodH26SkJLlcLnebUwNM7fzaeXXJzc1VbGys+5GamtrYbgMAgFai0SEmOztbX3/9tV5//XVf9qdO06ZNU2lpqfuxb98+v28TAAC0bF6Niak1YcIELV68WOvWrVOXLl3c051OpyorK1VSUuJxNKaoqEhOp9PdZuPGjR7rqz17qbbNL4WHhys8PLwxXQUAAK2UV0dijDGaMGGC3n33Xa1evVrnnnuux/wBAwYoNDRUq1atck/bsWOH9u7dq4yMDElSRkaGtmzZooMHD7rbrFixQjExMerdu3dTagEAAG2IV0disrOztXDhQr3//vvq0KGDewxLbGysIiMjFRsbq7FjxyonJ0dxcXGKiYnRb37zG2VkZOjSSy+VJF199dXq3bu3br/9dj3xxBNyuVx66KGHlJ2dzdEWAADQYF6FmOeff16SNHToUI/pL7/8su666y5J0pw5cxQcHKxRo0apoqJCWVlZmjdvnrttSEiIFi9erF//+tfKyMhQ+/btdeedd2rmzJlNqwQAALQpXoUYY8xZ20RERCgvL095eXn1tunWrZv++te/erNpAAAAD9w7CQAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAlQgxAADASoQYAABgJUIMAACwEiEGAABYiRADAACsRIgBAABWIsQAAAArEWIAAICVCDEAAMBKhBgAAGAlQgwAALASIQYAAFiJEAMAAKxEiAEAAFYixAAAACsRYgAAgJUIMQAAwEqEGAAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAlQgxAADASoQYAABgJUIMAACwEiEGAABYiRADAACsRIgBAABWIsQAAAArEWIAAICVCDEAAMBKhBgAAGAlQgwAALASIQYAAFiJEAMAAKxEiAEAAFYixAAAACsRYgAAgJUIMQAAwEqEGAAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAlQgxAADASoQYAABgJUIMAACwktchZt26dbruuuuUkpKioKAgvffeex7z77rrLgUFBXk8hg8f7tGmuLhYY8aMUUxMjBwOh8aOHasjR440qRAAANC2eB1ijh49qosvvlh5eXn1thk+fLgOHDjgfrz22mse88eMGaOtW7dqxYoVWrx4sdatW6fx48d733sAANBmtfN2gREjRmjEiBFnbBMeHi6n01nnvO3bt2vp0qX67LPPdMkll0iSnn32WV1zzTWaPXu2UlJSvO0SAABog7wOMQ3x0UcfKTExUR07dtRVV12lxx9/XPHx8ZKkgoICORwOd4CRpMzMTAUHB2vDhg268cYbT1tfRUWFKioq3M/Lysr80W0AsF/JPqn8UOOWjYqXHKm+7Q/gRz4PMcOHD9dNN92kc889V7t27dKDDz6oESNGqKCgQCEhIXK5XEpMTPTsRLt2iouLk8vlqnOdubm5mjFjhq+7CgCtS8k+KW+QVFXeuOVDo6TsjQQZWMPnIeaWW25x/7tv377q16+fzj//fH300UcaNmxYo9Y5bdo05eTkuJ+XlZUpNZU3GQB4KD90MsDc9ILU6QLvlv3pW+mdcSfXQYiBJfzyc9KpzjvvPHXq1EmFhYUaNmyYnE6nDh486NHmxIkTKi4urnccTXh4uMLDw/3dVQBoHTpdIKWkBboXgN/5/Tox33//vQ4dOqTk5GRJUkZGhkpKSrRp0yZ3m9WrV6umpkbp6en+7g4AAGglvD4Sc+TIERUWFrqf7969W5s3b1ZcXJzi4uI0Y8YMjRo1Sk6nU7t27dKUKVPUvXt3ZWVlSZJ69eql4cOHa9y4cZo/f76qqqo0YcIE3XLLLZyZBAAAGszrIzGff/65+vfvr/79+0uScnJy1L9/f02fPl0hISH66quv9C//8i+64IILNHbsWA0YMEB///vfPX4OevXVV9WzZ08NGzZM11xzjS6//HL9z//8j++qAgAArZ7XR2KGDh0qY0y985ctW3bWdcTFxWnhwoXebhoAAMCNeycBAAArEWIAAICVCDEAAMBKhBgAAGAlQgwAALASIQYAAFiJEAMAAKxEiAEAAFYixAAAACsRYgAAgJUIMQAAwEqEGAAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAlQgxAADASoQYAABgJUIMAACwEiEGAABYiRADAACsRIgBAABWIsQAAAArEWIAAICVCDEAAMBKhBgAAGAlQgwAALASIQYAAFiJEAMAAKxEiAEAAFYixAAAACsRYgAAgJUIMQAAwEqEGAAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAldoFugMAgF/46dvmXc5X64iKlxypTe8D0ECEGABoKaLipdAo6Z1xjV9HaNTJ9QRq29kbCTJoNoQYAGgpHKknQ0D5ocavo7FHQ5q67Z++PRmAyg8RYtBsCDEA0JI4UgMXAgK5baARGNgLAACsRIgBAABWIsQAAAArEWIAAICVCDEAAMBKhBgAAGAlQgwAALASIQYAAFiJEAMAAKxEiAEAAFbyOsSsW7dO1113nVJSUhQUFKT33nvPY74xRtOnT1dycrIiIyOVmZmpnTt3erQpLi7WmDFjFBMTI4fDobFjx+rIkSNNKgQAALQtXoeYo0eP6uKLL1ZeXl6d85944gk988wzmj9/vjZs2KD27dsrKytLx48fd7cZM2aMtm7dqhUrVmjx4sVat26dxo8f3/gqAABAm+P1DSBHjBihESNG1DnPGKO5c+fqoYce0vXXXy9J+tOf/qSkpCS99957uuWWW7R9+3YtXbpUn332mS655BJJ0rPPPqtrrrlGs2fPVkpKShPKAQAAbYVPx8Ts3r1bLpdLmZmZ7mmxsbFKT09XQUGBJKmgoEAOh8MdYCQpMzNTwcHB2rBhQ53rraioUFlZmccDAAC0bT4NMS6XS5KUlJTkMT0pKck9z+VyKTEx0WN+u3btFBcX527zS7m5uYqNjXU/UlO5VTwAAG2dFWcnTZs2TaWlpe7Hvn37At0lAAAQYD4NMU6nU5JUVFTkMb2oqMg9z+l06uDBgx7zT5w4oeLiYnebXwoPD1dMTIzHAwAAtG0+DTHnnnuunE6nVq1a5Z5WVlamDRs2KCMjQ5KUkZGhkpISbdq0yd1m9erVqqmpUXp6ui+7AwAAWjGvz046cuSICgsL3c93796tzZs3Ky4uTl27dtXEiRP1+OOPq0ePHjr33HP18MMPKyUlRTfccIMkqVevXho+fLjGjRun+fPnq6qqShMmTNAtt9zCmUkAAKDBvA4xn3/+uX71q1+5n+fk5EiS7rzzTuXn52vKlCk6evSoxo8fr5KSEl1++eVaunSpIiIi3Mu8+uqrmjBhgoYNG6bg4GCNGjVKzzzzjA/KAQAAbYXXIWbo0KEyxtQ7PygoSDNnztTMmTPrbRMXF6eFCxd6u2kAAAA3K85OAgAA+CVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAlQgxAADASoQYAABgJUIMAACwEiEGAABYiRADAACsRIgBAABWIsQAAAArEWIAAICVCDEAAMBKhBgAAGAlQgwAALASIQYAAFiJEAMAAKxEiAEAAFYixAAAACsRYgAAgJUIMQAAwEqEGAAAYCVCDAAAsFK7QHcAAFqckn1S+SHvl/vpW9/3BUC9CDEAcKqSfVLeIKmqvHHLh0ZJUfG+7ROAOhFiAOBU5YdOBpibXpA6XeD98lHxkiPV9/0CcBpCDADUpdMFUkpaoHsB4AwY2AsAAKxEiAEAAFYixAAAACsRYgAAgJUIMQAAwEqEGAAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAlQgxAADASoQYAABgJUIMAACwEiEGAABYiRADAACsRIgBAABWIsQAAAArEWIAAICVCDEAAMBKhBgAAGAlQgwAALASIQYAAFiJEAMAAKxEiAEAAFbyeYh59NFHFRQU5PHo2bOne/7x48eVnZ2t+Ph4RUdHa9SoUSoqKvJ1NwAAQCvXzh8r7dOnj1auXPn/N9Lu/2/mt7/9rZYsWaI333xTsbGxmjBhgm666SZ98skn/ugKAAvtLzmmn49WBmTbET8dUXdJhT8e0XFT6pdtdGwfps6OSL+sG2hL/BJi2rVrJ6fTedr00tJSvfTSS1q4cKGuuuoqSdLLL7+sXr16af369br00kv90R0AFtlfckyZT63VsarqgGy/T9BuLQmXHnh9s7b6KcREhoZo5aQhBBmgifwSYnbu3KmUlBRFREQoIyNDubm56tq1qzZt2qSqqiplZma62/bs2VNdu3ZVQUFBvSGmoqJCFRUV7udlZWX+6DaAFuDno5U6VlWtuTenqXtidLNvP+KnWOld6Y+3pOl4p74+X3/hwSOauGizfj5aSYgBmsjnISY9PV35+fm68MILdeDAAc2YMUNXXHGFvv76a7lcLoWFhcnhcHgsk5SUJJfLVe86c3NzNWPGDF93FUAL1j0xWhd1jm3+DQedDE7dE6KllABsH0CD+TzEjBgxwv3vfv36KT09Xd26ddMbb7yhyMjG/V/HtGnTlJOT435eVlam1NTUJvcVAADYy++nWDscDl1wwQUqLCyU0+lUZWWlSkpKPNoUFRXVOYamVnh4uGJiYjweAACgbfN7iDly5Ih27dql5ORkDRgwQKGhoVq1apV7/o4dO7R3715lZGT4uysAAKAV8fnPSf/xH/+h6667Tt26ddMPP/ygRx55RCEhIRo9erRiY2M1duxY5eTkKC4uTjExMfrNb36jjIwMzkwCAABe8XmI+f777zV69GgdOnRICQkJuvzyy7V+/XolJCRIkubMmaPg4GCNGjVKFRUVysrK0rx583zdDQAA0Mr5PMS8/vrrZ5wfERGhvLw85eXl+XrTAACgDeHeSQAAwEqEGAAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAK/n8OjEAgLMrPHgk0F3wqYifjqi7pMIfj+i4KVXH9mHq7GjcTX+BhiLEAEAz6tg+TJGhIZq4aHOgu+JTfYJ2a0m49MDrm7XVlCoyNEQrJw0hyMCvCDEA0Iw6OyK1ctIQ/Xy0MtBd8amIn2Kld6X5w6NVWBOm2ct36Nh30VJ5dMNXEhUvOVL910m0OoQYAGhmnR2Rre8IRftuUmiUUtc8oFRJvwqX9K6X6wiNkrI3EmTQYIQYAEDTOVJPBpDyQyr88YgeeH2z/nhLmronNPBIzE/fSu+Mk8oPEWLQYIQYAIBvOFIlR6qOm1JtNaU63qmvlBIb6F6hFeMUawAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKxFiAACAlQgxAADASlzsDgDgF97cqfuXd8Fu6bhLd8tAiAEA+FRj7tT9y7tgt3TcpbtlIMQAFtpfcqzV3QW5ljf/946WqTF36q69C/Yfb0k7ebuCFqzw4BFNXLRZPx+tJMQEGCEGsMz+kmPKfGqtjlVVB7orfhMZGqKO7cMC3Q00gdd36g46eaPI7gnR3G8JDUaIASzz89FKHauq1tyb09Q9sYF3CLYM4w0ANAQhBrBU98RoXdSZ/2MF0HZxijUAALASIQYAAFiJEAMAAKxEiAEAAFYixAAAACtxdhJandZ8ITiJi8EBQC1CDFqVtnAhOImLwQGARIhBK9MWLgQncTE4AJAIMWiluBAcALR+DOwFAABW4khMG9SaB74y6BUA2g5CTBvTFga+MugVANoGQkwb0xYGvjLoFQDaBkJMG8XAVwCA7RjYCwAArESIAQAAVuLnpDpw9g4A4Gxa++epDeMLCTG/wNk7AIAz6dg+TJGhIZq4aHOgu+JXkaEhWjlpSIsOMoSYX+DsHQDAmXR2RGrlpCGt9oi9dPIo08RFm/Xz0coW/X1BiKkHZ+8AAOrT2RHZor/c2woG9gIAACsRYgAAgJUIMQAAwEqEGAAAYCVCDAAAsBIhBgAAWIkQAwAArESIAQAAViLEAAAAKwU0xOTl5emcc85RRESE0tPTtXHjxkB2BwAAWCRgIWbRokXKycnRI488oi+++EIXX3yxsrKydPDgwUB1CQAAWCRgIebpp5/WuHHjdPfdd6t3796aP3++oqKitGDBgkB1CQAAWCQgN4CsrKzUpk2bNG3aNPe04OBgZWZmqqCg4LT2FRUVqqiocD8vLS2VJJWVlfm8b0cOl6mmolxHDpeprCzI5+sH0MIdPiJVmJP/9cNnDOrB696i+Ou7sPZ72xjjk/UFJMT89NNPqq6uVlJSksf0pKQkffPNN6e1z83N1YwZM06bnpqa6rc+Zsz126oB2OD3VwS6B20Tr3uL4q/vwsOHDys2NrbJ6wlIiPHWtGnTlJOT435eU1Oj4uJixcfHKyjozAmxrKxMqamp2rdvn2JiYvzd1YCi1tanrdQpUWtr1FbqlNpOrU2t0xijw4cPKyUlxSf9CUiI6dSpk0JCQlRUVOQxvaioSE6n87T24eHhCg8P95jmcDi82mZMTEyr/sM6FbW2Pm2lTolaW6O2UqfUdmptSp2+OAJTKyADe8PCwjRgwACtWrXKPa2mpkarVq1SRkZGILoEAAAsE7Cfk3JycnTnnXfqkksu0aBBgzR37lwdPXpUd999d6C6BAAALBKwEHPzzTfrxx9/1PTp0+VyuZSWlqalS5eeNti3qcLDw/XII4+c9nNUa0StrU9bqVOi1taordQptZ1aW1qdQcZX5zkBAAA0I+6dBAAArESIAQAAViLEAAAAKxFiAACAlawLMcXFxRozZoxiYmLkcDg0duxYHTly5IzLHD9+XNnZ2YqPj1d0dLRGjRrlcaG9/Px8BQUF1fmovav2Rx99VOd8l8vV6mqtrfef/umfFB4eru7duys/P99fZUryT6218vPz1a9fP0VERCgxMVHZ2dnueXv27KnztVi/fr3Pa5QCV6ckffXVV7riiisUERGh1NRUPfHEEz6t7Zf8VWtd++v11193z2/u92qg6pRa1/tUkg4dOqQuXbooKChIJSUl7umt5fP3VPXVKjXvfvVHnYcOHdLw4cOVkpKi8PBwpaamasKECR73O/TZPjWWGT58uLn44ovN+vXrzd///nfTvXt3M3r06DMuc++995rU1FSzatUq8/nnn5tLL73UXHbZZe755eXl5sCBAx6PrKwsM2TIEHebNWvWGElmx44dHu2qq6v9VWrAav3HP/5hoqKiTE5Ojtm2bZt59tlnTUhIiFm6dKm/SvVLrcYY89RTT5mUlBTz6quvmsLCQvPll1+a999/3z1/9+7dRpJZuXKlx2tSWVnZquosLS01SUlJZsyYMebrr782r732momMjDT//d//7Zc6jfFfrZLMyy+/7LG/jh075p7f3O/VQNXZmt6nta6//nozYsQII8n8/PPP7umt5fP3VPXV2tz71R91FhcXm3nz5pnPPvvM7Nmzx6xcudJceOGFHuv11T61KsRs27bNSDKfffaZe9rf/vY3ExQUZPbv31/nMiUlJSY0NNS8+eab7mnbt283kkxBQUGdyxw8eNCEhoaaP/3pT+5ptS/4qX9s/hTIWqdMmWL69Onj0e7mm282WVlZTSmpXv6qtbi42ERGRpqVK1fWu+3aEPO///u/vinmDAJZ57x580zHjh1NRUWFe9rUqVPNhRde2NSy6uTPv19J5t1336132835Xg1kna3lfVpr3rx5ZsiQIWbVqlX1hpjW8vl7plqbc7821/eMMcb88Y9/NF26dHE/99U+tSrEvPTSS8bhcHhMq6qqMiEhIeadd96pc5m6/kiMMaZr167m6aefrnOZ2bNnm9jYWFNeXu6eVvuCd+vWzTidTpOZmWk+/vjjphV0BoGs9YorrjAPPPCAR7sFCxaYmJgY7wtpAH/VumjRIhMeHm5eeeUV07NnT9O5c2fzr//6r2bv3r3u9rUhJjU11SQkJJjBgwd7HMHwpUDWefvtt5vrr7/eYx2rV682kkxxcXHTi/sFf/79SjIpKSkmPj7eDBw40Lz00kumpqbGPb8536uBrLO1vE+NMWbr1q3G6XSa7777rs4vt9b0+Xu2WptzvzbX98z+/fvNkCFDzJgxY9zTfLVPrRoT43K5lJiY6DGtXbt2iouLq/d3NJfLpbCwsNNuGJmUlFTvMi+99JJuvfVWRUZGuqclJydr/vz5evvtt/X2228rNTVVQ4cO1RdffNG0ouoRyFpdLtdpV05OSkpSWVmZjh071ohqzsxftf7jH/9QTU2Nfve732nu3Ll66623VFxcrH/+539WZWWlJCk6OlpPPfWU3nzzTS1ZskSXX365brjhBn3wwQetqs769mntPF/z59/vzJkz9cYbb2jFihUaNWqU7rvvPj377LPu+c35Xg1kna3lfVpRUaHRo0frySefVNeuXetcT2v5/G1Irc25X/39PTN69GhFRUWpc+fOiomJ0Ysvvuie56t92iJCzH/+53/WO9i09vHNN980S18KCgq0fft2jR071mP6hRdeqH//93/XgAEDdNlll2nBggW67LLLNGfOHK/Wb0OtvhLoWmtqalRVVaVnnnlGWVlZuvTSS/Xaa69p586dWrNmjaSTd1TPyclRenq6Bg4cqN///ve67bbb9OSTTzZ4OzbU6SuBrlWSHn74YQ0ePFj9+/fX1KlTNWXKFI/95Yv3qg11+kqga502bZp69eql2267rd42reXztyG1+kKg66w1Z84cffHFF3r//fe1a9cu5eTkuOf5ap8G7N5Jp5o0aZLuuuuuM7Y577zz5HQ6Pc6gkaQTJ06ouLhYTqezzuWcTqcqKytVUlLikRyLiorqXObFF19UWlqaBgwYcNZ+Dxo0SB9//PFZ253KhlqdTudpI+qLiooUExPjccTmbAJda3JysiSpd+/e7vkJCQnq1KmT9u7dW2+f0tPTtWLFijP2+1Q21FnfPq2d11CBrrUu6enpeuyxx1RRUVHv/Vy8fa/aUGdreZ+uXr1aW7Zs0VtvvSVJMv93J5xOnTrpv/7rvzRjxow6123j529DavXFfg10nae2dTqd6tmzp+Li4nTFFVfo4Ycfdn9m/VJj9qlVY2JqByF9/vnn7mnLli1r0CCkt956yz3tm2++qXMQ0uHDh010dLR59tlnG9SfzMxMc+ONNzaikrMLZK1TpkwxF110kce00aNH+33AoK9r3bFjh/vMo1qHDh0ywcHBZtmyZfX259/+7d9M//79m1rWaQJZZ+3A3lPPupo2bZrfB/b66+/3VI8//rjp2LHjGfvjr/dqIOtsLe/TwsJCs2XLFvdjwYIFRpL59NNPTVFRUb39sfHztyG1Nud+bc6/37Vr1xpJZvfu3fW2acw+tSrEGHPydLD+/fubDRs2mI8//tj06NHD47St77//3lx44YVmw4YN7mn33nuv6dq1q1m9erX5/PPPTUZGhsnIyDht3S+++KKJiIioc7T0nDlzzHvvvWd27txptmzZYh544AETHBx8xjNCmipQtdae4jd58mSzfft2k5eX1yynbvqj1uuvv9706dPHfPLJJ2bLli3m2muvNb1793Z/mefn55uFCxea7du3m+3bt5tZs2aZ4OBgs2DBglZVZ0lJiUlKSjK33367+frrr83rr79uoqKi/H6Kta9r/eCDD8wLL7xgtmzZYnbu3GnmzZtnoqKizPTp091tmvu9Gqg6W9P79FR1DXZtbZ+/teqqtbn3qz/qXLJkiVmwYIHZsmWL2b17t1m8eLHp1auXGTx4sLuNr/apdSHm0KFDZvTo0SY6OtrExMSYu+++2xw+fNg9v/ZskzVr1rinHTt2zNx3332mY8eOJioqytx4443mwIEDp607IyPD3HrrrXVu9w9/+IM5//zzTUREhImLizNDhw41q1ev9nl9pwpUrcacfHOlpaWZsLAwc95555mXX37Zl6Wdxl+1lpaWmnvuucc4HA4TFxdnbrzxRo+zdvLz802vXr1MVFSUiYmJMYMGDfI4dbC11GmMMV9++aW5/PLLTXh4uOncubP5/e9/77c6/VXr3/72N5OWlmaio6NN+/btzcUXX2zmz5/vcW2J5n6vBqpOY1rP+/RUdX2xt7bP31r1nWbcnPvVH3WuXr3aZGRkmNjYWBMREWF69Ohhpk6d6pd9GmTM//0oBwAAYJEWcXYSAACAtwgxAADASoQYAABgJUIMAACwEiEGAABYiRADAACsRIgBAABWIsQAAAArEWIAAICVCDEAAMBKhBgAAGAlQgwAALDS/wM5/cZN/eVHNAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.hist(density135_csiborg1, bins=\"auto\", histtype=\"step\", label=\"CB1\", density=True)\n",
    "plt.hist(density135_csiborg2, bins=\"auto\", histtype=\"step\", label=\"CB2\", density=True)\n",
    "\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.        , 0.        , 0.        ],\n",
       "       [0.11111111, 0.11111111, 0.11111111],\n",
       "       [0.22222222, 0.22222222, 0.22222222],\n",
       "       [0.33333333, 0.33333333, 0.33333333],\n",
       "       [0.44444444, 0.44444444, 0.44444444],\n",
       "       [0.55555556, 0.55555556, 0.55555556],\n",
       "       [0.66666667, 0.66666667, 0.66666667],\n",
       "       [0.77777778, 0.77777778, 0.77777778],\n",
       "       [0.88888889, 0.88888889, 0.88888889],\n",
       "       [1.        , 1.        , 1.        ]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = numpy.linspace(0, 1, 10)\n",
    "numpy.vstack([x, x, x]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_borg2.npz\")\n",
    "V = d[\"enclosed_volume\"]\n",
    "dist = d[\"distances\"]\n",
    "rho_mean = 0.3111 * 277.53662724583074\n",
    "density = d[\"enclosed_mass\"] / V / rho_mean - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bulk flow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_bulkflow_amplitude(cumulative_velocity, subtract_observer):\n",
    "    if isinstance(subtract_observer, bool):\n",
    "        if subtract_observer:\n",
    "            subtract_observer = 0\n",
    "        else:\n",
    "            return numpy.linalg.norm(cumulative_velocity, axis=-1)\n",
    "\n",
    "    if not isinstance(subtract_observer, int):\n",
    "        raise TypeError(\"Incorrect type for `subtract_observer`.\")\n",
    "\n",
    "    for i in range(len(cumulative_velocity)):\n",
    "        for j in range(3):\n",
    "            cumulative_velocity[i, :, j] -= cumulative_velocity[i, subtract_observer, j]\n",
    "\n",
    "    return numpy.linalg.norm(cumulative_velocity, axis=-1)\n",
    "\n",
    "\n",
    "def plot_bulkflow_amplitude(subtract_observer=False, to_save=True):\n",
    "    with plt.style.context(\"science\"):\n",
    "    # if True:\n",
    "        plt.figure()\n",
    "\n",
    "        # CSiBORG2 main\n",
    "        d = numpy.load(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_shells/enclosed_mass_csiborg2_main.npz\")\n",
    "        dist = d[\"distances\"]\n",
    "        cumulative_velocity = d[\"cumulative_velocity\"]\n",
    "        cumulative_velocity_amplitude = process_bulkflow_amplitude(cumulative_velocity, subtract_observer)\n",
    "\n",
    "        for i in range(len(cumulative_velocity_amplitude)):\n",
    "            plt.plot(dist, cumulative_velocity_amplitude[i], c=\"C0\", alpha=0.25, ls=\"dashed\")\n",
    "        plt.plot(dist, numpy.mean(cumulative_velocity_amplitude, axis=0), c=\"C0\", label=\"CSiBORG2\")\n",
    "\n",
    "\n",
    "        plt.axvline(135, c=\"k\", ls=\"--\")\n",
    "        plt.xlabel(r\"$r ~ [\\mathrm{Mpc} / h]$\")\n",
    "        plt.ylabel(r\"$\\langle U \\rangle ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "        plt.legend()\n",
    "\n",
    "        if to_save:\n",
    "            fout = join(plt_utils.fout, f\"enclosed_flow.png\")\n",
    "            print(f\"Saving to `{fout}`.\")\n",
    "            plt.savefig(fout, dpi=plt_utils.dpi, bbox_inches=\"tight\")\n",
    "\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving to `../plots/enclosed_flow.png`.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUYAAAD4CAYAAACQYE9BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3DUlEQVR4nO3de3xcdZ34/9eZWzK5TpLe09J20pZCKaVJqqyKIp2CixKUJrBZXdddaAIr+9vFS0Pd/e6qX4UmILreIIGv7n5VYptQBV1WybSICFRyaQttoZdMb6TXZHJynfuc3x9h5pvLNE3aJGeSeT8fDx4685kz5z2nM+98Pp/zuSiapmkIIYSIMugdgBBCxBtJjEIIMYwkRiGEGMakdwATZdWqVeTl5V3ydW1tbeTm5l7WOS73WD3OeSXHTsU5vV4vJ06cYPHixSQnJ0/ZeSfquOl4rMQbW2trKwcOHBj6pDZD3HHHHRP6uok8Vo9zXsmxU3HO5uZmDdCam5un9LwTddx0PFbiHftrpSk9DqWlpVN6nJ7H6nFOPY7V4xpd6Xmn2/dQj3NecbyXnZLjzFTUGBPFVFyjQ4cOaR/72Me0Q4cOTfq5Jot8l8Ym3q9TrPhmTB/jWOlVQ5hOpuIarVixgj/84Q+Tfp7JJN+lsZmO10nRtJkxwLuoqIgXXnhB7zDEGIXDYQKBAGazGYNBenSEfmLlDvlGCl3s3buX5ORk9u7dq3coQowwYxJjW1sbRUVF1NbW6h2KEGIaqK2tpaioiLa2thFlM6aPMTc3V5rSQogxKy0tpbS0lKKiohFlM6bGOB4ef5DjF3oJhcN6hyKEiEMJmRgDIY1T7X14/CG9QxFCxKEZ05QejxSLEQCPP0RaslnnaBLTddddx6lTp5gzZ47eoQgxQkImRpPRgNlkoN8X1DuUhGWxWFi4cKHeYQgRU0I2pQGsZqM0pXXkcrkoKSnB5XLpHYoQIyRkjRFgTmYyM2No+/Skqir19fVs2bJF71BmhKqqKgBsNhvZ2dkAuN1uHA4HdrudlpYWnE4nNpsNgOzsbFwuF5s3bwagoKCALVu2UFxcTEtLC9XV1dTU1LB58+boqlXNzc0UFBRQVlYW8/x2uz16XmDI64a/Z05ODh0dHaiqSnl5Ofn5+eP6PJH3dDqdADQ2NlJZWRktu2JTPjFxksT7fEwxVKzVdcTlsdvtWl1d3ZDnWltbNZvNFr2+DodjSHl1dbW2efPm6OOGhgats7Mz+rizs1OLlR6AIefq7OzU7Ha71traOuL9i4uLhzwXec/B54l8DwYfP5bP09nZqVVWVkbL6+rqNLvdPiLesZDVdd7X3uNl5/7TtHd5ZMiOmNZKSkpwOBwUFxcPed5ut0drbC6Xa0SXRVlZGTk5OdHHDocjWpscTX5+Po2NjdHHmzZtori4eERNraysDJfLRU1NzSXfDwZqf5f6PINbF01NTVRUVAyJP9bnvFwJ2ZTu94Z467iKxxviw9fMJSvVondIIg71+4IcPtM95eddMT+DlKSx/TTr6+tpaGiIWRZJJDabDZvNFm0u5+fnY7fbo83olpYWNm3aRHl5ecxmckRLSwvZ2dlDElR9fT11dXUxX19YWEhdXd2o71lfX09xcXE0EY72eQa/j8PhoLm5Ofq4qakJYMKa0gmZGBUF2tx9zMlMwuMLSmLUwYIFC3jkkUdYsGCB3qFc1OEz3dz0b7+b8vO++s1PcMOS7Eu+LlI7ivTBDTe4Btjc3ExFRQUVFRW4XC7y8/PZuXMnNpuN/Px87rnnHlRVHfEe9fX12Gw2XC4XDQ0NlJeXR983cv6LJaO8vDy2b98+4vlIX2ddXV30f8f7eYAh/ZKVlZVUV1fHPO5yJGRizEgxY1AU+n0hPAG5M62HefPmxf2NlxXzM3j1m5/Q5bxjMfxmx6VUVlZSWVmJqqps2rSJTZs2XbS2FzG4SVtWVkZBQQENDQ1UV1dHzx9JtMM1NjbGTJqRZrvD4WDDhg1UVVWxefPmi34ep9NJQ0MDLS0t5OfnU1lZOaS8pqaGkpKSUWum43ZZvZVxKD8/X7vjjju0Z5999pKvDYXC2r/Wtmjf/+8D2tsn3FMQnRius7NTe/7554d0xIvxKy4u1srKymKWRW5gNDc3j7jJFblpElFZWTnkZsbFbr5UVlZqNpst+tjhcFz0/Ha7Xauurh7xnoP/zTdv3jwkjtE+DzFu1jU0NIy4UTNWzz77rHbHHXdo+fn5I8pmzM2XyCISY1kU02BQsM9LJystaQoiE7G4XC7uvPNOGcd4hZ5++mmcTif19fUjygbXvB599NEhZU1NTTgcjnGfr7GxkcLCwujjuro6nE5n9OZJRGT4zqVqcXl5edHvgKqq0c8z/KZNrO9J5JyRWm1NTU3M7oCLKS0t5YUXXoi5aVZCNqUBblo5h67+ANddlaV3KEJcNpvNRmtrK1VVVVRVVQ0Zpzg4Ka1bty6abFRVpaOjI9on19LSwrZt24CBJKOqarSsoqKCvLw8VFWltbWV7Oxsnn766RHnr6iooKmpKTo+0mazDbmJoqpq9Pw1NTWUlZVhs9koKyujoaEhGntZWVn081RUVETvnLe2tlJZWTmkf7OgoGDEtZio5nRCruDt9Yc4cqYLk0lh+bxMjAYFRVEmOUIxWEtLCwUFBTQ3N8fsnxJiqsTKHQlZY3T3+djx5gmy05Jpc3tZsziL2RnJeoclhIgTM6aPcTxSLEaCQY2TF3rp6PHKnGkdJCcnc+2115KcLH+QRPxJyBpjutVMSpKJM11e+nxB/EFJjFPt2muv5cCBA3qHIURMCVljNBoMZKZa0LQwHl+IQFCmBQoh/p+ETIw9ngBef5D09xepDYRmxP2naWXv3r1kZGTILoEiLiVkU9pkVFgxP5O/WDGHPl+AJXNS9Q4p4YTDYXp6egjLIh4iDiVkjdGggC8YxmxUyM1OJdlk1DskIUQcScgaY1iD810eXjl4lusXZ+Hxh1i7NFvGMgohgBlUY2xra6OoqIja2tpLvtZqMWFLtXCqo4+TF3o5dLqLQEiadEIkktraWoqKimhraxtRNmNqjJG50mO1ckEmthQL7b1+znV58QfDWKRJPWVWrlxJc3MzK1eu1DsUkaBKS0spLS2lqKhoRNmMSYzjEQqHudDjAwVCoRD9vqDcmZ5iKSkpMhVQxK0Z05QeD6PBwFWzUggGw/R6gwRCYQIyyHtKnTx5ki9+8YucPHlS71CEGCEhE+PuIxd44rcHuWFpNrMzrGSmWEiSZvSUam9v58c//jHt7e16hyLECAmZGAPBMLWvHSfDaibZYkBRGPMeG0KImS8hs8HapVlkpVr41ZunQNMIhTQ6er3Ms6XoHZoQIg4kZI0x2WzkhiVZ7DvRSY8vwMmOPt4+2al3WEKIOJGQiTEUhkU5qXR7AqRYjJxq72PfCUmMU2nOnDk89NBDzJkzR+9QhBhhyprS9fX1uN3u6NaIkX0aXC4X9fX12O12XC5XdMnzS5VdCX8whDXJSGqSiQvdfpJNBs52eq74fcXYLVy4kCeeeELvMISIaUoSY2STms2bN+NyudiwYUM0MZaUlEQ3zlZVlZKSkuheEaOVXQmLycjKBZnkzU3j8JlulsxORe314QuESDLL3emp0Nvby9tvv83q1atJS0vTOxwhhpiSpnRFRQWbN28GBvbCjSS74TuL2Ww2mpqaUFV11LIrZTIqFNpzWLMkC28gTFqSGU2B9m7vFb+3GJvDhw/zoQ99iMOHD+sdihAjTHpijCQ4p9OJ0+mkoqIiuq1jZFexwbKzs2lqahq17Er1eoN09Pi5cflsstMsmIwGQiGN9h7fFb+3EGL6m/SmdKSWZ7fbsdvtFBYWUlBQQGtra8zan81mQ1XVUctiiSwiERGZBxlLWrIJi9nAVbNSWTI7jdOd/Vw1OwW13385H1EIMY3U1tYOWWxGl0Uk7HY7NpsNu90ODCQ3t9tNS0tLzBspqqpis9lGLYtlPItIGA0GjAaFrv4gmqZxutNDW0c/rjM9fOzaeWP9aEKIaWh4pSnWIhKT3pS22+0janmRx4WFhdFmdYTb7Y7WLC9WNhHCmsYZtZ+sFAvBsMbZLi9qvxd/QOZMTwWTycSsWbMwmRJyjoGIc1NSY8zPz8flckWH3USeA4YkzcFN7kuVXQlN0/D4ghgUhXlZA9t39vuCHD3fS483SI7cmZ50119/PRcuXNA7DCFimpI/13V1dVRXV5OXl0dzc/OQITd1dXVUVVVht9tpbGykrq5uTGVXQu3z4zrfS968dLy+ENveOInHH+R8p5fOPh856UkTch4hxPQ0JYnRbrdTWVkZsyw/Pz9ae4yMbRxL2ZUKhjR8/iBrl2aTk55EnzdAlydAR4+PZdLNOOkOHDjAnXfeyfPPP8+qVav0DkeIIRJySmBGipnl89Pp8gTx+EN8eMVs/MEw7T0eOnpkBsxU8Pl8tLa24vPJECkRfxKy59toMDA308rhM914/EFy0i0EQhpmk5GOnoDe4QkhdJaQNUZN0zh0upv2bh9oChfen/ESDGpYLXLjRYhEl5A1RkVRWDE/nfRkE0kWI28eNWM2Kpzv9tLjlaadEIkuIWuMAAaDgdycFMxGA0vnppFhNdPrDeJ866zeoSWEZcuW8bvf/Y5ly5bpHYoQIyRkYtQ0jV5vgPfc/eRmWTEqClkpFkJhjdMd/fgCQb1DnPEyMjK47bbbyMjI0DsUIUaYMYkxMld68BzIi/H4Q4TCGmnJZvwhjZz0JBbOSkEDznZ76eqXGzCT7cyZM3z961/nzJkzeociElRtbS1FRUUx50rPmMQYmSt9sYUjBlMUMBoUPN4A+467WTY/g1uumw9AIKjh7pXFJCbbmTNn+MY3viGJUeimtLSUF154gdzc3BFlMyYxjofVYmJ+VgrvdXpocnVwocuLpg0kTBQNtyw/JkRCS8i70gC52VYW5qSgaRrt3R7ePHqBtCQTYQ3MRk3v8IQQOkrIGiPAqY5+0pLMbPzgVayYnwkoWC1Ger0BDp3p0Ts8IYSOEjIx+oMhznR6WJiTQlZaMmvtOczPsmI1G/EGwuw+Iqu+TLasrCw++9nPkpWVpXcoQoyQkE1ps9HA1QsyyEq1AHCm08OFbi9GI2gaHJA9pifd0qVL+fnPf653GELElJA1RkVRmJ2RjMk48PFXLMjAajFGdwg81eGhR7Y5mFRer5ejR4/i9coGZCL+JGRiHM5oMHDXBxezYv7AYONAKEx7r/xgJ9PBgwdZvnw5Bw8e1DsUIUaQxPi+dctmcd3CLExGhbAGCoreIQkhdCKJ8X3vtfdzprOfVIuRUFhD7ZWxjEIkKkmM79PQ8IfDpCSZ8AfDvLT/LP6gbIwlRCKSxPi+WelJ5GYNDPj2BcMcO9czsF6jECLhzJjEOJ5FJGLJTE3CPi8z+njfCTdt7r6JCk8Mk5+fj6Zp0T19hJhqsojEGGRYzXxoxSxys60AnHF72H9KncAIhRDxZLRFJGIO8D5+/PiQPZ3HwmazsWTJksuJL24smZPOlz55LX/zo9fp8wU5eaGHcFjDYJA71BPt0KFDfOELX+A///M/ufrqq/UOR4ghYibGrVu3cuutt6JpY19MYfv27Wzbtm3CAtPLNYuyybCa6fYESLYYJSlOkr6+Pnbv3k1fn3RXiPgTMzFu2LCBu+66a1xv5Ha7JyQgvSmKQk66BbU/wLtt3fgCoeiMGCFEYojZx7hx48Yhj7ds2cKOHTvo7u7m1ltv5bbbbmPHjh1DXrNp06bJi3IKzc9MJtNqwaDAnuOdvHxAFlIVItGM6eZLYWEhd911F4888ggFBQX8/ve/p6OjY7Jj00VKsonMVAtWi5FT7f0cauvC45c9YIRIJGNKjJGloXbu3Mk999wDgN1un7yodKQoCh++ejYWowF/KMwLTac43enRO6wZZ8mSJfzsZz+b9jfsxMw0pmXHmpub0TSNzs5ObrjhBo4dO8axY8cmOzbd3POhJfzx4Dn++O4F3j6pcuCUm8WzUqOr8Ygrl52dzec+9zm9wxAipjH90svKymhpaaGhoYGuri6qq6vHPZxnOlkyJ4PvfH4dVrMBbyDMG+9ewBuQ6YET6cKFC/zoRz/iwgVZFFjEnzElxszMTL761a+ydOlSMjMz2bp1K1/5ylcmOzZdrVxo48MrZqEBzv2n6Zb1GSfUqVOnePDBBzl16pTeoQgxwojEuGfPHu6//362bNnC8ePHdQgpfnx45RxSk4wcPtPL1l+/zfHzsheMEIlgRGJcu3YtTz31FA8//DB1dXU88MADPPPMM3rENi5XOlc6FvvcTG5bPR8NeL7xFDvePEE4LDsICjETjDZX+qI3XyLNZxioRT788MN0dXVRUlLCLbfcMnnRXqbIXOmJtGRuGquXZPPakQuc6/Kx/fXjfOyauRTkzZ7Q8wghpl5paSmlpaUUFRWNKBtTH+PatWvZunUrTz75JJ2dnQnT1L4218aK+RkULM3GoMC7bd388MWD9HkDeoc27aWnp3PrrbeSnp6udyhCjDDuXQI3btzIxo0b6erqYvv27TQ3N1NYWMh99903GfHpKtli5Obr5qMo0Hquh0Nnenn5nfM07DvNpz+4WO/wprXly5fz+9//Xu8whIjpsrdPzczMZNOmTWzatImurq6JjCmuZFjNrF+9AJNB4Qs/fg13X4Ddh8/zqcJFMq7xCoRCIfr6+khNTcVolLnoIr5c9i+7u7s7+v8zMzNHeeX0l5JkYv318/mAPQdNg+1/PsGRMzP3j8FU2LdvH5mZmezbt0/vUIQYYcyJsbu7m+PHj0f/q6iomMy44o7FZOLv1y8nO9XMhW4/3/3NQUKhsN5hCSEmwZia0vfffz9OpxObzRZ9bs+ePTz55JOTFVdcunPdVby8/zQ/feU4L+5to6PXx5xMq95hCSEm2JhqjBs2bODo0aM0NTVF/3vqqacu64Tl5eVDphO6XC6qqqqor6+nqqpqzGV6MBgM3HvL1ViMCl2eID94UTaLF2ImuuybLxs2bBj3MU6nk+3bt1NZWRl9rqSkhObmZgBUVaWkpISGhoZLlull9eIsbr5mNi/tP8//efko5RtWsnBWqq4xCSEm1phqjPn5+Tz++OPs2LGDXbt2sWPHjnH3MUZqe9nZ2dHnWlpahrzGZrPR1NSEqqqjlunJYDDwv0sLsFoM9HhDPNNwWNd4pqvVq1dz/vx5Vq9erXcoQowwphpjRUUFqqqO6GMcj+3bt1NWVjbkuaampiGJEgYSZ1NTEy6X66JlDodjXOeeaNcstPHJG+ZT/2YbzzWe4O9vXcFVOVJrHA+z2czs2TKDSMSnMSXG8vJy1q9fP+S55557bswncTqd3H333SOej1X7s9lsqKo6apneFEXh7265mt/uOcPJjn5qXjrE//6rtSiKbJw1Vq2trTz00EN897vfJS8vT+9whBhiTIkx1g9+rF/mSCIbXNuMiPVcpGY6WlkskUUkIiLzICdLvj2Hv1g+i5cPnue/XmnlznWLWLdMakBj1dXVxW9+8xu+/vWv6x2KSDC1tbVDFpsZ1yISg7W2tlJdXc26desA0DSN7du309jYeMljnU4nbrebmpoaYOBOc01NDQ6Hg8LCQqqrq4e83u12Y7fbyc7OvmhZLJOxiMRo0pLN/OPt19Dc2o7aH+AbdXv5zcMOqTUKEeeGV5piLSIxpsRYXV3NPffcM2Sf6bHuOV1cXDzkcXl5OcXFxdEEN7hprKoqdrt9TGXx4GPXzuOTBQupff0kr757gV++2krpR5fpHZYQ4gqNKTFWVlaO6GMc7w0QVVWjtcbKykoqKiqw2+3U1dVRVVWF3W6nsbGRurq66DGjlcUDi8nIt/4qn2aXm8Nne9nyy718smARGalJeocmhLgS2hg899xzI557+OGHx3LolLnjjjt0O/cLjSe0OffWaul/8wvt0ef26RbHdHL27FntO9/5jnb27Fm9QxEJLlbuGNM4xuFjFh977DGqqqomJVFPRx++eg43XzcXDajZeZgT53v1DinuzZ07ly996UvMnTtX71CEGGFMiXH9+vU888wzHD9+nMLCQhobGy97SuBMlJ2ezL9vvGFggYkeP/+2rWXMfbCJqrOzk7q6Ojo7O/UORYgRxpQYn3rqKUpKSigpKaGqqort27df1pTAmeyahTYeuHUFCvB8cxuvHTitd0hx7dixY9x9990zen9yMX1d9ObLjh07Rjy3fv16WlpaUFWV6upqWYF5EEVRuOuDS/i/r7Ryyu3lcz9+g4NPFJGSbNE7NCHEOF20xrh582YaGxt58803o/8BtLe38+abb9LR0TFlQU4XKxZk8u3SfJJMCh29fu5/erfeIQkhLsNFa4zV1dUjhugMtnPnzkkJaLorWncV+0+pVL1wkN+2tLHvWAdrluboHZYQYhxi1hh37do1alIERpTv2rVr4qKaxowGAw/cdjXL56QSCGl89vt/pMfj0zusuGO1Wlm7di1Wqyz0K+KPosW4fXr//ffz8MMPj/lNNE2jqqpK1xW9CwoKyM3NnfQ50mO119XOJ7fuotsb5CNXz+KFilswmy57+UshxASLzJlua2uLrvsaETMx7ty5c9yr2NhstkvWMidTUVHRlM6VHovftZzgr77/GqEwPLBhOVV/s07vkIQQw8TKHTGrMHomuJnktrVXcc8HT/HsGyepdh7hE2vmccv1i/QOKy7s2bOHG2+8kd27d7N27Vq9wxFiCNkYeRIpisK3/7qAaxekEtbgr3/wOu1d/XqHFRc0TcPv98tAeBGXJDFOslmZVn7yDx9ldpqZPl+I4ideIRSWbVeFiGeSGKfAqquy+P4XCjEZFVqOdfLY82/rHZIQYhSSGKfIJ9ct4eE7rgVg6/MHcO45pXNEQoiLkcQ4RRRFoeKuNfzl2rmEwnDP9//E7kNnEraP7ZprrmH//v1cc801eocixAiSGKfYTx64idlpZvwhjTuqXqG64d2E7HO0Wq2sWrVKBniLuDTmxPj0009z//33s2XLFvbu3TuJIc1sqckWXvq3W8lJM+MNhHn4F3v40n/9mc7exJodc+LECe677z5OnDihdyhCjBBzgHfEnj172LZtG4qicM8993DDDTcAA0mypaWFgoIC7r77bjIyMqYq3ouKxwHeo3Gd6+GBmtd4/YgbgAW2ZP759pXcvHo+FqORrFQL2enJOkc5eSLfn+bmZvLz8/UORySwMQ/wfuyxx3C5XGzYsIGtW7eOKN+0aRMwsKbeI488QldXF+Xl5dHEKS7NPjedX33VwQ9/d4Dvvfgup1Uvm5/di8JeDAokmQ0sm5vGZ2/K4/MfzSMtRZYvE2KqxKwxdnV1kZmZOa43upxjJtJ0qzFGhMMavR4/j73wNj/5QyseX4gwEBrU7Tg73cLWz65l441LMRpmRrew1BhFvIiVOy7alN67d++0qgHG2yISl8PjC6J6/Ki9fvYca+d3e9p45d0LuHv9KMDKBen88p8/in2efn+AJookRqG3cS8iAXDPPfewbdu2KQlwIkzXGuOlBIJhyp96jfo3T6Ex0MT+t7tW8/99cpXeoV2RtrY2fvjDH/Lggw+Sm5urdzgigcXKHRdtlzU3N/P444/zzDPPyF1oHZlNBn7y4E384h8/RE6aCV8gzL9u28d9T/5xWo+BzM3N5dFHH5WkKOLSRRcIrKysZOPGjcDA3emnn34aRVGw2+3ccsstUxagGHDHuiXcdM08tu54iyd3HmXbG++B8hrf+9sPkGadfjdmenp6aG5upqCggPT0dL3DEWKIi9YYI0kRYO3atWzatIn77ruPrKwsHn/88SkJTgxlS0tm6+c/QM2mdRgNsO31k9z+iJOuvuk3BvLIkSN8/OMf58iRI3qHIsQI477FGUmSsXYRFFPjno8s5+nyGzEbFPacUPnEIw30efx6hyXEjHHRxDharTAzMxO32z0pAYmxKb5xKd/9QgEKsP9UN595bBf+QEjvsISYES6aGKurq/nVr35Fd3d3zPLxbn0gJpaiKHz+Y8v51t2rURR446ib+6v/NK1vyAgRLy6aGDVN46tf/So2m43ly5ezZcsWdu3aRXd3N8ePH+fo0aNTGaeIQVEUHrz9Op78uwIMCtS92cY//eQNvcMaE7PZTG5uLmazWe9QhBjhknelu7q6aGhoYPv27VRXV6OqKnl5eTQ0NExlnOIiDAaFz958NUajQvkzTfz0leOcVT1s/3J879uzevVq3nvvPb3DECKmUReRmE5m6gDv8fht43E+96PXCYXBPjuFlqoijMaZMYVQiMky5gHejz/++EX7FmPZuXMnu3bturLoxBX71LolPPdPN2FUwHWhn+u/8jz+QFDvsGJ6++23WbhwIW+/Lds8iPgTMzFu2rSJ6upqHnjggYsOyzl+/DiPPfYYDzzwAIqi6D7ou62tjaKiImpra3WNQ2/r1y7ij1/fgMWocLLDw9rNv6Gr16t3WCMEAgHa2toIBAJ6hyISVG1tLUVFRbS1tY0ou2RTeufOnTQ0NETXZGxtbaWhoYFly5axadMmXVfUGUya0kPtP36B9d/aSb8/zIKsZPZW3YE1KX5udMgiEiJejHk9xsHWr1/P+vUDHfnPPfccdrudp556anIiFBPmuiWzef1bn+DmrzdwutPLDZt/w57KT5KSnKR3aELEvXH1zG/cuJG1a9dOVixiguXNs7Fn6yexpZg53ell6YO/4sS5sfcdC5Go5JblDDfLlsKeytvJsJro94cp/NqL7D/WrndYLF++nJdffpnly5frHYoQI0hiTACzMlM5+MSdzMmw4A2E+eg3G9i1b2SH81RKT0/n5ptvlpV1RFySxJggMlOTOPidT7N6UQaBkMZd332F/2k5qVs8bW1tbNmyJeYdQSH0JokxgSQlmXj1m7fzqRvmEwrDX/3Hn/jeb/YTDk/9GP9z586xdetWzp07N+XnFuJSJDEmGKPRwLMP3UzZejuaBv+r7i0+9WgDwcG7bwmR4C45XGcitLS04HQ6AWhsbKSyshK73Q6Ay+Wivr4eu92Oy+WirKwMm812yTJx+RRF4Tt/eyO5WSn8e/1+Xj3UzvVf/jWvfuN2cjJn7l7WQozVpCdGVVVxOp1s3rwZgPr6ejZs2EBraysAJSUl0R26VFWlpKQkukDFaGXiyn2p6HpuWJxNyfde5ZTby6qvPM9vKz5O4bI5eocmhK4mvSnd1NRERUVF9LHD4cDlcuFyuWhpaRnyWpvNRlNTE6qqjlomJs4taxZy8Ik7WJSVTJ8vxIZv7+Rnf5j87QZycnK49957ycnJmfRzCTFek54YHQ7HkD1bm5qaALDb7TQ1NZGdnT3k9dnZ2TQ1NY1aJibW3Kw09n/303z82tkEQxr/8JNGPv3oS4RCk7ci+OLFi3nmmWdYvHjxpJ1DiMs1JX2Mg+fCVlZWUl1dDcReBdxms6Gq6qhlsUQWkYgoLS2ltLT0iuJOJAaDgRce3kDVr/fx7V8dYOc77eQ9uIPdj3ySeVkpE34+j8eDy+XCbrdjtVon/P2FuJja2tohi83EGjI2JYkxoqamhpKSEsrKygBi3khRVRWbzTZqWSy5ubmyiMQE2PzpNThWzeXWR3fR0Rdg9Zdf4L+/tp4PLJs9oed55513ZBEJoYvhlabBFaqIKRuu43Q6yc7OjiZFgMLCwhGbarndbux2+6hlYnLlL5/HyR9tJDcrGW8wzCe+7eTlt2S1bZE4piQxRm6kFBcXAwM1R1VVyc/PH9I0VlUVu92O3W4ftUxMvhRrEvu/82lutGe9P1PmVepec+kdlhBTYtKb0i6Xi4KCgiHP2Wy2aM2xrq6Oqqoq7HY7jY2N1NXVRV83WpmYfCaTgZf+/RP89ff+wG/3nOHemt0cO9/DVz99PYqi6B2eEJNG9nwRY1L25J+ofWNgbvXffXQp3/37D2I0XH6DY8+ePdx4443s3r1blrITurqshWqFAKh54CMoyis8+3obP/3jMc529fPLhz6O4TKT49q1a/H5fBMcpRATQ+ZKizGrvv9jPHznNRgU+J995/jCj/7EDGlwCDGEJEYxLv+ycS1PbfoABgV+1fge9z/1Gr7A+AeCv/POO+Tn5/POO+9MQpRCXBlJjGLcSj+yjCc+txaDAs++cZJ7f/RHfOPcptXj8bBnzx48Hs8kRSnE5ZPEKC7LvRuu4Vv3rEEBnm85Q/Hju+J2D2shxksSo7hs/3j7Kv7lM9cC8Id32ln/zZcIBiU5iulvxiTGyFzpwXMgxeSr+MwN/Nc/fBCjAntPqDi+2UBIFr0V00BtbS1FRUUx50rLOEYxIZoOn+W2rX/AHwxTsNTG85vXk5l68T2sOzs7cTqdOBwOsrKypjBSIYaKlTtmTI1R6KtwxTz+9PVbsZoNNB9TufFfXkTtvfg4xaysLEpKSiQpirgkiVFMmGuuyuaNb96G1WzgPbeHlQ/9GtdZNeZrz507xxNPPCGbYYm4JIlRTKi83Cz2Vn6KtCQDfb4QhQ//D9v+NHJF8La2Nr785S/L9qkiLkliFBNuwaw0XD/cyIp5aQTCGvfVNHLvk6/KLBkxbUhiFJPCmmSmqfIOtty5EoMC2984hf2Lz/H6u2f1Dk2IS5LEKCaNoih8bWM+f/w3BzlpFtp7/dz2yC5u+tcX6fXIAhIifkliFJNuTd4cjnz/MxR/YCEAe0+qfOo7u1ly/YcxmGQfaxF/JDGKKWE2Gfnpgx+l+dt/yZI5VrS0ubSv+Qcc33uLf61tJDiJOxIKMV6SGMWUWrEoi7ce+zT1/3QjudYAPr+f//ifI8y5bzt/+4NX8Pj8eocohCRGMfUURWG2onKo5gt8/zNzuWpWCoGQxo7GNnLvf46v/fzPsiCF0JUkRqGrNUtnc+CJT1P/zzcxO91MIKTxg5daWXh/PV/7eSNevyRIMfVmTGKURSSmt9vyF+H6UQn/9/4bWWBLxhMI84OXjrDw/no+/x+v0NXn1TtEMcPIIhIi7rS0tFBQUEBzczP5+fkjyl/e9x4P/ayZ1vN9ABgU+My6XJ74wgfJTpM72WLiyCISYtr4+JqF7H38Tv747w4WZCYR1uC5N9u4+p9+zeb/epNAUJY2E5NHEqPQxZo1a+jq6mLNmjWjvm5t3hwO/WAjz3/5o8zLTMIbCPPkzqMs/eJzPPbrtzhxoVemGooJJ4lR6MJoNJKRkYHRaBzT629Zs5AjP9jIb776MeZmJtHlCfDNHfu5/isvcOO//De1f2qlvVv6IcXEkMQodHHkyBFuu+02jhwZufLOaG5encuh732G/1O2jrkZFsIaHHyvm7KaP3Pdl5+n6NEGWlzt+AJBwmFpbovLY9I7AJGYenp6eOmll+jp6Rn3sUajgbs/spy7P7KcQ6fcPPzLPfz5SDs93hAvv3OBl7/+EhaTgsVoYJ4tmY9fN49brp2HNcnErAwrS+akkWQy0OMNEgiFSTYbsZqNJFuMKIoyCZ9WTDeSGMW0dvWibH711fUANB89R9ULB3j13Xb6fEH8wRBHz/Vx9FwrT+9sxaiA1WIkNcmEhkaKxcQCm5U5tmTCGqxalMHV87PwBYO09/ixmAxYjAbcvT6sFhNL56ZhSzGjeoKkJ5tYkJVCWrIJk2kgsZqNBowGBYNBkut0J4lRzBgFy+ay7UtzCYc1Ons97Dvh5rdNp9l95Dxnu7z0+UJ4AyF6fZF52X6Ot/dHj3+heeh4tuHpbXBl0qCA6f1ECJCVaiEr1YLJAItnpfGXa3NZmZuJ2u8j3ZrMnMxkMlIsJJsNJJulZhrvJDGKGcdgUMjJSOGW1Sncsnph9Pl+b4B3Tnfy9oluTrt7udDt5VR7Pz3eAP5gmI5eP519fgLBEIqioAFmo4GUJBOhcBiPP4QvGCYU1tA0Df/7Q4Y0Dfp8Ht5zewDYc6KLXze3oQBGg4LFpJBsNmE2KWRazXz46jncev0C0qxGMlKSyM1OwZaaRJJ5bDeixOSTxCh0sWjRIn74wx+yaNGiKTtnSrKZAvscCuxzYpYPJL8ggaBGKKzhD4ZIT7aQZDHQ6wlyocdDKAyhkMa5Lg++YJAMaxJqr5e3TnbyblsX73X00+0JcL7LSzCsEdIgGArT6fejAee6fBw+28tPX3FhVAaSuNVixGoxMDsjmU+syeW2NfNZMd9GSrKZZIskSz3IzBchJlg4rNHZ56O7389Z1Ut3v5+2zn5cZ7s5+F4377Sp9AfChEJhfMGB/wb/ChXAbFRISTKRm21l+fwMPrxiFh+/Lpe8eemYjDKYZCLFyh0zJjEWFBSQm5tLaWkppaWleocjLsHtdvPiiy9y++23k52drXc4U+qc6uFCtxd3j5cwCu4eH2c6+3jlwFneek9F7Rto2gdDGoN/nAYFbClmrpqVxprFNtYszsI+N5OVuelkpSVjlbvq41JbW0ttbS1tbW00NzcPKZsxiVFqjNPLpeZKJ5JwWENRBpZju9Dt5b2OPt7r6MPd4+Wtkyp7j7t5t62LPn+I0KChmdE+TKNC3rx0PrFmPutXL2BelhUFBVtaEunJZswmqWGOJlbukD5GIXQ2eHjP7IxkZmcks3ZpDuGwRr8/iPL+a94+0cmuA2doPdPF26e6aO/xofb56Q+EefvUwHOP//Zdks0GbKkWbKkWrl6Qzm3X57J6cTbXLbJhNBrwBUJYTAapXY5CEqMQccpgUEhLNkcff2D5bNYtm0W/P0SfN0BXvx+PL8TBtk7qdx/nrRMq3f0BvIEwZ1QvZ1Qv77R183xjGxaTgZw0Czcszeaa3AwyrBaWzcsgNyeFebYUZqVbUBQFk8Eg4zCRxCjEtKIoCqlJJlKTTMzJtAJw/ZJs7v4LO2q/n65+P4dPd/Num8oeVwevHT6Pu3egz/K06uX0ntO8uOc0BmWgGZ6ZYmZJTgqrrspi9VU2TCYjgWCIFIuJdKuFwrxZZKZa6Orz4/EHCYY1TEYFo2JgVoaFFIuZUFjDGwihoWFQFIyGgWFOJqMBw/u1Uk2DSI+p8v4I0Uj3QTg80J8aCmsYDKBpCv5ACKNRQVEUuvv99HgDKIqCQQGPP0Q4rDHPloIvEOBAWxdaGKxJJgwGuGnlvCu+zpIYhS5SU1O58cYbSU1N1TuUGcFgUMhOSyI7LYmlc9K57YZcYKD/sscb4OSFHnYfvsC2149xVvXQ2R+k1zsww6e9x0/TcRUYuBs+MINnYAB7isWI0WBgdkYSszOSSTYb8QaCgEJqkhFfMEy61YLFaMAXCuP1BzGiYDYrmIxGUizGgWFLIY2u/gBJZoUUi5HUJDPtvV4yUywoKJzr8uB/PyGHwnC+x0NOqpWcDAseXwh3r5dks4mwptHe48OgwDxbCtYkI+4eH2aLAYOmEAiGJDGK6evqq6/mjTfe0DuMGc9gUMhMsbB6cQ6rF+dw7/qr6fUFUXt97D5ynlcPnuVkex/nuryc6/bi8Q8MYPf6woS1EJ19AQBOdgzMELIYDSgG0DRtYBB8WHu/yW8iN8tKMKTR7Q3gC4QJhMIYgDSrmcWz00i1GDnf7edUez8oA++RYjExP9vKDYuzOdHez6EzPSjKwOo2/b4eur1JFNqzmZ1h4ejZXvr9ATJTTIBGd7+PebYMcrNSeKdN5XSXB1uqGde5Huxz06/oukliFCKBGAwKGVYzGVYzV81O4+4P2QHw+UPsO+nmbGcfvd4Qh053s+94B2EgK83CmU4PHd0+UKDfF8IfCGM2GwiFod8XxOMP8c7pnuhsoME6+gKcGDT1cjBFgWMX+mhq7cSg8H5zeqCxrWlgMPTx9kkVTQNvIIQ/GCasDTS7ARr2nxvxnrMyDvDUphuv6DpJYhS6kOE68SXJYuQDy2YDswmFw/R4gnT3+zEaFHzBMAff6+TPR9pJthhw9/jZf1LFYjawYn4GZpORUx292FIsWC0mTrv7sJgMzM9Kpat/YJqlxWyg3xui3xfEaBxo9l/o9uILhNE0jSSzgWAIDIaBeedGg4FzXR7Cmkay2YCiGAiHw2SlJpGeYqa9y0sw9H5/pxECQY1gWGNOehL/a+P1V3w94joxulwu6uvrsdvtuFwuysrKsNlseoclxIxmNPy/4T4Ri2enctM18+juD9DnC3Kh24snECQ3KwVfIMSu/Wfo84YIaxrWJDP93iChsMb61fM5erabM51ers210uMN8F5HPxlWE9fkziGsQU+fn6yMJNRePxd6vaRaLJiMCrMykuj3BclOs9DtCdLr8YMGaRYDyTkp9PmCdPb5CIUVer0B+rxBLCYDudkpV3wN4joxlpSUREekq6pKSUkJDQ0NOkclROIxGgxkpljITBlIlitzM6Nl/mCIubYUPP4gvsDAYhtef5BVV9nISk3ijcPnOXBSxRMIYXx/brgtNYk7ChbR6w2w6+0zdHkCpCZbCGkDd7Tzl2aTNzedRlc7be5+LCYjRkUhpGksmpXOpwoXcfCkysHTXaQnm0i3msm0mrjuqqwJ+bxxmxhbWlqGPLbZbDQ1NaGqqtQahYgjFpNx1Fra2iU5rFqUhS8QIhgaGJozOyOZZIuR810e1i2bRbcngDcQIhTWyEwxc92iLNKSTXR7A8zLtJJsMWE1G0i1mlkxLwNrkon5Niu3rV2AxTTxC23EbWJsamoaMYc2OzubpqYmHA6HTlEJIcYr2WIkGSNYzSPK5mRao+MxY/nIyrkXLUtJmrz0FbeJUVXVEc/ZbLaYzwO0tbVRVFQUfSyLScS3a6+9liNHjrBw4cJLv1iICRRZPCKira1txGviNjHGai6P1ozOzc2VRSSmkeTkZJYtW6Z3GCIBDa80Da5QRcTtshuFhYW43e4hz7ndbux2u04RiYl07NgxPve5z3Hs2DG9QxFihLhNjPn5+UOazaqqYrfbrzgxDq5Ci9im4hp1dnbyi1/8gs7Ozkk/12SR79LYTMfrFLeJEaCuro6qqirq6+t59NFHqauru+L3nI7/SFNNrtHYyHUam+l4neI6Mebn57N582aKi4uprKzUvRl9uf/AV/LF0OtYPc6px7F6/Win22dNpO8SxHlijDfT7QuZSF9mSYzxec4roWe8M2Zrg1WrVpGXl3fJ17W1tZGbm3tZ57jcY/U455UcO93ivZJjp1u8V3KsxBtba2srBw4cGPLcjEmMQggxUaQpLYQQw0hiFEKIYSQxCiHEMJIYhRBimLidKy0mT319PcXFxUOeG21RYFkwWCSahLgrLT/sAfX19bjdbsrLyxn+zx7ZZgBGLgo8WtlM1dLSgtPpBKCxsXHIBAP5IzLA6XSiqiput5uGhga2bNkS3aZi2l8jLQHk5+dH/39nZ6fmcDh0jEZ/w//Zm5ubh1wjTdM0m82mdXZ2jlo2U3V2dmqVlZXRx3V1dZrdbo8+Hu37lEjfNUBrbm7WNE3TqqurZ9Q1mvF9jKOtBC4GjLYo8GhlM1VTUxMVFRXRxw6HA5fLhcvlGvX7lGjftdbW1iEbmUVqfTPhGs34xJiIP+zxGm1R4PEuGDwTOByOaNcBEP2u2O12+SMyyOC1C1pbW6OLvMyEazTjb74k4g97vEZbFHi8CwbPFINrQpWVlVRXVwPyR2Q4l8tFdXU1LS0t0c85E67RjK8xJuoPezxGWxQ40RcMrqmpoaSkhLKyMkD+iAxnt9uprKxkw4YNrF+//qKfd7pdoxmfGBP9hz0Woy0KPFkLBk8HTqeT7OzsaFIE+SMS4XK5hvTDFhcXo6oqTU1NM+IazfimdCL/sIdzOp3Rzu+Kigo2bNgQ3XExsiiw3W6nsbFxyKLAo5XNVJHrFBnvWVNTw913333J71OifNdcLld0OFPksc1mo7CwcETTeDpeo4QYxxgZkxb5YZeXl8fdP4SIHy6Xa8QSdjabLboNw2jfp0T6rtXU1ET///BxjNP9GiVEYhRCiPGY8X2MQggxXpIYhRBiGEmMQggxjCRGIYQYRhKjEEIMI4lRCCGGkcQoxBgMXxVGzGySGEXccjqdlJeXD5lhoZfIIhKXUl9fT3l5OS6Xa5IjEpNJEqOIW6qqUlFRgcPhoKWlhfLychRFGTLjYrCSkhKysrKoqqqa0NVaYi1yUFFRMWJ2DAxMIdywYUPcrRYjxmfGz5UWM0N+fj6VlZXAwDJggxd2gIFpfJFFCjZv3jyh566pqaG8vHzIc+vWrZPkN4NJjVFMGzabjQ0bNgAj+/xizW+eKK2trSPm8jY0NFBSUjIp5xP6k8QopoyqqtHal9PpxOl0XlZyKS8v59FHH73k65xOJ3l5eZSXl1NfXx/t/xtc02tpaaGiooL6+npqampG9A22tLRQUFAQ871VVcXpdFJRUSF9ijOMJEYxZZxOJ2VlZWzfvp3s7GwKCwtZt27duN+nrKwsmpgi7xtZPm0wh8NBcXExeXl5FBcXU1xcTElJSTQZR3Y8rKyspLi4mNbWVurr64e8x7Zt22I2291uNw6HA4fDQV5e3ojjxPQmiVFMmchiptnZ2eTn52Oz2S6rP9Bms+FwOC56EybW6yMcDkc0qW7fvn1IE7mysnJEPLH6EZ1OJ1u2bIm+b0NDw5CtEMT0J4lRTKmL1e7Ga8uWLVRXV+NyuSgsLJyAyEaqr6+P2dRvaGgY8hkm6jOJ+CGJUUypbdu2XfZNi8G1t0gNrbq6+pL7hQw+rr6+HofDEa11Du8bHDxmcngCHPyayPkjSTHS3yhmBhmuI6aUqqqXVbtqaWnh0UcfjSay4uJiKioqok3h+vp66urqaGpqoqqqirKysmjCbG1tjTafB2/NYLfbqauro6KiItrXGYntYhs0uVyuIfFHluV3Op3RbRDE9CcreIu4VV9fT35+/hUtex8ZiD38Bsql1NTU4HA4Luvc9fX10Y3ExPQkTWkhYmhubo67fUjE1JHEKGasyFjJyIbwY6WqasyxiyJxSFNaxC2n00ldXR0lJSXT5q5vTU0NDQ0NVFZWSo1zGpPEKIQQw0hTWgghhpHEKIQQw0hiFEKIYf5/CIfq0vy9kY8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 350x262.5 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_bulkflow_amplitude(False, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv_csiborg",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}