csiborgtools/notebooks/field_velocity_fof_sph.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tests of velocities of haloes in CSiBORG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import csiborgtools\n",
    "\n",
    "from scipy.stats import spearmanr, binned_statistic\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FoF vs SPH velocity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "kind = \"main\"\n",
    "nsnap = 99\n",
    "nsim = 17417\n",
    "\n",
    "field_reader = csiborgtools.read.CSiBORG2Field(nsim, kind)\n",
    "catalogue = csiborgtools.read.CSiBORG2Catalogue(nsim, nsnap, kind)\n",
    "boxsize = csiborgtools.simname2boxsize(\"csiborg2_main\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity_field = field_reader.velocity_field(\"SPH\", 1024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos = catalogue[\"cartesian_pos\"] / boxsize\n",
    "vel = catalogue[\"cartesian_vel\"]\n",
    "mass = catalogue[\"totmass\"]\n",
    "\n",
    "spherical_pos = catalogue[\"spherical_pos\"]\n",
    "RA = np.deg2rad(spherical_pos[:, 1])\n",
    "dec = np.deg2rad(spherical_pos[:, 2])\n",
    "\n",
    "def project_radial(vx, vy, vz, RA, dec):\n",
    "    return vx * np.cos(dec) * np.cos(RA) + vy * np.cos(dec) * np.sin(RA) + vz * np.sin(dec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "vx, vy, vz = csiborgtools.field.evaluate_cartesian_cic(velocity_field[0], velocity_field[1], velocity_field[2], pos=pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABdEAAAGNCAYAAAD+VjiVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gVxdfA8e/ubemNFkIJTemgoiJKUUE6KoqKoBRBkCqCFSu8KApIkSoWQH8iiF2QKlVARAQpCtJ7Cun9lp33j5vc5KZAkBLK+fjkkczO7s7uTWY3Z2fPaEophRBCCCGEEEIIIYQQQgghCtBLugFCCCGEEEIIIYQQQgghxJVKguhCCCGEEEIIIYQQQgghRBEkiC6EEEIIIYQQQgghhBBCFEGC6EIIIYQQQgghhBBCCCFEESSILoQQQgghhBBCCCGEEEIUQYLoQgghhBBCCCGEEEIIIUQRJIguhBBCCCGEEEIIIYQQQhRBguhCCCGEEEIIIYQQQgghRBEkiC6EEEIIIYQQQgghhBBCFEGC6EIIIYS4rmiaxltvvVXSzbhgn3/+ObVq1cJisRASElLSzRGX2Ny5c9E0jSNHjlyS7b/11ltomnZFtOVKsHbtWjRN4+uvvy7ppoizOHLkCJqmMXfu3GLVv1b6fyGEEEJcfhJEF0IIIa4zBw8epH///lSrVg0fHx+CgoK46667mDJlChkZGSXdPFEMe/fupVevXlSvXp2PPvqI2bNnF1k3Jzha2NesWbPOa78Xsq0ZM2YUO9B1Pux2O1OmTOHmm28mKCiIkJAQ6tatS79+/di7d6+nXk7gN+fLx8eHG2+8kcGDBxMdHe2pd67gaa9evQgICLjox3Eleuedd/j+++9LuhnnZdOmTbz11lskJiYWq/78+fOZPHnyJW2TuLx+/vlnCZQLIYQQ4qIzl3QDhBBCCHH5LFmyhEceeQSbzUaPHj2oV68edrudX3/9lRdeeIE9e/acNSB7LcjIyMBsvrpvgdauXYthGEyZMoUaNWoUa52ZM2cWCP42btz4P+3/v2xrxowZlC5dml69ev2nfRbl4YcfZunSpTz++OM8/fTTOBwO9u7dy+LFi7nzzjupVauWV/3Ro0dTtWpVMjMz+fXXX5k5cyY///wzu3fvxs/P76K27Wry2muv8fLLL3uVvfPOO3Tp0oUHH3zQq/zJJ5+ka9eu2Gy2y9jC4tm0aROjRo2iV69exXpDY/78+ezevZthw4Zd8raJiy8yMpKMjAwsFoun7Oeff2b69OmFBtKvhf5fCCGEECVD7iCEEEKI68Thw4fp2rUrkZGRrF69mvLly3uWDRo0iAMHDrBkyZISbOGlYxgGdrsdHx8ffHx8Sro5FywmJgbgvNK4dOnShdKlS1+U/V/MbV2IrVu3snjxYt5++21GjhzptWzatGmFjkZu164dt956KwB9+/alVKlSTJw4kR9++IHHH3/8cjT7imQ2m4sdXDSZTJhMpkvcIiHOLeetkuK6Fvp/IYQQQpQMSecihBBCXCfGjRtHamoqn3zyiVcAPUeNGjV49tlnPd87nU7+7//+j+rVq2Oz2ahSpQojR44kKyvLa70qVarQsWNH1q5dy6233oqvry/169dn7dq1AHz77bfUr18fHx8fGjVqxPbt273Wz0mPcejQIdq0aYO/vz8RERGMHj0apZRX3QkTJnDnnXdSqlQpfH19adSoUaFpNzRNY/DgwXzxxRfUrVsXm83GsmXLPMvyjlBMSUlh2LBhVKlSBZvNRtmyZbnvvvv4888/vba5aNEiGjVqhK+vL6VLl+aJJ57g5MmThR7LyZMnefDBBwkICKBMmTI8//zzuFyuIj4ZbzNmzPC0OSIigkGDBnkFg6tUqcKbb74JQJkyZS5ajt/iHN9/VaVKFfbs2cO6des86VTuvvtuz/JDhw7xyCOPEBYWhp+fH3fccUexHugcPHgQgLvuuqvAMpPJRKlSpc65jXvvvRdwP2S6WCZMmICmaRw9erTAsldeeQWr1UpCQoKnbMuWLbRt25bg4GD8/Pxo0aIFGzduLNa+zvXzkncf7du3JzQ0FH9/fxo0aMCUKVM8y/PnRNc0jbS0NObNm+f5zHLeIigqJ/rSpUtp1qwZ/v7+BAYG0qFDB/bs2eNVJyoqit69e1OxYkVsNhvly5fngQceOGd+9Z07d9KrVy9PGqrw8HCeeuop4uLivI7hhRdeAKBq1aqedhe17bvvvpslS5Zw9OhRT90qVap41TEMg7fffpuKFSvi4+NDy5YtOXDgQKHn979+hlOnTqVu3br4+fkRGhrKrbfeyvz58z3Ljx49ysCBA6lZsya+vr6UKlWKRx55pMBx5Xwuv/76K0OHDqVMmTKEhITQv39/7HY7iYmJ9OjRg9DQUEJDQ3nxxRcL9LGGYTB58mTq1q2Lj48P5cqVo3///l4/r0U5n748LS2NESNGUKlSJWw2GzVr1mTChAkF6q1cuZKmTZsSEhJCQEAANWvW9Hpglj8neq9evZg+fTqAV/qmHIX1l9u3b6ddu3YEBQUREBBAy5Yt+e233wo9txs3bmT48OGUKVMGf39/OnfuTGxs7DnPjRBCCCGufjISXQghhLhO/PTTT1SrVo0777yzWPX79u3LvHnz6NKlCyNGjGDLli2MHTuWf/75h++++86r7oEDB+jWrRv9+/fniSeeYMKECXTq1IlZs2YxcuRIBg4cCMDYsWN59NFH2bdvH7qe+yzf5XLRtm1b7rjjDsaNG8eyZct48803cTqdjB492lNvypQp3H///XTv3h273c6CBQt45JFHWLx4MR06dPBq0+rVq/nqq68YPHgwpUuXLhAcy/HMM8/w9ddfM3jwYOrUqUNcXBy//vor//zzD7fccgvgDqD07t2b2267jbFjxxIdHc2UKVPYuHEj27dv9xoR7nK5aNOmDY0bN2bChAmsWrWK999/n+rVqzNgwICznvO33nqLUaNG0apVKwYMGMC+ffuYOXMmW7duZePGjVgsFiZPnsxnn33Gd99950mr0qBBg3N+nvHx8V7fm0wmQkNDz/v4zrWtwkyePJkhQ4YQEBDAq6++CkC5cuUAiI6O5s477yQ9PZ2hQ4dSqlQp5s2bx/3338/XX39N586di9xuZGQkAF988QV33XXXf0rTkBOIzx9wT0lJ4cyZMwXq53+IVJhHH32UF198ka+++soT1M3x1Vdf0bp1a8/5Wr16Ne3ataNRo0a8+eab6LrOnDlzuPfee9mwYQO33357kfspzs8LuAORHTt2pHz58jz77LOEh4fzzz//sHjxYq8HZ3l9/vnn9O3bl9tvv51+/foBUL169SLb8vnnn9OzZ0/atGnDe++9R3p6OjNnzqRp06Zs377d8/v38MMPs2fPHoYMGUKVKlWIiYlh5cqVHDt2rMjf0ZxjOHToEL179yY8PNyTemrPnj389ttvaJrGQw89xL///suXX37JpEmTPG9LlClTptBtvvrqqyQlJXHixAkmTZoEUCBN0bvvvouu6zz//PMkJSUxbtw4unfvzpYtWzx1LuQz/Oijjxg6dChdunTh2WefJTMzk507d7Jlyxa6desGuN+42LRpE127dqVixYocOXKEmTNncvfdd/P3338XSEM0ZMgQwsPDGTVqFL/99huzZ88mJCSETZs2UblyZd555x1+/vlnxo8fT7169ejRo4dn3f79+3v6g6FDh3L48GGmTZvG9u3bvX6milKcvlwpxf3338+aNWvo06cPN910E8uXL+eFF17g5MmTns9iz549dOzYkQYNGjB69GhsNhsHDhw468OJ/v37c+rUKVauXMnnn39+1rbm7KNZs2YEBQXx4osvYrFY+PDDD7n77rtZt25dgTRVQ4YMITQ0lDfffJMjR44wefJkBg8ezMKFC8+5LyGEEEJc5ZQQQgghrnlJSUkKUA888ECx6u/YsUMBqm/fvl7lzz//vALU6tWrPWWRkZEKUJs2bfKULV++XAHK19dXHT161FP+4YcfKkCtWbPGU9azZ08FqCFDhnjKDMNQHTp0UFarVcXGxnrK09PTvdpjt9tVvXr11L3
      "text/plain": [
       "<Figure size 1500x400 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 4, figsize=(15, 4))\n",
    "fig.suptitle(\"Comparison of FoF to SPH velocities at the same position\")\n",
    "\n",
    "axs[0].hexbin(vel[:, 0], vx, gridsize=50, bins=\"log\", mincnt=1)\n",
    "axs[0].set_xlabel(r\"FoF $v_x ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[0].set_ylabel(r\"SPH $v_x ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "axs[1].hexbin(vel[:, 1], vy, gridsize=50, bins=\"log\", mincnt=1)\n",
    "axs[1].set_xlabel(r\"FoF $v_y ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[1].set_ylabel(r\"SPH $v_y ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "axs[2].hexbin(vel[:, 2], vz, gridsize=50, bins=\"log\", mincnt=1)\n",
    "axs[2].set_xlabel(r\"FoF $v_z ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[2].set_ylabel(r\"SPH $v_z ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "\n",
    "vr_fof = project_radial(vel[:, 0], vel[:, 1], vel[:, 2], RA, dec)\n",
    "vr_sph = project_radial(vx, vy, vz, RA, dec)\n",
    "axs[3].hexbin(vr_fof, vr_sph, gridsize=50, bins=\"log\", mincnt=1)\n",
    "axs[3].set_xlabel(r\"FoF $v_r ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[3].set_ylabel(r\"SPH $v_r ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "for i in range(4):\n",
    "    axs[i].axline([0, 0], [1, 1], color=\"red\", ls=\"--\")\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"../plots/fof_to_sph_velocity_comparison.png\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlation of the peculiar velocity and total mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "kind = \"main\"\n",
    "nsnap = 99\n",
    "nsim = 17417"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FoF haloes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SignificanceResult(statistic=-0.011239061725954754, pvalue=1.637324987020833e-17)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "catalogue = csiborgtools.read.CSiBORG2Catalogue(nsim, nsnap, kind)\n",
    "\n",
    "vel = catalogue[\"cartesian_vel\"]\n",
    "mass = catalogue[\"totmass\"]\n",
    "velmag = np.linalg.norm(vel, axis=1)\n",
    "\n",
    "spearmanr(mass, velmag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAHLCAYAAADcEAx6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydd3hTVRvAf/cmaboXbSlQ9p4yZO+9ZQ9RlqCCoIK4B+LCLbhwC4iiMmTvjciQPQXZs4XumaZJ7vn+SJM2bVoKFAqf5/c8fWzuPefc99xcuW/fqQghBBKJRCKRSCQSF9SiFkAikUgkEonkbkQqSRKJRCKRSCRukEqSRCKRSCQSiRukkiSRSCQSiUTiBqkkSSQSiUQikbhBKkkSiUQikUgkbpBKkkQikUgkEokbpJIkkUgkEolE4gapJEkkEolEIpG4QSpJEsn/AVOmTEFRlNu2/ogRIyhXrtxtW/9ewWq18vzzz1O6dGlUVaV3795FLZJbZs2ahaIonDt3rqhFyZcPP/yQChUqoNPpqFu3blGLI5HkQipJkv8UjpeHu58XX3yxyNa610hLS2PKlCls3ry5qEW5o/z44498+OGH9O/fn9mzZzNx4sQilWfq1KksXry4SGW4WdauXcvzzz9P8+bNmTlzJlOnTi1qkSSSXOiLWgCJpCh48803KV++vMuxWrVqFfladyvfffcdmqY5P6elpfHGG28A0KZNmyKS6s6zceNGSpUqxbRp04paFMCuJPXv3z+XRWvo0KEMHjwYo9FYNIIVgI0bN6KqKj/88AMeHh5FLY5E4hapJEn+k3Tt2pX777//rlvrbsVgMBS1CHcF165dIzAwsKjFuC46nQ6dTlfUYuTLtWvX8PLykgqS5K5GutskEjds3LiRli1b4uPjQ2BgIL169eKff/4plLXHjx+Pr68vaWlpuc49+OCDhIeHY7PZnMdWrVrllMXPz4/u3btz9OjR617HarXy1ltvUbFiRYxGI+XKlePll1/GbDbnGrtq1Spat26Nn58f/v7+NGzYkLlz5zrPZ49JOnfuHKGhoQC88cYbThfjlClTmDlzJoqisH///lzXmDp1KjqdjsuXL+cp8/nz53niiSeoWrUqXl5eFCtWjAEDBuSKrbFYLLzxxhtUrlwZT09PihUrRosWLVi3bl2+9yQuLo5nn32W2rVr4+vri7+/P127duXgwYP5zjt37hyKorBp0yaOHj3q3PPmzZvZvHmz83d3c2bNmuU8NmLECHx9fbl8+TK9e/fG19eX0NBQnn32WZfvHEDTND799FNq166Np6cnoaGhdOnShT179gCgKAqpqanMnj3bKc+IESOAvGOSZsyYQc2aNTEajZQsWZJx48aRkJDgMqZNmzbUqlWLY8eO0bZtW7y9vSlVqhQffPBBvvfIQUGeO0VRmDlzJqmpqU7Zs9+nnDhkOnToEK1bt8bb25tKlSqxYMECALZs2ULjxo3x8vKiatWqrF+/3mV+YT5XUVFRjBw5koiICIxGIyVKlKBXr153ffyX5OaQSpLkP0liYiIxMTEuPw7Wr19P586duXbtGlOmTOGZZ55h+/btNG/e3O0/hPmt5Y5BgwaRmprKihUrXI6npaWxbNky+vfv77QCzJkzh+7du+Pr68v777/Pa6+9xrFjx2jRosV1/1EePXo0kydPpn79+kybNo3WrVvz7rvvMnjwYJdxs2bNonv37sTFxfHSSy/x3nvvUbduXVavXu123dDQUL766isA+vTpw5w5c5gzZw59+/alf//+eHl58csvv+Sa98svv9CmTRtKlSqVp8y7d+9m+/btDB48mM8++4wxY8awYcMG2rRp46JUTpkyhTfeeIO2bdvyxRdf8Morr1CmTBn27duX7z05c+YMixcvpkePHnzyySc899xzHD58mNatW3PlypU854WGhjJnzhyqVatGRESEc8/Vq1fP93rusNlsdO7cmWLFivHRRx/RunVrPv74Y7799luXcaNGjWLChAmULl2a999/nxdffBFPT0927twJ2J8No9FIy5YtnfI8/vjjeV53ypQpjBs3jpIlS/Lxxx/Tr18/vvnmGzp16oTFYnEZGx8fT5cuXbjvvvv4+OOPqVatGi+88AKrVq267v4K8tzNmTOHli1bYjQanbK3atUq33Xj4+Pp0aMHjRs35oMPPsBoNDJ48GB+//13Bg8eTLdu3XjvvfdITU2lf//+JCcnO+cW5nPVr18/Fi1axMiRI5kxYwZPPfUUycnJXLhw4br3RnIPIiSS/xAzZ84UgNsfB3Xr1hVhYWEiNjbWeezgwYNCVVUxbNiwG1rLHZqmiVKlSol+/fq5HJ83b54AxNatW4UQQiQnJ4vAwEDx6KOPuoyLiooSAQEBLsdff/11l+seOHBAAGL06NEuc5999lkBiI0bNwohhEhISBB+fn6icePGwmQy5ZLTwfDhw0XZsmWdn6OjowUgXn/99Vz7e/DBB0XJkiWFzWZzHtu3b58AxMyZM/O5M0KkpaXlOrZjxw4BiJ9++sl57L777hPdu3fPdy13pKenu8glhBBnz54VRqNRvPnmm9ed37p1a1GzZk2XY5s2bRKA2LRpU651c+55+PDhAsh1rXr16okGDRo4P2/cuFEA4qmnnsolQ/bvxcfHRwwfPjzXGMezefbsWSGEENeuXRMeHh6iU6dOLvv/4osvBCB+/PFHlz3mvN9ms1mEh4fnemZzUtDnznEvfHx88l0vp0xz5851Hjt+/LgAhKqqYufOnc7ja9asyXXfC+u5io+PF4D48MMPCyS35N5HWpIk/0m+/PJL1q1b5/IDEBkZyYEDBxgxYgTBwcHO8XXq1KFjx46sXLmywGvlhaIoDBgwgJUrV5KSkuI8/vvvv1OqVClatGgBwLp160hISODBBx90sVLpdDoaN27Mpk2b8ryGQ85nnnnG5fikSZMAnFasdevWkZyc7LRS5JTzZhg2bBhXrlxxke+XX37By8uLfv365TvXy8vL+bvFYiE2NpZKlSoRGBjo8td8YGAgR48e5eTJkzckm9FoRFXt/+zZbDZiY2Px9fWlatWq17VCFSZjxoxx+dyyZUvOnDnj/Lxw4UIUReH111/PNfdmvpf169eTkZHBhAkTnPsHePTRR/H3989l1fT19eXhhx92fvbw8KBRo0YuMrqjoM/dzeDr6+tijapatSqBgYFUr16dxo0bO487fs8ua2E9V44Yqs2bNxMfH3/Te5HcO0glSfKfpFGjRnTo0MHlB+yxC2D/Bzgn1atXJyYmhtTU1AKtlR+DBg3CZDKxdOlSAFJSUli5ciUDBgxwvgQd/1C3a9eO0NBQl5+1a9dy7dq1PNc/f/48qqpSqVIll+Ph4eEEBgY693n69GmgcLPxOnbsSIkSJZwuN03T+PXXX+nVqxd+fn75zjWZTEyePJnSpUtjNBoJCQkhNDSUhIQEEhMTnePefPNNEhISqFKlCrVr1+a5557j0KFD15VN0zSmTZtG5cqVXdY/dOiQy/q3E0d8UXaCgoJcXrqnT5+mZMmSLor6rZDXc+3h4UGFChWc5x1ERETkUsZyypjXdQry3N0M7mQKCAigdOnSuY4BLrIW1nNlNBp5//33WbVqFcWLF6dVq1Z88MEHREVF3fS+JHc3UkmSSIqAJk2aUK5cOebNmwfAsmXLMJlMDBo0yDnGkXI/Z86cXJaqdevWsWTJkute53YWmMwLnU7HkCFDWLhwIenp6WzatIkrV664WCby4sknn+Sdd95h4MCBzJs3j7Vr17Ju3TqKFSvmUoKgVatWnD59mh9//JFatWrx/fffU79+fb7//vt81586dSrPPPMMrVq14ueff2bNmjWsW7eOmjVruqx/I+R1j3MGYju427POIG8ZhRAFmn87nru8ZCqIrIX5XE2YMIF///2Xd999F09PT1577TWqV6/uNllBcu8jSwBIJNkoW7YsACdOnMh17vjx44SEhODj41Mo1xo4cCCffvopSUlJ/P7775QrV44mTZo4z1esWBGAsLCwAlmnslO2bFk0TePkyZMuwcVXr14lISHBuU/HNY4cOZLrr//
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.title(\"FoF velocity as a function of mass\")\n",
    "plt.hexbin(np.log10(mass), velmag, mincnt=1, bins=\"log\")\n",
    "\n",
    "\n",
    "plt.xlim(np.log10(np.min(mass)))\n",
    "plt.ylim(0)\n",
    "plt.xlabel(r\"$\\log M_{\\rm FoF} ~ [M_\\odot / h]$\")\n",
    "plt.ylabel(r\"$V_{\\rm pec} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "plt.savefig(\"../plots/fof_velocity_to_mass.png\", dpi=300, bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subfind subhaloes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalogue = csiborgtools.read.CSiBORG2SUBFINDCatalogue(nsim, nsnap, kind)\n",
    "\n",
    "mass = catalogue[\"totmass\"]\n",
    "vec = np.linalg.norm(catalogue[\"cartesian_vel\"], axis=1)\n",
    "is_central = catalogue[\"Central\"]\n",
    "\n",
    "def get_binned_trend(x, y, nbins=20):\n",
    "    median, bin_edges, __ = binned_statistic(x, y, bins=nbins, statistic=\"median\")\n",
    "    lower, __, __ = binned_statistic(x, y, bins=nbins, statistic=lambda x: np.percentile(x, 16))\n",
    "    upper, __, __ = binned_statistic(x, y, bins=nbins, statistic=lambda x: np.percentile(x, 84))\n",
    "    std, __, __ = binned_statistic(x, y, bins=nbins, statistic=\"std\")\n",
    "    xrange = (bin_edges[1:] + bin_edges[:-1]) / 2\n",
    "    std = (upper - lower) / 2\n",
    "    return xrange, median, lower, upper, std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAGGCAYAAABmGOKbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd1gUxxvA8e/e0aUJiIBgr9hil9h77zUae4k9ahKNiYlRf5bEGDXNaEw0RU2ssWDvDXvvvSJgpde7/f1xcHoCCngI6vt5Hp+Evd3Z2aW8Ozsz7yiqqqoIIYQQQgghhBDC7DRZXQEhhBBCCCGEEOJNJY1uIYQQQgghhBAik0ijWwghhBBCCCGEyCTS6BZCCCGEEEIIITKJNLqFEEIIIYQQQohMIo1uIYQQQgghhBAik0ijWwghhBBCCCGEyCTS6BZCCCGEEEIIITKJNLqFEEIIIYQQQohMIo1uIYQQQgghhBAik0ijWwghhBBCCCGEyCTS6BYijRYsWICiKFy/fj2rqyJEmnz11VcoisL9+/ezuipCiLdASnHy2W0SS4VIncTtN5c0ukWmuHLlCh988AEFCxbExsYGR0dHqlWrxqxZs4iOjs608+7bt4+vvvqKx48fZ9o53kYRERGMGzeOxo0b4+LigqIoLFiwIMV9jxw5QuPGjXF0dMTBwYGGDRty/Pjxl973WTt27EBRlBT/7d+/P2MXmkZnzpyhQ4cOFCxYEDs7O9zc3KhZsyZr1qxJtu/LXCPApUuX6Ny5M97e3tjZ2VG8eHEmTJhAVFSUGa/o5en1enLlysU333yT1VUR4q1x6tQp2rdvT758+bCxsSFPnjw0aNCAH374Id1lZWX8lNidNumJPYcOHWLIkCGULFmSHDlykDdvXjp27MjFixefe45JkyahKAqlSpVKc71eNs5lVHqeTQBiY2MZPXo0Xl5e2NraUqVKFTZv3myyjzmeLbLqfmSExO5XxyKrKyDePP7+/nTo0AFra2u6d+9OqVKliIuLY8+ePXzyySecOXOGuXPnZsq59+3bx/jx4+nZsyfOzs6Zco630f3795kwYQJ58+albNmy7NixI8X9jh49SvXq1fHx8WHcuHHo9Xp+/vlnatWqxcGDBylWrFiG9n2eYcOGUalSJZNthQsXzvC1psWNGzcIDw+nR48eeHl5ERUVxfLly2nZsiVz5syhf//+wMtf461bt6hcuTJOTk4MGTIEFxcXAgICGDduHEeOHGHVqlWZep3pcfDgQe7fv0+zZs2yuipCvBX27dtHnTp1yJs3L/369cPDw4Nbt26xf/9+Zs2axdChQ9Nd3quIn926daNz585YW1u/8nO/7tIaewC+/vpr9u7dS4cOHShTpgxBQUH8+OOPlC9fnv3796fYqL59+zaTJ08mR44caa6TuWJ5RqT12SRJz549WbZsGcOHD6dIkSIsWLCApk2bsn37dqpXr26yb0afLbLyfmSExO5XSBXCjK5evara29urxYsXVwMDA5N9funSJXXmzJmZdv5p06apgHrt2rUX7hsREZGusufPn5/mst80MTEx6t27d1VVVdVDhw6pgDp//vxk+zVt2lTNmTOnev/+feO2wMBA1d7eXm3btm2G903J9u3bVUBdunRpBq/KoFatWmqPHj1eqgxVVdWEhAS1bNmyarFixYzbXvYaJ02apALq6dOnTbZ3795dBdSHDx8+9/hx48apgHrv3r10Xk36ffHFF2q+fPky/TxCCIOmTZuquXLlUh89epTss+Dg4HSXl574mZqU4mRaYqc5zv22Sin2qKqq7t27V42NjTXZdvHiRdXa2lrt2rVrimV16tRJrVu3rlqrVi21ZMmSaTr/y8a5p6U3Hqf12URVVfXAgQMqoE6bNs24LTo6Wi1UqJDq5+dn3PayzxYvez9eZdxWVYndr5IMLxdm9c033xAREcFvv/2Gp6dnss8LFy7Mhx9+aPz6zp079O7dm9y5c2NtbU3JkiX5/fffTY5Jmt9y+fJl41twJycnevXqZTLE9quvvuKTTz4BoECBAsbhQNevXzeWcfbsWbp06ULOnDmNbzVv3LjBoEGDKFasGLa2tri6utKhQ4c0zTcLDw9n+PDh5M+fH2tra9zd3WnQoAFHjx7NyO0DYO7cuZQvXx47O7tkQ5sKFiyY4XJfhrW1NR4eHi/cb/fu3dSvXx9XV1fjNk9PT2rVqsXatWuJiIjI0L4vEh4eTkJCQpr3zwxarRYfHx+T4ZEve41hYWEA5M6d22S7p6cnGo0GKyurdNfzxo0bFC5cmFKlShEcHGz83bh48SLvv/8+Tk5O5MqViy+++AJVVbl16xatWrXC0dERDw8Ppk+fnmK5/v7+Jm/KM+N3QwjxxJUrVyhZsmSKPcPu7u7G/09LjHte/EySlnidFs/O6X7RudN63rclHj8rpdgD8O677yaLEUWKFKFkyZKcO3cuWTm7du1i2bJlzJw5M13nN2csT6+0PpsALFu2DK1WazIawMbGhj59+hAQEMCtW7eSHZORZ4vMuB/Pxm1AYvdrSIaXC7Nas2YNBQsW5N13333hvsHBwVStWhVFURgyZAi5cuVi/fr19OnTh7CwMIYPH26yf8eOHSlQoABTpkzh6NGjzJs3D3d3d77++msA2rZty8WLF1m8eDEzZszAzc0NgFy5chnL6NChA0WKFGHy5MmoqgoY5j3t27fPOG/2+vXrzJ49m9q1a3P27Fns7OxSvYYBAwawbNkyhgwZgq+vLw8ePGDPnj2cO3eO8uXLp/f2MWLECGbOnEnDhg3p1asXt2/fZsaMGcTHx9O8eXMqVKiQrvLi4+MJDQ1N074uLi5oNC/3Hi42NhZbW9tk2+3s7IiLi+P06dNUrVo13fs+T69evYiIiECr1VKjRg2mTZtGxYoVU9w3pfsRHx9PbGxssqQlabkfkZGRREdHExoayurVq1m/fj2dOnUyfv6y11i7dm2+/vpr+vTpw/jx43F1dWXfvn3Mnj2bYcOGpWsIIBge0uvWrYuLiwubN282/o4AdOrUiRIlSjB16lT8/f353//+h4uLC3PmzKFu3bp8/fXXLFy4kI8//phKlSpRs2ZN47FBQUEcO3aMCRMmGLeZ+3dDCGEqX758BAQEcPr06efOv01LjHtR/ExvvE6P5507Ped9m+Lxi2JPalRVJTg4mJIlS5ps1+l0DB06lL59+1K6dOk01TFJRuOcuePxixw7doyiRYvi6Ohosr1y5coAHD9+HB8fH+P29DxbPM1czzZJnhe3QWL3ayVrO9rFmyQ0NFQF1FatWqVp/z59+qienp4mQ3BUVVU7d+6sOjk5qVFRUaqqPhlq07t3b5P92rRpo7q6uppsS22IWlIZ7733XrJ6JJ3naQEBASqg/vnnn8ZtKQ2Rc3JyUgcPHpym632RXbt2qYA6cOBAk+3jx49XAfXgwYPpLjNpmFRa/qV1WN/zhnCVLl1aLVq0qJqQkGDcFhsbq+bNm1cF1GXLlmVo35Ts3btXbdeunfrbb7+pq1atUqdMmaK6urqqNjY26tGjR1/J/fjggw+M+2s0GrV9+/YmQ75f9hpVVVUnTpyo2tramtTt888/f+Fxqmo6TO3cuXOql5eXWqlSJZM6Ju3Tv39/47aEhATV29tbVRRFnTp1qnH7o0ePVFtb22TD/3777TfV1tbW5HfJnL8bQojkNm3apGq1WlWr1ap+fn7qqFGj1I0bN6pxcXEm+6U1xj1viHda43VahpentE9q507reVX17YrHL4o9qfnrr79UQP3tt99Mtv/444+qk5OTGhISoqqqmq7h5RmNc+aOxy8aXl6yZEm1bt26ybafOXNGBdRffvlFVdWMPVuY434kSUvcfno/id2vD+npFmaTNBTWwcHhhfuqqsry5cvp2LEjqqqavNVs1KgR//zzD0ePHqVatWrG7QMGDDApo0aNGqxcuZKwsLBkby5T82wZgMkbyfj4eMLCwihcuDDOzs4cPXqUbt26pVqes7MzBw4cIDAwEC8vrzTVITUzZszAxcWFadOmmWyvVasWABcvXkyW1ONFypYtmywzZ2rSOkT
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbins = 15\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(10, 4))\n",
    "\n",
    "# Centrals\n",
    "mask = is_central & (mass > 1e12)\n",
    "x, y = np.log10(mass[mask]), vec[mask]\n",
    "xrange, median, lower, upper, std = get_binned_trend(x, y, nbins=nbins)\n",
    "sigma_v, e_sigma_v = np.mean(std), np.std(std)\n",
    "axs[0].set_title(r\"Centrals, $\\sigma_v = {:.1f} \\pm {:.1f} ~ \\mathrm{{km}} / \\mathrm{{s}}$\".format(sigma_v, e_sigma_v))\n",
    "axs[0].hexbin(x, y, mincnt=1, bins=\"log\")\n",
    "\n",
    "axs[0].plot(xrange, median, color=\"red\")\n",
    "axs[0].fill_between(xrange, lower, upper, color=\"red\", alpha=0.3)\n",
    "\n",
    "# Satellites\n",
    "mask = ~is_central & (mass > 1e12)\n",
    "x, y = np.log10(mass[mask]), vec[mask]\n",
    "xrange, median, lower, upper, std = get_binned_trend(x, y, nbins=nbins)\n",
    "sigma_v, e_sigma_v = np.mean(std), np.std(std)\n",
    "axs[1].set_title(r\"Satellites, $\\sigma_v = {:.1f} \\pm {:.1f} ~ \\mathrm{{km}} / \\mathrm{{s}}$\".format(sigma_v, e_sigma_v))\n",
    "axs[1].hexbin(x, y, mincnt=1, bins=\"log\")\n",
    "axs[1].plot(xrange, median, color=\"red\")\n",
    "axs[1].fill_between(xrange, lower, upper, color=\"red\", alpha=0.3)\n",
    "\n",
    "for i in range(2):\n",
    "    axs[i].set_xlim(np.log10(1e12))\n",
    "    axs[i].set_ylim(0)\n",
    "    axs[i].set_xlabel(r\"$\\log M_{\\rm tot} ~ [M_\\odot / h]$\")\n",
    "    axs[i].set_ylabel(r\"$V_{\\rm pec} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"../plots/velocity_to_mass_subhaloes.png\", dpi=300, bbox_inches=\"tight\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SUBFIND haloes but subtracting the velocity of the velocity field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "field_reader = csiborgtools.read.CSiBORG2Field(nsim, kind)\n",
    "catalogue = csiborgtools.read.CSiBORG2SUBFINDCatalogue(nsim, nsnap, kind)\n",
    "boxsize = csiborgtools.simname2boxsize(\"csiborg2_main\")\n",
    "\n",
    "mass = catalogue[\"totmass\"]\n",
    "pos = catalogue[\"cartesian_pos\"] / boxsize\n",
    "subfind_vel = catalogue[\"cartesian_vel\"]\n",
    "is_central = catalogue[\"Central\"]\n",
    "\n",
    "parent_mass = catalogue[\"ParentMass\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity_field = field_reader.velocity_field(\"SPH\", 1024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "vx, vy, vz = csiborgtools.field.evaluate_cartesian_cic(velocity_field[0], velocity_field[1], velocity_field[2], pos=pos)\n",
    "\n",
    "residual_vel = np.copy(subfind_vel)\n",
    "residual_vel[:, 0] -= vx\n",
    "residual_vel[:, 1] -= vy\n",
    "residual_vel[:, 2] -= vz\n",
    "\n",
    "residual_vel_mag = np.linalg.norm(residual_vel, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAGMCAYAAADHg8H9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzddXgUVxfA4d9u3EMgikNwd4J7gOAuheAOBVqgtHhLKVRwKcVKkRZooTgELyVIcXcrkKBJCPHd+f7Ily1LbBN2SQLnfZ6U7sydO2cmmz17Z+7cq1IURUEIIYQQQgghhBBGp87oAIQQQgghhBBCiPeVNLqFEEIIIYQQQggTkUa3EEIIIYQQQghhItLoFkIIIYQQQgghTEQa3UIIIYQQQgghhIlIo1sIIYQQQgghhDARaXQLIYQQQgghhBAmIo1uIYQQQgghhBDCRKTRLYQQQgghhBBCmIg0uoUQ741JkyahUqkMKqtSqZg0aZJJ46lTpw516tQx6T4yYr/58uWjR48eJqtfvF/u3LmDSqVixYoVqZa9fv06jRo1wsnJCZVKxaZNm0wW14EDB1CpVGzYsMFk+0hJwnn57rvvMmT/mcUvv/xC0aJFsbCwwNnZGXi7zzBDP59WrFiBSqXizp076dqPEEKkhTS6hRBGl/BlJuHH3NycnDlz0qNHDx48eJDR4Qkju3TpEpMmTcr0X14jIiKYNGkSBw4cyOhQMqUFCxYY1DA2JX9/f86fP8/UqVP55ZdfqFixYobGYwzbt283+QW+1ISHhzNx4kRKliyJnZ0d2bNnp2zZsnz88cc8fPhQVy7hwmXCj62tLcWLF2fcuHGEhYXpyiV8xv/zzz9J7q9OnTqULFky1biuXLlCjx49KFiwID/99BOLFy9++4MVQohMyDyjAxBCvL+mTJlC/vz5iYqK4ujRo6xYsYLDhw9z4cIFrK2tjb6/cePG8dlnnxm9XqHv6tWrqNX/XbO9dOkSkydPpk6dOuTLly/jAktFREQEkydPBsiQHgiZ3YIFC8iRI0eG9WKIjIwkMDCQL774giFDhmRIDKawfft25s+fn2EN79jYWGrVqsWVK1fw9/dn6NChhIeHc/HiRdasWUPr1q3x8vLS22bhwoXY29sTHh7O7t27mTp1Kvv27ePvv/82uDeRIQ4cOIBWq2X27Nl4e3vrlu/evdto+xBCiMxAGt1CCJNp0qSJ7k5Vnz59yJEjB9OnT2fz5s106NDB6PszNzfH3Fw+1kzNysoqo0MQ76EnT54A6LoYG8OrV6+ws7MzWn1Z0aZNmzh9+jSrV6+mS5cueuuioqKIiYlJtE27du3IkSMHAAMGDKBt27b88ccfHD16FB8fH6PF9vjxYyDx79zS0tJo+xBCiMxAupcLId6ZmjVrAnDz5k295VeuXKFdu3a4uLhgbW1NxYoV2bx5s16Z2NhYJk+eTKFChbC2tiZ79uzUqFGDgIAAXZmknumOjo5mxIgRuLq64uDgQIsWLfj3338TxdajR48k79ImVefy5cupV68ebm5uWFlZUbx4cRYuXJimc5GgZMmS1K1bN9FyrVZLzpw5adeund6yWbNmUaJECaytrXF3d6d///68ePEi1f08fvyY3r174+7ujrW1NWXKlOHnn39Ocr+zZ8+mVKlSWFtb4+rqSuPGjfW6kb7+zOSKFSto3749AHXr1tV1Sz1w4AD+/v7kyJGD2NjYRPtp1KgRRYoUSTHmhC6qly5dom7dutja2pIzZ05mzJiR5uO7c+cOrq6uAEyePFkXZ0p3H58/f86nn35KqVKlsLe3x9HRkSZNmnD27NlEZaOiopg0aRKFCxfG2toaT09P2rRpo/deN+TcxsXF8eWXX1KwYEGsrKzIly8fn3/+OdHR0Xr7Sy72N59nTegG/PfffzNy5EhcXV2xs7OjdevWukZuwnYXL17k4MGDunPzem+AkJAQhg8fTu7cubGyssLb25vp06ej1Wr19h8SEkKPHj1wcnLC2dkZf39/QkJCkj3HCSZNmkTevHkBGDVqFCqVSu/v8fTp0zRp0gRHR0fs7e2pX78+R48e1asj4VgPHjzIoEGDcHNzI1euXKnuW6PR8Pnnn+Ph4YGdnR0tWrTg/v37uvUTJ07EwsJC73wl6NevH87OzkRFRSVZd48ePZg/fz6AXrftNy1evFj3O69UqRInTpxIVMaQz8mkJLwHq1evnmidtbU1jo6OqdZRr149AG7fvp1qWUPly5ePiRMnAuDq6qr3nk7qme7o6GgmTpyIt7c3VlZW5M6dm9GjRyf620jKxYsXqVevHjY2NuTKlYuvvvoq0XtXCCFMSW4JCSHemYRnfrNly6ZbdvHiRapXr07OnDn57LPPsLOzY926dbRq1Yrff/+d1q1bA/FfyqdNm0afPn2oXLkyYWFh/PPPP5w6dYqGDRsmu88+ffqwatUqunTpQrVq1di3bx9+fn5vdRwLFy6kRIkStGjRAnNzc7Zs2cKgQYPQarUMHjw4TXV17NiRSZMmERQUhIeHh2754cOHefjwIZ06ddIt69+/PytWrKBnz54MGzaM27dvM2/ePE6fPs3ff/+NhYVFkvuIjIykTp063LhxgyFDhpA/f37Wr19Pjx49CAkJ4eOPP9aV7d27NytWrKBJkyb06dOHuLg4/vrrL44ePZrk87W1atVi2LBhzJkzh88//5xixYoBUKxYMbp168bKlSvZtWsXzZo1020TFBTEvn37dF+4U/LixQsaN25MmzZt6NChAxs2bGDMmDGUKlWKJk2aGHx8rq6uLFy4kIEDB9K6dWvatGkDQOnSpZPd961bt9i0aRPt27cnf/78BAcH8+OPP1K7dm0uXbqk65Kr0Who1qwZe/fupVOnTnz88ce8fPmSgIAALly4QMGCBQ0+t3369OHnn3+mXbt2fPLJJxw7doxp06Zx+fJlNm7cmOr5Ss7QoUPJli0bEydO5M6dO8yaNYshQ4bw22+/ATBr1iyGDh2Kvb09X3zxBQDu7u5AfLf82rVr8+DBA/r370+ePHk4cuQIY8eO5dGjR8yaNQsARVFo2bIlhw8fZsCAARQrVoyNGzfi7++fanxt2rTB2dmZESNG0LlzZ5o2bYq9vT0Q/xlRs2ZNHB0dGT16NBYWFvz444/UqVOHgwcPUqVKFb26Bg0ahKurKxMmTODVq1ep7nvq1KmoVCrGjBnD48ePmTVrFg0aNODMmTPY2NjQrVs3pkyZwm+//abX7T0mJoYNGzbQtm3bZB+X6d+/Pw8fPiQgIIBffvklyTJr1qzh5cuX9O/fH5VKxYwZM2jTpg23bt3S/U0b+jmZlISLGStXrmTcuHHp6h6e0HDPnj273vLQ0FCePn2aqHxSF9reNGvWLFauXMnGjRt13dmT+3vUarW0aNGCw4cP069fP4oVK8b58+eZOXMm165dS3HAvaCgIOrWrUtcXJzu3C1evBgbG5tUYxRCCKNRhBDCyJYvX64Ayp49e5QnT54o9+/fVzZs2KC4uroqVlZWyv3793Vl69evr5QqVUqJiorSLdNqtUq1atWUQoUK6ZaVKVNG8fPzS3G/EydOVF7/WDtz5owCKIMGDdIr16VLFwVQJk6cqFvm7++v5M2bN9U6FUVRIiIiEpXz9fVVChQooLesdu3aSu3atVOM+erVqwqgzJ07V2/5oEGDFHt7e92+/vrrLwVQVq9erVdu586diZa/ud9Zs2YpgLJq1SrdspiYGMXHx0ext7dXwsLCFEVRlH379imAMmzYsERxarVa3f/nzZtX8ff3171ev369Aij79+/X20aj0Si5cuVSOnbsqLf8hx9+UFQqlXLr1q0Uzkz8cQDKypUrdcuio6MVDw8PpW3btmk+vidPniT6vackKipK0Wg0estu376tWFlZKVOmTNEtW7ZsmQIoP/zwQ6I6Es6bIec24f3ap08fvfWffvqpAij79u3TLUvuON783ST8LTZo0EDvdzhixAjFzMxMCQkJ0S0rUaJEku/XL7/8UrGzs1OuXbumt/yzzz5TzMzMlHv37imKoiibNm1SAGXGjBm6MnFxcUrNmjUVQFm+fHm
      "text/plain": [
       "<Figure size 1000x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbins = 15\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, figsize=(10, 4))\n",
    "fig.suptitle(\"Residual velocity not accounted for by the SPH field.\")\n",
    "# Centrals\n",
    "mask = is_central & (mass > 1e12)\n",
    "x, y = np.log10(mass[mask]), residual_vel_mag[mask]\n",
    "xrange, median, lower, upper, std = get_binned_trend(x, y, nbins=nbins)\n",
    "sigma_v, e_sigma_v = np.mean(std), np.std(std)\n",
    "axs[0].set_title(r\"Centrals, $\\sigma_v = {:.1f} \\pm {:.1f} ~ \\mathrm{{km}} / \\mathrm{{s}}$\".format(sigma_v, e_sigma_v))\n",
    "axs[0].hexbin(x, y, mincnt=1, bins=\"log\")\n",
    "\n",
    "axs[0].plot(xrange, median, color=\"red\")\n",
    "axs[0].fill_between(xrange, lower, upper, color=\"red\", alpha=0.3)\n",
    "\n",
    "# Satellites\n",
    "mask = ~is_central & (mass > 1e12)\n",
    "x, y = np.log10(mass[mask]), residual_vel_mag[mask]\n",
    "xrange, median, lower, upper, std = get_binned_trend(x, y, nbins=nbins)\n",
    "sigma_v, e_sigma_v = np.mean(std), np.std(std)\n",
    "axs[1].set_title(r\"Satellites, $\\sigma_v = {:.1f} \\pm {:.1f} ~ \\mathrm{{km}} / \\mathrm{{s}}$\".format(sigma_v, e_sigma_v))\n",
    "axs[1].hexbin(x, y, mincnt=1, bins=\"log\")\n",
    "axs[1].plot(xrange, median, color=\"red\")\n",
    "axs[1].fill_between(xrange, lower, upper, color=\"red\", alpha=0.3)\n",
    "\n",
    "for i in range(2):\n",
    "    axs[i].set_xlim(np.log10(1e12))\n",
    "    axs[i].set_ylim(0)\n",
    "    axs[i].set_xlabel(r\"$\\log M_{\\rm tot} ~ [M_\\odot / h]$\")\n",
    "    axs[i].set_ylabel(r\"$\\Delta V_{\\rm pec} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"../plots/residual_velocity_to_mass_subhaloes.png\", dpi=300, bbox_inches=\"tight\")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArUAAAHbCAYAAADVtmL2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gUVRfA4d/spnfSCCQhlNBC74Tee++odFSqgoqKfgpW7IANRemCShEB6b333nuHJLR00nbv90fYJUsKCQTC6nmfJ49y9+7Mmdnd2bN3ztzRlFIKIYQQQgghrJgurwMQQgghhBDicUlSK4QQQgghrJ4ktUIIIYQQwupJUiuEEEIIIayeJLVCCCGEEMLqSVIrhBBCCCGsniS1QgghhBDC6klSK4QQQgghrJ4ktUIIIYQQwupZRVI7duxYNE3LVl9N0xg7duwTjadBgwY0aNDgia4jL9ZbuHBh+vbt+8SWL/5dLly4gKZpTJ8+/aF9T58+TbNmzXB3d0fTNP7+++8nFteGDRvQNI358+c/sXVkxbRfvvrqqzxZ/7Ni1qxZlCpVCltbWzw8PIDHO4Zl9/g0ffp0NE3jwoULj7SeZ0VsbCwDBw7Ez88PTdMYMWJEXoeUoZx8P+eVlJQU3nzzTQIDA9HpdHTo0CGvQxJPSI6SWtPBwvRnY2ODv78/ffv25erVq08qRpFHjh07xtixY5/5L4f4+HjGjh3Lhg0b8jqUZ9KPP/6YrcTzSerTpw+HDx/mk08+YdasWVStWjVP48kNy5Yte+I/oB8mNjaWMWPGULZsWZydnfHy8qJixYq8+uqrXLt2zdzPlHiY/pycnAgJCeF///sf0dHR5n6mY/yePXsyXF+DBg0oW7bsQ+M6ceIEffv2pVixYvzyyy9Mnjz58Tf2P+bTTz9l+vTpDB48mFmzZtGrV688i8Xaj7FTp07lyy+/pEuXLsyYMYORI0fmdUhPzKeffvpEBw2edTaP8qQPP/yQIkWKkJCQwI4dO5g+fTpbtmzhyJEjODg45HaM/O9//+Ptt9/O9eUKSydPnkSnu/8759ixY3zwwQc0aNCAwoUL511gDxEfH88HH3wAkCcj6M+6H3/8EW9v7zwbhb979y7bt2/n3XffZdiwYXkSw5OwbNkyfvjhhzxLbJOTk6lXrx4nTpygT58+DB8+nNjYWI4ePcqcOXPo2LEjBQsWtHjOpEmTcHFxITY2llWrVvHJJ5+wbt06tm7dmqujbRs2bMBoNDJx4kSCg4PN7atWrcq1dfzbrVu3jpo1azJmzJi8DiXLY6w1fD+vW7cOf39/xo8fn9ehPHGffvopXbp0+c+ORj9SUtuyZUvzSMvAgQPx9vbm888/Z/HixXTr1i1XAwSwsbHBxuaRQhU5YG9vn9chiH+hGzduAJhPQeeGuLg4nJ2dc2151ujvv/9m//79zJ49m+eee87isYSEBJKSktI9p0uXLnh7ewMwaNAgOnfuzF9//cWOHTsIDQ3NtdgiIiKA9K+5nZ1drq3j3y4iIoKQkJC8DuOhrOH7OSIiIlePP48rISEBOzs7i0EkkTtyZY/WrVsXgLNnz1q0nzhxgi5duuDp6YmDgwNVq1Zl8eLFFn2Sk5P54IMPKF68OA4ODnh5eVGnTh1Wr15t7pNRzU5iYiIjR47Ex8cHV1dX2rVrx5UrV9LF1rdv3wxHGTNa5rRp02jUqBG+vr7Y29sTEhLCpEmTcrQvTMqWLUvDhg3TtRuNRvz9/enSpYtF24QJEyhTpgwODg7kz5+fl19+mTt37jx0PREREQwYMID8+fPj4OBAhQoVmDFjRobrnThxIuXKlcPBwQEfHx9atGhhcZoxbc3a9OnT6dq1KwANGzY0n7bcsGEDffr0wdvbm+Tk5HTradasGSVLlswyZtMpzGPHjtGwYUOcnJzw9/fniy++yPH2XbhwAR8fHwA++OADc5xZjZ7dvn2bN954g3LlyuHi4oKbmxstW7bk4MGD6fomJCQwduxYSpQogYODAwUKFKBTp04W7/Xs7NuUlBQ++ugjihUrhr29PYULF+add94hMTHRYn2Zxf5gPaHpNPHWrVt57bXX8PHxwdnZmY4dO5qTSNPzjh49ysaNG837Ju1IS2RkJCNGjCAwMBB7e3uCg4P5/PPPMRqNFuuPjIykb9++uLu74+HhQZ8+fYiMjMx0H5uMHTuWoKAgAEaNGoWmaRafx/3799OyZUvc3NxwcXGhcePG7Nixw2IZpm3duHEjQ4YMwdfXl4CAgIeu22Aw8M477+Dn54ezszPt2rXj8uXL5sfHjBmDra2txf4yeemll/Dw8CAhISHDZfft25cffvgBwOK0/oMmT55sfs2rVavG7t270/XJznEyI6b3YO3atdM95uDggJub20OX0ahRIwDOnz//0L7ZVbhwYfPooo+Pj8V7OqOa2sTERMaMGUNwcDD29vYEBgby5ptvpvtsZOTo0aM0atQIR0dHAgIC+Pjjj9O9dzNz6NAh+vbtS9GiRXFwcMDPz4/+/ftz69Yti34xMTGMGDGCwoULY29vj6+vL02bNmXfvn1ZLv/ixYsMGTKEkiVL4ujoiJeXF127dn1oOZepJvz8+fMsXbrU/N66cOFCpvXCpuekLQ/IyXE2q+Pcw46xGX2XZvd4V7hwYdq0acOWLVuoXr06Dg4OFC1alJkzZ2a5j0zi4uJ4/fXXzcevkiVL8tVXX6GUAu7Xt69fv56jR49afI9lxhTTqlWrqFixIg4ODoSEhPDXX39Z9Mvu94jptfnjjz/43//+h7+/P05OTuayn507d9KiRQvc3d1xcnKifv36bN261WIZpn185swZ+vbti4eHB+7u7vTr14/4+HhzP03TiIuLY8aMGeZtzeoMnSm2uXPn8sEHH+Dv74+rqytdunQhKiqKxMRERowYga+vLy4uLvTr1y/da5jdvGnPnj00b94cb29vHB0dKVKkCP3797fo88cff1ClShVcXV1xc3OjXLlyTJw4MdP4M5IrP69MH7B8+fKZ244ePUrt2rXx9/fn7bffxtnZmblz59KhQwcWLFhAx44dgdQXa9y4cQwcOJDq1asTHR3Nnj172LdvH02bNs10nQMHDuS3337jueeeo1atWqxbt47WrVs/1nZMmjSJMmXK0K5dO2xsbFiyZAlDhgzBaDQydOjQHC2re/fujB07lrCwMPz8/MztW7Zs4dq1a/To0cPc9vLLLzN9+nT69evHK6+8wvnz5/n+++/Zv38/W7duxdbWNsN13L17lwYNGnDmzBmGDRtGkSJFmDdvHn379iUyMpJXX33V3HfAgAFMnz6dli1bMnDgQFJSUti8eTM7duzIsL6xXr16vPLKK3z77be88847lC5dGoDSpUvTq1cvZs6cycqVK2nTpo35OWFhYaxbty5bp8vu3LlDixYt6NSpE926dWP+/Pm89dZblCtXjpYtW2Z7+3x8fJg0aRKDBw+mY8eOdOrUCYDy5ctnuu5z587x999/07VrV4oUKUJ4eDg///wz9evX59ixY+ZTtgaDgTZt2rB27Vp69OjBq6++SkxMDKtXr+bIkSMUK1Ys2/t24MCBzJgxgy5duvD666+zc+dOxo0bx/Hjx1m4cOFD91dmhg8fTr58+RgzZgwXLlxgwoQJDBs2jD///BOACRMmMHz4cFxcXHj33XcByJ8/P5B6SrF+/fpcvXqVl19+mUKFCrFt2zZGjx7N9evXmTBhAgBKKdq3b8+WLVsYNGgQpUuXZuHChfTp0+eh8XXq1AkPDw9GjhxJz549adWqFS4uLkDqMaJu3bq4ubnx5ptvYmtry88//0yDBg3YuHEjNWrUsFjWkCFD8PHx4f333ycuLu6h6/7kk0/QNI233nqLiIgIJkyYQJMmTThw4ACOjo706tWLDz/8kD///NOiLCIpKYn58+fTuXPnTMupXn75Za5du8bq1auZNWtWhn3mzJlDTEwML7/8Mpqm8cUXX9CpUyfOnTtn/kxn9ziZEdOPhZkzZ/K///3vkcoHTImxl5eXRXtUVBQ3b95M1z+jH7IPmjBhAjNnzmThwoXmcofMPo9Go5F
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbins = 15\n",
    "\n",
    "fig, axs = plt.subplots(1, 1)\n",
    "axs = [axs]\n",
    "fig.suptitle(\"Residual velocity not accounted for by the SPH field as a function of parent mass\")\n",
    "\n",
    "# Satellites\n",
    "mask = ~is_central & (parent_mass > 1e12)\n",
    "x, y = np.log10(parent_mass[mask]), residual_vel_mag[mask]\n",
    "xrange, median, lower, upper, std = get_binned_trend(x, y, nbins=nbins)\n",
    "sigma_v, e_sigma_v = np.mean(std), np.std(std)\n",
    "axs[0].set_title(r\"Satellites, $\\sigma_v = {:.1f} \\pm {:.1f} ~ \\mathrm{{km}} / \\mathrm{{s}}$\".format(sigma_v, e_sigma_v))\n",
    "axs[0].hexbin(x, y, mincnt=1, bins=\"log\")\n",
    "axs[0].plot(xrange, median, color=\"red\")\n",
    "axs[0].fill_between(xrange, lower, upper, color=\"red\", alpha=0.3)\n",
    "\n",
    "axs[0].set_xlim(np.log10(1e12))\n",
    "axs[0].set_ylim(0)\n",
    "axs[0].set_xlabel(r\"$\\log M_{\\rm parent} ~ [M_\\odot / h]$\")\n",
    "axs[0].set_ylabel(r\"$\\Delta V_{\\rm pec} ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(\"../plots/residual_velocity_to_mass_subhaloes_parent.png\", dpi=300, bbox_inches=\"tight\")\n",
    "fig.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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}