update

docs/source/Examples/Density_fields.ipynb

@@ -1 +1,267 @@
-{"cells":[{"metadata":{},"id":"b63aadc6","cell_type":"markdown","source":"# Creating density fields from snapshots\n\n[![Binder](https://mybinder.org/badge_logo.svg)](https://binder.flatironinstitute.org/v2/user/fvillaescusa/Quijote?filepath=/Tutorials/Density_fields.ipynb)"},{"metadata":{"trusted":true},"id":"bfa73080","cell_type":"code","source":"import numpy as np\nimport readgadget\nimport MAS_library as MASL","execution_count":1,"outputs":[]},{"metadata":{},"id":"4ed62239","cell_type":"markdown","source":"Define the value of the parameters"},{"metadata":{"trusted":true},"id":"d999bae9","cell_type":"code","source":"snapshot = '/home/jovyan/Data/Snapshots/fiducial/0/snapdir_004/snap_004' #location of the snapshot\ngrid     = 512    #the density field will have grid^3 voxels\nMAS      = 'CIC'  #Mass-assignment scheme:'NGP', 'CIC', 'TSC', 'PCS'\nverbose  = True   #whether to print information about the progress\nptype    = [1]    #[1](CDM), [2](neutrinos) or [1,2](CDM+neutrinos)","execution_count":2,"outputs":[]},{"metadata":{},"id":"9098a761","cell_type":"markdown","source":"Read the header and the particle positions"},{"metadata":{"trusted":true},"id":"8427b23a","cell_type":"code","source":"# read header\nheader   = readgadget.header(snapshot)\nBoxSize  = header.boxsize/1e3  #Mpc/h\nredshift = header.redshift     #redshift of the snapshot\nMasses   = header.massarr*1e10 #Masses of the particles in Msun/h\n\n# read positions, velocities and IDs of the particles\npos = readgadget.read_block(snapshot, \"POS \", ptype)/1e3 #positions in Mpc/h","execution_count":3,"outputs":[]},{"metadata":{},"id":"62506aef","cell_type":"markdown","source":"Print some information about the data"},{"metadata":{"trusted":true},"id":"d4808044","cell_type":"code","source":"print('BoxSize: %.3f Mpc/h'%BoxSize)\nprint('Redshift: %.3f'%redshift)\nprint('%.3f < X < %.3f'%(np.min(pos[:,0]), np.max(pos[:,0])))\nprint('%.3f < Y < %.3f'%(np.min(pos[:,1]), np.max(pos[:,1])))\nprint('%.3f < Z < %.3f'%(np.min(pos[:,2]), np.max(pos[:,2])))","execution_count":4,"outputs":[{"output_type":"stream","text":"BoxSize: 1000.000 Mpc/h\nRedshift: 0.000\n0.000 < X < 999.992\n0.000 < Y < 999.992\n0.000 < Z < 999.992\n","name":"stdout"}]},{"metadata":{},"id":"3d75aa92","cell_type":"markdown","source":"Define the matrix that will contain the value of the density / overdensity field"},{"metadata":{"trusted":true},"id":"951f674a","cell_type":"code","source":"delta = np.zeros((grid,grid,grid), dtype=np.float32)","execution_count":5,"outputs":[]},{"metadata":{},"id":"92a7fc19","cell_type":"markdown","source":"Now construct the 3D density field"},{"metadata":{"trusted":true},"id":"4c3b5a85","cell_type":"code","source":"# construct 3D density field\nMASL.MA(pos, delta, BoxSize, MAS, verbose=verbose)","execution_count":6,"outputs":[{"output_type":"stream","text":"\nUsing CIC mass assignment scheme\nTime taken = 5.884 seconds\n\n","name":"stdout"}]},{"metadata":{},"id":"05947788","cell_type":"markdown","source":"We can make some tests to make sure the density field has been computed properly"},{"metadata":{"trusted":true},"id":"ad8baf57","cell_type":"code","source":"# the sum of the values in all voxels should be equal to the number of particles\nprint('%.3f should be equal to\\n%.3f'%(np.sum(delta, dtype=np.float64), pos.shape[0]))","execution_count":7,"outputs":[{"output_type":"stream","text":"134217728.019 should be equal to\n134217728.000\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"As this point, delta contains the effective number of particles in each voxel. \nIf you want instead the effective mass in each voxel you can just do"},{"metadata":{"trusted":true},"cell_type":"code","source":"delta *= Masses[1]\n\n# now check that the mass in the density field is equal to the total mass in the simulation\nprint('%.3e should be equal to\\n%.3e'%(np.sum(delta, dtype=np.float64), pos.shape[0]*Masses[1]))","execution_count":8,"outputs":[{"output_type":"stream","text":"8.812e+19 should be equal to\n8.812e+19\n","name":"stdout"}]},{"metadata":{},"id":"77d940b3","cell_type":"markdown","source":"If needed, the overdensity is easy to calculate"},{"metadata":{"trusted":true},"id":"4593380c","cell_type":"code","source":"# at this point, delta contains the effective number of particles in each voxel\n# now compute overdensity and density constrast\ndelta /= np.mean(delta, dtype=np.float64);  delta -= 1.0\n\nprint('%.3f < delta < %.3f'%(np.min(delta), np.max(delta)))\nprint('<delta> = %.3f'%np.mean(delta))","execution_count":9,"outputs":[{"output_type":"stream","text":"-1.000 < delta < 1195.511\n<delta> = -0.000\n","name":"stdout"}]}],"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3 (ipykernel)","language":"python"},"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":5}
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "b63aadc6",
+   "metadata": {},
+   "source": [
+    "# Creating density fields from snapshots\n",
+    "\n",
+    "[![Binder](https://mybinder.org/badge_logo.svg)](https://binder.flatironinstitute.org/v2/user/fvillaescusa/Quijote?filepath=/Tutorials/Density_fields.ipynb)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "bfa73080",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import readgadget\n",
+    "import MAS_library as MASL"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4ed62239",
+   "metadata": {},
+   "source": [
+    "Define the value of the parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "d999bae9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "snapshot = '/home/jovyan/Data/Snapshots/fiducial/0/snapdir_004/snap_004' #location of the snapshot\n",
+    "grid     = 512    #the density field will have grid^3 voxels\n",
+    "MAS      = 'CIC'  #Mass-assignment scheme:'NGP', 'CIC', 'TSC', 'PCS'\n",
+    "verbose  = True   #whether to print information about the progress\n",
+    "ptype    = [1]    #[1](CDM), [2](neutrinos) or [1,2](CDM+neutrinos)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9098a761",
+   "metadata": {},
+   "source": [
+    "Read the header and the particle positions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "8427b23a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# read header\n",
+    "header   = readgadget.header(snapshot)\n",
+    "BoxSize  = header.boxsize/1e3  #Mpc/h\n",
+    "redshift = header.redshift     #redshift of the snapshot\n",
+    "Masses   = header.massarr*1e10 #Masses of the particles in Msun/h\n",
+    "\n",
+    "# read positions, velocities and IDs of the particles\n",
+    "pos = readgadget.read_block(snapshot, \"POS \", ptype)/1e3 #positions in Mpc/h"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "62506aef",
+   "metadata": {},
+   "source": [
+    "Print some information about the data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "d4808044",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "BoxSize: 1000.000 Mpc/h\n",
+      "Redshift: 0.000\n",
+      "0.000 < X < 999.992\n",
+      "0.000 < Y < 999.992\n",
+      "0.000 < Z < 999.992\n"
+     ]
+    }
+   ],
+   "source": [
+    "print('BoxSize: %.3f Mpc/h'%BoxSize)\n",
+    "print('Redshift: %.3f'%redshift)\n",
+    "print('%.3f < X < %.3f'%(np.min(pos[:,0]), np.max(pos[:,0])))\n",
+    "print('%.3f < Y < %.3f'%(np.min(pos[:,1]), np.max(pos[:,1])))\n",
+    "print('%.3f < Z < %.3f'%(np.min(pos[:,2]), np.max(pos[:,2])))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3d75aa92",
+   "metadata": {},
+   "source": [
+    "Define the matrix that will contain the value of the density / overdensity field"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "951f674a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "delta = np.zeros((grid,grid,grid), dtype=np.float32)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "92a7fc19",
+   "metadata": {},
+   "source": [
+    "Now construct the 3D density field"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "id": "4c3b5a85",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "Using CIC mass assignment scheme\n",
+      "Time taken = 5.884 seconds\n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "# construct 3D density field\n",
+    "MASL.MA(pos, delta, BoxSize, MAS, verbose=verbose)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "05947788",
+   "metadata": {},
+   "source": [
+    "We can make some tests to make sure the density field has been computed properly"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "id": "ad8baf57",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "134217728.019 should be equal to\n",
+      "134217728.000\n"
+     ]
+    }
+   ],
+   "source": [
+    "# the sum of the values in all voxels should be equal to the number of particles\n",
+    "print('%.3f should be equal to\\n%.3f'%(np.sum(delta, dtype=np.float64), pos.shape[0]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a6deccd3",
+   "metadata": {},
+   "source": [
+    "As this point, delta contains the effective number of particles in each voxel. \n",
+    "If you want instead the effective mass in each voxel you can just do"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "id": "1531aaa3",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "8.812e+19 should be equal to\n",
+      "8.812e+19\n"
+     ]
+    }
+   ],
+   "source": [
+    "delta *= Masses[1]\n",
+    "\n",
+    "# now check that the mass in the density field is equal to the total mass in the simulation\n",
+    "print('%.3e should be equal to\\n%.3e'%(np.sum(delta, dtype=np.float64), pos.shape[0]*Masses[1]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "77d940b3",
+   "metadata": {},
+   "source": [
+    "If needed, the overdensity is easy to calculate"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "id": "4593380c",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "-1.000 < delta < 1195.511\n",
+      "<delta> = -0.000\n"
+     ]
+    }
+   ],
+   "source": [
+    "# at this point, delta contains the effective number of particles in each voxel\n",
+    "# now compute overdensity and density constrast\n",
+    "delta /= np.mean(delta, dtype=np.float64);  delta -= 1.0\n",
+    "\n",
+    "print('%.3f < delta < %.3f'%(np.min(delta), np.max(delta)))\n",
+    "print('<delta> = %.3f'%np.mean(delta))"
+   ]
+  }
+ ],
+ "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.9.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/source/tutorials.rst

@@ -8,7 +8,8 @@ We provide multiple tutorials showing how to read and manipulate Quijote data. I
 
 
 .. nbgallery::
-
+
+   Examples/Reading_snapshots
    Examples/Density_fields
    Examples/Density_fields2