halo-units.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "71600f06-3607-415a-83dc-3c1644008dc7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "yt : [INFO     ] 2023-11-15 16:51:54,102 Loading plugins from /home/jwise/.config/yt/my_plugins.py\n"
     ]
    }
   ],
   "source": [
    "import yt\n",
    "import ytree\n",
    "yt.enable_plugins()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "0013e36b-474e-4959-9fcd-aea712a5ad50",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_76287/2857990963.py:2: UserWarning: This dataset appears to be of type EnzoDataset, but the following requirements are currently missing: libconf\n",
      "Please verify your installation.\n",
      "  ds = yt.load('DD0125/output_0125')\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,934 Parameters: current_time              = 75.89987395114\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,935 Parameters: domain_dimensions         = [64 64 64]\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,936 Parameters: domain_left_edge          = [0. 0. 0.]\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,936 Parameters: domain_right_edge         = [1. 1. 1.]\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,936 Parameters: cosmological_simulation   = 1\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,937 Parameters: current_redshift          = 11.181356517874\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,937 Parameters: omega_lambda              = 0.6889\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,938 Parameters: omega_matter              = 0.3111\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,938 Parameters: omega_radiation           = 0\n",
      "yt : [INFO     ] 2023-11-15 17:02:10,939 Parameters: hubble_constant           = 0.6766\n",
      "Parsing Hierarchy : 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 777/777 [00:00<00:00, 4655.97it/s]\n",
      "yt : [INFO     ] 2023-11-15 17:02:11,128 Gathering a field list (this may take a moment.)\n"
     ]
    }
   ],
   "source": [
    "arbor = ytree.load('rockstar_halos/trees/tree_0_0_0.dat')\n",
    "ds = yt.load('DD0125/output_0125')\n",
    "pfilters = ['dm', 'p2', 'p3']\n",
    "for pf in pfilters:\n",
    "    ds.add_particle_filter(pf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b5a2bd4e-d0e7-462f-8660-b83c1791fe48",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading tree roots: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 7934216/7934216 [00:00<00:00, 269986320.83it/s]\n"
     ]
    }
   ],
   "source": [
    "most_massive = arbor[0]  # most massive halo at last output\n",
    "\n",
    "# all length units are in comoving in ytree\n",
    "pos = most_massive['position']   # comoving\n",
    "rad = most_massive['virial_radius']  # comoving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "d75bc8ec-b779-4083-b490-978721507a43",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pos = [0.49441326 0.51073015 0.49368903] Mpc, radius = 13.925509452819824 kpc <== both comoving length\n"
     ]
    }
   ],
   "source": [
    "print(f\"pos = {pos.to('Mpc')}, radius = {rad.to('kpc')} <== both comoving length\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d2e65a1-53b4-47b5-bb6f-3da66faf2bdb",
   "metadata": {},
   "source": [
    "## Now let's try to inspect the dark matter, gas, stellar, and total mass within this halo (sphere)\n",
    "First we have to convert the comoving lengths to proper (without the universe's expansion factor).  The easiest way to do this is to create arrays / quantities associated with the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "760afc97-a80c-416a-84d9-6da514e40d4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_pos = ds.arr(pos.to('Mpc'), 'Mpccm')  # pos is comoving, and we have to specific 'Mpccm' for comoving length when creating this array\n",
    "ds_rad = ds.quan(rad.to('kpc'), 'kpccm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "b3b927fa-689c-4aac-bc0d-d9861fe54575",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ds values :: \tpos = [0.0405877  0.0419272  0.04052825] Mpc or [0.49441326 0.51073015 0.49368903] unitary\n",
      "\t\tradius = 1.1431821584701538 kpc <^ both proper length\n"
     ]
    }
   ],
   "source": [
    "print(f\"ds values :: \\tpos = {ds_pos.to('Mpc')} or {ds_pos.to('unitary')}\\n\\t\\tradius = {ds_rad.to('kpc')} <^ both proper length\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c289faa-043c-4c29-9b97-54fef7d99235",
   "metadata": {},
   "source": [
    "Now calculate the sphere masses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "661cc00e-01d7-4597-a92b-c8567e6c4def",
   "metadata": {},
   "outputs": [],
   "source": [
    "sp = ds.sphere(ds_pos, ds_rad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "37784460-6ed3-4b0d-b679-68136fc205d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "particle_mass = sp.quantities.total_quantity(('nbody', 'particle_mass'))\n",
    "star2_mass = sp.quantities.total_quantity(('p2', 'particle_mass'))\n",
    "star3_mass = sp.quantities.total_quantity(('p3', 'particle_mass'))\n",
    "gas_mass = sp.quantities.total_quantity(('gas', 'mass'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "fb4b3d04-6502-44d0-86db-be135f900fac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M_nbody = 80141058.29937911 Msun, M_p2 = 10498.532113547553 Msun, M_p3 = 403.959146873838 Msun, M_gas = 6685490.269615518 Msun\n"
     ]
    }
   ],
   "source": [
    "print(f\"M_nbody = {particle_mass.to('Msun')}, M_p2 = {star2_mass.to('Msun')}, M_p3 = {star3_mass.to('Msun')}, M_gas = {gas_mass.to('Msun')}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "51930339-d9cb-433d-85bf-5d6369f2b3e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_mass = particle_mass + gas_mass\n",
    "dm_mass = particle_mass - star2_mass - star3_mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "b9a8170b-dd05-4ec6-b5d3-3a387985e004",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M_DM = 80130155.80811867 Msun, M_total = 86826548.56899461 Msun\n"
     ]
    }
   ],
   "source": [
    "print(f\"M_DM = {dm_mass.to('Msun')}, M_total = {total_mass.to('Msun')}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ff9a8fe-3b57-4726-b7f5-c040d3e9cf2f",
   "metadata": {},
   "source": [
    "Check mass from ytree (rockstar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "a32562b1-8c30-46fa-9c4d-480bc6838353",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M_ytree = 79470880.0 Msun, difference = 0.830%\n"
     ]
    }
   ],
   "source": [
    "ytree_mass = most_massive['mass'].to('Msun')\n",
    "relative_diff = abs(ytree_mass - dm_mass) / ytree_mass\n",
    "print(f\"M_ytree = {ytree_mass.to('Msun')}, difference = {100*relative_diff.v:.3f}%\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b66da2f1-1ee9-4e8d-9457-8ac9ca696bee",
   "metadata": {},
   "source": [
    "Close!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "496d1c42-5369-4ed3-94bc-d59bc412fd17",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}