Metallicity and SSFR.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pynbody\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import numpy as np\n",
    "import pynbody.plot as pp\n",
    "import pickle\n",
    "import pandas as pd\n",
    "\n",
    "# set the config to prioritize the AHF catalog\n",
    "pynbody.config['halo-class-priority'] =  [pynbody.halo.ahf.AHFCatalogue,\n",
    "                                          pynbody.halo.GrpCatalogue,\n",
    "                                          pynbody.halo.AmigaGrpCatalogue,\n",
    "                                          pynbody.halo.legacy.RockstarIntermediateCatalogue,\n",
    "                                          pynbody.halo.rockstar.RockstarCatalogue,\n",
    "                                          pynbody.halo.subfind.SubfindCatalogue,\n",
    "                                          pynbody.halo.hop.HOPCatalogue]\n",
    "\n",
    "# we define the operation as a function, which takes in filepath and halo nums as input\n",
    "# you should rename 'operation' to correspond to whatever it is doing\n",
    "# i.e. if it computes gas fraction you should change the line below to `def gasFrac(filepath,halo_nums):`\n",
    "\n",
    "# prop is the property that you are computing\n",
    "prop = 'GasMetallicity'\n",
    "\n",
    "def operation(filepath,halo_nums):\n",
    "    # first, load in the simulation and halos\n",
    "    s = pynbody.load(filepath)\n",
    "    s.physical_units()\n",
    "    h = s.halos()\n",
    "    # then we open the .data file\n",
    "    with open(filepath+prop+'.data','wb') as f:\n",
    "        # we loop through all the halos, compute a value for that halo, and add it to the .data file\n",
    "        for halo_num in halo_nums:\n",
    "            # we load the copy of the halo to minimize computational stress\n",
    "            halo = h.load_copy(halo_num)\n",
    "            halo.physical_units()\n",
    "            \n",
    "            ###################################\n",
    "            # here is where you compute your stuff\n",
    "            # i.e. calculate the color, gas fraction, mass of each halo\n",
    "            # so in this example we will just compute the stellar mass and the gas mass\n",
    "            \n",
    "            gem = np.sum(h[halo_num].gas['HI']*halo.gas['metals']*halo.gas['mass'])\n",
    "            mzgas = gem/np.sum(halo.gas['mass']*h[halo_num].gas['HI'])\n",
    "            \n",
    "            \n",
    "            sfr, bins = pynbody.plot.stars.sfh(halo, filename=None, massform=False, clear=False, legend=False, subplot=False, trange=False, bins=69)\n",
    "            SFR = sfr[-1]\n",
    "            sSFR = SFR/np.sum(halo.star['mass'])\n",
    "            \n",
    "            #print(halo_num,nstar,ngas)\n",
    "            ###################################\n",
    "            \n",
    "            # now we add these to the .data file\n",
    "            pickle.dump({\n",
    "                'haloid': halo_num,  # always put the haloid here\n",
    "                'gas HI': h[halo_num].gas['HI'],\n",
    "                'gas Metals': halo.gas['metals'],\n",
    "                'gmas': halo.gas['mass'],\n",
    "                'MZgas': mzgas,\n",
    "                'SFR': SFR,\n",
    "                'sSFR': sSFR\n",
    "            \n",
    "            },f,pickle.HIGHEST_PROTOCOL)\n",
    "            \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/engelhar/.local/lib/python3.5/site-packages/pynbody/array.py:346: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return np.ndarray.__truediv__(self, rhs)\n",
      "/home/engelhar/.local/lib/python3.5/site-packages/pynbody/plot/stars.py:632: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  binnorm = 1e-9 * bins / (trange[1] - trange[0])\n"
     ]
    }
   ],
   "source": [
    "# for ease of use I have the filepaths for all of the simulations here\n",
    "#sim1 = '/home/akinshol/Data/Sims/h148.cosmo50PLK.3072g3HbwK1BH/h148.cosmo50PLK.3072g3HbwK1BH.004096/h148.cosmo50PLK.3072g3HbwK1BH.004096'\n",
    "#sim2 = '/home/akinshol/Data/Sims/h229.cosmo50PLK.3072gst5HbwK1BH/h229.cosmo50PLK.3072gst5HbwK1BH.004096/h229.cosmo50PLK.3072gst5HbwK1BH.004096'\n",
    "sim3 = '/home/akinshol/Data/Sims/h242.cosmo50PLK.3072gst5HbwK1BH/h242.cosmo50PLK.3072gst5HbwK1BH.004096/h242.cosmo50PLK.3072gst5HbwK1BH.004096'\n",
    "#sim4 = '/home/akinshol/Data/Sims/h329.cosmo50PLK.3072gst5HbwK1BH/h329.cosmo50PLK.3072gst5HbwK1BH.004096/h329.cosmo50PLK.3072gst5HbwK1BH.004096'\n",
    "\n",
    "# here is where we will put whatever halo numbers we decide are interesting and worth computing\n",
    "# i.e. all the halos with stars in them126, 129, 135, 170, 195, 444, 686, 942, 1418]\n",
    "\n",
    "nums1 = [1, 2, 3, 5, 6, 7, 9, 10, 11, 12, 13, 14, 17, 18, 20, 21, 22, 23, 26, 27, 28, 30, 31, 32, 34, 36, 37, 38, 40, 41, 42, 45, 46, 48, 55, 57, 58, 60, 64, 71, 73, 77, 80, 91, 94, 95, 99, 106, 115, 121, 125, 126, 131, 140, 143, 160, 223, 252, 264, 271, 304, 353, 372, 373, 435, 465, 590, 647, 677, 682, 738, 869, 961, 980, 1146, 1155, 1381, 2792, 4897, 10814]\n",
    "nums2 = [1, 2, 3, 4, 7, 14, 16, 17, 18, 19, 21, 22, 25, 27, 28, 29, 30, 33, 36, 41, 50, 51, 52, 56, 60, 62, 70, 73, 75, 95, 104, 108, 134, 203, 277, 553, 982, 1319, 1409, 1481, 4380, 5722]\n",
    "nums3 = [1, 4, 9, 10, 11, 12, 19, 24, 29, 30, 33, 36, 39, 40, 45, 46, 48, 53, 57, 62, 66, 69, 70, 72, 75, 76, 85, 89, 102, 133, 152, 185, 211, 302, 425, 457, 536, 1773, 2748, 2891, 9059, 9126, 10085, 11901]\n",
    "nums4 = [1, 8, 9, 13, 14, 19, 25, 31, 32, 40, 47, 63, 92, 99, 126, 129, 135, 170, 195, 444, 686, 942, 1418]\n",
    "#nums4 = [104, 108, 134, 203, 277, 553, 982, 1319, 1409, 1481, 4380, 5722]\n",
    "\n",
    "sim  = sim3\n",
    "nums = nums3\n",
    "\n",
    "# here we perform the operation, which in your code would be named by what it actually does\n",
    "operation(sim,nums)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ8AAAEQCAYAAABvBHmZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+8FPV97/HXOwgCakVJbMlBLxgM\noBJFj7+C9ZpEA4miaGyUNrdNJVJzY2PbW/rAxjyMbb2YktwmJOYHicSaeCXGGgpCY34YrmmjBogo\nGDSixnKIDUjEEIUE9HP/2DmwHHbP7p4zM7t75v18PPbBme/Oznx3GPic78xnPl9FBGZmZnl6XbM7\nYGZmxePgY2ZmuXPwMTOz3Dn4mJlZ7hx8zMwsdw4+ZmaWOwcfMzPLnYOPmZnl7qBmdyANkn4f+CNK\n3+f4iHhrk7tkZma9aPrIR9IiSVskre/RPk3Sk5I2Sprb2zYi4gcRcTVwL/DPWfbXzMz6T80uryPp\nHODXwO0RcWLSNgj4KXA+0AWsAmYCg4B5PTZxZURsST53FzArInbk1H0zM+uDpl92i4gHJI3p0Xw6\nsDEingGQtBi4OCLmARdW2o6kY4CXHHjMzFpf04NPFR3AprLlLuCMGp+ZBXyl2puSZgOzAQ455JBT\nJ0yY0N8+mpkVypo1a16IiDeksa1WDT4Ni4gbary/UNLzwPRRo0adunr16px6ZmY2MEh6Lq1tNT3h\noIrNwNFly6OTtn6JiGURMfvwww/v76bMzKwfWjX4rAKOkzRW0hDgCmBpfzcqabqkhS+99FK/O2hm\nZn3X9OAj6U7gQWC8pC5JsyJiD3ANcB+wAbgrIh7v77488jEzaw1Nv+cTETOrtK8AVqS5L0nTgenj\nxo1Lc7Nm1kZ2795NV1cXu3btanZXWtbQoUMZPXo0gwcPzmwfTX/Opxk6OzvDCQdmxfTss89y2GGH\nMXLkSCQ1uzstJyLYtm0bO3bsYOzYsfu9J2lNRHSmsZ+mj3wGgiWPbGb+fU/y8+07eeOIYcyZOp4Z\nkzua3S0zq2DXrl2MGTPGgacKSYwcOZKtW7dmup9CBZ8sLrsteWQz192zjp27XwVg8/adXHfPOgAH\nILMW5cDTuzyOT9MTDvKURcLB/Pue3Bt4uu3c/Srz73sytX2YWTGMGTOGF154AYBDDz205vof+9jH\n+MQnPtHrOkuWLOEnP/lJKv1LU6GCTxap1j/fvrOhdjOzPDn4tIAsRj5vHDGsoXYzay9LHtnMlJvv\nZ+zc5Uy5+X6WPNLv592ZMWMGp556KieccAILFy5s6LM33XQTb37zmzn77LN58sl9V1i+9KUvcdpp\np3HSSSfxnve8h1deeYUf/vCHLF26lDlz5nDyySfz9NNPV1yvGQoVfLIwZ+p4hg0etF/bsMGDmDN1\nfJN6ZGZp6b6nu3n7ToJ993T7G4AWLVrEmjVrWL16NQsWLGDbtm11fW7NmjUsXryYtWvXsmLFClat\nWrX3vUsvvZRVq1bx6KOPMnHiRG699Vbe+ta3ctFFFzF//nzWrl3Lm970porrNUOhEg6y0J1U4Gw3\ns4Gnt3u6/fk3vmDBAr75zW8CsGnTJp566ilGjhxZ83M/+MEPuOSSSxg+fDgAF1100d731q9fz/XX\nX8/27dv59a9/zdSpUytuo971slao4JPVQ6YzJnc42JgNQFnc0125ciXf/e53efDBBxk+fDjnnntu\nKg+8vv/972fJkiWcdNJJ3HbbbaxcubJf62WtUJfdXF7HzBqRxT3dl156iSOOOILhw4fzxBNP8NBD\nD9X92XPOOYclS5awc+dOduzYwbJly/a+t2PHDkaNGsXu3bu544479rYfdthh7Nixo+Z6eStU8DEz\na0QW93SnTZvGnj17mDhxInPnzuXMM8+s+7OnnHIKl19+OSeddBLvete7OO200/a+9/d///ecccYZ\nTJkyhfL5yq644grmz5/P5MmTefrpp6uulzeX1zGzQtmwYQMTJ06se/2iVjCpdJxcXqePXFjUzBrl\ne7rZKNRlN9/zMTNrDYUKPmZm1hoKddmtHRX1erNZliLCxUV7kUcugEc+LSyrp6vNimzo0KFs27Yt\nl/9g21H3fD5Dhw7NdD8e+bSwrJ6uNiuy0aNH09XVlfl8Ne2seybTLDn4tDBXzDZL3+DBgw+YodPy\nV6jLbllMqZAlV8w2s4GqUMGn3VKtXTHbzAYqX3ZrYa6YbWYDlYNPi/PT1WY2EBXqspuZmbUGBx8z\nM8udL7u1KVc+MLN21vbBR9IxwALgl8BPI+LmJncpc92VD7ofQO2ufAA4AJlZW2jqZTdJiyRtkbS+\nR/s0SU9K2ihpbo3NTALujogrgcmZdbaF9Fb5wMysHTR75HMb8Fng9u4GSYOAW4DzgS5glaSlwCBg\nXo/PXwk8BNwt6Urgqzn0uelc+cDM2l1TRz4R8QCly2XlTgc2RsQzEfFbYDFwcUSsi4gLe7y2AH8K\n3BARbwcuyPcbNIcrH5hZu2vFbLcOYFPZclfSVs23gA9L+gLws2orSZotabWk1e1eUNCVD8ys3TX7\nslu/RcR64LI61lso6Xlg+pAhQ07NvmfZceUDM2t3rRh8NgNHly2PTtr6LSKWAcs6OzuvSmN7zeTK\nB2bWzlrxstsq4DhJYyUNAa4Alqax4Xaram1mNlA1O9X6TuBBYLykLkmzImIPcA1wH7ABuCsiHk9j\nf+1W1drMbKBq6mW3iJhZpX0FsCLt/UmaDkwfN25c2ps2M7MGtOJlt8x45GNm1hoKFXzMzKw1FCr4\nOOHAzKw1FCr4+LKbmVlrKFTwMTOz1lCo4OPLbmZmraFQwceX3czMWkOhgo+ZmbWGVqztlpmiPmTq\nKbfNrNUUauRTxMtu3VNub96+k2DflNtLHkmlVquZWZ8UKvgUkafcNrNW5OAzwHnKbTNrRQ4+A5yn\n3DazVlSo4FPE53w85baZtaJCBZ8iJhzMmNzBvEsn0TFiGAI6Rgxj3qWTnO1mZk1VqFTr/mrXlGVP\nuW1mrcbBp07dKcvdmWPdKcuA/2O3vdr1FxSzvBXqslt/OGXZavEzVWb1qznykXRkHdt5LSK2p9Cf\nluWUZault19QPPox2189l91+nrzUyzqDgGNS6VGG+lNe540jhrG5QqBxyrJ18y8oZvWr57Lbhog4\nNiLGVnsB27LuaBr6k+3mlGWrxc9UmdWvnuBzVkrrtDWnLFst/gXFrH41L7tFxK401hkInLJsvek+\nN5ztZlZbn1KtJR0EnABsiohfptsls/blX1DM6tNwqrWkUcBtwLHA30q6JO1OmZnZwNaX53z+CvgI\n0BURfw28L90umZnZQNeX4DMceBk4OeW+9Imk4yXdJenzki5rdn/MzKy2vgSfTwP/CPxI0l8B3+nr\nziUtkrRF0voe7dMkPSlpo6S5NTbzLuAzEfFB4I/72hczM8tPwwkHEfFTSXOAc4BvR8T6Wp/pxW3A\nZ4HbuxskDQJuAc4HuoBVkpZSepB1Xo/PXwl8FbhB0kXAyH70xczMctKnbLeI2AZ8s787j4gHJI3p\n0Xw6sDEingGQtBi4OCLmARdW2dSHkqB1T3/7ZGZm2auntts7gMciYmsO/QHoADaVLXcBZ1RbOQle\nfwscAszvZb3ZwGyAY45p+UpAuXM1ZjPLUz0jn+8AWyS9BqwH1gGPJX8+HhG/ybB/NUXEz0iCSo31\nFgILATo7OyPjbrUVTxdhZnmrJ+HgzykVFl0A/APwBHAq8EnguQz6tBk4umx5dNLWb/VOo73kkc1M\nufl+xs5dzpSb7x/wJfE9XYSZ5a1m8ImIW4ApQACfAnYD10bE2yLi9zLo0yrgOEljJQ0BrgCWprHh\negqLFnFOFldjNrO81ZVqHRE7I+LjwNuAcZTSrKveh6mXpDuBB4HxkrokzYqIPcA1wH3ABuCuiHi8\nv/tK9ldz5FPEUYCrMZtZ3upJODgHmJC8JgJHATtIIa05ImZWaV8BrOjv9itsdxmwrLOz86pq6xRx\nFDBn6vj97vmAqzGbWbbqSThYCawFFgMLkhv8bameyeTSnjSuHbLIXI3ZzPKmiN4TvyT9GXAi8BZK\no59tlDLd1gHrI2JJ1p1MW2dnZ6xevbriez0zv6A0CujL3D1pbsvMrNkkrYmIzjS2Vc98Pl/ssfPR\nwCRKweg9QNsEn3pGPmmOAnq7f+TgY2ZFVnPkMxD1NvJJ09i5y6l0dAU8e/MFme/fzCxNaY58ama7\nSfpxGusUkbPIzMwqqyfhYKKkx3p5X0D1B2daSD2X3dLkLDIzs8rqCT4T6ljn1dqrNF89qdZpchaZ\nmVll9SQcZFFCpzBmTO5wsDEz66FPUyq0q6wvu7XDMz1mZq2gLzOZtq16arv1VRFrwpmZ9VXdwUfS\nH0g6LPn5ekn3SDolu661lyLWhDMz66tGRj4fjYgdks4GzgNuBT6fTbfaT7Xab5VK9ZiZFV0jwaf7\n1/oLgIURsRwYkn6X2lO1Z3cEvvTW5oo2v5NZHhoJPpslfRG4HFgh6eAGP9909U4m1xdzpo5HFdoD\nfOmtjflenlk2Ggke76U0x87UiNgOHAnMyaRXGcky4WDG5I6KpXRgYE/HMND5Xp5ZNupOtY6IV4B7\nypafB57PolPtqiPl6Ris+Yo4v5NZHpztlqI5U8czbPCg/dpcTqe9uT6fWTac7ZaiGZM7mHfpJDpG\nDEOURkKeu6e9+RcKs2w0UuHggGw3Sf+QQZ/amsvpDCyuz2eWjUaCT3e22/nAx9s1242y8jouh2P1\n8C8UZukrbLabU2jNzJqn7uCTZLt9HzhC0jnAccCurDqWNafQmpk1T92X3SR9ALgWGA2sBc4EHgTe\nnk3XsuUUWjOz5mnkstu1wGnAcxHxNmAysD2TXuXAKbRmZs3TSPDZFRG7ACQdHBFPAG2bb1ophVaU\n7v24fld+XDfNrJgayXbrkjQCWAJ8R9KLQNvOclqeQrt5+04Ee8vjdCcflK9n6etO+ui+9+bjblYc\njSQcXBIR2yPiY8BHgS8DF2fVsTzMmNzBf8x9Ox0jhh1Ql83JB9lz0odZcdUc+UhaWu0t4CrgolR7\nVIOkY4GPAIdHxGVJ2yHA54DfAisj4o5Gtunkg+bwcTcrrnouu50FbALuBB6GijMH1EXSIuBCYEtE\nnFjWPg34NDAI+HJE3FxtGxHxDDBL0t1lzZcCd0fEMklfBxoKPm/spSBotQdR/YBq//V23M1sYKsn\n+PwepaoGM4E/BJYDd0bE433Y323AZ4HbuxskDQJuSfbRBaxKRluDgHk9Pn9lRGypsN3RwLrk51cr\nvF9RdwDpec8HSvW73jbhDQfck/iLr6/lb+95jN2vBbtfjb3t192zjtXP/ZLvP7HVAalOc6aO3+/4\nguummRVFzeATEa8C3wK+lZTUmQmslHRjRHy2kZ1FxAOSxvRoPh3YmIxokLQYuDgi5lEaJdWji33P\nH1W8jyVpNjAb4JhjjjngZnfA3gDUkQSOSvckAF7Z/doBbTt3v8rXHvrPvcu+eV6b66aZFZciqk2B\nVrZSKehcQCnwjAGWAosiouG82CT43Nt92U3SZcC0iPhAsvw/gDMi4poqnx8J3ERppPTliJiX3PP5\nLKWKC/9e655PZ2dnHHzZP1a85NMxYhj/Mbf03OzYucurThBXr/LtmZm1M0lrIqIzjW3Vk3BwO3Ai\nsAK4MSLWp7HjvoqIbcDVPdpeBv601mfLC4vuqeNmd7V7Eo3wzXMzswPVk2r9Pkp13K4FfijpV8lr\nh6RfpdCHzcDRZcujk7bUlRcWrafCQaUHURvlm+dmZgeqGXwi4nURcVjy+p2y12ER8Tsp9GEVcJyk\nsZKGAFdQuqyXOknTJS186aWX6pokrHtyuBHDBh+wrcGvE6+rkffnm+dmZpXVdc8ntZ1JdwLnAq8H\nfgHcEBG3Sno38ClKGW6LIuKmLPvR2dkZq1evbihdutK6wAHZWj2TFiptz2naZtaO0rznk2vwabay\nez5XPfXUU6lss9FA0jPLDkojJE+3bWatzsGnn7pHPtVkOTKZcvP9NbPszMxaUa7ZbgNJz2m0K6lU\n7HLONx7lxmWPs/2V3f0ORi4pY2bWQGFRSX8g6bDk5+sl3SPplOy6lr7ybLdqKj1Yuvu14MVXdqcy\n3bbnETIza2w+n49GxA5JZwPnAbcCn8+mW81TzwikP5WX68myMzMb6BoJPt3DgQuAhRGxHBiSfpey\nU55qXU29I5C+XibrTt/uGDEMUbrX42QDMyuauhMOJN1L6eHPd1KaQnsn8KOIOCm77mWjs7Mzrv/S\nv1atVj3nG4+y+7Xej4sTBMysaNJMOGhk5PNe4D7gnRGxHTgSmJNGJ/K2/ZXdXHfPOjZv33nAfZwZ\nkzs4dGjveRi+TGZm1j+NZLv9dfLnidJ+j/Z/O73u5OO/frWL11eZQXPG5A62v7K76md7e3jUzMzq\n00jwebns56GUpjvYkG53stWdan3QiFEV3+++j1OtoKgvtZmZpaPu4BMRnyxflvQJSpfh2kZELAOW\nHTp6/FWV3u9ONkh7kjOX0zEz219/HjIdTqkCddvZ/eprFWcu7Q4uaU5yVumhVU8yZ2ZFV3fwkbSO\nff9fDwLeAPxdFp3KQ3ngOWL4YG6YfsJ+wWDG5I5UgkOlh1bL7y+ZmRVRIyOf8imt9wC/iIg9Kfcn\nU9Xu+eyqMC12WlxOx8zsQHWnWkfEc8AIYDpwCXB8Vp3KSnd5ndcNPWS/9v5ULKjF5XTMzA7USG23\na4E7gKOS1x2S/jyrjuUtq5GIy+mYmR2okctus4AzIuJlAEkfBx4EPpNFx/KW1UgkzeQFM7OBopHg\nI/bVdyP5ucZE0u0h65FIWskLZmYDRSPB5yvAw5K+mSzPoFTZuu0cfcRwfnfEMI9EzMyapKGZTJP5\ne85OFn8QEY9k0quM1ZrJtC/8IKnlzeec5a0phUVVKuh2PDAiIhYA2ySdnkYn8lLPlAp90f0gaaVC\npWZZ8Dln7a6RqtafA84CZibLO4BbUu9RhuqZybQvenuQ1CwLPues3TVyz+eMiDhF0iMAEfGipLaa\nTC4rfpDU8uZzztpdI8Fnt6RBJJVpJL0ByK40QA7SumZerQq2HyS1rPics3bXyGW3BcA3gaMk3QT8\nO/C/M+lVDtK8Zu4HSS1vPues3TUypcIdktYA76D0fM+MiGir+XzKpVnw0w+SWt58zlm7ayjVeqDo\n7OyMbefdSKVvLuDZmy/Iu0tmZi2vWanWfyDpsOTn6yXdkzz3kytJx0q6VdLdvbXV4oKfZmbN08g9\nn49GxA5JZwPnUapu8PlGdiZpkaQtktb3aJ8m6UlJGyXN7W0bEfFMRMyq1VZLf66ZL3lkM1Nuvp+x\nc5cz5eb7/WyFmVmDGgk+3TdILgAWRsRyoNFU69uAaeUNSQbdLcC7KD3EOlPS8ZImSbq3x+uoBvdX\n1YzJHcy7dBIdI4YhoGPEMOZdOqnmNXM/3Gdm1n+NpFpvlvRF4Hzg45IOprHgRUQ8IGlMj+bTgY0R\n8QyApMXAxRExj/0nsEtdXwp+emZSM7P+ayR4vBe4D5gaEduBI4E5KfShA9hUttyVtFUkaaSkLwCT\nJV1Xra3C52ZLWi1p9datW/vcWT/cZ2bWf42kWr8i6fvAcZLOSZp3ZdOtXvuxDbi6VluFzy2U9Dww\nfciQIaf2df9+uM/MrP8ayXb7APAApdHPjcmfH0uhD5uBo8uWRydtqUujtpsf7jMz679G7vlcC5wG\nPBQRb5M0gXQqHKyiNJoaSynoXAH8YQrbPYCk6cD0cePG9Xkb9Tzc51L3Zma9q/shU0mrIuI0SWsp\nFRn9jaTHI+KEuncm3QmcC7we+AVwQ0TcKundwKeAQcCiiLip0S/SiCzm8+nWnQ1XnpQwbPCgujLp\nzMxaWZoPmTYy8umSNAJYAnxH0ovAc43sLCJmVmlfAaxoZFt9kcbIpxZnw5mZ1Vb3PZ+IuCQitkfE\nx4CPAl+mNJV228hqPp9yzoYzM6utZvCRdLGkD5UtPwx8Bfg0cFGGfWtLLttjZlZbPSOfvwGWli0f\nTCnx4Fzggxn0KTNZTaNdztlwZma11XPPZ0hElD8E+u/JczXbJB2SUb8yERHLgGWdnZ1XVXq/Z5ba\n2ya8ge8/sbWhrDWXujczq61mtpukjRFR8Q69pKcj4k2Z9CwDZQkHVz311FP7vVcpS60nZ62ZWZHl\nPaXCw5IOGClI+jPgR2l0Ii+9JRxUylLrqTtrzczM+qeey25/CSyR9IfAj5O2Uynd+2mrbLfe1JuN\n5qw1M7P+qxl8ImIL8FZJbwe6HyhdHhH3Z9qznFWr2VZpPTMz659GnvO5PyI+k7zaMvD0lu1WKUut\nJ2etmZmlo6H5eNpdb/d8Kk0u974zj2l4sjkzM6utkfI6A15fJpczM7PGFWrkY2ZmraFQwSePCgdm\nZlZboYJPHoVFzcystkIFHzMzaw0OPmZmljsHHzMzy52Dj5mZ5c7Bx8zMcleo4ONUazOz1lCo4ONU\nazOz1lCo4GNmZq3BwcfMzHLn4GNmZrlz8DEzs9w5+JiZWe7aLvhIOlbSrZLuLmubIelLkr4u6Z3N\n7J+ZmdWWa/CRtEjSFknre7RPk/SkpI2S5va2jYh4JiJm9WhbEhFXAVcDl6ffczMzS1PeM5neBnwW\nuL27QdIg4BbgfKALWCVpKTAImNfj81dGxJZetn99si0zM2thuQafiHhA0pgezacDGyPiGQBJi4GL\nI2IecGE925Uk4Gbg3yLix+n12MzMstAK93w6gE1ly11JW0WSRkr6AjBZ0nVJ858D5wGXSbq6yudm\nS1otafXWrVtT6rqZmfVF3pfd+i0itlG6t1PetgBYUONzCyU9D0wfMmTIqRl20czMamiFkc9m4Oiy\n5dFJW+pc283MrDW0QvBZBRwnaaykIcAVwNIsduSq1mZmrSHvVOs7gQeB8ZK6JM2KiD3ANcB9wAbg\nroh4PIv9e+RjZtYa8s52m1mlfQWwIuv9S5oOTB83blzWuzIzs160wmW33HjkY2bWGgoVfMzMrDUU\nKvg44cDMrDUUKvj4spuZWWsoVPAxM7PWUKjg48tuZmatoVDBx5fdzMxaQ6GCj5mZtYZCBR9fdjMz\naw2FCj6+7GZm1hoKFXzMzKw1OPiYmVnuHHzMzCx3hQo+TjgwM2sNhQo+TjgwM2sNhQo+ZmbWGhx8\nzMwsdw4+ZmaWOwcfMzPLXaGCj7PdzMxaQ6GCj7PdzMxaQ6GCj5mZtQYHHzMzy52Dj5mZ5c7Bx8zM\ncufgY2ZmuWu74CPpWEm3Srq7rG2ipC9IulvSB5vZPzMzqy3X4CNpkaQtktb3aJ8m6UlJGyXN7W0b\nEfFMRMzq0bYhIq4G3gtMSb/nZmaWprxHPrcB08obJA0CbgHeBRwPzJR0vKRJku7t8Tqq2oYlXQQs\nB1Zk130zM0vDQXnuLCIekDSmR/PpwMaIeAZA0mLg4oiYB1zYwLaXAkslLQf+bzo9NjOzLOQafKro\nADaVLXcBZ1RbWdJI4CZgsqTrImKepHOBS4GDqTLykTQbmJ0s/qbnpb8Cez3wQrM70SJ8LPbxsdjH\nx2Kf8WltqBWCT0MiYhtwdY+2lcDKGp9bCCwEkLQ6Ijoz6mJb8bHYx8diHx+LfXws9pG0Oq1ttUK2\n22bg6LLl0UmbmZkNUK0QfFYBx0kaK2kIcAWwtMl9MjOzDOWdan0n8CAwXlKXpFkRsQe4BrgP2ADc\nFRGPZ9yVhRlvv534WOzjY7GPj8U+Phb7pHYsFBFpbcvMzKwurXDZzczMCsbBx8zMcle44NNIKZ92\nJ+loSd+X9BNJj0u6Nmk/UtJ3JD2V/HlE0i5JC5Jj85ikU5r7DdInaZCkRyTdmyyPlfRw8p2/niS9\nIOngZHlj8v6YZvY7bZJGJLUQn5C0QdJZRT0vJP1l8u9jvaQ7JQ0t0nlRqexZX84FSX+SrP+UpD+p\ntd9CBZ9qpXya26tM7QH+V0QcD5wJfCj5vnOB70XEccD3kmUoHZfjktds4PP5dzlz11JKbOn2ceCf\nImIc8CLQXTdwFvBi0v5PyXoDyaeBb0XEBOAkSsekcOeFpA7gw0BnRJwIDKKUcVuk8+I2epQ9o8Fz\nQdKRwA2UCgScDtzQHbCqiojCvICzgPvKlq8Drmt2v3L8/v8KnA88CYxK2kYBTyY/fxGYWbb+3vUG\nwovSM2TfA94O3AuI0pPrB/U8PyhlX56V/HxQsp6a/R1SOg6HA8/2/D5FPC/YV2HlyOTv+V5gatHO\nC2AMsL6v5wIwE/hiWft+61V6FWrkQ+VSPh1N6kuukssDk4GHgd+NiOeTt/4L+N3k54F+fD4F/A3w\nWrI8EtgepXR/2P/77j0WyfsvJesPBGOBrcBXkkuQX5Z0CAU8LyJiM/AJ4D+B5yn9Pa+hmOdFuUbP\nhYbPkaIFn0KSdCjwL8BfRMSvyt+L0q8pAz7fXtKFwJaIWNPsvrSAg4BTgM9HxGTgZfZdVgEKdV4c\nAVxMKSC/ETiEAy9BFVpW50LRgk/hSvlIGkwp8NwREfckzb+QNCp5fxSwJWkfyMdnCnCRpJ8Biyld\nevs0MEJSd43D8u+791gk7x8ObMuzwxnqAroi4uFk+W5KwaiI58V5wLMRsTUidgP3UDpXinhelGv0\nXGj4HCla8ClUKR9JAm4FNkTE/yl7aynQnY3yJ5TuBXW3/3GS0XIm8FLZ0LutRcR1ETE6IsZQ+nu/\nPyL+CPg+cFmyWs9j0X2MLkvWHxAjgYj4L2CTpO4Kxe8AfkIBzwtKl9vOlDQ8+ffSfSwKd1700Oi5\ncB/wTklHJKPJdyZt1TX7RlcTbqy9G/gp8DTwkWb3J+Pvejal4fJjwNrk9W5K16i/BzwFfBc4Mllf\nlLIBnwbWUcoAavr3yOC4nAvcm/x8LPAjYCPwDeDgpH1osrwxef/YZvc75WNwMrA6OTeWAEcU9bwA\nbgSeANYDX6U0NUthzgvgTkr3u3ZTGhXP6su5AFyZHJeNwJ/W2q/L65iZWe6KdtnNzMxagIOPmZnl\nzsHHzMxy5+BjZma5c/AxM7PcOfiYmVnuHHxswJP0kaRk/mOS1ko6I5lS4H/W+flfV/q5j30JSV8r\nWz5I0lYlUzykSdIYSTslrS3ylJpBAAADfUlEQVRr+7OkD+eWtX0oaTu/zu0OS47jbyW9Pu1+WzE4\n+NiAJuks4ELglIh4C6VyKpuAEUBdwacf+5aknv/GXgZOlDQsWT6fbEvVPB0RJ5ctTwIeBSYkfRwO\nfIBSodHH6tlgROxMtvnzlPtqBeLgYwPdKOCFiPgNQES8EBE/B24G3pT8Bj8fQNL7JP0oaftiMv9T\nRZKWSFqTjKhml7WPUWmywtspPTF/dIWPrwAuSH6eSekJ8163LekQScslParSpGeXV2qr43i8hVJt\nuwnJ8ocpPbH/WkT8oo7Pm6XCwccGum8DR0v6qaTPSfrvSftcklFBRMyRNBG4HJiS/Fb/KvBHvWz3\nyog4FegEPiypvKz+ccDnIuKEiHiuwmcXA1dIGkopGDzc4/1K254G/DwiTorSpGffqtJWy0TgLmCC\npBHJd/4hpUBplhsHHxvQIuLXwKmUZl3cCnxd0vsrrPqOZL1VyT2Sd1Cq71XNhyU9CjxEaXRzXNl7\nz0XEQ7306TFKk3fNpDQKqmfb64DzJX1c0u9HxEtV2qqSdDSwLSKeAY4C5gCfAd6cbKt83fFJsP6k\npDf2tl2zvjio9ipm7S0iXgVWAislraNUpXdlj9UE/HNEXFdre8nN+vMozWj5iqSVlApOdnu5jm4t\npTSJ2bmUTUZWbdsR8VNJp1AqDPsPkr4XEX9Xqa2XfU5iX5DZQWnkdDqlSfZ+XNaH11Oa5fdvgOHA\nJyS9PyJ+W8f3MquLRz42oCW/wZePSk4GnqP0n+9hZe3fAy6TdFTyuSMl/bcqmz0ceDEJDhOAM/vQ\ntUXAjRGxrkd7xW0no49XIuJrwHzglEptNfb5FvYFn/nANUlgLg9KAOdQmuvoEkozWP4LcEIfvqNZ\nVR752EB3KPCZ5P7GHkrl3mdHxDZJ/yFpPfBvyX2f64FvJxlqu4EPUQpUPX0LuFrSBkpz2Fe9xFZN\nRHQBCxrY9iRgvqTXkr59sEpbbyZRCiRERHlq9/HA42XLTwHvjIhPAkh6L/D/6v92ZrV5SgWzAUrS\nGErzFp3Yh89+kFIauICvRcS/VFjnZ5Tmc3mhfz21InLwMRugkgSDH1JKMji51voNbHcY8CDwBmBS\nRPwyrW1bcTj4mJlZ7pxwYGZmuXPwMTOz3Dn4mJlZ7hx8zMwsdw4+ZmaWOwcfMzPLnYOPmZnlzsHH\nzMxy9/8Bib9ex/9ZE6MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = []\n",
    "with open(sim3 + prop+'.data', 'rb') as f:\n",
    "    while True:\n",
    "        try:\n",
    "            data.append(pickle.load(f))\n",
    "        except EOFError:\n",
    "            break\n",
    "        \n",
    "    data = pd.DataFrame(data)\n",
    "\n",
    "\n",
    "### edit the code below to plot what you want\n",
    "\n",
    "plt.figure()\n",
    "plt.semilogy()\n",
    "plt.xlim(0,1000)\n",
    "plt.ylim(1e-13, 1e-7)\n",
    "# notice how you can retrieve the things we stored using the data array with data['key']\n",
    "plt.scatter(data['haloid'],data['sSFR'],label='all data') \n",
    "# notice also that you can filter the stored arrays based on the values of other stored arrays\n",
    "# using data['key'][True/False statement]\n",
    "#plt.scatter(data['mstar'][data['haloid']>=100], data['mgas'][data['haloid']>=100],label='haloid > 100')\n",
    "plt.legend()\n",
    "plt.xlabel('Stellar Mass [$M_\\odot$]')\n",
    "plt.ylabel('Gaseous Mass [$M_\\odot$]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0                      nan\n",
       "1                      nan\n",
       "2                      nan\n",
       "3                      nan\n",
       "4    0.0006731487161022437\n",
       "5                      nan\n",
       "6                      nan\n",
       "7                      nan\n",
       "8                      nan\n",
       "Name: MZgas, dtype: object"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data['MZgas']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "GMZ = []\n",
    "for i in range(5,13000):\n",
    "     if not h[i].properties['n_gas']==0 and not h[i].properties['n_star']==0:\n",
    "        halo = h.load_copy(i)\n",
    "        halo.physical_units\n",
    "        gem = numpy.sum(halo.gas['HI']*halo.gas['metals']*halo.gas['mass'])\n",
    "        mz = gem/numpy.sum(halo.gas['mass']*halo.gas['HI'])\n",
    "        print(i, mz)\n",
    "        GMZ = numpy.append(GMZ,mz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "SFH = []\n",
    "for i in range (5,13000):\n",
    "    if not h[i].properties['n_star']<100:\n",
    "        sfr, bins = pynbody.plot.stars.sfh(h[i], filename=None, massform=True, clear=False, legend=False, subplot=False, trange=False, bins=140)\n",
    "        SFR = numpy.sum(sfr [-10:])\n",
    "        print(SFR)\n",
    "        SFH = numpy.append(SFH,SFR)"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}