{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "importing packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import path as Path\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from astropy.cosmology import FlatwCDM\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "from astropy.cosmology import LambdaCDM\n",
    "from astropy.cosmology import wCDM\n",
    "from astropy.cosmology import Flatw0waCDM\n",
    "import scipy\n",
    "from escape_functions_noastropy import *\n",
    "from multiprocessing import Process\n",
    "from multiprocessing import Queue\n",
    "from multiprocessing import Pool\n",
    "import emcee\n",
    "import corner\n",
    "\n",
    "# from multiprocessing import cpu_count\n",
    "# ncpu = cpu_count()\n",
    "# print(\"{0} CPUs\".format(ncpu))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "initializing cosmology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cosmology(cosmology):\n",
    "    case = cosmology.name\n",
    "    if case == 'Flatw0waCDM':\n",
    "        return [cosmology.Om0, cosmology.w0, cosmology.wa, cosmology.h]\n",
    "    \n",
    "    elif case == 'FlatwCDM':\n",
    "        return [cosmology.Om0, cosmology.w0, cosmology.h]\n",
    "\n",
    "    elif case == 'wCDM':\n",
    "        return [cosmology.Om0, cosmology.Ode0, cosmology.w0,cosmology.h]\n",
    "        \n",
    "    elif case == 'LambdaCDM':\n",
    "        return [cosmology.Om0, cosmology.Ode0, cosmology.h]\n",
    "\n",
    "    elif case == 'FlatLambdaCDM':\n",
    "        return [cosmology.Om0, cosmology.h]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "creating fake data: true values, then adding error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEGCAYAAAB2EqL0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5zWc/7/8cers1FC5VRNUys2rCWjnNY6S5KzZYdCX21kZXexksP3i2h/whKxsW1hSBFinXJmNzQlbYQGHUbaklMMUb1+f7w/o6tpZvrMzHVdnzk877fbdbuu6/35fGZe8+lqXvP5vA8vc3dEREQ2pUnSAYiISP2ghCEiIrEoYYiISCxKGCIiEosShoiIxNIs6QAypX379p6Xl5d0GCIi9cqsWbM+c/cOFW1rsAkjLy+PoqKipMMQEalXzGxRZdt0S0pERGJRwhARkViUMEREJJYG24chIlJdP/74IyUlJXz//fdJh5JxrVq1olOnTjRv3jz2MUoYIiKRkpIS2rRpQ15eHmaWdDgZ4+6sXLmSkpISunbtGvs43ZISEYl8//33tGvXrkEnCwAzo127dtW+klLCEBFJ0dCTRZma/JxKGOUUFkJeHjRpEp4LC5OOSESkbshowjCzzmb2opnNN7N3zGxY1H6Dmb1nZnPN7BEz2zLlmOFmVmxm75vZkSntfaK2YjO7NBPxFhbC4HPWsWgRuMOiRTB4sJKGiGTHl19+ydixY6t93IQJE1i6dOlP7/Py8vjss8/SGRqQ+SuMNcCf3L0HsA8w1Mx2AaYDu7n77sAHwHCAaNupwK5AH2CsmTU1s6bA7cBRwC7AadG+aTViBJR+t+EpKS0N7SIimVZZwli7dm2Vx5VPGJmS0VFS7v4p8Gn0epWZzQc6uvuzKbu9DpwUvT4WmOTuq4GPzawY6BVtK3b3jwDMbFK077vpjHfx4uq1i4ik06WXXsqHH37IHnvsQfPmzWndujXbb789c+bM4cknn6Rfv37MmzcPgNGjR/PNN9+w2267UVRUREFBAZttthkzZswAYMyYMTz++OP8+OOPTJkyhZ///Oe1ji9rw2rNLA/YE3ij3KazgQej1x0JCaRMSdQGsKRce+90x5ibG25DVdQuIo3MhRfCnDnp/Zp77AF//Wulm0eNGsW8efOYM2cOL730EkcffTTz5s2ja9euLFy4sMJjTjrpJG677TZGjx5Nfn7+T+3t27dn9uzZjB07ltGjR3P33XfXOvysdHqbWWvgYeBCd/86pX0E4bZVWS9BRd32XkV7+e8z2MyKzKxoxYoV1Y5z5EjIydmwLYdvGXl6Wi9kRERi6dWrV7XmSaQ64YQTANhrr70qTTbVlfErDDNrTkgWhe4+NaV9INAPONTdy375lwCdUw7vBJTdmKus/SfuPg4YB5Cfn79RQtmUgoLwPGJEuA2V22kdI9dcRcFd98KQWdCpU3W/pIjUV1VcCWTL5ptv/tPrZs2asW7dup/eb2oORcuWLQFo2rQpa9asSUs8mR4lZcDfgfnuflNKex/gz0B/dy9NOWQacKqZtTSzrkB34E1gJtDdzLqaWQtCx/i0TMRcUAALF8K6dbBwcRMKnh8Uer5POglWr87EtxQRAaBNmzasWrWqwm3bbrsty5cvZ+XKlaxevZonnngi1nHplOkrjP2BM4D/mFnZzcDLgFuBlsD0aPLI6+4+xN3fMbPJhM7sNcBQd18LYGbnA88ATYHx7v5OhmMPevSAiRPhxBNh2DC4886sfFsRaXzatWvH/vvvz2677cZmm23Gtttu+9O25s2bc+WVV9K7d2+6du26QSf2mWeeyZAhQzbo9M4EW383qGHJz8/3tBZQGj4cRo2Cu++GQYPS93VFpM6YP38+PXr0SDqMrKno5zWzWe6eX9H+mukd17XXwmGHwdChMHNm0tGIiGSdEkZcTZvCAw/AdtuF21M1GIUlIlKfKWFUR/v2MHVqSBanngppGnkgInVHQ71NX15Nfk4ljOrq2TN0fL/wAlx2WdLRiEgatWrVipUrVzb4pFFWD6NVq1bVOk4FlGpi4EB480244QbYe284+eSkIxKRNOjUqRMlJSXUZOJvfVNWca86lDBq6uab4a234KyzYJddYNddk45IRGqpefPmNZ5Z3RjollRNtWgBDz0ErVvD8cfDV18lHZGISEYpYdTGDjvAlCnw8ccwYECYHi4i0kApYdTWr34FN90E06bBddclHY2ISMYoYaTD+efD6afDlVfCU08lHY2ISEYoYaSDGfztb7D77vDb38KHHyYdkYhI2ilhpEtOTpjUZxZmgpeWbvoYEZF6RAkjnbp1g/vvh7lzYfBgaOCTf0SkcVHCSLc+feCaa6CwEG67LeloRETSRgkjE4YPh/794Y9/hFdfTToaEZG0UMLIhCZN4J57oGvXsGzI0o2qyYqI1DtKGJnSti088gh8800o7/rDD0lHJCJSK0oYmbTrrvCPf8CMGfCHPyQdjYhIrShhZNrJJ8PFF8PYsTBhQtLRiIjUmBJGNlx3HRxyCAwZArNnJx2NiEiNKGFkQ7NmMGkSbLMNnHACfPZZ0hGJiFSbEka2dOgQZoIvWwannQZr1yYdkYhItShhZFN+fujLeO45uPzypKMREakWJYxsO/ts+N3vYNQoePjhpKMREYlNCSMJt9wCvXvDmWfC/PlJRyMiEosSRhJatgzlXXNyQnnXr79OOiIRkU1SwkhKp04weTIUF8PAgSrvKiJ1nhJGkn79axg9Gh59FP7yl6SjERGpkhJG0oYNC8NsL78cnn026WhERCqlhJE0M7jrrrDu1GmnwcKFSUckIlKhjCYMM+tsZi+a2Xwze8fMhkXtW5vZdDNbED1vFbWbmd1qZsVmNtfMeqZ8rYHR/gvMbGAm4866zTcPK9uuWwcnnEDhP34gLy+skp6XF2oxiYgkLdNXGGuAP7l7D2AfYKiZ7QJcCjzv7t2B56P3AEcB3aPHYOAOCAkGuAroDfQCripLMg3Gz34GhYUUvtWDwec4ixaFCq+LFoVqr0oaIpK0jCYMd//U3WdHr1cB84GOwLHAxGi3icBx0etjgXs8eB3Y0sy2B44Eprv75+7+BTAd6JPJ2BPRty8j2t5O6dqWGzSXlsKIEQnFJCISyVofhpnlAXsCbwDbuvunEJIKsE20W0dgScphJVFbZe3lv8dgMysys6IVK1ak+0fIisVft624fXGWAxERKScrCcPMWgMPAxe6e1Wz1KyCNq+ifcMG93Hunu/u+R06dKhZsAnLza3oR4Xc3CwHIiJSTsYThpk1JySLQnefGjX/N7rVRPS8PGovATqnHN4JWFpFe4MzcmSYAJ4qp9U6Ro5MJh4RkTKZHiVlwN+B+e5+U8qmaUDZSKeBwGMp7QOi0VL7AF9Ft6yeAY4ws62izu4jorYGp6AAxo2DLl3AzOnStIRxTYdQsHNR0qGJSCOX6SuM/YEzgEPMbE706AuMAg43swXA4dF7gCeBj4Bi4C7gPAB3/xy4BpgZPa6O2hqkgoIwHWPdOmPhBz9Q0GF6qNj36qtJhyYijZi5b9QV0CDk5+d7UVED+av8k0/gsMPCGNtHH4Ujjkg6IhFpoMxslrvnV7RNM73rg44d4eWXYeed4ZhjwiQ/EZEsU8KoL7bZBl54AXr2hJNPhvvuSzoiEWlklDDqk622gunTwyq3AwbAnXcmHZGINCJKGPVN69bwz3/C0UfDuefCDTckHZGINBJKGPVRq1YwdSr85jdwySVw1VVh4SkRkQxqlnQAUkPNm4cVCTffHK6+GlatghtvDMuli4hkgBJGfda0aail0aYN3HxzSBp33hnaRUTSTAmjvmvSJCSLNm3g2mvh229h4sRwBSIikkZKGA2BGVxzTUgaf/5zSBoPPhj6OkRE0kSd3g3JJZfA7bfDtGnQrx98803SEYlIA6KE0dCcd164JfXii3DkkfDll0lHJCINhBJGQzRgAEyeDDNnhkUL62kxKRGpW5QwGqoTTwy3pubPDzPDP/kk6YhEpJ5TwmjI+vSBZ56BkhI48ED4+OOkIxKReqzKhGFmW1f1yFaQUgsHHgjPPw9ffAG/+hW8917SEYlIPbWpK4xZQFH0vAL4AFgQvZ6V2dAkbfbeOyyPvmZNSCBz5iQdkYjUQ1UmDHfv6u7dCOVQj3H39u7eDugHTK3qWKljfvGLULGvVSs4+GCYMSPpiESknonbh7G3uz9Z9sbdnwJ+nZmQJGO6dw9Jo317OPzwUF9DRCSmuAnjMzO73MzyzKyLmY0AVmYyMMmQLl3glVega1fo2xeeeCLpiESknoibME4DOgCPRI8OUZvUR9tvDy+9FG5THX98WEZERGQTYq0l5e6fA8PMrLW7a72JhqBduzB6ql8/OO20sP7U2WcnHZWI1GGxrjDMbD8zexd4N3r/SzMbm9HIJPO22AKefhqOOAIGDYJbb006IhGpw+LekroZOJKo38Ld3wYOzFRQkkU5OfDYY3DCCTBsGIwcqep9IlKh2DO93X1Juaa1aY5FktKyZejHOOMMuPxyGD5cSUNENhK3HsYSM9sPcDNrAVwAzM9cWJJ1zZrBhAnQujX85S+het+YMaFAk4gI8RPGEOAWoCNQAjwLnJepoCQhTZqEehqtW8MNN4R6Gn//e0gmItLoxf1NsLO7F6Q2mNn+wL/SH5IkyixcYWyxBVxxRUga998fbluJSKMW937DmJht0hCYhb6Mm2+GqVPhmGNgpeZpijR2VV5hmNm+wH5ABzP7Y8qmLYCmmQxM6oALL4Qtt4Tf/Q723DMUZdpnn6SjEpGEbOoKowXQmpBY2qQ8vgZO2tQXN7PxZrbczOaltO1hZq+b2RwzKzKzXlG7mdmtZlZsZnPNrGfKMQPNbEH0GFj9H1NqorAQ8v73TJr8+D15S/9N4X63w1//qhFUIo2UeYz//GbWxd0XVfuLmx0IfAPc4+67RW3PAje7+1Nm1he4xN0Pil7/HugL9AZucffeUd2NIiAfcMKy6nu5+xdVfe/8/HwvKiqqbsgSKSyEwYOhtHR9W07T7xm39mwKjv8exo8PVx8i0qCY2Sx3z69oW9w+jLvN7KffDma2lZk9s6mD3P0V4PPyzYRbWgBtgaXR62MJicXd/XVgSzPbnjBhcLq7fx4lielAn5hxSw2NGLFhsgAoXduKEVvdCY8/Dj17wiyVRBFpTOImjPbu/mXZm+gX9zY1/J4XAjeY2RJgNDA8au8IpE4OLInaKmvfiJkNjm5zFa1YsaKG4QnA4sWVtH+5RVgifc0a2G8/uOMO3aISaSTiJox1ZpZb9sbMuhCuFGriXOAP7t4Z+APw97IvW8G+XkX7xo3u49w9393zO3ToUMPwBCA3t4r2ffaBt96Cww6D884LixeuWpXV+EQk++ImjBHAa2Z2r5ndC7zC+iuD6hrI+mp9U4Be0esSoHPKfp0It6sqa5cMGjkyLDOVKicntANhtdvHH4frr4eHHoL8fJg7N+txikj2xEoY7v400BN4EJhM6HTeZB9GJZayvlrfIYQa4QDTgAHRaKl9gK/c/VNCedgjon6TrYAjojbJoIICGDcu1FsyC8/jxoX2nzRpApdeGir3rVoFvXuHznDdohJpkGLN9DYzI3Q0d3P3q80s18x6ufubmzjuAeAgoL2ZlQBXAecAt5hZM+B7YHC0+5OEEVLFQClwFoRaHGZ2DTAz2u/qqD6HZFhBQbkEUZkDD4Q5c8LOgwbByy/D2LGw+eYZj1FEsifusNo7gHXAIe7eI/pL/1l33zvTAdaUhtUmYO1auPZa+L//gx49wq2qHj2SjkpEqiEdw2p7u/tQwhVB2SipFmmKTxqKpk3hqqtg+nT47LPQr3HffUlHJSJpEjdh/GhmTYlGJ5lZB8IVh8jGDj00jKLKzw81NgYPhu++SzoqEamluAnjVuARYBszGwm8BlyXsaik/tthh1AzfPhwuOsu2HdfWLBg08eJSJ0Vd5RUIXAJcD3wKXCcu0/JZGDSADRrBtddB08+CUuWwF57wRR9bETqqyoThpltET1vDSwHHgDuB/4bDXPVirWyaUcdFUZR7bYbnHIK/P73sHp10lGJSDVt6grj/uh5FmEBwFkpj9nAMjPTrSnZtM6dw3DbP/0JbrsNDjgAPv446ahEpBqqTBju3i967uru3aLnnx7AdsDx2QhUGoDmzWH0aHj0USguDgsYPvZY0lGJSExxO70xs/5mNjp6lCWSte6ugfZSPcceC7Nnw447wnHHhauOH39MOioR2YRYCcPMRgHDgHejxzAzuz6TgUkD17UrvPYanH8+3HQT/PrXoWNcROqsuFcYfYHD3X28u48nLBNydObCkkahZUsYMyaUfp03L5SBfeqppKMSkUrEviUFpJZXa5vuQKQRO/nkUIypUyfo2xcuuyzU2xCROiVuwrgeeMvMJpjZRMIoKY2OkvTp3h1mzIBzzglLph96KCzVKvYidckmE0a0Uu1rwD6EOhZTgX3dfVKGY5PGZrPNwhrq994LRUXhFtXzzycdlYhENpkwPCxn+6i7f+ru09z9MXdfloXYpLE6/XSYORPat4fDDw+r365dm3RUIo1e3FtSr5tZnV3KXBqgXXaBN9+E00+n8H/fJ6/1ZzRp4uTlQWFh0sGJNE6xCigBBwNDzGwh8C2hzra7++6ZCkyEzTen8MiJDH5wDaXfNwdg0SIYPDiUeo9V3ElE0iZuAaUuFbW7+6K0R5QmKqDUMOTlhSRRXpcdfmThJ82zHo9IQ1frAkpRYmgHHAv0B9rV5WQhDcfixZW0L20KI0fCDz9kNyCRRizuTO8rgYmEpNEe+IeZXZ7JwEQAcnMrac9ZCZdfHtaj+te/shuUSCMVt9P7NGBvd7/K3a8iDLHVHWTJuJEjISdnw7acHBg5rgM8/jisWhVWvh0yBL78MpkgRRqJuAljIdAq5X1L4MO0RyNSTkFBmJrRpQuYhedx40I7/frBO+/AH/4Qqvr16BGWGYnRLyci1Re30/tRYG9gOqGu9+GEyXzLAdz9ggzGWCPq9G5kZs0KtcNnz4ajj4bbbw/ZRUSqpapO77jDah+JHmVeqm1QImm1117wxhthMcMrrgjzOK65Bi64IJSKFZFai3WFsckvYvawu5+YhnjSRlcYjdjixTB0KDzxROgUHzcuJBQR2aRaD6uNoVuavo5I7eXmwrRpMGUKfPop9OoV+jm++SbpyETqtXQlDPUySt1iBiedBPPnw+9+B7fcEm5TPf540pGJ1FvpShgidVPbtjB2bJir0bYt9O8fEomWTheptnQlDEvT1xHJjH33DSOorrsO/vnPMAR37FhYty7pyETqjbgzvTc3syYp75uYWep0qj+nPTKRdGveHIYPD+Vge/UKHeP77w//+U/SkYnUC3GvMJ4HUhNEDvBc2Rt3f7aig8xsvJktN7N55dp/b2bvm9k7Zvb/UtqHm1lxtO3IlPY+UVuxmV0aM2aRiv3sZ/Dss6FQU3FxGEk1fDh8913SkYnUaXETRit3/2mISfQ6p4r9y0wA+qQ2mNnBhEUMd3f3XYHRUfsuwKnArtExY82sqZk1BW4HjgJ2AU6L9hWpObNQqOm99+CMM2DUKNhtN5g+PenIROqsuAnjWzPrWfbGzPYCNvnnmLu/AnxervlcYJS7r472WR61HwtMcvfV7v4xUAz0ih7F7v6Ru/8ATIr2Fam9du1g/Hh48cUwwe+II0IiWb5808eKNDJxE8aFwBQze9XMXgUeBM6v4ffcCfiVmb1hZi+nVPLrCCxJ2a8kaqusfSNmNtjMisysaMWKFTUMTxqlgw6Ct98Os8QnTw6d4uPHb7AuVWFhqM/RpAmq/CeNUtx6GDOBnxOuDs4Derj7rBp+z2bAVoQVby8GJpuZUfFIK6+ivaI4x7l7vrvnd+jQoYbhSaPVqhVcfTXMmRPmbAwaBAcfDO+/T2FhWKpq0aKQQ0LlPyUNaVzijpI6mdCPMY9wO+jB1FtU1VQCTPXgTWAdocZGCdA5Zb9OwNIq2kUyY5dd4OWXw5Iib78Nu+/OiKFfUlq64W6lpTBiRDIhiiQh7i2pK9x9lZkdABxJKKZ0Rw2/56PAIQBmthPQAvgMmAacamYtzawr0B14E5gJdDezrmbWgtAxPq2G31skniZN4JxzwkzxE05g8VdbVLhbZRUBRRqiuAljbfR8NHCHuz9G+EVfJTN7AJgB7GxmJWY2CBgPdIuG2k4CBkZXG+8Ak4F3gaeBoe6+1t3XEPpLngHmA5OjfUUyb7vt4IEHyN1mdYWbK6sIKNIQxa2H8QTwCXAYUDZC6k13/2Vmw6s5rVYr6VRYCIPPcUq/W9+lltPsB8aN+Z6CIRVffYjUR+lYrfYUwl/4fdz9S2BrQoe1SKNQUADj7rKo8p/TpfVnjFtzFgWXdIJrr4Vvv006RJGMiztKqpRQXe+AqGkNsCBTQYnURQUFsHAhrFtnLFzVnoJ3L4fDDgtDcX/2M7jjDvjxx6TDFMmYuKOkriKsFzU8amoO3JepoETqhR49YOpU+Pe/Yaed4LzzYNddQx0O1RWXBijuLanjgf7AtwDuvhRok6mgROqVffcNw3CfeAJatoRTToHeveGFF5KOTCSt4iaMHzz0jjuE1WszF5JIPWQGRx8dJv1NmADLlsGhh0KfPqFNpAGImzAmm9nfgC3N7BzCSrV3ZS4skXqqaVMYOBA++ABGj4aZM2HPPcP6VB9/nHR0IrUSt9N7NPAQ8DCwM3Clu4/JZGAi9VqrVvCnP8GHH4al06dOhZ13hmHDQOucST0Vu+Keu09394vd/SJ332ANaDObkf7QRBqALbcMVf4WLIAzz4Tbbgsjqq65Br75ZpOHi9Ql6SrR2ipNX0ekYerYMaxN9c47cPjhcOWVsOOOoUyshuJKPZGuhKExhCJx/Pzn8PDDMGNGGIo7dGhY7HDyZA3FlTovXQlDRKpjn33WD8Vt1Qp+85tQZ1xDcaUOiztx7w9m1qmqXdIUj0jjkToUd+LEUOXv0EPhyCPhrbeSjk5kI3GvMLYAnokq7g01s23LbT8jzXGJNB5Nm8KAAfD++3DjjVBUBD17hrVIPvoo6ehEfhJ3WO3/ufuuwFBgB+BlM3suZfu8DMUn0ni0agV//OP6obiPPBL6PMoNxVWpWElKdfswlgPLgJXANukPR0R+GopbXAxnnQW33w7dusHVV1M4/nuVipXExK2HcS7wG6ADYQLfg+7+boZjqxXVw5AG4733Qi3YqVPJa7KYRes6b7RLly5hJV2R2qqqHkazmF+jC3Chu2tRHJFsKxuK+/rrLN63Y4W7qFSsZEPcPoxLlSxEErbPPuTmVjwgMbfTuiwHI42R5mGI1CMjrzNycjZsy+FbRn5xbqj898UXyQQmjYIShkg9UlAQVhgJpWLD87grllDw609C5b/cXLjoIli6NOlQpQGK1eldH6nTWxqduXPhL3+BSZOgWbMwt+OSS6B796Qjk3qkqk5vXWGINBS77x7G1y5YAIMGwb33hg7z3/xGM8clLZQwRBqabt3CKriLFoUrjKefDjPH+/SBl17SIodSY0oYIg3VttvC9deHMbfXXx+uMg4+GPbbDx57DNZpZJVUjxKGSEPXti1cemmY2Td2LPz3v3DccfCLX8A996geh8SmhCHSWGy2GZx7bqg3Xli4vv74jjvCmDFQWpp0hFLHKWGINDbNmsFvfwtvvx3qcXTuDBdcEMboai6HVEEJQ6SxKqvH8dpr8Oqr0Lv3+rkcF1+suRyyESUMEYEDDghXG2+/Df37w003QdeuYSncBQuSjk7qiIwmDDMbb2bLzWyjehlmdpGZuZm1j96bmd1qZsVmNtfMeqbsO9DMFkSPgZmMWaRRKz+X4557NJdDfpLpK4wJQJ/yjWbWGTgcSF1j8yige/QYDNwR7bs1cBXQG+gFXGVmW2U0apHGrmwux8KFlc7lUCGnxiejCcPdXwE+r2DTzcAlQOoMomOBezx4HdjSzLYHjgSmu/vn7v4FMJ0KkpCIZMB221U4l6Ow+/8yeNAaFXJqZLLeh2Fm/YFP3P3tcps6AktS3pdEbZW1V/S1B5tZkZkVrUgpaSkitVRuLseIhedQunrDcjqlpaHOkzRcWU0YZpYDjACurGhzBW1eRfvGje7j3D3f3fM7dOhQ80BFpGLRXI7F6yor5KRlRxqybF9h/AzoCrxtZguBTsBsM9uOcOWQWnuyE7C0inYRSUilhZx8ERxxRFh6ZO3aLEclmZbVhOHu/3H3bdw9z93zCMmgp7svA6YBA6LRUvsAX7n7p8AzwBFmtlXU2X1E1CYiCRk5ko0LOW22jpEnvQXz54elR7p1C/0euj3cYGR6WO0DwAxgZzMrMbNBVez+JPARUAzcBZwH4O6fA9cAM6PH1VGbiCSkwkJOdzWhYMrx8PHHoQb5jjvCZZdBp05wxhnwxhtaKbeeUwElEcmc+fPD8NyJE2HVKthrLxg6FE49NfSHSJ2jAkoikowePcLChp98ArffDt99B2efHa46Lr4YPvoo6QilGpQwRCTz2rSB886DefPgxRfhkEPg5pvDbat+/eCpp1Sfox5QwhCR7DGDgw6CKVPCnI7LL4eiIujbF3baCW68ET5XF2VdpYQhIsno1AmuvjrMIn/gAdh+e7joIujYMaxjNXt20hFKOUoYIpKsFi1CJ/irr8KcOTBgAEyaFDrI99svrDeyenXSUQpKGCJSl/zyl/C3v4VO8ptvhs8+g9NPD0WeRoyAJUs22F0LIGaXEoaI1D1bbgkXXgjvvQfPPAP77gujRoWscPzx8PzzFBY6gwejBRCzSPMwRKR+WLgwXH3cdResXElesxIWrdl4TasuXcKuUjOahyEi9V9eXlhqpKQEJk5k8ZrtK9xt8eIKmyUNlDBEpH5p1QoGDCC3S8W/vnLbfgXLlmU5qMZBCUNE6qUKF0Bs8h0jvzw3DNk95hiYOhV++CGZABsgJQwRqZcqXADxns0omH9lmM8xaxaceCLssAMMGxaG7EqtqNNbRBqmNWtg+nSYMAEefTRcafzyl3DWWfDb34KKrFVInd4i0vg0awZHHQUPPgiffgq33QbNm4fhuh07wgknwOOPwyLBlfUAAAokSURBVI8/Jh1pvaGEISIN39Zbh2XVZ86EuXPhggvgX/+C/v3DpMCLLgoLI0qVlDBEpHH5xS9g9OgwPHfatLD8yC23hPa99w7LsGsBxAopYYhI49S8+fqRVEuXwl//Gm5PnX9+WAjxlFPCsuuqTf4TJQwRkQ4d1o+keustGDIEXnghLLuemwuXXhqWKWnklDBERFLtsUe4RbV0aahNvtde4RZWjx5hTatx4+CrrzY4pLEsgqhhtSIim7JsGdx3H/zjH/Duu2G2+QknwFlnUbjsEAb/rgmlpet3z8kJeaWgILmQa6qqYbVKGCIicbmHCoETJsD998OXX5LXdAmL1nbaaNf6ugii5mGIiKSD2fqRVJ9+Cg8+yOK1O1S4a0NcBFEJQ0SkJlq1glNOqXwRxJb/hb//Hb74IsuBZY4ShohILVS4CGLzHxi5xSj4n/+BbbcNw3cLC2HVqmSCTBMlDBGRWqhwEcR/tKBg2U2hv6NsuO7pp8M228DJJ4fRV999l3To1aZObxGRTFu3DmbMgEmTYPJkWL4cWreG446DU0+Fww+HFi2SjhJQp7eISLKaNIH994cxY+CTT+C550Ki+Oc/oV8/2G47OOcceP75Oj2zXAlDRCSbmjWDQw8NtcmXLQsr5vbtG64+DjssrKT7+9+HxRHXrUs62g0oYYiIJKVFi3CFcd998N//wpQpcMABcPfd4TkvDy6+OBSDqgPdBxlNGGY23syWm9m8lLYbzOw9M5trZo+Y2ZYp24abWbGZvW9mR6a094nais3s0kzGLCKSiJwcOOkkeOihkDzuvRd23z0sipifDzvtBFdcAe+8k1iImb7CmAD0Kdc2HdjN3XcHPgCGA5jZLsCpwK7RMWPNrKmZNQVuB44CdgFOi/YVEWmYttgijKp64omQPO66Kwy/uu462G23sBT7yJFQXJzVsDKaMNz9FeDzcm3Puvua6O3rQNmc+mOBSe6+2t0/BoqBXtGj2N0/cvcfgEnRviIiDd/WW4f5HM89FzrMx4yBtm3h8suhe/cw8/zGG2HJkowvgph0H8bZwFPR647AkpRtJVFbZe0bMbPBZlZkZkUrVqzIQLgiIgnabrtQr+O112DRIrjhhtC3cdFFFOb+mcEDvmPRotC0aBEMHpzepJFYwjCzEcAaoOzHsQp28yraN250H+fu+e6e30EF3kWkIcvNDaVli4rggw8Y0XYspes222CX0lIYMSJ937JZ+r5UfGY2EOgHHOrrZw6WAJ1TdusELI1eV9YuIiLdu7P464o3pXMRxKxfYZhZH+DPQH93T1lBnmnAqWbW0sy6At2BN4GZQHcz62pmLQgd49OyHbeISF2Wm1u99prI9LDaB4AZwM5mVmJmg4DbgDbAdDObY2Z3Arj7O8Bk4F3gaWCou6+NOsjPB54B5gOTo31FRCRS4SKIOaE9XbSWlIhIA1FYGPosFi8OVxYjR1a/6l9Va0kl0ochIiLpV1CQ2bKwSQ+rFRGRekIJQ0REYlHCEBGRWJQwREQkFiUMERGJRQlDRERiUcIQEZFYGuzEPTNbASxKOo4qtAc+SzqIKii+2lF8taP4aqc28XVx9wpXb22wCaOuM7OiymZT1gWKr3YUX+0ovtrJVHy6JSUiIrEoYYiISCxKGMkZl3QAm6D4akfx1Y7iq52MxKc+DBERiUVXGCIiEosShoiIxKKEkWZm1sfM3jezYjO7tILtfzSzd81srpk9b2ZdUratjaoQzjGzjJShjRHfmWa2IiWO/0nZNtDMFkSPgQnFd3NKbB+Y2Zcp27Jx/sab2XIzm1fJdjOzW6P455pZz5Rt2Th/m4qvIIprrpn928x+mbJtoZn9Jzp/Gak+FiO+g8zsq5R/xytTtlX52chSfBenxDYv+sxtHW3LxvnrbGYvmtl8M3vHzIZVsE/mPoPurkeaHkBT4EOgG9ACeBvYpdw+BwM50etzgQdTtn1TB+I7E7itgmO3Bj6KnreKXm+V7fjK7f97YHy2zl/0PQ4EegLzKtneF3gKMGAf4I1snb+Y8e1X9n2Bo8rii94vBNonfP4OAp6o7WcjU/GV2/cY4IUsn7/tgZ7R6zbABxX8H87YZ1BXGOnVCyh294/c/QdgEnBs6g7u/qK7l0ZvXwc61aX4qnAkMN3dP3f3L4DpQJ+E4zsNeCDNMVTJ3V8BPq9il2OBezx4HdjSzLYnO+dvk/G5+7+j7w/Z//zFOX+Vqc1nN7ZqxpfE5+9Td58dvV4FzAc6ltstY59BJYz06ggsSXlfwsb/mKkGEf4SKNPKzIrM7HUzOy7B+E6MLmUfMrPO1Tw2G/ER3crrCryQ0pzp8xdHZT9DNs5fdZX//DnwrJnNMrPBCcUEsK+ZvW1mT5nZrlFbnTp/ZpZD+GX7cEpzVs+fmeUBewJvlNuUsc+ganqnl1XQVuG4ZTM7HcgHfp3SnOvuS82sG/CCmf3H3T/McnyPAw+4+2ozGwJMBA6JeWw24itzKvCQu69Nacv0+Yujsp8hG+cvNjM7mJAwDkhp3j86f9sA083svegv7myaTVjL6Bsz6ws8CnSnjp0/wu2of7l76tVI1s6fmbUmJKsL3f3r8psrOCQtn0FdYaRXCdA55X0nYGn5nczsMGAE0N/dV5e1u/vS6Pkj4CXCXw9Zjc/dV6bEdBewV9xjsxFfilMpdzsgC+cvjsp+hmycv1jMbHfgbuBYd19Z1p5y/pYDjxBuA2WVu3/t7t9Er58EmptZe+rQ+YtU9fnL6Pkzs+aEZFHo7lMr2CVzn8FMdtA0tgfhiu0jwq2Sso65Xcvtsyeh8657ufatgJbR6/bAAtLcqRczvu1TXh8PvO7rO8w+juLcKnq9dbbji/bbmdDBaNk8fynfK4/KO22PZsMOxzezdf5ixpcLFAP7lWvfHGiT8vrfQJ8E4tuu7N+V8At3cXQuY302Mh1ftL0toZ9j82yfv+hc3AP8tYp9MvYZ1C2pNHL3NWZ2PvAMYVTHeHd/x8yuBorcfRpwA9AamGJmAIvdvT/QA/ibma0jXPmNcvd3E4jvAjPrD6wh/Kc4Mzr2czO7BpgZfbmrfcPL8WzFB6GzcZJH/wsiGT9/AGb2AGEkT3szKwGuAppH8d8JPEkYpVIMlAJnRdsyfv5ixncl0A4YG33+1nhY1XRb4JGorRlwv7s/nUB8JwHnmtka4Dvg1OjfucLPRgLxQfhD6ll3/zbl0KycP2B/4AzgP2Y2J2q7jPCHQMY/g1oaREREYlEfhoiIxKKEISIisShhiIhILEoYIiISixKGiIjEooQhIiKxKGGIJChailr/D6Ve0AdVJMvMLC+qZzCWsHZS500dI1IXaOKeSJZFq4x+RFie4/VkoxGJT1cYIslYpGQh9Y0Shkgyvt30LiJ1ixKGiIjEooQhIiKxqNNbRERi0RWGiIjEooQhIiKxKGGIiEgsShgiIhKLEoaIiMSihCEiIrEoYYiISCz/H6EMstm0om1oAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cosmo = FlatLambdaCDM(H0=70, Om0=0.2, name = 'FlatLambdaCDM')\n",
    "cosmo_params = cosmology(cosmo)\n",
    "\n",
    "z = 0.2\n",
    "radial_bins = 10\n",
    "y_err = 10\n",
    "N = 100\n",
    "\n",
    "Omega_M = 0.2\n",
    "little_h = 0.7\n",
    "\n",
    "M200 = 5e14*u.solMass\n",
    "# M200,R200,conc,rho_s, sigma_rho_s,r_s, sigma_r_s =  nfws_errors(M200_orig, 0.2, z,cosmo_params, cosmo.name)\n",
    "\n",
    "radius_array = np.linspace(0.1,2.0,radial_bins).round(3) #specify radius array for profiles. used in v_esc(r) funcs below.\n",
    "xdata = theta_data_array = radius_array / D_A(z, cosmo_params, cosmo.name).value\n",
    "\n",
    "ydata_err = np.repeat(y_err, len(ydata))\n",
    "\n",
    "r, truth = v_esc_NFW_M200(theta_data_array,z,M200.value,N,cosmo_params, cosmo.name)\n",
    "ydata = truth + np.random.normal(0,y_err,size=radial_bins) \n",
    "\n",
    "plt.plot(r, truth, \"r\", label = 'truth')\n",
    "plt.errorbar(r, ydata, yerr=ydata_err, fmt = \"bo\")\n",
    "plt.ylabel(\"v_esc_projected\")\n",
    "plt.xlabel(\"r\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "fitting M200 to v_esc_NFW_M200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-185765.7380385856"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def lnprior(theta):\n",
    "    p_M200 = theta[0]\n",
    "    \n",
    "    if not(1e13 < p_M200 < 1e17):\n",
    "          return -np.inf\n",
    "        \n",
    "    return 0.0\n",
    "\n",
    "\n",
    "def lnlike(theta, x, y, yerr): \n",
    "    p_M200 = theta[0]\n",
    "    p_theta_array = x\n",
    "    \n",
    "    ymodel = v_esc_NFW_M200(p_theta_array, z, p_M200, N, cosmo_params, 'FlatLambdaCDM')\n",
    "\n",
    "    inv_sigma2 = 1.0/(ydata_err**2)\n",
    "    return np.nan_to_num(-0.5*(np.sum((y-ymodel)**2*inv_sigma2)))\n",
    "\n",
    "def lnprob(theta, x, y, yerr):\n",
    "    lp = lnprior(theta)\n",
    "    ll = lnlike(theta, x, y, yerr)\n",
    "    \n",
    "    if not np.isfinite(lp):\n",
    "        return -np.inf\n",
    "    if not np.isfinite(ll):\n",
    "        return -np.inf    \n",
    "    \n",
    "    return lp + ll\n",
    "\n",
    "# lnprob([5e13], xdata, ydata, ydata_err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndim, nwalkers, nsteps = 1, 50, 1000\n",
    "p0 = np.transpose([np.random.uniform(1,10000,size=nwalkers)*1e13])#print np.shape(p0)\n",
    "\n",
    "pool = Pool(processes=20)     \n",
    "sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(xdata, ydata, ydata_err),pool=pool)\n",
    "sampler.run_mcmc(p0, nsteps)\n",
    "\n",
    "burn = 100\n",
    "samples = sampler.chain[:, burn:, :].reshape((-1, 1))\n",
    "print np.shape(sampler.chain)\n",
    "fig = corner.corner(samples[:,:], labels=[\"M200\"], truths = [M200.value])\n",
    "plt.show()\n",
    "percentile_array = np.arange(33-16.5,67+16.5, 1.0)\n",
    "print len(percentile_array)\n",
    "M200_fit =  np.percentile(sampler.chain[:,burn:,0],percentile_array)\n",
    "M200_fit_50 =  np.percentile(sampler.chain[:,burn:,0],50)\n",
    "M200_fit_33 =  np.percentile(sampler.chain[:,burn:,0],33-16.5)\n",
    "M200_fit_67 =  np.percentile(sampler.chain[:,burn:,0],67+16.5)\n",
    "print 'median(M200) = ', M200_fit_50, '+/-', M200_fit_67-M200_fit_50, M200_fit_50-M200_fit_33\n",
    "print 'median(logM200) = ', np.log10(M200_fit_50), '+/-', np.log10(M200_fit_67)-np.log10(M200_fit_50), np.log10(M200_fit_50) -np.log10(M200_fit_33)\n",
    "sigma_M200_fit = (M200_fit_67-M200_fit_50 + M200_fit_50-M200_fit_33)/2.0\n",
    "print 'Truth: ', np.log10(M200.value)"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}