CLASS vs. CAMB.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook compares TT spectra and P(k) between CLASS and CAMB, for reference spectra parameters (i.e. very high precision). \n",
    "\n",
    "Note that these plots below were generated with a CLASS input\n",
    "```\n",
    "\"k_max_tau0_over_l_max\":15.\n",
    "```\n",
    "but values of 10ish will work too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import sys, platform, os\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "# print('Using CAMB installed at %s'%(os.path.realpath(os.path.join(os.getcwd(),'..'))))\n",
    "#uncomment this if you are running remotely and want to keep in synch with repo changes\n",
    "#if platform.system()!='Windows':\n",
    "#    !cd $HOME/git/camb; git pull github master; git log -1\n",
    "sys.path.insert(0,os.path.realpath(os.path.join(os.getcwd(),'..')))\n",
    "import camb\n",
    "from camb import model, initialpower\n",
    "# print('CAMB version: %s '%camb.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:accuracy parameters are changed globally, not yet per parameter set\n"
     ]
    }
   ],
   "source": [
    "#Set up a new set of parameters for CAMB\n",
    "pars = camb.CAMBparams()\n",
    "camb.set_halofit_version('takahashi')\n",
    "#This function sets up CosmoMC-like settings, with one massive neutrino and helium set using BBN consistency\n",
    "pars.set_cosmology(H0=67.05, ombh2 = 0.02225, omch2 = 0.1198, \n",
    "                   tau=0.079, num_massive_neutrinos=0, mnu=0.0, \n",
    "                   standard_neutrino_neff=3.046)\n",
    "pars.InitPower.set_params(As=2.2065e-9, ns=0.9645)\n",
    "pars.set_for_lmax(4000, max_eta_k=12000, lens_potential_accuracy=4);\n",
    "pars.set_accuracy(AccuracyBoost=3, lAccuracyBoost=3, lSampleBoost=3, DoLateRadTruncation=False )\n",
    "pars.AccuratePolarization = True\n",
    "pars.AccurateReionization = True\n",
    "pars.YHe = 0.24\n",
    "# pars.omegan = 0.0006445\n",
    "pars.omegak = 0.\n",
    "pars.set_nonlinear_lensing(True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total\n",
      "lens_potential\n",
      "lensed_scalar\n",
      "unlensed_scalar\n",
      "unlensed_total\n",
      "tensor\n",
      "(4042, 4)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAECCAYAAADuGCyPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XtwVOX9BvBncyNcQjaJCbcY03Ct\nIpfNC6SFCMEglOvYidARncGREDuorVOh3qZWfxYKVWRqK1XHy0jHaSV2EMQGSQW5NEHeREEKgpMA\nQgK5LwQC5LL7+2PPriFssmeTPTl7zj6fmZ2cfc/Zs99Nsk/evOecdy1OpxNERGR8YXoXQEREgcFA\nJyIyCQY6EZFJMNCJiEyCgU5EZBIMdCIik2CgExGZRITeBRDpwWKxrACQDmAXgEkAcgD8FkA8gPsA\nbOlqvdPpnKVD2URdYqBTqKp3Op15AGCxWOwAbE6nM1+5Xw8AvtYTBRsGOoWq0m6uU7OeSBcMdApJ\nTqezvDvr1Kwn0gsPihIRmQQDnYjIJBjoREQmwUCnkGaxWLLhOk0xzWKx5FgsFqs/64mCiYXzoRMR\nmQN76EREJuHztEUhhA1AGgBIKfOVthwAdgA2KeV6f9qIiEgbanroTytBniaEsCkBDyllIQC7P20a\nvQYiIoKPQFd62IcAQEq5XkpZCmAJXL1uACgHkO1HGxERacTXkMskwDPskq0Mm1gBtJ/LIsGPthvc\nfffdzqFDh3ajbCKi0HT8+PFaKWWit3VqLv2vk1KWCiGylR57wERGRuL48eOdrs/NzUVeXl4gn5KI\nyNCEEGc6W+cr0OvgGi4BXMMnk5Sv8UqbVdkGfrR5JCYmoqCgwEcJRESkhq9Az4drHmjAFcqH4Ap4\nobSlAShUltW2ERGRBro8KCqlLIfrDJUcAAlSynzlwCiEENkA7FLKUrVtmr4SIqIQp+uVog8++KBz\n8+bNuj0/EZHRCCFKpJTC2zpeKUpEZBIMdCIik2CgExGZREgG+ieffIKMjAx88cUXepdCRBQwIRfo\nbW1tWLlyJQ4ePIiVK1eC0wcTkVmEXKAfPHgQ33//PaZPn47//e9/+O677/QuiYgoIEIu0Pfv3w8A\nWLt2LQDgP//5j57lEBEFTMgFenFxMUaMGIGMjAzEx8fjq6++0rskIqKACLlAP378OMaOHQuLxYLx\n48fjyJEjepdERBQQIRXora2tKCsrw+jRowEA48aNwzfffMMDo0RkCiEV6KdPn0ZLS4sn0EeOHImm\npiZUVVXpXBkRUc+FVKCfOHECADyB/qMf/QgAUF5e3uljiIiMIqQC3X2K4qhRowAAaWlpABjoRGQO\nIRXoZ8+eRb9+/ZCQ4Po0vNTUVADAqVOndKyKiCgwQirQz507h+TkZFgsFgBAdHQ0hg4dyh46EZlC\nSAV6RUUFhg0bdkPbrbfeioqKCp0qIiIKnJAKdHcPvb2hQ4eisrJSp4qIiAInZALd4XCgoqKCgU5E\nphUygV5TU4PW1tabhlyGDBmChoYGXL16VafKiIgCI2QC/dy5cwDgtYcOAOfPn+/1moiIAinkAr1j\nD90d6Bx2ISKjC5lAv3DhAgDXEEt7DHQiMouQCfSamhoAQGJi4g3t7kDnqYtEZHQhE+jV1dWIjY1F\nVFTUDe3x8fGIiIhAdXW1TpUREQVGyAR6TU0NkpKSbmq3WCy45ZZbPD14IiKj8hnoQoh1ytcV7dpy\nhBDZQojV/rbppbq6+qbhFrfExETU1tb2ckVERIGlpoe+QghRBqAcAIQQNgCQUhYCsAshbGrbNHkF\nKnXWQwdcgc4eOhEZnZpAz5VSDleCGQCWALAry+UAsv1o001XPXQOuRCRGagJ9LQOwyZWAPXt1if4\n0aYLh8OB2tpa9tCJyNQifG0gpVwPAEKIWUKIgPaya2pqIITodH1ubi7y8vJ6/DwNDQ1oa2vrcgzd\nbrejpaUFkZGRPX4+IiI9dBnoyoHQeillPoA6AGlwDaPEK5tYlXb40eaRmJiIgoKCbhevlrv33VUP\nHQDq6uowePBgzeshItKCrx66hHIwFMBwAG8obe5udRoA99i62rZe5z7HvKseOuAKfgY6ERlVl2Po\nUspSAIuFEDkAyqSUpUoblOEXuz9tmr6SLvjqod9yyy03bEdEZERqxtDfDGSbHtxB7Q7ujtw9dJ6L\nTkRGFhJXitbXu062iY+P97q+/ZALEZFRhUSgNzQ0IDo6GtHR0V7XJyS4zqhkoBORkYVMoMfFxXW6\nPiIiAlar1dOTJyIyIga6Ij4+noFORIbGQFfExcUx0InI0EIi0O12u6oeekNDQy9VREQUeCER6Bxy\nIaJQwEBXcMiFiIzO9IHe1taGixcvqh5ycTqdvVQZEVFgmT7QL168CACwWq1dbhcfH4+2tjY0Njb2\nRllERAFn+kB3H+hUM+QCgMMuRGRYDHSFe1oAnulCREbFQFe4A509dCIyKga6gkMuRGR0DHQFh1yI\nyOgY6AoOuRCR0YVEoEdFRaFv375dbte3b1/06dOHgU5EhhUSgR4XFweLxeJzW87nQkRGFjKBrgbn\ncyEiI2Ogt8P5XIjIyEwf6GqmznXjkAsRGZnpA51DLkQUKhjo7XDIhYiMzNSB7nA4/BpyiYuLw5Ur\nV9DS0qJxZUREgWfqQL906RKcTqfPqXPdeLUoERlZhNoNhRCrpZTrleUcAHYANn/bepPaq0Td3Ns1\nNDQgKSlJs7qIiLSgqocuhMgGMEtZtgGAlLIQgF0IYVPbpsUL6Ep3A53j6ERkRN0ZclkCV68bAMoB\nZPvR1qv8DXQOuRCRkfkMdCGETellu1kBtO/CJvjR1qt6MuRCRGQ0anro8ZpXoREOuRBRKOnyoKiX\n3jngGkZxh7wVQJ2yrLbNo6amBkKITp8/NzcXeXl5XZXYJX8D3X02DHvoRGREvs5ySRNCpMEVzPHK\ngc1/AnCncBoAd+CrbfNITExEQUFBN0v3raGhAREREejfv7+q7SMjIxETE8NAJyJD6nLIRUqZL6XM\nV+5albZSwHPmi11KWaq2TasX0Rl/ps51i4uLY6ATkSGpOg9dSvkmgDc73Pe2jc+23uTPZf9uvPyf\niIzK1FeK+nPZvxt76ERkVKYO9O700DmFLhEZFQO9A/bQicioGOgdcAydiIzKtIHudDq7PYZ+7do1\nXLt2TaPKiIi0YdpAb2xsRFtbm+qpc904nwsRGZVpA93fq0TdOJ8LERkVA70DzudCREbFQO+APXQi\nMioGegccQycio2Kgd8AeOhEZFQO9g9jYWAAcQyci4zF1oIeHhyMmJsavx4WHhyM2NpY9dCIyHFMH\nutVq9WvqXDfO50JERmTaQO/OVaJuvPyfiIzItIHenXlc3DhBFxEZEQPdCwY6ERkRA90LjqETkREx\n0L1wj6E7nc4AV0VEpB1TBrrT6UR9fX2PAr2lpQVNTU0BroyISDumDPTLly+jra2tR4EO8GpRIjIW\nUwZ6d68SdeN8LkRkRKYOdHcw+4tT6BKREZk60DnkQkShxJSB7u5ZM9CJKJSYMtA5hk5EoSjC1wZC\niGxlcZaU8rdKWw4AOwCblHK9P229oaeBHhMTg7CwMI6hE5GhdNlDV8L8PillIQCbEMImhLABgNJm\n96dN01fSTnenznULCwuD1WplD52IDKXLQJdSFkop85S7aVLKUgBL4Op1A0A5gGw/2nqF+yrR7kyd\n68bL/4nIaFSNoQshVgNwB7sVQPuxiAQ/2npFT64SdeMEXURkND7H0AFASrleCLFFCCED+eQ1NTUQ\nQnS6Pjc3F3l5eZ2u70xP5nFx45zoRGQ0XQZ6u3HwUriGTVbANYzivmLHCqBOWVbb5pGYmIiCgoLu\n1t6phoaGbl9U5BYXF4dTp04FqCIiIu356qFnAyhVlq0ADgEoBODuVqcp9+FHm+YaGhowYsSIHu2D\nY+hEZDS+xtDfBJAmhFgBAFLKfKW37j4Dxi6lLFXbptmr6CBQQy4NDQ2cQpeIDKPLHrqU0g5XqHds\n73ab1hwOR8ACva2tDY2NjRg4cGCAqiMi0o7prhRtbGyEw+EISKADvFqUiIzDdIHe05kW3Xj5PxEZ\njWkDPVA9dJ66SERGYbpA7+lMi24cciEiozFdoAeqh37LLbcAAGpra3tcExFRb2CgdyIxMRGA62pW\nIiIjMG2g9/SgaFRUFKxWK6qrqwNRFhGR5kwZ6JGRkejXr1+P95WYmMhAJyLDMF2gu2da7MnUuW5J\nSUl+B/q1a9ewZMkSJCcn49133+1xDUREapku0Ovq6pCQEJiZersT6M888ww+/PBDxMbGYvny5fjy\nyy8DUgsRkS+mC/Ta2lrPGSo95W+gV1VVYdOmTXjooYdQXFyM+Ph4/O53v+txHQ0NDbj//vsxZcoU\n7Nq1q8f7IyJzMl2g19XVBTTQ6+rq0NbWpmr7Dz74ANeuXcOqVasQExODxx9/HDt37uzRNLxOpxNL\nlixBfn4+qqqqsGDBApw8ebLb+yMi8zJdoAe6h+5wOFRfLbplyxaMHz8eP/7xjwEAy5Ytg8Viwfvv\nv9/tGrZt24Zdu3Zhw4YNKC4uRlRUFJ555plu74+IzMtUge50OlFbWxvQMXQAqoZdqqqqUFRUhPvu\nu8/Tduutt2LatGnYunVrt2vYuHEjUlNT8cgjj2Dw4MH45S9/ia1bt6KysrLb+yQiczJVoDc2NqK1\ntTVgPXT3xUVqAr2oqAgAMHPmzBva582bh6+//hoVFRV+P//JkyexZ88erFixAhERrpmOc3Nz0dbW\nhs2bN/u9PyIyN1MFuvsy/UAOuQDqAr24uBiRkZGYOHHiDe3z588HAOzYscPv59+8eTPCwsKwbNky\nT9uIESNgs9mwbds2v/dHROZmykDXY8iluLgYEydORHR09A3tt99+O1JSUrBz506/n3/btm2YNm0a\nhgwZckP7ggULUFRUxHlmiOgGpgz0QPXQ4+PjERYW5nM+l9bWVhw6dAgZGRk3rbNYLMjKysLevXvh\ncDhUP/eZM2dw5MgRLFy48KZ18+fPh9Pp7NYfCSIyL1MFel1dHYDABXp4eDiSkpJ8HoA8evQompqa\nvAY6AMyYMQO1tbU4duyY6ufevn07AFdvvKOJEyciNjYWe/fuVb0/IjI/UwV6oIdcACA5OdnnAc3i\n4mIA6DLQAWDPnj2qn3f79u0YPXo0Ro0addO68PBw/PSnP8W+fftU74+IzM90gR4eHo7Y2NiA7TM5\nORnnzp3rcpuioiIkJSUhNTXV6/rU1FTcdtttqgO9qakJe/bswbx58zrdJjMzE8ePH+c4OhF5mC7Q\nExISEBYWuJelJtCLi4uRkZHR5YRgM2bMwBdffKFqHH3v3r1obm7G7NmzO91m2rRpAID9+/f73B8R\nhQZTBXogL/t3S05Oht1ux+XLlzt9zpMnT3Y63OKWlZWF2tpaHD161OdzfvbZZ+jTpw8yMzM73WbS\npEmIiIjg5F9E5GGqQA/kVaJuycnJANDpOLo7UNUEOgDs3r3b53Pu2rULmZmZ6Nu3b6fbREdHY+zY\nsZBS+twfEYUG0wW6Fj10AJ0OuxQXFyMsLAxCiC73k5KSguHDh+Pzzz/vcrvKykocPXoU99xzj8/a\n0tPTUVJSAqfT6XNbIjI/UwW6VkMuQNeBPnbsWMTExPjc18yZM7Fnzx60trZ2uk1hYSEAYNasWT73\nl56ejvr6epw5c8bntkRkfhG+NhBCrFAWh0spf6u05QCwA7BJKdf706YVh8OhSQ992LBhAICzZ896\nfc6DBw/iF7/4hap9zZw5E2+99Ra++uorTJo0yes2n332GZKSkjBu3Dif+0tPTwcAlJSUdHqGDRGF\nji576EKIbACFUso3AaQJIbKFEDYAkFIWArALIWxq27R8IfX19WhtbcWgQYMCut/o6GgMGjQIp0+f\nvmndiRMncPHiRZ/j526+xtGdTicKCwsxa9YsVWfqjBs3DhERESgpKVH1/ERkbr5SIw1AtrJcrtxf\nAlev292W7UebZqqqqgAg4IEOACNHjvT6oRK+LijqaNCgQbjjjjs6HUcvKSlBVVWVqvFzgAdGiehG\nXQa6lPJNpXcOADYAEoAVQPtPfEjwo00zWgb6qFGjvAZ6UVERrFar16s5O5OVlYV9+/bh+vXrN637\n+OOPER4e3uUFRR2lp6ejtLSUB0aJSN1BUWW4pFRKWRrIJ6+pqYEQotPbG2+8oXpfWgd6VVUVLl26\ndEP7gQMH8JOf/MSvC5nmzp2LpqYmr58NunXrVmRmZvp16qXNZkNdXZ3XMX4iCi0+D4oqst0HROEa\nRolXlq0A6pRltW0eiYmJKCgo8Kvgzmg95AIA3333nedAZH19PY4dO4b777/fr33dfffdiI2NxUcf\nfeSZKx1wjccfPXoUr776ql/7s9lchya++uorpKSk+PVYIjIXn11LIcSKdmetZAP4J1xj6VC+FvrR\nppmqqipERkYiLi4u4Pt2D6l8++23njb3+PnUqVP92ldUVBQWLVqErVu34tq1a572d955B+Hh4arP\nmHEbN24cwsLCUFoa0H+eiMiA1Jzlsk4IUSaEaAAA97CLss4upSxV26blC6mqqkJSUlJA53FxGz16\nNKKjo284m+TAgQOIiIjA5MmT/d7fsmXLYLfb8fe//x2AazKu9957D/Pnz8fgwYP92le/fv0wZswY\nBjoRdT3kopxyeFOXt92BUr/btFJVVaXJcAsAz0fLHTp0yNP22WefYfLkyejXr5/f+5sxYwYmTpyI\nP/7xj7j//vvx6quvorq6Gk8++WS36rPZbKqmFCAiczPNlaJaBjoATJ48GSUlJWhpacH58+chpfTr\nbJT2LBYLXn75ZZSVlWHKlCn4/e9/j5ycHM8Miv6aOHEiKioqPMcRiCg0mSbQL1y4oGmgZ2Vl4erV\nq9i9eze2bt0KAN0OdMB11ei7776L1tZWLF68GG+//Xa399X+wCgRhS61Z7kENafTierqak0D/Z57\n7kH//v3x3nvv4fDhw5gwYYKqy/O7smzZMixbtqzHtU2YMAGAK9DnzJnT4/0RkTGZItAbGhrQ0tKi\naaD37dsXeXl52LBhAwBgy5YtXX6gRW+yWq1IS0vjgVGiEGeKQNfyHPT2XnrpJQwcOBApKSnIycnR\n9Ln8ZbPZGOhEIc4UgX7hwgUA8PuUP3/17dsXzz//vKbP0V02mw35+fmw2+2wWq16l0NEOjDFQVH3\npwm5p7oNRe4Do19//bXOlRCRXhjoJjFx4kQA4LALUQgzTaAPHDgQAwYM0LsU3SQlJWHYsGEMdKIQ\nZppAD+XeuZvNZuO56EQhjIFuIhMnTsS3336LK1eu6F0KEemAgW4iNpsNDocDR44c0bsUItKB4QO9\nra0N58+fZ6DjhzNdOI5OFJoMH+jV1dVoa2tjoANITk5GQkICx9GJQpThA919ymJycrLOlejPYrHw\nilGiEGaaQGcP3cVms+Ho0aNeP4SaiMyNgW4yGRkZaGlpwZdffql3KUTUy0wR6BEREUhKStK7lKAw\nffp0WCwWfP7553qXQkS9zPCBfvbsWQwdOlSTzxI1ori4ONhsNgY6UQgyfAqePn0aqampepcRVLKy\nslBcXIympibVj2lpaUFbW5uGVRGR1hjoJpSdnY3m5mZVvfSysjLMmzcPffr0gdVqxerVq9Hc3NwL\nVRJRoBk60Jubm1FRUcFA7yArKwuxsbHIz8/vcruDBw9CCIGioiI8+eSTWLRoEf70pz9h0aJFaGlp\n6aVqiShQDP0BF+fOnYPD4WCgdxAVFYVFixbh448/RnNzM6Kiom7a5vTp01i4cCHi4+Oxa9cupKWl\nAQAyMzPxyCOP4Omnn8bLL7/c26UTUQ8Yuod++vRpAGCge7F48WLY7XZs3br1pnUXL17E/Pnz0dzc\njB07dnjCHADy8vLwyCOPYMOGDSgqKurNkomoh1QFuhDC1uF+jhAiWwix2t+2QGKgd27OnDkYPnw4\nNmzYAKfT6WlvbW3FkiVLcOLECXz00UcYM2bMTY9dv349hg0bhsceewwOh6M3yyaiHvAZ6EKIbABb\n2t23AYCUshCAXQhhU9sWqKL37duH/fv348yZMwgLC+Nl/16Eh4dj1apVOHjwIP7xj38AABwOBx59\n9FHs3LkTmzZtwsyZM70+NiYmBn/4wx9QUlKCLVu2eN2GiIKPz0BXArm8XdMSAHZluRxAth9tAXHX\nXXchMzMTp0+fRnJyMiIjIwO1a1NZvnw5pkyZguXLl2PdunVYsGAB3njjDTz11FNYvnx5l49dunQp\n7rzzTjz77LM8QEpkEN0ZQ7cCqG93P8GPtoAqKyvjcEsXwsPD8fHHH2P8+PF46qmnsHfvXmzcuBFr\n1qxR9di1a9eirKwM77zzTi9US0Q9ZeiDoidPnsTo0aP1LiOoDRo0CAcOHMDZs2dRW1uLX/3qV7BY\nLKoeO3fuXEydOhUvvvgirl69qnGlRNRT3Tlt0Q4gXlm2AqhTltW2edTU1EAI0ekT5ebmIi8vr9P1\nNTU1GDVqlOrCQ5XFYunWcQaLxYI1a9Zg+vTp+Otf/4onn3xSg+qIKFC6E+j/BOBO4TQAhcqy2jaP\nxMREFBQUdKOEH7CHrq277roLc+bMwdq1a5Gbm4vY2Fi9SyKiTqg5yyXH9UXkAICUslRpzwZgl1KW\nqm3T4gWwh669P/zhD6ivr8eGDRv0LoWIumBpf45yb3vwwQedmzdv9vtx7ceAm5ubeZZLL1i8eDH+\n/e9/o7y8HImJiXqXQxSyhBAlUkqvY9WGPiialJTEMO8lL774IpqamrB27Vq9SyGiThgy0MPDwwEA\n8+bN07mS0DFmzBgsW7YMr7/+Os6ePavqMRUVFdiyZQs++OADnDhxQuMKiciQgT5p0iQAwKZNm3Su\nJLQ8//zzsFgsePjhh7ucO/3IkSOYO3cukpOTsXjxYixduhRjxozB7NmzcerUqV6smCi0GDLQAeCe\ne+5Bnz599C4jpKSkpODPf/4zdu3ahd/85jfoePzlwoULePjhhzFhwgQUFxfjhRdeQElJCb755hus\nW7cORUVFsNls+O9//6vTKyAyN0NOn6vngdxQt3z5chw7dgwbN25EeXk5nnjiCfTp0wfbtm3D66+/\njmvXruGJJ57Ac889h7i4OM/jxo4di5ycHMyZMwezZ8/Gvn37MGHCBB1fCZH5GDLQAai+2pECy2Kx\nYMOGDUhJScFzzz2H7du3AwDCwsJw7733Ys2aNZ2eSpqWloY9e/YgIyMDc+fORUlJCYYMGdKb5ROZ\nmiGHXNhD15fFYsETTzyB8+fP49NPP8Unn3yCiooK5Ofn+7wuYOjQodixYwfsdjuWLl3KzzElCiBD\nBjrAHnowGDhwIH72s59h3rx5GDx4sOrH3XnnnfjLX/6C3bt346WXXtKwQqLQYshAZw/d+B566CE8\n8MADePHFF7Fnzx69yyEyBUMGOsAeutFZLBZs2rQJw4cPx9KlS1FTU6N3SUSGZ8hAZw/dHAYMGIAP\nP/wQtbW1WLZsGT/ujqiHDBvo7KGbw4QJE/DKK6/g008/xauvvqp3OUSGZshABzjkYiYrV67Evffe\ni6eeegoHDhzQuxwiwzJkoHPIxVwsFgvefvttpKamYv78+fjmm2/0LonIkAwZ6AB76GYTFxeHXbt2\noV+/frj77rs5PQBRNxgy0NlDN6fU1FTs3r0bsbGxmDFjBp555hnU1d30yYVE1Ale+k9BZdSoUTh4\n8CB+/etfY+3atXjllVcwefJkpKamIjw8HPX19aiqqkJ1dTWuXr0Kh8OBxMREjB49GllZWbj33nu7\n9fmpRGbAHjoFnfj4eLz//vs4fPgwHnvsMTidTuzfvx+7d+/G999/D6vVimnTpmHRokX4+c9/jlGj\nRuHrr7/G448/jttuuw0LFy7kkA2FJPbQKWiNGzcOL7/8surtv/vuO7z77rt46623MHXqVMyePRsv\nvPACpkyZomGVRMGDPXQyjZEjR2LNmjU4ffo01q1bByklMjIyMG/ePBw6dEjv8og0Z8hAB9hDp871\n798fq1evxqlTp7BmzRoUFxdj8uTJmD9/Pj755BO0tLToXSKRJgw55MIeOqkRExODp59+GitXrsRr\nr72GjRs3YseOHUhISMD06dM9B1sHDx6MPn36ICoqCpGRkZ6vkZGRiI+PR//+/fV+KUSqGDLQAfbQ\nSb2BAwfi2WefxapVq1BQUID8/Hzs378f//rXv1Q9PikpCSNHjkRGRgamTp2Ku+66CwkJCRpXTeQ/\nQwY6e+jUHVFRUVi4cCEWLlwIALh06RLOnj2LqqoqNDc34/r162hpaUFzc7Pna21tLU6dOoVjx47h\ntddewyuvvIKwsDBMmzYNCxYswMKFC31+qAdRbzFsoLOHTj01cOBA3HHHHbjjjjtUbX/t2jVIKbFz\n505s27YNq1atwqpVqzBixAjMmjULs2bNQlZWFqxWq8aVE3lnyEAHOORCvS86OhrTpk3DtGnT8H//\n9384c+YMtm/fjoKCArz//vvYtGkTANdnp44fPx6jRo3CkCFDMGTIEMTHx2PAgAGeW0xMDGJjYxER\nYdi3IAUhQ/42cciFgsFtt92GRx99FI8++iiam5tx8OBB7N27F4cPH8bhw4dVnVHTv39/WK1WxMXF\nYdiwYUhOTvbchg0bhmHDhmHIkCFISEhAWJhhT0qjXqJpoAshcgDYAdiklOsDuW/20CmYREVFITMz\nE5mZmZ42p9OJ+vp6VFZWoqGhAVeuXMHly5dx5coVXLp0CRcvXoTdbofdbkddXR0qKytx5MgRXLhw\n4aZOS2RkJAYPHuzp8VutVk9P3/21f//+N90GDBhww/2+ffvyvWNimgW6EMIGAFLKQiFEmhDCJqUs\nDcS+2UMnI7BYLEhISPD7jJiWlhacP38eZ8+eRWVlJc6fP++5VVZWory8HBcvXsTly5dx+fJlNDc3\n+1XTgAEDMHDgwBtuMTExXtsGDBiAqKgoREVFeU7t9Hbzti4yMpJ/PHqZlj30JQB2KcvlALIBBCTQ\ny8rKcPvttwdiV0RBJzIyEikpKUhJSVG1fXNzMxobG3HlypUub+7/DhobG9HY2IhLly55bpWVlTe0\nBerjACMjIxEREYGIiAjPspqiM5t1AAAEg0lEQVS28PBwhIWFeb6qWe6tbS0WS49vERERWLRoUUC+\nx+1pGehWAPXt7gfsxN3r16+jX79+gdodkaFFRUV16z+BzjidTjQ1NaGxsdHzH0D72/Xr1/1qa21t\n9dxaWlpuWu7Y5r7f1tYGh8MBh8Phc7m76/X6HNt+/frhypUrAd+vrgdFa2pqIITodH1ubi7y8vJu\nav/000+Rnp6uZWlEIctisXjG3M3O6XTC6XT69cfB4XB4Htfdm2ZDUT0trLNbenr6uvT09GxlOSc9\nPX11x20eeOABZ3f97W9/6/ZjtRSMdQVjTU5ncNYVjDU5nazLH8FYk9MZuLrS09Ols5Pc1fI8qH8C\nSFOW0wAUBnLnb731ViB3FzDBWFcw1gQEZ13BWBPAuvwRjDUBvVOXZoHuPqNFCJENwB6oM1yIiMg7\nTcfQpZRvarl/IiL6AS89IyIyCQY6EZFJMNCJiEyCgU5EZBIWp47zogghagCc6ebDhwKoDGA5gRKM\ndQVjTUBw1hWMNQGsyx/BWBMQuLpuk1Imeluha6ATEVHgGHI+dLpRx5ksvU1brLZN47rWSSl/K4RY\n4T6lVY+6iLQihFjt73sukL/vhhtDF0LkCCGyhRCrdaxhnfJ1RVd19UatyoVbW9rd90xbDMAuhLCp\nbdOyLsUKIUQZXLNvqq41wHWtUG7r2rWp+tlp9fPspCbdf8eU/WYH2ffKW026f6/ctQGYpSzr8j40\nVKBr/Wb3g+7B5Kbsv7xd0xK4/toDP0xbrLZNy7oAIFdKOVxZ50+tAaG84QqV/w7SlDe3rn8AvdWk\nrNL1d0yp4z5l3zZ/vi8af69uqElZFTTvx3Z0eR8aKtChcQj5Qddg8sHbtMVq27TmDlF3D6m360rD\nDz+HcuW+3n8AvdUE6Pw7JqUslFK6pzpNU4bOdP1edVITEATvR+EaXmw/X5Uu70OjjaHrEULeuHtS\n7jGvYAnMoNZuzHBWu55obz5/+6kobHBNIJcOHf8AdlITECS/Y8ofX3eIBkVnoUNNQHB8r+IDvL9u\nMVoPPShIKdcrf40T9AgmH+z44ZfLCqDOjzbNKGPEOcrdOrh6orrUpfy7XRpME8Z1rClYfseUgMwT\nQlj1qqGjjjXp/b3y0jsHdHofGq2H3qsh5I1y4KVeSpmProMJnbRp7Z8A3J8a0n7aYrVtWpH4YUx9\nOIA3lDY96sqWUv5WWfbnZ6flz9NTUzD8jrUbcy6F6+e2ws8aAl6Xt5qEEHbo/35ME0KkKfuPV+rU\n5X1otB66pnOsqyTbPe9w5b63unqlVqXXK9y9X2/TFqtt64W6Fiv3y3Ssa0W7oZ9sqP/Zafbz9FJT\nMPyOZePGACz3owat6vJWk+7fKyllvvIHxV2Xbu9Dw11YpPReyuE6KKLL9LzuHpRSw/p2bTfUFQy1\n0g/ED6dS1sMVDPdJKQvV/uy0+Hn6qEm33zFlOGOxcjfdfTBS5+9VVzXx/QgDBjoREXlntCEXIiLq\nBAOdiMgkGOhERCbBQCciMgkGOhGRSTDQiYhMgoFORGQSDHQiIpP4f94doE7f4ARYAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc45801b950>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "#calculate results for these parameters\n",
    "results = camb.get_results(pars)\n",
    "\n",
    "#get dictionary of CAMB power spectra\n",
    "powers =results.get_cmb_power_spectra(pars, CMB_unit='muK')\n",
    "for name in powers: print(name)\n",
    "\n",
    "#plot the total lensed CMB power spectra versus unlensed, and fractional difference\n",
    "totCL=powers['total']\n",
    "print(totCL.shape)\n",
    "#Python CL arrays are all zero based (starting at L=0), Note L=0,1 entries will be zero by default.\n",
    "#The different CL are always in the order TT, EE, BB, TE (with BB=0 for unlensed scalar results).\n",
    "ls = np.arange(totCL.shape[0])\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(ls,totCL[:,0],'k-')\n",
    "ttCAMB = totCL[:,0]\n",
    "ax.set_title('TT');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from classy import Class\n",
    "\n",
    "# Define your cosmology (what is not specified will be set to CLASS default parameters)\n",
    "params = {\n",
    "    'output': 'tCl lCl mPk',\n",
    "    'l_max_scalars': 4000,\n",
    "    'lensing': 'yes',\n",
    "    'non linear': 'halofit',\n",
    "    'omega_b' : 0.02225,\n",
    "    'omega_cdm' : 0.1198,\n",
    "    'h' : 0.6705,\n",
    "    'ln10^{10}A_s' : 3.094,\n",
    "    'n_s' : 0.9645,\n",
    "    'tau_reio' : 0.079,\n",
    "    'N_eff': 3.046,\n",
    "    'YHe' : 0.24,\n",
    "    'N_ncdm' : 0,\n",
    "    'halofit_k_per_decade' : 3000.,\n",
    "    'l_switch_limber' : 40.,\n",
    "    'accurate_lensing':1,\n",
    "    'num_mu_minus_lmax' : 1000.,\n",
    "    'delta_l_max' : 1000.,\n",
    "\n",
    "    \"recfast_Nz0\":100000,\n",
    "    \"tol_thermo_integration\":1.e-5,\n",
    "    \"recfast_x_He0_trigger_delta\" : 0.01,\n",
    "    \"recfast_x_H0_trigger_delta\" : 0.01,\n",
    "    \"evolver\":0,\n",
    "    \"k_min_tau0\":0.002,\n",
    "    \"k_max_tau0_over_l_max\":10.,\n",
    "    \"k_step_sub\":0.015,\n",
    "    \"k_step_super\":0.0001,\n",
    "    \"k_step_super_reduction\":0.1,\n",
    "    \"start_small_k_at_tau_c_over_tau_h\" : 0.0004,\n",
    "    \"start_large_k_at_tau_h_over_tau_k\" : 0.05,\n",
    "    \"tight_coupling_trigger_tau_c_over_tau_h\":0.005,\n",
    "    \"tight_coupling_trigger_tau_c_over_tau_k\":0.008,\n",
    "    \"start_sources_at_tau_c_over_tau_h\" : 0.006,\n",
    "    \"l_max_g\":50,\n",
    "    \"l_max_pol_g\":25,\n",
    "    \"l_max_ur\":50,\n",
    "    \"tol_perturb_integration\":1.e-6,\n",
    "    \"perturb_sampling_stepsize\":0.01,\n",
    "    \"radiation_streaming_approximation\" : 2,\n",
    "    \"radiation_streaming_trigger_tau_over_tau_k\" : 240.,\n",
    "    \"radiation_streaming_trigger_tau_c_over_tau\" : 100.,\n",
    "    \"ur_fluid_approximation\" : 2,\n",
    "    \"ur_fluid_trigger_tau_over_tau_k\" : 50.,\n",
    "    \"ncdm_fluid_approximation\" : 3,\n",
    "    \"ncdm_fluid_trigger_tau_over_tau_k\" : 51.,\n",
    "    \"l_logstep\":1.026,\n",
    "    \"l_linstep\":25,\n",
    "    \"hyper_sampling_flat\" : 12.,\n",
    "    \"hyper_sampling_curved_low_nu\" : 10.,\n",
    "    \"hyper_sampling_curved_high_nu\" : 10.,\n",
    "    \"hyper_nu_sampling_step\" : 10.,\n",
    "    \"hyper_phi_min_abs\" : 1.e-10,\n",
    "    \"hyper_x_tol\" : 1.e-4,\n",
    "    \"hyper_flat_approximation_nu\" : 1.e6,\n",
    "    \"q_linstep\":0.20,\n",
    "    \"q_logstep_spline\": 20.,\n",
    "    \"q_logstep_trapzd\" : 0.5,\n",
    "    \"q_numstep_transition\" : 250,\n",
    "    \"transfer_neglect_delta_k_S_t0\" : 100.,\n",
    "    \"transfer_neglect_delta_k_S_t1\" : 100.,\n",
    "    \"transfer_neglect_delta_k_S_t2\" : 100.,\n",
    "    \"transfer_neglect_delta_k_S_e\" : 100.,\n",
    "    \"transfer_neglect_delta_k_V_t1\" : 100.,\n",
    "    \"transfer_neglect_delta_k_V_t2\" : 100.,\n",
    "    \"transfer_neglect_delta_k_V_e\" : 100.,\n",
    "    \"transfer_neglect_delta_k_V_b\" : 100.,\n",
    "    \"transfer_neglect_delta_k_T_t2\" : 100.,\n",
    "    \"transfer_neglect_delta_k_T_e\" : 100.,\n",
    "    \"transfer_neglect_delta_k_T_b\" : 100.,\n",
    "    \"neglect_CMB_sources_below_visibility\" : 1.e-30,\n",
    "    \"transfer_neglect_late_source\" : 3000.\n",
    "#     'm_ncdm' : 0.06\n",
    "}\n",
    "cosmo = Class()\n",
    "cosmo.set(params)\n",
    "\n",
    "cosmo.compute()\n",
    "cls = cosmo.lensed_cl(4000)\n",
    "ttCLASS = cls['ell'] * (cls['ell']+1) * cls['tt'] * (1e6 * 2.7255)**2 / (2*np.pi)\n",
    "# cosmo.struct_cleanup()\n",
    "# cosmo.empty()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# High ell TT Comparison\n",
    "Here, we compare the CAMB and CLASS high-ell TT. Pretty good!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/zequnl/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in divide\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEPCAYAAABV6CMBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XlwFGX6B/DvSyAQuSYhUQgsx+Cx\nuMJCaPB2VxNU8JYA64G1qBBwA1Z5RXdR94fXIrh44LoEOTxZVtTiKI0SVKBYBZooEFAwJHIEyUGY\nXEaO8P7+mO6hk/RcSXp6ZvL9VFH09Dvd/UxPpp9+j+4WUkoQEREFqp3dARARUWRh4iAioqAwcRAR\nUVCYOIiIKChMHEREFBQmDiIiCkp7uwOg0BJCTAEwHMBaACMApAPIApAAYByAD3yVSylH+Vh3GoAU\nAIUAXAAqAEyQUmZZ9XmCJYTYJ6UcGOJtpgBQAIySUo4zKQ9ovwkhPmi8vBDCCff3o0gph5uUbQOQ\nJaXM1l7PBuAE8IL2thEA1kopc1v+SZt8LieABdr287y8x/Sza8t5YpVSrvCxnSb7RZs/RZus0P5P\nkFJm+ysj/5g42p4KKWUGAAghXABS9B+lEKICAPyVm9F+iAONBzvtwJFm1QdpjlAnDc0TACYD+G/j\ngkD3mxAiHUC6EMIhpXTp86WUhUKID7T3pDQ6QDvg/r6zDe9dDncC0w/EK4QQx4QQA4zrbQ3a9vLg\nPulowtdnl1JmCSHWAhjuJ2mY7hctWTuklC8a1pvhr4wCw6aqtsf0zM9Q5q/cmwWNz5CllIUAlgcR\nm6WEEE7tQBNqDimly8uBOdD9lgBgBYApMLcA7gTVHBVwn9mHWmv8zXjbLwlw16bM1uurjALAxNHG\naD8Sr2X+ys3ma80NpmWGs7opQog07V+6vpwQYp9h/gIhRIr2b7Z2Zqi/75j2v17mNG6/8fxG6/5A\nCOGA+wA5Wwjh0N6fri9nWFeDOL3E6PCyH8w+YwqABG1dzkbv97vftPfpsS+AlzNjrabh1GMzqX2Y\n0uJz+WpK8rNvm+wT436Al4QU6Gf3E7vX/aI3vWlxLhBCpOmf0VcZBYaJgyynH0SllLnaj3aEdmDL\nhfvgUahN74O7rT4P7jPACfpy+nu0sgVw98HoMgzLZBmXgbupZpzhjF8/WOnr1pczjRPuA1OTGIP4\njHnaOnJ9JWU/xkspV2jrTdATqonlAMYHsL4ELRmkwX3Anezjvb72bZN9oq3TYdgPzf3MgfC5X7R+\nj1Fw9/MsMPRr+Cwj/5g4qMW0H66vM8tRaHgAOYqGB98KL9OmZ/baAdjYbp6lHbgbH9ATvJ1Jak0k\nTiHENpw5W/UVp7+4/H1Gsxj87TfdQK3Gkg4gF97b47MBZGhn4r76KyqklHnawT0DwEIfTXje9i1g\nvk8a7wdTQXx2X2Ve94ueRLRadLbWt5Xhr4wCw8RBrSVDCLHAOMPQfLENDQ8SAwGozd2Qtt5CbToN\n7lE7K+A+eOidnf7WMUVK+aI2EsklzoxAam6czV3W137TD3LLtTPrFXDXDkxrFYYa1fggazeFMLT5\nG7bdnH27tdG6TJO/xudn90JvvvS3X5zGJKOtVw2gjALAUVVtlPbDGQf3jygd7qYUV6DljWnDPdOE\nEI/hzNBKY3vyY3ozBoBtUso87cfvBDBeCPFfaGerQohcuJuSFCGEUz8IanEUwj3iSB9+WaGV6c0U\nCTjT1u/UEkS24T1OuA8wDm19LgD7tG1kN45TW6dpjI1GNzVZttFnnAKgSdu9r/2mLbsQ7qGzes1J\n72fQh6sC7n6bXK0WtUDbvr6/PPtAO+hPMHx2z/oad1IHsm+97JMVwj0IIQ1nOt0zhBBq478fP5/d\nqa1X//vTh4vnBbFfGn9O42f0VUZ+CMnbqlMEEEJsk42uUyAie7CpisKeftbuox2eiEKINQ4iIgoK\naxxERBQUJg4iIgoKEwcREQUl6ofjpqamyuTkZLvDICKKKN9//325qqpJZmVRnziSk5Pxzjvv2B0G\nEVFEURRlv7cyNlUREVFQmDiIiCgoTBxERBQUJg4iIgoKEwcREQWFiYOIiILCxEFEREFh4iAiikJZ\nWVmYN2+eJetm4iAiijJfffUV3n33XWzYsMGS9TNxEBFFmSVLlgAA0tK8Pr69RZg4iIiiTKdOndCj\nRw9MmzbNkvXbeq8qRVH0Zz6nqKra5HnMZuVe5unPdYaqqitCFD4RUViqrq5G37590a6dNXUD22oc\n2sEeqqrmAnDpr32V+1jmCS1hOBuvh4ioramurkbXrl0tW7+dNY4JANZq04UA0gDk+Snv0XieoihO\nAFsBwKzWQkTU1tTW1qJnz56Wrd/OPg4HgArD6x4BlJvNGwGgh1YjecyKQImIIkl1dTW6dOli2fqj\npXP8qKqqeYCnD4SIqM2qqamJ2qYqF4AEbdoB4GiA5WbzCg3LjADg6SAvKyuDoiheg5g8eTIyMjKa\nET4RUfiRUqKmpgadO3e2bBt2Jo7lAPQjuhNALgAoiuJQVdXlrdxknguAXstwQOvv0CUlJSEnJ8eK\n+ImIwk5dXR1Onz5taY3DtqYqQ9NSGgCX/hrAOm/lXuYVwj3CKh1ADw7HJaK2rLq6GgAs7eOw9ToO\nVVWzTeYN91Puax6TBhG1aTU1NQCsTRzR0jlORERg4iAioiCFoqmKiYOIKIqUlpYCcA8MsgoTBxFR\nlKitrcX06dMBAMnJyZZth4mDiChK7N692zNt5XUcTBxERFGiqqoKABAbG2vpdpg4iIiihMvlAgDk\n5ub6eWfLMHEQEUWJyspKAEB8fLyl22HiICKKEnqNo1u3bpZuh4mDiChKuFwudO3aFe3bW3tTECYO\nIqIoUVlZie7du1u+HSYOIqIoUVVVZXkzFcDEQUQUNax+1riOiYOIKEpUVVUxcRARUeCqq6vZVEVE\nRIFjjYOIiAK2ceNGz3BcqzFxEBFFgT/96U8ArL/4D2DiICKKKidOnLB8G0wcREQR7vjx457piRMn\nWr49Jg4ioginPy72ueeeQ2JiouXbY+IgIopw+l1xQ9G/ATBxEBFFPP0BTkwcREQUEL2piomDiIgC\nwqYqIiIKil7jCMXFfwATBxFRxNP7OELxLA6AiYOIKOJVVVWhXbt26Ny5c0i2Z+3zBYmIyFL5+fl4\n5ZVXAABCiJBskzUOIqII9uqrr4Z8m0wcREQRTEoZ8m0ycRARRTB9RNW8efNCtk0mDiKiCFZVVYVr\nrrkG48ePD9k2mTiIiCJYZWVlyIbh6pg4iIgiWGVlZciuGNcxcRARRSgpJaqqqkKeOGy9jkNRlHQA\nLgApqqq+GEi5r2UURXnMbD1ERNGotrYW9fX1baepSlGUFABQVTUXgEt/7avc1zKKoqQBGBWq+ImI\n7FRXV4c77rgDQOhubqizs6lqAtw1BwAoBJAWQLm/ZYiI2oRdu3YhLy8PANCjR4+QbtvOxOEAUGF4\n3fiTm5WbLqMoSopWCyEiahNqamo807179w7ptqOlczzB7gCIiELJ5XJ5pkOdOOzsHHfhzAHfAeBo\ngOUN5vmrbZSVlUFRFK9BTJ48GRkZGUGGTkRkL2PiiI+PD+m27UwcywHoR3QngFwAUBTFoaqqy1u5\nyTynoihOuBNKgpZI8vSNJCUlIScnx9IPQkQUavpT/woKCkJ2V1ydbU1V+sFdGw3lMhzs13kr9zJv\nhaqqK7RlHSH9EERENnG5XIiLi0NcXFzIt23rdRyqqmabzBvup7zJPMN80zIiomhTWVkJh8Oec+Vo\n6RwnImozioqKsG3bNiYOIiIKzNVXX42CggImDiIiCszJkycBwJb+DYCJg4goYpWXl9uyXSYOIqII\ncurUKc/0sGHDbImBiYOIKIIcPeq+FnrChAmYNWuWLTEwcRARRRC9eSo1NRXt29tzRQUTBxFRBNFr\nHImJibbFwMRBRBQhNm3a5HkGR6hvpW7ExEFEFCFef/11zzRrHERE5Fe/fv0806F+XKwREwcRUYSo\nrq72TIf6jrhGTBxERBGivLwcXbp0wX/+8x9b42DiICKKEOXl5bj88stx5ZVX2hoHEwcRUYQ4evSo\nrZ3iuiaJQwhxvxCif+hDISIiM3V1dbj88stRWlqKhIQE/wtYzKzGUSil/CnUgRARkbni4mL89NNP\nAIBzzjnH3mBgnjhGGWscQohuIYuGiIia0K8WB4CePXvaGImb2Y1OKgBMFUI4tdcSwITQhUREREZl\nZWWe6fPOO8/GSNzMEodLSjlHfyGEmBzCeIiIqBH9xoarVq3Cueeea3M0Xvo4hBBD9RdSyoUhjIeI\niBopLy+HEAJDhw71/+YQMEscTgCuUAdCRERN5efnY+XKlUhISEBMTIzd4QAwb6rax1FVRETh4brr\nrgMAnH/++TZHcgZHVRERRQCXK3wagjiqiogoTEkpPdPhcOGfLqhRVUKIsQCOQbtIUAgxVEr5XQji\nJCJqc4y1jH//+982RtJQsKOq8gAMB5AlhLgd7o50IiKygD4Md/78+WFx/YauSY1DSrlOCDFWCJEg\npfyiUVkRgDlCiO4AFLibtYiIyAL6hX/hcGNDI293x3XB3c9xvz6jUS2kUkq5Tkr5rdUBEhG1RQcP\nHsS4ceMAhF/iMOvjAIDuUsrxxhlSyu+EEKkAtkopq6wPjYio7XrjjTc800lJSTZG0pS3Godp972U\nch2ANOvCISIiAOjQoYNnOj4+3sZImvKWOOJ5/QYRkX30/o377rsvbK4Y13lLHNkAFgoh+pmUcSQV\nEZHFSktLcfHFF2PWrFl2h9KEaR+HlLJSCPE43MlDwj0M9yiAEQCWhzA+IqI26ciRIxg8eLDdYZjy\n+sxxKWWRlPJaAI8DKARQCeBxKeVHoQqOiKgt+vjjj1FUVBQWT/sz06TGIYQYAPeoqu8AQBtyy2G3\nREQhUFdXh8zMTADh8ZhYM01qHNpFfkIIMVkIcb/x+g0iIrLWkSNHPNPduoXnGCVvfRyeWoYQYph2\nvyoJ9z2qvjBbhoiIWq6kpMQzHa59HN4uAPRolEQGGJJIBYDcllwMqChKOtxXqaeoqvpiIOVe5k3R\nFhmoqmpWc+MhIrKbnjiWLVuGIUOG2ByNOa+d42a0DvOFUso3AbToYkBFUVIAQFXVXAAu/bWvci/z\n0gDkqqqaDcCpvSYiikh6U1W4Jg3AJHFo/Rr9/S2o3a+qJSOsJuDMI2oL0TQJmZWbzXMali0ErzMh\noghWUlKCjh07onv37naH4pVZU1VhiB4d60DDu+v2CKC8ybxGTVwp4HUmRBShDhw4gI8//hhnn302\nhBB2h+OVWeIYJYTwJA8hRLdIuKmh1oyVp6pqnnF+WVkZFEXxutzkyZORkZFhdXhERH7NnDkTpaWl\nYXdTw8bsfHSsC2dupuiA+8r0QMq9LZNm1jGelJSEnJycVgmYiMhK27ZtAwB07NjR5kh8M+scd0kp\nH5dSjtdurZ5r0baX40x/hFPfjqIoDh/l3paZYhhhxc5xIopIcXFxAIDFixfbHIlvwT46ttXoTUra\ngd5laGJa563cbJ42PVtRlH2KohyzIlYiIqudOnUKJSUlmDFjBn73u9/ZHY5PZk1VTgD7QrFxbQht\n43nD/ZRnN3qdCyC8blZPRBSkkpISnD59GsnJyXaH4pdZjWNfiEZVERGR5vDhwwAQsYljlPE6Dj7Q\niYjIWgcPHsStt94KAOjfv7+9wQTAzlFVREQEYOXKlZ7pvn372hhJYMwSh0tKOUd/od2bSp8eC+AY\ntIsEhRBD9duvExFR85SWlnqmjc8aD1fBjqrKAzAcQJYQ4nbw9h5ERC1WUlKCgQMHori42O5QAtKk\nxiGlXCeEGCuESGh8C3XtWR1zhBDdAShoePsPIiJqhp9//jlsH9pkxtvdcV1w93Pcr89oVAuplFKu\n0265TkRELVBUVIQBAwbYHUbAvD2Po7t21biHlPI7IUQqgK2RcO8qIqJIkJGRgYqKCjidkdPy763G\nkWA2U0rZomdwEBHRGSdOnMCaNWsAAH369LE5msB5SxzxvH6DiMha+kV/AHDZZZfZGElwvCWObAAL\nhRD9TMoipz5FRBTGDh06BAD473//i4QE04aesGTaxyGlrBRCPA538pBwD8M9CmAE+KAkIqJWoSeO\nSGqmArx3jutDb68VQgyDe+gtADyuzSciohY6dOgQhBDo1auX3aEExWvi0GlDbr8FACFEfyHEG1LK\naZZHRkQUxTZs2IB58+YhJiYGsbGxdocTFG99HA0IIboJIf4N94OTxvt7PxER+TZr1iwAQH19vc2R\nBM9n4jAkjG0AVCnluQDYVEVE1EK//vorAODDDz+0OZLgeU0cQoh/wJ0wPpdSnielfFMrkiGJjIgo\nSkkpcfjwYUydOhWXXHKJ3eEEzVcfx1oAW2HdM8eJiNqk8vJyHD9+HL1797Y7lGbxNapqHQAIIVK1\nmxrm8lYjREQtI6XE9ddfDyAynr1hJpBRVcYE4gAgLI+KiChKbdy4EUeOHAEADBs2zOZomsdv4tAZ\nEoh10RARRbkdO3YAAGbMmIEePXrYHE3zBDQc10hKGXlDAIiIwsT+/fuRmJiIrKwsu0NptqATBxER\nNd/evXsj6tkbZpg4iIhCZOzYsVBVFb///e/tDqVFmDiIiELgxIkT+OabbwBEbqe4jomDiCgECgoK\nAADjxo3DTTfdZHM0LcPEQUQUAt9//z0AYNq0aYiJibE5mpZh4iAistjBgwcxY8YMAIioZ4t7w8RB\nRGSx119/HQAQGxuLDh062BxNyzFxEBFZbMuWLejSpQvy8vLsDqVVMHEQEVmotLQUe/bswYMPPoj4\n+Hi7w2kVTBxERBapra31DL29/PLLbY6m9TBxEBFZ5PPPPwcA9OrVC0OGDLE5mtbDxEFEZJGcnBwk\nJSVhy5YtUXWDWCYOIiILvP/++1izZg0GDx6Mdu2i61AbXZ+GiCgM1NfX4/nnnwcA/O1vf7M5mtbH\nxEFE1Mruu+8+HDt2DNnZ2fjtb39rdzitLuAHOVlBUZR0AC4AKaqqvhhIeaDziIjssG7dOqxduxYA\ncO2119ocjTVsq3EoipICAKqq5gJw6a99lQc6L1Sf4amnnkLv3r0hpQzVJokojJ0+fRr33HMPAGDr\n1q1RcZW4GTubqibAXUsAgEIAaQGUBzovJBYtWgQAcLlcft5JRG1BZmYmACArKwvJyck2R2MdOxOH\nA0CF4XXjh++alQc6z1L5+fn44YcfPK+Li4ut3mSb8e6776J37944ceKE3aEQBWXLli1YuXIlLr/8\nckyfPt3ucCxlax9HKJSVlUFRFK/lkydPRkZGRsDrKy8vx3XXXddg3saNG7FgwQL885//jNqqaajo\nz2Heu3cvLrroIpujIQpMRUUFbrvtNgDAv/71r6i6ZsOMnYnDBSBBm3YAOBpgeaDzAABJSUnIyclp\npZCBbdu2NZn37LPPAgA++ugjLFu2DFdddVWrba+t6tixo90hEAXs5ZdfBgBMmjQJiYmJNkdjPTub\nqpYD0G9M7wSQCwCKojh8lAc6zzIlJSUNXvfs2bPB60WLFuGFF15o0JQVjU6ePImnn34an332GaZP\nn46ysjKcOnWqxeuNjY0F4B4H/+2336K+vr7F6ySy0tSpU7Fo0SKMGTPGcxIZ7WxLHKqq5gGAoihp\nAFz6awDrvJUHOs/KuA8ePNjg9YYNGzB27Fj8/e9/R8eOHbFx40bMnz8fqamprbK9uro65OTk4OTJ\nk62yvtaSmpqKN998E/feey8++ugjDB06FP369cPx48c97/nkk09QVFQU1Hr1RJGamoobb7yx1Wtv\nx44dw6233oply5bhk08+adV1U9uze/durF69GgDwf//3fzZHEzoi2oeSTpw4Ub7zzjutsq49e/bg\nmmuuAWDeIf7ss8/ijTfe8Lxevnw5PvzwQzz99NNwOBxN3u/NqVOn8N5772H16tX4+uuvPfN/+OEH\nFBQUID09HZs2bWpS2wHcQwAvvPBCtG/fHhkZGVAUBRkZGXjmmWdw1lln4f7770diYiLKy8vRrl07\nJCQkNFmHUXV1NYqKihrcoG3Dhg244447TN8/bNgwVFdXe56vHBsba5o8jh8/jlOnTqFz584AgNWr\nV2Pq1Kmm67zxxhsxYcIEHDp0yDPUMVilpaU4evQo0tLMB93NnTsXN910E7p06dKs9UspUVNTg65d\nuzZreYo8xt/Brl27gvqNRwJFUbapqmraQczEEYQ1a9Z4OtLNEseSJUswc+ZMz+tu3bqhqqoKDz/8\nMB566KGAtlFQUIA//OEPAb330ksvRfv27XHw4EGkpKQgKysLF198cUDL6t5++22ftaP09HR8/fXX\neOSRRzB37tyg1q0z7qsTJ05g1apVePDBB5u1rh9++KHBwXnJkiWorKxESUkJnnjiCXTr1g07duxA\n3759sW3bNlxwwQWYPHkyduzYEdD633zzTaSmpnqazHypqqpCZWUlLrnkkiZlf/3rXzFixAiMHDky\n8A9HEePdd9/1DOTIzMzEE088YXNErc9X4oj6UVWt6ZdffvFZrh/wb775ZqxatQpVVVUAgM2bNwNw\nn5X6Gm0xe/ZsvPrqqwHHY6yN/PTTT/joo48CXlann8GvXbsWP/74I/bs2YNXXnkFAPCPf/zDsw2z\npDFz5kyUlJRg4cKFPrfx8ccfIzMzE/Hx8UhMTMSPP/4YVIyjR4/Gp59+CgANbt/wl7/8xfNITgDY\nt28fNm3a5Hd9/fv3R3x8PL799tsmZffffz8AYN68eRg/fnyT8tLSUiQmJmL79u248cYbvW5Dv08R\n4N53o0ePbtEZaUVFBeLi4hAXF9ek7MCBAzh9+jT69+/f7PVT4KSUnqRx++23R2XS8Ic1jiAsXLgQ\nf//73zF37lyvTTW//vorYmNj0a9fP5w+fdr0Pfv370f79mdydl1dHTZv3oy77rrLM2/UqFHo3bs3\nHnnkEcTGxuL8888HADgcDtxzzz0BJZhx48Zh9+7d2LVrVzAfMyBz5szBnXfeCcB9tez+/ftxxRVX\ntMq6U1NTMWDAAEybNq1Bc9wNN9yA7777rlnrzMvLQ/fu3VFQUGA6zPfQoUNB19aC9dBDD+HOO+/E\n2WefjZiYGK/vq66uRufOnXHq1CkUFBRg//79noQGAHfffTe2bt2Kjh074tNPP0Xv3r0B8HqiUJg4\ncSK++OILz2vj7yDasKmqlRLHk08+icWLF+PgwYN+b5Os/5jvuusuvPfeew3KrrjiCixduhSxsbH4\n5ZdfmtwEbd++fejUqZPP9RcWFmLBggW4+eabcdlll6GqqgozZ85E9+7dsWTJEmRnZ+OGG24AACxd\nuhRXXnklBg4ciJMnT2L79u248MILoaqq1wRo1LNnTxw5cgQ9evTA0aPu0c5mByn9MzfHSy+9hE2b\nNqF79+4+R6ZIKbF27VpMmjQp4HUXFRUF1PQEADU1NbjgggsCXjfg/o5vu+029O3bF2eddRbq6+tR\nXV2NqVOnIj8/33SZXbt2IS4uDrW1tYiPj4eqqqisrITD4cAtt9wS1PZ1W7ZsadF3QL7l5+c3uIYr\nJSUFq1atitprNpg4WilxBHNmN3z4cBw5cgRffvklpk2bFvDw3IKCAtPmCCsdPXoUQ4YMQWZmJjp3\n7oz09HSMGDHCU37gwAHExMSguLgYI0eORFpaGt56660m6zEetHr16oWff/7ZdHurV69GQUEBxo0b\nh9ra2mZ1SM+bNy+gPpfrrrsOixcvDnr9xcXFuOWWW/Dzzz9j0KBB+OMf/4jBgwcjLy8PTqcTV111\nFQYMGOBzHbW1tcjNzcUDDzwQ9Pab65tvvsEHH3yA0aNHIyYmxlNT9ebHH39EQkICTpw4gcTERF7A\n6sXKlSubfI9bt26N6tuKMHHYkDh27NiB/Pz8JtXYxiOvdKmpqViyZInPJoxwoKoqhgwZYnoGP3/+\nfCxatAilpaUYNGgQbr75ZsyePRsAMH36dGRmZjZ71JKZSZMmeR7NabR+/XpPf9POnTv9jhwLhV9+\n+QXnnXde0Mv16dMHTz31FNavX49OnTp57o8WqO3bt+PPf/4z3njjDfzmN79pULZ3715cffXVDeal\np6dj7ty5TCCaQ4cOYezYsTh06JBpWbTWNgDfiQNSyqj+d/fdd8vW4HK5ZHJyskxOTm7Reurq6mRy\ncrKcP3++zM/Pl1u2bJHHjx9vlRjDwc6dO2VycrL88ssvLd/Wk08+6flO3n//ffnFF1/I+vp6y7fb\nXK+99ppMTk6Wu3btkqNGjfLEbvx39913y8rKSq/rOHz4sFy/fr2UUspjx46ZrsPbvxUrVsjk5GQ5\nb948v++tqakJ1W4Ja7NmzWqyb7Zv3253WCExfPhwVXo5rrLGEaADBw7g0ksvBcBOSH9OnjwZkjPW\nzz//HJMmTcILL7zQ7Os7Qun06dM4cOBAq45+OnnypKWjqe69915kZmbinHPOsWwb4er48eNwOp1N\n5utNt9HOV42DTwAMUE1NDQA0ucEhNRWqZo5rr70Wq1atwsSJE0OyvZZq165dqx/kO3TogOLiYhQU\nFKBbt26tum4AWLx4MVJSUtC7d29s3769TZ00mSWN4uLiNpE0/GHiCNBjjz0GABgzZozNkZDR8OHD\no7qdOVBxcXHYvXt3UMskJSVh06ZNDYY4z5kzB8XFxdi/f3+T948ZMwYjR47E9ddfj82bN+PXX3/1\nlBlviVNTU4Oamhrs2LEDhYWFKC4u9tz1IBIcOXLEdHRaW7xewxteABgg/WIx44+FKJz4S6CfffYZ\nKioq8Pzzz+PBBx/E6NGjPWXr169HTU0Nhg4dCgBo3769p3bx2muvYcmSJZ4bfO7cuRO33357s+Mc\nM2YMHn30Ub8jvuywbt0602bP++67L6jHL0Q79nEESD8D2bt3r+f+SkTh5tixY16fYxLI9Ue+SCkx\nc+ZMLF26tNnrMNP4glg77NmzB2+99ZbpMPN169Y1udaqLeAtR1pBTEwM6uvrmTQorMXHx+Occ85B\nSUkJNm7ciPj4eGzevBnLli1rUdIA3DWa5557Ds899xwAd+27tLQU+fn5SEhIwMmTJ9G9e3d06tQJ\n5513HoQQqKurQ1xcHE6fPo3TmNe7AAAIi0lEQVQJEybgf//7X5P19uvXD4C71uN0OlsU59dff42R\nI0d6+iFOnz6N3NxcpKWlma738OHDDa5ZMoqNjUVBQQH7NEywxhEAKSX69OkDgCOqiFpDXV0dXnjh\nBa/XpSxbtgznn3++55YzUrvPW1FREW655RY8+uijePvtt/H888/jwIEDmDFjBgYPHoydO3d63eag\nQYNQVFSEvn37Yu/evbjwwgu99gulp6d77tnWVvECwBYmjrKyMgwdOhS9evWCqqqtFBkR1dXV4cUX\nX0R2drbdoXiEQ9NZOOBw3BZ65plnALibAYio9cTFxeHpp59GcXExDhw4gGHDhtkWy5w5c/DVV18x\naQSAeygAx44dAwDT5y4QUeuIiYnBmjVrALjv2fb222/7vcWKEAJSSlx00UU4dOgQXC5XUNt0OBx4\n4IEHkJycjNtuu63Zsbc1TBwB0G/Wl5mZaXMkRG3Dueeei1mzZmHWrFkA3M8jqa2tRW1tLbp27er1\nLsAvvfQSPvjgAyxduhROpxOxsbF47733sGvXLlx11VVo166d587KF1xwQYNbpFPgmDj82LFjB/r0\n6YPvv/++Td52gSgcJCQkBHSzyocffhgPP/xwg3nG59wA7gEupaWl6NGjR6vG2JYwcfihXyQ1aNAg\nmyMhotZy9tln2x1CRGPneID69u1rdwhERGGBicMH46Nf9es4iIjaOiYOH4w3buON9IiI3Jg4fKiv\nr/dMc/QFEZEbE4cPxhqH/ghUIqK2jonDB2ONo/HzmomI2iomDh+MNQ6Hw2FjJERE4YOJw4dTp055\nprt06WJjJERE4YOJwwdj4uCoKiIiNyYOH4yJg4iI3Jg4fGDiICJqionDBz1x3HTTTTZHQkQUPpg4\nfNATx+23325zJERE4YOJwwc9cXTo0MHmSIiIwgcThw964oiJibE5EiKi8MHE4cPLL78MgDUOIiIj\nJg4vSkpKsGHDBgCscRARGTFxeFFeXu6Z5sV/RERn2ProWEVR0gG4AKSoqvpiIOVe5k3RFhmoqmpW\na8TWqVMnzzQTBxHRGbbVOBRFSQEAVVVzAbj0177KvcxLA5Crqmo2AKf2usWMiaNdO1bMiIh0dh4R\nJ8BdcwCAQgCND/hm5WbznIZlC7XXLdaxY0fPNGscRERn2NlU5QBQYXjdI4DyJvMaNXGlAFhuXElZ\nWRkURfEaxOTJk5GRkdFkvnEkFUdVERGdYWsfR2vSmrHyVFXNM85PSkpCTk5O0Ovr1q2bZzo5ObnF\n8RERRQtLE4eh09qoUO+jAJCgzXMAONrofd7KvS2T1lod40DD5qmEhAQf7yQialssTRxah7U3ywHo\nbUhOALkAoCiKQ1VVl7dyL8tMMYywStMSExERWcC2znG9SUkbBeUyNDGt81ZuNk+bnq0oyj5FUY6F\n/IMQEbUxQkppdwyWmjhxonznnXeatWzv3r0BAMXFxa0ZEhFR2FMUZZuqqqYji3iBAhERBSVqRlVZ\nIScnB3FxcXaHQUQUVpg4fBg8eLDdIRARhR02VfmwYMECu0MwFa5xAeEbG+MKDuMKTluLi4nDh4UL\nF9odgqlwjQsI39gYV3AYV3DaWlxMHEREFBQmDiIiCgoTBxERBYWJg4iIgsLEQUREQWHiICKioET9\nvaoURSkDsL+ZiycDONyK4bSWcI0LCN/YGFdwGFdwojGufqqqJpkVRH3ioNalKEqK8WFZiqKkw/3s\nlBTDre0DmmdxXLNVVc3SbrmfbVdcRFZTFOWxYH97Lf27Z1OVF4qipCuKkqYoymM2xjBb+3+KYV6T\nuEIVq3YL+w8Mr1MAQH8wl6IoKYHOszIuzRRFUfbB/Rz6gGNt5bimaP9mG+YF9P1Z+Z16icv2vzVt\nvWlhuL/M4rJ9f+mxARilTYfs98jEYcLqA0oQbD/4GWnbKDTMmgD3WQu0+WlBzLMyLgCYrKrqQMND\nvUIal/aDztVqO07t4GF7ojWLSyuy9W9Ni2Octu6UYPZNqOPSisLqt6kJ2e+RicOcpQe6INh68AuA\nA0CF4XWPIOZZTT9Y62d6oY7LiTPfRaH2OhwSrVlcgM1/a6qq5qqqmqHHqDU72r6/vMQFhMFvU3E3\nzxqfdhqy3yPvjmvOjgOdGf2MUG+HDJeDctgztOWOMpxVh3L7xscmp8D9KOThsDnReokLCJO/NS3R\n6wfqsDkxaRQXEB77K6GV1xcw1jjCmKqqL2pnFD3sOPgFwIUzf7wOAEeDmGcZrf0+XXt5FO6zalvi\n0pon8gxnqmGhcVzh8remHYQzFEVx2BWDmcZx2b2/TGobQAh/j6xxmAvpgc6M1ulWoarqCvg++MHL\nvFBYDkB/tKQTgP6HHOg8q6g40+cxEMACbZ4dcaWpqpqlTQfz/Vn9nXriCoe/NUN/QB7c392UIGMI\nWVyKorhg/2/TqSiKU1t/ghZnyH6PrHGYW44zbb+hONCZUQ3bHai9NosrZLFqZ/GKfjavn61qZ1wu\nVVXzAp0XgrjGa6/32RjXFEOTWRoC//4s/U5N4gqHv7U0NDzIFgYRQ6jjsn1/qaq6Qktcelwh/T3y\nOg4vtLOwQrg7xLL9vd/CGCq0GF40zGsQVzjESg0pZ4YIV8B94BmnqmpuoN+fVd+pn7hs+1vTmoDG\nay+H6x3SYbC/fMXVZn+bTBxERBQUNlUREVFQmDiIiCgoTBxERBQUJg4iIgoKEweRjYQQaUKIdP/v\nJAofTBxENhFCOOC+7YfL33uJwgkTB5F9xgNY4fddRGGGiYPIPsMBOKWUdtyZgKjZmDiI7ONE0+eI\nEIU9Jg4iGwghnADypJRMHBRxmDiI7FEB4KgQwqklEaKIwXtVEdlICOGQUnJUFUUUJg4iIgoKm6qI\niCgoTBxERBQUJg4iIgoKEwcREQWFiYOIiILCxEFEREFh4iAioqAwcRARUVD+H3Q5HrTOgxn1AAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc3b1494690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot( (ttCLASS[:4001]-ttCAMB[:4001]) / ttCAMB[:4001], '-'  )\n",
    "plt.xlabel(r'$\\ell$')\n",
    "plt.ylabel(r'$\\Delta C_{\\ell}^{TT} / C_{\\ell}^{TT}$')\n",
    "plt.title('TT Comparison of CAMB and CLASS');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# P(k) Comparison\n",
    "Next, we compare the linear P(k) and HALOFIT P(k). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars.set_matter_power(redshifts=[0.], kmax=2.0)\n",
    "pars.NonLinear = model.NonLinear_none\n",
    "results = camb.get_results(pars)\n",
    "kh, z, pk = results.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)\n",
    "s8 = np.array(results.get_sigma8())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "PCLASS  = np.array([cosmo.pk_lin(ki*cosmo.h(), 0) for ki in kh]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,u'LINEAR')"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAECCAYAAADelD2uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd4FFXDBfCzk00hpDdCUwi8wIuI\nECYIKjUREZAaKQooJVQBUT8QEAElIAhKiUDoSJciIJ0gEIqCAxZUFDGKUhLSSSFly/dHNrwRgbTN\n3p3d83uePGS2zJxck+Ps3dkZjdFoBBERqYckOgAREZUOi5uISGVY3EREKsPiJiJSGRY3EZHKsLiJ\niFSGxU1WT6PRBGk0mm0ajea8RqMJf8h9YablaI1G83uRx3iZbos2ff+P9Zm+Jmg0mtQHbH/bQzId\nMT1/mGkdweYfAaJ/0vA4blIDjUYzDEAzo9E4vLj7NBpNGIBgAHWK3BYEAEajMe5B6zM9L67wMabb\nwgFsA+BtNBrT7tluOIBni2zDC8AfRqPR23w/OdG/cY+bbJLRaJwLQDaVcbE0Go2X0WiMAeB1z10+\nALYDGFaC1fjc5/lEZsfiJlv2IoDo4h5k2hsv3CO/UOR2LwAppnX8a0/fxEej0QSb9uDnAGhW3tBE\nxWFxk80yTXlEazSaB5V3kGm640H39zYajdtNe+I+D5i/TjEajReMRuNyAEcApN3nMURmxeIm1SrJ\nNEjhlAmA+z02zmg0bkfBnnmaaZ1BRe6vU/jmJYAYPHivu6iJJXgMUbloRQcgKodgFBRqcV4EcB4P\nmMYwvelYuKccDCDOtHe9tXDqRKPRxAD4Aw8v7xSYplxMc+bc+6YKwT1usnqmveBnUfBmY+Ghd0cA\n+N5zX5hpeaJGo5lQ+HzTlMnEh6yvcJ3n8b/SXgFTCZsUFnK06VBAryLrKJwf3w4gzfRKoOhzicyK\nhwMSEakM97iJiFSGxU1EpDIsbiIilWFxExGpDIubiEhlLHIcd2hoqLFatWqW2BQRkU24dOlSkqIo\n/ve7zyLFXa1aNaxfv94SmyIisgmyLF990H2cKiEiUhkWNxGRyrC4iYhUxiJz3BkZGYiMjLy73KpV\nK7Ru3doSmyYisjkWKW53d3dMmTLFEpsiIrJ5nCohIlIZFjcRkcrwQgpExbgWn4SjZ3/Gt5ev47eb\nt5GcqUO2ToM8gxZ6SJBggAQjHCU9vJyBQC8n1K/hgx7tmiC4IU/LTebH4ia6R0JSGlbsPIGTP97A\n1dsScp18oNFIAJwAvQe0ukw4aXRwd8yHowOgNxR85eol3Mhxwd/JblBSjNj4w7eQco+hZuVc9GlT\nH4N7tIGTI//kqPz4W0QE4K8bifhoYwxif0lBquQHjaQFdJ7wcEhDU//baNGwBp56oi6aNqgFrdbh\noevKzM7B0bM/4sCZS/j2ag6u5nhh7tE0zDuwHo39dfhgVGfUq8VTQFDZsbjJbul0ekRtOYKtJ39H\ngt4XcHCCg94Fj3mmoVfrBuj3fEu4ODuVer1uri7o1k5Gt3YyACA1PROLNsdgj5KNb9MC0HHmCTzm\neRsfvtYFDWpXN/ePRXaAxU1251p8EqYt34fYOB30zl5Avjvqe6QhokswurfvAUky73v23p5umDai\nO6YBOPHNT5i6+gR+zPBHp8hYtK2ZhyVv9yvT/yDIfrG4yW5c/vMG/i9qLy6mugFaN1TW3EJ4Uwe8\nObA33FxdLJKhTchjiA15DLHKz3hj2TEcvxGIpsNXY/aAZugeGmKRDKR+LG6yeRcvX8XEpQfwS4Yn\njJI3qjnewoQ+T6Jbu17CMrWWG+Lc8gaYtWov1pxxwfgtf+Dz2J+wauqAYufQiXgcN9msuL8T0Gn8\nMnSb+xUuZfqgplMq1kY0xOllo+/OP4skSRLeieiKE7O7oKqUiJPx7mgxYhn+upEoOhpZORY32ZyM\nzGwMnrEOYdOO4FKWLx51TsGm0U1wYskotAl5THS8f6kR6IdTS0fi+do5SNb4of3kPThw8lvRsciK\nsbjJZhgMBkxf+jmCX9uM4zfc4KVJx5IBdXHsk1Fo2aSe6HgPJUkSPnn7Zczp9QiMGgeMXvszln12\nVHQsslI8OyDZhC/P/ojx0bHIcKwCJ+jwejsPjO4rbg67rHp3bIm6j1RBv7kxmHMkGVfjt2P22HDR\nscjKaIxGY4VvZMCAAUZeuowqQmp6JobO2owLKR7Q6PPRub6E+W/0Vv0nFK/FJ6Hz5M+Q4VgF7apn\nYtW7r4iORBYmy/J5RVHu+2YMp0pItRZtOoTm47bi2zQf1NAmYf/Udlg88SXVlzZQMO99OmowAozx\nOHbdDQPeXSM6ElkRFjepzrX4JLQZtQQLTmRCAwMmdfDGyaWjbe5TiG6uLji+OALVNPE4neCBvpNX\nwWAwiI5FVoDFTaqyYMNBtJm0F3/l+6GJZwrOR72MiF7tRceqMC7OTvhy8XDUdEjAuWQvhL+9kuVN\nLG5Sh8K97EUns6A15mFe70exc24E3N1cRUercE6OWhxdNBy1nRLwXbovur61guVt51jcZPU+3XMS\nbSftwV/5fmjqlQIlqj96Pfuk6FgWpdU64MjCEajnkoifs/zQ+Y3lLG87xuImq5WTm4fek1Zi2t54\naIxGzO5ZAzvm2Mde9v1IkoT9Hw9DA9dE/HrHH8+Pj2Z52ykWN1mlsz/8hmYjVkFJ8UZNh0TEftgT\nfZ9/SnQs4SRJwt75w9DILQm/5QTguXHR0On0omORhbG4yeq8E7UD/Raewx3JAwOaSIhdOgpV/b1F\nx7IakiRh14cReMIzBb/nBaDD6yxve8PiJquRkp6JNqOWYNNFwNVwG1tefxIzRvYQHcsqSZKEHR8M\nQbB3Cv7Mr4K2ry1DZnaO6FhkISxusgonvvkZT43biL90/gj2ToGybAiaP15XdCyrJkkStn8QgWeq\n3MYNYyCeeW01EpLSRMciC2Bxk3CRK/Zg0LJvkSe5YnxbD2z/IIJXhCmFT98bhBcbGpDu4Ie2b23F\nT1f+Fh2JKhiLm4TJzslF5zeisUrJh4shA9vfehpjX3pOdCxVmjPuRYxv54lcB3d0jTyKDXtPiY5E\nFYjFTUL88OufaD5yNS5l+aFepUSci3oFwQ2DRMdStbEvPYfFAxtAMuoxdfd1jP5gIw8XtFEsbrK4\nDXtPoccHx5EteWJAEwkHF4yAW+VKomPZhM5tgnHsg27wNSbhwB8ueGbUUlxPSBYdi8yMxU0WNX7e\nFkzdfR2SUYdlQxrxqJEKUCPQD18vG4GnAtJxU++PNhN3YdO+06JjkRmZtbhlWY425/rIdmTdyUGH\nsUux+zdHeBmScXRWN3R4qrHoWDZLq3XAhvcHY1aP6tAAmLLrOvpNWYW8fJ3oaGQGJSpuWZaD71kO\nl2U5TJblCQ96DFGhy3/eQItRq3ElNwCPuyfj7LJheKSav+hYdqFfp6dxfE53VHdIxNkkLzQbtgIX\nfo4THYvKqdjilmU5DMC2IsvBAKAoSgyANFmWg2VZDgKQZvoiumv3l9+g04yDyHLwwcuNNdg9b5hN\nXOhATapX8cWpZaPRtxGQJXkifN4ZvBe9W3QsKodii9tU0EX/F90H/yvoOABhAIIAeAEIMpU4EWav\n+gLjN10BAHzUrw7eH91TcCL7NmtML6x/rRlcjJlYe0GHdqOXIDU9U3QsKoOyzHF7AUgpsuxbpNy9\nzJKKVM1gMOCVaWux/GwOKunTsPfd59A9NER0LALwdNP6UJYOQhPPFPyZ74+24z/lpy1VyGxvTiqK\nkqYoyrOKonACzY5lZucgdMwynIx3R1XNLZxaMMDmLimmdq4uztg5NwL9Hpdw28EPof+3GSnc81aV\nskw2pgHwMX3vBaDYg0QTExMhy/e9WDEAICIiAsOHDy9DFLImf16/ha5TtyPTsQqaeadi88wR0God\nRMeiB5g1phcqLd+N1YofOry5Dqc+4akG1KIsxb0VQGELBwGIKe4J/v7+OHjwYBk2RWoRq/yMIVGn\noXP0Rb9GwKwxQ0VHohKYOqwbElI3Yn9cIPpMWYvd84aJjkQlUJKjSsIL/pHDAUBRlAum28MApBUu\nk/1atfMYXl1yHnrJGe91rYZZY3qJjkSlEDXxZdRxuoWLGb6Yv26/6DhUAhqj0VjhG+nevbsxJOR/\nb061atUKrVu3rvDtUsWbvHg7Nv+gh2P+bax7vQ1aNqknOhKVQXpGFlqM2YBcqRIOvhuGerWqiY5k\n92RZPq8oyn3nmC1yQK27uzumTJliiU2RhRgMBgycthZnbnnCQ5+MfR/0Ro1AP9GxqIw83Stj0fCW\nGLb6El6d/TnORI8WHYkegucqoVLLyc3Dc+OiceaWJ6pJCTgdNZilbQOebdkYLatkIR6BmLeWUybW\njMVNpZKcloFnRq/A73kBeMIzBSeiRsDN1UV0LDKTle+8DG1uCqKP30RGZrboOPQAFpkqycjIQGRk\n5N1lznGr05WrN9Ft+i5kOwagU+0cfPJ2hOhIZGauLs4Y+9yj+Oh4BsbM+wxrp78qOhLdB+e4qURO\nX/gFry46Bb3WGyNaumLioBdFR6IK8lq/Dlh79BPE/u2FK1dvou6jVUVHontwqoSK9dnBrzAw6hz0\nkhNm9qyJiYM6i45EFWzO0NYwOjhh7AKejMoasbjpoeav24+3t/8JyZCPNSOb4aVOT4uORBYQ2uJx\nBDkl4VKmN7767rLoOHQPFjc90Ph5WxB1Ogsu+nR8Ma0j2oQ8JjoSWdCCsZ2gMRrwf8uOiI5C9+Cb\nk/QvRY/R9jYk4ND8l+Hv4yk6FlnY4/UeRROfw/guPRC7jn7DMzxaEb45Sf+g0+nR9f9W4Jdsf1ST\nEnB44RC4ujiLjkWCLH6zJ1pN2ofpm1jc1oRTJXRXdk4u2o1Zhl+y/dHANQnHFw9nadu56lV8ERoE\n3NZWwZKtxZ5PjiyExU0AgJT0TLR6bSWuGwLxVEA69s6P4ClZCQDw0evhkHLTsPjA79Dp9KLjEFjc\nBODq9Vto/fp6pEgBeKFuHja8PxiSxF8NKuDu5oq+shdynf3w3nIeHmgNeHZAO/fDr1fx4uxDyHP0\nxJDmlfBORFfRkcgK6XR6NBq6EnpocSGqH9zdXEVHsnk8OyDdV6zyMwZHfQWDgxve7uCH4S+Gio5E\nVkqrdcCYjkGYf+w2/m/hDiybMkB0JLvG18N2atfRbzBoiQKD5Ih5/eqytKlYI3uHwiM/AYd/N+LK\n1Zui49g1FrcdWrHjS7yx6TI0Rh3WjGyGnmHNRUciFZAkCR8MfgpGB2e8NHMnDAaD6Eh2i8VtZ+au\n2YdZB27BUZ+NnW+356chqVQ6PtME7WrkIEkKxBsffSY6jt1icduRKVE7sPRMFlz1aTgU2Q2N69cS\nHYlUaNnkl+GpS8DuXyVEbzsqOo5d4kfe7cSYORux93cneOiTcGTeywjw5UfYqWycHLXYO6sP2k/c\ngTkH8+Hl/hX6dGwpOpZd4VElduCVaWtxMt4dvoZ4HF34KjzdK4uORCpXvYovNk8IQ5+5X2LStjjk\n5ukwsGsr0bHsBqdKbJjBYEDPCStwMt4dVTUJiI0aytImswluGIQdk8Kg1Wdj+p5rnDaxIBa3jdLp\n9Og0fjm+S/dBbcdbOLZ4GM87QmbXuH4t7JnWCY66DMw5lIgFGw6JjmQXWNw2KCc3D+3HLMPlHH88\n5paMI4uGw8nRIrNiZIca1K6OA+93g7MuDQtPpGHeOl4hvqKxuG1MRmY22oxZgWuGQLTwS8fuD4fy\nvCNU4YJqVsHh2b3gkp+KJbGpOHfxiuhINo1/0TYkOS0DrcetQaImEM8+ko1NkTxZFFlOzar+WPtG\nwSdwB390BDm5eYIT2S7+VduI6wnJaDN+PdIcAtCrgR7RPJcECfBk4/+gbxMXZDsFYNLinaLj2Cwe\nx20Drly9iRem7UaOox8GNXPEu8O7i45Edmzm6B7YMyQae36qhCmpt+Hn7SE6ks3hcdwq98OvfyJ8\n9hHotF4Y18YTr/d/TnQksnOSJOHtXo/j3S8SMGbeNmyOHCI6ks3hVImKffXdZfSafRQ6B1dM7VKV\npU1Wo3+XZ+BriMfZm05ISr0tOo7NYXGr1NGvL6L/wjPQS06Y26cuBnVvIzoS0T+MeaEx4OiKadFf\niI5ic1jcKrTvxAUMW/4tAA2WDH4c4R2eFB2J6F/6d3kaLnmJOPzLHeTl60THsSksbpXZeeQsxqz7\nCTDqsXqUjI7PNBEdiei+JElCnycDoXf25icqzYzFrSIb9p7CW1uuQDLkYePrT/Nc2mT1Jg7qBCkv\nHetj/xIdxaawuFVi5c5jePfzv+Cgv4NtE9ujxRP1REciKpaLsxPaBmmR5RSAbYe/Fh3HZrC4VSBq\n82HM2p8AR30Wdk/tiKb/rS06ElGJvT/iBUCXg493nhcdxWbwAzhW7sO1+7HkdAacdenY91431Hkk\nUHQkolKp6u+Nx7wy8eNtf1z4OQ7BDYNER1I9fgDHir0XvRtrlFxU0qXi4KxeeKSav+hIRGUyY8iz\n6PWxgmmrjuCL+cNFx1E9TpVYqcmLt2PN+Ty46ZJxdG4fljapWnDDIFR1SMRPaW64mZgqOo7qsbit\n0Ph5W7D5ohEeukQc++hlVPX3Fh2JqNzeCg8BtC4Y99EO0VFUj8VtZUbO2oDdvznCx3ALxz8eyBP0\nkM3oGdYcgZp4KIluUH78XXQcVWNxW5FXp6/FoauV4GeMx4mFg+Ht6SY6EpFZRY19HoARoxYehMFg\nEB1HtVjcVsBgMKDv5FWIvVlwUd/ji4bCrXIl0bGIzC64YRDa1MxDkhSIkbM3io6jWixuwQwGA3pN\nXIlzyV54RJvAi/qSzVsxpT98DfE4/KczVu44JjqOKrG4BTIYDOj8xnJ8f9sXdZxvIWYhL+pLtk+r\ndcCumX3hmJ+OWfuv49Dp70VHUh0WtyA6nR4dxi3Dr3f80bByEg4tGA6t1kF0LCKLqF7FF5veag+N\nQYfRK8/jrxuJoiOpCotbgJzcPLQfswxxeVXQ1CsFe+ZF8KK+ZHfkRnUwu0996LVu6DPjM75ZWQps\nCwvLzslFu7ErcM0QiJb+6dgxh6VN9qt3x5Z4OjATCQjEjOjdouOoBhvDgjIys9FmzEokIBDtqmdh\n48zBoiMRCbf63YFwyk3CprMpyM7JFR1HFVjcFpKanok249YgWQrE87VzsOrdgaIjEVkFJ0cthrWv\nCb2zFyZHfS46jirw7IAWkJKeifbj1yFdG4Be9XWY98bLoiMRWZXX+3fEqmPR2PuTC2bl5PKQ2GLw\n7IAVLDElHaFvbUSGNgD9Hpcwa8yLoiMRWR1JkjCkfS1EncrGh2sPYNqI7qIjWTVOlVSg+MRUtHtz\nIzIc/DCwqRazxvQSHYnIao3t1wEOuWnYdvam6ChWj8VdQa4nJCN0whZkaf0wOMQZM0b2EB2JyKpp\ntQ5oU8cR2U4B2HX0G9FxrBqLuwL8dSMRz07chmytD0Y8VRlTh3UTHYlIFaZHdAJ0ufh4B4v7YVjc\nZhb3dwKem7wDdxy98doz7pg4qLPoSESqUSPQD3Uqp+PvPB9cucopkwdhcZvRb1dvoNPUXcjVeuGN\ndt5485VOoiMRqc6Efk8DDo6YvvKA6ChWi8VtJpfirqHLtL3I03pgQgd/jHmpg+hIRKr0bMvG8MhP\nwFfXNMjIzBYdxyqxuM3g4uWr6PbeAeQ7uGFKp6oY0TtUdCQiVYvo8B8YnTwwcfFO0VGsEou7nL69\nFIees45A51AZ07vVxJCebUVHIlK9kb1D4ZafgEOX9UhMSRcdx+qwuMvhm4tX0HvOMegdKmFmr0cx\nsGsr0ZGIbIIkSZgU/gSMTm4YNXeb6DhWh8VdRl9/fxn95sdCLzljTp86eKnT06IjEdmUfp2eRjVN\nPJRkD+w7cUF0HKvC4i6Dk+cvof+C0zBIjvi4fwO82KGF6EhENunTKb0g6e7grTXneObAIljcpXTs\n3I94NeosDBoHLHrlMXRrL4uORGSzgmoGYsjTPsh19kfETF5cuBCLuxQOn/keQ5eeBzQaLB3yBLq0\nCRYdicjmTR76AgIRjzMJlXHw1Hei41gFFncJ7Y+9gJErvwdgxIrhzfDc00+IjkRkNz6d3BMafS4m\nrj4jOopVYHGXwJ5jCl5b+yNg1GPN6OZo/2Qj0ZGI7ErdR6uifW0gw7EK1u2OFR1HOLMUtyzLXrIs\nh8myPEGWZS9zrNNabD/8NV5ffwmSUYf1455Ca7mh6EhEdmnu2J7Q5GXio90/io4iXImKW5bl4HuW\nwwuLuvAmAAqANABB5o0ozpYDZzBh6xU4GPKw6Y3WeKpJfdGRiOyWj6cb2pr2ujftPy06jlDFFrcs\ny2EAthVZDgYARVFiAKTJshxs+h4AvBRFsYkDLtd/cQqTt/8BB30utk5oh+aP1xUdicjuffBadyD/\nDhbtsu83KYstblMpxxW5qQ8K9qxhuj1MluVhiqKkAbhQZC9ctVZ/fhzTdv0Frf4OdkwKRXBDm3kR\nQaRq/j6eaOybjQQE4OvvL4uOI0xZ5ri9AKQUWfYFoJj2zIMAbDdHMFGitx3FzH034ajPxO6pHdG4\nfi3RkYioiPeGdgCMRsxY86XoKMKY5WLBtjI9ErX5MD46mgInXQb2TO+CerWqiY5ERPdoXL8Wamj3\n4dcMb/x1IxGPVPMXHcniylLcaQB8TN97AUgu7gmJiYmQ5Qd/wjAiIgLDhw8vQxTz+Xj9ASyKvQ1n\nXTr2v98NQTUDheYhogeb2KcFxm7+E1Oj92HdjFdFx7G4shT3VhQcRQIUTI3EPOSxAAB/f38cPHiw\nDJuyjLlr9mHpmUy45KfiQGRP1KoeIDoSET1El7bN8M6nZ3HqL1dkZufAzdVFdCSLKslRJeEF/8jh\nwP+mRUxz2mlqnyaZuXw3lp7JgqsuBUc+CGdpE6nEkLA6MDq5Y8JCVb+tVibF7nErirId97zhqCjK\n8gpLZEHTl36OdRfyUVmXjCNz+6Kqv7foSERUQqP7PouVMctw8Dc3XItPQo1AP9GRLMYsb04WJyMj\nA5GRkXeXW7VqhdatW1ti0w80efEObL5ogLs+CUfnvwx/H0+heYiodCRJwvsDmmPc5j8w6sMd2DNf\n7PtklmSR4nZ3d8eUKVMssakSmbDgM2z7WQMPXSKOzu8PP28P0ZGIqAy6tpPx0bZzuHjbF199dxkt\nm9QTHcki7O4kU69/uBnbftbAS5+I4x8PZGkTqdyicZ2hMRrwxtIjoqNYjF0V9+gPNmLPFSf4GG7h\n+IJX4O3pJjoSEZVT4/qPIsQ/GwkIxJYD9nHaV43RaKzwjXTv3t0YEhJyd1nEHPewmZ8i5u/K8DXE\n48uPB8HdzdWi2yeiipOYko4n39oDX4fb+GbFaNFxzEKW5fOKotz3AzB2Mcc9aMZanLjhDn9jPI4u\nHGJ3x3wS2Tp/H0+EVMnDNymBOHjqO3R8ponoSBXK5qdK+k9djRM33BGIeBxbNJSlTWSjZo/qAuhy\nEbnR9k/5arPFbTQa0XfyKpy55YlqUgK+XBQBVxdn0bGIqIIE1ayC/3pm4JrOH2d/+E10nAplk8dx\nGwwG9J60ChfSfFDTIQFHFg6Dk6NFflQiEihyWAf0mH8OU1fG4PCi/4iOU2Fsbo7bYDCgx4QVuJjh\nh9qOt3BowXBotQ4W2TYRidWkQW3Ucj6AK9k++OWP62hQu7roSBXCpqZKDAYDXnhzOS5m+KGu8y0c\nWjCMpU1kZ6a/0hZGyQGTlu4XHaXC2Exx6/V6PD8+Gpey/dHANREHuadNZJfahDREoOYWvk92xfWE\nYs86rUo2Udw6nR7PvR6N33IC0MgtCXvnD4Mk2cSPRkRl8HbfJwGtC97+ZI/oKBVC9e2Wl69D2Lho\nxOVVQRPPFOz6MIKlTWTnurWT4aWLx+lrDkhMSRcdx+xUfVRJXr4OoWOjcd0QiBCfNGyOHMLSJiIA\nwITwppi86ybeXLATn743SHQcs1LtUSU5uXkIHbcCN42BaOmfjo0zh5h1/USkbn2ffwof7vgEp655\n4MatFFQL8Cn+SSqhyt3T7JxctB1TUNqtq2Zg48zBoiMRkRWa2q854OiKiYt3i45iVqor7oysO2g7\nZiVuaQLRvkYW1k5/VXQkIrJS3UND4KWLx5lrEtIzskTHMRtVFfftzCy0HbcaSVIgOta6g5VTB4qO\nRERWbkTH/8Lo5Ibpy78QHcVsVFPcqemZaDtuLVIdAtGlTh6WTOovOhIRqcDQXm3hnJuIfRczkK/T\ni45jFqo4qiQlPRPtx69DujYAvRroMG98v4qISUQ2SJIk9G0RiHXf6jH/0wN4e3AX0ZHKzeqPKklM\nSUfoWxuRoQ1A30YazB7bx8zpiMjWTRzUCRu+Xo/1J/Mw4VWD6g8btur08YmpaPfmRmQ4+KF/EwfM\nHhsuOhIRqZCLsxM6NHDBHacAbNir/vN1W21xJ6XeRuiELcjS+mGQ7Iz3R/UUHYmIVOz9kd2A/Cws\n2vOD6CjlZrXF7ePphkc8NRjWohLeHd5NdBwiUjkfTzc0D8xHikMgYr66KDpOuVhtcUuShAMLRmDS\nkBdERyEiGzFrZGdAn4eZG06KjlIuVlvcRETmFlQzEPUqp+Nqni8uXr4qOk6ZqeJwQCIic5kxJBR9\no77DO8sPYfe8YaLjlInVHw5IRGROTzb+D6pKh3Ex1QvXE5JRvYqv6EilxqkSIrI7E15sDmhd8M6y\nvaKjlAmLm4jsTvfQEHjkJyD2TyMys+6IjlNqLG4isksRHf4Do5M7Zqjw5FMsbiKySyN7h8IpNwl7\nvk+HTmUnn2JxE5FdkiQJL4b4I9/ZBws2HhIdp1RY3ERktyYN6QwpLx3rjqvrmG4WNxHZLVcXZ7Sv\n44gspwBs2qeek0+xuInIrkWO6grkZ2Ph7u9ERykxfnKSiOyav48nggNycD4lAMfP/YS2zR8THalY\n/OQkEdm9WSM6o+PMWLz/aawqiptTJURk9+rVqoY6lVLxR64PLv1+TXScYrG4iYgATBvUFkaNhCnR\nB0RHKRaLm4gIwDPB/0UAbuG7pEpISEoTHeehWNxERCbjezQFHCthytI9oqM8FIubiMik7/NPwS0/\nAcfj9MjOyRUd54FY3ERERQxhS30xAAAHuklEQVRuHwSDkwdmrrDeU76yuImIihj7Ugc45iZj5/lk\nGAwG0XHui8VNRFSEJEno0dQbec6+iNpyRHSc+2JxExHd452IF6DJu41VMb+LjnJfLG4ionu4ubqg\nTS0JGY5VsO3w16Lj/AuLm4joPiJHdQV0d/DRjvOio/wLi5uI6D6q+nvjCd87iDcG4PSFX0TH+Qee\nHZCI6AFmDuuILnPOYMbaYzgc3EB0nLt4dkAiogd4rG5N1HZOxpVsH1y5ehN1H60qOhIATpUQET3U\nOwNaAw6OmLp8v+god7G4iYgeov2TjeCtj8e5G1pkZGaLjgOAxU1EVKzBYfVgdHJD5Kp9oqMAYHET\nERVrZO9QOOYmY/e3KVbxMXgWNxFRMSRJQqdGHsh19sPa3bGi47C4iYhK4t2ILkB+Fpbt/0l0FBY3\nEVFJeHu6oal/LhI1ATj97a9Cs7C4iYhKaPqQDoDRgJnrjgnNweImIiqhx+s9ihraJPya4YHrCcnC\ncrC4iYhKYXyvEEDrjGnLxR0ayOImIiqFnmHNUTnvFo7H6ZCbly8kA4ubiKiUXnqmBgxOnpj36UEh\n22dxExGV0psDn4dDbhq2nL4mZPssbiKiUnJy1KJtXUdkOQVgx5GzFt8+i5uIqAzeG/4CoMvBxzst\nf4Ucs5yPW5blIABeAMIAbFcUJc4c6yUislZV/b3RwCMDlzL98MOvf6Jx/VoW23aJ9rhlWQ6+Zzlc\nluUwWZYnmG4KBhAHIAZAuHkjEhFZp3dfbQ9oNJi+6ohFt1tsccuyHAZgW5HlYABQFCUGQJosy8GK\nomxXFCUNpj3uigpLRGRNWjxRDwG4he+SXJCSnmmx7RZb3KaCLjr10QdAmun7OBSUdWHBx3CahIjs\nychOjQBHV7y3/AuLbbMsb056AUgpsuxrKu2JAIbLssypEiKyGwO7toJzbiL2/5RpsXN1m+XNSdNe\neYw51kVEpCaSJKFHMz9s+VGDT7bEYMxLHSp8m2Up7jQAPqbvvQAUe6aVxMREyLL8wPsjIiIwfPjw\nMkQhIhJv8pAu2Dp6C9Z+GW+1xb0VQGELB6EEe9r+/v44eFDMR0OJiCqam6sLQqrqcS45ELHKz2gt\nN6zQ7ZXkqJLwgn8K5q4VRblguj0MQFrhMhGRPZs29DkYDTrM3lDxlzYrdo9bUZTtuOcQP0VRlpdm\nIxkZGYiMjLy73KpVK7Ru3bo0qyAismr/DaphOle3F24mpqKqv3eFbcssb04Wx93dHVOmTLHEpoiI\nhBnboxkm7riO91bsw9LJ/StsOzxXCRGRmbzYoQUq5d1CzOUc6HT6CtsOi5uIyIzCQ6pA7+yNT7ZW\n3BHSGqPRWGErL9S9e3djSEjI3WXOcRORrcq6k4PGo7bAU8rChVWjy7weWZbPK4py3+OoOcdNRGRG\nlSu54MlqenydFIhj535Eu+aNzL4NTpUQEZnZtKEdYTToMGfjqQpZP4ubiMjM6teujpraJFzO9MDN\nxFSzr98iUyVERPZmxqA2+P3vW6ji62n2dVtkj7vwAziFX7GxJf9kUXR0dAUmsz0cr9LheJUOx6vk\n2jVvBH3Sb5Ak89esRY4qGTBggHH9+vVleq4sy1AUxcyJbBfHq3Q4XqXD8Sqd8ozXw44q4Rw3EZHK\nsLiJiFTG4sV97/x2ccuWzGLO5zzscQ+67363F3cbx6v0t1lqzDhepcPxKjmLvzm5aNGif/xgJ0+e\n/Mdj712uSGXZVkmf87DHPei++91e3G0cr9LfZqkx43iVDser5Cz+ycnIyEh+3J2IqBw4x01EpDIW\nORxQluVEAFdNi54A0ovcXdxyNQA3Kijavdsy53Me9rgH3Xe/24u7jeNV+tuKLnO8OF7WOl6PKori\nf787LFLcFUmW5WhFUXil4WLIsuyFgmuFBgNYrihKmuBIVk2W5SAUXAw7DMB2RVHiBEeyerIsBwPw\nURSl4s5nqmJF/ga9AMSU529Q+FSJ6T920eVwWZbDZFmeUNrn2oNyjJcMQAGQhoKLPNuFcoxXMIA4\nFFwMO7yi8lmb8vw9AuiDglKySyUYu0mm/6nFABhWnm0JLW7TBYe3FVkOBgDTD5f2sGI27RGlmb7s\nQnnGq8hekJe9XOC5nOO13bRHFIZ7rrlqq8ozXiZbKzCeVSvh2HmZbksDUKc82xN6kilFUWJkWS76\nErQPgCOm7+NQ8EdzofAK80XEoGCvMQVAkCzLQfbwUrac49VbUZTlsixfkGV5gqIocy0QWajyjJei\nKGmmP8YYe/jdAso/XpbIaK1KOHZppukSAPi9PNuztrMDeqGgjAv5AnevNH+vGNMg2O1LM5RivGRZ\nVkxFFAQ72YO8j9KMVxiAiQDiZFk+8oDfQVtXmr9HoKCcfGVZtvsix/3HLhoFYwQAy8uzcmsr7lIx\n/XI8KzqHGtjL9Ii5FJmLpBKyh1dx5WF65WaWV2/C35y8RxoAH9P3XgCSBWZRA45X6XC8SofjVXYV\nOnbWVtxb8b8jHoLAPZ7icLxKh+NVOhyvsqvQsRN9VEl4wT8Fb3YUvpw3zS+m8eX9P3G8SofjVToc\nr7Kz9Nip/gM4RET2xtqmSoiIqBgsbiIilWFxExGpDIubiEhlWNxERCrD4iYiUhkWNxGRyrC4iYhU\nhsVNRKQy/w9B7fbZlKaRAwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc3b181a410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plt.loglog(kh, Pmm, '-')\n",
    "plt.loglog(kh, PCLASS/(cosmo.h())**(-3), '-')\n",
    "plt.loglog(kh, pk[0,:], '-')\n",
    "plt.title('LINEAR')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear P(k) Difference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,u'$|\\\\Delta P(k)| / P(k)$')"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAECCAYAAADzStBRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8U3W+//HX6ZJu0JYuChUplFXc\nEA8iGihluei96r0qwnVcUEek4jiAjqAy6tWro4jyQB2ZCug4V9CLoleZQR2pzoDlx6jHMgoK2IVK\nWQRa2hS60Tbn90cW0zZNkzZpkp7P8/Hog+TkJPlyKHnnuyu6riOEEMK4IoJdACGEEMElQSCEEAYn\nQSCEEAYnQSCEEAYnQSCEEAYnQSCEEAYXFewCdMXUqVP1jIyMYBdDCCHCyp49eyo0TUtvezwsgyAj\nI4M33ngj2MUQQoiwoqrqj+6OS9OQEEIYnASBEEIYnASBEEIYnASBEEIYnASBEEIYnASBEEIYnASB\nEEKEiZaWloC8rgSBEEKEuPr6epYsWcKIESP48MMP/f76EgRCCBGiVq9ezZgxYxg9ejTr1q2joaGB\nV1991e/vE9JBoKrqK8EugxBC9LRdu3bxr//6rzz++OMcP34cRVEASElJYdGiRX5/v4AuMaGq6lhN\n0wpd7s8EqoGxmqY929lzA1k2IYQINUVFRTz00EPs2LGDyMhIAPr168evfvUrCgoKyM3NxWw2+/19\nAxYEqqpOA14BhtrvjwXQNC1fVdWstiHR5rlZ2AKjOlDlE0KIULFixQpWrVpFfX2989jo0aNJS0tz\nfvjn5uYG7P0DFgT2D/xSl0OzgS3226XANKDQXktwlQ9kASeALFVVszRNK0UIIXqh7du3s2LFCnRd\nJz4+nrq6OlJSUvjtb38bkG//7vTk6qPJ2D7cHVIBNE3b6ObcfFVVk+3PEUKIXmfFihWsXbuWkydP\nYjKZGD58ONdee21Am4A6ErLLUGuaVg1MD3Y5hBDCn3RdZ+HChWzcaPsOHBsbS0NDA+np6eTm5ga0\nCagjPRkE1UCK/XYyUNmD7y2EECHh5ZdfdoZA247gYOnJINgAqPbbWdj6Arrk+PHjqKra4eNz585l\n3rx5XX15IYQIiOXLl7Ny5UrOPfdc0tPTufvuuwPeEeyNQI4ammn7Q52padpGTdMKVZtpQHVHI4a8\nkZ6ezscff+y/wgohRIB98MEHrFy5EoDU1FTWr18f5BL9LJCjhjYCG9scWx2o9xNCiFD1+uuv89vf\n/paEhASysrK45557gl2kVkK2s1gIIcLdiRMnWLBgAZ999hmJiYnU1NSQlpbWoyOCvBGWQXDy5Eme\neuop5/2JEycyadKkIJZICCF+tm7dOlauXElVVRUNDQ0ADBkyhJSUlKD3B7gTlkHQt29fli5dGuxi\nCCFEO6WlpTz88MO0tLRgMpkA2xpBDz/8cMjVBBzCMgiEECKUNDQ0sGbNGt544w0qKiqIjo5m3Lhx\nTJ8+PSgTxHwlQSCEEN3w2muv8V//9V+tNo258MILeffddwFCsimorbAMAukjEEKEAl3Xee6552hp\naSErK4u4uDgiIyN5+OGHg100n4RlEEgfgRAiFOzYsQOLxcLw4cN58sknQ7r5x5OwDAIhhAiWU6dO\nERMTQ0NDAw8++CAmk4lHH300bEMAJAiEEMJr9fX1XHrppZw8eZLIyEgaGxsBWz/BlClTgly6rpMg\nEEIIL73//vtUVVUB0NzcTFZWFpmZmWHRIexJWAaBdBYLIXqaruu8+uqrDBo0iH79+qEoCg899FBY\nNwk5hGUQSGexECKQrFYr9957L4cPH+b+++/HbDazZs0a9uzZw7x583j00UeDXUS/igh2AYQQItRs\n3bqV999/ny+//JK8vDwA1qxZA8CePXuCWbSAkCAQQog21q1b57w9c+ZMdF3n9OnTpKWlhdzKof4g\nQSCEEC6OHDnCli1buOaaa1AUhR9++IEffviBiooKFi9e3Cv6BNoKyz4CIYQIFMdM4WnTpnHw4EFW\nrVrFDz/8AEBOTk6QSxcYYRkEXR01dO211/L1118TGxvL6dOniYiIoF+/flRXV2O1Wp2TRBRFoU+f\nPtTW1qIoCsnJyVgsFhRFIT09nYqKCqxWK2eeeSYnTpwAICMjg8rKStLS0vjpp5+IiIhg2LBh7N+/\nn4iICIYPH86+fftQFIUxY8bw/fffoygK48aNo7CwEEVRmDt3LpMnTyY5OZm+fftitVppaGggPj4e\nRVECdj2FED/Lz7ftovt///d/REZG0tTUxEcffURCQgKlpaVkZGQEuYT+p+i6Huwy+OyWW27R33jj\nDZ+fN2XKFPbt20dMTIxzIkh0dDRNTU0AmEwmTp8+3e626zlRUVE0NzcDEBkZ6VxoynE7IiICq9Xa\nvb9gG5GRkURERBAVFcUZZ5xBfX09w4cPp6SkBEVRuOCCC9i7dy/R0dEsWrSIq6++msjISL+WQYje\nZtGiRWzfvp0VK1a0au4ZP3489fX1rFq1CoDf/e53fPPNN4CtRuDafxBuVFX9WtO0dhu+GyoICgoK\nyMvLw2w28+c//xmAq6++2m+3CwoK2r32pk2bALjmmmv485//jK7rTJ06lY8//hir1Up2djZ/+9vf\nsFqttLS0sH//fgBiY2OdG1q43nYNIk9MJhOpqamcPHmSqKgozj33XCoqKpg5cyZXXnklGRkZxMTE\n+HwNhegNmpubGTp0KM3NzUyaNIm33noLsA0bHTZsGLfffjuPPPKI8/zNmzfz5ptvOjebD1cSBGGg\noKCAp59+GvAucD744AMApk2bxubNm2lqaqKpqYny8nKgdU3GnaioKDIzM7FYLERHR3PxxRdz4MAB\nbr/9dq6++mri4uIC+dcVwu/y8/N5/vnnWbp0qccP7O3btzNr1iwAli5dyvz58wE4ePAg48ePZ9my\nZdx88809UuaeJEFgEB2Fyb/927/x8ssvU11dTUJCArW1tUDr5q22TCYT6enp1NbWYjKZUFWVgwcP\nSlCIkJWdnU1xcTGDBg1ix44dHZ736KOPsm7dOnRd5/bbb3dOENu2bRs33ngj77zzDpdddllPFbvH\ndBQEYdlZLDpmNpvZvHmz877rGigXXHBBh01jjmaq119/nZMnTxIXF0d9fT1Hjx51NkV9+OGHgK1t\nddGiRURGRhITE8PgwYOpqKjAZDJx8803c8UVVzBw4EAJCtHjzjjjDIqLizlw4AA5OTn893//d7ua\nga7r/PWvf2XixInU1dWxbds252OlpaWAbX9hI5EagWjFtR/Ftc9D13UmT57MH//4R2pqapxBAXTY\nQR4fH4/VaiUyMpKRI0dy9OhRYmJiuO2227jkkkuIj49n4MCB0lch/Obyyy9n1KhRbNu2jbq6Ored\nu/v27WPKlCmMHDmSsWPH8tZbb/HPf/6T9PR0HnvsMdavX09RUVGvHKnXq2oEsuhc4JjNZuc3KEdt\nwrVWcdlll7mtVWzatImmpiZOnz5NcXExAC0tLc7RWYWFhc7XaLtOS1xcnHP4bXR0NJdeeimHDx9m\n4cKFYb20r+hZFouFsrIyZs+ezYkTJ9i1a5fbVUF//PFHwBYIffr0AeAXv/gFjz32GKWlpQwZMqRX\nhoAnYRkEsuhc8LgLCtfbnjq833vvPRobG7Farc4qeGxsLPX19ZSVlTn7KhzPueWWWxg1ahRWq5V7\n7rmH66+/3nD/QYX3du3aBdj2Cz558iQ7d+5kwoQJ7c47evQoABMmTODee+/l5ptv5vvvvycvL4+y\nsjJGjx7do+UOBWEZBCJ0eeqj6Cgs2jZBXX755fzP//wPp06dYt++fei6zoIFC3jkkUeIjIykb9++\nLFmyhP/4j//o2b+cCGnffvstAOeffz4HDx6kqamJI0eOMHDgwFbnHT9+HIA333wTk8nE0KFDOXbs\nGHPnzuXmm2+mpaXF+TtpFBIEosd1FBauoZGdne1sgnrxxRexWCzU19fT1NREVVUV99xzD4sWLaJv\n377MmTOHX/7ylyQnJ/f430WEjm+//Zazzz6blJQUMjMzAdi/f3+7IDh27BgpKSmYTCbA1rT8v//7\nv/Tp0wer1cqBAwecv3tGIYvOiZBkNptZt24dubm5rF69mpycHB588EEuvPBCRo4cSUZGBqdPn6ay\nspIVK1Zw7rnnMmzYMGbNmsVVV13F559/Huy/gugmxyRKb33xxRc0NDRQUFDgHPXj6A9wdezYMc44\n4wzn/QsuuIC6ujpeeOEFIiIiMJvNYb/jmK+kRiBCnrt+CUfzkq7r1NbWUlxcTH19Pdu3bwfgxhtv\ndH4YLFmyhKuuuio4hRdd8u2333LNNdfQv39/li9fzsSJEz2e39LSwrFjxwDIy8vjT3/6EyaTibKy\nsnbnHjt2jPT0dOf9MWPGAPDpp58yZcoUjDgiUWoEIiw5mpc+/PBDnnrqKcaMGcOYMWN44IEHSEhI\nQNd1fvzxR0pLS5k3bx45OTlMnTqVt956i3AcMm00zzzzDM3NzZSXlzNnzhwKCgo8nu/YR3j48OHk\n5uYSGRnJoEGD3NYIjh8/3qpGMHToUBISEgAM+4VBgkCEPUcobN68mYULF/Laa6+Rk5PDww8/TFJS\nEgBFRUXs3buX3/zmN4wcOZIJEyawfv16CYUQtGPHDrZu3cpNN92EyWSisbHRuQBcRxwdwI5tJQEy\nMzOda3c56LreLggiIiLIzMxEURTn74vRhGXTkMwjEJ64NiWdd9555OXlcdlll7FlyxYOHTrEoUOH\nqK2tZfHixTz66KOkp6fz+OOPM2PGjCCXXAA8/fTTmEwmZsyYwYABA1i+fHmn80kcQeD6AT948GB2\n7NiBruvOYccWi4XGxsZW54Ft3S1d11m3bh1XXHGFn/9GoS8sg0DmEQhvuYbC/PnzKSgo4KmnnqKx\nsZHa2loOHjxIeXk5d9xxBwMGDGDJkiXccMMNQS61cVmtVr755huam5t57bXXeOmll1i+fLlzbayO\nOIIgLS2t1WvV1dWxefNmZ5OPu8AA28JzeXl5huskdpCmIWEoZrOZjz76iM8++4znn3+eCy+8kOHD\nh5OWlsaRI0dYuHAhEydOZMaMGWzdujXYxTWcoqIimpubGTVqFLm5ufTr149zzjmHf/zjHx6f5+4D\n/quvvgJ+3nQecHYou3YWw8+j1Iw0ZNSVBIEwLLPZzIcffsjf//53Xn75Zec8hP3797N7925uuukm\nfvGLXzBr1qxOOyuFf7h+eDs+lMePH4+maR6XVD927BixsbHOJSMAZs+eDdiWaXc9D+DMM8/0e9nD\nmQSBENhC4ZVXXnF2Mvfp0wdd19m6dSvbt2/ntttuY8uWLcEuZq/31VdfkZqa2mr1z/Hjx1NXV8f1\n11/fYSAfP36c9PT0VkuQOPoVXD/0O6oRGJ0EgRB2juaB+fPn8+qrr5KTk8O8efOIjo6mvr6eO+64\ng0mTJvHZZ58Fu6i9lqZpjBs3rtUHuqraFsv8+uuvycvLc/s8RxC4cvQXVFZWtjovNjaWxMREfxc9\nrEkQCOGGIxQcG5hcfPHFJCQkUFJSwq233spjjz3msalC+O748eOUlZVRVFTU6pv/gAEDiI6O5uyz\nz+6wM7eioqJdECQkJBAbG0tFRYXz2NGjR9vVHIQEgRCdMpvNbNq0ibVr19K3b190XWft2rWcc845\nvPjiizIXwU8cq4eWlJS0+uavKApDhgzh3HPP7bAzt+1sYcfzUlNTWwVBUVERFotF+nzakCAQwktm\ns5m1a9cyefJkhgwZQn19PcuWLWPSpElce+218uHSTY59LC6//PJ23/zPPvtst7OEwbYR/YkTJ9oN\nCQVb85Br01BpaSk1NTUdNjEZlQSBED4wm82sX7+eZ555hn79+gG2UUZffvklt956K5988kmQSxi+\niouL6devH2+//Xa7b/6ZmZmUl5e7rX1VVlai63qrOQQOqampzqGlYJs4NmDAAMPOF+iIBIEQXWA2\nm8nLyyMnJ4f77ruPmJgYGhsbueuuu5x7LAjflJSUMGzYMLePDRo0iFOnTjnXFHLV0SQxsI0OcjQN\nWa1WTp48yaxZsww7X6AjYTmzWJaYEKHAddbyJZdcwrJly9i9ezeFhYX8+te/ZtSoUcyfP18+dLxU\nXFzM9OnT3T7m2F/gxx9/JCUlpdVj7mYVOziahnRdp7q6GqvV6vY8owvLIJAlJkSocYTCli1bWLJk\nCUePHuXo0aM0NjZKEHihqqqKioqKDmsEZ599NgAHDhzgoosuavWYpxpBamoqTU1N1NTUOM9LTU31\nZ9F7BWkaEsKPpk+fTmFhIQsWLCAiIoIvv/ySlStXBrtYIa+kpASwLQntzqBBgwBbELSlaRrwc2ez\nK8e3/4qKCmcTkdQI2pMgECIAFi9ezPjx47FarSxfvpwJEybI2kUeOIKgoxpBQkICaWlpboPAsRvd\na6+91u4x10llEgQdkyAQIkAWLlxIdnY2qampHDhwgFtvvZW//OUvwS5WSCouLsZkMjmbgNwZNGiQ\n2yAYPnw4JpOJu+++u91jjmagiooK5zBSWV6iPQkCIQLEbDbz5ptvsmrVKhISEmhubuaee+6RUUVu\nfPHFF0RHR3tcZTQuLg5N09pdu9jYWDIzM932xbRtGoqIiHAuLih+JkEgRICZzWZee+01xo0bR0RE\nBIWFhdxxxx0SBi727t1LbW2tx4leBw4coKGhgT/84Q+tjldXV3f44e4YYeRoGkpNTSUiQj722pIr\nIkQPMJvNvP/++6xatYqoqChqa2t55JFHgl2skKDrOs3NzQwcONDjRK/s7GwAbrvttlbHLRZLh1tM\nmkwmkpOTnTUC6R9wT4JAiB505ZVX8vrrr5OamsoPP/zA5MmTDV8zqKqqorGxkTvvvNPjUNtx48YB\n7TuUPQUB2JqU/vKXv1BaWipDRzsgQSBED8vJyUHTNNLS0igqKuL22283dBiUl5cDeOwohtbt/a48\nNQ2BbUZxRUUF5eXlUiPogASBEEFgMpl44YUXiI6Opq6ujvvvv5/m5uZgFysoHEEwcOBAj+e5C4Lm\n5mZOnjzpMQiuu+46ABoaGqRG0AEJAiGCZPLkyfzpT39i4MCBHDx4kEsuucSQNYODBw8CndcIXIeC\nOtTU1AB4bBq6/fbbnbdl6Kh7EgRCBFF2djZffPEFgwcP5ujRo9x2222GC4Py8nISExM9fpiD+yCo\nrq4G8FgjOOuss5wzlqVpyD0JAiFCwLJlyzCZTNTX1xtuHa3y8vJOm4Xg5xFArvsLWCwWwHONAGwL\nU4KsM9QRCQIhQoDZbOaPf/wj6enpFBcXk52dbZiawcGDBzttFnJou7+Ao0bQWRD0798faL1/sfhZ\nWAaBYxlqx8+2bduCXSQhum3y5MlomuYMg9/97nfBLlLA6brudY0A2u845qgRODYJ6sgXX3wBwObN\nm7tY0t5NlqEWIoRERUXx3HPPMW/ePHbv3s3bb7/NrFmzgl2sgKmqqqKurs6nGkFRUZHzvrc1AsdE\nNdmZzL1Og0BRlOuB6UA/4ASgADqwRdf19wJbPCGMZ9q0aYwZM4Z//OMfLF68mLfffpuFCxf2yn0N\nHIvwnTp1yqvz09LS2LFjh/O+t0HguomQaK/DpiFFUS5SFOVO4Gtd13N1XZ+t6/rd9tt3AzsVRble\nUZQxPVdcIYxh0aJFqKpKS0sLO3bsaLe+Tih54403mDJlChs2bHC7p7An69evB+Dvf/+7V+enp6dT\nVVXlnHNhsViIj4/HZDL59L6iNU99BNW6rq/Vdb3M3YO6ru/Xdf1dwBKQkglhYGazmQ8++IA777wT\ngPPOOy/IJWqvubmZuXPn8uCDD7Jv3z7uu+8+Ro4c6VNHt6qqAMybN8+r8x2jfk6cOAHYagSd1QZE\n5zoMAl3X9ztuK4qywfHNX1GUwa61ANfzhBD+9cgjjzBs2DBeeeUV3nsvdFpirVYr9913Hx9++CFg\na5o566yzqK2tpbi4mCVLlmC1Wjt9nbi4OGJiYrjyyiu9el/HPADHyCGLxSLLSvuBt6OGVuu6/k9F\nUYYAucDsAJZJCGHnWD+/qamJBQsWcNVVV4XEsNJ7772Xd999l+zsbHJycli9ejUrVqzg/PPPp1+/\nfpSVlXm1oN7hw4cZMGAAiqJ49b6uO45B5+sMCe946iP4RFGUPyiKch1QYv9zKpCn6/pDPVZCIQzu\ngQceYPz48URFRbFz505+//vfB7U83377Le+//z4AiqKwbt06Z2fsxx9/zK5du+jbty8lJSUe9xcA\nOHLkCAMGDPD6vdvOLu5s5VHhHU81giXAg9hGCeW6/MyUDmIheo7ZbOa9995z7l/g2GwlGJqbm7n7\n7ruJiorivPPOc7s9pKIoTJ06FYA5c+Z4fD1fg8CxVtDKlSspKCiQPgI/8dRHsFPXdYuu6+/quv6g\nruv/ouu6CnwKDO25IgohAO644w6ys7PZtGkTb7/9dlDKsGXLFsrKymhubiY9Pb3DIZk333wzgMd+\nAqvVyk8//eRTECQmJhIZGemsbVRWVvL//t//C4nmsnDmqWko0d1xe0C829l5Qgj/a25uRtd1nnji\nCZ+HavrDW2+9Rb9+/cjOzvY4Oeviiy8mPj6erVu3dnhOZWUlTU1NPgWBoigMGDCAtLQ0brnlFpqa\nmigvL++0CUp45qlpaJy9X6BD9slmqn+LJIToyK9//WuGDRtGVVUVV1xxRY9+Ez5y5Ah/+9vfuOmm\nm3jzzTc9TtAymUyMGjWKt99+u8MyHjlyBMCnIAAYMWIEAwYMICMjA7ANrZUZw93T4cxiXdc/VRQl\nSVGUB7A1BbX9+lGCbTRRTSALKIT4mdls5tNPP2X06NHs3r2bl156qcdmzD733HNYrVaGDx/u1fk1\nNTXU19ezcuVKt2V0BIHjA91bGRkZ7Ny5k7KyMgCef/75kJxnEU48Dh+19xEsd8wmbvPznISAED0v\nKiqKpUuXoihKjzUP6brOpk2bAJwjhjrzy1/+EoBLLrnE7eNdrREMHDiQqqoq9uzZA8DgwYN9er5o\nz2MQKIpyZ2fNQ0KInjdnzhx+85vfsH37dqZPnx7wJqK9e/dSV1fHiBEjvG6GufHGG4mLi+twHaHD\nhw8THR3t8x4BjpVKd+zYQWpqKn369PHp+aI9T53FecAw4F8URflrzxVJCOGNX/3qVyQmJvL999/z\nwgsvBPS9Nm3aREREBO+8847XTVHR0dFcdNFFfPXVV24fP3LkCP379yciwrfV8M866ywAdu7cSWZm\npk/PFe55+hf42j5sNBeYJTUDIUJLVFQUTzzxBIqisGvXroDVChzNQpdddpnPWz2qqsp3331HXV1d\nu8f27t1LTU2Nz+V2BEFTU5MEgZ94CgLn7g+6rluQxeWECDk33HADQ4cOdW7WFAjfffcdZWVlHD9+\n3OcP7XHjxtHS0sLOnTvbPVZSUoLFYvF56OeZZ55JZGQkgASBn3gKgmftS0zcaZ9J7Jy+J3MHhAgd\njz/+OHFxcVRUVNDU1OT313fsALhv3z6fP7THjh0LwJIlS1qFSF1dHY2NjQwZMsTnoZ9RUVHODmYJ\nAv/wFATLgNXYNqR5GFswFCmK8gdgTSALpapqlqqq76iqelcg30eI3mDy5MmsWrWKw4cPM2PGDL83\nERUUFDBw4EBycnJ8/tBOTk4mISGB/fv3twqR0tJSAB588MEuDX91NA/JiCH/8LTExBr7LOLluq7P\n0nV9GLbJY+9iC4dOqao6ts39maqqTlNVdbEXT5+radpqb95HCKObPn06ffr0Yd++fX7dxKaxsZEv\nvviCGTNmOBeX89Wll15KVFRUqz0HiouLARg2bFiXyhUdHQ3ATz/91KXni9Z86q63zyvIBzrdRUJV\n1WnAOy73xwJompYPVLcNiTZOAFn24JA1ZoXohKIoznH7juYYf/j6669paGjo1qS1q666iubmZs44\n4wznseLiYiIiIrr8jd4RABs2bOhyucTP3AaBoihDPK0w6s1mNPYP/FKXQ7OBavvtUmAaOGsJrj/J\ngKppWqH9PGkeEsIL999/P2eeeSZr1qzxW/PQ559/TmRkJBMmTOjyazh2IXMdRlpcXMygQYOIjY3t\n0ms++eSTZGdnu139VPjO7RITuq7vt+9ZPBfb0hKaruv/7OZ7JWP7pu+QCqBp2sa2J6qqWmqvMahA\nu8eFEO1FRkaSlJTE0aNHefrpp9m8eXO3X/OTTz4hPj6eb775psu1giFDhpCamoqmac5VSYuLixk6\ntOuLGE+cOJGJEyd2+fmiNU9rDe0EdoJzI3tHKJTquv5ZIAulaZqjJlEYyPcRord55JFHuOOOO/yy\n9ITVaqWoqIiWlhby8vK6HASKojBu3DhnjaClpYX9+/czadKkbpdR+IdXfQT2TuM1uq6vBfYrijLX\nsfyED0NJqwHHjhrJuMxTEEL4x5QpU7jvvvv45ptv+Pd///duNREdPHiQlpYWRo4c2e3VPdPS0igr\nK+P999/n0KFDNDQ0dLmjWPhfhzWCjtj7B9YAKIqShK2t35tdtTfw85LVWUC+r+/tcPz4cWe7oztz\n585tNUJBCCO59dZbefbZZ9E0rVvf5L///nvAtupodzugHaOEli9fzg033ADgdraxCA6fg8CVfcax\n2xBQVXWm7Q91pqZpGzVNK1RtpgHV9s7gLklPT+fjjz/u6tOF6NWSk5M5//zzKSoq6tY3+T179qAo\nCqNGjep2mRYtWkRubi6nTp1yzif49NNPufPOO7v92qL7OgwC+6YzHc0XUOx/Vuq67jYI7J3AG9sc\nk3kBQvSA2bNns3TpUvr379/l1/j+++8ZPHgw8fHx3S6P2Wzm6aefdgbT6NGjueeee7r9usI/PHUW\nv9vRY0KI0DZ9+nSWLl1Kfn5+l9vi9+zZwznnnOO3Ms2YMYPkZNu0oMcee6zHNtQRnfNt/dcQ4Vhg\ny/HjWAtFCGFz1llnkZmZyYsvvtilDuO6ujrKysoYPXq038pkMpkYPXo01dXVssdwiOm0j8DeRDQd\nWzPRCWzNQjqwpaNmoUDr27cvS5cuDcZbCxFWLBYLv//9733+9r1hw4aA7H62YMECYmJiZI/hEOOp\nj+Ai4GIg310zkX328fVAiR9gRxujAAAPSklEQVQmmwkh/GzOnDk88cQTXHbZZT4/d926dYBtZvGi\nRYv8Viaz2SxNQiHIU9NQta7ra3VdL3P3oK7r++0BIfsUCBGC5syZQ0xMDFVVVT4/d+TIkSiKwoIF\nCwJQMhFqPHUWO9cTUhRlA/C0ruv/VBRlMJDsqAV4s+6Qv7XdhGPixIkyS1GINmJjY1FVle3bt/v8\n3MbGRoYNG0Z2dnYASiZCjbfzCFbbQ2AItpVHdSBozUHSRyCEdy6//HKeffZZTpw4QUpKSudPsCsp\nKSErKyuAJROhxNPm9Z/Ydyi7Diix/zkVyNN1/aEeK6EQossuv/xyAP7zP//T69FDzc3NlJWVdWtR\nOBFePPURLAEexDZKKNflZ6anJaqFEKHjwgsvRFEUvvvuO6+HbJaXl9PU1CRBYCCdrT4Kth3JnKOG\n7KOJhhLEpiEhhHeio6MZMmQIVVVVXg/ZLCkpAZAgMBBPw0cTdV2vaXvcdXlqT+cFknQWC+G9CRMm\nsHnzZmczUWckCIzHU2fxOEVRkjxNGrPPI6gCAro/QVvSWSyE98477zzWr1/PoUOHGDhwYKfnl5SU\nkJyc7FPnsghvnpqGPlUUJUlRlAewNQW1nWZYgm00UY/WBoQQvjnvvPMA2L17t1dBsHPnTqxWKwUF\nBTL5yyA8Dh+1LzO9vIfKIoQIgHPOOYfIyEh2797NFVdc0en5JSUlNDY2dmsvAxFePC4659iFrKcK\nI4Twv7i4OIYNG8auXbs6Pbe5uZmmpiYyMzNlPSAD8TSPIA8YBvyLoih/7bkiCSH8LT09nW3btnU6\nl+Do0aNYrVbmz58vtQED8dQ09LWu684tKRVFuS5Yq422JaOGhPDNkSNHOH36NC+99JLHD/jDhw8D\nkJGR0VNFEyHAUxA4N5fXdd2iKErILC4no4aE8M0NN9zAM888w/Tp0z2e5wiCs846qyeKJUKEpz6C\nZ+1LTNxpn0mc5HhAUZTEwBdNCOEv111n6+ozmUwezzt06BAgNQKj8VQjWAZowDTgYWCsoijLgHwg\nBZgd+OIJIfwhIyODhIQEfvjhB4/nHTp0iKSkJPr27dtDJROhwNM8gjX2m66ziJOAccDiAJdLCOFH\niqIwYsQI9u3b5/G8Q4cOSW3AgHzas1jXdYuu6/nYlqIWQoSRESNGeFUjkCAwHp83r7dvTCM1AiHC\nzIgRI6ioqODEiRMdnnP48GHpKDYgbzemcXQQP4utz6AfcHegCtUZGT4qhO9GjhwJwL59+5gwYUK7\nx2tra6murpYgMKBOg8AlAKYCy3Rdz1UURQt4yTyQ4aNC+G7EiBEAPPTQQzz55JPt5hPI0FHj6myJ\niWeAr4FPdF0fruv6WvtDbRegE0KEuIyMDBRFoaioyO0mNZ988gkAlZWV7R4TvVtnfQRbsO1Slt8D\nZRFCBJCiKCQnJ9O/f3+36wi9955t4YCPPvqop4smgqyz1Uc/BVAUZap96Gi+LDstRPjKyMigf//+\nbpeZuOSSS9i7dy9z584NQslEMHnVWdwmEJKx7WMshAgzqampHY4aSk1NBeh0GQrR+3g9aghaBUJg\nSiOECKjU1FTKysrcPlZTU0NiYiKRkZE9WygRdD7PIwDQdf3dzs8SQoSalJSUDjuDq6urSUyUZcSM\nyNPm9ddjmy/gyYlgLE0t8wiE6JrU1FRqa2tpaGggNja21WM1NTUkJSV18EzRm3laayhkv/XLPAIh\nusbRD1BZWdluvoDFYpEagUF1qWlICBGeHEHgrsO4pqaG5OTkni6SCAESBEIYiGuNoC3pIzAuCQIh\nDCQlJQVwHwTSR2BcEgRCGEhHNYKmpibq6uqkRmBQEgRCGEhSUhKRkZHtgqCmxrZggPQRGJMEgRAG\nEhERQUpKSrvO4urqagCpERiUBIEQBpOWltZhjUD6CIxJgkAIg3E3u9hisQASBEYlQSCEwaSmpkoQ\niFZ8WnQuVMgSE0J0nbsVSB1BIH0ExhSWQSBLTAjRdbW1tVgsFrZu3Up2djYgfQRGJ01DQhhMYWEh\nAKtWrXIes1gsmEymdgvRCWOQIBDCYKZOnQrAjTfe6DxmsVhISkqSvUYMSoJACIO58MILARg9erTz\nmKw8amwSBEIYTEJCAmDrK3CQdYaMTYJACINxFwSOpiFhTBIEQhhMnz59AKirq3Me++mnn/juu+8o\nKCgIVrFEEEkQCGEw8fHxQOsaQVVVFceOHSMvLy9YxRJBFJbzCIQQXeeuaSg6Opr09HRyc3ODVSwR\nRFIjEMJgHEFw6tQp57GWlhauvvpqzGZzsIolgkiCQAiDcTQNOfoIdF2nsbGRmJiYYBZLBJEEgRAG\nExkZSWxsrLNpqLm5GV3XMZlMQS6ZCBYJAiEMKCEhwRkEjY2NALK8hIFJEAhhQH369GkXBNI0ZFxh\nOWpIlqEWonvi4+OdfQQSBCIsg0CWoRaie9w1DUkfgXFJ05AQBpSQkOAcPio1AiFBIIQBJSQkSNOQ\ncJIgEMKA4uPjnU1Dp0+fBiQIjEyCQAgDcu0jaGhoACQIjEyCQAgDkqYh4UqCQAgDSkhIoLGxkaam\nJmkaEhIEQhiR6wqkMnxUSBAIYUDugkCWmDAuCQIhDMgRBHV1ddJZLCQIhDAi113KHH0E0jRkXBIE\nQhiQu6YhqREYlwSBEAYkQSBcSRAIYUBtgyAmJgZFUYJcKhEsEgRCGFDbIJD+AWOTIBDCgNzVCIRx\nSRAIYUCuG9hLEAgJAiEMKCoqiujoaDZs2MChQ4ekacjgQnaHMlVV7wJKgROaphUGuzxC9Ebl5eXU\n19eTnp4e7KKIIApoEKiqOtb1Q1xV1ZlANTBW07RnPTzvLiBf07TSQJZPCCNLS0ujpaWFjIyMYBdF\nBFnAmoZUVZ0GvONyfyyApmn5QLXjfgcuBsaqqjpTVdWsQJVRCCPr378/o0ePJi4uTvoIDC5gNQJN\n0/JVVXX9Rj8b2GK/XQpMAwrttQRX+dhqDflACjAPWBKocgphVImJiVgsFiIiIpydx8KYerKPIBk4\n4XI/FUDTtI1tT1RV9WlsQQHwSuCLJoTxJCYmcvDgQeLi4ujXr1+wiyOCKCQ7izVNqwbaBYQQwn8S\nExOpqakhIiJCmoYMrieHj1Zja+oBW+2gsgffWwjRRlJSEjU1NTKPQPRojWADoNpvZ2HrA+iS48eP\no6pqh4/PnTuXefPmdfXlhTCExMREGhsbqampkSAwuIAFgb0TWFVVdaamaRs1TStUbaYB1d2ZG5Ce\nns7HH3/sv8IKYUB9+/YFoLq6WoLA4AI5amgjbdr5NU1bHaj3E0L4JikpyXlbgsDYQrKzuDMnT57k\nqaeect6fOHEikyZNCmKJhAg/iYmJztuyxISxhWUQ9O3bl6VLlwa7GEKENdcgkBqBscmic0IYlGvT\nUGxsbBBLIoJNgkAIg5IagXCQIBDCoFxrBNJHYGxh2UcgncVCdF9sbCzR0dE0NTVJjcDgwjIIpLNY\niO5TFIXExEQqKyslCAxOmoaEMDBHP4EEgbFJEAhhYI5+AgkCY5MgEMLApEYgIEz7CKSzWAj/cKw3\nJEFgbGEZBNJZLIR/OJqGZPiosUnTkBAG5mgakpnFxiZBIISBSR+BAAkCIQytpqYGgO+++y7IJRHB\nJEEghIHt2bMHgHXr1gW5JCKYwjIIHKOGHD/btm3z+rmvvPJKAEvW+8j18k24Xa/58+eTk5NDbm5u\nUN4/3K5XsAXqeim6rgfkhQPplltu0d94440uPVdVVTRN83OJei+5Xr6R6+UbuV6+6e71UlX1a03T\n2m34HpY1AiGEEP4jQSCEEAYX9kHgrn+g7TFf+hD8XRZ/PcfTeR095s21cXdMrpdvx0L5enn7PH9d\nL3fH5Xp5fiwUrlfYB8Hnn3/e6TF35/RUWfz1HE/ndfSYN9fG3TG5Xr4dC+Xr5e3z/HW93B2X6+X5\nsVC4XmG5xITrWkM7d+5k27ZtstaQEEJ0UVgGgetaQ0899ZSEgBBCdENYDh9VVfU48KP9bhJgaXNK\n22Ou9zOAwwEqmruy+Os5ns7r6DFvro27Y3K9fDsWytfL2+f563q5Oy7Xy/NjPXm9MjVNS297MCyD\nIJBUVX1F07R5wS5HqFNVNRlQgbHAak3TqoNcpJCmqmoWkAxMAzZqmlYa5CKFPFVVxwIpmqblB7ss\nocjl/2AykN+d/4Nh31nclv2Xx/X+TFVVp6mqutjX5xpBN66XCmhANZAVqPKFmm5cr7FAKZAPzAxU\n+UJNd/4/ArOxfcgZkhfX7iF7SOYDd3XnvXpVEKiqOg14x+X+WAD7xar29EFv/8ZWbf8xhO5cL5dv\nacmaphUGtKAhopvXa6P9G9s0YGOgyxoKunO97DYEsHghzctrl2w/Vg0M7c77hWVncUc0TctXVdW1\nyj0b2GK/XYrtP2Ghqqptv5HlY/tWewLIUlU1ywhV925er1mapq1WVbVQVdXFmqY92wNFDqruXC9N\n06rt/7nzjfC7Bd2/Xj1RxlDl5bWrtjcPAZR05/16VRC4kYztw90hFWzfztycm2+/qIatiuLD9VJV\nVbN/sGVhkG+4bvhyvaYBS4BSVVW3dPA72Nv58v8RbB92qaqqGj4YcH/tXsF2jQBWd+fFe3sQ+MT+\nyzY92OUIB0ZpDvIXl7Zc4SUj1DK7w16z9Evtslf1EbhRDaTYbycDlUEsSziQ6+UbuV6+kevVdQG9\ndr09CDbw84iWLOQbWWfkevlGrpdv5Hp1XUCvXa8KAnunk+rofHI0X9jbZ6ulOaM1uV6+kevlG7le\nXdfT104mlAkhhMH1qhqBEEII30kQCCGEwUkQCCGEwUkQCCGEwUkQCOFniqIYZhE+0TtIEAjhI0VR\nxiqKsqWjx7AvU+LpPCFCiQSBEL7zNLVf1XW90IvzhAgZEgRC+G4aP68E6Y/zhAgqCQIhfDcOKFQU\nZZqiKM414+3NQlqb81Ls53Vr4xAhAklWHxXCd2OBsbqub1QUxXVpYFXX9dVtznvafvuGHiudED6S\nGoEQvkvB9k0/y6U/wJ1qXdcdu5JJE5EIWRIEQvjA3vyzAdvqj/McQ0Xtf5a2Oe8r+93p2JqSZFip\nCEkSBEL4JgvI13Xd8aHv2NFumq7r+W3Ps98uwdaUJCOIREiS1UeF8ANFUe5q0z8gRNiQGoEQ3aQo\nStv9ZIUIKxIEQnTfNGS3LRHGpGlICCEMTmoEQghhcBIEQghhcBIEQghhcBIEQghhcBIEQghhcBIE\nQghhcP8fNjP71qaOwbkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc3b15e2510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog( kh, np.abs( (PCLASS/(cosmo.h())**(-3) - pk[0,:])/pk[0,:]   ) )\n",
    "plt.xlabel(r'$k/h$')\n",
    "plt.ylabel(r'$|\\Delta P(k)| / P(k)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear P(k) Difference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "pars.NonLinear = model.NonLinear_both\n",
    "results.calc_power_spectra(pars)\n",
    "kh_nonlin, z_nonlin, pk_nonlin = results.get_matter_power_spectrum(minkh=1e-4, maxkh=1, npoints = 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "PCLASS  = np.array([cosmo.pk(ki*cosmo.h(), 0) for ki in kh]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0,0.5,u'$|\\\\Delta P(k)| / P(k)$')"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAECCAYAAADzStBRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xt4FfWdx/H35EYCFEICgoJVAlK0\nCBhHqjVCY6PtrmWVluI+XrY+VgxVFG0fFeh217VFij5Vn8IuAfVxbRHlTlGESmxV2O6jjuBquYgQ\nuYtc4gkJIRdyZv84Fw8hOcnJuc6Zz+t58nhmzmV+jMl8z/f7u4xh2zYiIuJeGclugIiIJJcCgYiI\nyykQiIi4nAKBiIjLKRCIiLicAoGIiMtlJbsBXfHd737XPu+885LdDBERR9m+ffsxy7L6td7vyEDw\nta99jeHDhwe3r7nmGsaOHZvEFomIpD7TNPe2td+xgeCXv/xlspshIpIW1EcgIuJyjswIamtrmTVr\nVnBbpSERka5zZCBQaUhEJHYcGQiUEYiIxI4jA4EyAhGR2FFnsYiIyykQiOs8+uijXHPNNWzatCnZ\nTRFJCY4sDamPQKKxbNkyPB4PFRUVlJSUJLs5IknnyECgPgKJRmFhIc3NzUyZMiXZTRFJCSoNiesY\nhkH//v2VDYj4KRCI69TV1XH69OlkN0MkZTiyNKQ+AonGyZMnMQwj2c0QOcvp06f58ssvOX78OG+/\n/TYvvvgiLS0tDB8+nI8//hjbtnn88cf5h3/4h5ge15GBQH0E0lW2bXPy5Elyc3OT3RRxKdu2OXjw\nIMuXL2fFihX07t2bTz/9lKamJpqbm7Ft+6z3HDhwIPh44cKFCgQi0Th16hRer5fm5uZkN0Vc4tix\nY7z77rusXr2at99+m4aGBlpaWoLPG4YRvPjn5ubS0NAAwJAhQ8jKyiIzM5Mf/OAHrF+/noyMDH7x\ni1/EvI0KBOIqdXV1AOojkLipq6ujoqKCJUuW0NjYyPHjx4EzL/jdu3envr6ePn36MHXqVF599VUA\nxo8fH3w8Y8aMMwY0TJs2LW5tViAQV1EgkHg4evQo8+fP5+WXX6a2tjZ4wc/Ozg6+5tJLLyUjwzc+\nZ/z48WzatIkpU6ZQUlJyxlDmZAxrViAQVzl58iSASkMSterqap588klWrFgR/L0KGDJkCIMGDWLs\n2LHtfsNPpXksCgTiKoGMoKWlBa/XG/yGJtIZXq+XiooK5s2bR01NzRnPDRs2jLy8PAzDOOOin0oX\n/PY4MhBo+Kh0VSAQgC8r6NatWxJbI07xxRdf8MQTT7B69epgZy7AN77xDfLy8oCzv/E7iSMDgYaP\nSleFBoLTp08rEEhYa9eu5bHHHuPQoUN4vV4ALrjgAnr16kVmZqajL/6hHBkIRLqqdUYg0pYvvviC\nxx57jNWrVwPQrVs3GhsbKSgo4IknnkiLi38oBQJxldBOPY0cktY+/PBDZs+ezf/8z/8ER/707t2b\n+++//4xRPukmJQOBaZr5gAkUAwsty/IkuUmSJpQRSFveeOMNZs6cyeeffx7cN3z4cM4999w2h3im\nm7gGAtM0iy3L2hyyPRHwAMWWZT0R7q2ABRT5fzaHea1Ip7XuIxDZsGED5eXlNDU1BWf2FhQU8B//\n8R9p+e2/LXEbO2eaZhmwLGS7GMCyrErAE9hui/81APmhgUQkWqGlIWUE7lZdXc19993HHXfcQUZG\nBsOGDeOhhx6itLSU+fPnuyYIQBwzAsuyKk3TrArZdTOwwf+4CigDNvuzhFCVwCTLshaaprnZNM2H\nO8geRDpNGYEcPHiQRx99lHXr1mHbNgMGDODw4cMMHDiQKVOmpHUJqD2J7CPIB6pDtgsBLMta3vqF\npmla/oyiCDjreZGuUh+Bex07dowHH3yQv/71r2es8DlgwAAuvvhiVwaAgJTsLFY5SOJFo4bcae7c\nuTz55JPBVT8vueQScnJyAGdPBIuVRAYCD1Dgf5wPHO/qBx09ehTTNNt9fvLkyZSXl3f14yWN1dXV\nkZ2dTXNzszICl5g9ezbz5s0Dzh4JJD6JDARL8I0GAl/JpzLMa8Pq168f69evj0mjxF1OnjxJnz59\nOHLkiDICF3juueeCQaBPnz6uGgkUibgFAn8nsGma5kTLspZblrXZ9CkDPNGUf7TWkHRVXV0dBQUF\nHDlyRBlBGrNtm6lTp7J69WqGDx/OOeecw7333qsg0I54jhpaTquOXsuyFsbis7XWkHRVXV0dRUVF\ngPoI0tWpU6e46667eOuttwBfZ/BLL72U3EaluJTsLO6IMgLpitOnT9PQ0EB+fj6gUUPp6LnnnmPO\nnDnU19cDUFBQwM9+9rMktyr1OTIQKCOQrgiMGFIgSD/bt2/n4YcfZvNmX8W5qKiICy64QJ3CneTI\nQKCMQLoiMIcgEAhUGnK+U6dOMXXq1OCN3cHXKTx79mwFgAg4MhAoI5CuUEaQXg4ePMgPf/hDDhw4\nAPjmBvTr109ZQBc4MhAoI5CuCGQEq1atApQRONnSpUuZPn16MJgXFBTwq1/9SgGgixwZCJQRSFfU\n1tYC8OmnnwLKCJzq2LFjzJgxg8bGRr7xjW9w3nnnKQuIkiMDgUhXBG42Pnr0aD788ENlBA7U0NDA\nxIkTaWxsZNiwYTz22GMKADHgyECg0pB0RSAQPPXUU1x77bXKCBzGtm1+8YtfBDO6gQMHKgjEiCMD\ngUpD0hUnTpwAoLCwEFAfgdM8+OCDrF69mtLSUgBXrxYaa44MBCJdceLECbKzs/na174GqI/ASf70\npz+xbJnvPle2bWumcIwpEIhr1NTU0Lt37+Dyw8oInGHz5s1MmzaNHj16MGTIEM0UjgNHBgL1EUhX\n1NTU0KtXLwzDIDMzUxmBAxw8eJA777yTzMxMTp48SWFhofoF4sCRgUB9BNIVJ06coHfv3gBkZ2cr\nI0hxx44dY8KECVRXVzN58mQ++eQT9QvEiSMDgUhXBDICgKysLGUEKezYsWOMHz+egwcPAvDJJ5+w\naNGiJLcqfWUkuwEiiRKaEWRlZSkjSFG7du1i/PjxwaUjCgoKlAnEmSMzAvURSFeEZgSB21VKann/\n/fe59dZbaWpq4o477uCzzz7TrOEEcGQgUB+BRMq2bWUEKe61115j6tSptLS04PV6+eyzz1QOShCV\nhsQVGhoaaGpqUkaQombOnEl5eTnZ2dl4vV6VgxLMkRmBSKQCs4pDMwIFgtQwb948XnzxRQCGDBlC\n3759VQ5KMAUCcYXAOkOhGYFKQ8ll2zY///nPWbp0KSNGjKCwsJB77rlHASAJFAjEFQKBQBlB6pg2\nbRorVqwAoG/fvlo2IokcGQg0akgiFSgNKSNIDYsWLQoGAd1gPvkcGQg0akgipT6C1PFf//VfzJo1\niyFDhjBo0CCVg1KAIwOBSKRal4aUESTHwYMHmTNnDgDnnXceixcvTnKLBDR8VFwiEAgCS1ArI0i8\n5uZmbrnlFrxeL8OHD2fq1KnJbpL4KSMQVzhx4gS5ubl069YNUEaQDE899RS7du0C4Nxzz1U5KIWk\ndEZgmuaCZLdB0sOJEyfIz88PbisjSKznn3+e3//+94waNYrS0lJNFksxcc0ITNMstixrc8j2RMAD\nFFuW9URH741n28RdPB5PcMQQaImJRGpqauK3v/0t4Bu1pWUjUk/cMgLTNMuAZSHbxQCWZVUCnnAX\netM0i/AFDE+82ifusm/fPg4fPsymTZsAyMnJUUaQINOnT6e+vp7BgwerXyBFxS0jsCyr0jTNqpBd\nNwMb/I+rgDJgsz9LCFUJFAHVQJFpmkWWZVUhEoW9e/dSW1tLRUUFJSUlyggSpKqqiqVLlwJw4YUX\nql8gRSWyszgf38U9oBDAsqzlbby20jTNfP97RKLWs2dP8vLygrVpdRbHn23bzJgxg9zcXC677DL1\nC6SwlB01ZFmWB7gu2e2Q9OD1eikrKwt+I1VncfzNmTOHTZs2ceedd/LrX/862c2RMBIZCDxAgf9x\nPnC8qx909OhRTNNs9/nJkydTXl7e1Y+XNFRbW0vPnj2D28oI4quhoYGFCxcCvvKQpLZEBoIlQODq\nXYSvL6BL+vXrx/r162PSKEl/LS0t1NfXByeTgTKCeHvhhRdobGxk9OjRWkfIAeIWCPydwKZpmhMt\ny1puWdZm06cM8IQOK42UFp2TSJw8eRKAHj16BPcpI4if6upqnnrqKQoKCpgxY4Y6iB0gnqOGlgPL\nW+1bGIvP1qJzEom6ujqAszKCpqYmbNvGMIxkNS0tzZ07l/r6eurr64OjtCS1pWxncTjKCCQSgUDQ\nOiMAX9koK8uRfwYpaeXKlTz77LOMHj2aPn36aKSQQzjyL0AZgUSivYwAfAuhKRDEzqxZs7Btm+7d\nu2sGsYM48i9AGYFEIhAIWo8aAtRPEEP79u3jyJEjDBo0iGnTpiW7ORIBRwaCrmYEe/bsoa6uDsMw\nyMjICP6Ebrf3OLDdel97+wP7JPlqa2uBMwNBaEYgsTFv3jyysrJYvXo15557brKbIxFwZCDoakZw\n77338uGHH8azae3KyMjAtm3AdxFqaWkBfGveNDU1YRgGeXl5NDY2YhgGvXr1CgatwsJCvvzySzIy\nMujfvz8ej4dzzjmHzz//nMzMTC666CKqqqrIzMxk5MiRbN26lYyMDK666ip27NjBxIkTufLKK8nL\nywv+9OjRg8zMzKSci0QLVxpSRhAbhw4d4pVXXqF///7s3r1bgcBhHBkIupoRZGT41tjr2bPnGR2I\nocML23rcvXt36uvrz3qcl5fHqVOnzngcui83N5eGhgbAd8EPPM7MzAx+E7VtG6/XC/i+nQb219TU\nBB8fPnw4eMGqr6+npaWF6urq4PuOHTsW/DceOnQo+Hjv3r0A7Qa/jIwMevToQXNzM5mZmQwYMIDj\nx4+TmZnJkCFDOHLkCEOHDmXHjh1kZGRw5ZVXsnPnTiZNmsTVV19Njx49gj+pHFTCdRYrI4iN+fPn\n4/V6OXTokEYKOZAjA0FXPfLII8Ff0ldffRWA8ePHx+zxpk2b4vbZHR1nzZo1eL1errvuOtatW4fX\n66WkpISlS5cGZ9UGLoiBYOX1emlsbKSpqQnwlc4Cmcrx48eD+wL27dsHtB1YAqWxXr16cerUqWBg\nOXbsGJmZmQwdOpQ9e/aQmZnJqFGj2LZtGxkZGVx99dVs376dW265BdM0ycvLo7CwkLy8vEj+14bV\nVh+BMoLYOXLkCIsXL2bcuHEYhqGRQg5kBMoVTnLTTTfZV1xxRXBbncXt27RpU5eD38aNGxkzZgxr\n166lpaWFkpISlixZQm1t7RkZU2gWFCh1gS/zCQSWSGVkZJCVlcU555zDiRMnyMnJ4fLLL2fHjh1k\nZWVx//3384Mf/IDc3NwOP2vWrFk8//zzZyx1sHLlSu677z7efvtthg4d2qU2is8999zDn/70J555\n5hl+/OMfJ7s5EoZpmh9YlnXW+jyODAS33367/cc//jHZzXClaALLmjVrsG2bsrIyXn/9dbxeL1dd\ndRVbtmyhpqYmWMoKDSwdLRedlZVFv379qK+vJzs7m9GjR7Njxw4AfvrTn/Iv//IvPProo7z++ut8\n9NFHwfetWbOGn/3sZ7z55psMHz48xmfJPRobGxk2bBinT5+mtLRUQ0ZTXJcDgWEYP8K3CmgffMtI\nG4ANbLBte2Uc2tohBYL0s2nTJmbPng20H0xKS0t5/fXXOX36NE1NTezfvx8IHyz69+/PRRddxP79\n+/nb3/4W3L9u3Truuusu/vznPzNixIh4/tPS2qpVq5g6dSqjRo1i5syZ6htIcREHAsMwLgMuBypt\n297TxvODgWJgt23bCR2Ko9KQtBc4ysrKWLduHbZts3fvXk6ePElhYSEDBgzgjTfeCL5/w4YN3HHH\nHbz22mtcdtllSfk3pIPvfve77N27lxdeeIFrrrkm2c2RDrQXCMJ1Fnts236uvSdt2/4M+MwfEBJK\nM4ulpKSEtWvXBrdDOygffPBBAH7yk5/w1ltv0b9//zOGjoImlMXCJ598EizDLViwQIHAwdq9Z7H/\nQg+AYRhLDMMY7X98YeBx69eJpJIxY8Zw+vRpTp06dcbQUTh7+GhtbS1O7C9LpkWLFpGVlcXVV1+t\nkUIO19mb1y+0bftD/7f/KfjuPyyS0gYNGgT4hsCGywj27NnDpZdeyogRI4L315Xw6uvrWb58OePH\nj2fp0qXqG3C4dgOBYRhvGIYx3zCMHwK7/f/9LlBh2/aMhLVQpIvOP/98wDdpr3VGELrExCuvvEJz\nczMej4cHH3yQyy67jKuvvppFixYpS2jH7373O06cOMHo0aM7frGkvHB9BI8AVUAZviygGN+tJvMN\nw6hMdAdxKC06J50RCARAuxnBnDlzOHToEKNHj2bPnj14PB6qq6s5cuQIjzzyCLNmzSI/P5+HH36Y\nCRMmJLT9qezll18G4K9//St33XVXklsj0Wo3ENi2vcX/cIX/BwiOJhoCJC0QqLNYOqNv377BZT5C\nZxUDnHPOORiGwdatWwFfQOjdu3dwjsTcuXPxeDycOnWKEydOMHXqVJ555hny8vKYOXOmq7947N27\nl5qaGoqKinQbyjQRrjTUq639tm1vsW17RUevE0k2wzCC/QStA0H//v2ZOXMmGRkZZGZm0qNHD0pK\nSli0aBFTpkxhwYIFlJaWMn36dPLz8wHYvXs3H3/8MbfeeiszZswIrmrqNqtWrQLglVdeUd9AmgjX\nWXyFv1+gXf7JZmeNSRVJFYHyUOvSEPiWRrjqqqtoaWnhuefOHCndVlCYOXMmPXv2xOv18oc//IFL\nLrmE0aNHnzGMNd3Zts3KlSu56qqrGDhwYLKbIzESrjT0pmEYvQ3DeAhfKah1r9lufKOJTsSzgSLR\nCGQErTuLA+6//35ycnLCDn8sKSkJfvMdOXIkFRUVHDx4kJ07d3L06FGmTJnChRdeyG9+8xvGjRsX\n+39ECtm6dSu7d+8mKysruPihOF/Y1Udt264BnkxQW0RiLlxGAGde5Dsj8PrAzOaGhgY+//xzqqqq\nuP3225k+fTr33HNPTNqeitavXw/4JpNpuen0EXYegWEYd3VUHkqGwKihwM8777yT7CZJigqskBq6\n8mgsBGY2v/nmmyxcuJCePXvS0tLCrFmzGDlyJGvWrInp8VLF+vXrufjiiyktLdUksjQSrrO4AhgK\nXG8Yxp8T16SOBUYNBX7cPIJDwtuyxTf47c0334zbMUpKSnj++ecZN24c5513HsePH+fee++lrKyM\nTZs2xe24ibZ37162b9/OpEmTWLRokbKBNBIuI/jAtu3ptm1PASalYmYg0pF7772X0tLSuA9zLCkp\nYfHixTz99NN861vfIi8vj+3bt3PbbbexYcOGuB47USoqKgAoLCxMcksk1sIFguOBB/6+gpr4N0ck\ntgKjfxL17bWkpISVK1fy3HPPkZeXR3NzM3fffTcTJkxwfHYQKHcFho9K+ggXCJ7wLzFxl3+Rud6B\nJzR3QCS8sWPH8t///d9cdtlleL1e3nvvPX760586NhicOnWKuro6zj//fPUNpKFwgWAOsBDfDWlm\n4gsMnxqGMR94Np6NMk2zyDTNZaZp3h3P44jEU0lJCa+99hpz584lKyuLuro6x86If++99zh9+jS/\n/e1v1TeQhsItQ/2sfxbxk7ZtT7Jteyi+yWMr8AWHDpmmWdxqe6JpmmWmaT7cibdPtixrYWeOI5LK\n/umf/okXXniBvn37smvXLsaOHeu4zGDjxo3k5OQwZsyYZDdF4qCzy1ADvr4C27YrgfKOXmuaZhmw\nLGS7GMCyrErA0zpItFINFPkDR34kbRRJRddeey0ffPAB/fr1Y/fu3dx5552OCgYbN27k8ssvp3v3\n7sluisRBm4HAMIzBoTefaa0zN6PxX/BDB2/fDHj8jwOrmgayhNCffMC0LGuz/3UqD0layMrK4qmn\nniIrK4uTJ0/y+OOPJ7tJnVJdXc3f//53Dh8+7KjgJZ3XZiDwX+gNwzAmh3QWRysf3zf9gEIAy7KW\nt/rxAFX+jMEElsfg2CIp4dprr+XZZ5+lR48ebNu2jRtvvDHlL67/+7//C8Bnn30WHEIq6aWjZai3\ngG/pacMwJuNbb6jKtu2/xLNRlmUFMonN8TyOSDJcf/31vP3221x99dVYlsUzzzyT0h2w7777LtnZ\n2Xz729/WiKE0FXatoYBWQWFwSFCoBio7ufCcB9+NbcCXHRwP89qwjh49imm2v+jp5MmTKS/vsBtD\nJGnOPfdcfvOb3/DII4+wZcsW/vKXv3Dttdcmu1lteu+997jiiitYvHhxspsicdKpQBDKXzZ6FsAw\njN74av0rO/HWJXy1ZHURUBnpsQP69esXXPxKxKluueUW/vCHP/Dxxx8zffp03nvvvWQ36Sy1tbVs\n3bqVadOmJbspEkcRjRpqzT+KqM0gYJrmRN9/zIkA/s7fwGgiT2C7K7TonKSLf/3Xf2XQoEEcPHiQ\n73//+ynXX/DBBx/g9Xo1bDTNGe3dnNt/05n25gsY/v8eby8QxNPtt99u//GPf0z0YUXiorGxkZEj\nR1JXV8e3v/1tli1b1vGbEmTOnDn853/+J9u3b2/3ng7iHKZpfmBZ1ll19XCdxSvaey7ZdPN6SSfd\nunXj17/+NT//+c/56KOPUuqGL5WVlXTv3p0tW7akTJsk9trNCFKZMgJJR6WlpezcuZMRI0bw5z8n\nf+V327a54IILaGlpobS0lEWLFiW7SRKliDOCAH+J6Dp8ZaJqfGUhG9iQjLIQKCOQ9PToo49y1113\n4fF4aGlpITMzM6ntOXDgAC0tLQwbNkzDRtNcu4HAMIzLgMvxDQ89q0zkH0b6I2C3bdsfxrGNZwnc\nmEYknYwbN46nn36a8vJyvve97/Hoo48mtRyzbds2AH73u99RXBxuRRhxunAZgce27efae9I/jPQz\nwzAGx75Z4SkjkHT1j//4j/To0YPt27czf/78pAaCrVu3kpGRwcUXX5y0NkhihOssDq4nZBjGEmC2\nbdsfGoZxIZAfyAI6s+5QrCkjkHSVkZHBbbfdxoIFC/jOd76T1LZs3bqVwYMHk5eXl9R2SPx1dh7B\nQn8QGAxMwbeAnIjEwUMPPUR+fj7vv/9+UtuxdetWvvnNbya1DZIY4foI3gB2AxuAzf57FhcAFbZt\n70lM89qm0pCks7y8PK655hpeffVVVq9ezU033ZTwNtTU1LB//35uu+22hB9bEi9cH8EjfLVc9BSg\nGF8gyDcMozLRHcShVBqSdHf06FHAN6ErGYFg6dKlgG8IqaS/cHco2+JfQmKFbdvTbdu+3rZtE3gT\nGJK4Joq4z4MPPkjfvn358ssvaWhoSPjxX375ZQDeeuuthB9bEq/dQNDeDer9AWJFR6+LJ601JOmu\npKSEefPmUVtbyw033JDwNYgCfQP3339/Qo8ryRGuNHSFYRi9w00a888j+BKI6/0JWlNpSNygpKSE\nvLw8duzYQUVFRcKHkg4aNIhx48Yl9JiSHOGGj75pGEZvwzAewlcKal0s3I1vNFFn7kUgIhEyDIPv\nf//7rFq1iltvvTWhx963bx/nn39+Qo8pyRN2iQnbtmuAJxPUFhFppby8nFWrVlFTU5PQ4x44cEAj\n8Vwk7DwC//2Kf5ioxojImUaMGMGgQYNYt25dVJ+zZMkSrrzySsaMGcPw4cMZM2ZMu/0ODQ0NHD58\nmK9//etRHVOcI9w8ggp8t5fMNwyj3Lbt7yWuWeFpHoG4hWEYjBw5knXr1rFhwwauu+66iN7v9Xq5\n5557ePXVVwHIysri9OnT1NbWMnv2bNauXXvWew4cOACg0pCLhCsNfWDbdvCWlIZh/DBZq422ps5i\ncZMDBw5g2zbPPPNMxIHgxRdfDAaB/Px87rvvPlavXs22bdvanSMQCATKCNwjXGkoeHN5f19BYouU\nIgLA1KlTARg5cmRE71u2bBn/9m//xpAhQ/jOd77DggULmDJlCuvXr+ehhx7i//7v/7jxxhvPKhHt\n27cP8I0aEncIFwieMAxjvr+fYDTQO/BEMuYOiLjVDTfcwIABAzh58mSn31NVVcX06dPxer3079+f\nl1566Yzhpz/5yU/IzMzEsiwqKirOeO/+/fvJyclhwIABMfs3SGoLFwjmAAvx3ZBmJr7A8KlhGPOB\nZxPROBHxueSSS4L3B+jInj17uOGGG2hubmbYsGFMmzbtrNf06tWL8ePHA3DjjTee8dz+/fsZOHAg\nGRmdXZNSnC7cEhPP+mcRP2nb9iTbtocCJrCC9m9qLyJx8M1vfpOdO3d2uNzEiRMnuOOOO6ivr6el\npYWBAwe2OxHt3//93+nWrdtZq5xu374dj8eT8NnMkjwd3qoylL+voNIwjN1xak+naNSQuM2IESNo\naWlh586dYfsK7rjjDnbt2sVtt93GgQMHwt5i8pxzzmHs2LEsXryYsrIyrr/+esDXWdzQ0JCU2cyS\nHBEFAgD/jWkeBn4W68Z0lkYNidsE1v7ZunVru4Fg7dq1vPvuu4DvYt6Zm80fPnwY27aZO3duMBBk\nZGQwaNAg3afYRTpdBDQMo5d/bkElMCl+TRKR1i644AJyc3N5+umn2yzZ1NXV8dBDD5Gbm8uoUaM6\nfRGfPHkyAFdccQUALS0t1NfXM2nSJGUDLtJhIAgJAB8Alr+vIOG3pxRxs4yMDHJycjh48OBZo3wA\n/va3v1FTU0NDQwMFBQWdvojfdNNNdOvWDcMwAIJLWfTu3Tvc2yTNdLTExG/xBYA3bNu+KORm9rpb\nhUiCjRo1im7durX5bf/9998nKyuLsWPHRlTSyczMpKioiE8//RT4KhDk5+fHptHiCB31EWwA3sdX\nDkoo0zTvxneHtGrLsjYn+vgiqeZb3/oWGzdu5PLLLz/ruffee49Ro0YFbygTiaFDh/Lxxx8Dygjc\nKmxGYNv2m/6b0FxhGMYPI51IZppmcavtiaZplpmm+XAH77sbqLQsq1JBQMRn8ODBwFczfwMaGhr4\n6KOPgnX+SF100UXs27ePhoYGZQQu1anOYn9AWIkvIPwIMDp6j2maZcCykO1iAMuyKgFP6yDRyuVA\nsT9wFHWmjSLp7sILLwTgs89bzUpdAAAMNUlEQVTO7KL76KOPaGpq6nIgGDp0KF6vl6qqKjweD6CM\nwG0inUfwJhDsWArHsqxK0zSrQnbdjK/UBL6STxmw2TTNia3eWolv1dNKoAAoBx6JpJ0i6SgQCPbs\n2XPG/sCEsGgCAcCuXbuCgUAZgbtEPI8AIPSexRHIB6pDtgsBLMta3vqFpmnOxhcoABZ04VgiaSc/\nP58+ffqclRFs2LCBvLw8tm/f3qUhn0VFvqR7zpw5XHnllYAyArcJdz+CH9HxUhLV8Via2rIsD3BW\ngBBxuwsvvPCsQPD3v/+dU6dOdXkmcF5eHrm5uezZs4eWlhZyc3Pp1q1brJosDhDunsVd+dYfjgdf\nqQd82cHxMK8N6+jRo5im2e7zkydPpry8vKsfL5KyBg8eHJw9DNDc3ExjYyNf//rXo5oJPGTIEPbv\n38/QoUNpbm6ORVPFQbpUGuqiJfgWrQMoIoohqf369WP9+vUxaZSIkwwePJhVq1bR0NBAbm4u+/fv\nx+v18sADD0Q1E/jSSy/l2LFj5Obmqn/AheK2zqy/E9gMdAYHhoH6RxN5ohkWGlh0LvDzzjvvxKbR\nIimusbER27ZZvXo1ALt3+9Z/HDJkSFSfO2jQIL744gu++OIL9Q+4UNwyAn8n8PJW+xbG63gibhD4\n0vPSSy/xz//8z1RV+QbmRRsIAvcn3rFjh9YYcqFEloZiRquPilvddNNNfPTRR8GbyuzevZs+ffrQ\np090twgJBIL6+nqVhlzIkYFA9yMQtyou9s3DHD58OOC7JWW02QCceX9ilYbcx5GBQBmBuFV2djYA\nTU1NgC8QjBs3LurPHTBgANnZ2TQ3NysQuJBuSiriIIFA0NzcTF1dHV988UVMMoLMzEwGDhwIaFax\nGzkyI1BpSNwqJycH8GUEgaUm1qxZw+jRo6Pu5B00aBB79uxRIHAhRwYClYbErUIzgsBKoVu3bo3J\n/YUDHcYqDbmPSkMiDhIaCAIzgIuLi2Nyf+GWlhbg7EXtJP05MiNQaUjcKrSzONBh/Pjjj3PppZdG\n/dnbtm0D4PXXX+fOO++M+vPEORwZCFQaErcKzQgCgSCwL1q/+tWvePrpp3nggQdi8nniHI4MBCJu\nFegsDg0EgX3RKikp0axil3JkIFBpSNwqtDQU6COIVSAQ93JkIFBpSNwqtDTU2NgIKBBI9DRqSMRB\nDMMIzgAOZASx6iMQ91IgEHGYQCAI9BHobmISLQUCEYfJyck5ozSkjECi5cg+AnUWi5tlZ2cHO4sN\nwyAry5F/xpJCHPkbpM5icbPQ0lBOTg6GYSS7SeJwKg2JOEygNBQIBCLRUiAQcZhAaaipqUn9AxIT\nCgQiDhM6fFQZgcSCAoGIw4SOGtLQUYkFR3YWa9SQuJlKQxJrjgwEGjUkbqbSkMSaSkMiDpOTkxPM\nCBQIJBYUCEQcpvU8ApFopWRpyDTNicDNQDWw27KsJ5LcJJGUoUAgsRbXjMA0zeJW2xNN0ywzTfPh\nDt662bKsHwPLgIVxa6CIAwVKQ83NzRo1JDERt0BgmmYZvgt5YLsYwLKsSsDTOkiEsiyrKuSxJ15t\nFHGiQEbQ2NioUUMSE3ErDVmWVWmaZlXIrpuBDf7HVUAZsNlfBgpVaVmWx7+/ChE5g0YNSawlso8g\nH1/NP6AQwLKs5e28vgiojHejRJxG8wgk1lKysxhAHcQibcvOzub06dM0NTWpj0BiIpHDRz1Agf9x\nPnA8gccWSRtafVRiLZEZwRLA9D+Oquxz9OhRTNNs9/nJkydTXl7e1Y8XSWmB0lBWVpZKQxITcQsE\n/s5e0zTNiZZlLbcsa7PpUwZ4LMva3NXPzs3NZcKECcFtrTUkbpKdnY3X6+XUqVMqDUlMxHPU0HJg\neat9MZkToLWGxM0C5aCWlhZlBBITKdtZHI5WHxU3C734q49AYsGRgUAZgbhZ6MVfgUBiwZGBQBmB\nuJkyAok1RwYCZQTiZqGBQH0EEgtahlrEYUKzAI0aklhwZEag0pC4mTICiTVHBgKVhsTN1EcgsabS\nkIjDaNSQxJojMwKVhsTNlBFIrDkyEKg0JG6mjEBiTaUhEYdRZ7HEmiMzApWGxM1CL/4aPiqx4MhA\noNKQuFloOUgZgcSCSkMiDqPOYok1BQIRh1EgkFhzZGlIfQTiZgoEEmuODATqIxA306ghiTWVhkQc\nRovOSawpEIg4jDICiTUFAhGHUSCQWHNkH4GImxmGQXZ2NhkZGRiGkezmSBpQIBBxoEAgEIkFRwYC\nDR8Vt8vJyVEgkJhxZCDQ8FFxu+zsbDIzM5PdDEkTjgwEIm6XnZ1NVpb+fCU29Jsk4kA5OTnKCCRm\nUjIQmKaZDxQBBYBlWZYnyU0SSSnKCCSW4trbZJpmcavtiaZplpmm+XBHbwVuBgIBQURCZGdna50h\niZm4BQLTNMuAZSHbxQCWZVUCntZBohULKAPmAFXxaqOIU+Xk5CgQSMzELbe0LKvSNM3Qi/jNwAb/\n4yp8F/rNpmlObPXWSmCSZVmX+4PF3cAT8WqniBO1tLSwf/9+Nm3aRElJSbKbIw6XyCJjPlAdsl0I\nYFnW8tYvNE2zyp9RVANnPS/idj179uT48eNUVFQoEEjUUnJGimVZlf6fzZZlxbQ0tGDBglh+XNrT\n+YpMos7XAw88QGlpKVOmTEnI8eJFv1+Ridf5SmQg8OAbBQS+7OB4Ao8d9OyzzybjsI6l8xWZRJ2v\nkpISFi1a5PhsQL9fkYnX+UpkaWgJvtFA4BsJVNnVDzp69Cimabb7/OTJkykvL+/qx4uIuErcAoG/\nE9g0TXOiZVnLLcvabPqUAR7LsjZ39bNzc3OZMGFCcFtrDYmIdF3cSkP+i3+f0M5gy7IW+mv/C6P5\n7MBaQ7/85S/bDALvvPNO2O146cpxOvuecK9r77m29ndmn85XZPtS+Xx19n2xOl9t7df5Cv9cKpyv\nlOws7khg9dFZs2bx+9///qwTtXHjxrDb8dKV43T2PeFe195zbe3vzD6dr8j2pfL56uz7YnW+2tqv\n8xX+uVQ4X46cox66+uisWbNUFhIRiYIjA0Ho/Qi2bNnCO++8o2AgItJFhm3byW5DxEzTPArs9W/2\nBmpavaT1vtDt84BDcWpaW22J1XvCva695zpzbtrap/MV2b5UPl+dfV+szldb+3W+wj+XyPN1gWVZ\n/VrvdGQgiCfTNBdYlqWxpx3wrxBrAsXAQq0QG55pmkX45s+UActjPVEyHfmXmCnwr08mrYT8DeYD\nldH8DTqyszicKFY8Peu9bhDlCrEWvomCrlkhNorzVYxvja1KoPX6Wmkrmr9HvlqB2JU6ce5m+INk\nJb412bosrQJBNCue+r+xefw/rhDN+Qr5lpYfzZwQJ4nyfC33f2MrwyXrZ0W5AjH4JqG6UifPXb5/\nnwcYEs3xHNlZ3J4oVzwtwrfIXZFpmkVuSN1jsELsQtM0N5um+bBlWWm/Qmw058uyLI//j7vSDb9b\nEP35SkQbU1Unz53HXx4C2B3N8dIqELSh0yueApX+k+raVJTIVoi1/Be2IlzyDbcNkZyvMuARoMo0\nzQ3t/A6mu0j+HsF3sSs0TdP1gYG2z90CfOcIIKpJuukeCCLi/2W7LtntcAK3lINiJaSWK53khiwz\nGv7MMibZZVr1EbQhJVY8dRCdr8jofEVG56vr4nru0j0QLOGrES1RrXjqEjpfkdH5iozOV9fF9dyl\nVSAIXfEUvipfxGLF03Sk8xUZna/I6Hx1XaLPnSaUiYi4XFplBCIiEjkFAhERl1MgEBFxOQUCERGX\nUyAQiTHDMFyzCJ+kBwUCkQgZhlFsGMaG9p7Dv0xJuNeJpBIFApHIhZvab9q2vbkTrxNJGQoEIpEr\n46uVIGPxOpGkUiAQidwVwGbDMMoMwwiuGe8vC1mtXlfgf11UNw4RiSetPioSuWKg2Lbt5YZhhC4N\nbNq2vbDV62b7H/84Ya0TiZAyApHIFeD7pl8U0h/QFo9t24G7kqlEJClLgUAkAv7yzxJ8qz+WB4aK\n+v9b1ep17/s3r8NXStKwUklJCgQikSkCKm3bDlz0A3e0K7Ntu7L16/yPd+MrJWkEkaQkrT4qEgOG\nYdzdqn9AxDGUEYhEyTCM1veTFXEUBQKR6JWhu22Jg6k0JCLicsoIRERcToFARMTlFAhERFxOgUBE\nxOUUCEREXE6BQETE5f4f/YvmAJXDSqkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc3b1480e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog( kh, np.abs( (PCLASS/(cosmo.h())**(-3) - pk_nonlin[0,:])/pk_nonlin[0,:]   ) )\n",
    "plt.xlabel(r'$k/h$')\n",
    "plt.ylabel(r'$|\\Delta P(k)| / P(k)$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}