problem4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import cosmology as cos\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "import seaborn as sea\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# set plotting parameters\n",
    "plt.rc('text', usetex=True)\n",
    "sea.set(style='white')\n",
    "sea.set(style='ticks')\n",
    "sea.set_style({'xtick.direct'\n",
    "               'ion': 'in','xtick.top':True,'xtick.minor.visible': True,\n",
    "               'ytick.direction': \"in\",'ytick.right': True,'ytick.minor.visible': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4 (a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# given cosmological parameters\n",
    "z = 4.25\n",
    "h = 0.70\n",
    "omega_matter = 0.27\n",
    "omega_lambda = 0.73\n",
    "H0 = 100*h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model = FlatLambdaCDM(H0, omega_matter) \n",
    "\n",
    "def vel_now(redshift):\n",
    "    return model.H0 * model.comoving_distance(redshift)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$2.7880043 \\times 10^{8} \\; \\mathrm{\\frac{m}{s}}$"
      ],
      "text/plain": [
       "<Quantity 2.78800426e+08 m / s>"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cos.recession_velocity(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$531097.96 \\; \\mathrm{\\frac{km}{s}}$"
      ],
      "text/plain": [
       "<Quantity 531097.96263321 km / s>"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vel_now(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "redshifts = np.linspace(0.0000001, 11.09, num=1000) # can't start at 0, because it'll result in infinity\n",
    "ratio = vel_now(redshifts)/(cos.recession_velocity(redshifts)*10**(-3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Let's take a look at the discrepancy between my value and the astropy value as a function of redshift. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAJLCAYAAACi3DFcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3XlcVfeB//83Oy5cwH25qMFd3BdO\nJCYmcQnZd2lWTFM17TdOOv21sZ2Zdh7TZdqkTtqZZDo1tmmg2VxidoOJmhgTFFRccRcVrvsGF5T1\n3vv74xaXCAoIfO7yej4ePoR7Dpd3TqK88/l8zueEeDwejwAAAOATQk0HAAAAwAWUMwAAAB9COQMA\nAPAhlDMAAAAfQjkDAADwIeGmAzSXiooKbdu2TZ07d1ZYWJjpOAAAAPVyuVw6ceKEhg4dqujo6EuO\nBUw527Ztmx577DHTMQAAABrszTff1NixYy95LWDKWefOnSV5/yGXLVum9PT0Fvk+kyZN0ooVK1rk\nvTMyMlost7++tz9eb6715fz1mvjj9fbX6+GP11ry32vij9fbX69Hfe999OhRPfbYY+f7y8UCppzV\nTmV269ZNNptNdru9xb5XS713S+b21/eW/O96c60v56/XxB+vt79eD3+81pL/XhN/vN7+ej2u9t51\nLcUKyBsCkpOTTUdokpbM7a/v3ZJaKjfX+nL+ek388Xr76/Xwx2st+e818cfr7a/XoynvHRIoj29y\nOBznh1Jb8v8GBg4cqF27drXY++NSXO/Ww7VuXVzv1sO1bl1c74a5Um8JyJEzAAAAf0U5a6Rnn33W\ndISgwvVuPVzr1sX1bj1c69bF9b52TGsCAAC0MqY1AQAA/ATlDAAAwIcYL2cLFizQggUL9Pvf//6K\n582fP7+VEgEAAJhjtJxlZ2crJSVFaWlpKioqUnZ2dr3n1XcMAAAgkBgtZxcXsoSEBBUVFZmMAwAA\nYJzRcpaWlqa0tDRJ0vbt2zV06NDLzsnPz1dKSkprRwMAADDC+JozyVvAhgwZoqSkpMuOlZSUGEgE\nAABghk88+Dw7O1s/+clPLnu9KaNmkyZNuuy1Z599VrNnz25yPgAAgKZ6+eWX9corrzT4fOPlbMGC\nBZoxY4akCzcIOJ1O2Ww2FRUVqaioSCUlJSopKVF+fn6do2sXYxNaAADgS2bPnn3ZIFHtJrR1MX63\n5ty5czV58mSNGzfu/OvTp0+XJKWmpio1NVWS5HQ6TUQEAABoVUZHzlJSUrRu3brLXl+yZMkln198\n4wAAAEAg84kbAgAAAOBFOQMAAPAhlDMAAAAfQjkDAACQpIICaeJEKSLC+3tBgZEYlDMAAABJ7unT\npa++kmpqvL8/9ZSRHMb3OQMAADClorJGG3Ye1zdbDutH32RfOmr1j+d/tzbKGQAACCrnKqq1fscx\nfbPlsNbvOK6qapck6d5eQzTwwNYLJxp6tjflDAAABLyz5dXK3X5U32w+rLxdx1Vd45Yk9ezcTjeM\n6KkbhvfQdT94T/rud70jZikp0t/+ZiQr5QwAAASksnNVWrvtqL7Zclibdh9XjcsjSerVLUY3DO+h\nG4b3UK9uMQoJCfnHV8RKq1aZC/wPlDMAABAwzpZXKyf/iFZvurSQJfaIVcqI7koZ1kMJXWMMp7wy\nyhkAAPBr5ZU1ys0/qtWbDl0yZZnYM1YTRvTQDSN6qEen9oZTNhzlDAAA+J2Kqhpt2HFcqzcd0rod\nx84v6u/dLUY3juypCSN7qmdn/ylkF6OcAQAAv1Bd49KGnd5Clpt/VBVV3kLWs3M7TRjZUzeO7Kne\n3WyGU147yhkAAPBZ1TVubd5zQqs3HdLabUd0rqJGktS1Q1vd9Y9Cdl0P20WL+v1fwJWzjIwM2Ww2\nJScny7Is03EAAEAjud0ebSs4qVV5h5S95bDKyqslSZ3i2miq1Vs3juyp/glxflnIcnJylJubK6fT\nWe85AVfO0tPTZbfbTccAAACN4PF4VHCoRKs2HtJXGx06VVIhSepgi9LdYxN144ieGtg7XqGh/lfI\nLmZZlizLksPhUGZmZp3nBFw5AwAA/uPoqbNatdGhVXkOFR0rkyS1iw7XlOReunmMXUmJnRTm54Ws\nsShnAACgVZWUVerrTYf0ZZ5DOw+ekSRFhIcqZXh33TzarjGDuioyIsxwSnMoZwAAoMWVV9YoZ9sR\nfZnn0MbdJ+R2exQSIo3o30k3j7Zr/LAeatcmwnRMn0A5AwAALaLG5dbGXce1Ku+Q1uYfUeU/tr7o\nZ4/VxNEJunFkD3WMbWM4pe+hnAEAgGZTu7B/5foirdroUElZlSSpe8d2mjjaromje8rexbcfn2Qa\n5QwAAFyz084KfbnBoZXrC3XwaKkkydYuUndNuE43j7ZrQK94v9z6wgTKGQAAaJLKapfWbj2ilRuK\ntGnXcbk9UnhYiMYP665JYxM0elBXRYSHmo7pdyhnAACgwdxuj7bvP6WV64v0zZbD53fsH9grXreM\nTdCNI3vK1i7ScEr/RjkDAABXdeTkWa1cX6QvNhTp2Olzkrw79t95w3W6dWxCy60jW75ceuEFaf16\n7+djx3o/Hz3a+/nDD0uLF0vz5kkzZ176dQ8/LJ050/TvnZfn/b4Xv29Djl0jyhkAAKjTuYpqrd50\nWCvWFWrHgdOSpOjIMN06NkG3jk3QsL6dWnbH/hdflH77W2n+fOnzz6XiYmnhQm/xqi1nkpSYKM2Z\n07xFKS/PW+7mzGncsWZAOQMAAOd5PB5t339an+Uc1DdbDquyyqWQEGlk/866ZWyCxg/rrjZRrVAf\nCgq85WfDhgtFLC6u7gL20EPe0bPFi70fX6tZs7zf//Tpxh1rJpQzAACgUyXlWrm+SJ/nFurIybOS\npK4d2mpKci/dOraXOse38n5kc+Z4i9bFI2RXO7/2a67VvHne38eMadyxZkI5AwAgSFXXuLV+x1F9\nllOovJ3H5PZIkeGhunmMXVOSe2loYgtPW15JXl7jpg1nzvSev3y5NHny1c9fvFiaMePC58XF3iK4\nYUPjszYzyhkAAEGm8KhTn+cW6osNRec3ie2fEKcpyb104yi72vvCY5QKCrxryRpj5kzvzQINKWcP\nPXRhlK2gQOrb17u2zQdQzgAACALexf2H9HlOoXYVeu9gjGkbqXtuStSU5N7q091mOOG3JCZ6S1Ot\n5culKVO8H0+e7L1B4Nt+9jMpPt476tYYU6ZcegeoYZQzAAAClMfj0a7CM8pac0CrNx1WVbVLoSHS\nmEFdNCW5t5KTuioiPMx0zLqNHu0tYLU3AEyeLHk80quvXlj39W21Nwz89rfehfu1pky5sBWHdOn2\nGnPmeIvg8883/z9DEwVcOcvIyJDNZlNycrIsyzIdBwCAVldWXq1VG4qUtfagDhxxSpK6dWyrycm9\nNGlsL3WK84OHjb/wgneqsbHTm3PmeL+udpRNqnuUTfKOxr36qrR//7VlbYScnBzl5ubK6XTWe07A\nlbP09HTZ7XbTMQAAaFW1o2TL1hzUV5sOqarapbDQEN0wvIdSx/fW8H6dzS3ub4rERG9BGzNGWrTI\nO3JWXCzt23f1r3voIe/XXklxsXevsvnzvSNurcSyLFmWJYfDoczMzDrPCbhyBgBAMDlbXq0v6xgl\nm2r11uTkXoqPiTac8Bo8/7x3evOFF7xFSvI+IeBnP7vy1/3sZ95Sd6XS9dvfXihotRITveVvzhzv\niFpx8YUtOlas8Ga50rFmEuLxeDzN9m4GORwOTZo0SStWrGDkDAAQ0Dwej3YXntGytd5Rssoq7yiZ\nNbSbUq/voxH9/WyULAhdqbcwcgYAgJ84W16tL/McWrb2gPYf9o6Sde3QVrdd31uTx/VSvM2PR8lw\nHuUMAAAft/9wiT75Zr++zHOcHyVLGd5dt13fRyMZJQs4lDMAAHxQdY1b2VsO65Nv9p9/6HiX+Da6\nbVIfTU7upQ6MkgUsyhkAAD7kZHG5stYc0LKcgyourZQkjR7URXemXKcxg7sqjFGygEc5AwDAMI/H\noy17TuqT7P3KyT8qt9ujdm0idN/Evro9pY96dGpvOiJaEeUMAABDzpZXa+X6Ii3N3i/H8TJJUmLP\nWN15w3W6aVRPRUfyYzoY8W8dAIBWdvCIU598s19fbChSRZVL4WGhunmMXXemXKeBveMVEsLUZTCj\nnAEA0Apcbo9y84/qo9UF2rrvpCSpc3wbTZvcR1OSeysuJspwQvgKyhkAAC3obHm1Ps8t1MdfF+jY\n6XOSpJH9O+vOCddp3JBuLPDHZShnAAC0gMMnyvTR1wVasa5Q5ZUuRUaE6bbre+vuGxPVu5vNdDz4\nMMoZAADNxOPxaPOeE/rgqwJt2HlMHo/UKTZaD08aoNuu7yNbu0jTEeEHKGcAAFyjiqoafbnBoY++\nLlDh0VJJ0qDe8brnpr4aP6y7wsNCDSeEP6GcAQDQRCeLy/XJN/u1bO0BlZ6rVlhoiCaOsuuemxI1\noFe86XjwU5QzAAAaaZ+jWO+v2qfVmw7J5fbI1i5SaZMH6PaUPuoY28Z0PPi5gCtnGRkZstlsSk5O\nlmVZpuMAAAKEx+PRhp3H9d6Xe7Vlr3crjN7dYnTvTX01cbRdkRFhhhPCH+Tk5Cg3N1dOp7PecwKu\nnKWnp8tut5uOAQAIENU1Lq3Kc+i9VfvOrycb2b+z7r+5n0YN7MyGsWgUy7JkWZYcDocyMzPrPCfg\nyhkAAM2h7FyVPl1zQB+tLtCZ0kqFhYbo5jF23T+xnxJ7xpqOhwBGOQMA4CJHT53Vh6sL9HnOQVVU\nudQmKlz339xPd09IVOd41pOh5VHOAACQtLvwjN77cq+ytxyW+x/7kz162yBNtXqrXZsI0/EQRChn\nAICg5fF4tHH3CS1esef88y4Te8Tq/pv7asLInuxPBiMoZwCAoONye7Rm62EtXrlH+xwlkqRRAzrr\nwVv6a3j/Tizyh1GUMwBA0KiucemLDQ4t+WKPDp04q5AQ6YYRPfTQrf3Vzx5nOh4giXIGAAgC5ZU1\nWrb2gN5ftU+nSioUHhaiqVZvPXBLP/Xs3N50POASlDMAQMBynq3SR6sL9PHXBSorr1Z0ZJjum9hX\n903sy07+8FmUMwBAwDlZXK73Vu3VsrUHVVnlUkzbSD162yDdecN1srWLNB0PuCLKGQAgYBw9dVaL\nV+7RinWFqnF51Ck2Wk/cPli3Wb0VHcWPPPgH/ksFAPi9wyfLtGj5Hq3cUCS326Mendrp4Un9NXF0\ngiLC2Q4D/oVyBgDwW0XHSrVwxW59leeQ2yMldG2vaZMH6saRPRUWynYY8E/Gy9mCBQskSYWFhfrJ\nT37S6OMAgOBz8KhTCz/frdWbD8njkfp0t2na5AFKGd6DUga/Z7ScZWdnKyUlRQkJCfqnf/qn8583\n9DgAILgUHCrRguW7lL3liCQpsWesvjNlgKyk7gqllCFAGC1nRUVFKioqUlpamhISElRUVNSo4wCA\n4LCn6IwWfL5bOflHJUn9E+L0nakDNW5wV3bzR8AxWs7S0tLOf7x9+3bdcccdjToOAAhse4uK9eay\nnVq/45gkaXCfDvrOlIEaNbAzpQwBy/iaM0nKz8/XkCFDlJSU1KTjF5s0adJlrz377LOaPXv2NecE\nALSO/YdL9NaynVq7zTtSlpTYUY9MHajh/XjuJfzPyy+/rFdeeaXB5/tEOcvOzr7iYv+rHb/YihUr\nZLfbmysaAKAVFR516q3PdumbzYclSYN6x+vx1ME8jBx+bfbs2ZcNEjkcjjoHlCQfKGcLFizQjBkz\nJF24AcDpdMpms9V7HAAQWA6fKNPbn+3Sqo0OeTxSv4Q4PZ46SKMHdqGUIegY3ZkvOztbc+fO1eTJ\nkzVu3Ljzr0+fPv2KxwEAgeHoqbP673c26vsvrtSXeQ716W7Tvz2VrJeeu0ljBrHYH8HJ6MhZSkqK\n1q1bd9nrS5YsueJxAIB/O3GmXAuW79Ly3EK53B4ldI3RY7cN0vhhbIkBGJ/WBAAEjzOlFVq4fLey\n1hxUjcutnp3b6ZGpgzSBHf2B8yhnAIAWd7a8Wu+t2qsPVu1TRZVL3Tq21XemDNTNo+0KC+PZl8DF\nKGcAgBZTVe3S0uz9Wrh8j0rPVSk+JkpP3Z2kqVZvhVPKgDpRzgAAzc7lcmvl+iK9tWynTpZUqF10\nuJ68Y7DunpCo6Ch+9ABXwp8QAECz8Xg8WrP1iP7+6Q45jpcpMjxUD97STw/e2l8xbSNNxwP8AuUM\nANAstuw9oYxPtmt3YbFCQ0N02/W99Z0pA9Upro3paIBfoZwBAK7JPkexMj7Zro27T0iSbhjRQ4+n\nDpK9S4zhZIB/opwBAJrkxJly/f3T7fpig0OSNLJ/Zz1xx2AN6BVvOBng3yhnAIBGOVdRrcUr9+iD\nVftUVePWdT1seuquJI0a2MV0NCAgUM4AAA1S43Jr2dqDevuznSopq1LH2Gg9cftg3TwmgQ1kgWZE\nOQMAXJHH41FO/lG9/nG+Dp04qzZRYXri9sG656ZERUfyYwRobvypAgDUa3fhGb32Ub7yC04pNDRE\nt6f00SNTByo+Jtp0NCBgUc4AAJc5dvqc/r50h1Zt9C72t5K6Kf3OIUroyh2YQEsLuHKWkZEhm82m\n5ORkWZZlOg4A+JVzFdVatGKPPvhqn6pr3Opnj9V37x6qYf06mY4GBIScnBzl5ubK6XTWe07AlbP0\n9HTZ7XbTMQDAr7jdHq1cX6TMpdt1prRSneLaKP2OwbpplF2hLPYHmo1lWbIsSw6HQ5mZmXWeE3Dl\nDADQODv2n9arH2zV3qJiRUaE6dGpA3X/Lf1Y7A8Ywp88AAhSJ4vL9frH28+vK7tpVE9NvzNJneN5\n3BJgEuUMAIJMRVWN3vtyn979Yo8qq1zqZ4/VjPuGach1HU1HAyDKGQAEDY/Ho683HdbfPsnXiTPl\nio+J0jP3D9etYxNYVwb4EMoZAASBfY5ivfr+Vm3ff1rhYaF66Nb+enhSf7WNjjAdDcC3UM4AIICV\nnavSG1k79Wn2frk90vhh3fXUXUnq3qmd6WgA6kE5A4AA5N0ao1Cvf7JdJWVV6tm5vWbdP4yHkwN+\ngHIGAAFmr6NYf16yRbsOnlF0ZJim3zlE99zUVxHhoaajAWgAyhkABIiyc1X6+6c7lLXmgNweacKI\nHnr6nqHqFMfWGIA/oZwBgJ9zuz1asc47hek8WyV7F+8U5sgBTGEC/ohyBgB+bK+jWH9+d4t2FTKF\nCQQKyhkA+KFzFdV6I2unPvm6gClMIMBQzgDAj3g8Hq3ZekTz3tuq084K9ezcTs88MJwpTCCAUM4A\nwE8cP31Of35vi9ZtP6bwsFA9OnWgHprUXxHhYaajAWhGlDMA8HEul1sfri7Qm8t2qrLKpeH9Oun7\nDw6XvUuM6WgAWgDlDAB82K6Dp/W/izdr/2GnYtpG6gcPDtctYxIUEsKzMIFAFXDlLCMjQzabTcnJ\nybIsy3QcAGiSs+XV+vunO7Q0e788HmlKci9NvytJtnaRpqMBuAY5OTnKzc2V0+ms95yAK2fp6emy\n2+2mYwBAk3g8HmVvOaJX39+i085K2bu01w8eGqFhfTuZjgagGViWJcuy5HA4lJmZWec5AVfOAMBf\nnSop15+XbNHabUcVER6qx1MH6YFb+rHgHwgylDMAMMzj8ejz3EK99uE2na2oUVJiR82eNlI9O7c3\nHQ2AAZQzADDo6KmzennhJm3Ze1JtosL1g4dG6Dart0JDWfAPBCvKGQAY4HJ79NHqAr2RtUOVVS6N\nHdxV/++hEezwD4ByBgCt7eARp15euEm7Cs8opm2knn14pCaO6sn2GAAkUc4AoNVU17i1eMVuLVyx\nWzUujyaOsmvGfUMV2z7KdDQAPoRyBgCtYK+jWH98O08Hj5aqY2y0fvDQCCUP6WY6FgAfRDkDgBZU\nXePWwuXe0TK326Pbru+tp+5KUrs2EaajAfBRlDMAaCH7D5foj29vVMHhEnWKa6N/mjZSowZ2MR0L\ngI+jnAFAM3O53Fr8xR6989ku1bg8mpLcS0/fM5TRMgANQjkDgGZUdKxUf3g7T3uKitXBFq3Z00Zq\n7OCupmMB8COUMwBoBi63Rx+s2qs3snaqusatm8fYNeu+YWrflgeVA2gcyhkAXKPDJ8r0x3c2aseB\n04prH6UfPDRC44d1Nx0LgJ+inAFAE3k8Hn265oD++mG+qqpdmjCih555YDj7lgG4JpQzAGiCM84K\n/c/CTVq/45hi2kboh98ZpRtH9jQdC0AAoJwBQCOt3XZELy/cJOfZKo0c0Fk//M4odYzlmZgAmkfA\nlbOMjAzZbDYlJyfLsizTcQAEkPLKGv3lg236LOegIsJDNeO+obrrhkSFhvJMTAANk5OTo9zcXDmd\nznrPCbhylp6eLrvdbjoGgACz8+BpvfRWno6cPKvEHrH60WOj1bubzXQsAH7GsixZliWHw6HMzMw6\nzwm4cgYAzcnlcmvB8t1asHy3PB6PHrylnx5LHaSI8DDT0QAEKMoZANTj8MkyvfRmnnYVnlHn+Db6\n50dGa1jfTqZjAQhwlDMAqMPK9UX685LNKq906eYxdj1z/3AevwSgVVDOAOAi5yqq9X9LtujLDQ61\njQ7Xjx8bo4mjWccKoPVQzgDgH/YWFevFN9bryMmzGtgrXj9+fIy6dWxnOhaAIEM5AxD03G6PPvhq\nnzKXbleNy7vo//HbBys8LNR0NABBiHIGIKgVl1bqD+/kKW/nccXFROlHj4zWqIFdTMcCEMQoZwCC\n1qbdx/XSW3k6U1qp0QO76IePjFJ8TLTpWACCHOUMQNCpcbn1ZtZOvfvFHoWGhOi7dyfp3pv6stM/\nAJ9AOQMQVE6cKdeLf1+nnQfPqHvHdvrx42M0oFe86VgAcB7lDEDQyNt5XHPf3KDSc1W6cWRPPfvw\nCLWNZu8yAL6FcgYg4LncHr29bKcWrtitsNBQff/B4bp9fB+FhDCNCcD3GC9nCxYskCQVFhbqJz/5\nyWXHs7KyZLPZlJ+frxkzZrR2PAB+7oyzQnPf3KAte0+qa4e2+umT49QvIc50LACol9FNfLKzs5WS\nkqK0tDQVFRUpOzv7kuP5+fmSpJSUlPMFDQAaauvek3rupS+1Ze9JWUnd9Md/nkgxA+DzjJaziwtZ\nQkKCioqKLjm+dOlS2Wy288e/Xd4AoC5ut0cLl+/Wv/35GznPVunpe5L0r08lq33bSNPRAOCqjE5r\npqWlnf94+/btuuOOOy45XlpaqtjY2POfFxcXt1o2AP7JebZKL721QRt2HlfH2Gg9/8RYDbmuo+lY\nANBgxtecSd7pyyFDhigpKema32vSpEmXvfbss89q9uzZ1/zeAHzb3qJi/WdGrk6cKdfogV30o0dH\nK7Z9lOlYAILcyy+/rFdeeaXB5/tEOcvOzq7zZoCYmBiVlJRIkpxOp+Lirr5WZMWKFbLb7c2eEYBv\nW7GuUP+7eLNqXG49etsgpU0ewKayAHzC7NmzLxskcjgcdQ4oSYbXnEneuzVr78KsXVPmdDolSXfc\nccf5dWhFRUVKSUkxExKAz6qucev/3t2sP76zUZERYfrF09frkakDKWYA/JbxuzXnzp2ryZMna9y4\ncedfnz59uiSdn+bMzs6WzWZrlmlPAIHjVEm5/uVPX2tp9gH16W7TSz+8SWMHdzUdCwCuidFpzZSU\nFK1bt+6y15csWXL+44tvGgCAWvkFp/S7zHUqLq3UTaN6avbDIxUd5RMrNQDgmvA3GQC/4vF49PHX\n+/XXD7fJI+l79w7VPTcmsts/gIBBOQPgNyqqavSnxZv1xQaH4tpH6fknx2pY306mYwFAs6KcAfAL\nx8+c02/+lquCQyUa2CteP00fp05xbUzHAoBmRzkD4PN27D+t/3w9V8VllZpq9dYzDwxTRHiY6VgA\n0CIoZwB82uc5B/WndzfL7ZFm3T9Md95wHevLAAQ0yhkAn+RyufXax/n68KsCtW8ToZ8+OU4jBnQ2\nHQsAWhzlDIDPKTtXpRf+vl6bdp9QQtcY/fy7lrp3amc6FgC0CsoZAJ9SdKxUv3otR0dOntW4IV31\n48fGqG10hOlYANBqKGcAfMb6Hcf0+zfW61xFjR66tb8ev32wwngME4AgQzkDYJzH49H7q/bpbx/n\nKyIsVP/fY2N082i76VgAYATlDIBRNS63/rxki5atPagOtmj961PJGtAr3nQsADCGcgbAmLLyar2Q\nsU6b9pxQYs9Y/eJpSx1j2VgWQHCjnAEw4uips/rlX3NUdKxUyUO66cePj1EbHlwOAIFXzjIyMmSz\n2ZScnCzLskzHAVCHnQdP69ev5aikrEr33tRXT92dxMJ/AEEhJydHubm5cjqd9Z4TcOUsPT1ddjsL\niQFf9fXmQ/rDW3mqcbn1/QeH646U60xHAoBWY1mWLMuSw+FQZmZmnecEXDkD4Js8Ho8Wr9yjzKU7\n1CYqXP/yVLLGDOpqOhYA+BzKGYAWV13j1p8Wb9bydYXqFNdG//6969Wnu810LADwSZQzAC3qXEW1\nfvu6947Mfglx+vl3LXWwRZuOBQA+i3IGoMWcKinXf/xlrfYfdspK6qYfPzZG0dyRCQBXFGo6AIDA\nVHSsVD95ebX2H3bq9vF99LP0cS1XzAoKpDFjpPh47+8FBS3zfRqib19p+XJz3x+A36OcAWh2+QWn\n9PzLq3XiTLmeuH2wvv/gcIWFNfGvm7w8b+l68cX6z5kyRZo/XzpzRkpL8xa04mLvsTFjpJCQS3/V\nHgMAH0Q5A9Cssrcc1s/nZau8skbPpY3StMkDFBJyDXuYzZkjTZ5c//Hly73HR4/2fv7881KHDtLC\nhd7Pi4ulDRskj+fCr7i4pucBgBbG4g8Azebjrwv06vtbFRURpn97ytLoQV2u7Q1ffdU7KrZuXf3n\nTJ4sJSZe+tro0dK+fdf2vQHAEEbOAFwzt9uj1z/O17z3tmpyYa7e+dPjGj1+oHc6MiTEO7XYWMXF\n0gsveEfCrubb5Swvz7v2q9acOd4sfft6jzVE7XuEhEgPP+wtifH/eCD7iy96j8XHe1+vb5q0vvPm\nzPF+fvH3io/3rpWbNcv762KsYwOCCiNnAK5JdY1b/7Nwo77c4FDPzu007V/+RWEdf+M9WFDgLRbz\n53s/nzPn6uu9Hn7YOxo2Y4bdakeqAAAgAElEQVT3/MZavNj7+8yZF16Li/OuR1u+3FsUz5y5+tTm\nww9L8+ZJY8dKkyZJP/uZ9Pnn3mOJiRdG5h5+2Jtz3rzL36O+8154wVvOXnzRm/Phh6VFi7znz5rl\nzVj7fnl50unTV57aBRBQKGcAmqyiqka/y1inDTuPa2DveP38u5Zi20ddOGHKFG8RqV0P9sILDXvj\n5cu9pWTRosYFKijwFqDaEiV5P64dWaudAl248NLyVt971RaitDTv1OpDD3k/r/299lh9JfJK5y1a\n5C1hCxZ4C1nt9xo92vvr1Ve9GefNu3pWAAGFcgagScrOVemXf83RjgOnNWZQF/00fZyiIy/6K2XO\nHG8Rasi05LctWuQdLaqdRiwu9o6IzZtX/1qy4uJLR6BqfXvKMy6uYevRRo++cLPB559fOtVYUOAt\nmuvXX3kk8ErnxcV5r9GsWd4bFi42a5b3n2PmTG+R/PZxAAGNcgagcQoKVP1kutqsWaMnuw/Smh/9\nRk9911L4xVtlLF/uHfnZv//Sr23otOa8eZdOE06Z4v11paI3aZJ3+rR2lK5WcfGlU5i1U61Xk5jo\nLVa1a8BqR8GKi70jXosWeTMuX375GrGGnFdb3F544UKprDVzpvda5eV5c3y7YAIIaJQzAI1S9fiT\nilzzjSRp6KF8JX3wkkJ+dO+FE2pHsObPv3xdV0OnNa+mdpPZ2tLy7enTWsuXe1+fN8977quvel9v\nyDRhXp53xKxDh8vLneRdiybVP/V6tfNqS+jkyd6PX3zx0vI5bZp33V1dxQ9AQONuTQANduCIU6E5\nay95LSQ7+9KTfvvbCwWtdtPXhoxUNcYLL1xYv7V8uffXlCmXbjRbu47r4jstFy1q3BRh7Z2Wtf8M\neXneAjh27IWnEdTnSufNmuU9VrvObP78C6NrtWpHzqZNa3heAAEhxOPxeEyHaA4Oh0OTJk3SihUr\nZLfbTccBAs6O/af1H39dq397/XkNPZR/4cBNN0mrVpkL1hJefNF7A0DtaFdxsXfUbcGC1lv/tXy5\nt7A19qYIAH7hSr0l4KY1MzIyZLPZlJycLMuyTMcBAsKGncf0n6+vU43LLecrf5b+8K9SdraUkiL9\n7W+m4zW/uDhvIatdr1Y7rdmhQ+tlqF2PBiCg5OTkKDc3V06ns95zAq6cpaenM3IGNKNvNh/W3DfX\nKzQkRP86PVnJSd2k+wJspOzbZs68sKBf8t45Om1a64xivfqqd0pz2rTL19AB8HuWZcmyLDkcDmVm\nZtZ5TsCVMwDN54sNRfrj23mKigzTz5++XsP6djIdqfU8/3zTtgG5VjNnsq8ZEOQoZwDqtGztAf3v\n4s1qGx2hX84crwG94k1HAoCgQDkDcJkPV+/T/Pe3ydYuUr+alaLEnrGmIwFA0KCcAbjEohW7lbl0\nh+JjovTrZ1LUq5vNdCQACCqUMwCSJI/HozeX7dSCz3erU1wb/eaZFPXo3N50LAAIOpQzAPJ4PPrb\nx9v13pd71a1jW/3mmRvUpUNb07EAIChRzoAg53Z7NO+9LVqafUD2Lu3162dS1DG2jelYABC0KGdA\nEHO7PfrTu5u1bO1B9elu069mpSguJsp0LAAIapQzIEhdXMwSe8bqV7NSZGsXaToWAAQ9yhkQhL5d\nzH79TIpi2lLMAMAXhJoOAKB1UcwAwLdRzoAgQjEDAN9HOQOChNvt0f8t2UIxAwAfRzkDgkBtMcta\nc4BiBgA+jnIGBDiKGQD4F8oZEMA8Hu8GsxQzAPAfAbeVRkZGhmw2m5KTk2VZluk4gDEej0evfZSv\npdkHdF0PG8UMAHxATk6OcnNz5XQ66z0n4MpZenq67Ha76RiAcW8t26X3V+2TvUt7/XImxQwAfIFl\nWbIsSw6HQ5mZmXWew7QmEIAWrditdz7fpe4d2+nXz/BIJgDwJ5QzIMB8uHqfMpfuUKe4NjzEHAD8\nEOUMCCDL1h7Q/Pe3KT4mSr/5foq6dGhrOhIAoJEoZ0CA+GJDkf538WbZ2kXq18+kqEen9qYjAQCa\ngHIGBIBvthzWH9/ZqLbREfrlzPHq1c1mOhIAoIkoZ4Cfy9t5XHPfWK+oiFD9x4zr1dceZzoSAOAa\nUM4AP7bzwGn9Z0auQkNC9POnr9fA3h1MRwIAXCPKGeCnDh5x6j/+slbVNW49/8RYDevbyXQkAEAz\noJwBfujoqbP6xavZKiuv1nNpI2UN7W46EgCgmVDOAD9zxlmhX8xbo9POSn3v3qG6dWwv05EAAM2I\ncgb4kbLyav37/DU6cuqspk0eoHtv6ms6EgCgmflEOcvPz6/3WFZWlrKzs7VgwYJWTAT4noqqGv3q\nr2u1/7BTqeP76PHUQaYjAQBagPFylp2dreeee67OY/n5+UpISFBKSooSEhKuWOKAQFbjcuuFzPXa\nvv+0JozooWceGK6QkBDTsQAALcB4OastXvWZO3euJKmoqEhJSUmtFQvwGR6PR/+zYKPW7zimUQM6\n60ePjlFYKMUMAAKV8XJ2JUlJSUpISNC4ceMUGxtrOg5gRObSHfpig0MDesXpZ9OTFRHu039sAQDX\nyKf/lnc6nYqJidHMmTP185//XEVFRaYjAa3qk2/2a/HKPereqZ1+8fT1ahMVbjoSAKCF+fTf9AsW\nLNCsWbNks9mUkJCgrKwszZgx44pfM2nSpMtee/bZZzV79uyWigm0iDVbj2jee1sU1z5K/zFjvGLb\nR5mOBABogpdfflmvvPJKg8/3yXLmdDpls1364ObU1NQG3bG5YsUK2e32looGtIqdB05r7hvrFRkR\npp8/bal7p3amIwEAmmj27NmXDRI5HI46B5QkH5jWzMrK0rZt25SVlXX+tenTp0uSZsyYoQULFigr\nK0sLFixQWlqaoZRA6zl0oky//GuOatwe/fTJcRrQK950JABAKzI+cpaamqrU1NRLXluyZMn5j682\njQkEkjPOCv37q2tUeq5Ks6eN1NjBXU1HAgC0MuMjZwC8yitr9Mu/rtWx0+f0yNSBmmr1Nh0JAGAA\n5QzwATUut36XuU57HSWaktxLj0wdaDoSAMAQyhlgmMfj0Z+XbFHezuMaM6iLfvDQCHb/B4AgRjkD\nDHvvy31atvagEnvGas6T4xQexh9LAAhm/BQADFqz9bBe/yRfHWzR+sXTFpvMAgAoZ4Ape4uKNffN\nPEVFhOkXT1vqGNvGdCQAgA+gnAEGnDhTrl+9tlbVNS795PGx6muPMx0JAOAjKGdAKztXUa1fvbZW\np52VevqeoUpO6mY6EgDAh1DOgFbkcrn1+zc2aP9hp25P6aN7bkw0HQkA4GMoZ0Ar+utH+Vq/45hG\nDeisWfcNY8sMAMBlKGdAK/nkm/36aHWBenWL0ZwnxymMLTMAAHXgpwPQCjbvOaFX39+q2PaR+sXT\n16tdmwjTkQAAPopyBrSwo6fO6oXMdQoNkX6WnqyuHdqajgQA8GGUM6AFee/MzFHpuWo988AIJSV2\nNB0JAODjAm478oyMDNlsNiUnJ8uyLNNxEMTcbo9eeitPhUdLddeE63Tb9b1NRwIAGJaTk6Pc3Fw5\nnc56zwm4cpaeni673W46BqA3l+1UTv5RjejfSd+7Z6jpOAAAH2BZlizLksPhUGZmZp3nMK0JtIDV\nGw9p4fLd6taxrZ5/gjszAQANx08MoJntdRTrjws2qk1UmP7tu5Zs7SJNRwIA+BHKGdCMzpRW6Dev\n5ai6xqUfPzZWvbvZTEcCAPgZyhnQTGpcbv0uY51OllToidsH88xMAECTUM6AZvLaR/navv+0bhjR\nQw/d2t90HACAn6KcAc3giw1F+mh1gRK6xui5tFE8MxMA0GSUM+Aa7T9colcWbVbb6HD961PJahMV\ncDvUAABaEeUMuAZl56r0n6/nqqrapX9+ZLR6dm5vOhIAwM9RzoAmcrs9mvvmBh09dU7TJg/Q9UO7\nm44EAAgAlDOgid7+bJc27Dyu0QO76NHbBpmOAwAIEE0uZ4MHD27OHIBfyd1+VO98vktdO7TVjx8f\no7BQbgAAADSPJpczj8dzxeNlZWVNfWvApx0+WaaX3tygyPBQ/Sx9nGLa8gQAAEDzafJtZfVtFbBs\n2TI5nU7ZbDaVlJQoLi5OQ4YM4WHkCAiV1S799vV1OltRo39+ZLT62uNMRwIABJhmvee/rKxMt912\n22Wvl5aWNue3AYyZt2SLDhxx6vbxfXTr2ATTcQAAAahZbwho377ubQRiYmKa89sARqxcX6jPcwuV\n2DNW37t3qOk4AIAA1ahy1tB1ZIsWLZIkrVmzpvGJAB9UeNSpP727RW2jw/XTJ8cpMiLMdCQAQIBq\n0LTmxevIan+/Eo/Hox07djRLwMbKyMiQzWZTcnKyLMsykgGBpaKyRr/LXKfKKpd+mj5O3Tu1Mx0J\nAOCncnJylJubK6fTWe85Vy1n9a0ju5Jp06ZpzZo1Skho/TU56enp3HyAZuPxePSndzer6FiZ7rkx\nUTcM72E6EgDAj1mWJcuy5HA4lJmZWec5Vy1ntevIysrK1L59+/O/X8348eMbGRfwPZ/nFuqLDQ4N\n6BWn6XclmY4DAAgCDV5zNm/ePO3YsUPz5s2TVPc+Zw6HQ3PnztV//dd/yeFwNF9KwID9h0s0b8kW\ntW8ToeefGKeIcB6oAQBoeQ3eSiM2NlYej+eKm88uWLBAs2bNkiR9+umnSkhIYAQNfulcRbV+l7FO\nVTVuzXlynLp2aGs6EgAgSDS4nH3ve9+TdGFbjLo2oR02bNj549OmTdOyZcuaIyPQqjwej/530WYd\nPnlWD9zcT8lJ3UxHAgAEkQbdEHDxGrMrLfJ/5513VFJSoqFDh2rw4MH1PkUA8GXLcwv11aZDGtyn\ng564g2fIAgBaV4NuCFi2bJkcDoeGDBmi0tLSeqc2//u//1uStG3bNv3lL3/RmjVrNHXq1OZNDLSg\nomOlmvf+VrWLDtePHxuj8DDWmQEAWleDpjVrt9IoKipSfHz8VZ8EMH78eI0fP15paWnNFBNoeVXV\nLv3+jfWqrHLph0+OVRfWmQEADKiznM2dO1c33HCDJO80Zu2+YVfbt6y0tPSSRzXVftzQ7TcAkzI+\n2a79h5267fremjCip+k4AIAgVWc5mzVr1vlitWbNmgZv6hoTE3PJFKjT6VRISIiGDBlCOYNPW7f9\nqD5cXSB7l/b63j08NxMAYE6d5WzevHkqLS2VJA0d2rgfVBdPgcbFxfHQc/i8084K/fGdjQoPC9Xz\nT4xVdFSDb2IGAKDZ1flT6I477lBCQkKDi1Vd06AmHt0ENJbb7dFLb22Q82yVZt43TNf1iDUdCQAQ\n5OosZ0OGDGnUmzR1GhQwbcmXe7V5z0mNG9JVd024znQcAAAavgntlVzLNChgyu7CM3rj0x3qYIvS\nc2mj2JcPAOATmlzOLt7rrLHToIBpFZU1+q83N8jt8ehHj45RbPso05EAAJB0DeVs586d5z9u7DQo\nYNprH+Xr8Mmzum9iX43o39l0HAAAzmP7cwSddduP6tM1B9Snu01P8ngmAICPCbg9AzIyMmSz2ZSc\nnCzLskzHgY8pKavU/yzcpPCwUP3o0dGKCA8zHQkAEERycnKUm5srp9NZ7zkBV87S09O5WxR18ng8\nemXRJhWXVuqpu5LYNgMA0Oosy5JlWXI4HMrMzKzzHKY1ETSW5xZq7bajGta3k+6b2Nd0HAAA6kQ5\nQ1A4euqs5n+wVW2jw/XDR0YpNJRtMwAAvolyhoDncrn10lt5Kq906ZkHhqtLfFvTkQAAqBflDAHv\n3S/2aseB05owooduHs16RACAb6OcIaAVHCrRW8t2qoMtWj94aARPAQAA+DzKGQJWdY1bf3g7Ty63\nR8+ljVJM20jTkQAAuCrKGQLWwuW7deCIU1Ot3ho9qIvpOAAANAjlDAFpr6NYC1fsVqe4Nnr6niTT\ncQAAaDDKGQJOdY1Lf3w7T263R/80baTaRkeYjgQAQINRzhBw3vl8tw4eLVXq+D4aNZDpTACAf6Gc\nIaDsKTqjxSv3qEt8Gz111xDTcQAAaDSfKGf5+flXPJaVlaWsrKxWTAR/VF3j0h/e3uidzkwbxXQm\nAMAvGS9n2dnZeu655+o9Pm/ePKWmpqqoqOiKJQ54a9kuFR0r1R0pfTSif2fTcQAAaJJw0wFSUlKU\nkJBQ57GsrCwNGzZMkjRjxozWjAU/s+vgaS35Yo+6dmir6XdxdyYAwH8ZHzm7kq1bt6q4uFj5+fma\nP3++6TjwUdU1Lv33gk1ye6Tn0kapTZTx/+cAAKDJfLqcSVJcXJySkrwjIaw7Q10Wr9ijomOlun18\nHw3r18l0HAAArolPDzHExcWdn/K02WzaunWrUlNTr/g1kyZNuuy1Z599VrNnz26RjDDr4FGnFq7Y\nrY6x0ZrO3ZkAAB/08ssv65VXXmnw+T5ZzpxOp2w2m1JTU8+PljmdzvPrz65kxYoVstvtLR0RPsDl\n9ujlhZtU4/Lo+w8M5+5MAIBPmj179mWDRA6Ho84BJckHpjWzsrK0bdu2S6Ysp0+fLklKSEiQzWZT\nVlaWiouLrzpqhuCy9Jv92nXwjCaM6CFraHfTcQAAaBbGR85SU1MvK11Lliw5/3FaWtr584Bax0+f\nU+bS7WrfJkIz77/6iCoAAP7C+MgZ0Fgej0d/enezKqpc+t69QxUfE206EgAAzYZyBr+zKs+hDTuP\na+SAzrp1bN175AEA4K8oZ/ArJWWVevX9bYqKDNP/e2iEQkJCTEcCAKBZUc7gV/7ywTaVnqvS46mD\n1a1jO9NxAABodpQz+I28ncf1ZZ5D/RPidPeNiabjAADQIihn8AuV1S7935LNCg0N0expIxUWynQm\nACAwUc7gFxYt362jp87pnhsTdV2PWNNxAABoMZQz+LyiY6V694s96hQbrUdvG2Q6DgAALYpyBp/m\n8Xj0f+9uUY3Lo5n3D1ebKOP7JgMA0KIoZ/BpX2xwaOu+k0oe0k3XD+1mOg4AAC2OcgafVXquSq99\n5N3TbNb9w9jTDAAQFChn8FkZn2xXSVmVHpkyUF06tDUdBwCAVkE5g0/asf+0lq09qF7dYnTvxL6m\n4wAA0GooZ/A5NS63/vTuZknSDx4cofAw/jMFAASPgLv1LSMjQzabTcnJybIsy3QcNMFHqwt04IhT\nU5J7KSmxo+k4AAA0m5ycHOXm5srpdNZ7TsCVs/T0dNntdtMx0ESnSsr19mc7FdM2UtPvSjIdBwCA\nZmVZlizLksPhUGZmZp3nMF8En/K3j7arvNKl9DuHyNYu0nQcAABaHeUMPmPbvpNatdH7YPMpyb1M\nxwEAwAjKGXyCy+XWvPe2SpJm3T9MoTzYHAAQpChn8Amfrjlw/iaAgb07mI4DAIAxlDMYV1JWqTey\ndqpddLievGOI6TgAABhFOYNxmUt36Gx5tR5NHaS4mCjTcQAAMIpyBqN2F57R57kH1btbjO5Muc50\nHAAAjKOcwRi326N5722RxyPNemC4wngSAAAAlDOYs3J9oXYXFuumkT01rG8n03EAAPAJlDMYcba8\nWhmf7FBUZJieupsnAQAAUItyBiMWrdit4rJKPTypvzrFtTEdBwAAn0E5Q6s7cvKsPviqQJ3j2+i+\nif1MxwEAwKdQztDq/vZxvmpcbj11Z5KiIsJMxwEAwKdQztCqtu47qTVbj2hwnw6aMLKH6TgAAPgc\nyhlajcvt0V8+2CZJ+t69QxUSwvMzAQD4NsoZWs3KdYUqOFSiW8bYNaBXvOk4AAD4pHDTAZpbRkaG\nbDabkpOTZVmW6Tj4h3MV1fr7p96tM3h+JgAgWOXk5Cg3N1dOp7PecwKunKWnp8tut5uOgW9ZvHKP\nzpRW6tHbBrF1BgAgaFmWJcuy5HA4lJmZWec5TGuixR07fU7vr9qnTrHRuv/mvqbjAADg0yhnaHEZ\nn2xXdY1b6XcOUXRkwA3WAgDQrChnaFE7D57W6k2HNKBXnG4axXQzAABXQzlDi/F4PPrbR/mSpO/e\nPVShoWydAQDA1VDO0GJy8o9q+/7TspK6KSmxo+k4AAD4BcoZWoTL5dbrH29XaGiI0u9k6wwAABqK\ncoYW8VluoQ6dKNNUq7cSusaYjgMAgN+gnKHZlVfW6K1lOxUdGaZHpw40HQcAAL9COUOze//LvSou\nrdT9N/dTvC3adBwAAPwK5QzN6oyzQku+3Ku4mCjdN5ENZwEAaCzKGZrV25/vUkWVS49OHai20RGm\n4wAA4HcoZ2g2juOlWrb2oHp2bq8pVm/TcQAA8EuUMzSbzKU75HZ7lH7nEIWH8Z8WAABNwU9QNIud\nB05rzdYjGtyng64f2s10HAAA/BblDNfM4/Eoc+kOSVL6nUMUEsJjmgAAaCrKGa7Zpt0ntHXfSY0d\n3JXHNAEAcI3CTQdobhkZGbLZbEpOTpZlWabjBDyPx6PMT72jZo+nDjKcBgAA35aTk6Pc3Fw5nc56\nzwm4cpaeni673W46RtBYs/WI9hYVa8KIHuprjzMdBwAAn2ZZlizLksPhUGZmZp3nMK2JJnO5PXoj\na4dCQ0P0+O2DTccBACAgUM7QZKvyilR0rEyTxiaoZ+f2puMAABAQKGdokuoat95ctkvhYaH6Dg83\nBwCg2VDO0CSfrT2g46fP6Y4b+qhLfFvTcQAACBiUMzRaRWWN3lm+W9GRYXr41gGm4wAAEFAoZ2i0\nj7/Zr+LSSt17U1/FxUSZjgMAQEChnKFRysqr9e7KPWrfJkL339zPdBwAAAIO5QyN8uFX+1RWXq0H\nb+2vdm0iTMcBACDgUM7QYGXnqvTBV/sU2z5Sd91wnek4AAAEJMoZGuz9r/bpXEWNHri5v6KjAu7h\nEgAA+ATKGRqk9FyVPvyqQHHto3RHSh/TcQAACFg+Uc7y8/Oves78+fNbIQnq8/6qfSqvrNGDt/Zj\n1AwAgBZkvJxlZ2frueeeu+o52dnZrZQI3+Y8W6WPVu9TXEyUUsf3MR0HAICAZrycpaSkKCEhwXQM\nXMH7q/aqvNKlB2/pr+hIRs0AAGhJxsvZ1eTn5yslJcV0jKBVUlapj78uUHxMlG5nrRkAAC3O58tZ\nSUmJ6QhB7b0vvaNmD93aX1ERYabjAAAQ8Hx6jqopo2aTJk267LVnn31Ws2fPbq5YQaOkrFKffLNf\nHWxRuo21ZgAANMnLL7+sV155pcHn+2Q5czqdstlsKioqUlFRkUpKSlRSUqL8/HwlJSVd8WtXrFgh\nu93eSkkD25Iv9qqiyqUn7xjCqBkAAE00e/bsywaJHA5HnQNKkg9Ma2ZlZWnbtm3Kyso6/9r06dMl\nSampqUpNTZXkLWxoPcWllfoke786xkbrtut7m44DAEDQMD5ydnEBq7VkyZJLPk9LS1NaWlprxgp6\n76/aq8oql6bfOUSRjJoBANBqjI+cwfeUnqvS0uz9io+J0lSLUTMAAFoT5QyX+Wh1gcorXXrgln6M\nmgEA0MooZ7jEuYpqfbi6QDFtI5V6fR/TcQAACDqUM1xiafYBnS2v1r0TE3mGJgAABlDOcF5FVY3e\nX7VX7aLDddcNiabjAAAQlChnOO+ztQdVUlaluyYkql2bCNNxAAAISpQzSJKqa1x694u9io4M0903\nMmoGAIAplDNIkpavK9JpZ4VSx/dRbPso03EAAAhalDOoxuXW4pV7FBEeqvtv7mc6DgAAQY1yBn21\n0aHjp89pqtVbHWzRpuMAABDUKGdBzuX2aOHyPQoLDdEDtzBqBgCAaZSzILdm62EdOlGmW8cmqEt8\nW9NxAAAIepSzIObxeLR45R6FhEgP3drfdBwAACDKWVDbvOeE9jlKlDKsh3p0bm86DgAAkBRwz+fJ\nyMiQzWZTcnKyLMsyHcenvbtyrySx1gwAgFaSk5Oj3NxcOZ3Oes8JuHKWnp4uu91uOobP21tUrE17\nTmh4v04a0CvedBwAAIKCZVmyLEsOh0OZmZl1nsO0ZpB694s9kqQHb2GtGQAAvoRyFoSOnDyr7C2H\nldgjVqMGdjYdBwAAXIRyFoTe+3Kv3B7vWrOQkBDTcQAAwEUoZ0HmTGmFlq8rVNcObTVhRA/TcQAA\nwLdQzoLMR6sLVF3j1v0T+yosjH/9AAD4Gn46B5FzFdVamn1AtnaRmpTcy3QcAABQB8pZEFm29qDO\nllfr7hsTFR0ZcLuoAAAQEChnQaK6xq0Pvtqn6Mgw3XnDdabjAACAelDOgsSqPIdOlVRo6vW9FdM2\n0nQcAABQD8pZEPB4PHp/1V6Fhobo3pv6mo4DAACugHIWBDbtPqGDR0s1YUQPdYlvazoOAAC4AspZ\nEHh/1T5J0n0TGTUDAMDXUc4C3MEjTuXtOq6kxI7qn8ADzgEA8HWUswD3wVeMmgEA4E8oZwHsTGmF\nvtjgUPdO7TRuSDfTcQAAQANQzgLY0m8OqMbl1r039VVYKA84BwDAHwTcNvEZGRmy2WxKTk6WZVmm\n4xhTWe3S0uz9at8mQpPGJpiOAwAAJOXk5Cg3N1dOp7PecwKunKWnp8tut5uOYdwX64vkPFulhyf1\nV3RUwP1rBgDAL1mWJcuy5HA4lJmZWec5TGsGILfbow++2qfwsBAe1QQAgJ+hnAWgvF3H5ThepptG\n2dUxto3pOAAAoBEoZwHovS/3ShKPagIAwA9RzgJMwaESbdl7UsP7dVJiz1jTcQAAQCNRzgIMm84C\nAODfKGcBpLi0Ul9tPKSendtpzKCupuMAAIAmoJwFkGVrvZvO3jUhUaFsOgsAgF+inAWIGpdbS7P3\nq01UuG5l01kAAPwW5SxAZG85rNPOSk1J7qW20RGm4wAAgCainAWIj1YXKCREunMCm84CAODPKGcB\nYHfhGe08eEZjBnVVj07tTccBAADXgHIWAD7+ukCSdPeNiYaTAACAa0U583NnnBVavemQ7F3aa9SA\nzqbjAACAa0Q583NZaw+qxuXRXRMSFRLC9hkAAPg7ypkfq65x69Ps/WobzfYZAAAECsqZH/tmy2Gd\nKa3UlOTeahMVbjoOACJWF58AABJZSURBVABoBgH3Ez0jI0M2m03JycmyLMt0nBb18T+2z7iL7TMA\nAPALOTk5ys3NldPprPecgCtn6enpstvtpmO0uF0HT2tX4RlZSd3UrWM703EAAEADWJYly7LkcDiU\nmZlZ5zlMa/qpj7/eL0m6ewLbZwAAEEgoZ37otLNCX28+pISuMRrev5PpOAAAoBlRzvzQZzm122dc\nx/YZAAAEGMqZn3G53Mpac0BtosJ18+j/v727CYrqTvc4/hPfozTkJmZGwxnNy5RVdDPJJEprOy83\n0RlbKlW3yk1blQWyADciG90lLExW0Y0vm4QF0aosesNuTJvEO6+cCBrHGbrJNffORDiY+JJE+mB8\niULfhfZRAigip/+nm++naqqkuzn98AzV+fE/5zn/0r+2DgCAmYZwVmS6e8/r2+x1vbrK0mML5pou\nBwAATDPCWZE50nlWkrQptsJoHQAAwB+EsyIycHFIp//3kiLPPaHlPw2ZLgcAAPiAcFZEPvz0rCSp\nLsZNZwEAKFWEsyJx/cYtHevu1+Pl87UmstR0OQAAwCeEsyLxl9Pn9P31W/r9muWaO4f/2wAAKFX8\nV74I5HI5/aHzS5WVzVJ8zQrT5QAAAB8FYm/NTCajcDg87nPJZFKS1N/fr127dhWyrMD4ov+y/n0u\nq7U1S/Vk5ULT5QAAAB8ZXzmzbVstLS0TPheLxZRIJOQ4jmzbLnB1wXDEPitJquP2GQAAlDzj4SwW\ni8myrHGfuzeQWZYlx3EKWVogZK/c0F9Pn9PTSxbpF88vMV0OAADwWSBOa04kkUh4/+7t7VVdXZ3B\nasw4dqJfN2+NaFPsGZWVsY8mAAClzvjK2WRkMhlVV1dPeF1aqRoZyemIfVbz5s7W+lXjry4CAIDS\nEuiVszzbtic9DLB+/foxj23fvl3Nzc3TXZbvTp25qAvfXdXvan+mxY/NM10OAACYggMHDujgwYOT\nfn0gw5nrugqFbm9PlEwm1djYKOnugMD9HDt2TFVVVb7XWAhH7C8lSXXr2BEAAIBi1dzcPGaRaGBg\nYNwFJSkApzVTqZTS6bRSqZT32NatWyXdDmN79+7Vhg0btHr1akMVmnHhu6s6+fkFrfzZ43q+qtJ0\nOQAAoECMr5zF43HF4/FRj3V0dEi6Pcl54sQJE2UZ91FXn3I5aRO3zwAAYEYxvnKGsW4Nj+iT7j4t\nWjhXv3rxadPlAACAAiKcBdCJ3gv6zr2hV16q0vy5s02XAwAACohwFkBHj5+VJG1cu8JoHQAAoPAI\nZwFz8fJVnTpzUSuXP64VS0OmywEAAAVGOAuYj7v6lctJ8TXLTZcCAAAMIJwFyPDwiD7u7tNjC+bo\nVy8wCAAAwExEOAuQz/7nor7NXtdvX6rSgvnG73ICAAAMIJwFSOrOIEB8zQqjdQAAAHMIZwHxzeA1\nffb5Bf3cqtSzT1eYLgcAABhCOAuIj7v7NZKTNrJqBgDAjEY4C4DhkZw+6urTwvmz9ZtfMggAAMBM\nRjgLgL+fuahvBq/pN7+s0kIGAQAAmNEIZwFwlEEAAABwB+HMsG+z19Tde0HPVVXoeavSdDkAAMCw\nkjuHdujQIYVCIdXW1ioajZou54E+OdGvkZEcgwAAAMwAXV1d6u7uluu6E76m5MJZfX29qqqqTJcx\nKSMjOX10vE8L5s3WbxkEAACg5EWjUUWjUQ0MDOjw4cPjvobTmgad/uKSLl6+pl+/+LQeWzDXdDkA\nACAACGcGHe06K0mKr11htA4AABAchDNDLrvX1ZU+r2eWhfRzBgEAAMAdhDNDPjnRr+E7gwCzZs0y\nXQ4AAAgIwpkBI3d2BJg3d7b+86XiGF4AAACFQTgz4J//d0nnv72qX7+4TIsWMggAAADuIpwZcPR4\nnyR2BAAAAGMRzgpscOiGjqe/1vKflmvl8sdNlwMAAAKGcFZg/32yX7eGGQQAAADjI5wVUC6XU+p4\nn+bNKdMrLzMIAAAAxiKcFVDPv77R1998r3UvLNPix+aZLgcAAAQQ4ayA8oMAbHIOAAAmQjgrkOyV\nG7L/+bWsnyxW9TP/YbocAAAQUISzAvnjZ45uDY8wCAAAAO6LcFYAuVxOqU/7NHdOmV552TJdDgAA\nCDDCWQFk/v2tzl26oljNMoUWMQgAAAAmNsd0AdPt0KFDCoVCqq2tVTQaNV2OJOlo151BgLXLDVcC\nAABM6urqUnd3t1zXnfA1JRfO6uvrVVUVnHuIDV39QZ3/+EpPL1mkyLNPmC4HAAAYFI1GFY1GNTAw\noMOHD4/7Gk5r+uyPJx3dvMUgAAAAmBzCmY/yOwLMmV2mV1cxCAAAAB6McOajz89+J+fCkNbWLFXF\n4vmmywEAAEWAcOajuzsCMAgAAAAmh3DmkytXf9DfTp/T0icXqea5J02XAwAAigThzCd/OjWgH26N\naGN0ucrKGAQAAACTQzjzQS6X09HjfZpdNkuvrmYQAAAATB7hzAdn+i7r7Neu1tQs1ePlC0yXAwAA\nigjhzAep42clSXEGAQAAwEMinE2zK9du6q+nv9LSJxbpF88vMV0OAAAoMoSzafanzxz9cHNYv1/D\nIAAAAHh4hLNpdO8gwHoGAQAAwBQQzqbRmX4GAQAAwKMhnE2j1KdnJTEIAAAApo5wNk0YBAAAANOB\ncDZN/swgAAAAmAZzTBcw3Q4dOqRQKKTa2lpFo9GCvGcul1OKQQAAAPAAXV1d6u7uluu6E76m5MJZ\nfX29qqqqCvqe+UGAdS8sYxAAAABMKBqNKhqNamBgQIcPHx73NZzWnAZHP+2TxCAAAAB4dISzR3Tl\n2k395fQ5BgEAAMC0IJw9IgYBAADAdCKcPQIGAQAAwHQjnD2CM313dgSIsCMAAACYHoSzR/CHzi8l\nSXXrVpgtBAAAlAzC2RQNDt3Q3/7xlayfLFbNc0+aLgcAAJQIwtkUfdzdp1vDI6qLPaNZsxgEAAAA\n04NwNgXDIzl9+OlZLZg3W6+8zCAAAACYPoSzKTjZe16XLl/TKy9bWrRwrulyAABACSGcTcHdQYBn\nDFcCAABKTSDCWSaTmfC5VCol27bV1tZWwIomdu7SFf39i0sKP/uEViwNmS4HAACUGOPhzLZttbS0\njPtcPrTFYjGFQqH7hrhC+dA+K0mqi60wWgcAAChNxsNZLBaTZY1/Uf2RI0cUCt1enbIsS7ZtF7K0\nMa7/cEufnOhXZfl8ra1ZZrQWAABQmoyHs/sZGhpSRUWF9/Xg4KDBaqQ/nzqn76/d1Mbocs2dE+jW\nAQCAIjXHdAHTbf369WMe2759u5qbmx/52Ol/faM5s2cpvnbFIx8LAADMDAcOHNDBgwcn/fpAh7Py\n8nJls1lJkuu6qqysfOD3HDt2TFVVVb7Us/W1av3Xb5/Tk5ULfTk+AAAoPc3NzWMWiQYGBsZdUJIC\nelrTdV1JUl1dnRzHkSQ5jqNYLGayLD1RsVDPVz04IAIAAEyV8XCWSqWUTqeVSqW8x7Zu3SpJCofD\nkm5PdIZCIe9rAACAUmX8tGY8Hlc8Hh/1WEdHh/fvRCJR6JIAAACMMb5yBgAAgLsIZwAAAAFCOAMA\nAAgQwhkAAECAEM4AAAAChHAGAAAQIIQzAACAACGcAQAABAjhDAAAIEAIZwAAAAFCOAMAAAgQwhkA\nAECAEM4AAAACZI7pAqbboUOHFAqFVFtbq2g0arocAAAAT1dXl7q7u+W67oSvKbmVs/r6ejU3N/sW\nzA4cOODLcTE++l049Lqw6Hfh0OvCot/3F41G1dzcrPr6+glfU3LhzG8HDx40XcKMQr8Lh14XFv0u\nHHpdWPT70RHOAAAAAqQkw1lXV5fpEqbEz7qL9dh+8qtuej1WsfakGPtdrP0oxl5LxduTYux3sfZj\nKscuyXDW3d1tuoQp8bPuYj22n/yqm16PVaw9KcZ+F2s/irHXUvH2pBj7Xaz9mMqxS2Zac3h4WJJ0\n/vx5ua6rgYEB397Lr2P7WXexHlsqvn7T67GKtSfF2O9i7Ucx9loq3p4UY7+LtR8THfv8+fOS7uaX\ne83K5XI5X6opsJMnT+r11183XQYAAMCkffDBB1q1atWox0omnF2/fl3pdFpLlizR7NmzTZcDAAAw\noeHhYV26dEmRSEQLFiwY9VzJhDMAAIBSUJIDAQAAAMWKcAYAABAghDMAAIAAIZwB8LS1tZkuAUAR\nymQyo75OpVKybZvPlCkinD0EftkKJ5lMKplMas+ePaZLmTFs25Zt26bLKHmZTEapVEqpVMp0KTNC\n/nM7mUyaLqVk2batlpYW7+t8UIvFYgqFQmOCGx6McDZJ/LIVjm3bisViSiQSchyHwICS8u677yoe\nj8txHD5HfJbJZGRZlmKxmCzLot8+yfc378iRIwqFQpIky7L4DJ8Cwtkk8ctWOPcGMsuy5DiO4YpK\nXyaTUSwWM11GyUulUqqpqZEkNTY2KhwOG66o9O3du1fS7c8V+l0YQ0NDqqio8L4eHBw0WE1xIpxN\nEr9shZNIJJRIJCRJvb29ikQihisqfdls1nQJM0JPT48GBweVyWS4PKIAwuGwLMvS6tWrR31+A0FH\nOENgZTIZVVdX89euz1g1K6zKykrvd5rrzvzluq7Ky8vV1NSkN998k1X4AikvL/f+4HNdV5WVlYYr\nKj4ls/G53/hlKzzbtrVr1y7TZZQ8x3HkOI6y2ayy2awymQyB2CeVlZXetTmhUEg9PT2Kx+OGqypd\nyWRS27ZtUygUkmVZSqVSamxsNF1Wyaurq1M6nZZ0+/OFP/4eHitnk1RXV+f91cUvm/+SyaT3Icr1\nff6Kx+NeQHBd13A1pS0/CCDd7nX++jP4Lx6Pe9cNY3qlUiml02lvJTj/x51t2wqFQvyxNwXsrfkQ\nksmkd4F6/pooTL/8WHZFRYWy2az27dtHGEbJSCaTqqioUE9PDyvDBdDW1ibLspTNZvncRtEgnAEA\nAAQIpzUBAAAChHAGAAAQIIQzAACAACGcAQAABAjhDICvksmkGhoavP/5fSNQ13XV0NDwUN/jOM5D\nf8+D3mflypVqbW1Va2urGhoaJr3x9v2OO9FzP67ftm01NDSotbV1Uu8JIFi4CS0A37S1tamnp0ft\n7e2mS7kvy7J8qXH37t3evxsaGhSJRHy559OP629ra/O+TiaT3EICKDKsnAHwzXvvvae3337bdBmB\nUF1dXbA9TO99n8mu2AEIDlbOAPgik8nIsqwJ78re2trqneLcuXOnstmsd8NQy7KUSCS0detWRSIR\npdNpbdq0SY7jKJ1Oq6mpydvV4MfHyW+PNNH7hMNh77Hy8nJZlqXXXntN77zzjrfatGfPHvX29kqS\nGhsbFYlE9MYbb3jH3L9//0P1Ir/zQv5myve+f/5Y+dOS99Zv2/aYnmSzWbW2to7qg+u6amlpUXt7\nu9ra2uQ4jnbs2KGnnnrK+/eWLVu4mTNQJAhnAAouv9vG7t275bquF8ISicSoraQcx1FHR4ccx9GG\nDRt05swZOY6j1tZWxePxcY/z/vvv3/d93nrrLUlSe3u7GhoatG3btlG1pVIpDQ0NjTnNmQ9Rra2t\nsm1bkUjkgT/njh07NDQ0pGw2q46ODq+m8vJytbe3K5lMKpVKKZvNqrq6Wrt27fJ+vnwt9+tJvg/3\namxsVDKZ9Oo9derUQ4dJAGZxWhOAL8LhsBzHGXe/zs7OTm8VJ7+y1tTUJNu2tXnzZm+PvnwAqqio\n8K7Vym/rNdFx7n2/8Z4vKyvT4OCgpNun/368smfb9rgrTMlkUnv27FE6nZ70HqT79+9Xe3u7XNf1\n9ojNZDJesOrs7PQ2m1+3bp338+Xt3Lnzvj0p1GlSAIVFOAPgm6amplGnA6Xb4WndunVeWHFdV67r\nqqqqSrt371ZHR4f27t07qeOPd5x7w9Z4zy9evFhDQ0PasWOHmpqaxhwzFot535P/vmQyKdd1tWvX\nLq1du/bhmiBp37593mpYOBxWTU2Ndu/erf379yuRSCgcDquzs1OSRk2zhkKhh+4JgOLHaU0Avsmf\nYrv3Ng87d+5UIpFQa2urNm/eLOl2eMmf4pM06enC8Y4zmefT6bQsy1IymVQ2m9WmTZu874nH496t\nKPI/QyQSUUtLi3fNW01NzUP1IRwOa+PGjWpra1NjY6N3e43y8nJt2bJFiUTCu9XIvdecTaUnPxaJ\nRNTQ0KB4PM7UJlAk2PgcwIySv04rFot516HlrwcDgCAgnAGYUTKZzKjbS1iWpcbGRoMVAcBohDMA\nAIAAYSAAAAAgQAhnAAAAAUI4AwAACBDCGQAAQIAQzgAAAAKEcAYAABAg/w8pa47F0bVDQgAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1a1df4ef28>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (10,10))\n",
    "plt.plot(redshifts, ratio)\n",
    "gnz11 = np.where(redshifts == 11.09)\n",
    "ratio_gal = vel_now(z)/(cos.recession_velocity(z)*10**(-3))\n",
    "plt.plot(redshifts[gnz11], ratio[gnz11], 'o', c='red', ms=5)\n",
    "plt.text(redshifts[gnz11]-0.7, ratio[gnz11].value-0.08, 'GN-z11', color='red', fontsize = 15)\n",
    "plt.plot(z, ratio_gal, 'o', c='red', ms=5)\n",
    "plt.text(z+0.15, ratio_gal.value, 'z=4.25 galaxy', color='red', fontsize = 15)\n",
    "plt.rc('axes', labelsize=20) \n",
    "plt.xlabel('Cosmological Redshift')\n",
    "plt.ylabel(r\"$\\frac{v_{_{(z)}}}{v_{_{(Dc)}}}$\")\n",
    "plt.savefig('velocity_methods.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4 (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "term1 = (2/3)*(omega_lambda)**(-1/2)\n",
    "term2 = np.log((1+np.sqrt(omega_lambda))/np.sqrt(omega_matter))\n",
    "corr_fac = term1*term2\n",
    "t0 = (1/h)*9.78*corr_fac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "astropy_tl = model.lookback_time(z)\n",
    "astropy_age = t0 - astropy_tl.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "redshifts = np.linspace(0, 4.25, num=100)\n",
    "my_tl = cos.t_lookback(H0, omega_matter, omega_lambda, z)\n",
    "my_age = t0 - my_tl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Astropy age (in Gyr): 1.488174748352579\n",
      "Age from Gaussian quad (in Gyr): [ 0.76713638 13.86925375]\n"
     ]
    }
   ],
   "source": [
    "print('Astropy age (in Gyr):', astropy_age)\n",
    "print('Age from Gaussian quad (in Gyr):', my_age)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It appears there's something wrong with my code, because the uncertainty is large, and increasing the number of sampling points has no effect on the result. Don't have time to fix it, so I'll use the result from Mathematica, which does something similar to Gaussian quadrature.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "integral = 0.88636"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "my_tl = (100/H0)*9.78*integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Age from Mathematica (in Gyr): 1.4855383270249014\n"
     ]
    }
   ],
   "source": [
    "print('Age from Mathematica (in Gyr):', t0-my_tl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4 (c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D: 0.24977927467836344\n"
     ]
    }
   ],
   "source": [
    "from scipy.integrate import quad\n",
    "func = lambda z: (1+z) / model.H(z).value**3\n",
    "# find the constant of proportionality \n",
    "const = 1/quad(func, 0, np.inf)[0]\n",
    "integral = quad(func, z, np.inf)[0]\n",
    "print('D:', model.H(z)/model.H0 * const * integral)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "# repeat for Einstein-de-Sitter \n",
    "omega_m = 1\n",
    "omega_l = 0\n",
    "h = 0.5\n",
    "H0 = 100*h\n",
    "model1 = FlatLambdaCDM(H0, omega_m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D: 0.19047576679988082\n"
     ]
    }
   ],
   "source": [
    "func = lambda z: (1+z) / model1.H(z).value**3\n",
    "# find the constant of proportionality \n",
    "const = 1/quad(func, 0, np.inf)[0]\n",
    "integral = quad(func, z, np.inf)[0]\n",
    "print('D:', model1.H(z)/model1.H0 * const * integral)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4 (d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Comoving distance: 7587.113751902942 Mpc\n"
     ]
    }
   ],
   "source": [
    "print(\"Comoving distance:\", model.comoving_distance(z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4 (e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Light travel distance: 3796.056080219697 Mpc\n"
     ]
    }
   ],
   "source": [
    "print(\"Light travel distance:\", model.lookback_distance(z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4 (f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Luminosity distance: 39832.347197490446 Mpc\n"
     ]
    }
   ],
   "source": [
    "print(\"Luminosity distance:\", model.luminosity_distance(z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4 (g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Angular diameter distance: 1445.164524171989 Mpc\n",
      "Physical size (in kpc): 7.0063487463796035\n"
     ]
    }
   ],
   "source": [
    "da = model.angular_diameter_distance(z)\n",
    "print(\"Angular diameter distance:\", da)\n",
    "rad_to_arcsec = 206265\n",
    "actual_size = da.value*(10**3/rad_to_arcsec)\n",
    "print(\"Physical size (in kpc):\", actual_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4 (h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Angular diameter distance: 1786.7937786626903 Mpc\n",
      "Angular size at z = 1.5 (in arcseconds): 0.8088031990203196\n"
     ]
    }
   ],
   "source": [
    "da = model.angular_diameter_distance(1.5)\n",
    "print(\"Angular diameter distance:\", da)\n",
    "print(\"Angular size at z = 1.5 (in arcseconds):\", actual_size*rad_to_arcsec/(da.value*10**(3)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}