{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Examples of using velociraptor data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#load matplotlib as inline so figures display within the notebook\n", "%matplotlib inline\n", "from matplotlib.pylab import *" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#load other useful packages\n", "import numpy as np\n", "import astropy as ap\n", "import pynbody as pyn\n", "#for general python stuff\n", "import sys,os,string,time,re,struct\n", "import math,operator\n", "#for useful scipy stuff\n", "from scipy.stats.mstats import mquantiles\n", "from scipy.misc import comb\n", "import scipy.interpolate as scipyinterp\n", "import scipy.integrate as scipyint\n", "import scipy.optimize as scipyopt\n", "import scipy.special as scipysp\n", "import itertools\n", "#for useful mathematical tools\n", "from sklearn.neighbors import NearestNeighbors\n", "import scipy.spatial as spatial\n", "\n", "#to load specific functions defined in another python file\n", "basecodedir='/path/to/dir/'\n", "sys.path.append(basecodedir)\n", "import velociraptor_python_tools as vpt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#set some plotting parameters\n", "rcParams['figure.figsize'] = (10,6)\n", "rcParams['font.size'] = 18\n", "rc('text', usetex=True)\n", "rcParams['mathtext.fontset'] = 'custom'\n", "rcParams['mathtext.rm'] = 'Bitstream Vera Sans'\n", "rcParams['mathtext.it'] = 'Bitstream Vera Sans:italic'\n", "rcParams['mathtext.bf'] = 'Bitstream Vera Sans:bold'\n", "rcParams['xtick.direction']='in'\n", "rcParams['ytick.direction']='in'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading Halos" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#define formats\n", "ASCIIFORMAT=0\n", "HDFFORMAT=2\n", "\n", "#base filename \n", "inputfname='basevelociraptorfilename'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#now we try loading halo property data\n", "halopropdata,numhalos=vpt.ReadPropertyFile(inputfname,HDFFORMAT)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#then load the particles associated with halos\n", "haloparticledata=vpt.ReadParticleDataFile(inputfname,HDFFORMAT)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#list the keys available\n", "print(halopropdata.keys())\n", "print(halopropdata['SimulationInfo'].keys())\n", "print(halopropdata['UnitInfo'].keys())\n", "print(haloparticledata.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting halo mass function" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Set bin edges (log M)\n", "NBINS=50\n", "xlim=[7,15]\n", "#xlim=[0,5]\n", "deltax=(xlim[1]-xlim[0])/float(NBINS)\n", "xedges=np.arange(xlim[0],xlim[1]+0.5*deltax,deltax)\n", "xbins=np.arange(xlim[0]+0.5*deltax,xlim[1],deltax)\n", "\n", "#desired mass field\n", "massfield=['Mass_200crit','Mass_tot', 'Mass_200mean']\n", "massbindata=dict()\n", "#make histograms while converting units to Msun and kpc\n", "for key in massfield:\n", " massfac=np.log10(halopropdata['UnitInfo']['Mass_unit_to_solarmass'])\n", " massbindata[key],blah=np.histogram(np.log10(halopropdata[key])+massfac,xedges)\n", " massbindata[key]=np.float64(massbindata[key])/np.power(halopropdata['SimulationInfo']['Period']/halopropdata['SimulationInfo']['h_val'],3.0)\n", "#store mass value at which a halo is composed of 100 particles\n", "vlinedata=halopropdata[\"Mass_tot\"][np.where(halopropdata[\"npart\"]==100)][0]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "colors=['red','blue','green']\n", "vlabels=[r'$M_{200\\rho_c}$',r'$M_{\\rm{tot}}$',r'$M_{200\\rho_m}$']\n", "for i in range(len(massfield)):\n", " plot(np.power(10.0,xbins),massbindata[massfield[i]],color=colors[i],label=vlabels[i],linewidth=2)\n", "axvline(x=vlinedata,color='black',linewidth=2,linestyle='dashed')\n", "ylabel(r'$dn/d\\log M$ (kpc$^{-3}$)')\n", "xlabel(r'$M$ [M$_\\odot$]')\n", "yscale(\"log\")\n", "xscale(\"log\")\n", "legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting halo versus concentration" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#desired mass field\n", "massfield=['Mass_200crit','Mass_tot', 'Mass_200mean']\n", "massdata=dict()\n", "#get data\n", "sel=np.where(halopropdata['npart']>=1000)\n", "for key in massfield:\n", " massfac=np.log10(halopropdata['UnitInfo']['Mass_unit_to_solarmass'])\n", " massdata[key]=np.power(10.0,np.log10(halopropdata[key][sel])+massfac)\n", "cnfw=halopropdata['cNFW'][sel]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "markerstyle=['o','s','^']\n", "\n", "for i in range(len(massfield)):\n", " scatter(massdata[massfield[i]],cnfw,color=colors[i],label=vlabels[i],marker=markerstyle[i],alpha=0.5,s=50)\n", "axvline(x=vlinedata,color='black',linewidth=2,linestyle='dashed')\n", "ylabel(r'$c_{\\rm NFW}$')\n", "xlabel(r'$M$ [M$_\\odot$]')\n", "yscale(\"log\")\n", "xscale(\"log\")\n", "legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting Vmax distribution" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#Set bin edges (log M)\n", "NBINS=50\n", "xlim=[0.1,3.5]\n", "deltax=(xlim[1]-xlim[0])/float(NBINS)\n", "xedges=np.arange(xlim[0],xlim[1]+0.5*deltax,deltax)\n", "xbins=np.arange(xlim[0]+0.5*deltax,xlim[1],deltax)\n", "\n", "#desired mass field\n", "fields=['Vmax','sigV']\n", "bindata=dict()\n", "for key in fields:\n", " fac=np.log10(halopropdata['UnitInfo']['Velocity_unit_to_kms'])\n", " bindata[key],blah=np.histogram(np.log10(halopropdata[key])+fac,xedges)\n", " bindata[key]=np.float64(bindata[key])/np.power(halopropdata['SimulationInfo']['Period']/halopropdata['SimulationInfo']['h_val'],3.0)\n", "#store velocity value at which a halo is composed of 100 particles\n", "vlinedata=np.average(halopropdata[\"Vmax\"][np.where(halopropdata[\"npart\"]==100)])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "colors=['red','blue']\n", "vlabels=[r'$V_{\\rm max}$',r'$\\sigma_{v}$']\n", "fig1, ax = plt.subplots(figsize=(10,6))\n", "for i in range(len(fields)):\n", " ax.plot(np.power(10.0,xbins),bindata[fields[i]],color=colors[i],label=vlabels[i],linewidth=2)\n", "ax.axvline(x=vlinedata,color='black',linewidth=2,linestyle='dashed')\n", "ylabel(r'$dn/d\\log V$ (kpc$^{-3}$)')\n", "xlabel(r'$V$ [km/s]')\n", "yscale(\"log\")\n", "xscale(\"log\")\n", "ax.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting Vmax mass distribution" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#desired mass field\n", "massfield=['Mass_200crit','Mass_tot', 'Mass_200mean']\n", "massdata=dict()\n", "sel=np.where(halopropdata['npart']>=20)\n", "massfac=np.log10(halopropdata['UnitInfo']['Mass_unit_to_solarmass'])\n", "for key in massfield:\n", " massdata[key]=np.power(10.0,np.log10(halopropdata[key][sel])+massfac)\n", "vmax=halopropdata['Vmax'][sel]*halopropdata['UnitInfo']['Velocity_unit_to_kms']\n", "#store mass value at which a halo is composed of 100 particles\n", "massfac=halopropdata['UnitInfo']['Mass_unit_to_solarmass']\n", "vlinedata=halopropdata[\"Mass_tot\"][np.where(halopropdata[\"npart\"]==100)][0]*massfac\n", "#now get median and scatter in relation\n", "NBINS=50\n", "xlim=[7,15]\n", "deltax=(xlim[1]-xlim[0])/float(NBINS)\n", "xedges=np.arange(xlim[0],xlim[1]+0.5*deltax,deltax)\n", "xbins=np.power(10.0,np.arange(xlim[0]+0.5*deltax,xlim[1],deltax))\n", "xedges=np.power(10.0,xedges)\n", "nbindata=dict()\n", "bindata=dict()\n", "for key in massfield:\n", " nbindata[key]=np.zeros(NBINS)\n", " bindata[key]=np.zeros([NBINS,6])\n", " for ibin in range(NBINS):\n", " wdata=np.where((massdata[key]>=xedges[ibin])*(massdata[key]0)[0]\n", " bindata[key]=bindata[key][wdata].transpose()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "colors=['red','blue','green']\n", "vlabels=[r'$M_{200\\rho_c}$',r'$M_{\\rm{tot}}$',r'$M_{200\\rho_m}$']\n", "fig1, ax = plt.subplots(figsize=(10,6))\n", "for i in range(len(massfield)):\n", " print(massfield[i])\n", " ax.plot(bindata[massfield[i]][0],bindata[massfield[i]][3],color=colors[i],label=vlabels[i],linewidth=2,zorder=2)\n", " ax.fill_between(bindata[massfield[i]][0],bindata[massfield[i]][2],bindata[massfield[i]][4],\n", " facecolor=colors[i],alpha=0.35,interpolate=True,zorder=1)\n", " ax.fill_between(bindata[massfield[i]][0],bindata[massfield[i]][1],bindata[massfield[i]][5],\n", " facecolor=colors[i],alpha=0.15,interpolate=True,zorder=1)\n", "ax.axvline(x=vlinedata,color='black',linewidth=2,linestyle='dashed')\n", "xlabel(r'$M$ [M$_\\odot$]')\n", "ylabel(r'$V$ [km/s]')\n", "yscale(\"log\")\n", "xscale(\"log\")\n", "ax.legend()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.6.5" }, "widgets": { "state": {}, "version": "1.1.2" } }, "nbformat": 4, "nbformat_minor": 1 }