Lecture 1 Answers.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cosmological Parameter Estimation Lecture 1 Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Running These Exercises\n",
    "\n",
    "### MyBinder\n",
    "You can run these exercises by visting https://mybinder.org/v2/gh/joezuntz/michigan/master but your work will not automatically be saved at the end \n",
    "of each session.\n",
    "\n",
    "**To save your work:**\n",
    "- In the menu under the Jupyter logo at the top left of the page, select File > Download as > Notebook.  This will save your work to your machine.\n",
    "\n",
    "**To restore your work:**\n",
    "- Visit https://mybinder.org/v2/gh/joezuntz/michigan/master again\n",
    "- Click \"Upload\"\n",
    "- Find the file you downloaded and select it\n",
    "- Click the upload button that appears, and overwrite the existing file.\n",
    "\n",
    "### Local Jupyter\n",
    "\n",
    "You can also clone this repository from `git clone https://github.com/joezuntz/michigan` and run it locally if you already know how to do that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# Preamble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can run each cell by highlighting it and pressing Shift-Enter.\n",
    "\n",
    "We will first do some setup that lets us make plots inline in the notebook - run this next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import camb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "# Lecture 1: Essential Data Handling & Modelling\n",
    "\n",
    "## Exercise 1: Astropy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many public cosmology libraries that have been released to help calculate theory predictions or explore data.\n",
    "\n",
    "One of these is Astropy, a python library for astronomy research.  It contains tools for many common astronomy tasks, including basic (background) cosmology.\n",
    "\n",
    "In this first exercise we will explore how to create a cosmological model using Astropy, and use it to predict some observable values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.cosmology import LambdaCDM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you execute the cell below some documentation will pop up at the bottom.\n",
    "Read it and learn how to create a new Lambda CDM cosmology model with specified parameters.\n",
    "\n",
    "HINT: Some parameters are optional."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "LambdaCDM?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, create a test model called `cosmo` with Hubble parameter 72.0, matter density 0.3, dark energy density 0.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here\n",
    "cosmo = LambdaCDM(H0=72., Om0=0.3, Ode0=0.7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your `cosmo` object has functions attached to it (\"methods\") that calculate different cosmological quantities.\n",
    "Run the cell below and scroll down in the cell to find a listing of them.\n",
    "\n",
    "Find the method `distmod` and read about its meaning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on LambdaCDM in module astropy.cosmology.core object:\n",
      "\n",
      "class LambdaCDM(FLRW)\n",
      " |  LambdaCDM(H0, Om0, Ode0, Tcmb0=0, Neff=3.04, m_nu=<Quantity 0. eV>, Ob0=None, name=None)\n",
      " |  \n",
      " |  FLRW cosmology with a cosmological constant and curvature.\n",
      " |  \n",
      " |  This has no additional attributes beyond those of FLRW.\n",
      " |  \n",
      " |  Parameters\n",
      " |  ----------\n",
      " |  \n",
      " |  H0 : float or `~astropy.units.Quantity`\n",
      " |      Hubble constant at z = 0.  If a float, must be in [km/sec/Mpc]\n",
      " |  \n",
      " |  Om0 : float\n",
      " |      Omega matter: density of non-relativistic matter in units of the\n",
      " |      critical density at z=0.\n",
      " |  \n",
      " |  Ode0 : float\n",
      " |      Omega dark energy: density of the cosmological constant in units of\n",
      " |      the critical density at z=0.\n",
      " |  \n",
      " |  Tcmb0 : float or scalar `~astropy.units.Quantity`, optional\n",
      " |      Temperature of the CMB z=0. If a float, must be in [K].\n",
      " |      Default: 0 [K]. Setting this to zero will turn off both photons\n",
      " |      and neutrinos (even massive ones).\n",
      " |  \n",
      " |  Neff : float, optional\n",
      " |      Effective number of Neutrino species. Default 3.04.\n",
      " |  \n",
      " |  m_nu : `~astropy.units.Quantity`, optional\n",
      " |      Mass of each neutrino species. If this is a scalar Quantity, then all\n",
      " |      neutrino species are assumed to have that mass. Otherwise, the mass of\n",
      " |      each species. The actual number of neutrino species (and hence the\n",
      " |      number of elements of m_nu if it is not scalar) must be the floor of\n",
      " |      Neff. Typically this means you should provide three neutrino masses\n",
      " |      unless you are considering something like a sterile neutrino.\n",
      " |  \n",
      " |  Ob0 : float or None, optional\n",
      " |      Omega baryons: density of baryonic matter in units of the critical\n",
      " |      density at z=0.  If this is set to None (the default), any\n",
      " |      computation that requires its value will raise an exception.\n",
      " |  \n",
      " |  name : str, optional\n",
      " |      Name for this cosmological object.\n",
      " |  \n",
      " |  Examples\n",
      " |  --------\n",
      " |  >>> from astropy.cosmology import LambdaCDM\n",
      " |  >>> cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7)\n",
      " |  \n",
      " |  The comoving distance in Mpc at redshift z:\n",
      " |  \n",
      " |  >>> z = 0.5\n",
      " |  >>> dc = cosmo.comoving_distance(z)\n",
      " |  \n",
      " |  Method resolution order:\n",
      " |      LambdaCDM\n",
      " |      FLRW\n",
      " |      Cosmology\n",
      " |      builtins.object\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __init__(self, H0, Om0, Ode0, Tcmb0=0, Neff=3.04, m_nu=<Quantity 0. eV>, Ob0=None, name=None)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |  \n",
      " |  de_density_scale(self, z)\n",
      " |      Evaluates the redshift dependence of the dark energy density.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      I : ndarray, or float if input scalar\n",
      " |        The scaling of the energy density of dark energy with redshift.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      The scaling factor, I, is defined by :math:`\\rho(z) = \\rho_0 I`,\n",
      " |      and in this case is given by :math:`I = 1`.\n",
      " |  \n",
      " |  efunc(self, z)\n",
      " |      Function used to calculate H(z), the Hubble parameter.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      E : ndarray, or float if input scalar\n",
      " |        The redshift scaling of the Hubble constant.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      The return value, E, is defined such that :math:`H(z) = H_0 E`.\n",
      " |  \n",
      " |  inv_efunc(self, z)\n",
      " |      Function used to calculate :math:`\\frac{1}{H_z}`.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      E : ndarray, or float if input scalar\n",
      " |        The inverse redshift scaling of the Hubble constant.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      The return value, E, is defined such that :math:`H_z = H_0 /\n",
      " |      E`.\n",
      " |  \n",
      " |  w(self, z)\n",
      " |      Returns dark energy equation of state at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      w : ndarray, or float if input scalar\n",
      " |        The dark energy equation of state\n",
      " |      \n",
      " |      Notes\n",
      " |      ------\n",
      " |      The dark energy equation of state is defined as\n",
      " |      :math:`w(z) = P(z)/\\rho(z)`, where :math:`P(z)` is the\n",
      " |      pressure at redshift z and :math:`\\rho(z)` is the density\n",
      " |      at redshift z, both in units where c=1.  Here this is\n",
      " |      :math:`w(z) = -1`.\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data and other attributes defined here:\n",
      " |  \n",
      " |  __abstractmethods__ = frozenset()\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Methods inherited from FLRW:\n",
      " |  \n",
      " |  H(self, z)\n",
      " |      Hubble parameter (km/s/Mpc) at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      H : `~astropy.units.Quantity`\n",
      " |        Hubble parameter at each input redshift.\n",
      " |  \n",
      " |  Ob(self, z)\n",
      " |      Return the density parameter for baryonic matter at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      Ob : ndarray, or float if input scalar\n",
      " |        The density of baryonic matter relative to the critical density at\n",
      " |        each redshift.\n",
      " |      \n",
      " |      Raises\n",
      " |      ------\n",
      " |      ValueError\n",
      " |        If Ob0 is None.\n",
      " |  \n",
      " |  Ode(self, z)\n",
      " |      Return the density parameter for dark energy at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      Ode : ndarray, or float if input scalar\n",
      " |        The density of non-relativistic matter relative to the critical\n",
      " |        density at each redshift.\n",
      " |  \n",
      " |  Odm(self, z)\n",
      " |      Return the density parameter for dark matter at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      Odm : ndarray, or float if input scalar\n",
      " |        The density of non-relativistic dark matter relative to the critical\n",
      " |        density at each redshift.\n",
      " |      \n",
      " |      Raises\n",
      " |      ------\n",
      " |      ValueError\n",
      " |        If Ob0 is None.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      This does not include neutrinos, even if non-relativistic\n",
      " |      at the redshift of interest.\n",
      " |  \n",
      " |  Ogamma(self, z)\n",
      " |      Return the density parameter for photons at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      Ogamma : ndarray, or float if input scalar\n",
      " |        The energy density of photons relative to the critical\n",
      " |        density at each redshift.\n",
      " |  \n",
      " |  Ok(self, z)\n",
      " |      Return the equivalent density parameter for curvature\n",
      " |      at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      Ok : ndarray, or float if input scalar\n",
      " |        The equivalent density parameter for curvature at each redshift.\n",
      " |  \n",
      " |  Om(self, z)\n",
      " |      Return the density parameter for non-relativistic matter\n",
      " |      at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      Om : ndarray, or float if input scalar\n",
      " |        The density of non-relativistic matter relative to the critical\n",
      " |        density at each redshift.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      This does not include neutrinos, even if non-relativistic\n",
      " |      at the redshift of interest; see `Onu`.\n",
      " |  \n",
      " |  Onu(self, z)\n",
      " |      Return the density parameter for neutrinos at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      Onu : ndarray, or float if input scalar\n",
      " |        The energy density of neutrinos relative to the critical\n",
      " |        density at each redshift.  Note that this includes their\n",
      " |        kinetic energy (if they have mass), so it is not equal to\n",
      " |        the commonly used :math:`\\sum \\frac{m_{\\nu}}{94 eV}`,\n",
      " |        which does not include kinetic energy.\n",
      " |  \n",
      " |  Tcmb(self, z)\n",
      " |      Return the CMB temperature at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      Tcmb : `~astropy.units.Quantity`\n",
      " |        The temperature of the CMB in K.\n",
      " |  \n",
      " |  Tnu(self, z)\n",
      " |      Return the neutrino temperature at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      Tnu : `~astropy.units.Quantity`\n",
      " |        The temperature of the cosmic neutrino background in K.\n",
      " |  \n",
      " |  __repr__(self)\n",
      " |      Return repr(self).\n",
      " |  \n",
      " |  abs_distance_integrand(self, z)\n",
      " |      Integrand of the absorption distance.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : float or array\n",
      " |        Input redshift.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      X : float or array\n",
      " |        The integrand for the absorption distance\n",
      " |      \n",
      " |      References\n",
      " |      ----------\n",
      " |      See Hogg 1999 section 11.\n",
      " |  \n",
      " |  absorption_distance(self, z)\n",
      " |      Absorption distance at redshift ``z``.\n",
      " |      \n",
      " |      This is used to calculate the number of objects with some\n",
      " |      cross section of absorption and number density intersecting a\n",
      " |      sightline per unit redshift path.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      d : float or ndarray\n",
      " |        Absorption distance (dimensionless) at each input redshift.\n",
      " |      \n",
      " |      References\n",
      " |      ----------\n",
      " |      Hogg 1999 Section 11. (astro-ph/9905116)\n",
      " |      Bahcall, John N. and Peebles, P.J.E. 1969, ApJ, 156L, 7B\n",
      " |  \n",
      " |  age(self, z)\n",
      " |      Age of the universe in Gyr at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      t : `~astropy.units.Quantity`\n",
      " |        The age of the universe in Gyr at each input redshift.\n",
      " |      \n",
      " |      See Also\n",
      " |      --------\n",
      " |      z_at_value : Find the redshift corresponding to an age.\n",
      " |  \n",
      " |  angular_diameter_distance(self, z)\n",
      " |      Angular diameter distance in Mpc at a given redshift.\n",
      " |      \n",
      " |      This gives the proper (sometimes called 'physical') transverse\n",
      " |      distance corresponding to an angle of 1 radian for an object\n",
      " |      at redshift ``z``.\n",
      " |      \n",
      " |      Weinberg, 1972, pp 421-424; Weedman, 1986, pp 65-67; Peebles,\n",
      " |      1993, pp 325-327.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      d : `~astropy.units.Quantity`\n",
      " |        Angular diameter distance in Mpc at each input redshift.\n",
      " |  \n",
      " |  angular_diameter_distance_z1z2(self, z1, z2)\n",
      " |      Angular diameter distance between objects at 2 redshifts.\n",
      " |      Useful for gravitational lensing.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z1, z2 : array_like, shape (N,)\n",
      " |        Input redshifts. z2 must be large than z1.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      d : `~astropy.units.Quantity`, shape (N,) or single if input scalar\n",
      " |        The angular diameter distance between each input redshift\n",
      " |        pair.\n",
      " |  \n",
      " |  arcsec_per_kpc_comoving(self, z)\n",
      " |      Angular separation in arcsec corresponding to a comoving kpc\n",
      " |      at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      theta : `~astropy.units.Quantity`\n",
      " |        The angular separation in arcsec corresponding to a comoving kpc\n",
      " |        at each input redshift.\n",
      " |  \n",
      " |  arcsec_per_kpc_proper(self, z)\n",
      " |      Angular separation in arcsec corresponding to a proper kpc at\n",
      " |      redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      theta : `~astropy.units.Quantity`\n",
      " |        The angular separation in arcsec corresponding to a proper kpc\n",
      " |        at each input redshift.\n",
      " |  \n",
      " |  clone(self, **kwargs)\n",
      " |      Returns a copy of this object, potentially with some changes.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      newcos : Subclass of FLRW\n",
      " |      A new instance of this class with the specified changes.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      This assumes that the values of all constructor arguments\n",
      " |      are available as properties, which is true of all the provided\n",
      " |      subclasses but may not be true of user-provided ones.  You can't\n",
      " |      change the type of class, so this can't be used to change between\n",
      " |      flat and non-flat.  If no modifications are requested, then\n",
      " |      a reference to this object is returned.\n",
      " |      \n",
      " |      Examples\n",
      " |      --------\n",
      " |      To make a copy of the Planck13 cosmology with a different Omega_m\n",
      " |      and a new name:\n",
      " |      \n",
      " |      >>> from astropy.cosmology import Planck13\n",
      " |      >>> newcos = Planck13.clone(name=\"Modified Planck 2013\", Om0=0.35)\n",
      " |  \n",
      " |  comoving_distance(self, z)\n",
      " |      Comoving line-of-sight distance in Mpc at a given\n",
      " |      redshift.\n",
      " |      \n",
      " |      The comoving distance along the line-of-sight between two\n",
      " |      objects remains constant with time for objects in the Hubble\n",
      " |      flow.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      d : `~astropy.units.Quantity`\n",
      " |        Comoving distance in Mpc to each input redshift.\n",
      " |  \n",
      " |  comoving_transverse_distance(self, z)\n",
      " |      Comoving transverse distance in Mpc at a given redshift.\n",
      " |      \n",
      " |      This value is the transverse comoving distance at redshift ``z``\n",
      " |      corresponding to an angular separation of 1 radian. This is\n",
      " |      the same as the comoving distance if omega_k is zero (as in\n",
      " |      the current concordance lambda CDM model).\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      d : `~astropy.units.Quantity`\n",
      " |        Comoving transverse distance in Mpc at each input redshift.\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      This quantity also called the 'proper motion distance' in some\n",
      " |      texts.\n",
      " |  \n",
      " |  comoving_volume(self, z)\n",
      " |      Comoving volume in cubic Mpc at redshift ``z``.\n",
      " |      \n",
      " |      This is the volume of the universe encompassed by redshifts less\n",
      " |      than ``z``. For the case of omega_k = 0 it is a sphere of radius\n",
      " |      `comoving_distance` but it is less intuitive\n",
      " |      if omega_k is not 0.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      V : `~astropy.units.Quantity`\n",
      " |        Comoving volume in :math:`Mpc^3` at each input redshift.\n",
      " |  \n",
      " |  critical_density(self, z)\n",
      " |      Critical density in grams per cubic cm at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      rho : `~astropy.units.Quantity`\n",
      " |        Critical density in g/cm^3 at each input redshift.\n",
      " |  \n",
      " |  differential_comoving_volume(self, z)\n",
      " |      Differential comoving volume at redshift z.\n",
      " |      \n",
      " |      Useful for calculating the effective comoving volume.\n",
      " |      For example, allows for integration over a comoving volume\n",
      " |      that has a sensitivity function that changes with redshift.\n",
      " |      The total comoving volume is given by integrating\n",
      " |      differential_comoving_volume to redshift z\n",
      " |      and multiplying by a solid angle.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      dV : `~astropy.units.Quantity`\n",
      " |        Differential comoving volume per redshift per steradian at\n",
      " |        each input redshift.\n",
      " |  \n",
      " |  distmod(self, z)\n",
      " |      Distance modulus at redshift ``z``.\n",
      " |      \n",
      " |      The distance modulus is defined as the (apparent magnitude -\n",
      " |      absolute magnitude) for an object at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      distmod : `~astropy.units.Quantity`\n",
      " |        Distance modulus at each input redshift, in magnitudes\n",
      " |      \n",
      " |      See Also\n",
      " |      --------\n",
      " |      z_at_value : Find the redshift corresponding to a distance modulus.\n",
      " |  \n",
      " |  kpc_comoving_per_arcmin(self, z)\n",
      " |      Separation in transverse comoving kpc corresponding to an\n",
      " |      arcminute at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      d : `~astropy.units.Quantity`\n",
      " |        The distance in comoving kpc corresponding to an arcmin at each\n",
      " |        input redshift.\n",
      " |  \n",
      " |  kpc_proper_per_arcmin(self, z)\n",
      " |      Separation in transverse proper kpc corresponding to an\n",
      " |      arcminute at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      d : `~astropy.units.Quantity`\n",
      " |        The distance in proper kpc corresponding to an arcmin at each\n",
      " |        input redshift.\n",
      " |  \n",
      " |  lookback_distance(self, z)\n",
      " |      The lookback distance is the light travel time distance to a given\n",
      " |      redshift. It is simply c * lookback_time.  It may be used to calculate\n",
      " |      the proper distance between two redshifts, e.g. for the mean free path\n",
      " |      to ionizing radiation.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      d : `~astropy.units.Quantity`\n",
      " |        Lookback distance in Mpc\n",
      " |  \n",
      " |  lookback_time(self, z)\n",
      " |      Lookback time in Gyr to redshift ``z``.\n",
      " |      \n",
      " |      The lookback time is the difference between the age of the\n",
      " |      Universe now and the age at redshift ``z``.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      t : `~astropy.units.Quantity`\n",
      " |        Lookback time in Gyr to each input redshift.\n",
      " |      \n",
      " |      See Also\n",
      " |      --------\n",
      " |      z_at_value : Find the redshift corresponding to a lookback time.\n",
      " |  \n",
      " |  lookback_time_integrand(self, z)\n",
      " |      Integrand of the lookback time.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : float or array_like\n",
      " |        Input redshift.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      I : float or array\n",
      " |        The integrand for the lookback time\n",
      " |      \n",
      " |      References\n",
      " |      ----------\n",
      " |      Eqn 30 from Hogg 1999.\n",
      " |  \n",
      " |  luminosity_distance(self, z)\n",
      " |      Luminosity distance in Mpc at redshift ``z``.\n",
      " |      \n",
      " |      This is the distance to use when converting between the\n",
      " |      bolometric flux from an object at redshift ``z`` and its\n",
      " |      bolometric luminosity.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.  Must be 1D or scalar.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      d : `~astropy.units.Quantity`\n",
      " |        Luminosity distance in Mpc at each input redshift.\n",
      " |      \n",
      " |      See Also\n",
      " |      --------\n",
      " |      z_at_value : Find the redshift corresponding to a luminosity distance.\n",
      " |      \n",
      " |      References\n",
      " |      ----------\n",
      " |      Weinberg, 1972, pp 420-424; Weedman, 1986, pp 60-62.\n",
      " |  \n",
      " |  nu_relative_density(self, z)\n",
      " |      Neutrino density function relative to the energy density in\n",
      " |      photons.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array like\n",
      " |         Redshift\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |       f : ndarray, or float if z is scalar\n",
      " |         The neutrino density scaling factor relative to the density\n",
      " |         in photons at each redshift\n",
      " |      \n",
      " |      Notes\n",
      " |      -----\n",
      " |      The density in neutrinos is given by\n",
      " |      \n",
      " |      .. math::\n",
      " |      \n",
      " |        \\rho_{\\nu} \\left(a\\right) = 0.2271 \\, N_{eff} \\,\n",
      " |        f\\left(m_{\\nu} a / T_{\\nu 0} \\right) \\,\n",
      " |        \\rho_{\\gamma} \\left( a \\right)\n",
      " |      \n",
      " |      where\n",
      " |      \n",
      " |      .. math::\n",
      " |      \n",
      " |        f \\left(y\\right) = \\frac{120}{7 \\pi^4}\n",
      " |        \\int_0^{\\infty} \\, dx \\frac{x^2 \\sqrt{x^2 + y^2}}\n",
      " |        {e^x + 1}\n",
      " |      \n",
      " |      assuming that all neutrino species have the same mass.\n",
      " |      If they have different masses, a similar term is calculated\n",
      " |      for each one. Note that f has the asymptotic behavior :math:`f(0) = 1`.\n",
      " |      This method returns :math:`0.2271 f` using an\n",
      " |      analytical fitting formula given in Komatsu et al. 2011, ApJS 192, 18.\n",
      " |  \n",
      " |  scale_factor(self, z)\n",
      " |      Scale factor at redshift ``z``.\n",
      " |      \n",
      " |      The scale factor is defined as :math:`a = 1 / (1 + z)`.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      z : array_like\n",
      " |        Input redshifts.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      a : ndarray, or float if input scalar\n",
      " |        Scale factor at each input redshift.\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors inherited from FLRW:\n",
      " |  \n",
      " |  H0\n",
      " |      Return the Hubble constant as an `~astropy.units.Quantity` at z=0\n",
      " |  \n",
      " |  Neff\n",
      " |      Number of effective neutrino species\n",
      " |  \n",
      " |  Ob0\n",
      " |      Omega baryon; baryonic matter density/critical density at z=0\n",
      " |  \n",
      " |  Ode0\n",
      " |      Omega dark energy; dark energy density/critical density at z=0\n",
      " |  \n",
      " |  Odm0\n",
      " |      Omega dark matter; dark matter density/critical density at z=0\n",
      " |  \n",
      " |  Ogamma0\n",
      " |      Omega gamma; the density/critical density of photons at z=0\n",
      " |  \n",
      " |  Ok0\n",
      " |      Omega curvature; the effective curvature density/critical density\n",
      " |      at z=0\n",
      " |  \n",
      " |  Om0\n",
      " |      Omega matter; matter density/critical density at z=0\n",
      " |  \n",
      " |  Onu0\n",
      " |      Omega nu; the density/critical density of neutrinos at z=0\n",
      " |  \n",
      " |  Tcmb0\n",
      " |      Temperature of the CMB as `~astropy.units.Quantity` at z=0\n",
      " |  \n",
      " |  Tnu0\n",
      " |      Temperature of the neutrino background as `~astropy.units.Quantity` at z=0\n",
      " |  \n",
      " |  critical_density0\n",
      " |      Critical density as `~astropy.units.Quantity` at z=0\n",
      " |  \n",
      " |  h\n",
      " |      Dimensionless Hubble constant: h = H_0 / 100 [km/sec/Mpc]\n",
      " |  \n",
      " |  has_massive_nu\n",
      " |      Does this cosmology have at least one massive neutrino species?\n",
      " |  \n",
      " |  hubble_distance\n",
      " |      Hubble distance as `~astropy.units.Quantity`\n",
      " |  \n",
      " |  hubble_time\n",
      " |      Hubble time as `~astropy.units.Quantity`\n",
      " |  \n",
      " |  m_nu\n",
      " |      Mass of neutrino species\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors inherited from Cosmology:\n",
      " |  \n",
      " |  __dict__\n",
      " |      dictionary for instance variables (if defined)\n",
      " |  \n",
      " |  __weakref__\n",
      " |      list of weak references to the object (if defined)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(cosmo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, use `cosmo.distmod` to work out the distance modulus at the single redshift value $z=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$44.039065 \\; \\mathrm{mag}$"
      ],
      "text/plain": [
       "<Quantity 44.03906537 mag>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Write your code here\n",
    "cosmo.distmod(1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, create a numpy array of redshifts `z` = 0.1, 0.2, 0.3, ... 1.6 using numpy, and compute the distance modulus to it.\n",
    "Make a plot of distance modulus versus redshift in this model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'mu')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEKCAYAAAACS67iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl0leW5/vHvnYQQhjBkYgohQSYRBTEgglqKtU5otVqHOg/Q6tJ62qqt5/R09gy/9rRWTycrWq0D1tZa61DbgqggEBImgSgEEghjAiFhDhnu80e2/qIFswM7efdwfdbKWtk7786+CMmVJ89+3uc1d0dERGJfUtABREQkMlToIiJxQoUuIhInVOgiInFChS4iEidU6CIicUKFLiISJ1ToIiJxQoUuIhInUjrzybKysjw/P78zn1JEJOaVlJTsdPfsto7r1ELPz8+nuLi4M59SRCTmmdnGcI7TlIuISJxQoYuIxAkVuohInFChi4jECRW6iEicUKGLiMQJFbqISJxQoYuIdJBd++r566rt/PDlNRxqaOrw5+vUE4tEROLZ5t0HWFJRQ1H5bpZU1FBWtQ+A1JQkPj8+l9EDe3Xo86vQRUSOgbtTVrWPoooaisprWFJew9a6QwCkp6VQOKQvnx8/iIn5GZyc25uuKckdnkmFLiIShsamZlZv3UNReQ1FFTUUV9Sw+0ADANnpXZmYn8HM/L5MLMhkZP90kpOs0zOq0EVEjuDg4SaWVe5mSWj6ZOmm3Rw43DIPPiSzO+ec2I+JBRlMzM9gSGZ3zDq/wD9OhS4iAhw43EhxxW4WbtjF4g27eHdLHQ1NjhmM7JfOF07LZUJBBhPyM+jXKy3ouEcUdqGbWTJQDGxx9+lm9lvgU0Bd6JCb3H155COKiETeoYYmlm2qZeGGXSxcv5PllbU0NDkpScYpub259cyhTCzoy2lDMujdrUvQccPSnhH63UAp0Ppl2nvd/Q+RjSQiEnkNTc2s3FzLwvW7eGf9Lko27qa+sZkkg5MH9eaWMwuYfEIWhUP60qNrbE5ehJXazHKBi4AHgK91aCIRkQhoanZWb637sMCXVNR8OAd+4oBeXHv6ECafkMmEgtgZgbcl3F9DDwL3Aekfu/8BM/s2MAf4prvXRzKciEi4mpud93fs/bDAF5fvYu+hRgCG5fTk8vG5TD4hk9OHZpLRIzXgtB2jzUI3s+lAlbuXmNnUVh+6H9gOpAKPAN8Avn+Ex88EZgLk5eVFILKISItNuw7w1rpq3lm/k0UbaqjZfxhoWYUy/ZQBTBqayRlDM8mJ0hcxIy2cEfoU4BIzuxBIA3qZ2VPufl3o4/Vm9jhwz5Ee7O6P0FL4FBYWegQyi0iC2nuogYXrd/H2up28ta6ajbsOADCwdxqfHpnDGSdkcsYJmQzq0y3gpMFos9Dd/X5aRuOERuj3uPt1ZjbA3bdZy+LLS4FVHZpURBJOU7Ozaksdb62t5u11O1m6aTeNzU731GTOGJrJLVMKOGt4FgVZPaJiHXjQjuel3KfNLBswYDnw5chEEpFEtq3uIG+vbRmBLyjb+eHZmGMG9WLm2UM5a3g2pw3pS2qK9hb8uHYVurvPA+aF3p/WAXlEJMEcPNzE4vLQNMraataFNrTKSe/KtFH9OHtEFmcOyyKzZ9eAk0a/2FxsKSIxy915b/veD6dRiipqONzYTGpKEqcXZHBl4WDOGpHFyH7pmkZpJxW6iHS4/fWNzC/byZzSHbzxfjXVe1tWOI/sl84Nk4Zw9ohsJhZkkNal43ckjGcqdBHpEJU1B5j7XhVz3qti0fpdHG5qJj0thbNHZDN1RDZnDc+mf+/EWE7YWVToIhIRjU3NLKusZU5pFXPf28HaHS1z4UOze3Dj5CFMG9WPwvy+dEnWi5kdRYUuIses7kADb66rZm7pDuatrab2QAMpScbpQzO4akIe00blUJDVI+iYCUOFLiJhc3fWV+9n7ns7mFNaRfHG3TQ1Oxk9Upk2KodzRvXjrBFZ9EqLj71RYo0KXUQ+0eHGZpZU1Hw4lVIROjtzVP90vvypoUwb1Y9xg/sEcoUe+SgVuoj8k/31jcx7v5rXVm1j3vvV7KtvJDUlicknZHLrWUOZNionYU+vj2YqdBEBoO5gA3NKd/Daqu28tbaa+sZmMnukMv2UAZxzYj+mDMuke6oqI5rpf0ckge3aV8/f1uzgr6u28876nTQ0Of17pXHNxDzOH9OfCfkZmkqJISp0kQSzve4Qr6/ezmurtlFUXkOzQ15Gd26ZUsD5Y/ozNrcPSSrxmKRCF0kAlTUHeG3VNl5btZ1lm2oBGJ7Tkzs/PYzzxvRn9IBeOs0+DqjQReJUWdVeXnt3O39dvZ3VW/cALTsW3nveSM47qT/DcnoGnFAiTYUuEic+2PTq1XdbRuJloV0LTxvSl3+78ETOH9OfwRndA04pHUmFLhLjKmsO8NKKrfx5+RbW7thHksGkoZnccMYQzjupP/0S5PJrokIXiUk1+w/zysqtvLh8KyUbdwMwIb8vP7x0DBeM6a+9wxOUCl0kRhw43Mjf1+zgxWVbeHvdThqbnZH90rnv/JFcfMpATaeICl0kmjU0NTN/3U5eXL6Fv63ewcGGJgb2TuO2s4Zy6akDGdW/V9ARJYqo0EWijLtTsnE3f16+lVfe3UbN/sP06d6Fy8YP4tJxgygc0lfrxOWIwi50M0sGioEt7j691f0PAbe4u9ZAiRyHtTv28uKyLby0Yiubdx8krUsSnzmxH5eOG8TZI7J1UWRpU3tG6HcDpcCHf+OZWSHQN9KhRBLF1tqDoRUqWyndtofkJGPKsCy+du4IPntSf3p21R/REr6wvlvMLBe4CHgA+FrovmTgR8AXgcs6KqBIvDnc2Mw/Sncwe0klb6+rxh3GDe7Ddy8ezUWnDCQ7XStU5NiE++v/QeA+IL3VfXcCL7n7Np0yLNK2sqq9PLekkheWbmHX/sMM6J3GXZ8exuWn5TIkU1f1kePXZqGb2XSgyt1LzGxq6L6BwBeAqWE8fiYwEyAvL+94sorEnAOHG3ll5TaeW1JJ8cbdpCQZnzmxH1dNHMzZw7O1k6FElLn7Jx9g9p/A9UAjkEbLHHp96O1Q6LA8YIO7D/ukz1VYWOjFxcXHm1kkqrk7q7bsYfaSTby0fCt76xsZmtWDqyYM5vPjczWlIu1mZiXuXtjWcW2O0N39fuD+0CedCtzTepVL6P59bZW5SLyrO9DAn1dsYXZRJWu27SGtSxIXnjyAqyfkMSG/r3YzlA6nl9BFjoO7s7i8htlFm3ht1XbqG5sZM6gXP7h0DJeMHUjvbrpYsnSedhW6u88D5h3hfq1Bl4RStfcQfyzZwnNLNlGx6wDpaSlcWTiYqyYMZsyg3kHHkwSlEbpImJqanTfXVjG7qJI571XR1OxMLMjgK+cM54IxA+iWmhx0RElwKnSRNuw51MDvl1TyxMIKKmsOktUzldvOKuDKwsGckK0/TiV6qNBFjqJ8536eeKeC54sr2X+4iQn5fbn/ghM5d3Q/uiTrNHyJPip0kVbcnQVlu3h8QTlz368iJcm4eOxAbplSoLlxiXoqdBHgUEMTf1q2hccXlLN2xz6yeqbylWnDuXZSHjnpuuKPxAYVuiS0bXUH+d3CjTxbtIndBxoYPaAXP7riFC4eO5C0LnqRU2KLCl0S0tJNu3l8QQWvvbuNZnfOHd2PW6YUMLEgQycAScxSoUvCaGhq5tV3t/H4ggqWV9aS3jWFmybnc+PkfF2+TeKCCl3iXs3+wzxbtInfLdzI9j2HKMjqwfcuOYnLT8vVfuMSV/TdLHGrrGovj75dzp+WbaG+sZmzhmfxH58fw9QRObqEm8QlFbrEndJte/jfuWW8umobqclJfH58LjdPyWdEv/S2HywSw1ToEjdWbanj4bnreH31Dnp2TeGOqSdw65lDyeiRGnQ0kU6hQpeYt6KylofnruMfpVWkp6XwlXOGc8uUfPp0V5FLYlGhS8wq2bibh+euY9771fTu1oWvnzuCGybna8taSVgqdIk5ReU1PDx3HW+v20nf7l247/yRXD9pCOlpKnJJbCp0iQnuzsINu3hozjoWbaghq2cq/3rhKK49fQg9tPRQBFChS5Rzd+aX7eShOetYUrGbnPSu/Pv00XxxYp72Hxf5GBW6RCV3Z97aah6as45lm2oZ0DuN73/uJK4sHKw9VkSOQoUuUcXdmVNaxUNz17Fycx2D+nTjgcvGcMVpuXRNUZGLfJKwC93MkoFiYIu7TzezWUAhYMBa4CZ339cxMSXeuTt/W7ODn/1jHWu27WFwRjf++/KTuezUXFJTdDEJkXC0Z4R+N1AK9Ard/qq77wEws58AdwL/Fdl4kghWbq7lhy+XUlRRQ0FWD378hbF8btxAXRVIpJ3CKnQzywUuAh4AvgbQqswN6AZ4B2WUOLWt7iA/+uv7vLBsC1k9U/mPy07mysJcUlTkIsck3BH6g8B9wEc2wzCzx4ELgTXA1yMbTeLV/vpGfv3WBh55az3NDrdPPYE7pp6gdeQix6nNQjez6UCVu5eY2dTWH3P3m0Nz6w8DVwGPH+HxM4GZAHl5eZHILDGqudn5w9LN/Pj196naW8/0UwbwjfNHaS9ykQgx90+eKTGz/wSuBxqBNFrm0F9w9+taHXM2cJ+7T/+kz1VYWOjFxcXHHVpiz8L1u/jhK2tYvXUP4wb34d+nj+a0IX2DjiUSE8ysxN0L2zquzRG6u98P3B/6pFOBe4DrzWyYu5eF5tAvAd47vsgSj8p37uc/Xy3lb2t2MKhPN3529TguGTtQl3kT6QDHug7dgCfMrFfo/RXA7RFLJTGv7kADP5uzjicXVtA1JYl7zxvJrWcW6KQgkQ7UrkJ393nAvNDNKZEOI7GvoamZpxZt5Gdz1rHnYANXTRjMV88dQU56WtDRROKezhSViPjgDM//eLWUDTv3M2VYJt+6aDQnDujV9oNFJCJU6HLcVm+t44FXSnln/S6GZvdg1o2FTBuVo3lykU6mQpdjVrXnEP/zt7X8vqSS3t268N2LR3PtpCE6w1MkICp0abfmZufJhRX86PX3OdzUzK1TCrhr2nB6d9eJQSJBUqFLu5Tv3M83/rCSoooazh6RzfcuOYmCrB5BxxIRVOgSpqZm5/EF5fz4b+/TJTmJ/3fFKXzhtFzNk4tEERW6tKmsah/3/WEFSzfVcs6oHB647GT699YyRJFoo0KXo2psaubR+eX85O9r6dYlmZ9eNZZLxw3SqFwkSqnQ5YjW7tjLvc+vYMXmOs47qR8/uHSMTg4SiXIqdPmIhqZmfv3meh6aU0bPtBQevuZUpp8yQKNykRigQpcPlW7bwz3Pr2D11j1cdPIAvve5k8jq2TXoWCISJhW6cLixmZ+/UcbP3yijT/cu/PLa8Vxw8oCgY4lIO6nQE9yqLXXc8/wK3tu+l8+NG8h3Lj6JjB6pQccSkWOgQk9Q9Y1NPDynjF++uZ6MHqn85oZCzh3dL+hYInIcVOgJaEVlLff+YQVrd+zj8vG5fHv6aJ22LxIHVOgJ5FBDEz/9x1p+89YGctLTePymCXx6VE7QsUQkQlToCWLl5lr+5bnlbKjez1WFg/m36SfSK02jcpF4okJPAL8vruRbL64iq0cqT94ykbNHZAcdSUQ6gAo9jh1ubOYHL6/hd4s2MmVYJg9fM14rWETiWNiFbmbJQDGwxd2nm9nTQCHQABQBX3L3ho6JKe1VtfcQdzy1lOKNu/nS2UO597yRpOjCEyJxrT0/4XcDpa1uPw2MAk4GugG3RTCXHIelm3Zz8cPzWbW1joeuOZX7LzxRZS6SAML6KTezXOAi4NEP7nP3Vz2ElhF6bsdElPaYXbSJq3+9iNSUJF64fQqXjB0YdCQR6SThTrk8CNwHpH/8A2bWBbielhG8BKS+sYnv/WUNzyzexFnDs3j4mlPp013z5SKJpM1CN7PpQJW7l5jZ1CMc8gvgLXd/+yiPnwnMBMjLyzuOqHI0O/Yc4vanSli6qZbbp57APZ8dSXKSdkcUSTThjNCnAJeY2YVAGtDLzJ5y9+vM7DtANvCloz3Y3R8BHgEoLCz0CGSWVko21vDlp5ayv76Rn39xPBedok21RBJVm3Po7n6/u+e6ez5wNTA3VOa3AecB17h7cwfnlI9xd55evJGrH1lE99Rk/nTHFJW5SII7nnXovwI2AgtDFz94wd2/H5FU8onqG5v4zp9XM3tJJVNHZvOzq07VXiwi0r5Cd/d5wLzQ+zopKQDb6w7x5adKWF5Zy52fHsZXzx2h+XIRAXSmaEwpKq/hjqeXcvBwI7+6bjznj9EUi4j8fyr0GODu/G7RRr7/lzXkZXTn2RmnM7zfP60gFZEEp0KPcocamvj3F1fxfMlmzhmVw0+vHqddEkXkiFToUWxr7UG+/FQJKzfXcfc5w7n7nOEkab5cRI5ChR6lFm/YxR1PL6W+sVmXhxORsKjQo9Bba6u57clicvt245HrCxmW0zPoSCISA1ToUebtddXMeLKYE7J78sxtp9NX+5eLSJi0p2oUWVC2k9ueKKYgqwdPq8xFpJ1U6FHinbKd3PrEEgqyevDMjEm6spCItJsKPQosXL+LW55YwpCMlpG5ylxEjoUKPWCLNuzilt8uYXDf7jw943Qye3YNOpKIxCgVeoAWb9jFzY8vYVDfbjwzYxJZKnMROQ4q9IAsqajh5t8uYWCfNJ6ZcTrZ6SpzETk+KvQAFFfUcNNjRfTvncazMyaRk54WdCQRiQMq9E5WsrGGGx8rol+vNGbPmEROL5W5iESGCr0TLd20mxsfW0JOrzSenakyF5HIUqF3kmWbdnPjrCKyeqby7IxJ9FOZi0iEqdA7wfLKWm6YVURGz1SenTmJ/r1V5iISeSr0DrZycy3Xz1pM3x4tI/MBvbsFHUlE4lTYhW5myWa2zMxeDt2+08zKzMzNLKvjIsaudzfXcd2ji+nTvQvPzpzEwD4qcxHpOO0Zod8NlLa6vQD4DLAxoonixKotdVw3azG9unXh2RmTGKQyF5EOFlahm1kucBHw6Af3ufsyd6/ooFwxbdWWOq59dDE9u6bw7IxJ5PbtHnQkEUkA4Y7QHwTuA5o7MEtcWLN1D9fNWkyP1GRmz5zE4AyVuYh0jjYL3cymA1XuXnIsT2BmM82s2MyKq6urj+VTxIzSbXu49tFFdOuSzOyZZ6jMRaRThTNCnwJcYmYVwGxgmpk9Fe4TuPsj7l7o7oXZ2dnHGDP6vbd9D9c+upiuKS0j87xMlbmIdK42C93d73f3XHfPB64G5rr7dR2eLIZsqT3Itb9ZTJdkY/bMSQzJ7BF0JBFJQMe8Dt3MvmJmm4FcYKWZPdrWY+JRY1Mz/zJ7GYcamnhmxiTys1TmIhKMdl0k2t3nAfNC7z8EPBT5SLHlf98oY0nFbn561VhOyO4ZdBwRSWA6U/Q4FJXX8NCcdXz+1EFcdmpu0HFEJMGp0I9R3YEG/mX2MgZndOf7l44JOo6ISPumXKSFu/PNF1ZStbeeP94+mZ5d9WUUkeBphH4MZi+p5LVV2/n6Z0cydnCfoOOIiAAq9HYrq9rL9/6ymjOHZfGls4cGHUdE5EMq9HY41NDEXc8up3tqCj+5cixJSRZ0JBGRD2nytx3+67X3KN22h1k3FurycSISdTRCD9Pc93bw23cquGlyPuec2C/oOCIi/0SFHoaqPYe45/mVnDigF9+8YFTQcUREjkiF3obmZudrv1/BgcONPHzNONK6JAcdSUTkiFTobXjk7Q3ML9vJdy4+iWE56UHHERE5KhX6J1hRWcuPX3+fC8b05+oJg4OOIyLyiVToR7GvvpGvzF5GTnpX/uvzp2CmJYoiEt20bPEovv3iKiprDjB75hn07t4l6DgiIm3SCP0I/rRsMy8s28Jd04YzsSAj6DgiImFRoX/Mxl37+dafVlE4pC93TRsWdBwRkbCp0FtpaGrmK7OXk5xkPHj1OFKS9eURkdihOfRWfvL3tayorOUX144nt68u8iwisUVD0JAFZTv51ZvruXrCYC48eUDQcURE2i3sQjezZDNbZmYvh24XmNliMyszs+fMLLXjYnasXfvq+epzyxma1YNvXzw66DgiIsekPSP0u4HSVrf/G/ipuw8DdgO3RjJYZ3F3vvHHldQeaODha8bTPVWzUCISm8IqdDPLBS4CHg3dNmAa8IfQIU8Al3ZEwI725MKN/KO0im9eMIrRA3sFHUdE5JiFO0J/ELgPaA7dzgRq3b0xdHszMCjC2Tpc6bY9PPBqKdNG5XDzlPyg44iIHJc2C93MpgNV7l5yLE9gZjPNrNjMiqurq4/lU3SIg4ebuOvZZfTu1oUfXaFT+0Uk9oUzQp8CXGJmFcBsWqZafgb0MbMPJpxzgS1HerC7P+Luhe5emJ2dHYHIkfGDV9ZQVrWPn1w5lsyeXYOOIyJy3NosdHe/391z3T0fuBqY6+7XAm8AV4QOuxH4c4eljLA33q/imcWb+NKnhnLW8Oj5JSMicjyOZx36N4CvmVkZLXPqsyITqeP94o0ycvt24+vnjgw6iohIxLRrjZ67zwPmhd7fAEyMfKSOtXJzLUsqdvOti04kNUXnVYlI/Ei4Rps1v5yeXVO4ShesEJE4k1CFvr3uEK+s3MaVhYNJT9Me5yISXxKq0J9cWEGzu9aci0hcSphCP3i4iWeKNvHZ0f0ZnKGdFEUk/iRMof9x6WZqDzRwy5kFQUcREekQCVHozc3OYwvKOXlQbybk9w06johIh0iIQn9zbTUbqvdz65kFOsVfROJWQhT6YwvK6derqy5cISJxLe4L/f3te3l73U5uOCNfJxKJSFyL+4Z7bH45aV2SuPb0vKCjiIh0qLgu9J376vnT8i1cPj6XPt1j9gp5IiJhietCf3rRJg43NnPzFC1VFJH4F7eFXt/YxO8WbWTqyGyG5fQMOo6ISIeL20L/y4pt7NxXz606kUhEEkRcFrq7M2t+OSP7pXPmsKyg44iIdIq4LPSFG3ZRum0Pt5yZrxOJRCRhxGWhPza/nIweqXxu3KCgo4iIdJq4K/TynfuZ814V152eR1qX5KDjiIh0mrgr9McXlNMlKYnrzhgSdBQRkU7VZqGbWZqZFZnZCjNbbWbfC90/zcyWmtkqM3vCzNp1fdKOUHeggeeLN3Px2IHkpKcFHUdEpFOFM0KvB6a5+1hgHHC+mU0GngCudvcxwEbgxo6LGZ7ZSzZxsKFJSxVFJCG1WejeYl/oZpfQWxNw2N3Xhu7/O3B5x0QMT2NTM0+8U8EZQzMZPbBXkFFERAIR1hy6mSWb2XKgipbyLgJSzKwwdMgVwOCOiRie11ZtZ2vdIY3ORSRhhVXo7t7k7uOAXGAicBJwNfBTMysC9tIyav8nZjbTzIrNrLi6ujpCsf/ZrPnl5Gd2Z9qonA57DhGRaNauVS7uXgu8AZzv7gvd/Sx3nwi8Baw9ymMecfdCdy/Mzs4+/sRHULJxN8sra7l5SgFJSTqRSEQSUzirXLLNrE/o/W7AucB7ZpYTuq8r8A3gVx0Z9JM8Nr+cXmkpXHFablARREQCF85SwwHAE2aWTMsvgN+7+8tm9iMzmx6675fuPrcjgx7N5t0HeG3VNmacNZQeXQNfOSkiEpg2G9DdVwKnHuH+e4F7OyJUezy5cCNmxo2T84OOIiISqJg+U3R/fSPPFm3igjH9GdinW9BxREQCFdOF/nxxJXsPNWqpoogIMVzoTc3O4+9UcGpeH07N6xt0HBGRwMVsoc8p3cHGXQc0OhcRCYnZQn9sQTmD+nTj/JP6Bx1FRCQqxGShr95ax6INNdw4eQgpyTH5TxARibiYbMNZ88vpnprMVRPygo4iIhI1Yq7Qq/Yc4i8rtnJl4WB6d+sSdBwRkagRc4X+u0UbaWx2btKJRCIiHxFThX6ooYmnF2/inFH9yM/qEXQcEZGoElOF/uKyLdTsP6yliiIiRxAzhe7uPLagnNEDejFpaEbQcUREok7MFPrb63aydsc+bj2zADPteS4i8nExU+iz5peTnd6Vi8cODDqKiEhUiolCL6vay5trq7lh0hBSU2IisohIp4uJdpw1v4LUlCS+eLpOJBIROZqYKPQhmd259cwCMnt2DTqKiEjUiolrtn35UycEHUFEJOrFxAhdRETa1mahm1mamRWZ2QozW21m3wvdf46ZLTWz5WY238yGdXxcERE5mnBG6PXANHcfC4wDzjezScAvgWvdfRzwDPCtjospIiJtaXMO3d0d2Be62SX05qG3XqH7ewNbOyKgiIiEJ6wXRc0sGSgBhgE/d/fFZnYb8KqZHQT2AJM6LqaIiLQlrBdF3b0pNLWSC0w0szHAV4EL3T0XeBz4yZEea2YzzazYzIqrq6sjlVtERD6mXatc3L0WeAO4ABjr7otDH3oOmHyUxzzi7oXuXpidnX1cYUVE5OjCWeWSbWZ9Qu93A84FSoHeZjYidNgH94mISEDCmUMfADwRmkdPAn7v7i+b2Qzgj2bWDOwGbmnrE5WUlOw0s43HlThysoCdQYdoQ7RnjPZ8oIyREO35IPozHm++IeEcZC2LWBKPmRW7e2HQOT5JtGeM9nygjJEQ7fkg+jN2Vj6dKSoiEidU6CIicSKRC/2RoAOEIdozRns+UMZIiPZ8EP0ZOyVfws6hi4jEm0QeoYuIxJW4L3QzO9/M3jezMjP75hE+/jUzW2NmK81sjpmFtTyoMzO2Ou5yM3Mz69RX88PJZ2ZXhr6Oq83smc7MF05GM8szszfMbFno//rCTs73mJlVmdmqo3zczOyhUP6VZjY+yvJdG8r1rpm9Y2ZjOzNfOBlbHTfBzBrN7IrOytbqudvMaGZTQ7vUrjazNyMawN3j9g1IBtYDQ4FUYAUw+mPHfBroHnr/duC5aMsYOi4deAtYBBRGUz5gOLAM6Bu6nRNtX0Na5jBvD70/Gqjo5IxnA+OBVUf5+IXAa4DRsi/S4ijLN7nV/+8FnZ0vnIytvhfmAq8CV0RbRqAPsAbIC92O6M9KvI/QJwJl7r7B3Q8Ds4HPtT7A3d9w9wOhm4to2a8mqjKG/AD4b+BQZ4YjvHwzaNm0bTeAu1dFYcZAdwd197eAmk+RDl+pAAADqElEQVQ45HPAk95iEdDHzAZ0Trq287n7Ox/8/xLMz0k4X0OAu4A/Ap39PQiElfGLwAvuvil0fERzxnuhDwIqW93eHLrvaG6lZZTUmdrMGPrze7C7v9KZwULC+RqOAEaY2QIzW2Rm53dauhbhZPwucJ2ZbaZl9HZX50QLW3u/V4MUxM9Jm8xsEHAZLddqiFYjgL5mNs/MSszshkh+8pi4pmhnMLPrgELgU0Fnac3MkmjZyfKmgKN8khRapl2m0jJye8vMTvaWzdyixTXAb939f8zsDOB3ZjbG3ZuDDhZLzOzTtBT6mUFnOYIHgW+4e7OZBZ3laFKA04BzgG7AQjNb5O5rI/XJ49kWYHCr27mh+z7CzD4D/BvwKXev76RsH2grYzowBpgX+ibtD7xkZpe4e3EU5IOW0eRid28Ays1sLS0Fv6QT8kF4GW8Fzgdw94VmlkbL/hqB/Gl+BGF9rwbJzE4BHgUucPddQec5gkJgdujnJAu40Mwa3f3FYGN9xGZgl7vvB/ab2VvAWCAihR7vUy5LgOFmVmBmqcDVwEutDzCzU4FfA5cEMPfbZkZ3r3P3LHfPd/d8WuYvO6vM28wX8iIto3PMLIuWPys3dFK+cDNuomVUhJmdCKQB0bRB/0vADaHVLpOAOnffFnSoD5hZHvACcH2kRpOR5u4FrX5O/gDcEWVlDvBn4EwzSzGz7sDpRHCn2rgeobt7o5ndCbxOy6vfj7n7ajP7PlDs7i8BPwJ6As+HfrNvcvdLoixjYMLM9zrwWTNbAzQB93bmCC7MjF8HfmNmX6XlBdKbPLTMoDOY2bO0/NLLCs3jf4eWyzni7r+iZV7/QqAMOADc3FnZwsz3bSAT+EXo56TRO3kzrDAyBq6tjO5eamZ/BVYCzcCj7v6JyzDb9fyd+D0tIiIdKN6nXEREEoYKXUQkTqjQRUTihApdRCROqNBFROKECl1EJE6o0EVE4oQKXRKGmeWb2Xtm9lszW2tmT5vZZ0Kbiq0zs4lm9l0zu6fVY1aZWX5wqUXCp0KXRDMM+B9gVOjti7RsNHUP8K8B5hI5bip0STTl7v5uaJfF1cCc0BYA7wL5gSYTOU4qdEk0rXfTbG51u5mWvY0a+ejPRVon5RI5bip0kY+qoOUSYh9cWKQg0DQi7aBCF/moPwIZZrYauJMI7VMt0hm026KISJzQCF1EJE6o0EVE4oQKXUQkTqjQRUTihApdRCROqNBFROKECl1EJE6o0EVE4sT/AdxrRDxTqizPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Write your code here\n",
    "z = np.arange(0.1, 1.7, 0.1)\n",
    "mu = cosmo.distmod(z)\n",
    "plt.plot(z, mu)\n",
    "plt.xlabel(\"z\")\n",
    "plt.xlabel(\"mu\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's assume that the absolute magnitude of our supernovae is `M0 = -19.3` magnitudes\n",
    "\n",
    "Use the `distmod` method to calculate the apparent magnitude $m$ of these supernovae, as a function of redshift.\n",
    "\n",
    "Make a plot of $m$ vs $z$.\n",
    "\n",
    "HINT: Re-read the help on `distmod` for the relation between apparent and absolute magnitude.  Distmod returns values with units; to remove the units from a quantity `x` you can write `x = x.value`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'm')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VXed//HXJyTsaxb2hATKWmiBhq3QTau2VbtaLVas3dBxXKodtdatMy7j2rF1nBmpWLViFy22VWkrBSyFAiXQsIY9bCFkBRIIIcv9/P7IbX+xsgTIuecm9/18PHj05tx7c96F3LzP95zvOcfcHRERSVxJYQcQEZFwqQhERBKcikBEJMGpCEREEpyKQEQkwakIREQSnIpARCTBqQhERBKcikBEJMElhx2gJdLT0z07OzvsGCIibcqaNWvK3T3jTK9rE0WQnZ1NXl5e2DFERNoUM9vTktdp15CISIJTEYiIJDgVgYhIglMRiIgkOBWBiEiCC6wIzCzTzJaY2WYz22Rmn2/23GfNbEt0+Q+DyiAiImcW5PTRBuB+d19rZj2ANWa2EOgH3ABc7O4nzKxvgBlEROQMAisCdy8GiqOPq82sABgE3At8391PRJ8rDSqDiEhbVFVbz5o9h1hdWMmsaUMY0KtLoOuLyQllZpYNTABWAT8CLjOz7wK1wL+5++qTvGc2MBsgKysrFjFFREJRVn2C1bsreaOw6c+Wg1VEHJKTjEuG9Gn7RWBm3YFngfvcvcrMkoFUYCowCXjGzIa6uzd/n7vPAeYA5ObmOiIi7YC7s7eyhjcKK1m9u5LVuw9RWH4MgC4pHZiQ1ZvPvms4U3JSGZ/Vm64dg99eD3QNZpZCUwnMc/f50cX7gfnRX/xvmFkESAfKgswiIhKGSMTZWlL99hb/6t2VlFSdAKBXlxQmZfdh5uRMJmWnMnZQL1I6xH4yZ2BFYGYGzAUK3P3hZk89B1wFLDGzEUBHoDyoHCIisVTXEGFD0ZGmrf3oL/6q2gYA+vfszJScNCblpDI5O5XhfbuTlGQhJw52RDAdmAVsMLP86LIHgV8BvzKzjUAdcMc7dwuJiLQVJxoayd97mBW7Kli5q4L8fYeprY8AMDSjG9eNG8Ck7FQm56QyuE8XmraR40uQs4aWAaf6P/5YUOsVEQlSQ2OE9UVHWLGzghU7K8jbU0ltfQQzGDOgJzMnZzE5O5Xc7FQyenQKO26LtInLUIuIhKUx4hQUV7FiZwWv7yxn9e5DHD3RtKtnVP8e3DYpi2nD0piak0avrikhpz03KgIRkWbcne2lR3l9Rzmv76xgVWElR47XA027em4YP5BLh6UzZWgq6d3bxhb/magIRCShuTuF5cdYsaui6Rf/rgrKj9YBkJnahfdd2I9Lh6UzdWga/Xt1DjltMFQEIpJwSqpqeW17+dtb/QeraoGmWT2XDc9g2rA0pg1NIzO1a8hJY0NFICLtXm19I28UVvLa9jKWbitna0k1AGndOjJ1WBqXRn/x56R3i8tZPUFTEYhIu+PubCs5ytJtZSzdXsYbhZWcaIjQsUMSk3L6cPPEUVw2PIPRA3ok5C/+d1IRiEi7UHH0BMt2lPPa9nJe21729tm7w/t25/YpQ7h8RDpTctLo0rFDyEnjj4pARNqkuoYIa/ceYum2Ml7bXs7GA0dwh95dU5h+QTpXDM9gxvB0BvYO9oJt7YGKQETahLdm97y1xb9iZwXH6hpJTjImZvXhi1eP4LIRGYwb1IsOcXDZhrZERSAicau2vpFVhZUsLihh8dZS9lUeB2BIWldumjiIy6MzfHp0bpsncsULFYGIxJXSqlqWbC1lUUEpy3aUU1PXSOeUJGZckM7sy4dx+fB0hqR1Cztmu6IiEJFQuTsbi6pYtKWExVtKWb//CAADe3Xm5omDePeofkwblkbnFB3kDYqKQERirqaugeU7Klgc/eVfUnUCM5iQ2ZsvvW8k7xrVl1H9NbUzVlQEIhIT+w/VsGRLKYu2lPL6zgrqGiJ075TMFSMyeNeovlw5MoO0dnLtnrZGRSAigWiMOPn7DrGooJTFW0rZcrDpbN4haV352JQhXD26L7nZqXRMjv0dueQfqQhEpNWcaGhk+Y5yXtxwkEVbSqk8VkeHJGNSdh++dt1o3jW6L0MT9DIO8UxFICLnpaaugVe3lvHSpoMsLiil+kQDPTol867Rfbl6dD8uH5FBry6a3hnPVAQictaqautZsqWUFzcc5O/bSqmtj9CnawrXjRvANWP7c+kFaXRK1iyftkJFICItcuhYHQs3l/DixmKW76igrjFC3x6duPWSTK4d25/JOakkd9D+/rYosCIws0zgt0A/wIE57v6ImT0E3AuURV/6oLsvCCqHiJy70qpaXt50kJc2HWTlrkoaI87gPl34+LQhXDuuPxMy+5Ckyzm0eUGOCBqA+919rZn1ANaY2cLoc//l7j8OcN0ico72H6rhpY0HeWnjQdbsPYR70y0aP3XFUK4dO4ALB/bUwd52JrAicPdioDj6uNrMCoBBQa1PRM7d7vJj/HVDMS9tPMiGoqYze0cP6MkXrh7BtWP7M7xfj5ATSpBicozAzLKBCcAqYDrwGTP7OJBH06jh0EneMxuYDZCVlRWLmCIJpbS6lr+uL+a5/AOs23cYgPGZvfnqtaO4Zmx/Xc8ngZi7B7sCs+7Aq8B33X2+mfUDymk6bvBtYIC733W675Gbm+t5eXmB5hRJBNW19by8qYTn84tYvqOciMOFA3tyw/iBfOCigbp2fztjZmvcPfdMrwt0RGBmKcCzwDx3nw/g7iXNnn8M+EuQGUQSXV1DhL9vLeX5dQd4ZXMJJxoiZKZ24dNXXsCNEwZyQV/t9kl0Qc4aMmAuUODuDzdbPiB6/ADgJmBjUBlEElUk4ryxu5Ln84tYsOEgR47Xk9atI7dNyuT68YOYmNVbB3zlbUGOCKYDs4ANZpYfXfYgMNPMxtO0a2g38MkAM4gkDHenoLia5/OLeGHdAYqP1NK1Ywfed2F/rh8/kBkXpJOief5yEkHOGloGnGyTQ+cMiLSifZU1vLDuAM+9WcT20qMkJxlXjMjggWtH8Z4x/ejaUeeNyunpJ0SkDao4eoIFG5pm/KzZ0zTpblJ2H75z41iuGzeA1G4dQ04obYmKQKSNaIw4y3aU8/TqvSzcXEJ9ozOyXw++fM1IPnjRQDJTu4YdUdooFYFInCs6fJw/5O3jD3n7KTp8nD5dU7hjWja3XDKY0QN6hh1P2gEVgUgcqmuIsKighKdW72Pp9qbLcs24IJ0HrxvN1WP66sqe0qpUBCJxZEfpUZ7J28eza/ZTcayOAb0689l3DefWSwZr148ERkUgErKaugYWbDjI06v3snr3IZKTjKtH9+MjkzO5fHgGHXR1TwmYikAkBO7OxqIqnlq9lxfyD1B9ooGh6d346rWjuHniYDJ66CbuEjsqApEYOlJTz/PrinjqjX1sLq6iU3IS779oALdNymJSdh+d7SuhUBGIBMzdWVVYydOr97FgQzEnGiJcOLAn375xLNdfPFD385XQqQhEAlJb38jz+UU8vnw3Ww5W06NzMh/OzeQjkzIZO6hX2PFE3qYiEGllJVW1PLFiD/NW7eFQTT2jB/Tkh7dcxAcvHkiXjpr2KfFHRSDSSvL3Hebx5YX8dX0xje68Z3Q/7pyew9Shqdr3L3FNRSByHuobI7y08SCPLy9k7d7DdO+UzB2XZnPHtGyy0jTvX9oGFYHIOTh0rI4nV+/liRV7KD5Sy5C0rjz0wTF8KDeT7p30sZK2RT+xImdhW0k1jy/fzZ/e3E9tfYTpF6TxnRvHctXIviTpxC9po1QEImcQiTh/31bK48t389r2cjolJ3HThEF8Yno2o/rrom/S9qkIRE7h2IkG/rhmP79+fTeF5cfo17MTX3rfSGZOztL1/qVdURGIvEPR4eM8vqyQp/P2UV3bwMWZvXnktvFcN26AbvUo7ZKKQCRqX2UN//P3HfxxzX4iDteNG8Cd07OZmNUn7GgigQqsCMwsE/gt0I+mG9XPcfdHmj1/P/BjIMPdy4PKIXImu8uP8fMlO5j/ZhEdzLhtUhafunIYg3p3CTuaSEwEOSJoAO5397Vm1gNYY2YL3X1ztCTeC+wNcP0ip7Wz7Cg/X7yD5/KLSOmQxKypQ/jUFcPo36tz2NFEYiqwInD3YqA4+rjazAqAQcBm4L+ALwPPB7V+kVPZVlLNfy/ewZ/XH6BzcgfunpHDvZcPpW8PFYAkppgcIzCzbGACsMrMbgCK3H3d6U67N7PZwGyArKysGKSU9q6guIqfLd7OixsP0iWlA5+8fBj3XJZDendd+18SW+BFYGbdgWeB+2jaXfQgTbuFTsvd5wBzAHJzcz3IjNK+bSw6wqOLtvO3zSV075TMv155AXfNyNEUUJGoQIvAzFJoKoF57j7fzMYBOcBbo4HBwFozm+zuB4PMIoknf99hfrZoO4u2lNKzczL3XT2cOy/NoVdXXf9fpLkgZw0ZMBcocPeHAdx9A9C32Wt2A7maNSStac2eSh5dtINXt5XRu2sK//beEXz80mx6dlYBiJxMkCOC6cAsYIOZ5UeXPejuCwJcpySwVbsqeHTxdpbvqCC1W0e+cs0oZk0boovAiZxBkLOGlgGnvQqXu2cHtX5JHCt3VfBfC7exqrCS9O6d+Np1o7l9ahZdO6oARFpCnxRps3aVHeV7C7bwSkEJ/Xp24psfGMPMyVm6C5jIWVIRSJtzuKaORxZt54kVe+ic0oEvXzOSu6bn0DlFBSByLlQE0mbUNUT43co9PLJoO9W19XxkUhZffM8IMnroPACR86EikLjn7rxSUMr3FhRQWH6My4an87X3j9a9AERaiYpA4tqmA0f4zl8KWLGrgmEZ3Xj8E5O4cmSGbgYv0opUBBKXSqtq+fHftvKHNfvp3SWF/7jhQmZOztL9AEQCoCKQuHK8rpHHXtvF/726k/rGCPfMyOEz7xpOry46GUwkKCoCiQuRiPP8uiJ++NJWio/Ucu3Y/jxw7SiGpHULO5pIu6cikNCt3l3Jd/6ymXX7jzBuUC9++pHxTBmaFnYskYShIpDQ7K2o4fsvFbBgw0H69ezET269mJsmDCIpSQeCRWJJRSAxV1Vbz88X7+Dx5bvpkGTcd/VwZl8+VJeEEAmJPnkSM+7OX9YX89ALm6isqeOWiYP5t/eO1K0hRUKmIpCYKK2u5RvPbeTlTSVcPLgXv75zMuMG9wo7loigIpCAuTvP5Rfx0AubOV7fyAPXjuKeGTkk63wAkbihIpDAHDxSy9f+tIFFW0qZmNWbH37oYi7o2z3sWCLyDioCaXXuzh/W7Ofbf9lMfWOEr79/NHdOz6GDZgOJxCUVgbSqA4eP88D8DSzdVsbk7FR+8KGLyEnXSWEi8UxFIK3C3XnyjX18b0EBjRHn36+/kFlTh+icAJE2QEUg521fZQ0PzF/P8h0VTBuaxg9uuYistK5hxxKRFgqsCMwsE/gt0A9wYI67P2Jm3wZuACJAKfAJdz8QVA4JTiTi/G7VHr7/4hYM+M6NY/no5CyNAkTamCBHBA3A/e6+1sx6AGvMbCHwI3f/BoCZfQ74JvCpAHNIAPZUHOPLf1zPqsJKLhuezn/ePI7BfTQKEGmLAisCdy8GiqOPq82sABjk7pubvawbTaMFaSMiEefXr+/mRy9vJTnJ+MEt4/hwbqZuFCPShsXkGIGZZQMTgFXRr78LfBw4Alx1ivfMBmYDZGVlxSKmnMGusqN8+Y/rydtziKtGZvC9m8cxoFeXsGOJyHky92A3yM2sO/Aq8F13n/+O574KdHb3b53ue+Tm5npeXl6AKeV0GiPO3GW7+MnfttEpOYlvffBCbp44SKMAkThnZmvcPfdMrwt0RGBmKcCzwLx3lkDUPGABcNoikPDsLj/GF57J5829h7l6dD++e9NY+vXUReJE2pMgZw0ZMBcocPeHmy0f7u7bo1/eAGwJKoOcnyVbSvncU2+SZMYjt43n+osHahQg0g4FOSKYDswCNphZfnTZg8DdZjaSpumje9CMobgTiTj/8/cd/GThNkb378kvZl1CZqpmBIm0V0HOGloGnGzzcUFQ65Tzd/REA/c/k8/Lm0q4cfxA/vPmi+jSsUPYsUQkQDqzWN62s+won3xiDYXlx/jGB8Zw1/Rs7QoSSQAqAgHglc0lfOHpfFKSk3ji7slcOiw97EgiEiMqggQXiTiPLt7OT1/ZzthBPfnFrFwG9da5ASKJREWQwKpq6/ni0+t4paCEmycO4ns3jaNzio4HiCQaFUGC2lF6lNlP5LGnooaHPjiGOy7V8QCRRKUiSEAvbzrI/c+so1NyEvPumcLUoWlhRxKRELWoCMwsF/gaMCT6HgPc3S8KMJu0skjE+ekr23h08Q4uHtyL//3YJQzU8QCRhNfSEcE84EvABppOBJM25sjxer7wdD6Lt5Ry6yWD+faNY3U8QESAlhdBmbu/EGgSCcz2kmpmP7GGfZU1fPvGsXxsSpaOB4jI21paBN8ys18Ci4ATby08xYXkJI68tLGY+59ZR5eOyTw5eyqTslPDjiQicaalRXAnMApI4f/vGnJARRCnGiPOwwu38vMlO5mQ1Zv/vf0S+vfSVUNF5J+1tAgmufvIQJNIqzlSU8/nnnqTV7eVMXNyJg9dfyGdknU8QEROrqVF8LqZjXnHbSYlDm09WM3sJ/I4cPg437tpHB+doru7icjptbQIpgL5ZlZI0zECTR+NQxuLjnD7L1fRKTmJp2ZP45IhfcKOJCJtQEuL4JpAU8h523ygio/NXUW3jh14+pPTdP8AEWmxFhWBu+8JOoicu4LiKm7/5Uq6pHTgqdkqARE5O0lhB5Dzs+VgVXR3UAeemj2VrDSVgIicHRVBG7b1YDUffWwVKR2Mp2ZPZUhat7AjiUgbpCJoo7aVVPPRx1aSnGQ8NXsa2ekqARE5N4EVgZllmtkSM9tsZpvM7PPR5T8ysy1mtt7M/mRmvYPK0F5tj5ZAhyTjydlTyVEJiMh5CHJE0ADc7+5jaJp++q9mNgZYCIyNTj3dBnw1wAztzo7So8x8bBVmxu/vncqwjO5hRxKRNi6wInD3YndfG31cDRQAg9z9b+7eEH3ZSmBwUBnam51lR5n52EoAnrx3Chf0VQmIyPmLyTECM8sGJgCr3vHUXcCLp3jPbDPLM7O8srKyYAO2AbvKjjJzzkrcPVoCPcKOJCLtROBFYGbdgWeB+9y9qtnyr9G0+2jeyd7n7nPcPdfdczMyMoKOGdcKy48x87GVNEac3987leH9VAIi0noCvVWlmaXQVALzml+y2sw+AXwAeLe7e5AZ2rrd5ceYOWcl9Y3Ok/dOZYRKQERaWWBFYE13PpkLFLj7w82WXwN8GbjC3WuCWn97sKeiaSRwoqGRJ2dPZWR/lYCItL4gRwTTgVnABjPLjy57EHgU6AQsjN4la6W7fyrAHG3S3ooaZs5ZSW19I/Pumcqo/j3DjiQi7VRgReDuy2i6Suk7LQhqne3FvsoaZj62kpr6RubdM4UxA1UCIhIcnVkcZ/ZV1nDbnJUcPdHA7+6ewoUDe4UdSUTaORVBHNl/qGkkUF1bz7x7pjB2kEpARIKnIogTRYePM/OxlVQdr2fePVNVAiISM4FOH5WWOXD4ODPnrORwTT2/u3sK4warBEQkdjQiCFnxkaaRwKFjdTxx9xQuztQ1+EQktjQiCFFVbT23P7aKiqN1PHH3ZMarBEQkBCqCkLg7X//TRvZU1vD7e6YwIUs3mheRcGjXUEieXVvEC+sOcN+7hzNlaFrYcUQkgakIQrCr7CjffH4jU3JS+fRVF4QdR0QSnIogxuoaInz+qXw6Jifx09vG0yHpZCdfi4jEjo4RxNiP/7aVDUVH+MWsSxjQq0vYcURENCKIpVe3lTFn6S5un5LF+y7sH3YcERFARRAz5UdPcP8z6xjRrzvf+MCYsOOIiLxNu4ZiIBJx7n9mHVW19fzunsl0TukQdiQRkbdpRBADv1peyKvbyvj6+0frvgIiEndUBAHbWHSEH7y0hatH92PW1CFhxxER+ScqggAdO9HA5558k9RuHfnhhy4iekc2EZG4omMEAfr3P2+isOIY8+6ZQmq3jmHHERE5KY0IAvKX9Qd4Jm8/n75yGJcOSw87jojIKQVWBGaWaWZLzGyzmW0ys89Hl98a/TpiZrlBrT9M+ypr+Or8DYzP7M19V48IO46IyGkFuWuoAbjf3deaWQ9gjZktBDYCNwO/CHDdoWlojPD5p97EHR69bQIpHTToEpH4FlgRuHsxUBx9XG1mBcAgd18ItNsDp48u2s7avYd55LbxZKV1DTuOiMgZxWRz1cyygQnAqlisLywrd1Xw30t2cMvEwdwwflDYcUREWiTwIjCz7sCzwH3uXnUW75ttZnlmlldWVhZcwFZyuKaOLzydT1ZqV/79hgvDjiMi0mKBFoGZpdBUAvPcff7ZvNfd57h7rrvnZmRkBBOwlbg7X3l2PeVHT/CzmRPp3kmzckWk7Qhy1pABc4ECd384qPXEg9+/sZeXN5XwpfeNZNzgXmHHERE5K0Fuuk4HZgEbzCw/uuxBoBPwMyAD+KuZ5bv7+wLMEahtJdX8x583c9nwdO6ZMTTsOCIiZy3IWUPLgFNNDfpTUOuNpdr6Rj735Jt075TMTz58MUm625iItEHamX0e/nNBAVsOVvP4nZPo26Nz2HFERM6JznY6R69sLuE3K/Zw1/QcrhrZN+w4IiLnTEVwDg4eqeVLf1zHmAE9+cq1I8OOIyJyXlQEZ6kx4nzxmXxq6yP87KMT6JSsu42JSNumYwRn6RdLd/L6zgp+cMs4hmV0DzuOiMh504jgLBQUV/GTv23j/RcN4MO5mWHHERFpFSqCszBn6S46Jyfx3RvHttuL5olI4lERtFBJVS1/XneAW3Mz6d1VdxsTkfZDRdBCv12xm0Z37pyeHXYUEZFWpSJogeN1jfx+1V7eM7ofQ9K6hR1HRKRVqQha4E9vFnGopp67Z+SEHUVEpNWpCM7A3fnV8kLGDurJ5JzUsOOIiLQ6FcEZvLqtjB2lR7lreo5mColIu6QiOIO5ywrp26MTH7hoYNhRREQCoSI4jW0l1by2vZyPTxtCx2T9VYlI+6Tfbqfx+PJCOiUn8dEpQ8KOIiISGBXBKVQeq2P+2iJunjiY1G46gUxE2i8VwSnMW7mHEw0R7p6RHXYUEZFAqQhO4kRDI79duYcrRmRwQd8eYccREQlUYEVgZplmtsTMNpvZJjP7fHR5qpktNLPt0f/2CSrDufrLumLKqk9wl04gE5EEEOSIoAG4393HAFOBfzWzMcADwCJ3Hw4sin4dN9yducsKGd63O5cPTw87johI4AIrAncvdve10cfVQAEwCLgB+E30Zb8Bbgwqw7lYVVjJ5uIq7pqhE8hEJDHE5BiBmWUDE4BVQD93L44+dRDoF4sMLTV3WSGp3Tpy04RBYUcREYmJwIvAzLoDzwL3uXtV8+fc3QE/xftmm1memeWVlZUFHROA3eXHeKWghNunZNE5RfciFpHEEGgRmFkKTSUwz93nRxeXmNmA6PMDgNKTvdfd57h7rrvnZmRkBBnzbb9+fTfJScasqTqBTEQSR5CzhgyYCxS4+8PNnnoBuCP6+A7g+aAynI0jx+t5Jm8fH7xoIH17dg47johIzCQH+L2nA7OADWaWH132IPB94BkzuxvYA3w4wAwt9vTqvdTUNWrKqIgknMCKwN2XAaeadvPuoNZ7LhoaI/zm9T1MyUll7KBeYccREYkpnVkMvLyphKLDx3UHMhFJSCoCYO6yXQxJ68q7R8fVTFYRkZhI+CJ4c+8h1u49zJ2XZtMhSSeQiUjiSfgimLuskB6dk7k1NzPsKCIioUjoIig6fJwXNx7ktkmZdOsU5AQqEZH4ldBF8NvXd+Pu3HFpdthRRERCk7BFcOxEA0++sZdrxw5gcJ+uYccREQlNwhbBs2v3U1XboBPIRCThJWQRRCLO48t3Mz6zN5cMibv74oiIxFRCFsHiLaUUlh/TCWQiIiRoEcxdVsiAXp25Zmz/sKOIiIQu4Ypg04EjrNhVwR2XZpPSIeH+90VE/knC/SZ8fPluuqR0YOakrLCjiIjEhYQqgtLqWl7IP8CtuYPp1TUl7DgiInEhoYrgdyv3Uh+JcOd0HSQWEXlLwhRBbX0j81bu4d2j+pKT3i3sOCIicSNhiuD5/CIqjtVxl0YDIiL/ICGKwN2Zu6yQUf17MG1YWthxRETiSkIUwfIdFWwrOcrdM3Iw0z0HRESaC6wIzOxXZlZqZhubLbvYzFaY2QYz+7OZ9Qxq/c3NXbaL9O6duH78wFisTkSkTQlyRPBr4Jp3LPsl8IC7jwP+BHwpwPUDsKP0KEu2ljFr6hA6JXcIenUiIm1OYEXg7kuByncsHgEsjT5eCNwS1Prf8vjyQjomJ3H7VJ1AJiJyMrE+RrAJuCH6+FYg0PtDHjpWx7Nr93PT+EGkd+8U5KpERNqsWBfBXcCnzWwN0AOoO9ULzWy2meWZWV5ZWdk5rez3b+yltj7CnTOyz+n9IiKJIKZF4O5b3P297n4J8CSw8zSvnePuue6em5GRcU7r69ujEx/OHcyo/jE5Ji0i0ibF9I7tZtbX3UvNLAn4OvB/Qa7v1txMbs0NdO+TiEibF+T00SeBFcBIM9tvZncDM81sG7AFOAA8HtT6RUSkZQIbEbj7zFM89UhQ6xQRkbOXEGcWi4jIqakIREQSnIpARCTBqQhERBKcikBEJMGpCEREEpy5e9gZzsjMyoA9YeeISgfKww5xBsp4/uI9H8R/xnjPB+0/4xB3P+OlGdpEEcQTM8tz99ywc5yOMp6/eM8H8Z8x3vOBMr5Fu4ZERBKcikBEJMGpCM7enLADtIAynr94zwfxnzHe84EyAjpGICKS8DQiEBFJcCqCUzCza8xsq5ntMLMHTvL8F81ss5mtN7NFZjYk3jI2e90tZuZmFtPZES3JZ2Yfjv49bjKz38cyX0symlmWmS0xszej/9bXxTjfr8ys1Mw2nuJ5M7NHo/nXm9nEOMt3ezTXBjN73cwujmW+lmRs9rqwZsWaAAAEVUlEQVRJZtZgZh+KVbZm6z5jRjO70szyo5+VV1s1gLvrzzv+AB1ounvaUKAjsA4Y847XXAV0jT7+F+DpeMsYfV0PYCmwEsiNp3zAcOBNoE/0677x9ndI0/7Zf4k+HgPsjnHGy4GJwMZTPH8d8CJgwFRgVZzlu7TZv++1sc7XkozNfhYWAwuAD8VbRqA3sBnIin7dqp8VjQhObjKww913uXsd8BRwQ/MXuPsSd6+JfrkSGBxvGaO+DfwAqI1lOFqW717g5+5+CMDdS+MwowNv3eu0F003VIoZd18KVJ7mJTcAv/UmK4HeZjYgNunOnM/dX3/r35dwPict+TsE+CzwLBDrn0GgRRk/Csx3973R17dqThXByQ0C9jX7en902ancTdNWWSydMWN0N0Gmu/81lsGiWvJ3OAIYYWbLzWylmV0Ts3RNWpLxIeBjZrafpq3Fz8YmWoud7c9qmML4nJyRmQ0CbgL+N+wspzEC6GNmfzezNWb28db85jG9Z3F7ZGYfA3KBK8LO0lz0vtAPA58IOcrpJNO0e+hKmrYUl5rZOHc/HGqqfzQT+LW7/8TMpgFPmNlYd4+EHawtMbOraCqCGWFnOYmfAl9x94iZhZ3lVJKBS4B3A12AFWa20t23tdY3l39WBDS/6/3g6LJ/YGZXA18DrnD3EzHK9pYzZewBjAX+Hv3h7g+8YGbXu3teHOSDpq3XVe5eDxRG72c9HFgdg3zQsox3A9cAuPsKM+tM07VfQtmFcBIt+lkNk5ldBPwSuNbdK8LOcxK5wFPRz0k6cJ2ZNbj7c+HG+gf7gQp3PwYcM7OlwMVAqxSBdg2d3GpguJnlmFlH4DbgheYvMLMJwC+A60PYt33GjO5+xN3T3T3b3bNp2j8bqxI4Y76o52gaDWBm6TQNf3fFKF9LM+6laSsMMxsNdAbKYpjxTF4APh6dPTQVOOLuxWGHeouZZQHzgVmttfXa2tw9p9nn5I/Ap+OsBACeB2aYWbKZdQWmAAWt9c01IjgJd28ws88AL9M0m+BX7r7JzP4DyHP3F4AfAd2BP0S3JPa6+/VxljE0Lcz3MvBeM9sMNAJfiuUWYwsz3g88ZmZfoOnA8Sc8Om0jFszsSZrKMj16nOJbQEo0///RdNziOmAHUAPcGatsLcz3TSAN+J/o56TBY3yRtxZkDN2ZMrp7gZm9BKwHIsAv3f2002HPav0x/JkWEZE4pF1DIiIJTkUgIpLgVAQiIglORSAikuBUBCIiCU5FICKS4FQEIiIJTkUgcg7M7FPRa8Pnm1mhmS0JO5PIudIJZSLnwcxSaLqO/Q/d/c9h5xE5FxoRiJyfR4DFKgFpy3StIZFzZGafAIYAnwk5ish50a4hkXNgZpcAvwEua3YHLpE2SbuGRM7NZ4BUYEn0gPEvww4kcq40IhARSXAaEYiIJDgVgYhIglMRiIgkOBWBiEiCUxGIiCQ4FYGISIJTEYiIJDgVgYhIgvt/e1j9/MadyhEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# write your code here\n",
    "# m = M + mu\n",
    "M = -19.3\n",
    "m = M + mu.value\n",
    "plt.plot(z, m)\n",
    "plt.xlabel(\"z\")\n",
    "plt.ylabel(\"m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the final step in this exercise, write a `model` function to predict the apparent magnitude for any set of parameters, for any set of redshifts.\n",
    "We will use this tomorrow!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Complete this function\n",
    "def model(H0, Omega_matter, Omega_lambda, M0, z):\n",
    "    cosmo = LambdaCDM(H0=H0, Om0=Omega_matter, Ode0=Omega_lambda)\n",
    "    mu = cosmo.distmod(z)\n",
    "    m = mu.value + M0\n",
    "    return m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2: CAMB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "CAMB, the Code for Anisotropies for the Microwave Background, is Boltzmann code:  it integrates the Boltzmann equation, $df/dt = C[f]$, which describes how collisional processes affect densities, momenta, and higher order moments of the distribution of light, matter, and other cosmic components like neutrinos.\n",
    "\n",
    "First, create an array or redshift values from 0 .. 2 with spacing 0.1, at which CAMB will evaluate the matter power spectrum.\n",
    "\n",
    "HINT: `np.arange`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Complete this code:\n",
    "\n",
    "redshifts = np.arange(0., 2., 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the following parameters to CAMB, which describe the cosmology and tell it that we want the late Universe power spectrum, which is calculated from the transfer functions.  Fill in the `...` below to use all these and your `redshifts` variable from above as named arguments to `camb.set_params`, which will create the settings CAMB needs.\n",
    "\n",
    "$\\Omega_c h^2 = 0.1193$\n",
    "\n",
    "$\\Omega_b h^2 = 0.02242$\n",
    "\n",
    "$\\Omega_k = 0.0$\n",
    "\n",
    "$H_0 = 67.66$\n",
    "\n",
    "HINT: Have a look at the camb documentation on https://camb.readthedocs.io/en/latest/model.html for parameter names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: redshifts have been re-sorted (earliest first)\n"
     ]
    }
   ],
   "source": [
    "# Complete this code:\n",
    "\n",
    "p = camb.set_params(redshifts=redshifts, WantTransfer = True, omch2=0.1193, ombh2=0.022442, omk=0., H0=67.66)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run CAMB as described in the lecture using your settings `p`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Complete this code\n",
    "results = camb.get_results(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CAMB can generate an interpolator for the matter power, an object which stores $P(k, z)$ evaluated at grid points in $k$ and $z$, and knows how to interpolate between them to other values.\n",
    "\n",
    "Run the following cell to get a list of all the attributes on the results that start with `get_`.  Find the one that will compute this estimator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "get_BAO\n",
      "get_Omega\n",
      "get_all_fields\n",
      "get_background_densities\n",
      "get_background_outputs\n",
      "get_background_redshift_evolution\n",
      "get_background_time_evolution\n",
      "get_cmb_correlation_functions\n",
      "get_cmb_power_spectra\n",
      "get_cmb_transfer_data\n",
      "get_cmb_unlensed_scalar_array_dict\n",
      "get_dark_energy_rho_w\n",
      "get_derived_params\n",
      "get_fsigma8\n",
      "get_lens_potential_cls\n",
      "get_lensed_gradient_cls\n",
      "get_lensed_scalar_cls\n",
      "get_linear_matter_power_spectrum\n",
      "get_matter_power_interpolator\n",
      "get_matter_power_spectrum\n",
      "get_matter_transfer_data\n",
      "get_nonlinear_matter_power_spectrum\n",
      "get_redshift_evolution\n",
      "get_sigma8\n",
      "get_sigma8_0\n",
      "get_sigmaR\n",
      "get_source_cls_dict\n",
      "get_tensor_cls\n",
      "get_time_evolution\n",
      "get_total_cls\n",
      "get_unlensed_scalar_array_cls\n",
      "get_unlensed_scalar_cls\n",
      "get_unlensed_total_cls\n"
     ]
    }
   ],
   "source": [
    "# Run this cell\n",
    "for d in dir(results):\n",
    "    if d.startswith('get_'):\n",
    "        print(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and call it to create an interpolator.  It takes no arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Complete this code\n",
    "P = results.get_matter_power_interpolator()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use your interpolator to plot P(k) at varying redshift values 0, 0.1, and 0.2, using the method `P.P(z, k)`.\n",
    "For `k`, use log-spaced values from $10^{-3}$ to $10^{0}$\n",
    "\n",
    "HINT: use np.logspace for the k values.\n",
    "\n",
    "What are the units of these quantities?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'P(k) / (Mpc/h)^3')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3WdUVFcbhuF700QUkKKigKCiIqDYu8bee48Re40tGltsMUZjiRpN7B17b1E0sXdFsfeKAlZAQPows78fY+pn0ERgKPtaa9YKM2cOzwHDO2dXIaVEURRFUf7OyNABFEVRlPRJFQhFURTlnVSBUBRFUd5JFQhFURTlnVSBUBRFUd5JFQhFURTlnVSBUBRFUd5JFQhFURTlnVSBUBRFUd5JFQhFURTlnUwMHeBj2NvbS1dXV0PHUBRFyVACAgJCpZS533dchi4Qrq6uXLhwwdAxFEVRMhQhxOMPOU41MSmKoijvpAqEoiiK8k6qQCiKoijvlKH7IBRFUf6JRqMhODiY+Ph4Q0cxGHNzc5ycnDA1Nf1P71cFQlGUTCk4OBhLS0tcXV0RQhg6TpqTUhIWFkZwcDAFCxb8T+dQTUyKomRK8fHx2NnZZcniACCEwM7O7qPuoNQdhJJpaHWS8JhEQqMTiIrTkKjVodHqkBIszEzIkc0YS3NT8lmbY25qbOi4ShrIqsXhNx97/apAKBmORqvjzvM3XA2O5M7zKB6GxvDwVQzPIuPQfeAW63kss+FiZ4FnfmtKOllTyjkXBe1zZPk/KIrhJSQk0KVLFwICArCzs2PTpk28a0Lw/v37GTJkCFqtll69ejF69OgUz6IKhJLuJSRpufQkgtMPwjjzIJQrwZEkJukAyJnNhEK5c1De1QZnW0dyW2bDPmc2rLObYmZihKmxvhU1NjGJ2AQtEXEankbEEfw6lkehMWw6H8Sq04EAONlkp7Z7Hmq756Gqm/3v71WUtLR8+XJsbGy4f/8+GzduZNSoUWzatOkvx2i1WgYMGMCBAwdwcnKifPnyNG/eHA8PjxTNogqEki5Fxmo4fOcFB26+4NidV8QkajESUMLRmi6VXPB2zoW3Uy6cbbN/1Kf+JK2O+6+iuRD4mqN3XrL5QhCrzzzGPmc22pR1pEM5ZwrlzpmCV6ZkJYsWLWLRokUAREZG4urqypEjR5J9z65du5g4cSIAbdu2ZeDAgUgp//Lv3N/fHzc3NwoVKgRAx44d2bVrlyoQSuYVl6jlwK0X7L4cwrG7r9BoJbkts9G8lCO1iuWmYiE7rLP/t+F6/8TE2Ah3ByvcHazoXMmFeI2W43dfsflCMMtOPGLxsYfUcc/DgNpulClgk6LfW0k73/x8g5tPo1L0nB75rfi6mWeyx/Tr149+/fqh0WioXbs2w4YNo0OHDty5c+f/jh02bBhdunQhJCQEZ2dnAExMTLC2tiYsLAx7e/vfj/3zMQBOTk6cO3cuha7sD6pAKAYlpeRaSCQb/IP4+cpTohOScLAyp3vVgjTycsDbKRdGRsncIWjiIfwhhN2HyCCIfgHRryA+ApISQJuoP87UAsxygLk1WDtCLhewcYW8nvrn3zI3Naa+pwP1PR14GRXPxvNBrDz1iNYLTlPNzZ4RDYrh7ZwrdX8oSqYzZMgQateuTbNmzWjWrJmh43wwVSAUg4jXaNl9+SmrTgdy81kU5qZGNCmRn7ZlnahY0PbdRUGnhWdXIMgfnl6EkIv6wsCfeqaNzSBnXjDPBSZmYJxN/3z0C0iMgbjXEBv6x/HCCHK7g2MZKFwHCteG7PoCkMfKnMF1itCzWkHWn3vC4uMPaDH/FO3KOjGiYTHyWJqn3g9ISVHv+6SfmlatWsXjx4+ZN28ewHvvIBwdHQkKCsLJyYmkpCQiIyOxs7P7y7G/HfOb4OBgHB0dUzy7KhBKmnr1JoHVZwJZd+4J4TGJFMtrybctvWhRKj9W5u9oPooMgbv74MERCDwB8ZH653M6gGNZ8GoN9kXBzg1sXPSF4X19EomxEBkM4Q/g6SV9obm1By6tBWEMLlWgRDv9ubNZkiObCb1rFKJjBWfmHbnPipOP2Hf9OSMbFqNzRZfk73CULC0gIICZM2dy4sQJjIz0gx7+3uH8d82bN8fX15fKlSuzdetWateu/X/9bOXLl+fevXs8evQIR0dHNm7cyPr161M8v5DyA8cFpiNCiGZAMzc3t9737t0zdBzlAzwJi2XJiQdsvhCMRqujjnteelRzpXKhd0xkCn8I17bB7Z/1dwygbxIq9AkU/AQKVNY3E6UknRaCz8PdX+D2Hgi9C6Y5wLMVVOoPDl6/H/rwVTRf777BiXuhVHOzZ3rbkjjmyp6yeZSPduvWLYoXL27QDN27d+eXX34hT548AJQrV45ly5Yl+574+Hh8fHy4dOkStra2bNy4kUKFCvH06VN69eqFn58fAH5+fnzxxRdotVp69OjB2LFj33m+d/0chBABUspy78ufIQvEb8qVKyfVfhDpW2BoDD8dvs/OyyEYC0HrMo70qVHo/0cGxb2Gq1vg6iYIefs7dSoP7k2gWBPIXfSDvp+UktikWKISokjUJZKoTUQiyW6Snewm2bE0syTbb81O/3wSfbG4uBqubwdNDBRrDNWHg1PZ37/PBv8gJu+9ibEQTG7lRYtSKX+Lr/x36aFApAcfUyBUE5OSKoLCY5l76B47LoVgYiToVsWVPjUKkdfqT+32UsLj0xCwCm7uAm0C5PWCut9AibZg7fSP53+T+IYbYTe49/oegZGBBEYF8jT6KWHxYcQlxSWbzdbclvw58uNk6YS7rTte9l542HlgaWapP0AIcK6gf9SbBP5L4OxCuFMb3JtCgykIG1c6VSxANTd7hm2+zJCNlwl4/JqxTYqTzUTN0lYyB3UHoaSo0OgE5h2+z7pzjzESgs6VXOj7SaG/duhq4uHaFji3CF5ch2xWULI9lOkC+bzfed6XsS859+wc/s/9ufzyMoFRgb+/ZmVmhau1K045ncidPTf22e2xymaFmbEZpkamCATx2nhiNbFEJkTyLOYZz2Ke8TjqMSHRIQAIBCVyl6CGYw1qONXA3db9r01fCW/g7CI4OVvfHFVlEFQfBmY50Gh1fP/LHZYcf0gp51ws+KwM+VWTk8GpOwg91cSkGFxcopalJx6y6NgDEpJ0tC/nzJA6RXCw/lNhiHsN/sv0hSE2FPJ4QMV++g5hM4u/nE8nddwIvcGRoCMcCTrC/Yj7AOTKlotSeUpRwr4EXvZeFLMphq257X+eLBcRH8HNsJtcfnWZE8EnuB52HYDC1oVpVaQVTQs1xS77n0aQRD2FA1/Dtc1gUxBaLYIClQDYd+0ZI7ZeJbuZMcu7lqOkkxoOa0iqQOipAqEYjE4n2Xk5hBn77/A8Kp6Gng6MaFiMwn/uY3jzAk7/qG9KSoyGIvWh8kAoWOMvI46klNwOv43fIz/2PdrHi9gXGAtjyuYtSzXHalTKV4litsUwEqm3BEZoXChHgo6w895OroZexUSY0KRQE3qV6IWrtesfBwaehJ39ISIIqg6GmmPA1Jx7L97QfdV5QqMTmNuxNA08HVItq5I8VSD0VIFQDOJKUARf777B5aAISjpZM66JBxUK2v5xwJsXcGouXFgOWg14tYGqQ/4yIgjgdfxrdj/YzY57O3gQ+QATIxOqOVajvkt9ajjVwDqbdRpfmd791/fZfHcz2+9tR6PTUN+lPgNKDfijUCS8gV/HQ8BKfdNY+9Vg48qrNwn0Xn2BK8ERjG/iQY9q/20tfuXjqAKhpwqEkqbCohOYsf8OmwOCsMuRjdGN3Gld2vGP+QBxr+HkHDi3WN/xXLIj1BgOdoV/P4eUkiuvrrD+1noOPjmIRqehZO6StCjcgvou9cllnn6aZ0LjQllzcw0bb28kUZvIZ8U/o6933z86te/shx199P/deikUbUC8RssXGy+z/8ZzBtZy48v6RdVKsWlMFQg9VSCUNKHTSTZdCGLavtvEJCTRvaorg+sUwfK3CW6aOP1on1NzID5K37dQc/RfCoNGp+GXwF9Ye3MtN8JuYGlqSXO35rQp0oYiNkUMdGUfJjQulB8v/sjO+zuxMbdhZPmRNC7YWP+H/3UgbPKB51fhk9FQczRaCWN3XGPj+SA6VyrApOZealJdGsqoBeJDl/vu0aMHe/bsIU+ePFy/fv0fz6eGuSqp7vbzKMZsv8bFJxFUKGjL5JZeFM379hO0TgfXt8LBbyAqGIo0gDoT/tKUFJcUx/Z72/G94cuzmGe4WrkyruI4mhVuhoWpxT981/TFPrs9k6pOokOxDkw5N4XRJ0Zz6MkhxlUah62NK/T8FfYMg2PTIPwBxi3mM7V1CXJZmLHo2AOi4pKY3d4bE7WMuJKMD1nuG6Bbt24MHDiQLl26pFoWVSCUZMVrtMw7fJ9Fxx5gld2Ume28aVPG8Y/mkqDzsH8UhATo2+FbLwbXar+/P1YTy4bbG/C94cvrhNeUyVOGsRXHUt2peqp2NqcmT3tP1jRaw6obq5h/eT4BLwL4tuq31HCqAS0XgL0bHJoEkcGIDusY3cgdq+wmzNh/B41Wx9yOpTEzyZjXrvw7qbXcN0CNGjUIDAxMjdi/UwVC+Uf+j8IZve0qD0NjaF3GkXFNPLDNYaZ/8c0LODgRrqwHy3zQcqG+r+HtejNxSXFsuL2BVddX8TrhNVUdq9KnRB/K5C1juAtKQcZGxvQs0ZMaTjUYc3IMAw4NoHeJ3gwoNQDj6l/qV4rd0R9W1AefnXxe041sJsZ8u+cmmnUBzP+sjJpQl5b2jYbn11L2nA4loNG0ZA9JreW+04oqEMr/iU1MYsb+O6w6HYizbXbW9KxA9SK59S9qk+D8Ujjynb7PodpQ/RIU2fTDWjU6DTvu7WDhlYWExoVSNX9V+pfqj3fud0+Ay+iK2BRhbeO1TD03laXXlnI19CrTq0/HzqsNWOaH9R1gRQPw2UnPakUxMxaM33WDPqsDWOxTVu2NnUWo5b6VTOHMgzBGbrtCUHgc3aq4MrJhMSzM3v4zCb4Ae77QfxJzqwsNp+ubU9CPSjr45CBzAubw5M0TSucpzaxPZmWaO4bkZDPOxsQqE/HO7c2Uc1PotLcTC+ouoLBLZei+F9a01heJztvwqVwGU2Mjvtpxjd6rL7C0SzlVJNLCez7pp6bUWO47ragCoQD6mdDT999m1elAXOws2NSnEhULvf1HGR8Fh76B88vB0gHa+YJHi98nud0IvcGM8zO4+PIibrncmFd7HjWcamS5YZ2tirSiqE1RBh4eiI+fD7NqzqJy/srQYz+sbgm+zaHzNjpWqIiREIzafpVevvoikd1MFYnMKLWW+04zUsoM+yhbtqxUPt6lJ69lrZlHpMuoPXLCzmsyJkHzx4u3/aSc6S7l19ZS7h0hZVzk7y+9in0lx54YK71WeckaG2vIzXc2S41W847vkLU8ffNUttrVSpbyLSW3392ufzIiWMq5paWcnE/KRyellFJuuRAkXUfvkZ8uOfPXn7mSIm7evGnoCLJbt24yX7580tvbW3p7e8uePXu+9z1xcXGybdu2snDhwrJ8+fLywYMHUkopQ0JCZKNGjX4/rmPHjtLBwUGamJhIR0dHuWzZsnee710/B+CC/IC/sWoeRBaWpNUx/8gDfjx8j7yW2fi+nTdV3d52hMWEgt8IuLFdv2ZS85/AST9sWqPTsPH2RhZcXkC8Nh4fDx/6lOhDTrOcyXy3rCU6MZrhx4Zz6ukphpYdSg+vHvDmOfg20y/P0WkjFKrJ9ovBDN9yhXKutqzoVp6c2dRNfUrJqPMgUpqaB6H8a4/DYvhi02UuPYmgZan8fNPCC+vsbye83dgBe4frd2+rNRaqfqHfvhO49PISk85M4n7Efao6VmVU+VEUtFZLSfxdTrOc/FT7J8aeHMsPAT8QkRDB0DJDEd38YHULfed1x/W0LlMHE2Mjhm66TJfl51jVo8K7d9ZTFANQBSKLkVKy7WIIX++6jrGR4MdPS9PcO7/+xZhQ2DtMvzdD/tLQ4mfI6wHoVz2dc3EO2+5twyGHA3NqzaG2swHbRjMAU2NTplafilU2K1ZeX0lUQhQTKk/AqOvP+iKx4VPouJ7m3nUxNRIM2nCJzsvO4du9Aja/DSdWFANSBSILiYrXMG7HdXZfeUqFgrb80KHUH1tl3vaDnwfr7xrqTIAqQ8DYBCkl+wP3M81/GpEJkXT37E4/734ZZvazoRkbGTO24liszKxYem0pGp2GSVUmYdx1t75IbPwUOqyjUYn6LDI24vP1F2m/+Axrelb861LpimIAqkBkEZeDIhi04SJPI+L5sl5RPq/lhrGR0I9Q2v8VXF4LeUtAl12Q1xOAZ9HPmHxuMseDj+Nl58WSeksoZlvMwFeS8QghGFxmMGbGZsy/PB+NVsN31b/DpMsuWNMSNn0G7Xyp69EY3+4V6L36Am0WnmZtr4oUtM9h6PhKFqbm+2dyOp1kyfEHtF14Gp0ONvetxKA6RfTF4clZWFRNPxu6+pfQ+zDk9URKyZa7W2i1uxXnn59nZPmRrG28VhWHj9TPux9Dyw5lX+A+Rh4ficbcUl+QHUrAZh+4uYvKhe3Y0LsScRot7Rad5uKT14aOrWRhqkBkYuExifT0Pc93frepWzwvfoOrU9bFVr83w+EpsLKR/sDu+/TNSiZmhESH0OdAHyadmYSXnRfbm2/Hx8MHYyM1Tj8l9PDqwcjyIznw+ADDjw5HY5YTfHaAY1nY0h2ubaWEkzWb+1bGwsyEjkvOsutyiKFjK1mUKhCZVMDjcJr8eIJT98OY1MKThZ3LYG1hql+WemUjOD5Dv3ZSv5NQoJK+8/ruNlrvas3VV1eZUHkCS+svxcnSydCXkun4ePgwusJoDgcdZtixYWhMLaDzNv3Wpdt7Q8Aq3PLkZOeAqpRyzsWQjZeZ9esddLqMOyRd+XAJCQl06NABNzc3Klas+M4F+YKCgqhVqxYeHh54enoyd+7cVMmiCkQmI6Vk6fGHdFh8FlNjI7Z/XoUulV31o42ubYVF1eHVXWizHFotBHMrXsW+YuDhgUw8MxEvey92tNhBu6Lt1AilVPRZ8c8YU3EMR4OOMvToUBJNssFnW6Fwbfh5CJycg20OM9b2rEiHcs78dPg+XVf68+pNgqGjK6nsz8t9Dx06lFGjRv3fMSYmJsyaNYubN29y9uxZ5s+fz82bN1M+zIfMpkuvDzWT+q8i4xJlb9/z0mXUHtl39QUZGZeofyEhWsqdn0v5tZWUy+pJGR74+3sOBB6Q1TZUk2XXlJVrb66VWp3WQOmzpk23N0mvVV6y74G+Mk4TJ6UmQcrN3fS/q1/HS6nTSZ1OJzeceyyLjvWTZb89IE/cfWXo2BlCephJvXDhwt9nUbu6usqaNWu+9z3169eXp0+fllJKqdFopJ2dndTpdMm+p3nz5vLXX39952sfM5NajWLKJG4+jaL/ugBCXscxrklxelYrqL8DeHETtnSD0Lv6VVdrfgXGJsRoYpjuP50d93fgYefB1OpTKWRdyNCXkeW0L9YeEyMTJp6eyMDDA/mp9k9kb7MMsufS7+cd8QTRciEdKxSgdAEbBq6/iM+Kc/SoWpARDYqphf4+0HT/6dwOv52i53S3dWdUhf//dP9nabHcd2BgIJcuXaJixYoff1F/k64KhBAiB3AMmCil3GPoPBnFtoBgxuy4Ri4LUzb2qUQ5V1uQEi6uAb/hkM1K3xFauBYAV19dZdTxUTyNeUrvEr3p790fU2M1e9dQWhdpjYmRCeNPjaf/wf7Mqz2PnE1mQy4X/Z4brx9Dx/UUc8jHroFV+c7vFstPPuLI7ZfMbO9NmQI2hr4E5T1Sa7nv6Oho2rRpw5w5c7Cyskqx8/7uQ24z/usDWAG8BK7/7fmGwB3gPjD6T89PAkYCTT/k/Fm9iSlBo5XjdlyTLqP2yPaLTsuXUfFvX4iWcntffTPFqmZSRj2XUkqZpE2SS68ulaV8S8l6W+rJC88vGDC98nd+D/1kKd9Ssv3P7WVYXJj+yVt79Qv8zXSX8vGZ3489cfeVrPzdQVlw9B45Ze9NGZuQZKDU6Vd6aGKSUsqVK1fKxo0bS61W33zbvn3735ud/vzw9fWVUn54E1NiYqKsX7++nDVrVrLf/2OamFK7QNQAyvy5QADGwAOgEGAGXAE8gHpAR6CbKhDv9zwyTraaf1K6jNL/gdAkve07eHlHynkV9KuvHpkqpVb/h+NFzAvZY38P6bXKS3559EsZmRCZzNkVQzkWdEyWXVNWNtvRTD6LfqZ/8tlVKeeUlHKijZTHZ0n59g9NZFyiHL3tinQZtUfWmHFYnrqv+ib+LD0UiAsXLkhPT08ZHh7+we+ZN2+e7Nu3r5RSyg0bNsh27dr93zE6nU76+PjIIUOGvPd86bZA6HPg+rcCURn45U9ff/X2MQWYA/wK7AKM3nfurFogzj8Kk+UmH5DFx++Te648/eOFa1v1nzanF5Ly/uHfnz4RfEJW31Bdll9bXm6/u/29HV6KYZ1/dl5WWldJ1ttST94Nv6t/Mi5Cys1d9XeFq1tJGfXs9+NP3X8lP5lxWLqM2iPHbL8qo+PV0uFSpo8CkVrLfZ84cUICskSJEr+fe+/eve88X7pe7lsI4QrskVJ6vf26LdBQStnr7dc+QEUp5cC3X3cDQuU/9EEIIfoAfQAKFChQ9vHjx6maPz2RUrLu3BO++fkGjrmys6RLOYrmtYSkRDgwHs4tAueK0G4VWOVHo9Pw06WfWHl9JUVsijCzxkwK5VId0RnBrbBbDDg0gBhNDN9/8j01nGro+5UCVuqXRjHNDk1mgVcbAOI1Wmb9eodlJx/hbGPBzHbeVChoa+CrMCy13Lfexyz3ne7mQUgpV/1TcXj7+hIpZTkpZbncuXOnZTSDSkjS8tX2a4zbeZ1qbvbsGlhNXxyinoFvU31xqPQ5dNsLVvl5HvOcHvt7sPL6StoVbcf6xutVcchAitsVZ0OTDbhYuTDo8CBW31iNBCjXQz+50bYwbO0Bm7tCTCjmpsaMbeLBpj6VAeiw5Aw/HLirJtcpH8UQBSIEcP7T105vn1P+wcs38XRaeo6N54MYWMuN5V3L6/dueHwaFteA59eh7QpoOBWMTTkZcpJ2P7fj7uu7zKgxgwmVJ2BuolYGzWjy5sjLqoarqO1cm+8vfM+I4yOITowG+yLQ4xf98ii398K88nB1C0hJhYK27BtSndalnZh76B791wUQnZBk6EtRMihDFIjzQBEhREEhhBn6jundBsiRIVwNjqD5T6e4+TSK+Z3KMLxBMYwEcG6xfneybJbQ+xB4tUGr0zLv0jz6H+xPbovcbGy6kUYFGxn6EpSPYGFqwayas/iizBccfHyQ9nvacyPsBhib6BdY7HcCbAvB9l76TYiinpIjmwkz25VkQlMPDt56SZsFpwkKjzX0pSgZUKoWCCHEBuAMUEwIESyE6CmlTAIGAr8At4DNUsob//K8zYQQSyIjI1M+dDqy63II7RadwcRYsK1/FZqUzAeaONjZH/aNhCL1oc8RyFOc8Phw+h/sz+Kri2np1pJ1jdepnd4yCSNhRM8SPVnZcCWJ2kR8/HxYfm05Wp0W8hSHnr9Cg+/g0XFYUAmubkYAPaoVxLd7BZ5HxdNxyVmeRsQZ+lLSXGr3saZ3H3v9ak/qdEink8z89Q4Ljj6gQkFbFn5WBruc2fR7GW/qDM8uQ80xUGMEGBlx9dVVhh0dxuv414ytNJbWRVob+hKUVBIRH8E3Z77h4JODeOf2Zkq1KbhYuehfDHsAO/pBsD94tICmc8DClmvBkXRaepbcltnY3K8y9jmzGfYi0sijR4+wtLTEzs4uS64rJqUkLCyMN2/eULDgXz8sfmgntSoQ6Ux0QhJfbLzMwVsv6FSxABObeWJmYqTvb9jkA9pEaL0UijVESv2+DVP9p5LXIi+za87Gw87D0JegpDIpJX6P/JhybgoarYbBZQbTyb2Tfkl2nRZO/whHvgPLfPDpRsjrgf+jcLqsOEdB+5xs7F1Jv7JvJqfRaAgODiY+Pt7QUQzG3NwcJycnTE3/+vtWBSIDCn4dSy/fC9x7Gc2Eph50qeyi/+Rzfrm+ScmmIHy6AeyLEJ8Uz+Szk9n1YBdVHasyvfp0rLNZG/oSlDT0IuYFk85O4njwcUrmLsmkKpMonKuw/sXgANjYCRKjofUScG/Csbuv6OV7nooF7VjdowJGRlnvU7Wil2GHuWZVAY/DaTHvFCERcazqXp6uVVwRuiTYMxT2DtMvA937ENgX4XnMc7ru78quB7voW7Iv82vPV8UhC8qbIy/zas9javWpPIl6Qruf27Hs2jJ934RTWX3/lH1RfaE4u4hPiuZmUgsvTt4PZeGxB4aOr2QAGbJAZLZO6h2Xgvl0yTkszU3Y8XlVqhfJDTFhsLolXFgBVb/QNxWYW3P++Xk67OnAk6gn/FjrRwaWHqh2e8vChBA0LdSUnS12UtO5JnMvzqXb/m4ERQWBVX7o7gfuTWH/KLi4ho7lnWlSMh+zD9wl4HG4oeMr6ZxqYjIgnU7yw8G7/HT4PpUL2bGwcxlyWZjpl+je0BHePIcW86Bke6SUrL+9nu/Pf08BqwLMrTVXjVJS/kJKyd5He/nu7HckySTGVxpPs8LNICkBNnwKD49A2xVEFW5K0x9PotVJ/AZXzxL9EcpfqSamdC5eo2XQxkv8dPg+Hcs749ujgr443NkPy+vp/6fuvg9KtidRm8jXp79mmv80qjtVZ33j9ao4KP/nt7uJ7S2242HnwZiTY5gdMButkQl0WANOFWBbb6yCj/PTp6V5ERXP2J3XDB1bScdUgTCAV28S6LjkLH7XnjGmsTtTW5fAzFjoN4jZ0BHs3PTtx05leRX7iu6/dGfH/R308+7H3FpzyWmW09CXoKRjDjkcWFp/Ke2Ltmfl9ZUMPjKYaCR02qTvk9jWE2+rGIbUKcKeq884eueloSMr6ZQqEGns7os3tJx/ijvP37Coc1n61CiM0Gpg1wA4MEE/fr37PrDKz43QG3Tc05F7r+/xQ80fGFBqAEZC/cqU9zM1MmV85fGMqziOUyGn6HewH7EmZtB+tX6DHkVcAAAgAElEQVRxx+296VPdhcK5czB+13XiErWGjqykQ+qvTRo6fvcVbRacRqPVsblvZRp4Oug7o9e0hMvr4JPR+pVYzSzwe+hH1/1dMTEyYU2jNdR1qWvo+EoG1MG9A7M+mcW10Gt8ceQLEm0K6FeBfXyKbKdnM7llCYLC45h35J6hoyrpUIYsEBlxFNMG/yd0X3UeR5vs7BhQlRJO1vDqLiyrA8EXoM1yqPUVOiQ/XvyRUSdG4WXvxYamGyhmW8zQ8ZUMrI5LHb6p8g1nnp1h9InRJJVsByU7wrHpVDa+Resyjiw5/pB7L94YOqqSzmTIAiGl/FlK2cfaOv2P/dfpJNP23ear7deo6mbPln6VccyVHR4eheV19ROZuu2BEm2J1cQy7Ogwll5bSpsibVhabym25ll7TX8lZbR0a8nI8iM58PgAM87PgCYzwcYVdn7O2HquWJiZMG7n9Sy/dpHyVxmyQGQU8RotgzdeYtGxB3SqWIAVXcthaW4KF1fD2jZgmR96HQLnCjyLfkaXfV04EnSE0RVG83XlrzE1VsMPlZTj4+GDj4cPG25v4NBzf2g2FyIeY3fxJ0Y0KMa5R+Hsu/7c0DGVdEQViFTyOiYRn+Xn2HP1GaMbuTOlpRcmAjjwNeweBAVrQM9fwMaFq6+u8uneTwmJDmF+nfl8VvyzLLm4mJL6hpYZioedBxNOT+B5nqJQsgOcmkvHgnG4O1jynd8t4jWqw1rRUwUiFTwJi6XNwtNcCY5kXqfS9PukMCIpHrZ2g1NzoGx36LQFzK3Z/2g/PX7pQXaT7KxtvJZqjtUMHV/JxEyNTZlRYwZJuiRGHR9FUt1vwMwCk33DmdCkOMGv41h+8pGhYyrphCoQKexyUAStFpwiPDaR9b0q0rRkfoh+Bauaws3dUH8yNP0BaWTMwisLGXF8BJ52nqxvsv6PhdYUJRW5WLkwrtI4Lr68yLJHu6DO1xB4giqxh6nvkZf5R+7zMirrroCq/EEViBR04OYLOi45Q45sJmzvX4VyrrZ/jFR6cUM/m7XKIBJ1GsacHMOCywtoVqgZS+svxcbcxtDxlSykWeFmNHRtyJKrSwgsUgscy8KB8Yyr54xGq2P6/juGjqikA/9YIIQQDkKIhUKI+UIIOyHERCHENSHEZiFEvrQM+Y5s6W6Y65ozgfRdc4FieS3Z/nkVCuXOCYGn9MtmaGKh214o3ozX8a/p/Wtv9jzcw6DSg5hSbQpmxmaGjq9kQaMqjMLc2Jwp/lORDadD9AsK3FxMz2qF2HYxmIDHrw0dUTGw5O4gVgE3gSDgCBAHNAZOAItSPVky0tMw19+GsY7fdYPa7nnY0KeSfseua1v1E+By5oFeB8GpLIGRgXT268z10Ot8X+N7+pTsozqjFYOxz27PoDKDOPvsLPuTQqFEezg9j8FlzHCwMmf8zutodWrYa1aWXIHIK6X8SUo5DcglpZwupQySUv4EuKRRvnQtIUnL0M2XWXTsAZ9VLMCizmWxMDWGE7NhW09wKq/fL9jGlYAXAXTe15k3iW9Y3mA5DQs2NHR8RaF90fZ42Hkw4/wM3nwyHIQRFscnMb6pBzefRbHu3GNDR1QMKLkC8efXVv+L92UJUfEauq88z67LTxnRoBiTW3phgk6/uc+hb8CrLfjsgOw27Hm4h96/9sYmmw3rmqyjVJ5Sho6vKAAYGxkzodIEwuLCmP9wB1QdAjd20Ng6kGpu9nz/yx1CoxMMHVMxkOT+0O8SQuQEkFKO++1JIYQbcDe1g6VnzyPjab/oDP6Pwpnd3psBtdwQiTH6nbsurIBqQ6H1UqSxGYuvLOarE1/hndubtY3X4mzpbOj4ivIXnvaetCvajo23N/LQqxlY5kfsH83Epu7Ea7R853fL0BEVA/nHAiGlnCCljH7H8/ellG1TN1b6de/FG1ovOEXw6zhWdi9P6zJOEP0SfJvC/QPQZDbUnYgGLRPPTGTe5Xk0LdSUxfUWq21BlXRrQOkBWJhYMOvqQqg3CZ5dxi14O31qFGL7xRD8rj0zdETFAN7bVCSEsEyLIBnB+cBw2i46g0Yn2dinkn5r0ND7sKwuvLwNHddD+Z5EJ0Yz8NBAtt/bTt+Sffmu2ndqpJKSrtma29KnZB+OBx/ntJ0juFSFQ98wpLI9pZxzMWrrVR6HxRg6ppLGki0QQghHYE8aZUnX9l9/Tudl57DLYcb2/lXwcrSGIH/9MNbEGP0w1mKNeBHzgm77u+H/zJ9JVSYxsPRANVJJyRA6Fe+EU04nvr/wPUkNp0F8FGbHpzCvU2mEgAHrL5KQpJbhyEqSmwfhCewHRqVdnA+T1vMg1p59zOfrAiiez4qt/avgbGsBt/eCbzPInks/UsmpLHdf3+Uzv88IehPE/DrzaVWkVZrkU5SUYGZsxrByw7gfcZ/tUbehQh+4sBKnuDvMal+K6yFRTN5zS634moUkdwdxBOgrpTybVmE+VFrNg5BSMvvXO4zbeZ2axfKwvndFbHOYwfnlsKkz5PWEngfArjDnnp2j676uSCnxbeRLFccqqZpNUVJD3QJ1KZu3LPMuzSOqygDIYQ97v6ReMTt6Vy/ImrOPmbz3Fjo1PyJLSK5AnAey7EfgJK2OMTuu8ePh+7Qv58QSn7dzHA59qx/K6lYPuv4MOezZ83AP/Q72wyGHA+uarMPd1t3Q8RXlPxFCMLrCaCITI1lwZx00mg4hAXD0O75qVJxuVVxZfvIRw7dcQaPVGTquksqSKxDNAWshxIy0CpNexGu09Ft7kQ3+QQys5cb0NiUxQQu7B8KJmVDaBzquR5pasOzaMr468RWl85TGt5EvDjkcDB1fUT6Ku6077Yq2Y8PtDdzJ7wVlusCJ2Rg9PMzXzTz4sl5Rtl8Koc/qC0TGaQwdV0lFyQ1z1Uop+wD/N9Q1M4uM1eCz/ByHbr9gUgtPhjcohtDE6uc4XFqr3ze6+U9ohWDKuSnMvTiXRgUbsajuIqzMrAwdX1FSxKDSg7Ays2Kq/1Rkg2mQ2x129EVEv2RQnSJMaeXFiXuhNJ57gguB4YaOq6SS9w5zlVJOSosg6cGzyDjaLT7NlaBI5ncqQ5fKrhATpu+Mvn8Qms6BWl8Rr01g2NFhbLqzie6e3ZlWfZoaxqpkKtbZrBlcZjABLwLY//Q4tFsJCdH6JWQ08XxW0YWt/atgbCRov/gMPx66p/olMqEPWjJDCGEjhPAUQhQSQmTKZTbuv3xDmwWneRoRz6oe5WlcIh+8DoQV9d8u1b0WynUnIj6C3r/2/n1r0GHlhmGUOX8kShbX2q01xW2L8/3573lt5QDN5kDgCVjXFuKjKOWci72Dq9HcOz+zD9xl+NYrJKl+iUwluWGu1kKIMUKIa8BZYDGwGXgshNgihKiVViFT26Unr2m76AyJWsmmvpWoUtgenl+D5fUhJhS67AL3JjyNfkqX/V24GXaTmZ/M5LPinxk6uqKkGmMjYyZVnUREQgTjTo1DluwArZfC49OwujnEhGFpbsoPHUoxvH5Rtl8Mod/ai2rL0kwkuY++W9Ev9V1dSllMSllNSllOSukMTANaCCF6pknKVHT0zks6LT2HdXZTtvWvjGd+a3h0AlY2BiMT6PELFKjEnfA7dPbrTGhsKIvrLaa+a31DR1eUVOdu686X5b7kePBx1t5aCyXb61cMeHkLVjSAl7cQQjCwdhEmtfDk4K0XdF95nrhEVSQyg+Q6qetJKddIKSPe8VqAlPILKeXy1I33bik1UW7npRB6+V6gUO4cbO1XBRe7HHBjJ6xtDVb59RPg8rjj/8yfbvu7YSSM8G3kSzmHcil0JYqS/nVy70RN55rMDpjNjbAbUKyhfqXi+AhYUgsCfEFKulR2ZU6HUpx5GMbUfWqBv8zgQ/sgHIUQVYQQNX57pHaw5KTERLnlJx/xxabLlHe1ZWOfSuS2zAbnl8GWbpC/NHTfB9ZO7A/cT7+D/chrkZe1jddSxKZIyl2IomQAQgi+rfItduZ2DD86nOcxz8GlCvQ7BQUqws+DYWsPSIylZWlHelUryOozjzly+6Whoysf6UMW65sOnALGASPePoancq5UtfDoA77dc5OGng6s7F4ey2wmcGQq7P0SijYAn51gYcu6W+sYeWwkJexLqDkOSpaWyzwXs2vO5nXCa7ru68qTqCdgmRc674A6E+DGDljbBuKjGN6gGO4OlozYepUwtZdEhibet66KEOIOUFJKme5+0+XKlZMXLlz41++79SyKrQHBjGlcHGN04Ddcv49Dqc7QbC7SyJgfL/3IsmvLqO1cm+k1pmNuYp4KV6AoGcuNsBv0P9AfYyNjFtdbTFGbovoXrm+D7X0grxd03s7tN6Y0/+kUnxTLzRKfsmrBynRGCBEgpXxvW/mHNDE9BEw/PlL6UTyfFeObemCsTdA3Kf22yU+LeSQJ+Pr01yy7toy2Rdsyu+ZsVRwU5S1PO09WNVyFEUZ029+NI0+O6F/waqPvvH51G1Y2wt0qiZENi3Hg5gt+vqr2ksio/vEOQgjxEyABR8AbOAT8fhchpRycFgGT81/vIACIj9LPjg48AQ2mQuXPiUuKY8SxERwLPkZ/7/709+6vPvkoyjuERIcw9MhQboXfonPxzgwrOwxTY1P9CMA1raBYI3RtfWky7xQxCUkcHPYJZiZqvlB6kRJ3EBeAAGA38C1w+u3Xvz0yruiXsKoJPDmjH9dd+XMiEyLp86t+w5RxFcfxeanPVXFQlH/gmNORtY3X0sm9E2tvrcVnnw+hcaFQsDrUHgu3dmN0bROjGhbjSXgsG/yfGDqy8h8kdwexBNgHHJRSvknTVB/oP99BHJgA/kuh/WooUo/nMc/pd6AfT948YXqN6dRzqZfyYRUlkzr0+BBfnfyK3Nlzs6z+MvJZ5IFVTeH5NWT/k3Ta8py7L95wbGQtcmYzMXRchZS5g1iOvmnJTwhxSAgxSgjhnWIJDan2eOh1EIrU42HkQ3z2+fAi9gWL6y1WxUFR/qU6LnVYUm8Jr+Nf03V/Vx5HB0OrRQCInZ8zukERwmISWXr8oYGTKv9WchPlzkkpJ0opqwPtgSfAl0KIS0KIFUKI9mmWMqUZm0JeT66+ukqXfV3QaDWsbLiS8g7lDZ1MUTKkUnlKsazBMuKS4ui2vxvPTM30e0k8PoV36B4al3Bg6YmHvHqT7gZDKsn4oF4jKWWYlHKDlLKLlLI0MB/I0DPGToWcotevvbAys2JNozVqkx9F+Ugedh6sbLCSWE0s406NQ+fdEZzKw9FpjKhdgIQkHXMP3TV0TOVf+JCJct8JIXL96WsboJWUckqqJktFfg/9GHhoIC5WLqxutBpnK2dDR1KUTMHNxo1RFUbh/9yfdbfXQ92J8OYpBR+up1OFAmzwD+L+yyy1xUyG9iF3EI3+vB6TlPI10Dj1IqW+HKY5KOdQjhUNVmCf3d7QcRQlU2nl1oqaTjWZEzCHB7ny6bfnPTGbL6rlIbupMdP33zZ0ROUDfUiBMBZCZPvtCyFEdiBbMsenuo9drO8T509YUm8JlmaWKZxMURQhBF9X+ZocpjkYc3IMmlpjID4CuyuL6PdJIQ7cfIH/I7ULXUbwIQViHXBICNHz7fLeBwDf1I2VvJRYrE/NcVCU1GOf3Z4JlSdwM+wmm6NugVdbOLOAnqUscLAyZ4rfLd63zI9ieB+y5eh0YDJQHHAHvpVSzkjtYIqiZGx1CtShgkMFllxdQmyNYaDTkP3U9wyrX5QrQRFsCQg2dETlPT507vtl4ARwHLiSenEURckshBAMLjOY8Phw1jw7CeV6wsXVtCkQS4WCtozfeZ3rIR+3p4uSuj5kFFMvwB9oBbQFzgoheqR2MEVRMj7v3N7Ucq7FqhuriKjUF0wtMD48iQWflcE2hxl91wQQHpNo6JjKP/iQO4gRQGkpZTcpZVegLDAqdWMpipJZDCo9iBhNDCse7oJqQ+D2HuzDL7HYpyyvohMYuP4iSVqdoWMq7/AhBSIM+PNaTG/ePqcoivJeRWyK0LRQU9bfXs+Lku3AMh/8Op6SjtZMaenF6QdhDN9yhcQkVSTSmw8pEPeBc0KIiUKIr4GzwF0hxDAhxLDUjacoSmbweanP0Uoti275Qs2vINgfrm2lXTlnRjQoxs7LT+m20p+oeI2hoyp/8iEF4gGwE/3eEAC7gEeA5duHoihKspwsnehQrAM77u3gYcEq+iU49nwBL28zoJYbs9p54/8onHYLz/A0Is7QcZW33rvlaHr2URsGKYqSpsLjw2m8vTGV81Xmh7IjYHENMLeG3ofB3JpT90PptyaAPFbZ2N6/KtYWmWojy3Tlo5f7FkLsTu6RsnEVRcnsbM1t6ebZjYNPDnIl4RW084XXgbCjH+h0VHWzZ2nXcjwJj6X/ugA0quPa4JJrYqoMOKGf/zATmPW3h6Ioyr/SxaMLduZ2zL4wG+lSBepPgTt+sH806HRUKmTHtNYlOf0gjHE7rqvZ1gaWXIFwAMYAXsBcoB4QKqU8JqU8lhbhFEXJXCxMLejn3Y+LLy9y6MkhqNgXKn0O/otha3fQxNOmrBODarux6UIQy08+MnTkLC25DYO0Usr9b+c+VEI/mumoEGJgmqVTFCXTaVO0De627kw4PYGg6GBo8B3Unww3d8KalhAbztC6RalbPA+zfr3Ls0jVaW0oyY5iEkJkE0K0BtYCA4AfgR1pEUxRlMzJ1MiU2TVnAzD0yFDitPFQZRC0XQEhAbCuHUbaeL5u5olWSqb6qeXBDSW5TurVwBmgDPCNlLK8lPJbKWVImqVTFCVTcrZ0Zlr1adx9fZfJZyfr+xq82vxRJHYNwNkmO31rFGL3laecD1TLgxtCcncQndFvKzoEOC2EiHr7eCOEiEqbeIqiZFY1nGrQ37s/ux/sZtm1ZfoiUbyZfhe669vg6DT61yxMPmtzvvn5Blqd6rBOa8n1QRhJKS3fPqz+9LCUUlqlZUhFUTKnvt59aVywMT9e+pFp/tPQ6rRQdQiU6gzHpmFxdxejG7lzPSSKLReCDB03y0muiSnn+978Icekho/dUU5RlPTBSBgxtfpUunh0Yf3t9Yw4PoIEXSI0/QGcKoDfCJoXs6C8qw0zfrlDRKxa+TUtJdfEtEsIMUsIUUMIkeO3J4UQhd7uLvcL0DD1I/6/lNhRTlGU9MFIGDGi/AhGlBvBgccHGHZ0GFojY2gyC2LDEcdnMqmFF5FxGmb8csfQcbOU5JqY6gCHgL7ADSFEpBAiDP2IJgegq5Rya9rEVBQls+vi2YWxFcdyPPg4PwT8APlKQhkfOLeI4qYv6VbFlQ3+T7gcFGHoqFlGssNcpZR+UsrPpJSuUkprKaWdlLKKlHKKlPJ5WoVUFCVr6OjekU/dP8X3pi/b722H2uPBJDv8Oo4v6hYhd85sjNt5TXVYp5EP3XJUURQlTYwsP5Iq+avw7ZlvOR/9GGoMh7v7sAw5wfimHlwPiWLduceGjpklqAKhKEq6YmJkwveffI+jpSMTTk0gsXxPsCkI+7+iqac91YvYM2XvLQIeq7kRqU0VCEVR0h0rMyvGVBhDcHQw6+9tg4ZT4dVthP9i5nQoRf5c2enpe4H7L6MNHTVTS26Ya4AQYq4QoqEQwjwtQymKolRxrEJ1x+osvrqYcJeKUKQBHJ2GnS4c3+4VMDESdF3hz8uoeENHzbSSu4OoiH7dpZrAMSGEnxBiiBCiaJokUxQlyxtebjhxSXEsuLwAGk0DrQYOjKeAnQUru1UgIjaRbivPE5OQZOiomVJyw1yTpJRHpZSjpZQVgV7AG2Dy27uLBWmWUlGULKlQrkK0K9qOrXe38sBI6mdZX9sCgScp4WTN/M/KcOfFGwZvuKRGNqWCD+6DkFI+lVKukFK2B8oD61IvlqIoit7npT7HwsSC789/j6z6BVgXgN2DIDacmsXyMLG5J4duv2Ty3puGjprp/KdOaimlTkp5KqXDKIqi/J2NuQ39vPtx6ukpjr7wh9ZLIDIENnwKmjh8KrnQo2pBVp4KZPWZQEPHzVTUKCZFUdK9T4t/SmHrwsw4P4MEpzLQahEEnYUdfUGnY2yT4tRxz8Okn29y65labDqlqAKhKEq6Z2pkyqgKowiODsb3hi94tX67C90uODAeYyPBzHbeWGc3ZdS2qyRpdYaOnCm8b0c5cyFE27fDXbcIIVYLIUYKITzTKqCiKApA5fyVqVugLsuuLeN5zHOoPBAq9IEz8yBgFTY5zPimhSdXgyNZcUrtZZ0SkpsH8Q1wCqgMnAMWA5uBJGCaEOKAEKJkmqRUFEUBhpcfjk7q+O7cd0iABlOhcB3Y+yU8OkGTEvmoWzwvsw/cJTA0xtBxM7zk7iD8pZRlpZRfSinXSykPSin3SClnSymbAZ8BZmmUU1EUBcecjgwqPYgjQUeYd3keGJtAu5VgWxg2+yDCHzK5pRemRkaM3n4VnRr6+lGSmwexF0AIUfDvrwkhykspX0opL6RmOEVRlL/r4tGF1kVas+TqEnbc2wHm1tBpo/7FDZ/iYAHjmhbn7MNwFh57YNiwGdyHdFJvE0I4/vaFEOITYEXqRVIURflnQgjGVRpH5XyVmXRmEmeengHbQtBmGYTegROzaF/OmaYl8zHr1zv4P1KL+v1XH1Ig+gI7hRAOQojGwI9A49SNpSiK8s9MjUyZVXMWrtaufHn0Sx5FPgK3ulCyA5z8AfHqDlNbl8DFLgeDNlwkNDrB0JEzpPcWCCnleWAw8CswEagrpVS7hyuKYlCWZpbMrzMfU2NTBh8eTFRiFDT4DrLlhJ+HYGlmzLxOpXkdq2HopstqKY7/ILlRTD8LIXYLIXYDXwEWQAKw/O1ziqIoBpU/Z35m15xN8JtgRh4biTa7DdSfop9Ed3EVnvmtmdjMkxP3Qvnm5xtIqYrEv2GSzGsz0yyFoijKf1Q2b1nGVhrLN2e+4YeAHxhe7ku4sgEOTAS3unSqWIBHodEsPfGIfNbZ6V+zsKEjZxjJFYjj8j3lVggh3neMoihKamtbtC23w2/je9OXWgVqUbb5T7CoOmzrBd38+KpRcZ5HJTB9/23yWZvTsrTj+0+qJNsHcUQIMUgIUeDPTwohzIQQtYUQvkDX1I2nKIryYYaVHUb+HPmZfHYyGmsnaDYHgs7B0akYGQlmtitJpUK2jNh6hYDHrw0dN0NIrkA0BLTABiHEMyHETSHEI+Ae8CkwR0q5Kg0yKoqivJeFqQVjKo7hfsR9Vt9cDSXaQqnOcGIWPDxGNhNjFvuUw8HanMEbLhERm2joyOlechPl4qWUC6SUVYECQB2gtJTSRUrZW0p5KSWDCCGKCyEWCSG2CiH6p+S5FUXJGj5x/oTazrVZdGURIdEh0HgG2LnB9j4Q9Qzr7KbM+7QML9/EM3zLVdVp/R7JjWIyF0J8IYSYB3QHXkkpI/7NyYUQK4QQL4UQ1//2fEMhxB0hxH0hxGgAKeUtKWU/oD1Q9d9fiqIoCoyuMBohBFPPTUWaWkC7VZDwBjZ0gIRovJ1zMbpRcQ7eesHyk2pRv+Qk18TkC5QDrqGfGDfrP5x/Ffqmqt8JIYyB+UAjwAP4VAjh8fa15sBewO8/fC9FURTy5czHgFIDOBZ8jJ33d4KDl369pufX9J3WOi09qrpSzyMv0/bd5kKgmmn9T5IrEB5Sys5SysVAW6D6vz25lPI48PeffgXgvpTyoZQyEdgItHh7/G4pZSP0CwEqiqL8J52Ld6a8Q3mm+k/lcdRjKNoAGs2Au/vglzEIIZjZ1hsnm+z0XRNAUHisoSOnS8kVCM1v/yGlTErB7+kI/HkmdjDgKISoKYT4UQixmGTuIIQQfYQQF4QQF169epWCsRRFySyMjYz5rtp3mBqZMvr4aDQ6DVToDZUG/K+9e4+zqd7/OP767MtcMCmGJCGHJLcwhCgJcVBJF3GKonR+bvWrHCXdE0VRyq0G1SmdUr+ocLopDskQuZQ4Ti7dXHM3GN/fH2uXmTkbM8yePXvm/Xw85jEz3/Vda31mvo8971lr7fVdsHAcLJtKyWJBXurekIMZR+g1JY3dBw6deMNFzPECoq6Z7Qp97Abq/P61meX5M/2cc3Occ/2dc72dcy8cp98E51yKcy6lTJkyeV2GiBQS5YqX4+GmD7Ni2wpeXPqi19jmMajY1Ht+xLZ/U7VsCcZ2a8DaLXsYMFXTcWR3vHcx+Z1zp4U+kpxzgUxfn3YK+/wROCfT9xVCbSIieap1pdZcU+0aXl7+Mp9v/Bx8frhmgvd5Wi/IOESzask8fGVNPv1uM49qOo4sovFM6kVANTM718zigC6A5nYSkYgY1GgQNUrXYOAXA/l+x/dw+jnQ8Tn4aQl89gQANzWuRK9m5zJlwXomzl0X5YoLjogGhJm9ASwAqpvZJjPrGbqe0ReYDXwL/MM5tzKX2+1oZhN27tyZ90WLSKGSGEjkucueo3iwOP0+6ce2/dug5tVQ/2aYNwrWfAzA/X+uQfvaZzH0w++YseynKFddMFgsH06lpKS4tDQ91E5ETmzl1pV0n9WdGqVqkHpFKsGMg/BSa9i1CW77DEr/iQOHMrjp5YUs27iTV3s24qIqpaNddkSY2WLnXMqJ+kXjFJOISL6rmVyTx5s9ztItS73nWccVhxtfB/PB1K6QvpuEoJ+JN6dQoVQivV9bzLote6JddlQpIESkyGhbuS2dq3Vm0opJpP2SBmdU9u603roG3ukNR45werE4JvdohN+MWycvYvveojtnkwJCRIqUgQ0HUiGpAoPnDWb3wd1QpQVc8QSs/uCPi9YVSxdjws0p/LTzALe/ksaBQxlRrTlaYjIgdJFaRE5WsWAxhjYbyi/7fmHYV8O8xovugHo3wdwRsPQNABpUOoNnrq9L2vod3Dl1KYczjkSx6uiIyYBwzs1wzt1esmTJaJciIjHowrIXclvt25j+7+m8v+59MIP2z0Dl5jC9H6yfD3CxVagAABDaSURBVECHOuUZ0uECZq38hb9NW86RInYjXUwGhIjIqepdtzcNzmzAI/MfYfX21RCIgxte9a5LTO0G2737IXo2O5c7W1Vj2pJNPPr+qiJ1I50CQkSKpKAvyIhLR5AUl8Rdc+5i18FdkHgGdH0TcPB6FzjgzSo04PJq9Gp2LpPn/8AzH30f3cLzkQJCRIqs5MRkRrYYyc97fmbw3MEccUeg9J/g+ldg+7//mB7czBjcvgZdGp7D85+uZeIXReNuawWEiBRp9crW456G9zBn0xzGLxvvNZ57CbQbDmtmw8cPA2BmPNGpNu1rn8UTH37Lm4s2RK/ofBKIdgEnw8w6Ah2rVq0a7VJEpBDoen5XVm1bxYvLXqTaGdVoVakVNOwFm7+F+c9BmepQ7y/4fcazN1zInvTD3PfOcorHB+hQp3y0y4+YmDyC0LuYRCQvmRkPNnmQOsl1uH/e/d5Fa4C2w6DKZTC9P3w/G4C4gI9xf2lASqVSDJi6lJnLf45i5ZEVkwEhIpLX4v3xjLpsFEnBJAZ8NoAdB3aAP+i9s6lcbfhHd9i4CIDEOD+ptzSk3jmn0++Nr5m1onCGhAJCRCSkTLEyjG45mi37tnDv5/dy+MhhiE+Cbm9DUjl4/TrY4h1dlIgPMPnWRtQ953T6vl44Q0IBISKSSa3kWgxpMoSFvyxk1OJRXmOJMnDTO+ALwmudYZcXBiXiA0y+pSF1KpSkz+tf82EhO92kgBARyebqqlfTpXoXpqyawsz/zPQaS1WBbm/B/h3w92vhgDfVT1JCkCm3NuLC0OmmwvQsCQWEiEgYAxsOpF7Zejw0/6GjF63LX+hdk9jynXe39eF04GhI1K94OgOmfs17SwvHU5RjMiA0WZ+IRFrQH+SZFs+QFJdEn0/6sHnfZm/Bn1rCVS/CD3Phndsg4zDw++mmRjSsXIo731xaKO6TiMmA0NtcRSQ/JCcm88LlL7D74G76fNKHvYf2egvq3gBXDIVV78GM/nDEm+m1eCgkmlcrw9+mLWfSv/4TxepPXUwGhIhIfjm/1PmMuHQEa3asOfrOJoAmfaDF/bD07zBzIIQm8UuM8zPx5gZcUfNMHpmxihc+WxuzE/wpIERETqB5heYMbjyYuT/O5cmFTx79g3/pQGjaHxZNhI8e/CMk4gN+Xuhan071zubp2asZNvO7mAyJmJxqQ0Qkv1133nVs2r2J1BWpVDytIt1rdveeI9H6UTi415uSI1gMLrsPgIDfx8jr6pKUEGD8F+v4bd8hhl5TG7/PovyT5JwCQkQkhwbUH8DG3RsZmTaSCiUqcHmly72Q+PMI7x1Nnw/znivR/G4AfD7jkStrUjIxyPOfrmV3+iGeveFC4gP+KP8kOaNTTCIiOeQzH0ObDaV2cm0GzR3E8i3LQwt8cOVzUPs6+ORRWPDCH+uYGXe3qc4D7Wvw4fJfuGXSInYfOBSlnyB3FBAiIrmQEEhgdMvRlE4sTZ9P+vDDzh+8BT4/XD0OalwJs+/PEhIAvZpXYeR1dVn4n+3cOPFLtu5Jz//icykmA0L3QYhINCUnJjOu1TgAen/U++g9Ev4AXJt6NCTmj8myXucGFZh4cwPWbt7DtWPns2HbvvwuPVdiMiB0H4SIRFvlkpUZ22osv6X/Ru+PerMzPfQPqz/ohcQFV8E/B8P857Os1/L8M/l7r4vYse8Q14ydz4ofC+4/ujEZECIiBUHN5JqMbjma9bvW0+/Tfuw/vN9b4A9C55ehZif45wMwd2SW9RpUKsW0vzYhzm90mfAl89dujUL1J6aAEBE5BY3Pasyw5sNYtmUZ/zvnfzmUEboA7Q/CNS8dvXA9Z9gf90kAVC2bxLT/aUr50xPoPumrAjnJnwJCROQUtanchiGNhzDvx3kMnjeYjCMZ3gJ/ADqNhwu7wZwnvaDIFBJnlUzkrd5NqXfOGfSf+nWBm5pD90GIiOSBa8+7lp3pOxm1ZBRJcUk80PgBzMx7d9OVY7wjinnPwKF9cMWT3ltjgZLFgrzSsxEDpn7NIzNW8cuuAwxqe763bpQpIERE8kjP2j3ZdXAXqStSSQgkcE/KPaGQ8EGHURBXAhaMgfQ93n0TPu+GuYSgnxe7NeCh6SsY//k6tuxOZ3jnOgT90T3Jo4AQEclDd9a/k/2H9/PKqldIDCTSt15fb4EZtHncC4nPh8GhvdBpgnfnNeD3GY9dVYtypyUw4p/fs3XPQcZ2q0/x+Oj9mVZAiIjkITNjUKNBHMw4yPhvxpMQSKBX7V6/L/TmaoorDh8N8Y4krn8F4or9sW7fltUokxTPfe8sp+vEL0nt0ZDSJeKj8rPoIrWISB7zmY8hjYfQvkp7Ri8ZzeQVk7N2uLg/dBwNaz/2nnF9IOu9EDc0rMiEm1JY/eturh23gI3bo3NDXUwGhO6kFpGCzu/z8/jFj9O2cltGLh7JKytfydqhQQ+49mXY9BVM6Qh7s94L0eoC74a67XsP0nnsfL79eVf+FR8SkwGhO6lFJBYEfAGebP4krSu15um0p3lt1WtZO9TqDF3egC2rIbUt7NyUZXGDSqV4644m+My4fvwCFq7blo/Vx2hAiIjEioAvwPBLhnN5xcsZvmg4r656NWuH89rATe/Cnl/h5Stg65qsi8/0bqgrmxTPzalf8fGqX/OtdgWEiEiEBX1Bnr70aVpVbMVTi55iysopWTtUago93ofDB7wjiZ+WZll89umJvHVHU6qXS6L3a4uZtjjrkUakKCBERPJB0BfkqUufok2lNoxIG0HqitSsHc6qC7fOhmAiTO4AP8zLsrhU8Thev60xF51birvfWpYvIaGAEBHJJ0FfkOGXDKdd5XY8u/hZxi8bn7VDclUvJE4r7727afXMLItLxAeYdEtDujepxMVVkyNerwJCRCQfBXwBhjYfSscqHRmzdAzPf/08LtP8TJQ8G26ZCWUvgKndYNmbWdaPD/h55KpalCuZEPlaI74HERHJIuAL8NjFjxH0B5nwzQQOZRzirgZ3HZ1/qXhp6D4dpnaFd2/37pO46Pb8rzPf9ygiIvh9fh5q8hBBX5BJKydxIOMAgxoNwmehEzvxSdD1LZjWE2beCwd+g0vu9e7GzicKCBGRKPGZj8EXDSbBn8CUVVNIz0jnwcYP4g9N4kcwAa6bAtP7wmdPeEcSbR7Pt5BQQIiIRJGZcXfK3SQEEhj/zXgOHD7AE82eIOAL/Xn2B+CqFyH+tNBMsLu8mWF/D5EIUkCIiESZmdG3Xl8SAgmMXjKa9Ix0nrrkKeL83kyv+HzQbjgknAZfPO1N8tdp/B8zwUZKTL6LSXMxiUhh1Kt2LwY1GsQnGz6h/6f9jz7jGrzTSi0fgNaPwsp3YNFLEa/Hsry9KsakpKS4tLS0aJchIpKn3l3zLg/Nf4j6Z9ZnTMsxlIgrkbXDmo+hSgvv9NNJMLPFzrmUE/WLySMIEZHCrFO1Tgy/ZDhLNy/l9o9uZ2d6trMl1VqddDjkhgJCRKQAanduO55p8Qzfbf+OW2ffyrb9+TuTKyggREQKrJYVWzLm8jFs2LWBHrN68Ove/JvJFRQQIiIFWtPyTRnXehxb9m+hx6we/Ljnx3zbtwJCRKSAa3BmA15q8xK7Du6i+8zurN+1Pl/2q4AQEYkBtZJrkXpFKgczDtJjVg/W7lgb8X0qIEREYkT1UtWZ3HYyQV+QDbs3RHx/upNaRCSGVDm9CjM6zSDeHx/xfekIQkQkxuRHOIACQkREjkEBISIiYSkgREQkLAWEiIiEpYAQEZGwFBAiIhKWAkJERMKKyRvlzKwj0BHYZ2bfZltcEgj3qLns7cnA1shUeELHqjHS28lp/xP1O97ynP7+j9UWrXGJ1pjkZp2THZdTbddr5eT7FdTXSqUc9XLOxewHMCEnbeHagbSCVHd+bCen/U/U73jLc/r7P05bVMYlWmOSH+Nyqu16reT9mOR2XKL1Won1U0wzcth2vPZoyKtacrudnPY/Ub/jLc/N719jkrt1TnZc8qo9GvRaydl+IiKmn0l9KswszeXgmaySvzQuBY/GpGDKj3GJ9SOIUzEh2gVIWBqXgkdjUjBFfFyK7BGEiIgcX1E+ghARkeNQQIiISFgKCBERCUsBEYaZ1TCzcWb2tpn9Ndr1iMfMrjaziWb2ppm1iXY9AmZWxcxeNrO3o11LUWZmxc1sSuj10S2vtlvoAsLMUs1ss5mtyNbe1sxWm9laMxt0vG045751zt0BXA9cHMl6i4o8Gpf/c87dBtwB3BDJeouCPBqTdc65npGttGjK5fhcA7wden1cmVc1FLqAACYDbTM3mJkfeAFoB1wA3GhmF5hZbTN7P9tH2dA6VwIfAB/mb/mF1mTyYFxCHgitJ6dmMnk3JpL3JpPD8QEqABtD3TLyqoCYnIvpeJxzX5hZ5WzNjYC1zrl1AGY2FbjKOfck0OEY25kOTDezD4DXI1dx0ZAX42JmBgwDZjrnlkS24sIvr14rEhm5GR9gE15ILCUP//EvjEcQ4ZzN0XQF75d59rE6m1kLM3vOzMajI4hIytW4AP2AVsC1ZnZHJAsrwnL7WiltZuOAemZ2X6SLk2OOzztAZzMbSx5Oy1HojiDygnNuDjAnymVINs6554Dnol2HHOWc24Z3TUiiyDm3F7glr7dbVI4gfgTOyfR9hVCbRJfGpeDRmBRs+To+RSUgFgHVzOxcM4sDugDTo1yTaFwKIo1JwZav41PoAsLM3gAWANXNbJOZ9XTOHQb6ArOBb4F/OOdWRrPOokbjUvBoTAq2gjA+mqxPRETCKnRHECIikjcUECIiEpYCQkREwlJAiIhIWAoIEREJSwEhIiJhKSCkSDGzytmnTz5O35lmViFM+ygzuyT09Q9mlnyM9YNmtiT0tTOz1zItC5jZFjN7/yR/jhFm1vJk1hXJKQWESBhmlgiUds5tytZeGmjsnPsiB5tpBvwr9PVeoFZouwCtObUpEp4HjvusBpFTpYCQIiv0NLSvzaxhmMUtCD9hY2dgVra2fma2xMyWm9n5mdrbAjMzff8h0D709Y3AG5lqedjMXjWzBWa2xsxuy7Tsb6FtLzOzYQDOufVAaTMrl7OfViT3FBBSJJlZdWAa0MM5tyhMl3b8dxCA94TBxdnatjrn6gNjgXsytV9G1pCZCnQxswSgDrAw23bqAC2BJsCDZlbezNrhzfd/kXOuLvBUpv5L0BMPJYIUEFIUlQHeA7o555Ydo8/FwLww7WcBW7K1vRP6vBioDGBmZwPbnXP7fu/knPsmtPxGwj9n5D3n3H7n3FbgM7yHw7QCJv2+Hefc9kz9NwPlj1G/yClTQEhRtBPYgHeN4L+YWRVgo3PuYJjF+4GEbG3poc8ZHH3GSlu8CdWymw6MINPppUyyT4x2oonSEkL1iESEAkKKooNAJ+BmM+saZvmxTi+BN4Nm1RzsI/v1h9+lAo8455aHWXaVmSWELoS3wJva+SPgFjMrBmBmpTL1Pw/I0TuyRE6GAkKKpNATuDoAd5nZldkWt+XYAfEB3h/vYwo9WL6qc+67MPvdFHoyXjjf4J1a+hJ4zDn3k3NuFt5RR5qZLSV0jcPMgnhBlXa8WkROhab7FsnEzOKBfznnUo7TZx7QwTn32zGWNwP+4pzL8aM4zexhYI9zbkQO+3cC6jvnhuR0HyK5pWdSi2TinEsHjhkOIXcDFYGwAeGcm0f4C9x5KQCMjPA+pIjTEYSIiISlaxAiIhKWAkJERMJSQIiISFgKCBERCUsBISIiYSkgREQkrP8He34mEDsXKH8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Write your code here\n",
    "k = np.logspace(-3, 0, 100)\n",
    "for z in [0.0, 0.1, 0.2]:\n",
    "    p = P.P(z, k)\n",
    "    plt.loglog(k, p, label=f'z={z}')\n",
    "plt.legend()\n",
    "plt.xlabel(\"k / (h/Mpc)\")\n",
    "plt.ylabel(\"P(k) / (Mpc/h)^3\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More Activities\n",
    "\n",
    "If you complete the activities above, try one of these activities:\n",
    "- Plot some other CAMB outputs, like the CMB power spectra\n",
    "- Compare the angular diameter distance calculated by CAMB and by Astropy.  They should agree to 10^{-6} if all the parameters are the same.  Beware of neutrinos!\n",
    "- Dragan described the $\\sigma_8$ parameter, an integral over P(k).  Calculate it, and compare it to the camb internal version"
   ]
  },
  {
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}