PBH_massfraction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Cython extension is already loaded. To reload it, use:\n",
      "  %reload_ext Cython\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "from importlib import reload \n",
    "import scipy.interpolate as interp\n",
    "import scipy.constants as const\n",
    "import scipy.integrate as integ\n",
    "\n",
    "import cython\n",
    "%load_ext Cython\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "mpl.rc('xtick', direction='in', top=True)\n",
    "mpl.rc('ytick', direction='in', right=True)\n",
    "mpl.rc('xtick.minor', visible=True)\n",
    "mpl.rc('ytick.minor', visible=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "# importing MS and HM modules. Reload modules if needed\n",
    "try :\n",
    "    reload(ms_hm.MS)\n",
    "except :\n",
    "    print(\"Did not reload MS module\")\n",
    "\n",
    "try :\n",
    "    reload(ms_hm.HM)\n",
    "except :\n",
    "    print(\"Did not reload HM module\")\n",
    "    \n",
    "try :\n",
    "    reload(ms_hm.QCD_EOS)\n",
    "except :\n",
    "    print(\"Did not reload QCD_EOS module\")\n",
    "    \n",
    "import ms_hm\n",
    "from ms_hm.QCD_EOS import *\n",
    "from ms_hm.MS import *\n",
    "from ms_hm.HM import *\n",
    "\n",
    "qcd = QCD_EOS()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEKCAYAAABHZsElAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAliklEQVR4nO3dfZxcZX338c83z2xIFhISIAQMeTAYQROJFFREBS0gC1pFk9K7CClbqKjctq9W2r7KrbU31ba2PlC4F0Fs9SYSBGExYn2Am6oIeSA8xAjsJsQsJNmQkE12N8lmd3/3HzMLw2R2MrM7s3Nm9/t+vc4rOdc55zq/2ezML9d1rrkuRQRmZmZJM6rSAZiZmeXiBGVmZonkBGVmZonkBGVmZonkBGVmZok0ptIBDCfHHHNMzJo1a0jutWPHDqZNmzYk9yqWYxu4JMfn2AYuyfElIbY1a9a8HBGHBOEEVUKzZs1i9erVQ3KvxYsXD9m9iuXYBi7J8Tm2gUtyfEmITdLmXOXu4jMzs0RygjIzs0RygqpS9fX1lQ6hX45t4JIcn2MbuCTHl+TY5KmOSmfx4sVR6b5cM7NqI2lNRCzOLncLyszMEskJyszMEskJyszMEsnfg7Ky+lXzy9y/7qVKh2FmZfZn75nLSVNrSlqnE5SV1b8/1Mzjm3Zx9MSxlQ7FzMroj858Q8nrdIKysmre0c5Fbz2er3xsYaVDMbMq42dQVjbtB7rZ2rafOdOOrHQoZlaFnKCsbDbt6ABwgjKzAXGCsrJp2rEXgLnTJ1Y4EjOrRk5QVjbNrR2MHiVOmuIEZWbFc4Kysmne0c4bptQwbox/zcyseP7ksLJp3tHObD9/MrMBcoKysuju6eWFlzuZ4+dPZjZATlBWFi2v7KOrp9cj+MxswJygCiRptqTbJN1d6ViqQfOOdsBDzM1s4MqaoCTdLqlV0jM5js2XtC5j2yPpuozjR0m6W9JvJW2QdFYZYjhf0rOSmiR9Ll89EbExIpYNJIaR6LUE5S4+MxuYcreg7gDOz3UgIp6NiIURsRA4HegE7s045avAgxFxCvBWYEPm9ZKmS5qUVTa30BgkjQZuAi4AFgBLJS2QdJqkB7K26QW9WnvVppc7mDJxHEfVjKt0KGZWpco6F19EPCJpVgGnngs0R8RmAEmTgXcDn0jX0wV0ZV1zDnCNpAsjYr+kq4APAxcWGMMZQFNEbEzfczlwSUTcCFxU2Cu0/uzuPMiUiU5OZjZwSXkGtQS4M2N/NrAD+JakJyR9U9Lr+ooiYgXwILBc0mXAlcDHirjnCcCWjP2WdFlOkqZKugVYJOn6XOe0tbVRX19PY2NjEWEMT51dPdSMG13pMMwswRobG6mvrweozXW84glK0jjgYmBFRvEY4G3AzRGxCOgADnlGFBFfBvYDNwMXR0R7MbfOURb9nRwROyPi6oiYk25lHaK2tpaGhgbq6uqKCGN46uzqdoIys7zq6upoaGgAaMt1vOIJitQzoLURsT2jrAVoiYjH0vt3k0pYryPpbOBUUs+ubijyvi3AiRn7MwGvrFciqRaUV3Mxs4FLQoJayuu794iIbcAWSfPTRecCv8k8R9Ii4FbgEuAKYIqkLxZx31XAPEknp1txS4D7B/YSLJu7+MxssMo9zPxO4FFgvqQWScvS5SslzZBUA7wfuCfH5Z8CvivpKWAh8L+zjtcAl0ZEc0T0ApcDmwuNISK6gWuBH5MaIXhXRKwf9Is2INXFN9EtKDMbhHKP4lvaT3nmSLup/ZyzDlicp+5fZu0fJNWiKiiG9LGVwMr+jtvAdR7o4Qi3oMxsEJLQxWfDTETQebCHieOdoMxs4JygrOQOdPfS0xseJGFmg+IEZSW3r6sHwIMkzGxQnKCs5Dq6ugEnKDMbHCcoK7nXWlDu4jOzgXOCspLrcBefmZWAE5SVXOerXXxuQZnZwDlBWcl1Hki1oDzM3MwGwwnKSq7zoLv4zGzwnKCs5DoPuIvPzAbPCcpKrtODJMysBJygrOQ8SMLMSsEJykqus6uHMaPEuDH+9TKzgfMniJWc14Iys1JwgrKSSy337u49MxscJygruY6uHmr8HSgzGyQnKCu5zgPd7uIzs0FzgrKSSz2DchefmQ2OE5SVXGdXDxPdgjKzQXKCspLzIAkzKwUnKCs5DzM3s1JwgrKSc4Iys1JwgrKS6+zqpma8u/jMbHCcoKykurp7OdgT1Ix1C8rMBscJykpqX99M5m5BmdkgOUFZSXUe7JvJ3C0oMxscJygrqY4DXgvKzErDCapAkmZLuk3S3ZWOJcn6uvgm+ntQZjZIZU1Qkm6X1CrpmRzH5ktal7HtkXRdxvEXJD2dPra6TDGcL+lZSU2SPpevnojYGBHLBhrHSNHR5S4+MyuNcv839w7gG8B/ZB+IiGeBhQCSRgMvAvdmnfbeiHg5V8WSpgP7ImJvRtnciGgqJIb0PW8C3g+0AKsk3Q+MBm7MquPKiGjt70XaazxIwsxKpawtqIh4BNhVwKnnAs0RsbmI6s8B7pM0AUDSVcDXiojhDKAp3TLqApYDl0TE0xFxUdZWUHJqa2ujvr6exsbGIl7G8OIWlJkVqrGxkfr6eoDaXMeT8gxqCXBnVlkA/yVpjaT67AsiYgXwILBc0mXAlcDHirjnCcCWjP2WdFlOkqZKugVYJOn6XOfU1tbS0NBAXV1dEWEML51dHiRhZoWpq6ujoaEBoC3X8Yr3w0gaB1wMZH/ovzMiXkp35f1E0m/TraFXRcSXJS0HbgbmRER7MbfOURb9nRwRO4Gri6h/ROo80NeCqvivlplVuSS0oC4A1kbE9szCiHgp/WcrqWdTZ2RfKOls4NT08RuKvG8LcGLG/kzgpSLrsCydB92CMrPSSEKCWkpW956kiZIm9f0d+ADwTNY5i4BbgUuAK4Apkr5YxH1XAfMknZxuxS0B7h/wqzAAOg/0MEowfkwSfrXMrJqVe5j5ncCjwHxJLZKWpctXSpohqYbUKLp7si49FviFpCeBx4EfRsSDWefUAJdGRHNE9AKXA4cMsugvhojoBq4FfgxsAO6KiPWleeUjV2qxwjFIuXpQzcwKV9YHBRGxtJ/yCzN2p+Y4vhF462Hq/mXW/kFSLaqCYkgfWwmszHcfK077gYNM9BBzMysB98NYSbXtO0jtEWMrHYaZDQNOUFZSTlBmVir99sVI+my+CyPiK6UPx6rd7s6DnDilptJhmNkwkO9hwaT0n/OBt/PaCLc64JGcV9iIt8ctKDMrkX4TVER8HkDSfwFv65vzTtL/AlYMSXRWdXbvO8hRTlBmVgKFPIM6CejK2O8CZpUlGqtqXd29dHb1uAVlZiVRyHjg/wQel3QvqamAPkyO2cnN2vYdBKC2xgnKzAbvsAkqIv5B0o+As9NFV0TEE+UNy6rRqwnKLSgzK4F8o/imZOy+kN5ePRYRhSyjYSOIE5SZlVK+FtQaUl16/c36PbssEVnVatuXelTpBGVmpZBvFN/JQxmIVb++FtRRNeMqHImZDQcFTZom6WLg3endhyPigfKFZNVqd6e7+MysdA47zFzSPwKfAX6T3j4j6cZyB2bVp68FNXmCJ4s1s8Er5JPkQmBhekkLJH0beIJDV8C1Ea5t30EmjR/DmNGe4tHMBq/QT5KjMv5eW4Y4bBho6zzIZHfvmVmJFNKCuhF4QtJDpEb0vRu3niyHtn0HOcpf0jWzEinki7p3SnqY1ISxAv4qIraVOzCrPrs9UayZlVChXXyjgJeBV4A3Snr3Yc63EcgtKDMrpcO2oCR9Cfg4sB7oTRcHXnLDsnixQjMrpUKeQX0ImB8RB8oci1WxiPAgCTMrqUK6+DYC/tSxvPYf7KWrp5ejjvAsEmZWGvkmi/06qa68TmCdpJ8Br7aiIuLT5Q/PqsVuz8NnZiWWr4tvdfrPNby23LtZTp7J3MxKLd9ksd/OLpN0NHBiRDxV1qis6vTNw+dRfGZWKoXMxfewpMnp9aGeBL4l6SvlD82qiVtQZlZqhQySqI2IPcAfAN+KiNOB88obllUbJygzK7VCEtQYSccDHwO8zIbl1Na31Ia7+MysRApJUF8Afgw0R8QqSbOB58sbVvJImi3pNkl3VzqWJGrbd5BRgiPHeakNMyuNwyaoiFgREW+JiGvS+xsj4iOFVC7pdkmtkp7JcWy+pHUZ2x5J12WdM1rSE5IG3HI7TAznS3pWUpOkz+WrJ/26lw00juFuZ0cXR9eMY9QoVToUMxsmChkk8UZJP+v7gJf0Fkl/W2D9dwDn5zoQEc9GxMKIWAicTur7VvdmnfYZYEM/cU2XNCmrbG6hMUgaDdwEXAAsAJZKWiDpNEkPZG3T87xGA7bv2c+xkydUOgwzG0YK6eK7ldTyGgcB0kPMlxRSeUQ8Auwq4NRzSXUhbu4rkDQT+CDwzX6uOQe4T9KE9PlXAV8rIoYzgKZ0y6gLWA5cEhFPR8RFWVtrAa9hRNvWtp/jap2gzKx0CklQNRHxeFZZd4njWALcmVX2b8Bf8toEta8TESuAB4Hlki4DriQ1kKNQJwBbMvZb0mU5SZoq6RZgkaSc62G1tbVRX19PY2NjEWEMD25BmVmxGhsbqa+vh34Wwi3kifbLkuaQmvYISR8FtpYqQEnjgIvJWARR0kVAa0SskfSe/q6NiC9LWg7cDMyJiPZibp2ryjz32glcna/C2tpaGhoaighheDjQ3cPOji6Oc4IysyLU1dVRV1fHrbfe2pbreCEtqE8C/wc4RdKLwHUc5oO6SBcAayNie0bZO4GLJb1AquvtfZK+k32hpLOBU0k9u7qhyPu2ACdm7M8EXiqyDgNa96SmaDyudnyFIzGz4SRvgkoPJLgmIs4DpgGnRMS7Mp8VlcBSsrr3IuL6iJgZEbNIdf/9PCL+KCu2RaSej10CXAFMkfTFIu67Cpgn6eR0K24JnnNwQLbv2Q/gLj4zK6m8CSoiekiNsCMiOiJibzGVS7oTeBSYL6lF0rJ0+UpJMyTVAO8H7hlA7DXApRHRHBG9wOXAIYmzvxgiohu4ltR3vDYAd0XE+gHEMeJtc4IyszIo5BnUE5LuB1YAHX2FEXHYpBIRS/spvzBjd+ph6ngYeDhH+S+z9g+SalEVFEP62EpgZb772+Fta0slKD+DMrNSKiRBTQF2Au/LKAsG1uqxYWj7nv2MGzPKM5mbWUkdNkFFxBVDEYhVr+17DnDc5AlInkXCzErnsAlK0jTgKmBW5vkRcWX5wrJqsm3PfnfvmVnJFdLFdx/w38BPgZ7yhmPVaPue/bxl5lGVDsPMhplCElRNRPxV2SOxqhQRbGvbzwcW+DtQZlZahXxR9wFJFx7+NBuJ2vYd5EB3r4eYm1nJ9duCkrSX1Gg9AX8t6QCpCWMFRERMHpoQLcn6vgPliWLNrNT6TVARMam/Y2Z9/B0oMyuXQtaD+lkhZTYyeZojMyuXfF18E4CJwDGSjua12b8nAzOGIDarAtvaUhPFTp/sQRJmVlr5RvH9KamZy2cAa3gtQe0htRKtGS/s7OD42gmMHzO60qGY2TCT7xnUV4GvSvpURHx9CGOyKtLU2s7c6UdWOgwzG4YO+wzKycn609sbNO9oZ840JygzK71CvgdlltPWPfvp7OpxC8rMysIJygasqbUdwAnKzMoi3yi+t+W7MCLWlj4cqyZOUGZWTvlG8f1LnmPB69eHshGoqbWdo2rGMnXiuEqHYmbDUL5RfO8dykCs+jS3tjN32pFeB8rMyqKQ2cyRdCqwAHh1uoCI+I9yBWXVoWlHOx9YcGylwzCzYaqQBQtvAN5DKkGtBC4AfgE4QY1guzq62NXR5edPZlY2hYzi+yhwLrAtvfz7WwHPazPC9Q2QmOMEZWZlUkiC2hcRvUC3pMlAKzC7vGFZ0r06gs9f0jWzMinkGdRqSUcBt5Kak68deLycQVnyPdWym0kTxnDCUUdUOhQzG6YOm6Ai4s/Sf71F0oPA5Ih4qrxhWdI9tmkXZ8yawqhRHsFnZuVRyHpQH5ZUCxARLwC/k/ShMsdlCda6Zz+bXu7gzNlTKx2KmQ1jhTyDuiEi2vp2ImI3cEPZIrLE+/WmXQD83uwpFY7EzIazQhJUrnMK+v6UDU+/3riTI8ePYcHxkysdipkNY4UkqNWSviJpjqTZkv6V1GAJG6Ee27iTxbOOZsxozzVsZuVTyCfMp4Au4HvACmA/8MlyBpVE6eR8m6S7Kx1LJe3Ye4DmHX7+ZGblV8iChR0R8bmIWBwRp0fE9RHRUUjlkm6X1CrpmRzH5ktal7HtkXRd+tgESY9LelLSekmfL/qVFRbD+ZKeldQk6XP56omIjRGxbKBxDBePbdoJwO+d7OdPZlZe+Zbb+LeIuE5SI6nZy18nIi4uoP47gG+QY1qkiHgWWJi+12jgReDe9OEDwPsiol3SWOAXkn4UEb/OiG86qS8R780omxsRTYXEkL7nTcD7gRZglaT7gdHAjVl1XBkRrQW83mHvgSe3MmXiOE49obbSoZjZMJdvsMN/pv/854FWHhGPSJpVwKnnAs0RsTl9XZD6QjDA2PSWnSTPAa6RdGFE7Jd0FfBh4MICYzgDaIqIjQCSlgOXRMSNwEWFvL5sbW1t1NfXU1dXR11d3UCqSLQdew/w0w3bueKdsxjr509mNkiNjY00NjYC5Pwfb7+fMhHRNxBiYUT8v8yNdMunhJYAd2YWSBotaR2pqZV+EhGPZcW3AngQWC7pMuBK4GNF3PMEYEvGfku6LCdJUyXdAiySdH2uc2pra2loaBiWyQng+2tb6O4NPv72kyodipkNA3V1dTQ0NAC05TpeyH+DL89R9olBxPQ6ksYBF5MagPGqiOiJiIXATOCM9JIfZJ3zZVKDNm4GLo6I9uxz8t06R9khXZkZ99oZEVdHxJx0K2tEiQi+t2oLb591tGcwN7Mhke8Z1FLgD4GT089m+kwCdpYwhguAtRGxPdfBiNgt6WHgfOB1Ax0knQ2cSurZ1Q3AtUXctwU4MWN/JvBSEdePKI9u3Mmmlzu49r1zKx2KmY0Q+Z5B/QrYChzD65d/3wuUci6+pRzavTcNOJhOTkcA5wFfyjpnEakJbD8IbAK+I+mLEfG3Bd53FTBP0smkBmgsIZWQLUtPb/CPP/ot0yaN58LTjq90OGY2QuR7BrU5Ih6OiLOynkGtjYjuQiqXdCfwKDBfUoukZenylZJmSKohNYrunqxLjwcekvQUqUTyk4h4IOucGuDSiGhOLwdyObC50BjSr+Fa4MfABuCuiFhfyOsaab63agtPtbTxNxe+iSPGja50OGY2Qig1YC7HAekXEfEuSXt5/bMZkRpo53lusixevDhWr15d6TBK6pWOLt77Lw/zxmMn8b36M5E8e7mZlZakNRGxOLu83y6+iHhX+s9J5QzMkqv9QDdXfnsVHQe6+cIlb3ZyMrMhlXcUn6RRuWZgsOGv/UA3V35rFU+1tPH1pW/jlOPcYDazoZV3VvKI6E1PN3RSRPxuqIKyylr9wi7+513rePGVfXx1ySLOP/W4SodkZiNQIctmHA+sl/Q48OocfAVOdWRVYlvbfn68fhs/emYrj23axcyjj+CuPz2LxbM8556ZVUYhCWrAE7VaMu3q6GLD1j1s2LqH5h0d/GbrHp7cshuAedOP5NPvm8dV757NkeO97JeZVU6+L+rOBY5NT22UWf5uUt8bsoSLCH67bS8//20rv9vZydY9+3lu21627dn/6jlH14xlzrQj+YsPvJHzTz3es0SYWWLk+y/yvwF/naO8M31seE44N0z89/M7+PsHfsNz29uR4Jgjx3Ps5PGcNWcqC46fzIIZkznluElMPXJ8pUM1M8spX4KaFRGHzBgREasLnKHcKqDjQDfX3/M09z/5EicfM5F/+PCp/P6bj+MYJyIzqzL5EtSEPMeOKHUgNngv7d7Hsm+v5rnte7nuvHlc8545jB/jmR/MrDrlS1CrJF0VEbdmFqanClrTzzVWIa90dPHRm3/F3v3d3P6Jt3POG6dVOiQzs0HJl6CuA+5Nr7XUl5AWA+NILQxoCRERXH/P0+xoP8D3r3kHb5l5VKVDMjMbtHxTHW0H3iHpvaSWtAD4YUT8fEgis4KtWN3Cg+u38dcXnuLkZGbDxmG/6BIRDwEPDUEsNgA72w/whQd+w1mzp/In75pd6XDMzEqmkBV1LcFueqiZzq5u/v5Db2bUKE/mambDhxNUFduyq5Pv/Hozl55+InOne9J5MxtenKCq2L/+9DkkuO798yodiplZyTlBValNL3fwgyde5I/PegPH1/praWY2/DhBVambH25i7OhRXPVuD4wws+HJCaoKtbzSyT1rX2TpGScxfVK+CT/MzKqXE1QVanhkIxLUu/VkZsOYE1SVeXH3PpY/voWPvG0mM47ysyczG76coKrM1376PACfOtcj98xseHOCqiLNO9pZsWYLl515Eie49WRmw5wTVBX51588x4Sxo/nke+dWOhQzs7JzgqoSm3d28MOnt/LHZ83y4oNmNiI4QVWJb/3yBcaMEle8c1alQzEzGxJOUFWgbd9B7lq9hbq3zuDYyf7ek5mNDE5QBZI0W9Jtku4e6nvf+fjv6OzqYdm7Th7qW5uZVUxZE5Sk2yW1Snomx7H5ktZlbHskXZc+dqKkhyRtkLRe0mfKFMP5kp6V1CTpc/nqiYiNEbFsoHEMVETw3cc2c+bsKbx5Ru1Q397MrGLK3YK6Azg/14GIeDYiFkbEQuB0oBO4N324G/jziHgTcCbwSUkLMq+XNF3SpKyyXMPbcsYgaTRwE3ABsABYKmmBpNMkPZC1TS/4FZfY2t/tZsuufVx6+omVCsHMrCLKmqAi4hFgVwGnngs0R8Tm9HVbI2Jt+u97gQ3ACVnXnAPcJ2kCgKSrgK8VEcMZQFO6ZdQFLAcuiYinI+KirK21kNfb1tZGfX09jY2NhZxekPvXvcj4MaP4wJuPLVmdZmZJ0NjYSH19PUDO7qGkPINaAtyZ64CkWcAi4LHM8ohYATwILJd0GXAl8LEi7nkCsCVjv4VDk2BmHFMl3QIsknR9rnNqa2tpaGigrq6uiDD6193TywNPbeW8Nx3LpAljS1KnmVlS1NXV0dDQANCW6/iYoQ3nUJLGARcDh3zoSzoS+D5wXUTsyT4eEV+WtBy4GZgTEe3F3DpHWfR3ckTsBK4uov5B+2XzTnZ2dHHxwhlDeVszs0RIQgvqAmBtRGzPLJQ0llRy+m5E3JPrQklnA6eSenZ1Q5H3bQEyH+zMBF4qso6yum/di0yaMIb3zJ9W6VDMzIZcEhLUUrK69yQJuA3YEBFfyXWRpEXArcAlwBXAFElfLOK+q4B5kk5Ot+KWAPcPIP6yiAgeeW4H7ztlOuPHjK50OGZmQ67cw8zvBB4F5ktqkbQsXb5S0gxJNcD7gewW0juB/wG8L2MY+oVZ59QAl0ZEc0T0ApcDmwuNISK6gWuBH5MahHFXRKwv0UsftOdb23m5vYt3zjmm0qGYmVVEWZ9BRcTSfsozk83UHMd/Qe5nRJnn/DJr/yCpFlVBMaSPrQRW5rtPpTzavBOAs+Yc8uMxMxsRktDFZzn8qvllZh59BCdOqal0KGZmFeEElUC9vcFjm3Zx1my3nsxs5HKCSqAN2/awu/Mg75jrBGVmI5cTVAK9+vxptgdImNnI5QSVQL/euIuTj5nIcbVeWsPMRi4nqAR6dvse3jxjcqXDMDOrKCeohNnX1UPLK/uYN33S4U82MxvGnKASpnlHOxEw79gjKx2KmVlFOUElTFNrar7budOdoMxsZHOCSpjnW/cyepSYNXVipUMxM6soJ6iEaWptZ9bUGsaN8T+NmY1s/hRMmOdb2z1AwswMJ6hE6eruZfPOTj9/MjPDCSpRXtjZQU9veASfmRlOUIny/HaP4DMz6+MElSDPt+5FgjnTnKDMzJygEqSptZ0Tj65hwlgv8W5m5gSVIE2t7cxz956ZGeAElRjdPb1s3NHh509mZmlOUAnxu12ddPX0OkGZmaU5QSVE3xx88471l3TNzMAJKjGeTyeoOdM8B5+ZGThBJUZTazvH105g0oSxlQ7FzCwRnKASoqm13c+fzMwyOEElQG9vpIeY+/mTmVkfJ6gEeHH3PvYd7HELyswsgxNUAjTt6BvB5wRlZtbHCSoBmvomifUcfGZmr3KCKpCk2ZJuk3R3qet+vnUvxxw5jqMnjit11WZmVausCUrS7ZJaJT2T49h8Sesytj2Srivk2hLGcL6kZyU1SfpcvnoiYmNELBtMLP35vZOncvlZs8pRtZlZ1RpT5vrvAL4B/Ef2gYh4FlgIIGk08CJwbyHXpq+ZDuyLiL0ZZXMjoqmQGNL3vAl4P9ACrJJ0PzAauDGrjisjorXfVzlIHzl9ZrmqNjOrWmVtQUXEI8CuAk49F2iOiM1FXHsOcJ+kCQCSrgK+VkQMZwBN6ZZRF7AcuCQino6Ii7K2gpJTW1sb9fX1NDY2FnK6mdmI1tjYSH19PUBtruNJeQa1BLizmAsiYgXwILBc0mXAlcDHiqjiBGBLxn5LuiwnSVMl3QIsknR9rnNqa2tpaGigrq6uiDDMzEamuro6GhoaANpyHS93F99hSRoHXAzk/NDPJyK+LGk5cDMwJyLai7l1rirz3GsncHWRIZqZ2QAloQV1AbA2IrYXe6Gks4FTST27uqHIy1uAEzP2ZwIvFRuDmZmVRxIS1FKK7N4DkLQIuBW4BLgCmCLpi0VUsQqYJ+nkdCtuCXB/sXGYmVl5lHuY+Z3Ao8B8SS2SlqXLV0qaIamG1Ci6ewq9NkMNcGlENEdEL3A5sLnQeiKiG7gW+DGwAbgrItaX5pWXX7rfNpEc28AlOT7HNnBJji/JsSmi38cuVqTFixfH6tWrh+peDNW9iuXYBi7J8Tm2gUtyfEmITdKaiFicXZ6ELj4zM7NDuAVVQpJ2kKObsUyOAV4eonsVy7ENXJLjc2wDl+T4khDbGyJiWnahE5SZmSWSu/jMzCyRnKDMzCyRnKDMzCyRnKDMzCyRnKCGCUkfknSrpPskfaDS8WSTNFHSGkkXVTqWTJJGSfoHSV+XdHml48kk6SRJ96fXNMu7XtlQybVwZ/rf9tvp37/LEhhfIt4b/S16moT3Rj8/t4q/N5ygEqy/xRZzLbQYET+IiKuATwAfT1JsaX8F3FXuuAYQ2yWkZrE/SGp+xiTF9kbghxFxJbAgCTH1s3DnHwB3p3//Lk5afOV8b5TgZwdlem+UILYhfW/k4gSVbHcA52cW6LWFFi8g9aG1VFLmh9ffpo8nJjZJ5wG/AYqeELjcsQHzgUcj4rPANQmL7QlgiaSfAw8lJKZcZvLa0jU9CYyvTzneG4OKrczvjUHFxtC/Nw7hBJVg/Sy2mHOhRaV8CfhRRKxNUmzAe4EzgT8ErpJUiYUy+4utBXglfU45PlwHE9sVwA0R8T7ggwmJKZcWUkkKyvCZMtj4yvneKMHPrmzvjRL9uw7ZeyMXJ6jq099Ci58CzgM+KqlS61bljC0i/iYirgP+L3BrenLfRMRGaqLi35f0deCRCsQF/cf2IPBppRbKfCEJMSn3wp33AB+RdDMwVMtJFxPfUL83Co6tAu+NYv9dK/reqPiChVa0nAstRsTXyLHk/RDLuwhkRNwxdKEcor+fWyeQ67nAUOovtmeAjw51MGn9xXTIwp0R0UGqtTeUiolvqN8bBceWcfCOskb0mmJ+bhV/b7gFVX2SvNCiYxuYJMaWxJgyJTk+x1YiTlDVJ8kLLTq2gUlibEmMKVOS43NspRIR3hK6kVppeCuvDfNcli6/EHgOaAb+xrE5tuEWU7XE59jKu3k2czMzSyR38ZmZWSI5QZmZWSI5QZmZWSI5QZmZWSI5QZmZWSI5QZmZWSI5QZmZWSI5QZmZWSI5QZkNkKT2rP1PSPpGkXX8qrRRJY+kP5UUks7JKLs2XXZeJWOzZHOCMquA9BpFoyLiHZWOZQi8BXgKeBOApBpSs2TvAJ6uYFyWcE5QZmUg6bOSnklv16XLZknaIOnfgbXAiX2tMElXS1qX3jZJeqiAem6VtF7Sf0k6IkcMsyT9VtI309d/V9J5kn4p6XlJZ2Sc+wNJa9L11afLJkr6oaQn09d/PFdZAT+O00jNC3dKev/TwAqgNyKGapVlq0aVngzQm7dq3UitMrouY/sd8A3gdFItg4nAkcB6YBEwC+gFzsyooz2rzrHAfwN1h6mnG1iYvuYu4I9yxNd33mmk/jO6Brid1JpAlwA/yDh3SvrPI4BngKnAR0gtotd3Tm2usgJ+TjuA2aQWYKwltZT9e4CfVvrf0FuyN7egzAZuX0Qs7NuAv0uXvwu4NyI6IqKd1MqkZ6ePbY6IX+ep86vAzyOi8TD1bIqIdem/ryGVjHLZFBFPR2ql1vXAzyIiSCW+zGs+LelJ4Nek1gualz7nPElfknR2RLT1U9YvSScCOyNiIzAd+Evg68AbSXX7Iem+jPNXSBqdr04bOZygzEov16qlfTr6vUj6BPAG4PMF1HMg4+899L86duZ5vRn7vX3XSHoPqSXRz4qIt5Jq4UyIiOd4rRV3o6S/y1WWJ0ZIPX/qe860Fzgf+DapVt3T6QS2NeP8URHRc5g6bYRwgjIrvUeAD0mqkTQR+DCpbrt+STod+AtSXXW9A61ngGqBVyKiU9IpwJnpmGYAnRHxHeCfgbflKkuf+zNJJ+So+zReS1D/BFybTkCnkWpBnQ4skHSLpG+T4NVdbej1978uMxugiFgr6Q7g8XTRNyPiCUmz8lx2LTAFeEgSwOqI+JMB1DMQDwJXS3oKeJZUNx+kksg/SeoltejdNbnKJI0C5gK7ctR9GvB9gIh4IKN8AfAb4EPAn0fEKkkfBKaV+LVZFfOChWY2KJJOBa6MiM8O4Np7gY9HRJekzwN3R4SHnhvgBGVmZgnlZ1BmZpZITlBmZpZITlBmZpZITlBmZpZITlBmZpZITlBmZpZITlBmZpZITlBmZpZI/x/VOHXTyN7PpgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#obtaining deltac and density values\n",
    "combineddata = np.loadtxt('initial_critical_data_qcdw.dat')\n",
    "\n",
    "combineddata = combineddata[combineddata[:, 0].argsort()]\n",
    "deltac = (combineddata[:,1] + combineddata[:,2]) / 2\n",
    "deltac_up = combineddata[:,1]\n",
    "deltac_low = combineddata[:,2]\n",
    "rho = combineddata[:,0]\n",
    "\n",
    "T = qcd.TofRho(rho)\n",
    "Mhorizon = qcd.MH(T)\n",
    "\n",
    "Order = Mhorizon.argsort()\n",
    "Mhorizon = Mhorizon[Order]\n",
    "deltac = deltac[Order]\n",
    "deltac = interp.InterpolatedUnivariateSpline(Mhorizon, deltac, ext = 3)\n",
    "deltac_up = deltac_up[Order]\n",
    "deltac_up = interp.InterpolatedUnivariateSpline(Mhorizon, deltac_up, ext = 3)\n",
    "deltac_low = deltac_low[Order]\n",
    "deltac_low = interp.InterpolatedUnivariateSpline(Mhorizon, deltac_low, ext = 3)\n",
    "\n",
    "MH_grid = np.logspace(-10,10,500)\n",
    "\n",
    "# plt.loglog(MH_grid, deltac(MH_grid))\n",
    "# plt.loglog(MH_grid, deltac_up(MH_grid))\n",
    "# plt.loglog(MH_grid, deltac_low(MH_grid))\n",
    "plt.loglog(rho,deltac(rho))\n",
    "plt.xlabel(\"Horizon mass, $M_H$\");\n",
    "plt.ylabel(\"Critical threshold\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fc5931c3b80>]"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAD5CAYAAADV5tWYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAj40lEQVR4nO3dd3RVVd7/8ff3ptNCCWAIJZSA0hSM9KooiEYGdRSwC2JDdJgGz4yjU3x0dMQKYhgZrGBDJSIiNkCKQkSkQ0CBUEPvAsn+/WH0l+EheEPuzbn35vNaK2t5DvecfPbiri/bffbZ25xziIhIePN5HUBEREpPxVxEJAKomIuIRAAVcxGRCKBiLiISAVTMRUQiQLRXvzgpKcmlpqZ69ev9lpeXR82aNb2OETSR3D61LXxFcvtK07bs7OydzrlTXuxZMU9NTWXRokVe/Xq/paenh0XOMxXJ7VPbwlckt680bTOzDcX9mYZZREQigIq5iEgE8KyY79u3j6FDh5KVleVVBL8MHTrU6whBFcntU9vCVyS370zalpWV9dN1icV9xrxamyU9Pd1F6piYiEgwmFm2cy79VH+mYRYRkQigYi4iEgHCrpjnF2jJXhGRk4VVMT92ooD+Y+cy5rMcjp0o8DqOiEjICKvZLEeO51OvWgUem7GajGe+4OuNe4KYUEQkNETsbJaPV2znL+8tY+v+o9zQoQG/792MyvExAU4oIhJaIm42S6/mtfloRHdu6pjKyws2cPHo2Xy0fJvXsUREPBOWxRygUlw0D17Rgnfu6kzVCjEMfTmbO17OZvv+o15HExEpc2FbzH9yXr2qZN3ThT/0acZnq3fQ6/FZvLxgAwWa9SIi5UjYF3OAmCgfd/Vowoz7utG6XiL3v7uMq8fNY832A15HExEpEwEv5mbW3MzeMLPnzOzqQN//dFKTKvLK4Pb869fn8t3OQ1z29BxGf7Sao8fzyzKGiEiZ86uYm9kEM9thZstOOt/HzFabWY6ZjSw8fSnwjHPuTuDGAOf1JytXn1+Xj0d0J6N1HZ7+NIe+T81hwfpdZR1FRKTM+Nsznwj0KXrCzKKAMfxYvJsDA82sOfAyMMDMHgNqBC5qydSoFMfoa8/j5cHtOFHgGJC5gD++9S17Dx/zKpKISND4Vcydc7OB3SedbgfkOOfWO+eOAZOBfs65Hc65u4GRwM6Apj0DXdNqMuO+btzRvTFvfZ1Lr9GzmLpkC17NrxcRCYbSjJmnAJuKHOcCKWaWamaZwEvAY8VdnJeXR3p6+s8/mZmZpYhyegmxUYy89GymDutMnaoJDJ+0mFsnLiR3z+Gg/U4RkUDIzMz8uU4CScV9zu83QM0sFXjfOdey8PjXQG/n3JDC4xuAds65e/y5n1frmecXOF6c9z3/+mg1zsFvL2nKzZ1SiY6KiIk9IhLBgvUGaC5Qr8hxXWCLvxd7tdNQlM+4tUtDZo7oTsfGNfjHtJX0HzuPZZv3lWkOERF/BXRtllP0zKOBNcBFwGZgITDIObfcn/uFwk5DzjmmLd3Kg1NXsOfwMYZ0aci9vdKoEBvtaS4RkVMpdc/czCYB84FmZpZrZoOdcyeAYcAMYCXwhr+FHEJjD1Az4/LWdfhkRHeuSa/L87PXc8kTs5m1Js+zTCIiJ4vYVROD5cv1uxj1zlLW5x2i33l1uP/y5iRVivM6logIEKKrJoZCz/xk7RvV4IPhXRl+URofLN1Kr9GzeHPRJk1jFBFPqWdeCmu3H2DUlKUs2rCHTo1r8FD/VjRMquh1LBEpx9QzPwNptSvzxu0deah/S5bm7qP3k7O1XZ2IeEI98wDZvv8oD05dzvRl22hWuzIPX9WKtvWreR1LRMqZkOyZh5PaVeJ57vrzybzhfPYdOc5Vz83jgfeWceDoca+jiYgAGmYpkUtanMXMEd24qWMqLxVuVzdzxXavY4lIhNMwSxAt3riHUVOWsmrbAS5teRYPXtGC2lXivY4lIhFMwyxB0KZ+tZ+3q/t01Y/b1b2i7epExCMaZimFn7ar+/C+brSqm8if313Gr5+fz1ptVyciAaRhljLknOPtrzfzj2krOPTDCe7s0YS7ejQmPibK62giEiE0zFIGim5Xd1mrZJ7+ZC19n57Dl9quTkTKgIp5gCVViuPJAW148dZ2HM8v4NrMBYx8+1v2HdY0RhEJHhXzIOne9Mft6m7v1og3s3O5aPQssrRdnYgEiR6ABlGF2GhG9T2H9+7uTHJiPPdMWszgFxdpuzoRKRE9AA0hJ/ILmDjvex7/aA1m8NtLmnFzp1SifOZ1NBEJE3oAGgKio3wM6dqImSO60b5hdf7+/gr6j53L8i3ark5ESk/FvIzVrVaBCTdfwDMD27Bl7xGueHYuD09fyZFj+V5HE5EwpmLuATMj49w6fDyiO1e3rcvzs9bT+8nZzFmr7epE5MyomHuoaoVY/nl1aybd1oFon3HDC18x4vVv2H3omNfRRCTMaDZLCOjYuAYf3NuVey5swtQlW7jo8c+Z8nWupjGKCKDZLGFp9bYDjJryLV9v3EuXJkk81L8lDWpouzoR0WyWsNLsrMq8dUcn/t6vBd9s2kvvJ2czbtY6TuRruzoRKZ6KeQjy+YwbOqYyc0Q3uqbV5JHpq7ji2bkszdU0RhE5NRXzEJacmMD4G9MZd31bdh78gX5jvuChaSs4fOyE19FEJMSomIeBPi2TmTmiO9deUJ/xc76j95Oz+WLtTq9jiUgICXgxN7P6ZjbVzCaY2chA37+8SkyI4eErWzF5aAeifT6uf+FLfv/mEq3GKCKAn8W8sDDvMLNlJ53vY2arzSynSOFuCkxzzt0KNA9w3nKvQ6MaTL+3K3f2aMyUxZu5aPQsPly21etYIuIxf3vmE4E+RU+YWRQwBriUH4v2QDNrDiwGBpjZp8BngYsqP4mPieKPfc7mvbs7U7tKHHe88jV3vZpN3oEfvI4mIh7xq5g752YDu0863Q7Icc6td84dAyYD/YBbgAeccxcClwUyrPy3limJvHt3Z37fuxkfr9jBxU/M4u1svWwkUh6VZsw8BdhU5Di38NyHwHAzGwd8X9zFeXl5pKen//yTmZlZiijlV0yUj7t7NuGDe7vSuGYlfvvmEm6c8BWbdmvNdJFIkJmZ+XOdBJKK+5zfb4CaWSrwvnOuZeHxr4Hezrkhhcc3AO2cc/f4cz+9ARp4BQWOlxds4NEPV1Hg4L5eadzapSExUZq0JBIJgvUGaC5Qr8hxXWCLvxdrbZbA8/mMmzql8tGI7nRqXIOHp6/i8qe/YOH3J4+QiUg4CejaLKfomUcDa4CLgM3AQmCQc265P/dTzzz4Plq+jb9mrWDz3iP0O68Ov+/djLrVKngdS0TOUKl75mY2CZgPNDOzXDMb7Jw7AQwDZgArgTf8LeSgnnlZuKTFWcwc0Y1hPZvw4bJtXPj4LB6evlJz00XCjFZNlJ9t3XeEf81Yw5TFuVSMjebGjg0Y3KUhNSrFeR1NRPx0up65Z8U8LS3N9ezZk4yMDDIyMjzJUB6t2rafZz7N4YOlW4mPjuKa9Lrc1CmVRjUreR1NRIqRlZVFVlYW48ePz3HOpZ3qM+qZl1M5Ow4w9vN1ZC3ZwvF8R49mNbmpYyrdmtYkymdexxORU1DPXIq148BRJn25iVe+3EDegR9IqZrAwHb1uCa9HrWqxHsdT0RQz1xK4NiJAj5euZ1Xv9zA3JxdRPuM3i3O4roO9enYqAZm6q2LeO10PfPosg4joSk22kffVsn0bZXM+ryDvPblRt7MzmXa0q2k1arErV0a0r9NCvExUV5HFZFT0DCLFOvo8XyylmzhP3O/Z8XW/VSvGMsNHRpwc6dUqlWM9TqeSLmhYRYJCOccC9bv5oUv1vPxyh1UiI3i+g4NGNK1IbUqa1xdpKyE5ANQFfPwtHrbAcZ+nkPWki3ERPm4qVMqt3drpPnqImUgJIu5hlnC2/c7D/H0p2t5d/FmEmKiuKVzQ27r1ojEhBivo4lEHA2zSNDl7DjIEx+vYdq3W0lMiOGO7o25qVMDKsTq2bpIoIVkz1zFPLIs37KPxz9aw6erdpBUKY7bujbkug4NqBSnoi4SKCrmUmYWfb+bpz5Zy5y1O0lMiOHmTqlcc0E9UqomeB1NJOypmEuZ+2bTXp79NIePV24HoEOj6vQ7L4Vz61alca2KxEX/3/nqx04UsPvQMY4cz+fYiQKO5xdQq0ocNSvF6aUlEUK0mOsBaPmwcddh3v1mM+8s3sx3Ow8BEOUz6lVLICbKR4Fz5Bc4dh86xv6jJ055j8SEGJrUqkSPpjUZ0K4+NStr5oyUL3oAKiHDOce6vIOs3HqANdsPsH7nIZxzGIbPZ1SrEENSpTiqV4ylYlwUsVFRRPmMrfuOkLPjIKu2HSB7wx5io3xc1jqZu3o0Jq12Za+bJVKm9Dq/eM7MaFKrMk1qnXkBXpd3kJfnb+CtwmUG/tT3HG7s2EBDMCKUbg9QkTLVuGYlHryiBZ//vgddmiTxwNTlDHlxEbsO/uB1NBHPqZhL2EmqFMcLN6XzYEZz5uTspP/YeWzee8TrWCKeUjGXsGRm3Ny5IZOHdmDPoWMMyJxP7p7DXscS8YxnxVwbOksgtK1fjZeHtGfv4eMMyFyggi4RSRs6S7mxZNNern/hS2pUjOWduzpriV6JSKebzaJhFokI59arysRbLmDL3qPc8Uo2x04UeB1JpEypmEvEOL9BdR69ujVffrebP7+7FK/+r1PEC5pnLhHlV21SWJ93kKc/zaFxzUrc3r2x15FEyoSKuUSc+3o1Zd3OQzzy4Sqa1q5Mz7NreR1JJOgCPsxiZl3NbJyZ/dvM5gX6/iK/xOcz/nX1uTRPrsLwyYtZl3fQ60giQedXMTezCWa2w8yWnXS+j5mtNrMcMxsJ4Jyb45y7A3gfeDHwkUV+WUJsFJk3phMb5eO2lxax/+hxryOJBJW/PfOJQJ+iJ8wsChgDXAo0BwaaWfMiHxkETApARpEzklI1gbHXtWXjrsMMn7SYE/ma4SKRy69i7pybDew+6XQ7IMc5t945dwyYDPQDMLP6wD7n3P5AhhUpqfaNavC3fi35fHUef353mWa4SMQqzQPQFGBTkeNcoH3hfw8G/nO6i/Py8khP//9z34cOHfrTG04iATWofX227D3Cs5/lULtKPL+5uKnXkUT8lpmZSWZm5k+HScV9rjTF/FTrjjoA59wDv3RxzZo10RugUlZ+e0lTtu8/ylOfrKV2lXgGta/vdSQRvxTt6JrZzuI+V5rZLLlAvSLHdYEt/l6stVmkLJkZ/3tlK3o2q8mf3l3KS/O/9zqSiN8CujaLmaUC7zvnWhYeRwNrgIuAzcBCYJBzbrk/99PaLOKFo8fzGfbaYj5euZ17LmzCiIubanMLCRulXpvFzCYB84FmZpZrZoOdcyeAYcAMYCXwhr+FHNQzF2/Ex0Qx7vq2XJtej2c+zWHUlKVax0VCnlZNFCmGc47HP1rDs5/l0DKlCk8NaEPjmpW8jiVyWiG5aqJ65uIlM+N3vZuRecP55O45wuVPf8HrCzdq6qKEJPXMRfywbd9RRrzxDfPW7aJb05o89KuW1KtewetYIv+HeuYip3FWYjyvDG7P3/q14OsNe7j4iVlkzl6nN0YlZKhnLlJCW/cd4f53l/Pxyu20rpvIo1e35uyzqngdSwQI0Z65SChKTkxg/I3n8+ygNmzec4SMZ77giZlrNONFQp6GWUROYmZc3roOM0d0p2+rZJ76ZC2/fn4+G3dps2jxhoZZRALgg6Vb+ePb34KDh65sxRXn1vE6kpRTGmYRKYW+rZKZfm9Xmp5VmeGTFvOX95ZxXA9HJcRomEXED3WrVeD1oR0Y2q0RL83fwI0vfMWeQ8e8jiXlhIZZRIJgyte5jJyylNpV4phw0wWk1a7sdSQpJzTMIhJAV7aty+tDO3D0eAFXPTePr747ed8WkbKnYi5yBtrUr8aUOzuRVDmO61/4kulLt3odSco5FXORM1SvegXevqMTLetU4a7XvuZlrZEuHtIDUJFSqFYxlleHdOCis2tz/3vLeWLmGi3WJQGnB6AiZeREfgGjpizlzexcru9Qn79e0ZIonza9kMA63QPQ0uwBKiKFoqN8PHp1a2pUimPcrHXsOniMJ649j/iYKK+jSTmhYi4SIGbGyEvPJqlSLP+YtpJdh75i/I3pJCbEeB1NygE9ABUJsCFdG/HUgPNYvHEPvx43j637jngdScoBFXORIOh3Xgov3tKOLXuP0n/MPJZv2ed1JIlwms0iEiSdmiTx5h0d8Rn8etx8Plm53etIEqY0m0UkBOzYf5TBLy5i+ZZ9/Omy5tzaORUzzXSRktPr/CIeqlUlntdv70Cvc2rz9/dX8Pu3vuXo8XyvY0mEUTEXKQMVYqMZd/35DL8ojbeyc7k2cwHb9h31OpZEEBVzkTLi8xkjLm7KuOvbsnb7AS5/5gvmr9vldSyJECrmImWsT8tk3r27M1USornu3wt47vN1FBRoCQApHRVzEQ80rV2ZqcO6cGmrZP754SqGvryIvYe12YWcuYAXczPzmdlDZvaMmd0U6PuLRIpKcdE8O7AND2Q0Z9aaPC57+guyN+zxOpaEKb+KuZlNMLMdZrbspPN9zGy1meWY2cjC0/2AFOA4kBvYuCKRxcy4pXND3rqjEz4fXPP8fMbN0rCLlJy/PfOJQJ+iJ8wsChgDXAo0BwaaWXOgGTDfOTcCuDNwUUUi17n1qjJteFd6t6jNI9NXccvEhew8+IPXsSSM+FXMnXOzgZP3xmoH5Djn1jvnjgGT+bFXngv89P+Kmkwr4qcq8TGMGdSWv/+qJfPX76LvU3OYl7PT61gSJkozZp4CbCpynFt4bgrQ28yeAWYXd3FeXh7p6ek//2RmZpYiikhkMDNu6NCAd+/qTKX4aK574UtGf7SaE/kFXkcTj2RmZv5cJ4Gk4j5XmiVwT/U+snPOHQYG/9LFsbGxtG3bloyMDDIyMkoRQyTyNK9Thffv6cJf3lvO05/msOC73Tw9oA1nJcZ7HU3K2NChQ0lOTiYrK4vs7OzjxX3O77VZzCwVeN8517LwuCPwoHOud+HxKADn3MP+3E9rs4j45+3sXO5/bxnxMVE8fs259GxWy+tI4pFgrc2yEEgzs4ZmFgsMAKb6e7FWTRTxz1Xn12XqsC7UqhzHLf9ZyD8/XKVhl3ImYKsmmtkkoAc/jtdsBx5wzr1gZn2BJ4EoYIJz7iF/w6lnLlIyR4/n89es5Uz6ahMXpFbjmYFtNexSzpyuZ+7ZErhpaWmuZ8+eGjMXKaH3vtnMqClLiY+J4qkB59E1rabXkSTIsrKyyMrKYvz48TnOubRTfUbrmYuEoZwdB7nr1WzW7jjI8AvTGH5RGlE+rZEe6UJyPXONmYucuSa1KvHu3Z3p3yaFpz5Zy83/+YpdeskoYmmnIZEI55xj8sJNPDB1OTUqxvLsoLac36Ca17EkSEKyZy4ipWdmDGxXnyl3diI6yrj2+fn8Z+53eNVJE+9omEUkArRMSeT9YV3p0awmf81awbDXFnPgaLHvl0iY0TCLSDlTUODInLOex2aspn71Coy9ri3nJFfxOpYEiIZZRMoJn8+4o3tjXhvSnkM/nOBXY+by2pcbNexSDmiYRSQCtW9Ugw/u7Uq7htX5n3eWMnzyNxp2CWMaZhEp5woKHM/NWsfjH62mbrUKPD2wDefVq+p1LDlDGmYRKad8PuPunk144/aO5Bc4rn5uHmM+yyFfOxlFHBVzkXIgPbU6H9zbld4tz+KxGasZNH4BuXsOex1LAkjFXKScSEyI4dmBbXj06tYs27yPS5+cwzuLc/VwNELoAahIOWJmXJNej+n3dqPZWZX5zetLuPu1r9l96JjX0eQ09ABURIqVX+B4fvY6npi5hsSEWB65shW9mtf2Opachh6Aisj/EeUz7urRhKnDulCzchxDXlrE795cwn5NYQxLKuYi5dw5yVV47+7ODOvZhHcWb6b3E7OZtSbP61hSQirmIkJstI/f9W7GlDs7UTEumpsmfMWoKd/qRaMwomIuIj87t15V3r+nC7d3a8TrCzfR58k5zFmrXno40GwWEfkv8TFRjOp7Dm/d2Ym4GB83vKBeutc0m0VESuXo8XxGz1zDv+esJzkxgUeuaqU9Rz2k2SwickbiY6L4n8JeenxhL33k299qxksIUjEXkV/Utn41pg3vyh3dG/PGok1cMno2n63a4XUsKULFXET8Eh8TxchLz+aduzpTJSGaWyYuZMTr37BHb4+GBBVzESmRc+tVJeueLtxzYROmLtnCxU/MYtq3W7XGi8dUzEWkxOKio/jtJc2YOqwLyYkJ3P3a19z+cjbb9x/1Olq5FfBibmY9zGyOmY0zsx6Bvr+IhI7mdarwzl2d+GOfs5m1Jo9eo2cx6auNFGi99DLnVzE3swlmtsPMlp10vo+ZrTazHDMbWXjaAQeBeCA3sHFFJNRER/m4s0djPryvGy3qVGHUlKUM+vcCvtt5yOto5Yq/PfOJQJ+iJ8wsChgDXAo0BwaaWXNgjnPuUuCPwF8DF1VEQlnDpIpMuq0Dj1zZiuVb9tP7ydmM+SyH4/kFXkcrF/wq5s652cDuk063A3Kcc+udc8eAyUA/59xPf3N7gLiAJRWRkGdmDGhXn09GdOeis2vx2IzVZDzzBUs27fU6WsQrzZh5CrCpyHEukGJmV5rZ88DLwLPFXZyXl0d6evrPP5mZmaWIIiKhpFaVeJ67/nyev+F89hw+Rv+xc/lb1goO/XDC62hhJzMz8+c6CSQV97noUvwOO8U555ybAkz5pYtjY2Np27YtGRkZZGRklCKGiISq3i3OomPjGvxz+iomzP2OGcu38Y/+LenZrJbX0cLG0KFDSU5OJisri+zs7GJfvfV7bRYzSwXed861LDzuCDzonOtdeDwKwDn3sD/309osIuXLwu93M/Ltb1mXd4h+59Xh/subk1RJI7ElEay1WRYCaWbW0MxigQHAVH8v1qqJIuXLBanV+eDertzXK40Plm6l1+hZvJWtDaX9EbBVE81sEtCDH8drtgMPOOdeMLO+wJNAFDDBOfeQv+HUMxcpv3J2HGDk20tZtGEPXZok8VD/ljSoUdHrWCHvdD1zz5bATUtLcz179tSYuUg5VVDgePWrjfxz+ipOFBTwm15NGdylIdFRejH9ZFlZWWRlZTF+/Pgc51zaqT6j9cxFxFNb9x3h/neX8/HK7bRMqcIjV7amZUqxownlWkiuZ64xcxEBSE5MYPyN5zP2urZs2/cD/cbM5ZHpqzh6PN/raCFDOw2JSFjZd/g4D32wgjcW5ZJaowIPX9majo1reB0rZIRkz1xE5GSJFWJ49OpzeXVIewocDBy/gFFTlmpnIz9omEVEQk7nJknMuK8bt3VtyOsLN3LJ6Nl8umq717E8o2EWEQl732zayx/eWsKa7Qfp3yaFBzKaU7VCrNexPKFhFhEJW+fVq8r793Rl+EVpZC3ZQq/Rs5mxfJvXsUKOirmIhLzYaB8jLm7Ke8M6U7NyHLe/nM29kxez97D2H/2JxsxFJGy0qJPI1GGdua9XGtO+3UrvJ2fz2eodXscKOo2Zi0jEWrZ5HyPe+IY12w8yqH19/nzZOVSILc1CsKFPY+YiEnFapiQydVgXbu/WiElfbeTyp79gae4+r2N5RsVcRMJWfEwUo/qew6tD2nP4WD79x85l3Kx15XJDaY2Zi0jY69Q4iQ/v68olLWrzyPRV3PriQnYfipyHoxozF5FyxTnHK19u5O9ZK6heMZZnBrXhgtTqXscKGI2Zi0i5YGbc0KEBU+7qRFyMj4GZC3hp/vflYgMMFXMRiTgtUxLJuqcL3ZvW5C/vLWfk20v54URkr8KoYi4iEalKfAzjb0xnWM8mvL5oEwMzF7Dz4A9exwoaFXMRiVg+n/G73s0Ye11bVmzdT/+xc8nZccDrWEGh2SwiEvH6tkpm8tCOHDmWz5Vj5zFv3U6vI5WIZrOIiBSxafdhbp24kO93HeLRq1vTv01dryOViGaziIgA9apX4K07O5HeoDq/eX0JYz7LiZiZLirmIlKuJCbEMPHWC7ji3Do8NmM1f353GSfyC7yOVWqRvSqNiMgpxEVH8eS155FcNZ7nZ63n+12HeHZgW6pVDN9NL9QzF5FyyeczRl16Do9e3ZqF3+3hijFfsGrbfq9jnTEVcxEp165Jr8fk2zvww/ECrhw7jzcWbgrLcfSgFHMzq2hm2WZ2eTDuLyISSG3rVyPrni60rpvIH97+ltteWsSOA0e9jlUifhVzM5tgZjvMbNlJ5/uY2WozyzGzkUX+6I/AG4EMKiISTLWrxPPakA7cf3lzZq/dSe8nZvPeN5vDppfub898ItCn6AkziwLGAJcCzYGBZtbczHoBK4DtAcwpIhJ0Pp8xuEtDPhjehfo1KnLv5G+4ZeJCcvcc9jraL/KrmDvnZgO7TzrdDshxzq13zh0DJgP9gJ5AB2AQcJuZaVxeRMJKk1qVmXJnJ/5yeXO++m43F4+ezb/nrCc/hDe9KE2hTQE2FTnOBVKcc39yzt0HvAaMd86dcgJnXl4e6enpP/9kZmaWIoqISGBF+YxbuzRk5ojudGhUnX9MW0n/sXNZvqVst6bLzMz8uU4CScV9rjTzzO0U537+Z8s5N/F0F8fGxtK2bVsyMjLIyMgoRQwRkeBJqZrAhJsvYNrSrTw4dQVXPDuXIV0acl+vpiTERgX99w8dOpTk5GSysrLIzs4+Xtzn/F6bxcxSgfedcy0LjzsCDzrnehcejwJwzj3sz/20NouIhJt9h4/z8PSVTF64iXrVE3joV63o1rRmmf3+YK3NshBIM7OGZhYLDACm+nuxVk0UkXCTWCGGR65qzeShHYjx+bhxwlfcN3kxu4K8TnrAVk00s0lAD34cr9kOPOCce8HM+gJPAlHABOfcQ/6GU89cRMLZ0eP5jP0sh+dmraNiXDR/vqw5V7VNwexUI9CBcbqeuWdL4KalpbmePXtqzFxEwtra7QcYNWUpizbsoVPjGvxv/1akJlUM6O/IysoiKyuL8ePH5zjn0k71Ga1nLiJSSgUFjte+2sg/p6/iWH4Bwy9KY2i3RsREBXZmdkiuZ64xcxGJFD6fcX2HBnz82+5ceHYtHpuxmoxnvmDxxj0Bub92GhIR8cBHy7fxl/eWs/3AUW7qmMrvejejUlzpVxwPyZ65iEikuqTFWcwc0Y2bOqby4vzvuXj0LGauCO4KJxpmEREJgsrxMTx4RQum3NmJxIQYbntpEXe+ks2O/SVfjVHDLCIiIeB4fgGZs9fz1CdrqZoQw+w/9CQ+puRvj55umEXbxomIBFlMlI+7ezahb6tkVm7df0aF/JeomIuIlJGGSRVpGOA56D/RmLmISIjTmLmISATR1EQRkQinYi4iEgE0Zi4iEuL8GTP3rJgnJiaSmZkZ8ismRvp2dpHcPrUtfEVy+86kbRkZGT9dV+yedRpm+QWR/KWCyG6f2ha+Irl9wWqbirmISATwbGqimeUBGzz55SWTBOz0OkQQRXL71LbwFcntK03bGjjnTrnpqGfFXEREAkfDLCIiEUDFXEQkAqiYi4hEABXzUjCzHmY2x8zGmVkPr/MEkpmdU9iut8zsTq/zBJqZNTKzF8zsLa+zBEKktaeocvBdDEgdKbfF3MwmmNkOM1t20vk+ZrbazHLMbOQv3MYBB4F4IDdYWUsqEG1zzq10zt0BXAOccmEfrwSofeudc4ODm7R0StLOcGhPUSVsW8h+F4tTwu9oYOqIc65c/gDdgLbAsiLnooB1QCMgFlgCNAdaAe+f9FML8BVeVxt41es2BbJthddcAcwDBnndpmC0r/C6t7xuTyDaGQ7tKU3bQvW7GIj2BaqOlNvNKZxzs80s9aTT7YAc59x6ADObDPRzzj0MXH6a2+0B4oIS9AwEqm3OuanAVDObBrwWxMglEuC/u5BVknYCK8o4XqmUtG2h+l0sTgm/oz/93ZWqjpTbYl6MFGBTkeNcoH1xHzazK4HeQFXg2aAmK72Stq0HcCU/frk+CGawAClp+2oADwFtzGxUYdEPB6dsZxi3p6ji2taD8PouFqe49gWkjqiY/zc7xbli36pyzk0BpgQvTkCVtG2fA58HK0wQlLR9u4A7ghcnaE7ZzjBuT1HFte1zwuu7WJzi2heQOlJuH4AWIxeoV+S4LrDFoyyBFsltg8hv308iuZ2R3DYIcvtUzP/bQiDNzBqaWSwwAJjqcaZAieS2QeS37yeR3M5IbhsEu31eP/X18GnzJGArcJwf/8UcXHi+L7CGH586/8nrnGpb+WtfeWhnJLfNq/ZpoS0RkQigYRYRkQigYi4iEgFUzEVEIoCKuYhIBFAxFxGJACrmIiIRQMVcRCQCqJiLiEQAFXMRkQjw/wCs3qJtYVv7UQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Sigma and gamma will later be made more rigerous\n",
    "\n",
    "omegaCDM = 0.2589\n",
    "sigmasquared = 0.003\n",
    "\n",
    "k = 3.3\n",
    "gamma = 0.36\n",
    "Meq = 2.8e17\n",
    "\n",
    "def mu(M, MH) :\n",
    "    return M / (k*MH)\n",
    "\n",
    "# compute integral for a specific M value\n",
    "def fraction(M):\n",
    "    return (1/omegaCDM)*integ.simps( 2 / np.sqrt(2*np.pi*sigmasquared) * \n",
    "                np.exp(-( mu(M, MH_grid)**(1/gamma) + deltac(MH_grid))**2 / (2*sigmasquared) ) *\n",
    "                M/(gamma*MH_grid) * mu(M, MH_grid)**(1/gamma) * np.sqrt(Meq/MH_grid),\n",
    "            np.log(MH_grid)\n",
    "        )\n",
    "\n",
    "Ms = np.logspace(-6,5,100)\n",
    "points = np.array([fraction(m) for m in Ms])\n",
    "plt.loglog(Ms, points)\n",
    "# OmegaCDM = integ.simps(points, Ms)\n",
    "# print(\"$Omega_{CDM}$ is\", OmegaCDM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}