joezuntz/wz_distance_modulus

## wz_distance_modulus
{
 "metadata": {
  "kernelspec": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "display_name": "IPython (Python 2)",
   "language": "python",
   "name": "python2"
  },
  "name": "",
  "signature": "sha256:9eddd098cb45009744eca0c194d2b0b94f7ff8e4c30907dd1925c325bee5d786"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.integrate"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#This all assumes omega_k = 0\n",
      "\n",
      "lcdm = {\n",
      "\"Omega_m\":  0.3,\n",
      "\"Omega_v\":  0.7,\n",
      "     \"H0\": 72.0,\n",
      "     \"w0\": -1.0,\n",
      "     \"wz\":  0.0,\n",
      "}\n",
      "\n",
      "\n",
      "wcdm = {\n",
      "\"Omega_m\":  0.3,\n",
      "\"Omega_v\":  0.7,\n",
      "     \"H0\": 72.0,\n",
      "     \"w0\": -1.0,\n",
      "     \"wz\":  0.5,\n",
      "}\n",
      "\n",
      "\n",
      "c = 3.0e5 #same units as H0\n",
      "\n",
      "\n",
      "def Omega_Lambda(z, cosmo):\n",
      "    return cosmo['Omega_v'] * (1+z)**(3*(1+cosmo['w0']-cosmo['wz'])) * exp(-3*cosmo['wz']*z)\n",
      "\n",
      "def Omega_Matter(z, cosmo):\n",
      "    return cosmo['Omega_m'] * (1+z)**3\n",
      "\n",
      "def H(z, cosmo):\n",
      "    return cosmo['H0'] * sqrt(Omega_Matter(z, cosmo) + Omega_Lambda(z, cosmo))\n",
      "\n",
      "def iH(z, cosmo):\n",
      "    return 1.0/H(z, cosmo)\n",
      "\n",
      "def D_c(z, cosmo):\n",
      "    return c*np.array([scipy.integrate.quad(iH, 0.0, zi, args=(cosmo,))[0] for zi in z])\n",
      "\n",
      "def D_m(z, cosmo):\n",
      "    return D_c(z, cosmo)\n",
      "\n",
      "def D_L(z, cosmo):\n",
      "    return (1+z)*D_m(z, cosmo)\n",
      "\n",
      "def mu(z, cosmo):\n",
      "    return 5*np.log10(D_L(z, cosmo)) - 25"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "z = np.linspace(0, 2.0, 201)\n",
      "plot(z, mu(z, lcdm))\n",
      "plot(z, mu(z, wcdm))\n",
      "#We get a warning about divide by zero because at z=0 mu is undefined"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "-c:46: RuntimeWarning: divide by zero encountered in log10\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "[<matplotlib.lines.Line2D at 0x1065d6bd0>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH6xJREFUeJzt3XmcVNWZ//HPAzTQiNKiCCqbAu4SFUXcWzMYdRINUUfj\nOEaTMY4mMTNxiUsy1ZXMZPk5Y2Ic/Y0zMRNHk6gxanQQCSrtNlFcEVlUBARFFBdAaVCWZ/44t+jq\n7qqm1r61fN+v131VVdetex/K8lunzj33XHN3RESktvWKuwARESk/hb2ISB1Q2IuI1AGFvYhIHVDY\ni4jUAYW9iEgdKEvYm9k1ZjbfzGab2d1mNqgc+xERkdyUq2X/J2Bfd/8M8CpwZZn2IyIiOShL2Lv7\nDHffHD18Ghhejv2IiEhueqLP/qvAAz2wHxERyaJPoS80sxnAsAxPXeXu90frXA186u6/LXQ/IiJS\nPCvX3Dhmdi5wPvBZd1+fZR1NzCMikid3t3xfU3DLvjtmdgJwGXBMtqBPKaRo6crMWty9Je46aoXe\nz9LS+1k6hTaSy9Vnfz0wEJhhZi+Y2Y1l2o+IiOSgLC17dx9Xju2KiEhhdAZt7WiNu4Aa0xp3ATWm\nNe4C6l3ZDtDmtHMzV5+9iEjuCs1NtexFROqAwl5EpA4o7EVE6oDCXkSkDijsRUTqgMJeRKQOKOxF\nROqAwl5EpA4o7EVE6oDCXkSkDijsRUTqgMJeRKQOKOxFROqAwl5EpA4o7EVE6oDCXkSkDijsRUTq\nQFnD3swuMbPNZja4nPsREZHuleWC4wBmNgKYDLxRrn2IiNQLS5oBOxT6+rKFPXAtcDnwxzLuQ0Sk\n6lnSegFDgOFZlhHArsD6QvdRlrA3s1OAN939JTNdT1xE6lfUIt8eGBktI9Lup8J8F2A18GanZUba\n/bc84WutxbyQOgoOezObAQzL8NTVwJXA8emrd7OdlrSHre7eWmhNIiI9zZLWnxDg6SHe+f5GYGna\nsgx4iRDiy4DlnvCMrXYzawbOArCWwhvP5l7Ql0T2DZrtBzwMtEV/Gg68BUx093c7revurqa/iFSk\nqHtlKNlDfCTQRAjtVIgv7XR/mSd8TUH7N/oSum9S+xoB9qNCcrPk3Tju/jLhzQHAzBYDE9z9g1Lv\nS0SkGFGYDwN2A0ZnWEYAq+gY4kuBJ2gP9Hc84Zvz3reR+iJJC/Iu93cE3k7b17K8/5Gp/ZW6Zd9l\nB2aLgIMzhb1a9iJSTmkt89FkDvSRwIfAkizLG9m6V7a6b6OJjr8AOgf5rnT8Ikm/Td1f4c6mjtst\nLDfLHvbd7lxhLyJFiA5+7giMoT3A00N9JOHA55Isyxue8HV579dIDYMcFe0ndZt+3+ga4ulh/qZ7\n/qNrFPYiUpMsaX0IoT0G2D26HZP2eAOwKFqW0DXM28hTFOY7kTnEU7cb0veT4f4qd0oesAp7Eala\nlrSBZA/zEcAK4PVoWZR2/3VP+Kq89xfCfCgdfwWkB/tIwiCTJWQO8zfcWZ3vfktBYS8iFc2S1gSM\nA/aIbsfSHu7bAotJC3Hag32JJ/yTvPdn9CcEd+pLI7WMIYR8W7TPxWQO848L+XeWm8JeRGJnSRtA\nCPFUqO+Rdn8A8CrwWrQspD3UV+Q7oiVqnQ+hPcA7B/oQQh956ktjUdr9xe4UNBwybgp7EekRlrQG\nQsu4c5jvQThYupgQ6qlgT92+7Yn8AseMPoTW+Vi6BvruwKd0DPH0+292HslSCxT2IlJSlrTtgb2B\nvdJu9yL0ob9FxzBP3V/qCc8rYM3oHW1zXNqS+hIZSeivf42ugb7Ynbz766udwl5E8hYNXRxBx0Df\nO1oagQXA/LTbVwgHRT/Naz/hBKJd6Brm4wi/Et6n46+A1LKokOGJtUxhLyJZRV0vY2kP8lSw7wl8\nRMdAT90uL6DbZYe07aYH+hhgDR2DfEvfvTt5D4+sVwp7EUmdMToa2K/TMo4wf8s8OrXU8x26GHW7\njKK9Wyd96Uf7L4AOoe7OR8X96wQU9iJ1Jep+2Zmuob4PoUvk5bRlLjA/3zNFzRhAaKF3DvRxwHt0\n/CWQWlaU40QiaaewF6lRlrRBwGfoGuwb6BjqLwPzPOF5newTzeGyb7TN9AOxQwnDI1NBngr2Vyt1\nDHo9UNiLVLmotT4KOIAQ7qnbnQhBPpu01ronOk4ZvtXtG9sSWv6pYN83WpoIrf+5hEBPLUtqcehi\ntVPYi1QRS1o/QtCmQv0AYDzhrM4Xo2V2dPt6PsMZo+6XvWgP9NTtToQQTwV7qotnqTt5T9Er8VDY\ni1QoS9p2wIHAwbQH+xjCWPH0YJ/tCV+Z83bDcMYxhC+MzwD7E4J9V8JB0VSYp4J9sVrq1U9hL1IB\nogm9UsE+IbodQQjz54AXCOE+L5950qMumPG0B/t4QrC/H217NjCHEOoL3dlQon+SVBiFvUgPs6Rt\nQwj2VKhPIPS5zwGeJYT7s4SRMBtz2mZorY+mY6h/hnA1pbm0B/tsYE49nkFa7xT2ImVkSetN6Pee\nBBwaLbsRAjg92Od5wnNqVUfj1fckfElMAA4iBPsaOob6bEJrXV0worAXKSVL2lBCoE+KloOB5cBT\n0TILeDnXaQOiCb32IQR6KtjHE64v+jzhy+J54AV3dL1myUphL1KgaGTMAbQH+yRgECHQt4S7J7pe\nRznj9oy+hF8BqWCfQOhfX0rXYI/lAhhSvSou7M3sW8BFwCZgqrt/N8M6CnvpcZa0wcARwJHAUYSu\nk1dpD/angNdymV896mMfS/gVMDG63Y8wM2Mq1J8DZmu6ACmFQnOzT5mKORY4GRjv7hvMbEg59iOS\nC0vaSEKop8J9JPA08DjwPUKrPaczQs3YifZQnwgcQrig9axom3cQWuxrS/zPEClKWVr2ZnYn8O/u\n/shW1lPLXkoqmghsHzqGez9CsD8R3c7OZXRMdHLSgbQfkJ1IONv0GUKwzwJmufNO6f8lIplVVDeO\nmb0A/BE4AVgPXOruz2ZYT2EvRYmmGNgTOBY4DmgmtLQfoz3cF+YyVa8ZOxO6d1LLvoRZIlPB/jRh\n9kadbSqx6fFuHDObQRj729nV0Xa3d/dJZnYIcCfhEmKZttOS9rDV3VsLrUnqgyVtN0KwH0cI+Q3A\nI8B9wD94wt/c6jZCX/u+dAz3JuB/gSeBS4Fn3MlrpkiRUjOzZkIjprjtlKllPw34ibs/Gj1eCBzq\n7u93Wk8te9kqS9qutAf7cYRumZmEgH8EWLy1lnvUJTOR9mA/jDBN75NpywK12qXSVdQBWuBewv+U\nj5rZHkDfzkEvko0lrZHQ134C8DnCL8iZ0XINsCCHcB9ICPXmaBkPvEQI9f8EzlNfu9STcrXsG4Bf\nEcYufwpckql7Ri17gS397nsTgv1zhJB+CZgOPAg8t7VZHzOE+/6EIY+t0fKUumSkFlTUAdqcd66w\nr1vRTJCTaW+9OyHcpwMPb+1SeQp3qVcKe6l4lrTRwBei5TDCwdBphIDvtmsmmm7gEMIXxGTCr8b0\ncH9a4S71QGEvFSeaPGwi7QE/FJgK3A/M8IRnPaPUDCPM1Z4K92MJ0w3MAP4EPOFOW1n/ASIVSGEv\nFcGSNgA4nnAG9V8CKwnhfh/hTNWsfe9mDCYc2D+eEPB9CeE+A3hIB1RFFPYSI0vatsBJwGmEoH6W\ncFLd/Z7wxVlfF8a6H0j4UjiJcObr47QH/Dx34vuAilQghb30KEvaIELXzGmE1vgTwB+AP3rC38v6\nOiN1YPYvgRMJZ7tOjZYn3fmkzKWLVDWFvZRdNILmi8AZhHlnWoG7CC34jKNnor73PWlvvR9CODA7\nFXjAndfLX7lI7VDYS1lEc72fCJxFGCLZCtwO/E+2A6xR98wkYArhy6E/8AAh4B/WjJAihVPYS8lE\nM0ceQwj4LxGuqfpb4K5sF/CILthxLCHgTyFMRXBPtLyovneR0qi06RKkClnS9gHOA75MGEXzW+AA\nT/iyjOuHE5tOIAT8icACQrgf7c5rPVK0iOREYV/non74M4CvAqOA/waO94TPy7h+CPgvRK85jnBV\np3uAy9xZ3iNFi0je1I1Th6K5aI4mBPwpwMOEuYymZ7qohxmNhIOrZxD67Z8kXJHpPnc+7Km6RUR9\n9pIDS9rOhG6arwKfADcDt3nC3+2yrtGPMGb+DODzhLHzdwB3u6MZTEViorCXjKJW/FHANwjh/Xvg\nl8AzneeiiYZJHgZ8hTB+fh4h4O9yZ0VP1i0imekArXQQndV6NnAR4b/zjcDXPeGru6xrjAb+BjgH\n2ATcAhzoztIeK1hEykphX2MsaWOBbwN/TbjYx7eBmRla8dsRWu9fIVye7/boNc9omKRI7VHY14Co\nq+YI4BLCma3/AYzvfC3WqJvmCODrhInKWoGfE85k1TQFIjVMffZVzJLWBzgV+A4wGPgZcIsnvMMZ\nqmZsT+imuYDwBX8TcKs7K3u2YhEplvrs60h0jdbzCSG/DPgxYX6aLdMHpx1svYAwvPIBwkHaR9VN\nI1J/FPZVxJI2EPg7QnfN08AZnvCnO6xjbEvoh78A6Efo0rlUrXiR+laWsDezicC/AQ3ARuAid3+m\nHPuqB9F0wt8kOtgKfM4T/lKHdYzdgW8RRtQ8AlwMtKoVLyJQvpb9/wO+7+7TzezE6PGxZdpXzbKk\nNQH/QOh+mQYc4wmfv+X50FXTTPgSOIJwFqyGTIpIF+UK+7eBQdH9JuCtMu2nJkV98t8CLiNc0m+S\nJ3zhlufDDJNnEb4IGoDrgL/W1MEikk1ZRuOY2SjClYsc6AUc5t515kSNxukoGl1zLpAAZgFXe8IX\nbHk+TEL2t4Q++/nAvwAz1FUjUj96fDSOmc0AhmV46mpCf/HF7n6PmZ1O6F6YnGU7LWkPW929tdCa\nqlU0Tv5LwD8Dy4FTPeGztjxv7EBo6V8EPAp80Z3n4qhVRHqWmTUTumuL206ZWvZr3H276L4Bq9x9\nUIb16r5lb0k7nDA+vgG4ApiROtvVjF2BSwmja/4AXOPOq3HVKiLxq7Rx9gvN7Bh3f5Qw57kCqhNL\n2nDgp4Sphq8AfucJ3wxgxjDgu4SRNb8G9nfXcQ8RKVy5wv7rwA1m1g9YFz0Wthx8vRT4e9onJ1sL\nYMaOwOWEfvlbgX0126SIlEJZwt7dnwUOLce2q5kl7STgBsLc8BM84Utgy3QGlwAXAncC4915M9t2\nRETypTNoe0B00ZDrgIOA8z3hDwGY0Z9w4PVy4I/ABHeWxFWniNQuhX0ZWdJ6EaYt+AFh8rGveMLX\nmdELOBP4EfAicKQ7r8RXqYjUOoV9mVjS9iIMOd0MNHvC5wKYcQxhfLwD57jzWHxViki9UNiXWNSa\nv5hwvkEC+HdP+GYzhgP/SjiWcQVwpzub46tUROqJ5rMvIUvabsB/Eb5Ez/WEL4ymNvgOYQTOjcBP\n3GmLsUwRqWKVNs6+rkRnwP4toQ/+p8DPPOGbzDgeuB54DTjUnddjLFNE6pha9kWKph/+JTAWONsT\nPjcaL/9z4HDgYnf+J84aRaR2FJqbvcpRTL2wpB0CPA+8AxxGi88z48vAy9Hf9lfQi0glUDdOAaJu\nm9RB2Is84XeZMQL4/8Ao4GR3ZnW3DRGRnqSWfZ6ibpu7gbOBSbT4H8w4h9DCn0U4MUpBLyIVRS37\nPFjSxgL3EaYZPpMWHwDcDuwH/IU7s+OsT0QkG7Xsc2RJ+yzwJHC9J/xCWvxwYDahb/5gBb2IVDKN\nxtmKqH/+G8D3gC/T4k8APwT+BviaOw/GWZ+I1BeNsy8DS1oDYZz8kcDhtPg64BFgDXCAOyvjrE9E\nJFfqxsnCkrYNcC8wkhD0uxKmJv4T8AUFvYhUE4V9Bpa0HYGHgZWs3PsUWvw84C5Ct80PNaeNiFQb\ndeN0YkkbBUwH7uGat3/I2mH/DewBTHJncbzViYgURmGfxpK2O6FP/jpa/HbCEMsFhPnm18VanIhI\nEdSNE7GkjQNagZ/S4g8BTxGuHnWOgl5Eql3BYW9mp5vZXDPbZGYHdXruSjN7zcwWmNnxxZdZXpa0\nvYGZQJIWX0jor7/CnX9yJ76xqSIiJVJMN84cYArhcntbmNk+wBnAPsCuwENmtoe7V+RBzeis2IeA\nq2jxNuA24FR3Ho+3MhGR0im4Ze/uC9z91QxPnQL8zt03uPsSYCEwsdD9lJMlbTgwg9Ci70eYlvh4\nBb2I1Jpy9NnvAryZ9vhNQgu/oljShhCC/kZafHvCpQKP0bQHIlKLuu3GMbMZwLAMT13l7vfnsZ+s\n/d5m1pL2sNXdW/PYbkEsadsCDwJ30+K9gfOAo9x5q9z7FhHJh5k1A83FbqfbsHf3yQVs8y1gRNrj\n4dHfsu2jpYB9FMyS1ge4A3iO5MZVwPlAszvLe7IOEZFcRA3g1tRjM0sUsp1SdeOkT8pzH3CmmfU1\ns92AcVBR87tfCzTwzx+/hve+ADhWQS8ita6YoZdTzGwZMAmYambTANx9HnAnMA+YBlzkcU6tmcaS\n9i3gL/jFq/ezYZsLCUGvrhsRqXl1M8WxJe0E4L+44/ffZ/5pPwSOdue1nti3iEipaIrjbkTz3dzC\nrG+0MP+0HwMnKehFpJ7UfMvektYPeIyVe7Vyw/yvAOfqgiMiUq0Kzc16mBvnX9jYbyU3zJ0CXKWg\nF5F6VNNhb0k7E+dErl3aF3o96M6v4q5JRCQONRv2UT/99dxxTyttO/UFLom7JhGRuNRkn70lrRcw\ng9cnv82tfzoKOMSdd0u9HxGRnqbROB1dyIbGwfxm6njgBAW9iNS7muvGsaSNxUly85PO5oYfufNc\n3DWJiMStplr2lrTewK954avPsOLA3sB1cdckIlIJairsga+xrmk77v+PscBB7lTkBVNERHpazYS9\nJW0wbv/ErdM/xXtfqMnNRETa1UzYAz/g9eOXs3zibHfuibsYEZFKUhNDLy1p49nUp5V/Xb6JtiH7\nuLOyFPWJiFSauh16aUkznOuZ+YMPaRvyUwW9iEhXtTD08q9oG7IbT172LvDLuIsREalEVR32lrQG\nNvf+MXffui3e50KNvhERyayqwx44l5X79OX1z93izotxFyMiUqmq9gCtJa0/mxoW8euZjSw7Yow7\nH5S6PhGRSlOP89mfy/IJxrIjfqSgFxHpXlFhb2anm9lcM9tkZhPS/j7ZzJ41s5ei22OLLzVtv0nr\nw8Z+3+ehH/cC/q2U2xYRqUXFDr2cA0wBbgLS+4NWAp939xVmti8wHRhe5L7SncnKvfvzRvP33FlX\nwu2KiNSkosLe3RcAmFnnv6cfLJ0HNJpZg7tvKGZ/EI2r/3RACzN/sBl05SkRkVz0xElVpwLPlSLo\nI8exdqedeO2kq935pETbFBGpaVsNezObAQzL8NRV7n7/Vl67L/ATYHI367SkPWx199ZuC1rXdCV/\n/g54b51AJSI1z8yageait1OKoZdmNhO4xN2fT/vbcOBh4Fx3/3OW1+U1hMiSthufDJzHdYt+7muH\nXFl04SIiVaYS5sbZsnMzawKmAt/NFvQFWT/o73nxXKNtiC5KIiKSh2KHXk4xs2XAJGCqmU2Lnvom\nMAZImNkL0bJjUftKWgO26Txe/qup7qwoZlsiIvWmas6gtX9sOJkVB9zJfz5zlDvPlLs2EZFKVPtn\n0K4dejlzzloJPBt3KSIi1aYqwt6SthP9P5zI/C9d5058P0VERKpUdVy85KOdz+GNI2H1qJvjLkVE\npBpVRcueTQ3ns+TYp935MO5SRESqUcWHvSVtKI0fjOGVk38Rdy0iItWq8rtxPtztPJYfvImPdr0v\n7lJERKpVxbfsga+x5OgnNQ+OiEjhKrplb0lrYpsBu/Pq5y+LuxYRkWpW0WHPe3ueyuqRm1k9+oG4\nSxERqWaV3Y2zofEclk+Y686ncZciIlLNKjbsLWnGoKWH8Oak2+OuRUSk2lVuN87aHQ/i02378cop\nt8RdiohItavcsP9g7Pm8v+f77rwddykiItWuYrtx6Lt2Mu/t0Rp3GSIitaAiw96S1pumxaNYfoj6\n60VESqAyu3HeH3cw0ItFkx+MuxQRkVpQmWG/euRZrGta6U5b3KWIiNSCiuzGod+a4/hgnC5SIiJS\nIhUX9pY0o2nJON7d7964axERqRUFh72ZnW5mc81sk5kdlOH5kWb2sZldkteGV43aHbe+vPr5uwqt\nTUREOiqmz34OMAW4Kcvz1wJT897qO/udRsP6j339IF2oRESkRAoOe3dfAGDW9SLnZvZFYBGwNu8N\n995wFGt2XVRoXSIi0lXJ++zNbCBwOdBS0AYaP9iXj3Z+vpQ1iYjUu25b9mY2AxiW4amr3P3+LC9r\nAX7m7m2WqdnfdR8taQ9buWTYzqza/aGtvU5EpB6YWTPQXPR23L3YQmYCl7j789Hjx4AR0dNNwGbg\n++5+Y4bXurtv+UKwi8ftwMAV73HDvMG+eoT67EVEOumcm7kq1UlVW3bs7kenFZUAPsoU9Bmt3OdE\n1g/6REEvIlJaxQy9nGJmy4BJwFQzm1Z8Of5ZVo9cXvx2REQkXTGjce4B7tnKOsm8Ntp/9UG8s//L\nhdYkIiKZVdYZtANXjObjnR+LuwwRkVpTMROhhcsQNm7L+3tMj7sWEZFaUzkt+0XHjWTDAGPe6fPi\nLkVEpNZUTMueNcMPY8D7693ZFHcpIiK1pnJa9r02TmDtThpyKSJSBpUT9g1te9O2gy4uLiJSBpUT\n9v1W78Yng16PuwwRkVpUOWHf+OFQ1jdpjL2ISBlUTtgPfGc71u74TNxliIjUoooIe/v22L40vt/A\nW4f+Oe5aRERqUUWEPcsOP5iPdtnkbxy1Ku5SRERqUWWE/ebek/h46MdxlyEiUqsqI+z7rB9P25CV\ncZchIlKrKiPs+348jvVNy+IuQ0SkVlVG2DesG8KGAQp7EZEyqZCwX7sdG/u/FXcZIiK1qjLCvu9H\nA9jYf2ncZYiI1KrKCPt+a/qzoXFJ3GWIiNSqYq5Be7qZzTWzTWZ2UKfnxpvZn83sZTN7ycz6dbux\nxlV9WN+0sNBaRESke8W07OcAU4AOlxE0sz7ArcDX3X0/4BhgQ7aN2Bmn9afPOmPpkUuKqEVERLpR\nzAXHFwCYWeenjgdecvc50Xrdz1H/8dAxrBvsvvzgjYXWIiIi3StHn/04wM3sQTN7zswu63Ztt91Z\n35S15S8iIsXrtmVvZjOAYRmeusrd78/ysgbgSOBgYB3wsJk95+6PZFy718bRfLLd+txLFhGRfHUb\n9u4+uYBtLgMec/cPAMzsAeAgIHPYvzjjTN7Y2MvMWoBWd28tYJ8iIjXJzJqB5mK3U6oLjqd33E8H\nLjezRsKB2WOAa7O+cuLhi2n8cLDPeaOlRLWIiNSMqAHcmnpsZolCtlPM0MspZrYMmARMNbNpUWGr\nCOH+DPAC8Jy7T8u6oT7rd2JD4weF1iEiIltXzGice4B7sjz3G+A3OW2oYd2OtA3WvDgiImUU/xm0\nDW1NbBzwdtxliIjUsgoI+7XbahI0EZHyij/s+300gI39NAmaiEgZVUDYr+nLxsZFcZchIlLL4g/7\nxg/6sG57TYImIlJG8Yc9Do0froi7ChGRWhZ/2K8bvNlntnjcZYiI1LL4w/6TQZ/GXYKISK2LP+wX\nHzsr7hJERGqducfXg2Jm7u5dJsQXEZHMCs3N+Fv2IiJSdgp7EZE6oLAXEakDCnsRkTqgsBcRqQMK\nexGROqCwFxGpAwp7EZE6oLAXEakDxVxw/HQzm2tmm8zsoLS/9zez35nZS2Y2z8yuKE2pIiJSqGJa\n9nOAKcBjnf5+JoC7jwcmABeY2cgi9iM5MLPmuGuoJXo/S0vvZ/wKDnt3X+Dur2Z46m1gGzPrDWwD\nfAqsKXQ/krPmuAuoMc1xF1BjmuMuoN6VvM/e3acTwv1tYAlwjbuvKvV+REQkd326e9LMZgDDMjx1\nlbvfn+U1ZwONwM7AYOBxM3vY3RcXW6yIiBSm6CmOzWwmcIm7Px89vhH4X3e/LXp8M/Cgu/8+w2t1\nhSoRkTwVMsVxty37PKTveAFwHHCbmW0DTAJ+lulFmsteRKRnFDP0coqZLSOE+VQzmxY9dRPQ18zm\nALOAX7n7y8WXKiIihYr1SlUiItIzyn4GrZmdYGYLzOw1M/tulnV+ET0/28wOLHdN1Wxr76eZNZvZ\najN7IVq+F0ed1cDMfmVm70S/QrOto89mjrb2fuqzmTszG2FmM6MTV182s4uzrJf759Pdy7YAvYGF\nwGigAXgR2LvTOicBD0T3DwWeKmdN1bzk+H42A/fFXWs1LMBRwIHAnCzP67NZ2vdTn83c38thwAHR\n/YHAK8VmZ7lb9hOBhe6+xN03ALcDp3Ra52TgFgB3fxpoMrOhZa6rWuXyfkLHA+aShbs/DnzYzSr6\nbOYhh/cT9NnMibuvcPcXo/sfA/OBXTqtltfns9xhvyuwLO3xm9HftrbO8DLXVa1yeT8dODz6WfeA\nme3TY9XVHn02S0ufzQKY2WjCL6anOz2V1+ezVEMvs8n16G/nb3sdNc4sl/fleWCEu7eZ2YnAvcAe\n5S2rpumzWTr6bObJzAYCdwHfjlr4XVbp9Djr57PcLfu3gBFpj0cQvn26W2d49Dfpaqvvp7t/5O5t\n0f1pQIOZDe65EmuKPpslpM9mfsysAfgDcJu735thlbw+n+UO+2eBcWY22sz6AmcA93Va5z7gHAAz\nmwSscvd3ylxXtdrq+2lmQ83MovsTCcNrP+j5UmuCPpslpM9m7qL36WZgnrv/PMtqeX0+y9qN4+4b\nzeybwHTCSJKb3X2+mV0QPX+Tuz9gZieZ2UJgLXBeOWuqZrm8n8BpwIVmthFoI5pyWroys98BxwA7\nRicIJgijnPTZLMDW3k/02czHEcDZwEtm9kL0t6uAkVDY51MnVYmI1AFdllBEpA4o7EVE6oDCXkSk\nDijsRUTqgMJeRKQOKOxFROqAwl5EpA4o7EVE6sD/AVEfKOsZTcI2AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1065d6b50>"
       ]
      }
     ],
     "prompt_number": 4
    }
   ],
   "metadata": {}
  }
 ]
}
	{
	"metadata": {
	"kernelspec": {
	"codemirror_mode": {
	"name": "ipython",
	"version": 2
	},
	"display_name": "IPython (Python 2)",
	"language": "python",
	"name": "python2"
	},
	"name": "",
	"signature": "sha256:9eddd098cb45009744eca0c194d2b0b94f7ff8e4c30907dd1925c325bee5d786"
	},
	"nbformat": 3,
	"nbformat_minor": 0,
	"worksheets": [
	{
	"cells": [
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"pylab inline"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "stream",
	"stream": "stdout",
	"text": [
	"Populating the interactive namespace from numpy and matplotlib\n"
	]
	}
	],
	"prompt_number": 1
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"import scipy.integrate"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 2
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"#This all assumes omega_k = 0\n",
	"\n",
	"lcdm = {\n",
	"\"Omega_m\": 0.3,\n",
	"\"Omega_v\": 0.7,\n",
	" \"H0\": 72.0,\n",
	" \"w0\": -1.0,\n",
	" \"wz\": 0.0,\n",
	"}\n",
	"\n",
	"\n",
	"wcdm = {\n",
	"\"Omega_m\": 0.3,\n",
	"\"Omega_v\": 0.7,\n",
	" \"H0\": 72.0,\n",
	" \"w0\": -1.0,\n",
	" \"wz\": 0.5,\n",
	"}\n",
	"\n",
	"\n",
	"c = 3.0e5 #same units as H0\n",
	"\n",
	"\n",
	"def Omega_Lambda(z, cosmo):\n",
	" return cosmo['Omega_v'] * (1+z)*(3(1+cosmo['w0']-cosmo['wz'])) * exp(-3cosmo['wz']z)\n",
	"\n",
	"def Omega_Matter(z, cosmo):\n",
	" return cosmo['Omega_m'] * (1+z)**3\n",
	"\n",
	"def H(z, cosmo):\n",
	" return cosmo['H0'] * sqrt(Omega_Matter(z, cosmo) + Omega_Lambda(z, cosmo))\n",
	"\n",
	"def iH(z, cosmo):\n",
	" return 1.0/H(z, cosmo)\n",
	"\n",
	"def D_c(z, cosmo):\n",
	" return c*np.array([scipy.integrate.quad(iH, 0.0, zi, args=(cosmo,))[0] for zi in z])\n",
	"\n",
	"def D_m(z, cosmo):\n",
	" return D_c(z, cosmo)\n",
	"\n",
	"def D_L(z, cosmo):\n",
	" return (1+z)*D_m(z, cosmo)\n",
	"\n",
	"def mu(z, cosmo):\n",
	" return 5*np.log10(D_L(z, cosmo)) - 25"
	],
	"language": "python",
	"metadata": {},
	"outputs": [],
	"prompt_number": 3
	},
	{
	"cell_type": "code",
	"collapsed": false,
	"input": [
	"z = np.linspace(0, 2.0, 201)\n",
	"plot(z, mu(z, lcdm))\n",
	"plot(z, mu(z, wcdm))\n",
	"#We get a warning about divide by zero because at z=0 mu is undefined"
	],
	"language": "python",
	"metadata": {},
	"outputs": [
	{
	"output_type": "stream",
	"stream": "stderr",
	"text": [
	"-c:46: RuntimeWarning: divide by zero encountered in log10\n"
	]
	},
	{
	"metadata": {},
	"output_type": "pyout",
	"prompt_number": 4,
	"text": [
	"[<matplotlib.lines.Line2D at 0x1065d6bd0>]"
	]
	},
	{
	"metadata": {},
	"output_type": "display_data",
	"png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEACAYAAABS29YJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH6xJREFUeJzt3XmcVNWZ//HPAzTQiNKiCCqbAu4SFUXcWzMYdRINUUfj\nOEaTMY4mMTNxiUsy1ZXMZPk5Y2Ic/Y0zMRNHk6gxanQQCSrtNlFcEVlUBARFFBdAaVCWZ/44t+jq\n7qqm1r61fN+v131VVdetex/K8lunzj33XHN3RESktvWKuwARESk/hb2ISB1Q2IuI1AGFvYhIHVDY\ni4jUAYW9iEgdKEvYm9k1ZjbfzGab2d1mNqgc+xERkdyUq2X/J2Bfd/8M8CpwZZn2IyIiOShL2Lv7\nDHffHD18Ghhejv2IiEhueqLP/qvAAz2wHxERyaJPoS80sxnAsAxPXeXu90frXA186u6/LXQ/IiJS\nPCvX3Dhmdi5wPvBZd1+fZR1NzCMikid3t3xfU3DLvjtmdgJwGXBMtqBPKaRo6crMWty9Je46aoXe\nz9LS+1k6hTaSy9Vnfz0wEJhhZi+Y2Y1l2o+IiOSgLC17dx9Xju2KiEhhdAZt7WiNu4Aa0xp3ATWm\nNe4C6l3ZDtDmtHMzV5+9iEjuCs1NtexFROqAwl5EpA4o7EVE6oDCXkSkDijsRUTqgMJeRKQOKOxF\nROqAwl5EpA4o7EVE6oDCXkSkDijsRUTqgMJeRKQOKOxFROqAwl5EpA4o7EVE6oDCXkSkDijsRUTq\nQFnD3swuMbPNZja4nPsREZHuleWC4wBmNgKYDLxRrn2IiNQLS5oBOxT6+rKFPXAtcDnwxzLuQ0Sk\n6lnSegFDgOFZlhHArsD6QvdRlrA3s1OAN939JTNdT1xE6lfUIt8eGBktI9Lup8J8F2A18GanZUba\n/bc84WutxbyQOgoOezObAQzL8NTVwJXA8emrd7OdlrSHre7eWmhNIiI9zZLWnxDg6SHe+f5GYGna\nsgx4iRDiy4DlnvCMrXYzawbOArCWwhvP5l7Ql0T2DZrtBzwMtEV/Gg68BUx093c7revurqa/iFSk\nqHtlKNlDfCTQRAjtVIgv7XR/mSd8TUH7N/oSum9S+xoB9qNCcrPk3Tju/jLhzQHAzBYDE9z9g1Lv\nS0SkGFGYDwN2A0ZnWEYAq+gY4kuBJ2gP9Hc84Zvz3reR+iJJC/Iu93cE3k7b17K8/5Gp/ZW6Zd9l\nB2aLgIMzhb1a9iJSTmkt89FkDvSRwIfAkizLG9m6V7a6b6OJjr8AOgf5rnT8Ikm/Td1f4c6mjtst\nLDfLHvbd7lxhLyJFiA5+7giMoT3A00N9JOHA55Isyxue8HV579dIDYMcFe0ndZt+3+ga4ulh/qZ7\n/qNrFPYiUpMsaX0IoT0G2D26HZP2eAOwKFqW0DXM28hTFOY7kTnEU7cb0veT4f4qd0oesAp7Eala\nlrSBZA/zEcAK4PVoWZR2/3VP+Kq89xfCfCgdfwWkB/tIwiCTJWQO8zfcWZ3vfktBYS8iFc2S1gSM\nA/aIbsfSHu7bAotJC3Hag32JJ/yTvPdn9CcEd+pLI7WMIYR8W7TPxWQO848L+XeWm8JeRGJnSRtA\nCPFUqO+Rdn8A8CrwWrQspD3UV+Q7oiVqnQ+hPcA7B/oQQh956ktjUdr9xe4UNBwybgp7EekRlrQG\nQsu4c5jvQThYupgQ6qlgT92+7Yn8AseMPoTW+Vi6BvruwKd0DPH0+292HslSCxT2IlJSlrTtgb2B\nvdJu9yL0ob9FxzBP3V/qCc8rYM3oHW1zXNqS+hIZSeivf42ugb7Ynbz766udwl5E8hYNXRxBx0Df\nO1oagQXA/LTbVwgHRT/Naz/hBKJd6Brm4wi/Et6n46+A1LKokOGJtUxhLyJZRV0vY2kP8lSw7wl8\nRMdAT90uL6DbZYe07aYH+hhgDR2DfEvfvTt5D4+sVwp7EUmdMToa2K/TMo4wf8s8OrXU8x26GHW7\njKK9Wyd96Uf7L4AOoe7OR8X96wQU9iJ1Jep+2Zmuob4PoUvk5bRlLjA/3zNFzRhAaKF3DvRxwHt0\n/CWQWlaU40QiaaewF6lRlrRBwGfoGuwb6BjqLwPzPOF5newTzeGyb7TN9AOxQwnDI1NBngr2Vyt1\nDHo9UNiLVLmotT4KOIAQ7qnbnQhBPpu01ronOk4ZvtXtG9sSWv6pYN83WpoIrf+5hEBPLUtqcehi\ntVPYi1QRS1o/QtCmQv0AYDzhrM4Xo2V2dPt6PsMZo+6XvWgP9NTtToQQTwV7qotnqTt5T9Er8VDY\ni1QoS9p2wIHAwbQH+xjCWPH0YJ/tCV+Z83bDcMYxhC+MzwD7E4J9V8JB0VSYp4J9sVrq1U9hL1IB\nogm9UsE+IbodQQjz54AXCOE+L5950qMumPG0B/t4QrC/H217NjCHEOoL3dlQon+SVBiFvUgPs6Rt\nQwj2VKhPIPS5zwGeJYT7s4SRMBtz2mZorY+mY6h/hnA1pbm0B/tsYE49nkFa7xT2ImVkSetN6Pee\nBBwaLbsRAjg92Od5wnNqVUfj1fckfElMAA4iBPsaOob6bEJrXV0worAXKSVL2lBCoE+KloOB5cBT\n0TILeDnXaQOiCb32IQR6KtjHE64v+jzhy+J54AV3dL1myUphL1KgaGTMAbQH+yRgECHQt4S7J7pe\nRznj9oy+hF8BqWCfQOhfX0rXYI/lAhhSvSou7M3sW8BFwCZgqrt/N8M6CnvpcZa0wcARwJHAUYSu\nk1dpD/angNdymV896mMfS/gVMDG63Y8wM2Mq1J8DZmu6ACmFQnOzT5mKORY4GRjv7hvMbEg59iOS\nC0vaSEKop8J9JPA08DjwPUKrPaczQs3YifZQnwgcQrig9axom3cQWuxrS/zPEClKWVr2ZnYn8O/u\n/shW1lPLXkoqmghsHzqGez9CsD8R3c7OZXRMdHLSgbQfkJ1IONv0GUKwzwJmufNO6f8lIplVVDeO\nmb0A/BE4AVgPXOruz2ZYT2EvRYmmGNgTOBY4DmgmtLQfoz3cF+YyVa8ZOxO6d1LLvoRZIlPB/jRh\n9kadbSqx6fFuHDObQRj729nV0Xa3d/dJZnYIcCfhEmKZttOS9rDV3VsLrUnqgyVtN0KwH0cI+Q3A\nI8B9wD94wt/c6jZCX/u+dAz3JuB/gSeBS4Fn3MlrpkiRUjOzZkIjprjtlKllPw34ibs/Gj1eCBzq\n7u93Wk8te9kqS9qutAf7cYRumZmEgH8EWLy1lnvUJTOR9mA/jDBN75NpywK12qXSVdQBWuBewv+U\nj5rZHkDfzkEvko0lrZHQ134C8DnCL8iZ0XINsCCHcB9ICPXmaBkPvEQI9f8EzlNfu9STcrXsG4Bf\nEcYufwpckql7Ri17gS397nsTgv1zhJB+CZgOPAg8t7VZHzOE+/6EIY+t0fKUumSkFlTUAdqcd66w\nr1vRTJCTaW+9OyHcpwMPb+1SeQp3qVcKe6l4lrTRwBei5TDCwdBphIDvtmsmmm7gEMIXxGTCr8b0\ncH9a4S71QGEvFSeaPGwi7QE/FJgK3A/M8IRnPaPUDCPM1Z4K92MJ0w3MAP4EPOFOW1n/ASIVSGEv\nFcGSNgA4nnAG9V8CKwnhfh/hTNWsfe9mDCYc2D+eEPB9CeE+A3hIB1RFFPYSI0vatsBJwGmEoH6W\ncFLd/Z7wxVlfF8a6H0j4UjiJcObr47QH/Dx34vuAilQghb30KEvaIELXzGmE1vgTwB+AP3rC38v6\nOiN1YPYvgRMJZ7tOjZYn3fmkzKWLVDWFvZRdNILmi8AZhHlnWoG7CC34jKNnor73PWlvvR9CODA7\nFXjAndfLX7lI7VDYS1lEc72fCJxFGCLZCtwO/E+2A6xR98wkYArhy6E/8AAh4B/WjJAihVPYS8lE\nM0ceQwj4LxGuqfpb4K5sF/CILthxLCHgTyFMRXBPtLyovneR0qi06RKkClnS9gHOA75MGEXzW+AA\nT/iyjOuHE5tOIAT8icACQrgf7c5rPVK0iOREYV/non74M4CvAqOA/waO94TPy7h+CPgvRK85jnBV\np3uAy9xZ3iNFi0je1I1Th6K5aI4mBPwpwMOEuYymZ7qohxmNhIOrZxD67Z8kXJHpPnc+7Km6RUR9\n9pIDS9rOhG6arwKfADcDt3nC3+2yrtGPMGb+DODzhLHzdwB3u6MZTEViorCXjKJW/FHANwjh/Xvg\nl8AzneeiiYZJHgZ8hTB+fh4h4O9yZ0VP1i0imekArXQQndV6NnAR4b/zjcDXPeGru6xrjAb+BjgH\n2ATcAhzoztIeK1hEykphX2MsaWOBbwN/TbjYx7eBmRla8dsRWu9fIVye7/boNc9omKRI7VHY14Co\nq+YI4BLCma3/AYzvfC3WqJvmCODrhInKWoGfE85k1TQFIjVMffZVzJLWBzgV+A4wGPgZcIsnvMMZ\nqmZsT+imuYDwBX8TcKs7K3u2YhEplvrs60h0jdbzCSG/DPgxYX6aLdMHpx1svYAwvPIBwkHaR9VN\nI1J/FPZVxJI2EPg7QnfN08AZnvCnO6xjbEvoh78A6Efo0rlUrXiR+laWsDezicC/AQ3ARuAid3+m\nHPuqB9F0wt8kOtgKfM4T/lKHdYzdgW8RRtQ8AlwMtKoVLyJQvpb9/wO+7+7TzezE6PGxZdpXzbKk\nNQH/QOh+mQYc4wmfv+X50FXTTPgSOIJwFqyGTIpIF+UK+7eBQdH9JuCtMu2nJkV98t8CLiNc0m+S\nJ3zhlufDDJNnEb4IGoDrgL/W1MEikk1ZRuOY2SjClYsc6AUc5t515kSNxukoGl1zLpAAZgFXe8IX\nbHk+TEL2t4Q++/nAvwAz1FUjUj96fDSOmc0AhmV46mpCf/HF7n6PmZ1O6F6YnGU7LWkPW929tdCa\nqlU0Tv5LwD8Dy4FTPeGztjxv7EBo6V8EPAp80Z3n4qhVRHqWmTUTumuL206ZWvZr3H276L4Bq9x9\nUIb16r5lb0k7nDA+vgG4ApiROtvVjF2BSwmja/4AXOPOq3HVKiLxq7Rx9gvN7Bh3f5Qw57kCqhNL\n2nDgp4Sphq8AfucJ3wxgxjDgu4SRNb8G9nfXcQ8RKVy5wv7rwA1m1g9YFz0Wthx8vRT4e9onJ1sL\nYMaOwOWEfvlbgX0126SIlEJZwt7dnwUOLce2q5kl7STgBsLc8BM84Utgy3QGlwAXAncC4915M9t2\nRETypTNoe0B00ZDrgIOA8z3hDwGY0Z9w4PVy4I/ABHeWxFWniNQuhX0ZWdJ6EaYt+AFh8rGveMLX\nmdELOBP4EfAicKQ7r8RXqYjUOoV9mVjS9iIMOd0MNHvC5wKYcQxhfLwD57jzWHxViki9UNiXWNSa\nv5hwvkEC+HdP+GYzhgP/SjiWcQVwpzub46tUROqJ5rMvIUvabsB/Eb5Ez/WEL4ymNvgOYQTOjcBP\n3GmLsUwRqWKVNs6+rkRnwP4toQ/+p8DPPOGbzDgeuB54DTjUnddjLFNE6pha9kWKph/+JTAWONsT\nPjcaL/9z4HDgYnf+J84aRaR2FJqbvcpRTL2wpB0CPA+8AxxGi88z48vAy9Hf9lfQi0glUDdOAaJu\nm9RB2Is84XeZMQL4/8Ao4GR3ZnW3DRGRnqSWfZ6ibpu7gbOBSbT4H8w4h9DCn0U4MUpBLyIVRS37\nPFjSxgL3EaYZPpMWHwDcDuwH/IU7s+OsT0QkG7Xsc2RJ+yzwJHC9J/xCWvxwYDahb/5gBb2IVDKN\nxtmKqH/+G8D3gC/T4k8APwT+BviaOw/GWZ+I1BeNsy8DS1oDYZz8kcDhtPg64BFgDXCAOyvjrE9E\nJFfqxsnCkrYNcC8wkhD0uxKmJv4T8AUFvYhUE4V9Bpa0HYGHgZWs3PsUWvw84C5Ct80PNaeNiFQb\ndeN0YkkbBUwH7uGat3/I2mH/DewBTHJncbzViYgURmGfxpK2O6FP/jpa/HbCEMsFhPnm18VanIhI\nEdSNE7GkjQNagZ/S4g8BTxGuHnWOgl5Eql3BYW9mp5vZXDPbZGYHdXruSjN7zcwWmNnxxZdZXpa0\nvYGZQJIWX0jor7/CnX9yJ76xqSIiJVJMN84cYArhcntbmNk+wBnAPsCuwENmtoe7V+RBzeis2IeA\nq2jxNuA24FR3Ho+3MhGR0im4Ze/uC9z91QxPnQL8zt03uPsSYCEwsdD9lJMlbTgwg9Ci70eYlvh4\nBb2I1Jpy9NnvAryZ9vhNQgu/oljShhCC/kZafHvCpQKP0bQHIlKLuu3GMbMZwLAMT13l7vfnsZ+s\n/d5m1pL2sNXdW/PYbkEsadsCDwJ30+K9gfOAo9x5q9z7FhHJh5k1A83FbqfbsHf3yQVs8y1gRNrj\n4dHfsu2jpYB9FMyS1ge4A3iO5MZVwPlAszvLe7IOEZFcRA3g1tRjM0sUsp1SdeOkT8pzH3CmmfU1\ns92AcVBR87tfCzTwzx+/hve+ADhWQS8ita6YoZdTzGwZMAmYambTANx9HnAnMA+YBlzkcU6tmcaS\n9i3gL/jFq/ezYZsLCUGvrhsRqXl1M8WxJe0E4L+44/ffZ/5pPwSOdue1nti3iEipaIrjbkTz3dzC\nrG+0MP+0HwMnKehFpJ7UfMvektYPeIyVe7Vyw/yvAOfqgiMiUq0Kzc16mBvnX9jYbyU3zJ0CXKWg\nF5F6VNNhb0k7E+dErl3aF3o96M6v4q5JRCQONRv2UT/99dxxTyttO/UFLom7JhGRuNRkn70lrRcw\ng9cnv82tfzoKOMSdd0u9HxGRnqbROB1dyIbGwfxm6njgBAW9iNS7muvGsaSNxUly85PO5oYfufNc\n3DWJiMStplr2lrTewK954avPsOLA3sB1cdckIlIJairsga+xrmk77v+PscBB7lTkBVNERHpazYS9\nJW0wbv/ErdM/xXtfqMnNRETa1UzYAz/g9eOXs3zibHfuibsYEZFKUhNDLy1p49nUp5V/Xb6JtiH7\nuLOyFPWJiFSauh16aUkznOuZ+YMPaRvyUwW9iEhXtTD08q9oG7IbT172LvDLuIsREalEVR32lrQG\nNvf+MXffui3e50KNvhERyayqwx44l5X79OX1z93izotxFyMiUqmq9gCtJa0/mxoW8euZjSw7Yow7\nH5S6PhGRSlOP89mfy/IJxrIjfqSgFxHpXlFhb2anm9lcM9tkZhPS/j7ZzJ41s5ei22OLLzVtv0nr\nw8Z+3+ehH/cC/q2U2xYRqUXFDr2cA0wBbgLS+4NWAp939xVmti8wHRhe5L7SncnKvfvzRvP33FlX\nwu2KiNSkosLe3RcAmFnnv6cfLJ0HNJpZg7tvKGZ/EI2r/3RACzN/sBl05SkRkVz0xElVpwLPlSLo\nI8exdqedeO2kq935pETbFBGpaVsNezObAQzL8NRV7n7/Vl67L/ATYHI367SkPWx199ZuC1rXdCV/\n/g54b51AJSI1z8yageait1OKoZdmNhO4xN2fT/vbcOBh4Fx3/3OW1+U1hMiSthufDJzHdYt+7muH\nXFl04SIiVaYS5sbZsnMzawKmAt/NFvQFWT/o73nxXKNtiC5KIiKSh2KHXk4xs2XAJGCqmU2Lnvom\nMAZImNkL0bJjUftKWgO26Txe/qup7qwoZlsiIvWmas6gtX9sOJkVB9zJfz5zlDvPlLs2EZFKVPtn\n0K4dejlzzloJPBt3KSIi1aYqwt6SthP9P5zI/C9d5058P0VERKpUdVy85KOdz+GNI2H1qJvjLkVE\npBpVRcueTQ3ns+TYp935MO5SRESqUcWHvSVtKI0fjOGVk38Rdy0iItWq8rtxPtztPJYfvImPdr0v\n7lJERKpVxbfsga+x5OgnNQ+OiEjhKrplb0lrYpsBu/Pq5y+LuxYRkWpW0WHPe3ueyuqRm1k9+oG4\nSxERqWaV3Y2zofEclk+Y686ncZciIlLNKjbsLWnGoKWH8Oak2+OuRUSk2lVuN87aHQ/i02378cop\nt8RdiohItavcsP9g7Pm8v+f77rwddykiItWuYrtx6Lt2Mu/t0Rp3GSIitaAiw96S1pumxaNYfoj6\n60VESqAyu3HeH3cw0ItFkx+MuxQRkVpQmWG/euRZrGta6U5b3KWIiNSCiuzGod+a4/hgnC5SIiJS\nIhUX9pY0o2nJON7d7964axERqRUFh72ZnW5mc81sk5kdlOH5kWb2sZldkteGV43aHbe+vPr5uwqt\nTUREOiqmz34OMAW4Kcvz1wJT897qO/udRsP6j339IF2oRESkRAoOe3dfAGDW9SLnZvZFYBGwNu8N\n995wFGt2XVRoXSIi0lXJ++zNbCBwOdBS0AYaP9iXj3Z+vpQ1iYjUu25b9mY2AxiW4amr3P3+LC9r\nAX7m7m2WqdnfdR8taQ9buWTYzqza/aGtvU5EpB6YWTPQXPR23L3YQmYCl7j789Hjx4AR0dNNwGbg\n++5+Y4bXurtv+UKwi8ftwMAV73HDvMG+eoT67EVEOumcm7kq1UlVW3bs7kenFZUAPsoU9Bmt3OdE\n1g/6REEvIlJaxQy9nGJmy4BJwFQzm1Z8Of5ZVo9cXvx2REQkXTGjce4B7tnKOsm8Ntp/9UG8s//L\nhdYkIiKZVdYZtANXjObjnR+LuwwRkVpTMROhhcsQNm7L+3tMj7sWEZFaUzkt+0XHjWTDAGPe6fPi\nLkVEpNZUTMueNcMPY8D7693ZFHcpIiK1pnJa9r02TmDtThpyKSJSBpUT9g1te9O2gy4uLiJSBpUT\n9v1W78Yng16PuwwRkVpUOWHf+OFQ1jdpjL2ISBlUTtgPfGc71u74TNxliIjUoooIe/v22L40vt/A\nW4f+Oe5aRERqUUWEPcsOP5iPdtnkbxy1Ku5SRERqUWWE/ebek/h46MdxlyEiUqsqI+z7rB9P25CV\ncZchIlKrKiPs+348jvVNy+IuQ0SkVlVG2DesG8KGAQp7EZEyqZCwX7sdG/u/FXcZIiK1qjLCvu9H\nA9jYf2ncZYiI1KrKCPt+a/qzoXFJ3GWIiNSqYq5Be7qZzTWzTWZ2UKfnxpvZn83sZTN7ycz6dbux\nxlV9WN+0sNBaRESke8W07OcAU4AOlxE0sz7ArcDX3X0/4BhgQ7aN2Bmn9afPOmPpkUuKqEVERLpR\nzAXHFwCYWeenjgdecvc50Xrdz1H/8dAxrBvsvvzgjYXWIiIi3StHn/04wM3sQTN7zswu63Ztt91Z\n35S15S8iIsXrtmVvZjOAYRmeusrd78/ysgbgSOBgYB3wsJk95+6PZFy718bRfLLd+txLFhGRfHUb\n9u4+uYBtLgMec/cPAMzsAeAgIHPYvzjjTN7Y2MvMWoBWd28tYJ8iIjXJzJqB5mK3U6oLjqd33E8H\nLjezRsKB2WOAa7O+cuLhi2n8cLDPeaOlRLWIiNSMqAHcmnpsZolCtlPM0MspZrYMmARMNbNpUWGr\nCOH+DPAC8Jy7T8u6oT7rd2JD4weF1iEiIltXzGice4B7sjz3G+A3OW2oYd2OtA3WvDgiImUU/xm0\nDW1NbBzwdtxliIjUsgoI+7XbahI0EZHyij/s+300gI39NAmaiEgZVUDYr+nLxsZFcZchIlLL4g/7\nxg/6sG57TYImIlJG8Yc9Do0froi7ChGRWhZ/2K8bvNlntnjcZYiI1LL4w/6TQZ/GXYKISK2LP+wX\nHzsr7hJERGqducfXg2Jm7u5dJsQXEZHMCs3N+Fv2IiJSdgp7EZE6oLAXEakDCnsRkTqgsBcRqQMK\nexGROqCwFxGpAwp7EZE6oLAXEakDxVxw/HQzm2tmm8zsoLS/9zez35nZS2Y2z8yuKE2pIiJSqGJa\n9nOAKcBjnf5+JoC7jwcmABeY2cgi9iM5MLPmuGuoJXo/S0vvZ/wKDnt3X+Dur2Z46m1gGzPrDWwD\nfAqsKXQ/krPmuAuoMc1xF1BjmuMuoN6VvM/e3acTwv1tYAlwjbuvKvV+REQkd326e9LMZgDDMjx1\nlbvfn+U1ZwONwM7AYOBxM3vY3RcXW6yIiBSm6CmOzWwmcIm7Px89vhH4X3e/LXp8M/Cgu/8+w2t1\nhSoRkTwVMsVxty37PKTveAFwHHCbmW0DTAJ+lulFmsteRKRnFDP0coqZLSOE+VQzmxY9dRPQ18zm\nALOAX7n7y8WXKiIihYr1SlUiItIzyn4GrZmdYGYLzOw1M/tulnV+ET0/28wOLHdN1Wxr76eZNZvZ\najN7IVq+F0ed1cDMfmVm70S/QrOto89mjrb2fuqzmTszG2FmM6MTV182s4uzrJf759Pdy7YAvYGF\nwGigAXgR2LvTOicBD0T3DwWeKmdN1bzk+H42A/fFXWs1LMBRwIHAnCzP67NZ2vdTn83c38thwAHR\n/YHAK8VmZ7lb9hOBhe6+xN03ALcDp3Ra52TgFgB3fxpoMrOhZa6rWuXyfkLHA+aShbs/DnzYzSr6\nbOYhh/cT9NnMibuvcPcXo/sfA/OBXTqtltfns9xhvyuwLO3xm9HftrbO8DLXVa1yeT8dODz6WfeA\nme3TY9XVHn02S0ufzQKY2WjCL6anOz2V1+ezVEMvs8n16G/nb3sdNc4sl/fleWCEu7eZ2YnAvcAe\n5S2rpumzWTr6bObJzAYCdwHfjlr4XVbp9Djr57PcLfu3gBFpj0cQvn26W2d49Dfpaqvvp7t/5O5t\n0f1pQIOZDe65EmuKPpslpM9mfsysAfgDcJu735thlbw+n+UO+2eBcWY22sz6AmcA93Va5z7gHAAz\nmwSscvd3ylxXtdrq+2lmQ83MovsTCcNrP+j5UmuCPpslpM9m7qL36WZgnrv/PMtqeX0+y9qN4+4b\nzeybwHTCSJKb3X2+mV0QPX+Tuz9gZieZ2UJgLXBeOWuqZrm8n8BpwIVmthFoI5pyWroys98BxwA7\nRicIJgijnPTZLMDW3k/02czHEcDZwEtm9kL0t6uAkVDY51MnVYmI1AFdllBEpA4o7EVE6oDCXkSk\nDijsRUTqgMJeRKQOKOxFROqAwl5EpA4o7EVE6sD/AVEfKOsZTcI2AAAAAElFTkSuQmCC\n",
	"text": [
	"<matplotlib.figure.Figure at 0x1065d6b50>"
	]
	}
	],
	"prompt_number": 4
	}
	],
	"metadata": {}
	}
	]
	}