csiborgtools/notebooks/flow/flow_bulk.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calibrating the velocity field against observations "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copyright (C) 2024 Richard Stiskalek\n",
    "# This program is free software; you can redistribute it and/or modify it\n",
    "# under the terms of the GNU General Public License as published by the\n",
    "# Free Software Foundation; either version 3 of the License, or (at your\n",
    "# option) any later version.\n",
    "#\n",
    "# This program is distributed in the hope that it will be useful, but\n",
    "# WITHOUT ANY WARRANTY; without even the implied warranty of\n",
    "# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General\n",
    "# Public License for more details.\n",
    "#\n",
    "# You should have received a copy of the GNU General Public License along\n",
    "# with this program; if not, write to the Free Software Foundation, Inc.,\n",
    "# 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from flow_bulk import *\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enclosed monopole"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABiQElEQVR4nO3dd3hT9f4H8PdJ0qR7791C2QXKLksUBBQFHKg4AOWnoriuE66iV72KolfFhXq9ol4FlSsigrJEkD0KZRVKoaV7r3SnSc7vjzaRKqMjyUlO3q/n6QPJOTn59FDad79TEEVRBBERERE5PIXUBRARERGRZTDYEREREckEgx0RERGRTDDYEREREckEgx0RERGRTDDYEREREckEgx0RERGRTDDYEREREcmESuoC7JXRaERBQQG8vLwgCILU5RAREZGTEkURNTU1CA8Ph0Jx6TY5BruLKCgoQFRUlNRlEBEREQEAcnNzERkZeclzGOwuwsvLC0DLTfT29pa4GiIiInJWWq0WUVFR5mxyKQx2F2HqfvX29mawIyIiIsm1Z2gYJ08QERERyQSDHREREZFMMNgRERERyQSDHREREZFMMNgRERERyQSDHREREZFMMNgRERERyQSDHREREZFMMNgRERERyQSDHREREZFMMNgRERERyQSDHREREZFMMNgRERERyQSDHREREZFMMNhdxrG8KqlLICIiImoXBrvLOJxbJXUJRERERO3CYHcZqTmVUpdARERE1C4MdpeRmlsFURSlLoOIiIjoshjsLqOkRof8qgapyyAiIiK6LAa7dkjJZncsERER2T8Gu3Y4xGBHREREDoDBrh1SOIGCiIiIHACDXTucLKxBXZNe6jKIiIiILonB7jLCfFxhMIo4wvXsiIiIyM4x2F3GwChfAJxAQURERPaPwe4yBkb5AOA4OyIiIrJ/DHaXMTDKD0DLzFijkQsVExERkf1isLuMHqFecHNRQtuox9nSWqnLISIiIrooBrvLcFEqzOPsDnKcHREREdkxBrt2GBzT0h3LCRRERERkzxjs2sEU7LgDBREREdkzBrt2SIr2BQBkltWhok4nbTFEREREF8Fg1w6+7mp0D/YEwFY7IiIisl8Mdu00OLp1nB3XsyMiIiI7xWDXToNjW4PdOQY7IiIisk8Mdu1kmkBxJK8KOr1R4mqIiIiI/orBrp3iAz3g6+6CJr0RaYVaqcshIiIi+gsGu3YSBOGPcXacQEFERER2iMGuAwZxPTsiIiKyYwx2HTCkNdgdzK6AKIoSV0NERETUFoNdB/SP9IVKIaBY24T8qgapyyEiIiJqg8GuA9zUSvQN9wbAcXZERERkfxjsOojj7IiIiMheMdh1kGk9O+5AQURERPaGwa6DTMHuZGEN6pr0EldDRERE9AcGuw4K83FDhK8bDEYRR/KqpC6HiIiIyIzBrhNM4+y4bywRERHZEwa7ThgU7QsAOMRxdkRERGRHGOw6YVDr1mKHc6u4UDERERHZDQa7Tugd5g2NSoGq+mZkldVJXQ4RERERAAa7TlGrFEiM8AEAHMqpkrYYIiIiolYMdp1kXqiY4+yIiIjITjDYdZJpAsVhttgRERGRnWCw66Sk1gkU6UVa1HKhYiIiIrIDDHadFOLtighfNxhF4GhuldTlEBERETHYdUWSqTuWwY6IiIjsAINdF5i6Yw9lcwIFERERSY/BrgsGnddix4WKiYiISGoMdl3QN9wHapUCFXU6ZJfXS10OEREROTkGuy5QqxToF+4NgOvZERERkfQY7LrItG8sgx0RERFJjcGui0wTKLhQMREREUmNwa6LBsX4AgBOFdWgXseFiomIiEg6DHZdFObjhjAfVxiMIo7kVktdDhERETkxBjsL+GOhYo6zIyIiIukw2FmAeQJFdpW0hRAREZFTY7CzgD8mUFRyoWIiIiKSDIOdBfQN94aLUkB5nQ65FQ1Sl0NEREROisHOAlxdlOgb7gOA69kRERGRdBjsLIQLFRMREZHUGOwsxDwzlgsVExERkUQY7CxkUExLi93JQi0adAaJqyEiIiJnxGBnIeE+rgjx1kBvFHE0r0rqcoiIiMgJOVywO3fuHObOnYu4uDi4ubmhW7dueOGFF6DT6dqcd/ToUYwZMwaurq6IiorCkiVLrFqXIAhIimpd9iS3yqrvRURERHQhKqkL6KhTp07BaDTi448/Rvfu3XH8+HHce++9qKurw5tvvgkA0Gq1mDhxIiZMmICPPvoIx44dwz333ANfX1/cd999VqttUIwvNpwowqFsTqAgIiIi23O4YDd58mRMnjzZ/Dg+Ph7p6elYtmyZOdh9/fXX0Ol0+Oyzz6BWq9G3b1+kpqbirbfesmqwSzLPjK2CKIoQBMFq70VERET0Zw7XFXsh1dXV8Pf3Nz/es2cPxo4dC7VabX5u0qRJSE9PR2Wl9VrTEiN8oFIIKKttQl4lFyomIiIi23L4YHfmzBm89957uP/++83PFRUVISQkpM15psdFRUUXvE5TUxO0Wm2bj45qWajYGwDXsyMiIiLbs5tgt2DBAgiCcMmPU6dOtXlNfn4+Jk+ejBkzZuDee+/t0vsvXrwYPj4+5o+oqKhOXeePfWOrulQPERERUUfZzRi7J554AnPmzLnkOfHx8ea/FxQU4Morr8TIkSPxySeftDkvNDQUxcXFbZ4zPQ4NDb3gtRcuXIjHH3/c/Fir1XYq3CVF++Lz3ZwZS0RERLZnN8EuKCgIQUFB7To3Pz8fV155JQYPHozly5dDoWjb8JicnIxnn30Wzc3NcHFxAQBs3rwZPXv2hJ+f3wWvqdFooNFouvZJAOY9Y08X1cBoFKFQcAIFERER2YbddMW2V35+PsaNG4fo6Gi8+eabKC0tRVFRUZuxc7fffjvUajXmzp2LEydO4Ntvv8XSpUvbtMhZS2yAO9RKBRqaDciv4gQKIiIish27abFrr82bN+PMmTM4c+YMIiMj2xwTRREA4OPjg02bNmH+/PkYPHgwAgMD8fzzz1t1qRMTlVKB+CAPnCqqQXpRDaL83a3+nkREREQAIIimNERtaLVa+Pj4oLq6Gt7e3h167aPfHMaPqQV4alJPzL+yu5UqJCIiImfQkUzicF2xjqBHiBcAIKO4RuJKiIiIyJkw2FlBz9Zgl15cK3ElRERE5EwY7KzA1GJ3tqQWeoNR4mqIiIjIWTDYWUGknxvcXJTQGYzIrqiXuhwiIiJyEgx2VqBQCOgR4gmgZT07IiIiIltgsLOSBPM4OwY7IiIisg0GOysxTaA4zWBHRERENsJgZyU9Qk3BjjNjiYiIyDYY7KzENMYuq6wOTXqDxNUQERGRM2Cws5JQb1d4uapgMIrILK2TuhwiIiJyAgx2ViIIAsfZERERkU0x2FnRH+PsGOyIiIjI+hjsrKhHcMs4u/QiTqAgIiIi62OwsyJTi11GCVvsiIiIyPoY7KzINMYup6Ie9Tq9xNUQERGR3DHYWVGApwYBHmqIInCmhN2xREREZF0MdlbWI4QLFRMREZFtMNhZWU/OjCUiIiIbYbCzsoQQ08xYBjsiIiKyLgY7KzNNoMhgix0RERFZGYOdlSW0BruC6kZoG5slroaIiIjkjMHOynzcXBDm4wqArXZERERkXQx2NmBqteMOFERERGRNDHY20LN1AgVnxhIREZE1MdjZwB9r2THYERERkfUw2NkAgx0RERHZAoOdDZjWsiur1aG8tkniaoiIiEiuGOxswF2tQrS/OwBuLUZERETWw2BnI+yOJSIiImtjsLORHpwZS0RERFbGYGcjPUPZYkdERETWxWBnIz3MixTXQBRFiashIiIiOWKws5H4IA8oFQK0jXqU1HBmLBEREVkeg52NaFRKxAa0zIxNL2J3LBEREVkeg50NcZwdERERWRODnQ0lBDPYERERkfUw2NmQqcUunYsUExERkRUw2NmQaWZsRnENjEbOjCUiIiLLYrCzoZgAdygVAup1BpRyz1giIiKyMAY7G3JRKhDu6woAyKmol7gaIiIikhsGOxuL9m9Z8iSnnMGOiIiILIvBzsbMwY4tdkRERGRhDHY2FtUa7HI
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "cols = plt.rcParams['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "# Carrick+2015 must be multiplied by beta = 0.43\n",
    "r, Vmono = read_enclosed_monopole(\"Carrick2015\")\n",
    "Vmono *= 0.4\n",
    "plt.plot(r, Vmono[0], label=\"Carrick+2015\", color=cols[0])\n",
    "\n",
    "# r, Vmono = read_enclosed_monopole(\"Lilow2024\")\n",
    "# plt.plot(r, Vmono[0], label=\"Lilow+2024\", color=cols[1])\n",
    "\n",
    "# r, Vmono = read_enclosed_monopole(\"csiborg2X\")\n",
    "# ylow, ymed, yhigh = np.percentile(Vmono, [16, 50, 84], axis=0)\n",
    "# plt.fill_between(r, ylow, yhigh, color=cols[2], alpha=0.5, label=\"Manticore\")\n",
    "\n",
    "\n",
    "\n",
    "plt.xlabel(r\"$R ~ [\\mathrm{Mpc} / h]$\")\n",
    "plt.ylabel(r\"$V_{\\mathrm{mono}}(r \\leq R) ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "# plt.legend(loc=\"upper left\")\n",
    "\n",
    "plt.xlim(r.min(), r.max())\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"../../plots/enclosed_monopole.png\", dpi=450)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enclosed overdensity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9eZxkdX3vjz8/n7PU3vs6PTPMAgwgywSQEcTgQgJqYjCKirvXoDf5kqgk3oC/xPVGvPeqD0004RpDojciXjTBjYsiiguMDNuwDLMye+/7Uvs55/P741NV3dVV3dPVXdXLzHn6aKeps9SpU9XnvOq9vN5CKaXw8fHx8fHx8fFZ88iVPgAfHx8fHx8fH5/q4As7Hx8fHx8fH5/TBF/Y+fj4+Pj4+PicJvjCzsfHx8fHx8fnNMEXdj4+Pj4+Pj4+pwm+sPPx8fHx8fHxOU3whZ2Pj4+Pj4+Pz2mCL+x8fHx8fHx8fE4TfGHn4+Pj4+Pj43Oa4As7Hx8fHx8fH5/ThDUv7L761a+yadMmgsEgO3bsYNeuXfOuf++993LeeecRDAa56KKLuP/++4uWCyHK/vyv//W/Cuts2rSpZPnnPve5mrw+Hx8fHx8fH5+FsqaF3Xe+8x1uvfVWPvGJT/DUU09xySWXcN111zEwMFB2/UcffZSbbrqJ97///Tz99NPccMMN3HDDDTz//POFdXp7e4t+7rrrLoQQvOlNbyra16c//emi9f78z/+8pq/Vx8fHx8fHx+dUCKWUWumDWCw7duzgpS99KV/5ylcA8DyPDRs28Od//ufcdtttJeu/9a1vJR6P86Mf/ajw2Mte9jK2b9/OnXfeWfY5brjhBiYnJ3nooYcKj23atIkPf/jDfPjDH67uC/Lx8fHx8fHxWQLmSh/AYslkMjz55JPcfvvthceklFx77bXs3Lmz7DY7d+7k1ltvLXrsuuuu47777iu7fn9/Pz/+8Y/5xje+UbLsc5/7HJ/5zGfYuHEjb3/72/nIRz6CaZY/nel0mnQ6Xfhvz/MYGRmhubkZIcSpXqqPj4+Pj4/PaYJSisnJSdatW4eU1U+crllhNzQ0hOu6tLe3Fz3e3t7Ovn37ym7T19dXdv2+vr6y63/jG98gFovxx3/8x0WP/8Vf/AWXXnopTU1NPProo9x+++309vbyxS9+sex+7rjjDj71qU8t9KX5+Pj4+Pj4nOacOHGC9evXV32/a1bYLQd33XUX73jHOwgGg0WPz4z6XXzxxdi2zQc/+EHuuOMOAoFAyX5uv/32om3Gx8fZuHEjJ/753XiBGLt6HeKj5zPYv4XWtiO8JPwAg3vqSUTacILheY9xJJ7hdzY20hYrfV6AvskJHnlxgKjpERsbYu+rX89ka0clp2Feukdcrj4vyPZNdtX2WUAp1j35LdSeY6SeHkVtaCPRfn7RKmJsEpQi/YG3QMAqPD7yrMORH6eo22Rwzk2h6h/bHBx/IM3g01k6rrTpeuXSz4kYGcf+xn0gJKqlYf6VlUIc68F9+aU4r7piyc+9VLyU4sC3k8T7PKyo4Lx3hbEb1naEeiQ5gilNXtrxUqJ2dKUPx8fHZw0yMTHBhg0biMViNdn/mhV2LS0tGIZBf39/0eP9/f10dJQXLh0dHQte/9e//jX79+/nO9/5zimPZceOHTiOw9GjR9m2bVvJ8kAgUFbw1YVtvECAcFCi7DAhO0I4ECIatEhYNtgBnECwZLuZJLOCSCBALFRevExmM4Rsm7DhErIsAuEomXDklK9poZgJBysUJBqb/zgXi7VuC+Fjx4jbJspUyOCs8xh1EfEklmlAbPpG661zCNkGRloQjS3fDVgmJCHbpml9kGjMOvUGp9rfwAiWp1DrW8AwTrm+aGlEHT5B5torIbx8grYsMbjkT6Ls+9cEyUGPE/8hOP99Yey6tduzFYlGODl1kgFvgM5Yp19K4ePjs2hqdf1Ys1dY27a57LLLipoaPM/joYce4sorryy7zZVXXlm0PsCDDz5Ydv1/+Zd/4bLLLuOSSy455bHs3r0bKSVtbW0VvgqA4t4VQeW9LPNtoZRebkmJYZhQ7Xy+EDhudXc5k0y0DWkLMAQyky1dwTLBcRCJVNHDdr1+nelxxXL2B6VGPACCzdX5gxWjE6C8BYk6ANVYhxgZRx48VpXnXypWWLDtXSECjYL0qGLvvyVIj3krfViLRghBa6iV4xPH6YuXL+Hw8fHxWUnWrLADnRL953/+Z77xjW+wd+9e/vRP/5R4PM773vc+AN797ncXNVd86EMf4oEHHuALX/gC+/bt45Of/CRPPPEEt9xyS9F+JyYmuPfee/mTP/mTkufcuXMnX/rSl3jmmWc4fPgw3/rWt/jIRz7CO9/5ThobGxf9WhSLFwLz6RaVU3YCgSFMVJW/IUgUqWzthFM20gzhAEIKyDqlL9YyEdkywq5Ov07lgJNYHmHnZRWZcf1cwabq/GmJ0Qmo5LNhGGAaGM8eAG91CCi7TnLeu8PY9YL0iGLvvyZIDa+OY1sMQTOIFJJDY4fIemW+bPj4+PisIGta2L31rW/l85//PB//+MfZvn07u3fv5oEHHig0SBw/fpze3t7C+ldddRV33303X/va17jkkkv47ne/y3333ceFF15YtN977rkHpRQ33XRTyXMGAgHuuecerrnmGl7ykpfwd3/3d3zkIx/ha1/72tJeTE576IhdZeJr3ohd7kfkdlt1YSch7dROOGUiLahIFGF44LgINSs8KKUWr4lk8cOmwIrq15oXW7UmNarFihEAM1KliF3fIFiVpXRVUwPyeA+iZ7Aqx1ANAo2S8/9LmGCzIDOuxV2iv4ah3hrTHGpmKDnkR+18fHxWHWu2xi7PLbfcUhJxy/Pwww+XPHbjjTdy4403zrvPD3zgA3zgAx8ou+zSSy/lt7/9bcXHORclkmMRemC+VKOaIRiVEFBtYScEmRoKO88KkWloQZqHcR0P4bkoWfqxFclUyWN2vSA7pciMe0TWLSyVuRTSuShUoElWp3Yi6yCHRiFUYRNGOAj9Qxj7DuOsbz/1+stEoF5y/vvC7Ps/SZL9Hnv/LcG2d4aJdtXuvZk66dK3M8PkUZdwp6T+bJOGc0yCzUv7TmtKk4AR4PDYYTrCHVjG0uspfXx8fKrBmhd2pwM6qiZm/Ffl28+5LK/sFIDAq3KNnRTUVNgBZBo2EAjswk264JWP8sxOxYKus4t3e6SXK2I3kkvDLlE05BFjE5BKo+oqbP4QAmJR5PMH4KrtK99EMQMrKjn/PWH2fytBvNtj3zcSrH91gGCzJNAgsOslhr0wUexmFJkJhRUVGIHpQmTlKUb3OfTtzDB1YjrlO37IZfyQy/EH0gQaBXVbTKQJblrhZfT+vKxO44daJcFWSajVINgoEEbpMTUFm+iOd3N84jhbG7eilMJVLmaZLx4+Pj4+y4V/BVotFCJri992vkVSF9pVP2IndelbLcmGmwhGTNRQtjQVC2BIGJ8seThQn0/FLk89V75urGr1dWMTkMpAW3krm/lQTXWI473IQ8fxLi7t1F5JzLDgvHeHOfDtJJNHtdAqWh4RRNZJmi+yaDzPLBF6yQGX/sezDD2Txcvox6QJVkxgRSXZSY/0mP7kCwnNF5s0X2SR6PMYO+gwddwlPaoYfHJh9XHCgHCHpG6TSWyTQWyjgREQGNKg3q5n/9h+onaU4eQwE5kJXtrxUgxZ+wixj4+PTzl8YbdqyN28hKo4ZjffFtNpWoVCoET1I3ZpR3ee1qp12zNtZNhEqAyiXMTOshDjUyUP5201lq3GrtARW73GCaEUajFR1kITxX68i86tuqBfKkZAsO0dIXofyRDv9ciMeaTHPNw0OHHF+EGX8YMu0oLG87Qw8xxF/64sk0enPwPSBM/RP+lRRXpULzNDgraXWrS91MKO6fNXvxU6X27jphUTRxwdzRO6JlLaAsMWSBPSYx7JQf2TGvTwHIh3e8S7M/Q+AgiIrJM0bjNpe2kdcbeX54aeYzwzji1tJjOTNAQbVuCs+vj4+PjCbpVQKjwq6ZKdtyt25i9CVL95QghcT+EpKJOtqgqeYUP
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "cols = plt.rcParamsDefault['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "r, overdensity = read_enclosed_density(\"csiborg1\")\n",
    "c = cols[0]\n",
    "ymin = np.percentile(overdensity, 16, axis=0)\n",
    "ymed = np.median(overdensity, axis=0)\n",
    "ymax = np.percentile(overdensity, 84, axis=0)\n",
    "plt.fill_between(r, ymin, ymax, color=c, alpha=0.5, label=\"CB1\")\n",
    "\n",
    "\n",
    "r, overdensity = read_enclosed_density(\"csiborg2_main\")\n",
    "c = cols[1]\n",
    "ymin = np.percentile(overdensity, 16, axis=0)\n",
    "ymed = np.median(overdensity, axis=0)\n",
    "ymax = np.percentile(overdensity, 84, axis=0)\n",
    "plt.fill_between(r, ymin, ymax, color=c, alpha=0.5, label=\"CB2\")\n",
    "\n",
    "\n",
    "r, overdensity = read_enclosed_density(\"csiborg2_random\")\n",
    "c = cols[2]\n",
    "ymin = np.percentile(overdensity, 16, axis=0)\n",
    "ymed = np.median(overdensity, axis=0)\n",
    "ymax = np.percentile(overdensity, 84, axis=0)\n",
    "plt.fill_between(r, ymin, ymax, color=c, alpha=0.25, label=\"Random\", zorder=-1)\n",
    "\n",
    "r, overdensity = read_enclosed_density(\"csiborg2X\")\n",
    "c = cols[3]\n",
    "ymin = np.percentile(overdensity, 16, axis=0)\n",
    "ymed = np.median(overdensity, axis=0)\n",
    "ymax = np.percentile(overdensity, 84, axis=0)\n",
    "plt.fill_between(r, ymin, ymax, color=c, alpha=0.5, label=\"Manticore\")\n",
    "\n",
    "plt.xlabel(r\"$R ~ [\\mathrm{Mpc} / h]$\")\n",
    "plt.ylabel(r\"$\\delta_{\\rm enc} (r \\leq R)$\")\n",
    "plt.legend()\n",
    "\n",
    "plt.xlim(0, r.max())\n",
    "plt.ylim(-0.15, 0.075)\n",
    "plt.axhline(0, c='k', ls='--', lw=1)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"../../plots/enclosed_density.png\", dpi=450)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enclosed bulk flow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(101, 500, 3)\n"
     ]
    }
   ],
   "source": [
    "r, Vmag, l, b = read_enclosed_flow(\"csiborg1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 100, 3)\n"
     ]
    }
   ],
   "source": [
    "r, Vmag, l, b = read_enclosed_flow(\"csiborg2X\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABW0AAAHqCAYAAAB/bWzAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdeXwV9bn48c/MnDX7QlYIhCUoqwtuiNVWbdGi195yu3ordrPXaq21t7W21ra01m7XWm8tvfVn1S5qF5e27jsqICCrErYEkrBkT86+z8zvj28IBAIkIclJwvP2dV54zpkz8z2QnPOdZ57v82i2bdsIIYQQQgghhBBCCCGEGBH0dA9ACCGEEEIIIYQQQgghxEEStBVCCCGEEEIIIYQQQogRRIK2QgghhBBCCCGEEEIIMYJI0FYIIYQQQgghhBBCCCFGEAnaCiGEEEIIIYQQQgghxAgiQVshhBBCCCGEEEIIIYQYQSRoK4QQQgghhBBCCCGEECOIBG2FEEIIIYQQQgghhBBiBHGkewAjiWVZ7N+/n+zsbDRNS/dwhBBCDDPbtgkGg5SXl6Prcl2zL+S7UwghTm7y3Tkw8v0phBAnr75+d0rQ9hD79++noqIi3cMQQgiRZnv27GHChAnpHsaoIN+dQgghQL47+0u+P4UQQhzvu1OCtofIzs4G1F9aTk5OmkcjhBBiuAUCASoqKrq/D8TxyXenEEKc3OS7c2Dk+1MIIU5eff3ulKDtIQ4sS8nJyZEvTiGEOInJMsW+k+9OIYQQIN+d/SXfn0IIIY733SlFh4QQQgghhBBCCCGEEGIEkaCtEEIIIYQQQgghhBBCjCAStBVCCCGEEEIIIYQQQogRRGraCiHEcZimSTKZTPcwxCBwOp0YhpHuYQghhBBCCCGEEMckQVshhDgK27ZpamrC5/OleyhiEOXl5VFaWjomGqYsW7aMZcuWUVdXB8CsWbO44447uPzyywF4//vfz/Lly3u85ktf+hK//e1vu+83NDRw/fXX89prr5GVlcWSJUu46667cDhkiiCEEEIIIYQQ6SJnZEIIcRQHArbFxcVkZGSMiSDfycy2bSKRCC0tLQCUlZWleUQnbsKECfzkJz+hqqoK27Z5+OGHueqqq9iwYQOzZs0C4Itf/CJLly7tfk1GRkb3/5umyaJFiygtLWXlypU0NjZyzTXX4HQ6+fGPfzzs70cIIYRIN9M0+f73v8+f/vQnmpqaKC8v59prr+X222/vngvats33vvc97r//fnw+HwsWLGDZsmVUVVWlefRCCCHGEgnaCiFEL0zT7A7YFhYWpns4YpB4vV4AWlpaKC4uHvWlEq688soe9++8806WLVvG22+/3R20zcjIoLS0tNfXv/jii1RXV/Pyyy9TUlLC6aefzg9/+ENuvfVWvv/97+NyuYb8PQghhBAjyU9/+lOWLVvGww8/zKxZs3jnnXf47Gc/S25uLjfddBMAP/vZz7j33nt5+OGHmTx5Mt/97ndZuHAh1dXVeDyeNL8DIYQQY4U0IhNCiF4cqGF7aFaiGBsO/JuOtTrFpmny2GOPEQ6HmT9/fvfjf/7znxk3bhyzZ8/mtttuIxKJdD+3atUq5syZQ0lJSfdjCxcuJBAIsGXLlmEdvxBCCDESrFy5kquuuopFixZRWVnJf/zHf/ChD32INWvWACrL9p577uH222/nqquuYu7cufzhD39g//79PPXUU+kdvBBCiDFFMm2FEOIYpCTC2DPW/k3fffdd5s+fTywWIysriyeffJKZM2cC8OlPf5pJkyZRXl7O5s2bufXWW9m+fTtPPPEEoEqAHBqwBbrvNzU19Xq8eDxOPB7vvh8IBIbibQkhhBBpcf755/O73/2OHTt2MH36dDZt2sRbb73F3XffDcDu3btpamri0ksv7X5Nbm4u5557LqtWreKTn/xkr/uV708hhBD9JUFbIYQQYhQ75ZRT2LhxI36/n7///e8sWbKE5cuXM3PmTK677rru7ebMmUNZWRmXXHIJtbW1TJ06dUDHu+uuu/jBD34wWMMXQgghRpRvfetbBAIBTj31VAzDwDRN7rzzTq6++mrg4EXN3i56Hu2CJ8j3pxBCiP6ToK0QQvRTIJYkljCH5Vgel0GOxzksxxKjk8vlYtq0aQDMmzePtWvX8qtf/Yr/+7//O2Lbc889F4CamhqmTp1KaWlp93LPA5qbmwGOWgf3tttu45Zbbum+HwgEqKioGJT3IoQQQqTbX//6V/785z/zyCOPMGvWLDZu3MjNN99MeXk5S5YsGfB+5ftTCCFEf0nQVggh+iEQS/K/r+ykI5wYluMVZLr4yiVV/Q7cNjU1ceedd/LMM8+wb98+iouLOf3007n55pu55JJLqKyspL6+HgBd1ykpKeHyyy/nF7/4Bfn5+QDEYjH+67/+i3Xr1rF161auuOIKqdU2CliW1WP55aE2btwIQFlZGQDz58/nzjvv7G7MBvDSSy+Rk5PTXWLhcG63G7fbPfgDF0IIIUaAb3zjG3zrW9/qLnMwZ84c6uvrueuuu1iyZEn3Rc3m5ubu79MD908//fSj7le+P4UQQvSXBG2FEKIfYgmTjnACt8Mgw2UM6bEiXceKJcx+BW3r6upYsGABeXl5/PznP2fOnDkkk0leeOEFbrjhBrZt2wbA0qVL+eIXv4hpmuzYsYPrrruOm266iT/+8Y+Aamzl9Xq56aabePzxx4fkPYoTc9ttt3H55ZczceJEgsEgjzzyCK+//jovvPACtbW1PPLII3z4wx+msLCQzZs387WvfY0LL7yQuXPnAvChD32ImTNn8pnPfIaf/exnNDU1cfvtt3PDDTfIiaUQQoiTUiQSQdd79us2DAPLsgCYPHkypaWlvPLKK91B2kAgwOrVq7n++uuHe7hCCCHGMAnaCiHEAGS4DDLdQ/8RGk/1vwzDl7/8ZTRNY82aNWRmZnY/PmvWLD73uc9138/Ozu7OFhk/fjxLlizh0Ucf7X4+MzOTZcuWAbBixQp8Pt8A34UYKi0tLVxzzTU0NjaSm5vL3LlzeeGFF/jgBz/Inj17ePnll7nnnnsIh8NUVFSwePFibr/99u7XG4bB008/zfXXX8/8+fPJzMxkyZIlLF26NI3vSgghhEifK6+8kjvvvJOJEycya9YsNmzYwN133909h9I0jZtvvpkf/ehHVFVVMXnyZL773e9SXl7ORz7ykfQOXgghxJgiQVshhBhDOjo6eP7557nzzjt7BGwPyMvL6/V1+/bt41//+ld3zVMxOjzwwANHfa6iooLly5cfdx+TJk3i2WefHcxhCSGEEKPW//7v//Ld736XL3/5y7S0tFBeXs6XvvQl7rjjju5tvvnNbxIOh7nuuuvw+XxccMEFPP/883g8njSOXAghxFijH3+T4VFZWYmmaUfcbrjhBkDVVrzhhhsoLCwkKyuLxYsXdzdLOaChoYFFixaRkZFBcXEx3/jGN0ilUul4O0IIkRY1NTXYts2pp5563G1vvfVWsrKy8Hq9TJgwAU3TuPvuu4dhlEIIIYQQI1N2djb33HMP9fX1RKNRamtr+dGPfoTL5ereRtM0li5dSlNTE7FYjJdffpnp06encdRCCCHGohETtF27di2NjY3dt5deegmAj33sYwB87Wtf41//+hd/+9vfWL58Ofv37+ejH/1o9+tN02TRokUkEglWrlzJww8/zEMPPdTjiuiAhFogGTuxfQghxDCxbbvP237jG99g48aNbN68mVdeeQWARYsWYZr9L8kghBBHY4XD2PK5IoQQQgyL2PbtdPzpT1iJ4WmcLIQYOiMmaFtUVERpaWn37emnn2bq1KlcdNFF+P1+HnjgAe6++24uvvhi5s2bx4MPPsjKlSt5++23AXjxxReprq7mT3/6E6effjqXX345P/zhD7nvvvtInMiH1erfwqZHj7+dEEKMAFVVVWia1t1s7FjGjRvHtGnTqKqq4uKLL+aee+5h5cqVvPbaa8MwUiHEWGeGwgReeYWWX91LdNPmdA9HCCGEGPOSzc34Hn+cyLr1JGpr0z0cIcQJGjFB20MlEgn+9Kc/8bnPfQ5
      "text/plain": [
       "<Figure size 1400x500 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(14, 5))\n",
    "cols = plt.rcParamsDefault['axes.prop_cycle'].by_key()['color']\n",
    "thin_line_kwargs = {\"lw\": 0.3, \"zorder\": 0, \"ls\": \"dotted\"}\n",
    "\n",
    "# CSiBORG1\n",
    "r, Vmag, l, b = read_enclosed_flow(\"csiborg1\")\n",
    "c = cols[0]\n",
    "for i, y in enumerate([Vmag, l, b]):\n",
    "    ymin, ymed, ymax = np.percentile(y, [16, 50, 84], axis=0)\n",
    "    axs[i].fill_between(r, ymin, ymax, color=c, alpha=0.5, label=\"CB1\")\n",
    "\n",
    "# CSiBORG2\n",
    "r, Vmag, l, b = read_enclosed_flow(\"csiborg2_main\")\n",
    "c = cols[1]\n",
    "for i, y in enumerate([Vmag, l, b]):\n",
    "    ymin, ymed, ymax = np.percentile(y, [16, 50, 84], axis=0)\n",
    "    axs[i].fill_between(r, ymin, ymax, color=c, alpha=0.5, label=\"CB2\")\n",
    "\n",
    "# Manticore\n",
    "r, Vmag, l, b = read_enclosed_flow(\"csiborg2X\")\n",
    "c = cols[3]\n",
    "for i, y in enumerate([Vmag, l, b]):\n",
    "    ymin, ymed, ymax = np.percentile(y, [16, 50, 84], axis=0)\n",
    "    axs[i].fill_between(r, ymin, ymax, color=c, alpha=0.5, label=\"Manticore\")\n",
    "\n",
    "\n",
    "# Random\n",
    "r, Vmag, l, b = read_enclosed_flow(\"csiborg2_random\")\n",
    "c = cols[2]\n",
    "ymin = np.percentile(Vmag, 16, axis=0)\n",
    "ymax = np.percentile(Vmag, 84, axis=0)\n",
    "ymed = np.median(Vmag, axis=0)\n",
    "axs[0].fill_between(r, ymin, ymax, color=c, alpha=0.35, zorder=-1,\n",
    "                    label=\"Random\")\n",
    "\n",
    "\n",
    "axs[0].legend()\n",
    "for i in range(3):\n",
    "    axs[i].set_xlabel(r\"$R ~ [\\mathrm{Mpc} / h]$\")\n",
    "    axs[i].set_xlim(0, r.max())\n",
    "\n",
    "axs[0].set_ylim(0)\n",
    "axs[0].set_ylabel(r\"$V_{\\rm enc} (r \\leq R) ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[1].set_ylabel(r\"$l_{\\rm enc} (r \\leq R) ~ [\\mathrm{deg}]$\")\n",
    "axs[2].set_ylabel(r\"$b_{\\rm enc} (r \\leq R) ~ [\\mathrm{deg}]$\")\n",
    "axs[1].set_ylim(0, 360)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"../../plots/enclosed_flow.png\", dpi=450)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enclosed overdensity $\\texttt{BORG2}$ all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)\n",
    "nsims = paths.get_ics(\"borg2_all\")\n",
    "\n",
    "r, overdensity = read_enclosed_density(\"borg2_all\")\n",
    "overdensity_A = overdensity[:80]\n",
    "overdensity_B = overdensity[80:][::30]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 7500,  7600,  7700,  7800,  7900,  8000,  8100,  8200,  8300,\n",
       "        8400,  8500,  8600,  8700,  8800,  8900,  9000,  9100,  9200,\n",
       "        9300,  9400,  9500,  9600,  9700,  9800,  9900, 10000, 10100,\n",
       "       10200, 10300, 10400, 10500, 10600, 10700, 10800, 10900, 11000,\n",
       "       11100, 11200, 11300, 11400, 11500, 11600, 11700, 11800, 11900,\n",
       "       12000, 12100, 12200, 12300, 12400, 12500, 12600, 12700, 12800,\n",
       "       12900, 13000, 13100, 13200, 13300, 13400, 13500, 13600, 13700,\n",
       "       13800, 13900, 14000, 14100, 14200, 14300, 14400, 14500, 14600,\n",
       "       14700, 14800, 14900, 15000, 15100, 15200, 15300, 15400])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nsims[:80]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOz9d7gc+XnfiX4qdnXO3SdH5AFmgAnkDIccckiapCKtsFrJWlmyZXPvXdGWTdnX0tqWvOukZ4Nsybas62d9fVcPLUuWZF1LlsQgDkmRnByAQQ4nnz59TudUOd0/GnMAEMAQlEhOYH2eB88D9Klq/Kq6+tS33vB9hTAMQyIiIiIiIiIiIt7yiG/0AiIiIiIiIiIiIr4xRMIuIiIiIiIiIuJtQiTsIiIiIiIiIiLeJkTCLiIiIiIiIiLibUIk7CIiIiIiIiIi3iZEwi4iIiIiIiIi4m1CJOwiIiIiIiIiIt4mRMIuIiIiIiIiIuJtgvxGL+DtRhAE7OzskE6nEQThjV5ORERERERExFucMAwZDodMTU0hiq8fk4uE3TeYnZ0dZmdn3+hlRERERERERLzN2NraYmZm5nW3iYTdN5h0Og2MT34mk3mDV/Nn5zPndvhf/uAcT4UfIyZ4OGIC9e9deaOXFRERERER8U3Hcn0s1yeXUG/7mW57tHUbWRTRFAnCEEEQSGsya60RubhKOaPdss9Kc8RSKcnl3QHZuMpkLs5GWycdkwmBTFyhb7hYnkcQwlwhecv+g8GA2dnZfY3xekTC7hvMa+nXTCbzlhZ2hYJJLJ4i7iukBB9HElHfwscTERERERFxr3TbOvMTydtebw0tGlbAfKVILqHSHtmIgkA+qXJ6q8fcRJlyOnbLPlsdgwcWJ3hlq8fBmQr5ZIxaz2SmWiQIQoqpGDs9k1hcYWS43D+Tu+u67qXEKxJ2EXdkvaWj2x6uPL5EBN9+g1cUERERERHxzcd0fDKacstrPcPBdH0Mx+dQNU1clW4RdSvNIdM57TZR99p7vbTZ5ehEhpSm0Bha5OMKfctlMhunZzg4fkB75PDQfP6Oa+qM7v0eHHXFRtyREzM5cskYHhIAoRg9A0REREREvP1pjWzyyRsp2O2ugSgKJBSZQkIlrkr0DAfhuqir9008D8rpW9OvYRhS65ns9A0WiilSmkJ7ZBOTJVojh2IyRhiG2F5Ao2/dVdRtdgy2usY9rz8SdhF3xHR8XD/AuR7UFQL/DV5RRERERETENxfL9UnGxve9IAjZaOuUUjHSMRnb88knVQzHw/VDCkmV1tBir29xePL22rftrommiMRkiXI6th/hI4RsXEGVRWo9k1e3ejyyWLjjei7vDthojZjNJ+75GKIwTMQdGZgOjhfihjIIIBC80UuKiIiIiIj4ptIc2szk4xiOR89wmc7FkSVxv/nB8wO6118fWi7du9TEdXUHVRbYG9jcP5OjPbIRBIFcQmFvYDOR1bA9n/Pbfd5/tLpfO+cHIQPTpW86nK8NSMZkpnIahM49H0Mk7CLuyGI5SSEVwxtK11+JPPkiIiIiIt6+OF5AXJXoGi5hGDKZ1RAEAccLKFzvjt3qmiyWkgRByGpLZzoXRxSF296nrTs0BhbvOlCiObRRJIFcQuVaY8SBSgqAs9t9Ts3n8MOQ7sACAWw3gDDgWmPE/bNZZq93x67XW/d8HFEqNuKO6LaHYXu4r6ViiVKxERERERFvX5ojmzAc/72Yiu1H0XZ6JrmEwkpTZ7E0Flq1nkk+rlBK3WiW8PyAxsDiXK3HxXqfQxNpnltt4wcBGU2h3jeZK4xTqh3dZqurgyCwN7DwghDLCVBlkQu7Q548XGG2kES3PTbbBnFF4l6JhF3EHYnJEhDiXm+eEF+72iMiIiIiIt5m+EGIbrtIokDhpsaJoeWST6psdoz9SFtHd3CDgJnrdW9BENIYWHT0cbrUdHw+cKSKbnucmsuTS6jsDS32+jYj22Oro/OfntvkxHQOwx7X7U3l4lQyKpfqAz583yQhY+87PwyZKyaoD8x7PpYoFRtxR8ppjUJSxXWuR+wEIAjga4wyiYiIiIiIeKux0hiR0pRbRB1A33QpJNX9LlnD8djtW/sp2Nb1KF9Kk4krEi+sd5gpJKgPLJZKyf2on277nJzLYTgeV/dMHlnMc6Byo+GiObRZaQw5OZtnpTmikFBZLo+FZK1rYLn3njWLhF3EHbFcn7bu7qdiAQhcEGN33ykiIiIiIuItRs9w8MOQqVz8ltcbQ4tCUqXet24SWSbltEpak2kMLVRJ3J9O0RqOU6quH+5H9wAaA4vFUpKO7tDRbXqGy4ePTwLjlOyL610s12ehlKQ5sigkY5iuz6htMLJdzu30Wcree1AlEnYRdyQZk0nFJFz7pry+74AcCbuIiIiIiLcHA8tlt2/t177dTBiC5QYsFsd1dY2BhSQJpDWFvaFFPqGOR4oBrh/w0kaPyZzGcvnGxIogCHH8gL7p4voBuwObQxMZ1lsjWiOHoeUyndW476s6a10/YLOtszu0+K77p/Cse/exi4RdxB1JqBL5uII7uOkS8d03bkERERH3RN9wcfyxPdHI9sjGFRwvwPWDO3bwRUR8u2J743mwhaS67133Glsdg0p6PPqrkFTpGQ5t3aGUUmnrDpMZ7Zbv0mpjRDGpcGwyc8vYr7W2TlKVMRyPxsDGtL39ZopialzP91pDxmuYjs+1xgDLDXjiYAU/CFnrRMIu4s+J54fU+9b+5AlgHLGLiIh40zKwxqLutbFGXz3eaG9g4XiRwIuICMOQ3b5FRlNu+y6EYUhclegYDkvXU7DtkYOmiAQhTH9VyrZnOJzf6fO+I1Vk6UbKtDW06OoOYcg4dTswmSulmM0naAxt4opEKaXe4mHXHFp0RjZBCA8vFPYbMmbvEFG8G1ElfMQdiSsSU7n4/uQJALxoXmxExJsVw/GwHJ9SSr3rNtWMxmwhwe7AotYzCYKo2z3i25OVps58MYnrB2Tjt86FXW3pJBQJ6brga49seoZDLq5SzWi3vdeL6x2OTWX3Gy88P2C7Y/DU5San5vIslZKsNIbMl1IcnUjTMRwUWSCuSPtCsD2yeXW7hyyJKLLE4WqGnb6FJAoUU19fCVQk7CLuiCwJSKJwa/OEF0XsIiLejDheQM9wKadjt6SB7sZULs50Ls5aWyeMrIwivs1Yb+kcqKQYWu5t0bogCMknVBpDm8r12a+1nkklG7tlfuxrXKz3EQWBI5MZYNx4eLE+pDm0+aGHZxEFOFvrMbI9ZgsJmiMbSRAgFMgmFPqmy+mtLn4QMl8cGx8nYzIDy6WYVPcbMxzv3rtiI2EXcUdEQaAxtPHCm1Kx3r3n+CMiIr41BME4pTSR0e5J1N3McjnFSlP/Jq0sIuLNR2NoMZ0fp1JN17/FYBhgva0jEDKZG4u6nZ5Jd+Qwm0/e9l4Dw+HMVp+HF/LAuOHhzFaPiazGYnk8fmylOSIIQk7O5jEdH1kU8IKQfFJhpTmkZzjM5BNUMhqeHzCyPWKySCWj7TdmtEc2n7u4d8/HGAm7iDsiigKHqqmvsjuJ5sVGRLzZWG/rzBb+7DVz88UE7VFUZhHx9qdvusiiiCKJ2J6P598arfaDkGxcoa2710364cxWjycOl297rzAMeWmzy6NLeTJxFT8IeWalxUNzeYaWi6ZItEYOmiKhO+NomyqJeH6IH4TUuiapmMJEVqOUinF1b8jQ9qhkNIqpGGEY0h7ZbHcN6j2LdywW7/k4I2EXcVdGlr8/eQIA996dryMiIr65eP54nuTiTSaofxYUScRw/KjeLuJtjen4DC13vw6ub7q3+dZtdgyCMNy3K3lls8tyOXXH79fLm10EYK6QwvUCnl5p8cShCo2RTSmljkWkJFw3M9aQJZG2blPrmWQ0hYQqU0nHkASBi/UBqiyyXE4hiwKNgUW9b5GMySiSCAJfV51dJOwi7orlB7dG7NwoZRMR8WagObSp9y0OVO580/l6mS0kWG9H3++Ity/d6ylPGJcvmM6tNWuOF5BQJUwnQBAEXM+nb7gcmkjf9l7na336hssDs3kGpsvp7R7vOVi
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "cols = plt.rcParamsDefault['axes.prop_cycle'].by_key()['color']\n",
    "\n",
    "# r, overdensity = read_enclosed_density(\"csiborg1\")\n",
    "c = cols[0]\n",
    "for i in range(len(overdensity_A)):\n",
    "    plt.plot(r, overdensity_A[i], c=c, lw=0.2, zorder=0, ls=\"dotted\")\n",
    "plt.plot(r, np.median(overdensity_A, axis=0), c=c, lw=2, label=r\"Chain index: $7500$ to $15000$\")\n",
    "\n",
    "c = cols[1]\n",
    "for i in range(len(overdensity_B)):\n",
    "    plt.plot(r, overdensity_B[i], c=c, lw=0.2, zorder=0, ls=\"dotted\")\n",
    "plt.plot(r, np.median(overdensity_B, axis=0), c=c, lw=2, label=r\"Chain index: $15000$ to $17000$\")\n",
    "\n",
    "\n",
    "\n",
    "plt.xlabel(r\"$R ~ [\\mathrm{Mpc} / h]$\")\n",
    "plt.ylabel(r\"$\\delta_{\\rm enc} (r \\leq R)$\")\n",
    "plt.legend()\n",
    "\n",
    "plt.xlim(0, r.max())\n",
    "plt.ylim(-0.14, -0.005)\n",
    "plt.axhline(0, c='k', ls='--', lw=1)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"../../plots/enclosed_density_BORG2.png\", dpi=450)\n",
    "plt.show()"
   ]
  },
  {
   "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
}