CAMB-matter-power.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: Alkistis Pourtsidou, QMUL\n",
    "\n",
    "Using part of http://camb.readthedocs.io/en/latest/CAMBdemo.html \n",
    "\n",
    "To run this Jupyter notebook you need to have CAMB and the CAMB python package \n",
    "installed. In order to install the CAMB python package on your computer follow\n",
    "the instructions in http://camb.readthedocs.io/en/latest/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Import numpy, matplotlib, and camb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, platform, os\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import camb\n",
    "from camb import model, initialpower"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import scipy and interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy\n",
    "from scipy.interpolate import interp1d\n",
    "from __future__ import division"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pi=np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Set plotting and printing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "font = {'size'   : 16, 'family':'STIXGeneral'}\n",
    "axislabelfontsize='x-large'\n",
    "plt.rc('font', **font)\n",
    "plt.rcParams['text.usetex'] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Cosmological parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "c=3e5\n",
    "hubble=0.678\n",
    "omegab=0.022*pow(hubble,-2)\n",
    "omegac=0.119*pow(hubble,-2)\n",
    "om0=omegac+omegab\n",
    "H00=100*hubble\n",
    "Ass=2.14e-9\n",
    "nss = 0.968\n",
    "\n",
    "gamma=0.545"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set up the fiducial cosmology CAMB will use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set up the fiducial cosmology\n",
    "pars = camb.CAMBparams()\n",
    "#Set cosmology\n",
    "pars.set_cosmology(H0=H00, ombh2=omegab*pow(hubble,2), omch2=omegac*pow(hubble,2),omk=0,mnu=0)\n",
    "pars.set_dark_energy() #LCDM (default)\n",
    "pars.InitPower.set_params(ns=nss, r=0, As=Ass)\n",
    "pars.set_for_lmax(2500, lens_potential_accuracy=0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculate results for these parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "results = camb.get_results(pars)\n",
    "#print pars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get background quantities, e.g. the angular diameter distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEiCAYAAADAnxR7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xlwm/d95/H3j7d4gqDug6Kow7Jk\nyxJF2YntXDXVtEmaZBvJbtM0bZqYctuZzra7kZI27bQ7O+uVOpu2M7tpJGe7PZLdtaU01yZNItmJ\nE9utLIqSdVgnKYmiLkoiwUO8yd/+8TygIRAkCBLAAwKf1wxGwvN7AHzx2Hq++N3GWouIiEiiZHkd\ngIiIpDclGhERSSglGhERSSglGhERSSglGpmVjDE+Y0yNR59d7dVnzxa6RhJKiUaiMsbs9zqGUMaY\namA/8LIHn13jfvb+sOP1xpijyY4nFekaSTglGpmUMcYH1BhjtnkdS5C1thnY7tFnN07w2Q3Ai0kO\nZ1Je1SjieY1UK0oPSjQSTTWwG3jG60BCWWsDXscQylrbaK3d43UcYeq8DiDUNK9RSn0HmR4lGpmK\nl4CUqdHIlKXUj4NpSofvkPGUaGRCbl9Is1t7OJSM5rNInfxedvxPlTGmxhiTEr++3eu1F6c2mjJi\nuUap+h1kenK8DkBSms/tDwGnY3cHcCD0BDcBvAD4gM2809SxFdgd8vrg+buBJvf8ZsAfUtyA00xX\nB5iQ97/v2GRCkqHfjX9PSFkwVnD6EGrc77R9sqY49+ZY48YLEAgrrwa+6MZYPo14wq/dM8AuoJ3o\n17PePQ9gC/A88LT73Odeb4CD1tpDk73OWhtI9jVyr0/wXB+w1Vq7Y4rfYTrXdtx1dN+nOiSOQ2Hl\nEa/VRNdCIrDW6qFHxAdQE/J3H2Bx/kGHn1cNdAB1IcfqcG4MoeftBnaGPD8Y/Izg+wY/J8L72wif\n2xH2fFtofEB96OcFvxNOoqtzn+8Fqie5BtuA/WHHdgJNYcd804xnwms3heu5E9gW9l5HQ993gu80\n4euSeY3c120LO+dgWFwTfYcZXdtosU/1WukxtYfnAeiRmg/3H5Qv7NhBoD7CueOS0AQ33qawf/Q7\ngb2R3ivaMfd4+PvvD32/SDeFSLFGuQ4d4TdZ932bEhXPBMfuu+lGur7u8aNuooh4k472umReI/e/\nf31YeX3Ye0yUaGZybTsmij30NVO5VnpM7aGmM5lIHbDSmHGtVduBfZFeYMc3J/jCnoeXV0w7usie\nhbFmHB9uk0qE8wIRYh3HbX4JbT5MVDyRrl34sfaw4jqgPUK/2SEmN9XXJeMaHQCOGmO24ySOl6y1\nEf/fimAm19Y3UezWaT7c7P65jeldYwmjRCMTaQ7/R+/OqekwxvimchOKYC9Oogr+Q60jvvNh/Dg3\nrF3W2gNuv8DuCOeF37QnEtpun8h4pqMaJxkcCDt+AMb6RGJ+XYiEXyNrbbMxZgVOs9cOYK8xZped\n2hDoeFzbiLFbZx7QWPkUrpVEoVFnMk5wtFn4cTe5NPJOR22smoH97gzxepwO5vs+Z5oJLOggzo1n\n3C9ON0nGqpkJfiV7FE+omGILGe010+80ozhCGWNqrLUBa+0ea+1mnEECOyZKkmEj1uJxbaPFHu9r\nlbGUaCSSmkmaQl5k+rWQGmvtIWvtPvcx1eaWqENc3ZtTddiNJ/R1MQ89dn/ZBiLc+KLefBIRT5hD\nQHV4bCFDwcNrJNVTfF1MZnKNgGdCk4L7I2N3SKwRv0O8rm1I7OO+t/sZcb1WmUyJRiLZMknZIaBu\nmr/IK6b4D7Q57P1rw09wy0PPaQ85HhT6d3/Yn1O1y32Eeib8fWYQTyzCb8o7GN9c9LR1ZuAHInx+\n1NdNM7bpXiNwhjyH8uEMc2ei70B8r+2u8BiMMdustc1TvFYyFV6PRtAjdR44Qz2bcEbdbItQHlzM\n0uKMvKnHGeEUPBacYBd+zBfyeosz0qfD/az9jB+xVIc7rNR9VId8ZnVYHKHvv819vs19Dx/OTWK3\nG1MNTpNL8HVTGjnkvl+9+57BP637Xr7weMJeFy2eqNcuQty+sM/YHXKtfGHXMVgWfo0jvi6Z1yjk\n/59tIY+6sPeN+B3idW0nuBZTulZ6TP1h3AspklBu88MOQia7uceCEzI3W02CE0lLnow6M8bsdP+6\nBThiw0aZuOVjs8bt+NFPMyoXT+wAXgxNJtbpo2k2xmzB+VWq0TwiaSjpfTTGmL3WGWWyx1q7HadD\ncGdI+W6cobUH3ASxMnQc+0zLxTMHmXiBxGDHq4ikoaQ2nbmdd1+01u4KOVaPs/ZQufu8I/h393kd\nzjDGrfEoF++ErIcVwOnQDU6yO2CnPylSRFJcspvO/MBOt1YTemMJnakbbmxhwZmWi7esMxxVNReR\nDJPUpjM3uWwOSzJbeefm42f82PnQIY4zLRcRkSRL+mAAGzL+3L351+Es4w3vrFkUKpg4/HEoH+uI\ndpvs6gGKioo2r127NtavIiKS0Y4ePXrHWjsv2nler3W2H3gqpIYTaXhrMHG0x6F8jDtQYB9AbW2t\nbWhoiCFsERExxlyZynmerQzgjg7bbe+fYdvO+JnDPhibJTzTchERSTJPEo073Hhst7xgJ76bdMIT\ngh+3D2em5SIiknxezKOpw7n5N7iL01Vz//yKfWHzXrbiLBcRr3IREUkiL+bRdEQoOuBO3gyeF5zZ\nH9wPYqKZ/9MqD6c+GhGR2Bljjlprxy16Gy6pgwHcfpJxWzZGOG/SjY9mWi4iIsmjbQJERCShlGhE\nRCShlGhERCShlGhERCShlGhERCShlGhERCShlGhERCShlGhERCShlGhERCShlGhERCShlGhERCSh\nlGhERCShlGhERCShlGhERCShkrpNgEgqGR21dA8M09k7RFf/EPcGhukdGqFvcIShkVGGRizDI6P3\nvSY7y5CbnUVOtiE/J5s5udkU5GZRlJ9DcX4OJQU5lBTkkp0VdTcMkYyhRCNpq2dgmObbPVy+20tr\nRy+tHX1cD/TR1jVAW/cA7fcGGE3Qvn8lBTn4CnPxF+bhL8qjojifucX5zC/JZ35pPgtKC1hUVsCC\n0gJys9WwIOlNiUZmPWstLe29nGjt5PT1Lt6+0cX5m93c7Oq/77yKojwW+Zwb/IalZcwtzsdXmIuv\nMI/SghyK8nMozMtmTl42udlZ5GVnkZ1lMCb4OTAyahkedWo6A8Oj9A2N0Ds4Qu/AMN0Dw/T0D9PV\nP0Sgd4hA7yDtvUO0dQ9w5kY3d+8NMDRyf2YzBuYV57OkfA6LfXNYWj6HSn/h2GOJbw45SkQyyynR\nyKxjreXszW5ev3iHNy+109jSwZ2eQQBysw2r55fw+MoKVs4vZuW8YlbMLWJp+RyK8r3933101BLo\nG6Ktu59bXQPc7OzjeqCfG519XAv0cfpaJz8+ffO+ZJSTZVjmL2R5RSEr5hZRPa+Y6rlFrJpfzPyS\nfIxRE52kPiUamRV6BoZ59dxtDp25xc8v3B5LLJX+Qt67eh6bq8p5ZKmPNQtKyMtJzRpAVpbBX+Q0\npa1dGPmckVHLra5+Wtp7abnby+W797hyt5fmO/c43NxO39DI2Lkl+TlUzy9m9fxi1iwoZvX8EtYs\nLGFxWYESkKQUJRpJWb2Dw/z49C2++9Z1Xrtwh8GRUcoLc3nvmnk8uWouT6yay2LfHK/DjKvsLMNi\nn9OM9q7qivvKrLXc7Oqn+fY9mm73cLHNebx6/jYHjraOnVeSn8OahSU8sLCEBxeW8OCiUh5YWEJJ\nQW6yv44IAMbaBPWGziK1tbW2oaHB6zAE52Z65HIH/+fNFn50+ia9gyMsLivglx9exC+uW8Dm5eXq\ns4gg0DvI+Vs9nLvVzfmb3Zy72c3Zm1109Q+PnVPpL+TBRSWsW1TG+sWlrF9SysJS1X5k+owxR621\ntdHOU41GUsK9gWH2N1zlG4dbuNDWQ0lBDh/buISPb1zMlio/WRouPClfYR6PrvDz6Ar/2DFrLTc6\n+zlzo8t9dPP2jS5+/PYtgr8vK4ryWL+kjIcWl/LQkjIeXlLG0vI5Sj4SV6rRoBqNl253D/APb1zm\nn/7tCp19QzyyzMdvPFbJr2xYzJy8bK/DS0v3BoY5c6OL09e7OH29k1PXujh/q5thd6y3rzCXh5eU\nsWFpGQ8v8bFhaRmL1O8jEahGIyntTs8AX/lJE18/fIWhkVE+uG4h9e+rpqay3OvQ0l5Rfg61VX5q\nq96p/fQPjXDuZjcnr3Vy6lonJ1o7+eqrzYy4yWducT6PLC3jkWU+57G0DF9hnldfQWYZJRpJqu7+\nIfa+2szfvX6J/qERfrVmKb/3/pVUzyv2OrSMVpCbPZZEgvqHRnj7RhcnrgY40drJW60BXjnXNtbs\nVlVRyMZlPudRWc66RaUpO+JPvKVEI0kxOmr5ZmMru394jjs9A3xkwyL+cOsaVirBpKyC3GxqKsvv\nq2V29Q9xqrWT460BjrcEeKPpLt8+fh2AvOws1i8pZdOycjZV+thU6WOJT/09oj4aQH00iXbqWid/\n8u1TvHU1QE2ljz//6Ho2LPVFf6GkvOCAg+NXA86jJcCJawH6h5w14uaX5LtJx0lYG5aWUZCrvrd0\noT4a8Vz/0Ah/8/IF9v2smfLCPL789CN8fOMSjSBLI8a8M+/nQw8vAmBoZJRzN7tpbOmg8UoHx64G\n+NHpW4Cz0sG6xaXUVJazebnzSLe5UDKeajSoRpMIx1o6+A/736L59j22b17Klz68jrJCTRjMVHd7\nBjjWEqCxpYOjVzo40do5tsrBwtICNi8vp8ZNPOrrmT1UoxFPjI5avvqzJr784/MsKC3gH3/nUd67\nZp7XYYnHKorzqVu3gLp1CwCn1nP2RvdY4jl6pYPvn7wBQH5OFo8s8zk1HrfmU16kEW6zmWo0qEYT\nL23d/fzRi2/x2sU7fPjhRfyXX32YsjmqxcjU3OzsH0s6R1s6OH2tc2xuz8p5RWxeXk7tcj+bq8qp\nnlukQQYpYKo1GiUalGji4VhLB899/SidfUP8+a+s55kty3QjkBnpHxrhrasBjrZ0cPSyk3wCvUMA\n+IvyqKksp7aqnC1V5Ty0pIz8HA0ySDY1nUnS7G+4yp986xQLyvL59u8/wdqFpV6HJGmgIDebx6or\neMxdXHR01NJ8p4eGyx00uDWfQ2ecQQZ5OVlsWFLG5qpytiz3q7ktxahGg2o00zU6ann+X87wws8v\n8cSqCv77r9foH7ck1Z2egbHmtiOX2zl1rXNsP59V84upXV7urIKwvJzlFYWqZceZms5ioEQTu4Hh\nET6//wTffes6n373cv7sI+u0qrJ4LtjcFqzxNFxuH1vBem5xvpt4nOSzfnGpttGeITWdScJ09w/x\n3NeP8vrFu+z6pbU8975q/VKUlBCpue1CWw9HLrc7iedKOz88fROAObnZPLKsjC3uum+bKn2Uas+e\nhFCNBtVoYtHZO8Rv/t1h3r7exe5PbOATm5d6HZJITG519dNw2Wlqa7jSztvXuxi1YAysXVh6X61n\niSaTTkpNZzFQopmaQO8gn/qfhzl/s4e//VQNTz24wOuQRGbs3sAwx68Gxmo9jVc6uDfoTCZdXFbA\nZrePp7aqnLULS8nWyhZj1HQmcRXoHeQ3vnaYC2097P3NzXxg7XyvQxKJi6L8HJ5wtwYHGB4Z5ezN\nbhout9NwpYM3L93le285C4cW5+ewqdJH7XI/tVXlbFzmoyhft9FoVKNBNZpouvuH+OQLhzl3q5t9\nv7mZ9z+gJCOZw1rLtUCfO6y6nYbLHZy71Y21kJ1lWLeo1JlMWuVMKF1YVuB1yEmjprMYKNFMbGB4\nhM/8ryO8eamdFz5dq5qMCNDZN8SxlneGVR+/+s6K1Ut8c6itemfR0HRubkvppjNjzDZgi7V2V4Tj\n1cABoB2oBw5Ya5tDztkJNAN+AGvtvrD3mLRcpm5k1PKHLx7njaa7/NUzjyjJiLjK5uTy/gfmj9Xu\nh0ZGeft6lzusup1/bbrLd47f39wWXEJnY6WP4gxrbkvqtzXG1AE1wFacZBDOD+x2HwHg2bAksxs4\nYq09EHxujNkW+nyycpk6ay1/8b3T/ODkTb704Qf5d5s0ukxkIrnZWWM7lH72yRVYa2nt6BsbUt1w\nuYO/efkC1kJWcHSbW+upqSxnaXl6bxDnSdOZmxB81todYcfrgZcAf2iCCSnvsNaWhzyvA3ZZa7dO\npXwiajob7x//9TJ/9p3T1L+3mj/+0INehyMy63X1D3G8JTBW6zneEhgb3bagNH8s6WxeXs76xWWz\nYquElG46m4y1NoBTm7mPMaYmwuntQN1UymXqXr94h7/43tvUPTifL/zSWq/DEUkLpQW5vHfNvLFt\nM4ZHRjl3q9tdwcDp7/nBSWcyaX5OFhuWljl79FQ6e/XMLc73MvwZSblE49Zq2nGa0XzW2j1ukd89\nHirgvsYXrdxNYBLFpTv3+L1vNLJqXjF//WubtBumSILkZGexfnEZ6xeX8el3VwHvbJXQ2OI8/u61\nS+wdcRp3qioKqXGTTk1lOQ8sLJk1gwxSLdEcAtqDScEYs9cYU+926AeTSahgYvFPoVyJJop7A8M8\n+48NZBn42m/VZlyHpYjXFpYV8OENi/jwBmdb7P6hEU5d6xxbOPRnF+7wz8euAVCUl83GSp+TfCrL\n2VTpw1eYmovaptSdJEK/zEGcgQH7iJwogomlfQrl93FrTvUAlZWV0wk3rVhr+dK3T9F8u4evf/Yx\nlvkLvQ5JJOMV5GY7q09XObcyay1X2/vGajyNLR185adNjLgbxFXPKxpLOjWV5axZkBq1npRJNG7z\nVwdQHtLMFcAZ7gxOsvCFvcwHTr+OMWbS8vDPc2tJ+8AZDBCP7zCb7W9o5VvHrvFHW9fwuDtDWkRS\nizGGyopCKisK+fimJQD0Dg7z1tVOGls6ONbSwStn2zhwtBVwaj2PLPONJZ9NleX4PdjKI2USjWtP\nWFKoxh0Gba1tNMaEJww/TnNb1HKZ2NmbXfzpd07x5Kq5/P4HVnkdjojEoDAvh3evrODdK50Vq621\ntLT3uoknQGNLB3/76ju1nqqKQjYFE8+yctYuKkn4dgkpk2jcWsndsMPbgdBJnfvC5sVsBfbGUC5h\n+gZH+P1vNFI6J5e/emZjSlSzRWT6jDEsryhieUXR2Py33sFhTrZ2cuxqgMYrHbx28Q7fcvt63rN6\nLv/02ccSGlOyJ2zW4Aw33gb4jTFNwCFrbaN7yj53Zn8AWAnsDZ1saa3dZYzZGbKCQFMs5TLenh+d\npen2Pb7xuceYVzJ7h0+KyMQK83Lu26fHWsv1zn6OtXRQmJed8M/XWmdk7oTNN5ru8MkXDvPbj1fx\n5x9d73U4IjLLTHXCZupPPZWE6O4f4vP7T7BibhG7NClTRBIoZfpoJLn+8/87w43OPvY/9zhzklB1\nFpHMpRpNBnrtwh1ebLhK/XtXsnl5efQXiIjMgBJNhukfGuFPv3OK5RWF/Pu61V6HIyIZQE1nGeYr\nP23i0p17/NNnH6UgV01mIpJ4qtFkkKbbPXz1p0189JHFvGf1PK/DEZEMoUSTIay1fOlbp8jPzeJL\nH9H+MiKSPEo0GeIHJ2/yr8132flLa5lfUuB1OCKSQZRoMkD/0AjP/8sZ1i4s4ZOPaqVqEUkuJZoM\n8PdvXKa1o48vfXid1jITkaRToklzd3oG+B+vXOSptfN5crWW/xeR5FOiSXN/feg8fUMjfPFDGgAg\nIt5QokljF9u6+d+HW/jUu5azan6x1+GISIZSokljf3XwAnNys/mDp7QCgIh4R4kmTb19vYvvn7zB\n7zy5wpOtW0VEgpRo0tSXD56ntCCHz72n2utQRCTDKdGkoeNXAxw6c4tn31NN2Zxcr8MRkQynRJOG\nvnzwPOWFuXzmyRVehyIiokSTbo5eaedn52/zu+9fSXG+FucWEe/FnGiMMRsney7e+spPmigvzOVT\n71rudSgiIkCMicYY8wngZWNMVcjhTmPML8QzKJmesze7ePlsG595YgWFearNiEhqiLVGE7DWVlhr\nLwcPWGsvAb64RiXTsvfVZgrzsvn0u1WbEZHUEWuiKZvguH+mgcjMXG3v5btvXefXH63EV6h5MyKS\nOmJNNI9G6KP5BWBz/EKS6fjaz5vJMvC592ikmYiklpga8q21XzDGNBhjVgDNQLX751OJCE6m5m7P\nAC82XOXjG5ewqGyO1+GIiNwn5h5ja22tMaYGpxbTbK19Of5hSSy+/m8t9A+NsuN9WgVARFLPtIYm\nWWsbgcY4xyLTMDg8yjcOX+F9a+axan6J1+GIiIwznXk0nzfGtBtjRowxd40xn01EYDI1Pzx9k7bu\nAX778SqvQxERiSimGo0x5llgC7Adp29mJVBvjDHW2q8lID6J4h/euExVRSHvWzPP61BERCKKtels\npbX26ZDnl4BDxpjPxzEmmaKTrZ0cvdLBn31kHVlZxutwREQiirXp7MgEx9Vf44G/f+MyhXnZbKtd\n6nUoIiITijXRlE9wfGwipzHmc9MPR6bqbs8A3ztxnU/ULKW0QFsBiEjqirXp7BeNMbtw+meCqoFm\nY8wOwOAMe1Z/TYK91NDK4PAov/W4lpsRkdQWa6KpBr4Q5Zwd04xFpshay4tHWnh0hV9DmkUk5cWa\naHZFm6BpjGmerFxm7vCldi7f7eUPnlrtdSgiIlFN2kcTvvz/REkm9Dxr7bH4hCYTefHIVUoKcvjl\nhxZ5HYqISFTRajQ7jDFHcPpeJrMdeCU+IclkOvuG+MHJG2yvXcqcvGyvwxERiSpaotkObHP/Hkw2\nNuwc4x773TjGJRP47vFrDAyP8mtbKr0ORURkSqLWaNw/G7h/pFkoA+yLW0Qyqf975CrrFpXy0JKJ\ntgYSEUktkyYaa+0LAMaYTbh7zlhrxzWRGWP2JiQ6uc+pa52cvt7Ff/rYeq9DERGZsimNOgvt4DfG\nPIXTVNZurT3ulmurgCQ4cLSVvJwsPvbIEq9DERGZsunsR/MygDGmzB1t5gMOWWu74h2cvGN4ZJTv\nvXWdugfnU1aolQBEZPaIeZuAENuBPcAB4IvxCUcm8trFO9y9N8jHNqo2IyKzS6zbBGwEngPqgaPA\n3mA/Tozvsw3YYq3dFaFsJ87AAz+AtXZfPMtnq+8cv05pQQ7vf0DbAYjI7BK1RmOMqTLG/EdjzEXg\nZSCAs13AlliTjDGmzk0EO3Ca3MLLd+NsD33ATRAr3aQUl/LZqndwmB+dvsmHHl5Efo7mzojI7BJt\nZYAfA004m53tsNZWWGu/YK29FHbe81P5MGvtIWvtHibeVqDeWnsg5PlB7l87babls9KhM230Do6o\n2UxEZqVoTWe1OItoBoBqY8yKkLLgRM1yYCcz7KcxxtREONwO1MWjfDb7zrFrLCor4LEVfq9DERGJ\nWbREs89a+5fR3sQYszIOsfhxEkOogPv+vpmWW2sDcYgx6drvDfLq+dt89skV2kVTRGalaH00U52I\nuXumgeD02YT/ZA8mDn8cymel75+8wfCo5aMbF3sdiojItEyaaML7YmZ6XhSRahzBBNEeh/L7GGPq\njTENxpiG27dvxxpr0nz/xHVWziti3aJSr0MREZmWCRONMea/GmOeN8ZUJSmWdsaPRPMBuM1eMy2/\nj7V2n7W21lpbO29eag4ZvtszwJuX2vnQw4swRs1mIjI7TdhHY639gjGmDKg3xtQCR3D6bBKyAoC1\nttEYE54Q/MCheJTPRgffvsWohQ+uX+h1KCIi0xat6azTWvuX1tpngG8Cf2yMedEY8zljTCLacvaF\nzXvZyv39RDMtn1V+ePomy/xzWL9YzWYiMntNeWUAtx/mCzC2mvMeY4wFDlpr/3kq7+EOQa7D2ePG\nb4xpwlknrdH9jF3GmJ1usqgGmkLnxcy0fDbp7Bvi9Yt3+MwTK9RsJiKzWsyLasLYas7PgbOaszHm\nqzhzavZH2kYg5HWNOJM190xyzoRl8SifLV45e4uhEatmMxGZ9aaVaEK5qzkHV3R+1hjzHE5NQgtt\nzsAPT91kQWk+m5aNW6lHRGRWmXGiCeWufRbzIptyv97BYV49f5tnapdpkqaIzHoz2SZAEuSn527T\nPzTKBx9Ss5mIzH5xTTTGmF+N5/tlqh+dvom/KI9Hq2btggYiImNibjozxnwCWAFU4Izs8rl/VgBl\ngNaxn4HhkVF+eu42dQ8uICdbFU4Rmf2meyf7NeAu8BLOkOd91lo/7kg0mb7GlgCdfUM89eB8r0MR\nEYmLmBONtfab1tpaoNN5ao/h7FkTHAwgM/DK2TZysgxPrp7rdSgiInEx7VFn1toXjDErjDGfx5lD\nI3Hwk7NtbKnyU1qQ63UoIiJxMaNOAGvtJXe/mmPuHJrPxSmujNTa0cu5W91qNhORtBKX3mZr7ctu\ns9nT8Xi/TPWTs20AfGCtEo2IpI94D2vaFef3yygvn22jqqKQ6rlFXociIhI3U+qjMcb8AlCDM4T5\nIs6aZuO2C3AHBsg09A4O80bTXX7jsUotoikiaSVqojHGNODMlWnGmS9TjbMc/05r7X9LcHwZ442L\ndxkcHuWptQu8DkVEJK4mbTozxvwt8Ly1dpW19hfdP7OADwKPGWMuGGNKkhJpmnvlXBtFedk8ukKr\nAYhIeonWR2Ostd8MP2itPWStfRr4XeBrCYksg1hrefXcbR5fNZe8HK0GICLpJdpdrWmyQmvtIeAl\nDWuemSt3e7kW6OO9mqQpImkoWqKJOhHTrfHUxieczPTzi3cAeGKVEo2IpJ9oiaZiiu8zac1HJvf6\nhTss8c1hhYY1i0gaipZodhljjhhjnneHOEucjYxa3mi6w5Or5mpYs4ikpWiJ5gCwD1gJHDDGjEyQ\neLTW2TSdaA3Q1T/ME+qfEZE0FS3R7LXWvmCtfdrdBmA1ERIPsCPRgaar14P9Myun2kopIjK7TDph\n01r7ctjzZpyJmy8AGGNWAFuB5xMVYLr7+YU7rF9cSkVxvtehiIgkxLS3CQBn9WacVQLK4hRPRrk3\nMExjSwe/88QKr0MREUmYeM0O3Ben98kob15uZ2jEapMzEUlr8domoDMe75NpXrtwh7ycLLZUadkZ\nEUlfWu/EQ69fvMOWqnIKcrOzeH8VAAAKEElEQVS9DkVEJGGUaDwS6B3k7M1u3rVCo81EJL0p0Xjk\nzUvtADxWrUQjIulNicYjb15qJy8niw1LNWBPRNKbEo1HDl9qZ9Myn/pnRCTtKdF4oLt/iNPXO3lM\nm5yJSAZQovFAw5UORq36Z0QkMyjReODNS+3kZBk2Vfq8DkVEJOGUaDxwuPkuG5aWUZg3oxWARERm\nBSWaJOsbHOFEayePav6MiGQIJZoka2zpYHjU8li1BgKISGZQokmyw5fayTJQu7zc61BERJJCiSbJ\nDjffZf3iMkoKcr0ORUQkKZRokmhoZJS3WgNarVlEMooSTRKdudFF/9AoNcs1rFlEMocSTRI1XukA\noKZS/TMikjmUaJKosSXAwtICFvvmeB2KiEjSKNEkUWNLh5rNRCTjpNzUdGPMNqAaOAC0A/XAAWtt\nc8g5O4FmwA9grd0X9h6Tlnuhrbuf1o4+fuvdVV6HIiKSVKlYo/EDu4Em4BLQHJZkdrvHDrgJZKWb\nnKZU7pXGKwEA1WhEJOOkYqIBKAdWWmvLrbUHwsrqw44dBHbEUO6JYy0d5GYb1i/WRmcikllSrukM\nwFobAALhx40xNRFObwfqplLupWMtAdYvLtNGZyKScVIy0Rhj6nEShB/wWWv3uEV+93iogPsaX7Ry\nN4El3dDIKCeuBfjko8u9+HgREU+lYqI5BLQHk4IxZq8xpt7tbwkmk1DBxOKfQvlYonGTWT1AZWVl\nXL9AOE3UFJFMlnJ9NNba5rCax0Fgl/v3SDWSYGJpn0J56Ofss9bWWmtr582bN5OQo9JETRHJZCmV\naIwxPmOMdZvBggI4w53BSRbh1QIfjPXrRCv3RGNLgAWl+SwqK/AqBBERz6RUonHtCUsK1ThzYrDW\nNjK+1uLHaW6LWu6VxpYOairLMcZ4GYaIiCdSKtG4CeZu2OHtvNN0BrAvbF7MVmBvDOVJdbdngNaO\nPjYuU/+MiGSmVBwMsM+d2R8AVgJ7Q+fFWGt3GWN2hqwg0BRLebKdvNYJwMNLNX9GRDJTyiUat1az\nJ8o5MypPplNuonloiRKNiGSmlGo6S0cnWjupnltEqXbUFJEMpUSTYCevdao2IyIZTYkmgW53D3Cj\ns58N6p8RkQymRJNAwf6Zh1WjEZEMpkSTQCevdWIMrFeiEZEMpkSTQMGBAMX5KTe4T0QkaZRoEujk\ntQAblmqipohkNiWaBGnr6udW14BGnIlIxlOiSZDgigAacSYimU6JJkFOtHaSZWDdolKvQxER8ZQS\nTYKcutbJynnFFGkggIhkOCWaBDlxrVMLaYqIoESTEG1d/dzuHuChxUo0IiJKNAnw9o0uANYtVv+M\niIgSTQKcudENwIMLlWhERJRoEuDMjS4WlxVQVqitAURElGgS4MyNLh7UsGYREUCJJu76h0ZovnNP\niUZExKVEE2cX23oYGbVKNCIiLiWaOAuOOHtwUYnHkYiIpAYlmjg7c6OLObnZLK8o8joUEZGUoEQT\nZ2dudPHAwhKys4zXoYiIpAQlmjiy1nLmRrf6Z0REQijRxNGNzn46+4ZYp/4ZEZExSjRxdGZsIIBq\nNCIiQUo0cRRMNGuVaERExijRxNGZG91U+gsp1h40IiJjlGjiyFl6Rv0zIiKhlGjipG9whEt3tfSM\niEg4JZo4udDWjbWwdqFqNCIioZRo4uTCrR4AVi9QohERCaVEEyfn27rJy85iub/Q61BERFKKEk2c\nXLjVQ/W8InKydUlFRELprhgnF9q61WwmIhKBEk0c9A4Oc7W9j9Xzi70ORUQk5SjRxEFT2z0A1ixQ\nohERCadEEwfnb3UDsGq+ms5ERMIp0cTBhbYecrMNVRUacSYiEk6JJg4u3Oqmem6xRpyJiESgO2Mc\nnG/rZrX6Z0REIlKimaHewWFaO/pYo6HNIiIRKdHMUFPbPaxFQ5tFRCaQlhunGGN2As2AH8Bauy9R\nn3WhzRlxpsmaIiKRpV2NxhizG2i21h5wE8xKY8y2RH3e+VvOiLPlGnEmIhJR2iUaoN5aeyDk+UFg\nR6I+7GKbM+IsVyPOREQiSqu7ozGmJsLhdqAuUZ95/lYPqzTiTERkQmmVaHD6ZNrDjgUAjDG+eH9Y\n3+AIVzt6WaMVAUREJmSstV7HEDduX8wL1trykGM+oANYaa1tDjleD9S7Tx8Azs3go+cCd2bw+kyj\n6xUbXa/Y6HrFZibXa7m1dl60k9Jt1FkgwjG/++d9NR13oEBcRqMZYxqstbXxeK9MoOsVG12v2Oh6\nxSYZ1yvdms7agfAmMh+AtTZSEhIRkQRLq0RjrW1kfK3GDxzyIBwRESHNEo1rX9i8ma3A3kR/ZoLf\nP93oesVG1ys2ul6xSfj1SqvBAEEhKwNUA4FErgwgIiKTS8tEIyIiqSMdm84kRRljtrlLBIlICjHG\nJLR7Id2GNydVMhfvnM2MMXVADU5/WXOU08Xl/v8FsAU4Yq3d42U8qcydL1ePMxhoJYC1dpenQc0S\n7o+/hA5vVqKZJvc/zpHgumrGmN3GmG1h66wJYK09BBwyxlQwfvi5RGCM2Wut3RHy/KgxBiWbCX0x\nNLG416teP/4mZ4ypTsbnqOls+pK6eKdkDvfXefgw/b3AFz0IZ7bY5q72EdSMU4OWydXh3LsSSolm\nGrxYvFMyih/YGeHXpmqDE9saVnupBo54Fcxs4DZpv5SMz1KimZ6kLt4pmcVdk29z6Np8OL/ONfF4\nAmHrGNa4x9TMODlfslZMUR/N9Ph4Zw21oGDi8RN5zTWRKXNXuQDGfrzUAZu9iyj1udfpaWA78KzH\n4aS0ZPcnq0YzPVNevFMkDvYDT4XVcCSMtTZgrd1nrd0KvBDWZyMut0k2qT+GVaOZHi3eKUnhjm7c\nHVrDkfGMMeHNQHvdh0adjVcDVIf0NW8BfO5w+gOJ+EGjRDMN1tpGY4wW75SEctfsO+gOD8cYU6OE\nM57bqX3QGFMe/kMvQgLKeOFNZm7NrzqRfVpqOps+LxbvnJWMMTXur6VtwNPGmJ0TjNwTl3vz9AMN\nxhif29zxjMdhpaoGYF9YQtmK8+tcSWYSbpLZjlPD2ZmowUxa62wGtHinJELIrrDhDlhrtyc7ntnA\n/eESnF5QAVoZIJUo0YiISEKp6UxERBJKiUZERBJKiUZERBJKiUZERBJKiUZERBJKiUZERBJKiUZE\nRBJKiUZERBJKiUZERBJKiUZERBJKiUYkRbmLkXYYY/aHPNcOrjLrKNGIpK46YAWwN7iJl1YjltlI\ni2qKpDi3FuPXDpsyWynRiKSwYFOZajIym6npTCRFuZudjSWZ4HOR2UZbOYukIHcjr1qg3RhzCKe/\nJoCz0Z7IrKIajUiKCWku24eze+slYIu19pCngYlMk/poREQkoVSjERGRhFKiERGRhFKiERGRhFKi\nERGRhFKiERGRhFKiERGRhFKiERGRhFKiERGRhFKiERGRhPr/iav2m0iBVAEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11048fdd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ztab = np.linspace(0,4,100)\n",
    "DA = results.angular_diameter_distance(ztab)\n",
    "plt.plot(ztab, DA)\n",
    "plt.xlabel('$z$')\n",
    "plt.ylabel(r'$D_A /\\rm{Mpc}$')\n",
    "plt.title('Angular diameter distance')\n",
    "plt.ylim([0,2000]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Get matter power spectrum using linear model or Halofit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get matter power spectrum at some redshift\n",
    "pars.set_matter_power(redshifts=[0.5], kmax=2.0)\n",
    "\n",
    "#Linear spectra\n",
    "pars.NonLinear = model.NonLinear_none\n",
    "results.calc_power_spectra(pars)\n",
    "kh, z, pk = \\\n",
    "         results.get_matter_power_spectrum(minkh=1e-4, maxkh=2, npoints = 200)\n",
    "\n",
    "#Non-Linear spectra (Halofit)\n",
    "pars.NonLinear = model.NonLinear_both\n",
    "results.calc_power_spectra(pars)\n",
    "kh_nonlin, z_nonlin, pk_nonlin = \\\n",
    "        results.get_matter_power_spectrum(minkh=1e-4, maxkh=2, npoints = 200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Construct interpolating function for the linear case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Pk = interp1d(kh, pk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEUCAYAAAAWSSD2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xmc3fP1x/HXIaHWjAm1VImZRK1B\nMinV30+bzVbRLUtbraWYQVFVEtRaiiS0FI1MUEr7I4mtuiAJparIJJSiTSSqpVWJZKJILcn5/XG+\nI2Myy507997vXd7Px+M+Zu693/u9J/nO3DOf7XzM3REREcm1ddIOQEREypMSjIiI5IUSjIiI5IUS\njIiI5IUSjIiI5IUSjIiI5IUSjIiI5IUSjIiI5IUSjIiI5EWvtAPIt80339z79euXdhgiIiVl3rx5\nS919i56co2wTjJmNAkb179+fpqamtMMRESkpZvZST89Rtl1k7n6Pu9f36dMn7VBERCpS2SYYERFJ\nlxKMiIjkhRKMiIjkhRKMiIjkhRKMiIjkhRKMiIjkRdmug5HKtno1vPoq/OMf8NprsGRJfP3Pf+Ct\nt+DNN+MYgHXWgQ02gI03hk03hW22gW23XXPbaKN0/y0ipUoJRkraypXw/PPw5z/H7dlnYdEi+Nvf\n4J131j5+3XUjkWy0UXzvHolm5cpIPO++u/Zrqqth4EAYNAgGD46vO+4YiUlEOqYEIyXDPZLJ44+v\nuT3zDKxaFc+vvz7stBPsthuMGgU77ADbbQdbbgkf/ShssUW0VMw6fo+VK+Gf/4yWz8svx+3FF+Gp\np+Caa9YkraoqGDkSDjoIDjggWj0i8mFKMFK03GHBAnjwwTW3JUviuT59YMgQOOMM2GuvSCq1tdCr\nhz/RG2wQ56mtXfu5996LBDdvHjzyCNx7L8yYEc/tsQd87Wtw2GHwsY/1LAaRcmHunnYMeVVXV+eq\nRVY63n47EsmvfgW//nW0JCA+tIcOjdu++xZHF5V7tKDuvRfuugv++MdoHQ0bBkceCWPGRKtKpBSZ\n2Tx3r+vROZRgJG1LlsCdd8Ivfwlz5sB//xvjJPvvH91PQ4dC//6dd20Vg4UL4ZZb4Oabo1ttyy3h\nW9+C446L7jmRUlLRCcbMxgON7t7c2XFKMMVp6dJIKtOnR4tl1SqoqYFDDonbfvuV7l//7jBrFvzo\nR9G6WX99+OY34ayzYlaaSCmo2ARjZlXAPGCwEkzpePNNuP12+MUvoqWyahUMGABjx8Zt992Lv5XS\nXc8/H4nmxhujS6+hIcaNtt467chEOpeLBFOqEy3rgMVpByFdc4ff/z7+gt966xibWLQIJkyAJ5+E\nv/4VLroopgGXW3IB2HlnaGyMyQpf/3rMROvfP/7NK1emHZ1IfqWWYMxstJlN7OC58cnz9WZW3+a5\nGpRcit4//gE/+EEMxu+3X8y2GjcuZl8tXBjP7blneSaV9vTrB9ddB3/5S0xtPucc2GWXaNGVYCeC\nSEYKnmDMbEQyftIAVLXz/ERgsbvPdPdGoNbMRrc6pMbdlWCK0OrVMfbwhS/EB+rZZ8eYw003xar6\n666DT3+6cpJKe/r3h5kz4YEHYJNNYPTomHX25z+nHZlI7hU8wbj7bHefBMzv4JB6d5/Z6v4sIhlh\nZoMADagUmRUr4Morozto//3h0UejC2zRohjAP/xwlVtpa+hQmD8fpkyJqc6DBsH557dfSUCkVBXV\nGEySQNpaBoxIvq8BRiQtmhpgbDLgLyl45pmYgvuxj8Epp0RJlZtvju6xiy+OWWHSsV694v/vL3+J\nSQ4XXBCJ5vHH045MJDeKKsEA1URCaa0ZYuZY0m02s+UxoKqrWWSSW+5w330wfHgMzN90U4ytzJsX\nCw2//vXSnV6cls03j/Uzv/51tAY/9Sk49dSojSZSyootwVQRSaa1loTzweNJN1tt0tUmBfDee/Eh\nuOeecOCBMftr4sSo1XX99fGXt/TMwQdHsc7jj4+pzbvvHmM1IqWq2BJMe62RlsTStmXToWT2WZOZ\nNS1pKV4lWXnzTbjiiqjN9Y1vxNqVG2+ExYth/Hjo2zftCMvLppvGVOaHHooutOHDoxvtjTfSjkyk\n+4otwSxj7ZllVQDd6Qpz90Z3r3P3ui1UoyMr//53zALbbjv4zndiPOVXv4Knn4YjjoD11ks7wvK2\n335Rwfm002DatCjmee+9aUcl0j1FlWDcfT5rt2KqgdkphFORXnkFTjoJtt8+BuqHDYPHHoPf/Q4+\n97n0C0xWkg03hMmTY1bexhvH+pmjjoLly9OOTCQzxfhx0dhm3ctIYGpawVSKf/4TTj45usKuvTa6\nw/7611izsffeaUdX2fbeO6oenHVWzNLbddcoDCpS7NJYaDkoWWg5mphmPL719GR3nwDUJCv5xwOL\n2qyLyfR9RplZ44oVK3IXfBn617/g29+OLrApUyKxLFgQ3TIDBqQdnbRYf/2ofvDEE1GZ+fOfj/1n\nli5NOzKRjpVkscvuULHL9v3rXzELbOrUmCF25JHxF7LWrhS/d9+FSy+NemZVVTEpYPToyq6QILlX\nycUuJUuvvrpm0P7qq+Ov4AULooyLkktpWG89OPfcWHu0/faxSHP06Li2IsVECaZCLFkC3/1u7FN/\n1VXw1a/GGMv11yuxlKrdd4/FrZdeGos0d9wxJgW8807akYmEsk0wGoMJb78ds8Fqa2M9y7hxUZrk\nhhva33deSkuvXlH37emnY2rz+PExCWDmzCg+KpKmsk0w7n6Pu9f36dMn7VBSsWpVtE4GDIDvfW9N\nxd4bb4yKvlJedtwx1im17KA5ZkxUV7jrLm0HIOkp2wRTqdzjg2aPPeCYY2Kh5O9/Hx80O++cdnSS\nbwccEK2Zm2+O1usXvwiDB8M99yjRSOEpwZSRJ56Az34WRo2KmUYzZ8Yivf/5n7Qjk0Jad90oOvrc\nc9FiXbECDj00Es0dd6jrTAqnbBNMJY3BvPBCzCTae+8YX/nJT6Jo4pe/rKmrlaxXryjr0zLm9p//\nxM/EHnvArbdGN6pIPpVtgqmEMZjXXouyLjvvDL/5DZx3XiSb44+H3r3Tjk6KRe/eUWLm+eejIvaq\nVTGLcNdd4Wc/g/ffTztCKVdlm2DK2VtvxSK7/v1j9f0xx0RiOf/82IZXpD29esFhh8Vkjxkz4CMf\niRbOJz4R66C0m6bkmhJMCXGPwdsBA+Ccc2DEiOgKmzIFttoq7eikVKyzTizMfPLJqGnWty8ce2z8\nXDU2qutMckcJpkTMmwef/nTsb7/ttvDIIzFg+4lPpB2ZlCqzmBDy+OMxvXnbbaGhIX7Onn467eik\nHJRtgimXQf4lS6C+HoYMgUWLYrD2scfiQ0AkF8xievMjj8QYzaJFUFcHP/yhZpxJz5Rtgin1Qf73\n34+SLjvuCD/9adQPW7AgBmu1J4vkg1mM0fzlL7H3z3e/G9ObS/xvNEmRPqqK0IMPwl57xf4sdXXR\nXXH55VCiuVJKTN++0f161VVw332w776xRbZIdynBFJG//z3WswwbBm++Gb/k99+vFfhSeGZw4onx\n8/evf8E++8TsM5HuUIIpAitXwoUXwk47RZmX738/VmF/8YtaKCnpGjo0Kjb37h3fP/NM2hFJKVGC\nSZE73Hkn7LJL7O9xyCGxGO6cc2CDDdKOTiR84hPwu99FEU0lGekOJZiUPP98zNz50pdg443hgQdg\n+vTYQEqk2AwYAA89FIszDzwQ/vGPtCOSUlC2CaZYpymvWAGnngoDB8LcufDjH8eCt6FD045MpHO1\ntfDb38b44EEHQXNz2hFJsSvbBFNs05RXr47pxjvuGBt/ffObMe34pJOihIdIKdh995h8smBBVANQ\nHTPpTNkmmGIydy586lORVGpr4/7UqbDFFmlHJtJ9w4fHz++cOXDGGWlHI8VMCSaP3ngjWih77x1T\nkH/2M/jDH2JfDpFSdtRRcMIJsT7r1lvTjkaKVaedM2b2JaAmy3Mvdvc7snxtybvzzkgu//wnfOtb\n8IMfwKabph2VSO786Efwpz/B0UfHTMiBA9OOSIpNV73/+wMTszz3eKDiEszLL8cCtbvvjl+422+P\nFoxIuVlvvSj7P3hwrNmaNw+qqtKOSopJVwlmnru/mM2JzWx+Nq8rVatWwTXXwPe+F99PnBj1w7Tx\nl5SzrbeOrbk/85noNrvjDi0OljU6HYNx92nZnrgnr82FQk5TfuqpGMT/9rejyvGzz8L48UouUhn2\n3RcmTYK77opuM5EW7SYYM9vLzG4zs6OT+5ea2WozW2hmexY2xOwUYpryW2/B6adHQcqXXoJf/CLW\nCeywQ97eUqQonXJKdJNNmACPPpp2NFIsOmrBjAOmAwPMbAowF6hNHm8oUGxF7be/hd12g8suW7Pf\n+Ve/qu4BqUxmsVfRdtvBuHGxj5FIRwlmlrvf7u5nALOT71909/nA7ALGV3RefRW+8hU4+OAom/HQ\nQzBtGlRXpx2ZSLqqqmLQf+nS2CnzrbfSjkjS1lGCqTKz+8xsU3e/HcDMhpvZacBmhQuveKxeHfuV\n77xzTEG+4IIYe9lvv7QjEykegwbB//1fLCb+yle00r/StZtgkqRynLu/0erhZcAKd7+uIJEVkeee\ni1kyDQ2wxx6xAdi550Z1WRH5sC98Aa6+OraeOOGEqBoulanDacptpye7+5PAkx0db2b93P1vuQst\nff/9byyQnDgxKh7fcAMceaTGWUS6cvzxsSbs4oth223jDzKpPFmXWTSztuvSJwDH9yyc4vHgg9Fi\nWbgw9in/4Q/hox9NOyqR0nHRRZFkzjsPNtwQTjst7Yik0Lpdi8zMjjWz1cByoLnV1/ocx5aK11+P\nWWHDhsWCyfvvh1tuUXIR6S4zuO662Ab89NPj9t57aUclhZRNsctaYDN3X9fd12n5CkzOcWwF5Q43\n3xzbFt9yS1SJfeYZGDky7chESlfv3vDzn0eX2WWXxVjmSy+lHZUUSjYJZpa7t7c8/pKeBpNL3VnJ\n/8ILkUgOPxz694f58+GSS6JZLyI906sX/OQnUXX5z3+GPfeM8jJS/rJJMG5m/dp5/NiehZJbmazk\nf/fdGITcffeYVnnNNVFOf/fdCxioSIUYNy52b+3fH8aM0YLMStDlIL+ZHQuMbv0QMNjMlhFjLy2P\n7QVclvMI8+TRR6G+PuqGffnLsXXxNtukHZVIeautjd+9yZPh/PPhgQeidTNmTNqRST5kMouslli9\n31l1ZCPK8xe95mY480y49lr4+Mfhl7+MVcciUhi9e8NZZ8Ghh8aEmrFj44+8a66BLbdMOzrJpUwS\nzG3JGphOmdnrOYgnb9yj3/fkk+G116I434UXxvoWESm83XaDP/4xBv/POw9+9zuYPj1mcEp56HIM\nJpPk0p3j0vDSS9FKGTs29q94/PEoK67kIpKuXr1ixuaTT8JWW8EBB2gCQDnpasvkSwEHppbqKv1/\n/xt23TVaMJdfHi2YXlkvLxWRfNhll5hg87nPRVXyDTaI76W0mXdRKMjM+hCLKOuIsv2NbWqUFTWz\nOj/44CauuQb69Us7GhHpzIoV0UX27LOxyFnFZNNjZvPcva5H5+gqwbR5wx2I/WB2AGYB04s92fTv\nX+cLFzapfphIiVi6FP73f2HZsqhYvvXWaUdUmXKRYLq1DibZE+YMdx8HzAMmmdkUM/tST4LIp6oq\nFacUKSWbbw633w5vvglf+5pK/peybBZaAjGo7+7HufvxwAozuzZJNpoDIiI9sssuMGVKzCy76KK0\no5FsZZ1gWnP3Oa2STa2ZTTezoiodIyKl5fDD4etfjwTzZNHOUZXOdGsMphTV1dV5U1NT2mGISBaW\nLYtZoFtuGeWcevdOO6LKUfAxmFLSnWKXIlKcqquj6saf/gSXXpp2NNJd7SYYM9vLzG4zs2OS+5ea\n2WozW2hmexY2xOxkUuxSRIrf5z8fhTEvugj++te0o5Hu6KgFMw6YDvQ3synE+pfa5PGGAsUmIgLA\nlVfCRhtFgdrVq9OORjLVUYKZ5e63u/sZwOzk+xfdfT5R+FJEpGC23DJqlj38MNxwQ9rRSKY6SjBV\nZnafmW3q7rcDmNlwMzsN2Kxw4YmIhKOOgs9+NrZefvXVtKORTLSbYJKkclybVfrLgBXufl1BIhMR\nacUMpk6FlSvhsMO0ALMUdDiLzN1fbHP/SXefBtB2oL9UBv5FpLTtuGMkmQcegNNOSzsa6Uq3pymb\n2ZeBOW22TV6hFfwiUghHHBH7OV15Jfz0p2lHI53JZh1Ms7v3bV2+P2ntVOUsKhGRTkyeDMOHw3HH\nwWOPpR2NdCSbBNPRwpLqngQiIpKpXr3gtttg223hoIPgkUfSjkjak02C+WQ7YzDDgMG5CUlEpGt9\n+8ZYzJZbwsiR8Mtfph2RtNXtBJOsjbnOzF43s7lm9jowEZiQ8+hERDqx/fbRehk4EL74xRiXKfPy\niiUlq1pkSQG0kUAjMNbdhxT7xmMiUp423xzmzIFRo2Lw/+ij4Z130o5KoGf7wcx392nuPieXAYmI\ndNfGG8Mdd8C558bMsuHDobk57aikbKspi0hlWWcduOCCGPx/4olY9b9kSdpRVbasEoyZnW5my8xs\nVTIWc3SuA+vkvUeb2Qgzm2hmgwr1viJSGsaOhV/9Kiovjx4N772XdkSVK5uFlscCQ4AxQH+iwvIB\nLaX988nMaoAGd59NVHg+M9/vKSKlZ//9oyjmww/DqaemHU3l6pXFa2rdfWyr+y8Cs83s9BzF1CF3\nX2xmY5K7Q4BZ+X5PESlNX/0qzJsHl18OQ4bEFsxSWNl0kc3t4PH53TlJ0tU1sYPnxifP15tZfevn\n3L3ZzEYDVe7e2J33FJHKcumlsN9+cNJJ8MoraUdTebJJMB2V6/9ghX9n3WXJ+Ml4YuOytcrLJEln\nsbvPTBJIbZJQPuDuM4FZHSUoERGIFf/XXx/jMMcfrzUyhZZNgtk/2Tr5vla3hUBD8v39xMLLdrn7\nbHefRMctnvokgbSYRbKLppkNMrMRyePzgfFZxC8iFaR/f7jwQrjnHpg+Pe1oKks2YzA1wBldHJPV\ntsodzApbBrQklTqi1TM7+bo4m/cRkcpyyilw663w7W/HBIDNtG1iQWSTYCZ0tbjSzLL94K8mEkpr\nzck5q9y9MRmXGU1UEhjT9gQiIm2tuy40NkJdHZx5Jlx7bdoRVYYuu8ja7vPSUXJpfZy7P5llPFWs\nXZW5JeFUJ+duTMZnGty9WxMLRKRy7bVXtGSmToU//CHtaCpDJmMwDWa2iZlt2tmN3LQm2ivu0JJw\n2rZsOpS0cprMrGmJlvKKSOKCC2C77aC+Ht59N+1oyl8mCWYM8cG/PPna8n3rWzNQ39EJumEZa88s\nq4KYnpzpSZJWTp27122xxRY5CEtEysHGG8M118Bzz8HZZ6cdTfnLZAymZcC+iY4H1Y2orNwj7j7f\nzNomkmpiUF9EpMcOOSR2wpw8GT75ySgnI/nRZQsmqZg8Lbk7GBjs7iva3JqBqTmKqbHNupeR2Zzb\nzEaZWeOKFStyFJaIlIsrroB99oGjjoJnn007mvJlnsXKIzMbDjiwzN2f6uZrBxHTjhuI1sklwOzW\nA/bJQszFxJTo5p6s2K+rq/OmpqZsXy4iZeqVV2DwYFh/ffj972NsRtYws3nJ3l/ZnyObBNMqgD5E\nq6aKSBJFt+mYEoyIdGTePBg2LDYte+gh2HbbtCMqHrlIMD3dD2YMMAmYiSobi0iJGTwY7r8fli6N\n/WNUryy3sinXv6eZXWtmq4lurqnuvo67F1WC0RiMiGRi773hvvvgtddgxAhYlvGCCOlKRgnGzPqZ\n2Wlm9gIwh5iWXOvuQ1pNACgq7n6Pu9f36dOn64NFpKLts09sUrZ4MRx6KKxcmXZE5SGTlfz3A4uI\n/Vca3L2vu5/h7i+2Oe6SPMUoIpJ3++0Ht9wCjz4KRx6pysu5kMk6mDqiuGUzUGNmO7R6zojZZJsR\nlY2LqptMRKQ7xoyBRYuiXtmBB8Y0ZsleJgmm0d0nd3WQmdXmIB4RkVSdfjrcey+cfDJ85jNQU5N2\nRKUrkzGYTBc5FtXmXxrkF5FsrLsu3HQTrLNObLO8alXaEZWuThOMmV3SdqylI+2MyUzpSWA9pUF+\nEcnW9tvDVVdF1eUrrkg7mtLVVRdZrZmdlsV5DeibxetERIrCN74Bt98O3/seHHww7Lxz2hGVnq4S\nzIQenHtm14eIiBQns9g7ZtddY1bZH/4AvbLZorGCdfrflWn3mIhIOdpqK/jJT+ArX4nqy2dqnmy3\n9LRUTNHSIL+I5MLYsVHS/7zz4Jln0o6mtJRtgtEgv4jkglm0YjbbLJLNf/6TdkSlo2wTjIhIrmyx\nBdx6KyxYEIsvtco/M0owIiIZGDoUJk2KmWWTu1x6LpDZSv4PMbN+wCBis7AqYmOwZnd/IKeRiYgU\nmVNPhSeeiMH+gQOjnIx0LOMWjJmdnhS+nAjUEmtdViTf729m95vZFDPbMz+hioikywxuuAF22y1m\nli1YkHZExa3LHS2T4pYNwK1dbY+c7HBZDyx39+tyFmUWzGwUMKp///7HLly4MM1QRKTM/O1vUFcX\nYzNz58LGG6cdUe7lfUfLJLkMT8rzd5pcANx9RVIYc46ZHdOTwHpKs8hEJF/69YMZM6IFc/LJaUdT\nvDpNMO7+YjYtkWxfJyJSKoYOhbPOgp/+FG67Le1oilM2WyZf27qQpZntYGbDchuWiEjxO/fc2A2z\noQFeeintaIpPNtOUZwGNLXeS1soDZval3IUlIlL8eveGn/8cVq9Waf/2ZJNgmoktlDGzTc3skqTi\n8pCcRiYiUgJqauDHP4aHH4bLL087muLSZYIxs9NaTz129zlAXTJjbBIxXdlo1aoREakkRxwBX/oS\nnH02PNXldKjKkck05WXAMmAHYDbRRTYbqAFq3P2yfAeZDU1TFpFCWroUdt8d+vaFpib4yEfSjqhn\n8j5NOTHB3fsTG4g1Av2JvV5mAiPN7JhkdX9R0TRlESmkzTePGWXPPquy/i26TDDuPi352uzut7v7\ncUnC2YxIMvsD881sbn5DFREpbgceCCecENssz56ddjTp67KLLOMTmfVx96LbfKWurs6bmprSDkNE\nKsTbb8OgQfDmm/D001BdnXZE2SlUF1lGijG5iIgU2oYbwi23wGuvwTe/Wdml/bssFZNNyZdsXyci\nUg7q6mDiRLj7brjqqrSjSU+XpWKIumLXZrJaP1kXczpRv0ylYkSkYp1yChxyCJx+Osybl3Y06ch4\nDMbMjgXGAA7MB15PnupL7AtTSyzAnJQkpqKgMRgRScvrr8Oee8JGG8GTT8IGG6QdUeZyMQaT8YZj\nyWyyackCyzoiqVQDLwKL3f3JngQiIlJu+vaFG2+EESNi6vIVV6QdUWF1e0fLZDB/Th5iEREpO8OH\nw4knwpVXwqGHwrAKKg2cs1lkxcbMRplZ44oVmtwmIumaOBEGDICjjoJK+kjKWYJJBvg3LZbZY1rJ\nLyLFYsMN4Wc/g5dfhu98J+1oCqfHCSaZYXYfMINY2T+hx1GJiJSZffaBM86IcjJ33ZV2NIXR7TGY\ndsxIKiwDYGZ75eCcIiJl57zz4N574ZhjYO+9Yeut044ov3LRReZmtmmr+zvk4JwiImVnvfVig7K3\n366MVf65SDCNxGLMuWbWBEzLwTlFRMrSTjvB5MnRkrn22rSjya9cJJgGdx+S3OqAETk4p4hI2Trh\nBDjgAPjud2HBgrSjyZ8eJ5jW4y+JzXp6ThGRcmYGN9wQm5J94xvw/vtpR5QfWSUYM1uYTEkebmYv\nmNl9ye1+YjaZiIh0YpttYOpUeOIJuPjitKPJj2xnkdW5+xtmthgY2br2mGaRiYhkZswYOOww+P73\n4eCDowpzOcmqBdOy94u7v9hOYUvNIhMRydDVV8NWW8Hhh8N//5t2NLmVVQvGzF4gKid/8BBRZdmA\nwcAdPQ9NRKT8VVXBddfBQQfBuefCpElpR5Q72XaRNbQzuA+oi0xEpLsOPBCOPRYuuwy+8AXYd9+0\nI8qNjLrIzOySZBD/aGh35tgHVLZfRKT7Lr8ctt8ejjgC3nor7Whyo8sEY2bXEt1eK4DJZlYS8x1U\nTVlESskmm0SdshdeiJpl5SCTFky1u+/v7mPdvRqoNrN++Q2r51RNWURKzWc/G1stX301zJqVdjQ9\nl0mCWdz6jrsfB4zOTzgiIpXt4oujnMxRR0Fzc9rR9EwmCeb1dh6zXAciIiKwwQaxd8yrr0ZrppRl\n1EXWzmMlnldFRIrXkCFw5plw001wzz1pR5O9TBLMhGQG2WlmtkfyWJkXmRYRSdc558DAgVBfD8uW\npR1NdjJJMDOB2cAngQfNbBWRdC4xs2Ete8EUy1bJIiLlYL314MYbYelSOOmktKPJTiYJZqq7T241\ni2wAMAmoJZLPcjNbiLZKFhHJqb32grPPhl/8Au68M+1ous+8h1uqmdkOwEjgEnfvm5Oocqiurs6b\nmprSDkNEJCvvvRfbK7/yCjz7LGy+eWHe18zmJXt8ZS0X+8G86O6NaCdLEZGc6907usqWL4cTT0w7\nmu7JxY6WLS7J4blERCQxcCCcdx7cdhtMn552NJnLWYJpKeEvIiK5N2FCTF8+4YRYI1MKctmCERGR\nPOnVKxZgvvUWHHcc9HD4vCCUYERESsROO8Xul3ffDTNKYHP6kkswZlaf3KaaWU3a8YiIFNJ3vhNb\nK594IrzeXiGvIlJSCcbMBgFNyay1GclNRKRi9OoF118fs8qKvVZZSSUYoAZoSL5vSu6LiFSUgQOj\nVtktt8BvfpN2NB3r8ULLrN/YbDQwxN3XqgBgZuOJbQKqAZIWS8tzVe7enLx+nLuP6ex9tNBSRMrR\nO+/AoEHwxhuxAHPTTXN7/qJYaNldZjYiSSANQFU7z08EFrv7zCSx1CbJBAB3b6nkPA44thAxi4gU\nm/XXj66yV14p3h0wC55g3H22u08C5ndwSL27z2x1fxZrusWAD1o4x7ZKNiIiFWeffWIcZsoUePjh\ntKNZW1GNwSSD+G0tA0a0OmY00Jh0k41o53gRkYpx4YVQUwPHHAMrV6YdzYcVVYIhxlza7nzQDDH2\nkiSgacA8M1uOKjiLSIXbaCO4RrnVAAAKFUlEQVSYNg0WLoTzz087mg8rtgRTxdo7aLYknGp3n+/u\nm7l7bfJ1ZIHjExEpOsOGRQvmssugmOY0FVuCaW9MpSXhZLynW7IQs8nMmpYsWZKbyEREitjkybDV\nVnD00fDuu2lHE4otwSxj7ZllVfCh2WNdcvdGd69z97otttgil/GJiBSlqqoY7H/6aZg0Ke1oQlEl\nGHefz9qtmGpiy2YREenEoYfCuHEx8P/cc2lHU2QJJtHYet0LsVvm1O6exMxGmVnjihXaRUBEKseP\nfwybbBJdZatWpRtLGgstByXrWEYDY81sfOvpycnK/hozG50ct6jNupiMuPs97l7fp0+f3AUvIlLk\nPvpRuPJKeOwxuPrqdGNJrVRMoahUjIhUGncYNQoefBCeeSbWyXRXSZaKERGR/DKLAf9114X6+vQ2\nJyvbBKMxGBGpZB//eMwmmzMHbrghnRjURSYiUqZWr45FmE89FbPKttkm89eqi0xERDq0zjpRRuad\nd+CEEwrfVaYEIyJSxgYMiHUxd98NMwq8B7ASjIhImTvlFKirgxNPhNdfL9z7lm2C0SC/iEjo1Ss2\nJ1u+PJJNoZRtgtFCSxGRNQYOhLPOgltugd/8pjDvWbYJRkREPuyss2CXXeC44+CNN/L/fkowIiIV\nYv31Y03Myy/DGWfk//2UYEREKsjee8c4zJQp8NBD+X2vsk0wGuQXEWnfhRdGfbJjjoG3387f+5Rt\ngtEgv4hI+zbaCK67Dl54Ac49N3/vU7YJRkREOjZ0KDQ0wI9+BPPn5+c9lGBERCrUpZfG/jH19fD+\n+7k/vxKMiEiFqqqKzcnmzcvP5mRKMCIiFWzMGDj4YDj7bPj733N77rJNMJpFJiLSNTO45pqotHzi\nibmtuFy2CUazyEREMtOvH1xwAdxzD9xxR+7OW7YJRkREMnfKKbDnnnDSSZCrjh8lGBERoVcvaGyE\nf/87apblghKMiIgAMGRIjMNMmZKb8ynBiIjIBy66CD72sdycSwlGREQ+sMkmasF0SdOURUSyc8gh\nuTlP2SYYTVMWEUlX2SYYERFJlxKMiIjkhRKMiIjkhRKMiIjkhRKMiIjkhRKMiIjkhRKMiIjkhXku\ni/8XITP7D/DXtOMA+gC5XPWZ7fkyfV0mx3V1TEfPd+fxzYGlXcRRCLm8fvm+dpke29kx3X2uo+OL\n4fqV2u9eJscW4nfvE+6+SWdBdsndy/oGNKUdQxJHYzGcL9PXZXJcV8d09Hx3Hi/H65fva5eL69fd\n5zq5pqlfv1L73cvk2FL53VMXWeHcUyTny/R1mRzX1TEdPd/dx4tBLmPL97XL9NjOjunuc5Vy7Xpy\nvlxev5L43auELrImd6/L4LjRQDMwErjN3efnPTjpUjeuX33y7WBgorsvzm9kkolMr19ybA2wzN2b\n8xyWZKCja9fqs7IGmN3Z71oltGAauzog+cFucPfZwFzgzLxHJZnK5PoNIprzjcCM5CbFocvrB2Bm\nI4CpQEbJSApirWuXfFYOcffZye/bxM5OUPYJJvlP6OqYxcCY5O4QYFZeg5KMZXL9iL+kGpLvm5L7\nUgQyvH4kf9yp1VlEOrh2o4FFre4P6uwcvXIaUZ4lTbMh7j6hnefGEz+g1ZD5D3YLd29Ozl/V3ddK\nZvJ1/dx9ppnNTu6OAGZ3drxkJ5+/f1J4WV7PvrT5Q8DMqjrq1iyJFoyZjUj+wQ1AVTvPTwQWu/vM\n5D+iNvnP6xZ3nwnMSs4nOVKI69fqB3wccGxPY5Y1CvX7J4WRg+tZnfF7ldIgf/IPr3L3hjaPL3f3\nzVrdHwFMcPeRyf162vmPBOa7++ykD786+b4GWOTulr9/SWXK1/Vr9brxxHRLDRLnQQGu31RgRuvH\nJH+yuZ7J71hzS4um7bFtlVQXWXuS5NDWMqKrBMiouV5H/ALMTr6qL7hAcnT9Wpr7jUlX5wh9SBVG\nrq6fFIcMrudMYEJybMtnZodKPsEQzbVlbR5rhs77Bltz90Yzq08+pEayZsBf8q/H1y/5pZgGLDOz\namKgXwmmMHp8/ZJjR5PMIDOzZVomkJqurudiM5uXtGoGkSSbjpRDgqli7T7Blv+gapL/nK60+itr\nZo7iksz0+PolH0YdNtMlr3L1+zcT/e4Vgy6vZ6vPyi7/iCuJQf4utPcD3PIf1DYTS/HR9Sttun7l\nJafXsxwSzDLWHkCsgg/NLJLipetX2nT9yktOr2fJJ5ike6TtP7wa9cGXBF2/0qbrV15yfT1LPsEk\nGtvM0x5JlJ2Q0qDrV9p0/cpLzq5nSayDSWYJjSAWBlUDlxBF1ua3OqZl5WkNHx6IkpTp+pU2Xb/y\nUsjrWRIJRkRESk+5dJGJiEiRUYIREZG8UIIREZG8UIIREZG8UIIREZG8UIIREZG8UIIREZG8UIIR\nyYCZ1ZjZjJ7udpr2bqnJthRTteOkFIISjFS05AN3ebJyuStz3X2CmQ0ys4nJTo3dea+qVt/XJOdY\nbmazuojPk+Q2oqPjMpWsyJ5IrNAWySslGKloyQfuMrpRzC8pqfE63d+Yrp6kppO7L3b3CUSZjhHJ\nVt2dOVa7dEqpUYKRipa0Kmqy2EExmx0Xh7h72+24m4mNttbqskpqRjVl8T4iRUEJRirdCApQWj5p\noczt4OmpROHBtqqJgoMiJakctkwW6YmRwCyAZIyjChjp7u194K+l1bjISGBRJ1VnG4jusLW4+2wz\nqzazQV21pJL3m0gknllE996Q9t47mVAwl2glVSXbEosUjBKMVLoRwNTkg7uJ7g1+1wGL3X2xmTUB\ny4GOEkxNFzsCNhJJqAEikSSJ50O7CyaPXQJMdPeWMaCZySSAKneflLx+HjFuMz+5v9zMFmfRFSiS\nNXWRScVqGX9Jbk3u3uzu8zNtvQDLWsZUWpJH24SQPDaIpJXUianEJIBMte06mwqc2er92o4rDVZy\nkUJTgpFKNpYYf6kBJraXHLqQ6R7lDcD0zg5IEtX8ZFpyFdH11R2LWbOXeh1tElA7kwtE8k4JRirZ\nSGBG0q3UzJoWQM7XiHTRPdaiZbC/LovWRg1rEl7LToQiqVKCkUrWegbZ68kNcvjhnIztdNY99kGr\nKRmkH9T6sU60jbEBmJCcZzawuPXCzGRh56BM4xbJBSUYqUhJN9TiVl1HM4G+Zja6qwWNyQf1BGCQ\nmY03s6pWJWAmtmkBjWlv9lbygT81Ob51+ZhJJEkvSRDTWp239Ur+xWY22sxGJFUIZrWZRTYcGJN0\nuY0mu7U+Ij1i7p52DCJFL0kao1tmaXXjdVO7MWkg03OOBhrcfWSWr8/q3yLSXWrBiORJUqtsRtpx\niKRFCUYkf0bmun5Y0j3XANR1t9imSKFpoaVI5oZ0s8sr51ODk3GUrLrG4INW1WC6Xpcj0mMagxER\nkbxQF5mIiOSFEoyIiOSFEoyIiOSFEoyIiOSFEoyIiOSFEoyIiOTF/wMx25ymfYEiJQAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x181eeb1588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "kminn=1e-3\n",
    "kmaxx=1\n",
    "\n",
    "k1 = np.linspace(kminn,kmaxx,1000)\n",
    "\n",
    "plin = [Pk(k) for k in k1]\n",
    "\n",
    "plt.loglog(k1,plin,color='blue')\n",
    "\n",
    "plt.xlim(kminn,kmaxx)\n",
    "\n",
    "plt.xlabel(\"$k \\, \\mathrm{[h/Mpc]}$\")\n",
    "plt.ylabel(\"$P_{\\, \\mathrm{lin}}(k) \\, \\mathrm{[Mpc^3/h^3]}$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise: \n",
    "\n",
    "Plot the non-linear (HALOFIT) power spectrum and compare it to the linear one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Tonale 2018"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}