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": null,
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity_field = field_reader.velocity_field(\"SPH\", 1024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kind = \"main\"\n",
    "nsnap = 99\n",
    "nsim = 17417"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FoF haloes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity_field = field_reader.velocity_field(\"SPH\", 1024)\n",
    "vx, vy, vz = csiborgtools.field.evaluate_cartesian_cic(velocity_field[0], velocity_field[1], velocity_field[2], pos=pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "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": "markdown",
   "metadata": {},
   "source": [
    "### Calibration of $\\sigma_{\\rm sat}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "kind = \"main\"\n",
    "nsnap = 99\n",
    "nsim = 17417\n",
    "boxsize = csiborgtools.simname2boxsize(\"csiborg2_main\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "field_reader = csiborgtools.read.CSiBORG2Field(nsim, kind)\n",
    "catalogue = csiborgtools.read.CSiBORG2SUBFINDCatalogue(nsim, nsnap, kind)\n",
    "\n",
    "mass = catalogue[\"totmass\"]\n",
    "parent_mass = catalogue[\"ParentMass\"]\n",
    "pos = catalogue[\"cartesian_pos\"] / boxsize\n",
    "subfind_vel = catalogue[\"cartesian_vel\"]\n",
    "is_central = catalogue[\"Central\"]\n",
    "\n",
    "velocity_field = field_reader.velocity_field(\"SPH\", 1024)\n",
    "vx, vy, vz = csiborgtools.field.evaluate_cartesian_cic(velocity_field[0], velocity_field[1], velocity_field[2], pos=pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "RA = np.deg2rad(catalogue[\"spherical_pos\"][:, 1])\n",
    "dec = np.deg2rad(catalogue[\"spherical_pos\"][:, 2])\n",
    "\n",
    "\n",
    "dvx = subfind_vel[:, 0] - vx\n",
    "dvy = subfind_vel[:, 1] - vy\n",
    "dvz = subfind_vel[:, 2] - vz\n",
    "\n",
    "dvx_projected = csiborgtools.flow.project_Vext(dvx, dvy, dvz, RA, dec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigma_vir(Mh):\n",
    "    \"\"\"Suvodip's sigma_vir at z = 0\"\"\"\n",
    "    Omega_m = 0.3111\n",
    "    gv = 0.9\n",
    "    Delta_nl = 18 * np.pi**2 - 60 * (1 - Omega_m) - 32 * (1 - Omega_m)**2\n",
    "    return 476 * gv * Delta_nl**(1 / 6) * (Mh / 1e15)**(1/3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.925519739350855 0.3163983444613365\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAMWCAYAAADs4eXxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdfXyT9b0//leS5qZp0yb0jt4BhdaCshZh0zFL2RRRjwqbMOdRj2zimWc4p0yc85wdb3bc2dHJmGdnOPcdTPebzk3QgTqx4g2lylAptIJQChR6S2+TNm3uk+v3R3JdvZImbdomhZbX8/HwEUmuXNeVNMl1va/P+/N+KwRBEEBEREREREREMac81ztARERERERENFUx6CYiIiIiIiKKEwbdRERERERERHHCoJuIiIiIiIgoThh0ExEREREREcUJg24iIiIiIiKiOGHQTURERERERBQnDLqJiIiIiIiI4oRBNxEREREREVGcMOgmIhqlb3/725g1a9a53g0iIqIhFAoFHnvsMenfzz//PBQKBU6fPi3d99WvfhVf/epXJ3zfiC5UDLqJKK5OnjyJu+++G7Nnz4ZOp0NKSgquuOIKPPPMM7Db7XHZZmtrKx577DEcOnQoLusnIiISffbZZ1i9ejVmzpwJnU6H3NxcXH311fj1r3896nVt3rwZzz//fOx3cgQ8bhLFV8K53gEimrrefPNNfPOb34RWq8Udd9yB+fPnw+VyoaqqCg8++CCOHDmC3/3udzHfbmtrKx5//HHMmjULCxYsiPn6iYiIAOCjjz7C1772NcyYMQP/+q//iunTp6OpqQn/+Mc/8Mwzz+Dee+8d1fo2b96M9PR0fPvb347PDgdUVFQE/ZvHTaL4YtBNRHHR0NCAW265BTNnzsR7772H7Oxs6bF77rkHJ06cwJtvvnkO93CQzWaDXq8/17tBRESTzM9+9jOkpqbik08+gdFoDHqso6Pj3OxUFDQazbneBaILCtPLiSgunnrqKfT392PLli1BAbeosLAQ9913n/TvP/3pT1i0aBESExMxbdo03HLLLWhqagp6zle/+lXMnz8fn3/+Ob72ta9Br9cjNzcXTz31lLTMBx98gC996UsAgO985ztQKBRQKBRSup64jgMHDqC8vBx6vR7//u//DgDYsWMHrr/+euTk5ECr1WLOnDn4r//6L3i93hFf78svv4xFixbBYDAgJSUFX/jCF/DMM8+M+n0jIqLJ4+TJk7jkkkuGBNwAkJmZKf3/H/7wB1x55ZXIzMyEVqvFxRdfjGeffTZo+VmzZuHIkSPYs2ePdOySz7u2WCy4//77kZ+fD61Wi8LCQjz55JPw+Xyj3m/5nO6RjpsAsH//flx77bVITU2FXq/H0qVL8eGHHwat02q14v7778esWbOg1WqRmZmJq6++GtXV1aPeP6KphiPdRBQXr7/+OmbPno2vfOUrIy77s5/9DP/5n/+Jm2++GXfddRc6Ozvx61//GuXl5Th48GDQyYzZbMa1116Lm266CTfffDO2bduGhx56CF/4whdw3XXXYd68efjpT3+KRx55BN/97nexZMkSAAjaj+7ublx33XW45ZZbcPvttyMrKwuAv9hMcnIyfvjDHyI5ORnvvfceHnnkEfT19eEXv/hFxP1/55138M///M+46qqr8OSTTwIAjh49ig8//DDowgIREU0tM2fOxL59+3D48GHMnz8/4nLPPvssLrnkEqxYsQIJCQl4/fXXsW7dOvh8Ptxzzz0AgF/96le49957kZycjP/4j/8AAOn4ZLPZsHTpUrS0tODuu+/GjBkz8NFHH+Hhhx9GW1sbfvWrX435NYx03Hzvvfdw3XXXYdGiRXj00UehVCqliwh79+7FZZddBgD4t3/7N2zbtg3f//73cfHFF6O7uxtVVVU4evQoFi5cOOb9I5oSBCKiGOvt7RUACCtXrhxx2dOnTwsqlUr42c9+FnT/Z599JiQkJATdv3TpUgGA8Mc//lG6z+l0CtOnTxdWrVol3ffJJ58IAIQ//OEPQ7YnruO3v/3tkMdsNtuQ++6++25Br9cLDodDum/NmjXCzJkzpX/fd999QkpKiuDxeEZ8vURENHVUVFQIKpVKUKlUwuLFi4Uf/ehHwttvvy24XK6g5cIdX6655hph9uzZQfddcsklwtKlS4cs+1//9V9CUlKScPz48aD7f/zjHwsqlUpobGyU7gMgPProo9K///CHPwgAhIaGBum+pUuXBm0n0nHT5/MJRUVFwjXXXCP4fL6g11NQUCBcffXV0n2pqanCPffcM2TfiUgQmF5ORDHX19cHADAYDCMu++qrr8Ln8+Hmm29GV1eX9N/06dNRVFSE999/P2j55ORk3H777dK/NRoNLrvsMpw6dSrq/dNqtfjOd74z5P7ExETp/61WK7q6urBkyRLYbDYcO3Ys4vqMRiMGBgbwzjvvRL0PREQ0+V199dXYt28fVqxYgZqaGjz11FO45pprkJubi507d0rLyY8vvb296OrqwtKlS3Hq1Cn09vaOuJ1XXnkFS5YsgclkCjpWLlu2DF6vF5WVlXF5fYcOHUJ9fT1uvfVWdHd3S9sdGBjAVVddhcrKSim93Wg0Yv/+/WhtbY3LvhBNZkwvJ6KYS0lJAeAPXEdSX18PQRBQVFQU9nG1Wh3077y8PCgUiqD7TCYTamtro96/3NzcsEVkjhw5gp/85Cd47733pAsHouFOitatW4e//vWvuO6665Cbm4vly5fj5ptvxrXXXhv1PhER0eT0pS99Ca+++ipcLhdqamrw2muvYdOmTVi9ejUOHTqEiy++GB9++CEeffRR7Nu3DzabLej5vb29SE1NHXYb9fX1qK2tRUZGRtjH41W0rb6+HgCwZs2aiMv09vbCZDLhqaeewpo1a5Cfn49Fixbhn/7pn3DHHXdg9uzZcdk3osmEQTcRxVxKSgpycnJw+PDhEZf1+XxQKBR46623oFKphjyenJwc9O9wywCAIAhR7598xEFksViwdOlSpKSk4Kc//SnmzJkDnU6H6upqPPTQQ8MWqsnMzMShQ4fw9ttv46233sJbb72FP/zhD7jjjjvwwgsvRL1fREQ0eWk0GnzpS1/Cl770JVx00UX4zne+g1deeQW33347rrrqKsydOxe//OUvkZ+fD41Gg7///e/YtGlTVIXQfD4frr76avzoRz8K+/hFF10U65cjbRcAfvGLX0RsJSYep2+++WYsWbIEr732GioqKvCLX/wCTz75JF599VVcd911cdk/osmCQTcRxcUNN9yA3/3ud9i3bx8WL14ccbk5c+ZAEAQUFBTE7KQhdCQ8Gh988AG6u7vx6quvory8XLq/oaEhqudrNBrceOONuPHGG+Hz+bBu3To899xz+M///E8UFhaOen+IiGjy+uIXvwgAaGtrw+uvvw6n04mdO3dixowZ0jKh06eAyMevOXPmoL+/H8uWLYvL/g63XcB/MT2abWdnZ2PdunVYt24dOjo6sHDhQvzsZz9j0E0XPM7pJqK4+NGPfoSkpCTcddddaG9vH/L4yZMn8cwzz+Cmm26CSqXC448/PmS0WhAEdHd3j3rbSUlJAPyj19ESR9Dl++ByubB58+YRnxu6j0qlEiUlJQAAp9MZ9T4QEdHk8v7774fNtPr73/8OACguLg57fOnt7cUf/vCHIc9LSkoKe+y6+eabsW/fPrz99ttDHrNYLPB4PGN9CdJ2xXXJLVq0CHPmzMHTTz+N/v7+Ic/r7OwEAHi93iHTsDIzM5GTk8PjIBE40k1EcTJnzhy89NJL+Na3voV58+bhjjvuwPz58+FyufDRRx/hlVdewbe//W3cd999eOKJJ/Dwww/j9OnT+PrXvw6DwYCGhga89tpr+O53v4sNGzaMettGoxG//e1vYTAYkJSUhMsvvxwFBQURn/OVr3wFJpMJa9aswQ9+8AMoFAr8f//f/xdV2vpdd92Fnp4eXHnllcjLy8OZM2fw61//GgsWLMC8efNGte9ERDR53HvvvbDZbPjGN76BuXPnSse4v/zlL5g1axa+853voL29XcqGuvvuu9Hf34//9//
      "text/plain": [
       "<Figure size 1000x800 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = np.arange(np.log10(np.min(parent_mass)), np.log10(np.max(parent_mass)), 0.2)\n",
    "xbins = (bins[1:] + bins[:-1]) / 2\n",
    "bins = 10**bins\n",
    "xbins = 10**xbins\n",
    "\n",
    "\n",
    "fig, axs = plt.subplots(\n",
    "    nrows=2, ncols=2, figsize=(10, 8), sharey=\"row\", sharex=\"col\", gridspec_kw={\"height_ratios\": [3, 1]})\n",
    "cols = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
    "fig.subplots_adjust(wspace=0.0)\n",
    "\n",
    "# Plot the projected velocity difference\n",
    "axs[0, 0].scatter(parent_mass[is_central], dvx_projected[is_central], s=0.5)\n",
    "axs[0, 1].scatter(parent_mass[~is_central], dvx_projected[~ is_central], s=0.5)\n",
    "plt.plot(xrange, sigma_vir(xrange), color=\"black\")\n",
    "axs[1, 0].plot(xrange, sigma_vir(xrange), color=\"black\")\n",
    "\n",
    "# Plot the standard deviation in a bin\n",
    "xrange = np.logspace(np.log10(np.min(parent_mass)), np.log10(np.max(parent_mass)), 100)\n",
    "stat, __, __ = binned_statistic(parent_mass[is_central], dvx_projected[is_central], statistic=\"std\", bins=bins)\n",
    "axs[1, 0].plot(xbins, stat, color=cols[0], label=r\"$\\sigma_{(\\mathbf{V} - \\mathbf{V}_{\\rm SPH}) \\cdot \\widehat{\\mathbf{r}}}$\")\n",
    "axs[1, 0].plot(xrange, sigma_vir(xrange), color=\"black\", label=r\"$\\sigma_{\\rm vir}$\")\n",
    "\n",
    "# Plot the standard deviation in a bin\n",
    "stat, __, __ = binned_statistic(parent_mass[~is_central], dvx_projected[~is_central], statistic=\"std\", bins=bins)\n",
    "axs[1, 1].plot(xbins, stat, color=cols[0])\n",
    "axs[1, 1].plot(xrange, sigma_vir(xrange), color=\"black\")\n",
    "m = np.isfinite(stat)\n",
    "beta, alpha = np.polyfit(np.log10(xbins[m] / 1e15), np.log10(stat[m]), 1)\n",
    "print(alpha, beta)\n",
    "axs[1, 1].plot(xbins, 10**alpha * (xbins / 1e15)**beta, color=\"red\", ls=\"--\", label=r\"$\\beta = {:.2f}$\".format(beta))\n",
    "\n",
    "\n",
    "axs[0, 0].set_ylabel(r\"$\\mathbf{V} \\cdot \\widehat{\\mathbf{r}}  ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[0, 0].set_ylabel(r\"$(\\mathbf{V} - \\mathbf{V}_{\\rm SPH}) \\cdot \\widehat{\\mathbf{r}}  ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[1, 0].set_ylabel(r\"$\\sigma  ~ [\\mathrm{km} / \\mathrm{s}]$\")\n",
    "axs[0, 0].set_title(\"Centrals\")\n",
    "axs[0, 1].set_title(\"Satellites\")\n",
    "for i in range(2):\n",
    "    axs[0, i].axhline(0, color=\"red\", ls=\"--\", alpha=0.5)\n",
    "    axs[1, i].set_xscale(\"log\")\n",
    "    # axs[1, i].set_yscale(\"log\")\n",
    "    axs[1, i].set_xlabel(r\"Parent mass $[M_\\odot / h]$\")\n",
    "\n",
    "axs[1, 0].legend()\n",
    "axs[1, 1].legend()\n",
    "\n",
    "\n",
    "fig.tight_layout(w_pad=0)\n",
    "fig.savefig(\"../plots/projected_residual_velocity_to_parent_mass.png\", dpi=300, bbox_inches=\"tight\")\n",
    "# fig.savefig(\"../plots/projected_velocity_to_parent_mass.png\", dpi=300, bbox_inches=\"tight\")\n",
    "\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
}