compare cosmosis and pycamb phi-phi

compare.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ahead-enhancement",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "christian-consistency",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cosmosis.runtime\n",
    "import camb\n",
    "import fitsio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "accepting-reproduction",
   "metadata": {},
   "outputs": [],
   "source": [
    "cosmo = 'cosmological_parameters'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "martial-demographic",
   "metadata": {},
   "outputs": [],
   "source": [
    "ini = cosmosis.runtime.Inifile(\"./om_params.ini\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "everyday-university",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load planck data\n",
    "def load_data():\n",
    "    clik = \"../cosmosis-standard-library/likelihood/planck2018/baseline/plc_3.0/lensing/smicadx12_Dec5_ftl_mv2_ndclpp_p_teb_consext8_CMBmarged.clik_lensing/clik_lensing/\"\n",
    "    fn=f'{clik}/pp_hat'\n",
    "    f=fitsio.FITS(fn)\n",
    "    vals = f[0][:]\n",
    "    fn=f'{clik}/bins'\n",
    "    f=fitsio.FITS(fn)\n",
    "    d=f[0][:]\n",
    "    d2 = d.reshape((9, 2501))\n",
    "    bin_ell = np.arange(2501)\n",
    "    ell_mids = [where(b>0)[0].mean() for b in d2]\n",
    "    fn=f'{clik}/siginv'\n",
    "    f=fitsio.FITS(fn)\n",
    "    sigma = linalg.inv(f[0][:].reshape((9,9))).diagonal()**0.5\n",
    "    return ell_mids, vals, sigma\n",
    "\n",
    "planck_ell, planck_cl, planck_sigma = load_data()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "spread-membrane",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# set up consistency module\n",
    "consistency = cosmosis.runtime.Module.from_options('consistency', ini, root_directory=os.environ['COSMOSIS_SRC_DIR'])\n",
    "config_block = cosmosis.runtime.pipeline.config_to_block([\"consistency\"], ini)\n",
    "consistency.setup(config_block)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "hired-summary",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make initial data block\n",
    "block = cosmosis.datablock.DataBlock()\n",
    "H0 = 89.68525000000001\n",
    "h = H0 / 100\n",
    "block['cosmological_parameters','hubble'] = H0\n",
    "block['cosmological_parameters','ombh2'] = 0.03620155501741413\n",
    "block['cosmological_parameters','omch2'] = 0.15716533383444936\n",
    "# neutrino stuff\n",
    "block['cosmological_parameters','omega_nu'] = 0.01\n",
    "block['cosmological_parameters','massive_nu'] = 3\n",
    "block['cosmological_parameters','massless_nu'] = 0.046\n",
    "block['cosmological_parameters','omnuh2'] = block['cosmological_parameters','omega_nu'] * h**2\n",
    "# initial power\n",
    "block['cosmological_parameters','a_s'] = 2.597861e-09\n",
    "block['cosmological_parameters','n_s'] = 1.014815\n",
    "# other params\n",
    "block['cosmological_parameters','omega_k'] = 0.0\n",
    "block['cosmological_parameters','tau'] = 0.055\n",
    "block['cosmological_parameters','w'] = -1.0\n",
    "block['cosmological_parameters','yhe'] = 0.245341\n",
    "block['planck', 'a_planck'] = 1.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "stylish-permission",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency relation input parameters: omega_nu, ombh2, omch2, omnuh2, hubble, omega_k\n",
      "Calculated ommh2 = 0.20141 from omch2+ombh2+omnuh2\n",
      "Calculated baryon_fraction = 0.17974 from ombh2/ommh2\n",
      "Calculated h0 = 0.896853 from hubble/100\n",
      "Calculated K = -0 from -hubble*hubble*omega_k/299792.458/299792.458\n",
      "Calculated omega_m = 0.250403 from ommh2/h0/h0\n",
      "Calculated omega_b = 0.0450075 from ombh2/h0/h0\n",
      "Calculated omega_c = 0.195396 from omch2/h0/h0\n",
      "Calculated omega_lambda = 0.749597 from 1-omega_m-omega_k\n",
      "Calculated omlamh2 = 0.602934 from omega_lambda*h0*h0\n",
      "Model okay\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# fill in other parameters\n",
    "consistency.execute(block)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "opposite-threshold",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "**** WARNING: Parameter 'halofit_version' in the [camb_planck_lensing] section never used!\n",
      "\n",
      "**** WARNING: Parameter 'use_nonlinear_lensing' in the [camb_planck_lensing] section never used!\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Prep the CAMB module from cosmosis\n",
    "module = cosmosis.runtime.Module.from_options('camb_planck_lensing', ini, root_directory=os.environ['COSMOSIS_SRC_DIR'])\n",
    "config_block = cosmosis.runtime.pipeline.config_to_block([\"camb_planck_lensing\"], ini)\n",
    "module.setup(config_block)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "apart-folks",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run on cloned datablock to avoid messing with first and get ell, cl\n",
    "b = block.clone()\n",
    "module.execute(b)\n",
    "\n",
    "ell = b['cmb_cl', 'ell']\n",
    "\n",
    "# different scaling from what is saveed - this is accounted for in the likelihood\n",
    "cl = b['cmb_cl', 'PP'] * (ell * (ell + 1) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "finished-poison",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do the same in camb\n",
    "\n",
    "pars = camb.CAMBparams()\n",
    "\n",
    "# fill parameters from same data block as above\n",
    "pars.set_cosmology(\n",
    "    H0=block[cosmo, 'hubble'], \n",
    "    ombh2=block[cosmo, 'ombh2'],\n",
    "    omch2=block[cosmo, 'omch2'],\n",
    "    mnu=block[cosmo, 'omnuh2'] * 93.14,\n",
    "    omk=block[cosmo, 'omega_k'],\n",
    "    tau=block[cosmo, 'tau',],\n",
    "    num_massive_neutrinos=block[cosmo, 'massive_nu'],\n",
    "    standard_neutrino_neff=3.046,\n",
    "    YHe=block[cosmo, 'yhe']\n",
    ")\n",
    "\n",
    "pars.InitPower.set_params(As=block[cosmo, 'A_s'], ns=block[cosmo, 'n_s'])\n",
    "pars.NonLinearModel.set_params(halofit_version='takahashi')\n",
    "pars.set_for_lmax(2500, lens_potential_accuracy=1.1, max_eta_k=500000);\n",
    "pars.set_nonlinear_lensing(True)\n",
    "results_tak = camb.get_results(pars)\n",
    "\n",
    "cl_camb = results_tak.get_lens_potential_cls()[2:, 0]\n",
    "ell_camb = np.arange(cl_camb.size) + 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "norman-output",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.008055266653726856, 0.008043444067562501)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check omnuh2 close\n",
    "pars.omnuh2, block['cosmological_parameters', 'omnuh2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cosmetic-trustee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x133e51af0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxoklEQVR4nO3deVyU5f7/8dc1wAACoqIiQogaorIroJlLLplmuVtu5ZKaS3VaTrvfU9/T6VunziltO6aVtqplmmnlSUszl1RU3HBXTFxAAdlkn+v3B8rPfQBnuGH4PB+PeRyZuZc3MId399z3fV1Ka40QQghxIyajAwghhKj+pCyEEEJYJWUhhBDCKikLIYQQVklZCCGEsErKQgghhFXORgewNaXUvcC9Xl5eE1u1amV0HCGEqFG2bt16Vmvd6MrnlaPeZxETE6Pj4+ONjiGEEDWKUmqr1jrmyuflYyghhBBWSVkIIYSwSspCCCGEVQ53glsIYT9FRUUkJyeTn59vdBRxk9zc3AgICMDFxaVcy0tZCCHKLTk5GS8vL4KCglBKGR1HVJLWmrS0NJKTk2nevHm51pGPoYQQ5Zafn4+Pj48URQ2nlMLHx6dCR4hSFkKICpGicAwV/T1KWQghapTTp08zfPhwWrZsSfv27bn77rs5cOBAlWY4f/48o0aNIjw8nLCwMDp37kxOTg4Anp6eAFgsFh577DHCwsIIDw8nNjaWo0ePVkm+sWPHsmjRIptuU85ZCCFqDK01gwYNYsyYMSxYsACAHTt2kJKSQlWO2DBz5kx8fX3ZtWsXAPv377/qRPHChQs5efIkO3fuxGQykZycjIeHR5VltDU5shBC1BirV6/GxcWFyZMnlz0XGRlJly5d0Frz9NNPl/2X/MKFCwE4deoUXbt2JSoqirCwMH7//Xeg9Ajg6aefJjQ0lF69erF582buuOMOWrRowffffw+UnqMZN24c4eHhREdHs3r16rJt+vv7l2UICQnB1dX1sqynTp3Cz88Pk6n0z2xAQAD169e/6nvasmULnTp1IjIykri4OLKzs0lKSqJLly60a9eOdu3asWHDBgDWrFlDt27dGDBgAC1atOC5557jyy+/JC4ujvDwcA4fPly23VWrVhETE0OrVq1Yvnz5Tf/s5chCCFEp/7tsD4kns2y6zbZN6/LSvaHXfX337t20b9/+mq8tXryYhIQEduzYwdmzZ4mNjaVr16589dVX3HXXXbz44ouUlJRw/vx5AHJzc+nRowdvvvkmgwYNYvr06axcuZLExETGjBlD//79ef/991FKsWvXLvbt20fv3r05cOAA48ePp3fv3ixatIiePXsyZswYgoODL8tz33330blzZ37//Xd69uzJ6NGjiY6OvmyZwsJC7r//fhYuXEhsbCxZWVm4u7vTuHFjVq5ciZubGwcPHmTEiBFcHL5ox44d7N27lwYNGtCiRQsmTJjA5s2bmTlzJu+++y4zZswAICkpic2bN3P48GG6d+/OoUOHcHNzq+yvRo4shBCOYd26dYwYMQInJyd8fX3p1q0bW7ZsITY2lrlz5/Lyyy+za9cuvLy8ADCbzfTp0weA8PBwunXrhouLC+Hh4SQlJZVtc/To0QC0bt2aZs2aceDAAaKiojhy5AhPP/006enpxMbGsnfv3svyBAQEsH//fl577TVMJhM9e/bkl19+uWyZ/fv34+fnR2xsLAB169bF2dmZoqIiJk6cSHh4OMOGDSMxMbFsndjYWPz8/HB1daVly5b07t277Hu4mBtKy8pkMhEcHEyLFi3Yt2/fTf185chCCFEpNzoCsJfQ0NAKn7jt2rUra9eu5YcffmDs2LE8+eSTPPjgg7i4uJRdEWQymco+RjKZTBQXF1vdrqenJ4MHD2bw4MGYTCZ+/PFH2rRpc9kyrq6u9O3bl759++Lr68t3331Hz549rW777bffxtfXlx07dmCxWC47Irj0464b5b7yaqebvYpNykJYpbXm/Pnz5J47S15mKvnZaRTlpFGYl0thfh6WojxUSQGquKD0f0sKsGiwACgTSpkwmRQmkxMWJ1csznUwuXri7O6B2c0Tc526uNZtiLt3Y+o39MN8xWe/QlzUo0cPXnjhBWbPns2kSZMA2LlzJ5mZmXTp0oUPP/yQMWPGkJ6eztq1a3nzzTc5duwYAQEBTJw4kYKCArZt28aDDz5Yrv116dKFL7/8kh49enDgwAH+/PNPQkJCWL9+PW3btqV+/foUFhaSmJjIHXfccdm627Zto0mTJjRt2hSLxcLOnTuJiIi4bJmQkBBOnTpVdgSUnZ2Nu7s7mZmZBAQEYDKZ+PTTTykpKanwz+qbb75hzJgxHD16lCNHjhASElLhbVxKysLBjVsxDoC5feZe8/WCvBzO/HmAc6ePkHv2OCUZyTjnnsY9PwWPwjTqlGThpbPxUAWU5zoOi1YU4YwGTGhAY7r4UOUbDj8LD7JNdcly9iHHzY8irwCcGjTDo3FzfJq2pFFAS5xd65TvByAcilKKJUuW8Pjjj/PPf/4TNzc3goKCmDFjBp07d2bjxo1ERkailOKNN96gSZMmfPrpp7z55pu4uLjg6enJZ599Vu79TZ06lSlTphAeHo6zszPz5s3D1dWVw4cPM2XKFLTWWCwW+vXrx5AhQy5bNzU1taygAOLi4njkkUcuW8ZsNrNw4UIeffRR8vLycHd3Z9WqVUydOpUhQ4bw2Wef0adPn0pdRRUYGEhcXBxZWVnMmjXrps5XgMxn4fDGrRiHKoaXWkwl49guClIPYso4ikfucRoWnaQx6Zctb9GKdOVNulMjcl18KHarh3arD+6lD5OHD86eDXD1bICbhzceHnVwdauDk2sdnM3uuLi44uR09akwi0VTXGKhuCiP4vxc8nKzyM/J4nxuNoXnsyjOTaMk+wyWnLPo82k456fhUXCGBsUpNLKcxemKojmlGpPiGkRevWCcm7ShQVAk/sGRuHl42/XnWdvt3bv3qo9aRM11rd/n9eazkCMLB6ItJZxO2kvKgc0UnEjEJX0/L+YepJlOxSV+Mc0uLHeG+px1acqxeh047N0Mk08LPBsH0cAviEZ+zWhodqOhjbOZTAqzyQmziyfU8aRuA99yr1tSXMTpk0c5m3yY3JTDFKUl4ZZ5hAbnD9P21DbMp4shoXTZ46oppz3bUNwkmvrBHQgKuw23Ol42/m6EqH2kLGoobSnhxOHdnN73B4XHt+GdsYfAwkP4qTz8gBKtOGVqwiGnBmx0a0GbkP7UbRaOX/NQGtWtx1VzJlZjTs4uNAlsRZPAq2+6Ki4q5NjRvZw5spOCE7twT9tDs+wEGmf/Ageh+AcTR5ybkdqgPS4tutK8/Z00aNzUgO9CiJpNyqKGyEhP49iONeQd2YjXma0E5e8lgDwCgHztwp/mFiQ27IP2i8S7RQyBraII8PTiie+HkV2UzesdO9KqcZTR34bNObuYadYqkmatIi97/uzJY/y5Zx15R7fgdTaBiNRl1DmzCDbBUVMzzvjE4hrSk5Yd7sbTq54x4YWoQaQsqqm0E4f5c/svFCVtwCc9gaCSJKKUpkQrjjo3Z2eDPpgC2tEouAPNQqJoZb76CqKE1AT2Z+xHo5n480Tm9J5DlAMWxrU0bNqMhk2bAaMAKCrMZ/+OdaQn/orHyT8ITV2Gx5lFFP7uzC63CPKDehLYcSC+zcOMDS5ENSVlUU2cSTnJ0fgV6CNr8E/fTIA+hQ+Qo9056taGrbf0xiu4E4ERXbm1bgNuLcc241Pi0ZSeGC6yFBGfEl9ryuJKLmY3QmJ7QWwvAPLz89gdv5Kc3T/RJGUt4fvfhP1v8qfpFk40vYuGHe7j1tBYlEnuWxUCpCwMU5CXzaEtq8jau4qGqRtpWXyERkqTo905WCeSY/6j8AntTsuwDoSXcyarK8X4xqBQaDQuJhdifK+6wKHWcnNzJ6xzf+jcH4BjhxNJ3rQE76QVxB3/GKfkj/hzsT+n/O/Cr9NwAtvEgQzNLWoxKYsqdOLIXk5sXkKdY6todX4HoaqYQu3EIddQ4gMn4xPWm6CIzkS7mG2yv6jGUYTUDyk9Z9Hl9Vp7VFEezVq2pVnLtsCLnEtN5uBvC3A/uIyY43Nx+voTjjoFcSpoEMG9xtPIL9DouLWak5MT4eHhFBcX06ZNGz799FPq1Ln2fTdJSUls2LCBkSNHAjBv3jzi4+N57733qjJypaxZs4Z//etfNhkE0BakLOyoqKiQfVt+IXPHMvxT19JcH8cfOKYC2NJ4KO5t7iQk9k7aetnv3gBPsyeeZk8pigqo1ziA2GF/Bf7KmdPJHFz9BQ0Pf0unw29TfGgmO+rEQtRIQu+4H2dXd6Pj1jru7u4kJCQAMGrUKGbNmsWTTz55zWWTkpL46quvyspCVJ58IGtjeZlp7PjxI7a+NYTzrwYR/vNwOpxeQJ5rQzaHPE3y6PU0e2kPt0/7kHY9huJhx6IQ5TNuxbiyO92v1KhJAJ1GPEer6Vs4MfI34v1H0yTvIJEb/0LWa63YOPsxTibtr+LE4qIuXbpw6NAh/va3v5WNtgrw4osvMnPmTJ577jl+//13oqKiePvttwE4efIkffr0ITg4mGeeeaZsnfnz55dNZvTss8+WPe/p6cmLL75IZGQkHTt2JCUl5aocOTk5ZUOZR0RE8O233wIwZcoUYmJiCA0N5aWXXipbPigoiOeff56oqChiYmLYtm0bd911Fy1btmTWrFlly2VlZdGvXz9CQkKYPHkyFovFZj+7ipIjCxs4dzaFQ78vxPXA97Q+v41IVUI6Xuz37oK5bV9CbhtA27pXj2Mvahb/VlH4t3qX4qK32Pb7UlT8x8Sd+AzmfkaCx2243jaJ1rf3R5mcjI5aNX56Dk7vsu02m4RD39fLtWhxcTE//fQTffr0oW/fvgwePJjHH38ci8XCggUL2Lx5M5GRkZd9lDNv3jwSEhLYvn07rq6uhISE8Oijj+Lk5MSzzz7L1q1bqV+/Pr179+a7775j4MCB5Obm0rFjR1599VWeeeYZ5syZw/Tp0y/L8sorr+Dt7V02GVJGRgYAr776Kg0aNKCkpISePXteNj5UYGAgCQkJPPHEE4wdO5b169eTn59PWFhY2XwdmzdvJjExkWbNmtGnTx8WL17M0KFDbfKjrigpi0o6k3KK/b8twPPwckLztxOjSjhJYzb53o939EBax/QgrpInpkXVyinMIbsom4TUhHJ9XOfs4kK7HkOhx1BSjh/kyIr3aXViMT6/jOXP1beQEjaR8L4TcXOX8avsIS8vj6ioKKD0yOKhhx7CbDbj4+PD9u3bSUlJITo6Gh8fn2uu37NnT7y9S4/o27Zty7Fjx0hLS+OOO+6gUaPS21VHjRrF2rVrGThwIGazmXvuuQeA9u3bs3Llyqu2uWrVqrKZ+4CySY6+/vprZs+eTXFxMadOnSIxMbGsLPr3L724Ijw8nJycHLy8vPDy8sLV1ZVz584BpeNJtWjRAoARI0awbt06KYuaICs9lf1r5mPev5S2+Ql0ViWcUr5s9x9J/dj7uDWyM02r2aWW1xtAUJS62XtRfG8JxnfiDPLz/o+NK+bRaNdsYnf+jTM732Zb89GEDXicuvVsPXhKNVHOIwBbu/ScxaUmTJjAvHnzOH36NOPHj7/u+pcO8e3k5GR1OPJLhzIvz/IXHT16lH/9619s2bKF+vXrM3bsWPLz86/Kcekw4xe/vrgPWw8zfjOq11+2aig/J4Ody95jzz974T6zNbE7/0ajwmQSAkZxfOiP+P1tP3GT3iM4uqtck18DXetelMpwc6/DbYOm0nL6Nnb3/JQUt+Z0OvouphlhxM+eSkZqsi1ji2sYNGgQK1asYMuWLdx1110AeHl5kZ2dbXXduLg4fvvtN86ePUtJSQnz58+nW7du5d73nXfeyfvvv1/2dUZGBllZWXh4eODt7U1KSgo//fRThb+nzZs3c/ToUSwWCwsXLqRz584V3oatyJHFNRQXFZG4YRlF276i7bnfiFCFJOPL5iYjadTxfoIjb692RxCicmx9L4oymQjrMhC6DOTQjvWkr/w37U98ReH737CxyRBuHfQijZrcYpvw4jJms5nu3btTr149nJxKzxtFRETg5OREZGQkY8eOveYc2AB+fn68/vrrdO/eHa01/fr1Y8CAAeXe9/Tp05k2bRphYWE4OTnx0ksvMXjwYKKjo2ndujW33HILt99+e4W/p9jYWB555BEOHTpE9+7dGTRoUIW3YSsyRPkFWmsSd24hY/08WqX+RGPSycKDPQ164RH7AKFxPa859Lao+YZdHD/LTveiJO3fwZkfX6XduZ8pwMyupkNpPfhFvBv523xf9ladhyi3WCy0a9eOb7755qr5sMW1yRDlFZCacoL9K+fS6MgSQi2HKNYm9np24GTEcNp0u4/b3OQkpaOz970oQSGRBIV8zYlDOzn5/d+JOfkVBe8vYnPgA4QOnY6HXCl30xITE7nnnnsYNGiQFIWd1Moji4KCPHb9+g1q5wLCz/+BWZVwxLklma2GEtxzLJ4+MoS1sJ8j+7ZzdtnLxOWuIZ26HGg9lehBT+DqenMzmVWF6nxkISpOjiyuI2lvPKd/nUXImRXEkM1Z6rGj6XCadhtPi9YybpKoGi1aR9Oi9VL2xq/G8vP/0HHf6xx//TPSOj5H5J0PyIUSolpy+LLIP5/Drp8/xWvPF7QuSqSpdmaXV2fMMQ/QtvMAGjrLvRDCGG1iuqPbrWXnmq+pu+4fRG18jANbZ+Fy75s0DzfuqhchrsVhy6IgP5dN7z9EmzM/EUsufyp//rj1SVr1nkh7mSlNVBPKZCKix3CKugxm/eJ3ab13JvUX3cPmNf0JHvEG9Rs2MTqiEIADl4Vr+gGiU1PZWfcO3G97iLYd+xAoh/eimnJxMXP7/U9xLmMMm+a/QGzKN+S+t4rNbZ8gZvATmJxr7v9VL467JTeI1mwO+9cz182X3Ed2E/PUt4R2uls+BxY1Qr36Dblt6mxODF/JcfOtxCX+g6TXYzm64zejo1UbTk5OREVFERYWxrBhwzh//jxQOuCfLa1Zs6ZsmA/hwGXh0aAp9Rv5GR1DiEpp1iaG0OfW8Ef7t/AoPkezxQPY/J9J5GZlGB3NcBeH+9i9ezdms/myUVqF/ThsWQhR0ymTiY73PoTrX+LZ1HAQMae/JvutGLavWmB95WokpzCHU7mnSEhNsPm2Lw5Rftn+cnLo2bMn7dq1Izw8nKVLlwKlc1u0adOGiRMnEhoaSu/evcnLywPg0KFD9OrVi8jISNq1a8fhw4cv2+aWLVuIjo6+6vnaRMpCiGquXn0fbnt0Lvv7fUO+qQ7R6x5m+78Hkpl63OhoVl0cqPFEzgkm/jzRpoVxcYjy8PDwy553c3NjyZIlbNu2jdWrV/PUU09x8X6ygwcPMm3aNPbs2UO9evXK5p0YNWoU06ZNY8eOHWzYsAE/v///qcSGDRuYPHkyS5cupWXLljbLX9NIWQhRQ7SJu5Omz25hQ+DDhGb9jv6gIztXfGJ0rBuy1UCNl7o4RHlMTAyBgYE89NBDl72uteaFF14gIiKCXr16ceLEibIJi5o3b142vHn79u1JSkoiOzubEydOlI275ObmVjZN6969e5k0aRLLli0jMLB2T6dbcy+xEKIWMru60Wn8GxxMvJ+Sb6cQ8ccTbN/zPc0f/A/1quE5OlsP1AjXH6L8oi+//JIzZ86wdetWXFxcCAoKKhsa/MrhyS9+DHU9fn5+5Ofns337dpo2rd2X3MuRhRA1UHDb9rR4dh3rm00lNGstxe93ZPuq+UbHukpU4yhC6ofg7+lf4blCKiszM5PGjRvj4uLC6tWrOXbs2A2X9/LyIiAggO+++w6AgoKCsius6tWrxw8//MDzzz/PmjVr7Jy8eqsRZaGU6qKUmqWU+kgptcHoPEJUtWvNE242m7l93Gv8OfRHMp3qE71uMvEzR3A+K92glNfmafbEz8OvSooCSs8/xMfHEx4ezmeffUbr1q2trvP555/zzjvvEBERQadOnTh9+nTZa76+vixfvpxp06axadMme0av1uw+kKBS6hPgHiBVax12yfN9gJmAE/CR1trqtFtKqYGAr9b6Q2vLVnSIciGqM2s3thUW5LPl0+focOJTzpoaknvvLFq262nzHJUZSFBuyqu+KjKQYFUcWcwD+lwRxgl4H+gLtAVGKKXaKqXClVLLr3g0vmTVkcBXVZBZiGrF2uWnZlc3bp80g8S+31CiFc2WDmXzp89hKecUoPY0t89cKQoHYPey0FqvBa48Lo4DDmmtj2itC4EFwACt9S6t9T1XPFIBlFKBQKbW2vociUI4kIpcfhrRsRcef9nINq/uxB39D/ve6E5K8qHrLi9EeRl1zsIfuPQi8eQLz93IQ8AN//NEKTVJKRWvlIo/c+bMTUYUonqo6OWn9er7EPvkIjZFvUqzggO4fdSVnSs/r4qowoHViBPcAFrrl7TWNzy5rbWerbWO0VrHNGrUqKqiCWFXFy8/Bcp9+akymegw8BHSH1hFqpMfEesfYdsH4ygquPGlouXhqBOm1TYV/T0aVRYngEtnrQ+48JwQ4go3c/npLbeGc8vT6/i90XDapS7m2JtdOPPngUpncXNzIy0tTQqjhtNak5aWhptb+WdnrJJpVZVSQcDyi1dDKaWcgQNAT0pLYgswUmu9x1b7lKuhhLjcph8/pe2mZ7EoJ453e5uw7vdVeBtFRUUkJyeX3eQmai43NzcCAgJwcbl8AjjDplVVSs0H7gAaKqWSgZe01h8rpR4B/kvppbOf2LIohBBX63D3GJKC21M8/wHCfpvIH4fWEzvu3zhVYK4MFxcXmjdvbseUorqqkiMLI8iRhRDXlpebw845D9Ph3HL2uEYRMOFLvBsFGB1LVBNG3mchhKhG3D086fD4l2wMf4WW+Xso/KALSTK5krBCykKIWuq2IY9xZOBSirQzfouHsH3ZB0ZHEtWYw5WFUupepdTszMxMo6MIUe21jb4dlylr2O/aluitz7Nl1sOUFBcZHUtUQw5XFlrrZVrrSd7e3kZHEaJGaOTrT+u/rmJ9w2HEnl7A3jfvJDPttPUVRa3icGUhhKg4s9nM7Y98xKaIV2iVv4vc97pyfN8Wo2OJakTKQghRpsPgxzh8z9c46yJ8FvRj769fGh1JVBNSFkKIy7SJ7UnRQ79yzKkZIb9NY9v8v4ODXmIvyk/KQghxFf9bmhPw+C9s8ehKu/3/Ztv7D1JSVGh0LGEgKQshxDV5edWl/ZOLWdtkDO3Ofs++f/cm+9xZo2MJg0hZCCGuy9nZma6T32FD+CsE5+0k451unD62z+hYwgAOVxZyn4UQttdpyGPsu/NzvC0ZuM69kyPbfjE6kqhiDlcWcp+FEPYR0bkfacN/Ilt54L/0fvaunGd0JFGFHK4shBD206J1JK4P/8oB51a0Wf8XEr7+h9GRRBWRshBCVIhvk6YEPfFfNrl3ISrxTbbNnoK2lBgdS9iZlIUQosK8PL2IfmIJaxsMod3Jr9g5Y4hNpmwV1ZeUhRDiMuNWjGPcinFWlzObXejyyEf81uwxIrNWc/Ctu8jJTKuChMIIUhZCiEpTJhPdxr3CxsjXuTV/N2ff6U76qaNGxxJ2IGUhhLhMTmEOp3JPkZCaUO51bhs0hT3dP6ZhcSrFs3ty+uA2+wUUhpCyEEKUSUhNYH/Gfk7knGDizxMrVBjRdwzi2IBFoC14fNmPPxPkXgxHImUhhCgTnxKPpnTQwCJLEfEpFZvHPrRdZ7JG/Ug69Wj03XAOrltsj5jCAA5XFnIHtxCVF+Mbg0IB4GJyIcY3psLbuDW4Lc4T/sufpgCCVk5gz38/tnVMYQCHKwu5g1uIyotqHEVI/RD8Pf2Z03sOUY2jKrUd/4BAfKauZK9LG9pseIodi/9l26CiyjlcWQghbo6n2RM/D79KF8VFDRs2JOgvP7HVrQORO19h2+fPy7wYNZiUhRDCbup61SX8ye/Z4Hkn7Q5/wLbZk+Vu7xpKaQdt+piYGB0fX7GTc0II+yguLmbDB5Ppmv4NO3z6EjHlM5Sz2ehY4hqUUlu11ledrJIjCyGE3Tk7O9N52mx+9ZtIZNpP7H1nECWFMjxITSJlIYSoEiYnE90nvcnPQU/TOnM9h2b0ozgv2+hYopykLIQQVUYpRe+x0/m19cvcmruNozP7UJCbYXQsUQ5SFkKIKtdrxOP8FvFPgvL2kjzzLvIyZW7v6k7KQghhiB5DHuaP2BkEFBwm5d07yUk7aXQkcQNSFkIIw3S550G23j4L36JkMj7oTWbKn0ZHEtfhcGUhw30IUbN06j2M3d0/oX7xGXI/7E36iUNGRxLXcMP7LJRS7wLXXUBr/Zg9QtmC3GchRM2SsHElzVeMId9UB+dxy/AJbGN0pFrpevdZOFtZT/7aCiGqRNRtd7LDaT6BP4yiZO7dpI1Zik9QhNGxxAVyB7cQwq4uTtE6t8/cci2/Y9tGmi4djrNJYxm9BJ+W7e0ZT1yhUkcWSqll3PhjqP42yCaEEGUi293GDqfFNF48jDqfD+Ds6KU0vFUKw2jWPoaScYWFEFUuMrI9u9S3NPx2CHW+GEDqyCU0bhVrdKxa7YZlobX+7eK/lVLuQKDWer/dUwkhar3wiGh2m5agvxmM51cDSRmxGN+QDkbHqrXKdemsUupeIAFYceHrKKXU93bMJYRwEDmFOZzKPVWh+bwvCguL5Nz9S8jBDff5gzm97w/bBxTlUt77LF4G4oBzAFrrBKC5XRIJIRxGQmoC+zP2cyLnBBN/nlipwmjbNoLM+5eQjTt1Fgzm1N6Ntg8qrCpvWRRpra+8y80xL6MSQthMfEo8+sKfiiJLEfEplbtCsU2bCHKGf0c2dfBYOIRTiRtsGVOUQ3nLYo9SaiTgpJQKvnCznvy2hBA3FOMbg0IB4GJyIcb3qisyyy2kdRi5I5aWFsbXQ0nZK3+CqlJ5y+JRIBQoAL4CMoHH7ZRJCOEgohpHEVI/BH9Pf+b0nnPT83q3CgklZ/hSsvCgzsIhnNm33jZBhVXlKgut9Xmt9Yta69gLj+la6/yLr1840hBCiKt4mj3x8/C76aK4KKR1KFn3f8c57YX7gqGkHZBzGFXBVgMJ3m6j7dw0GUhQCMfXtk0oGfctIUN74vLVUNIPy2gN9uZwo85qrZdprSd5e3sbHUUIYUcRoaGkDf2WHO2K6YtBnEvaYXQkh+ZwZSGEqD2iwiM4PfBrCi1OWD7tT+bxvUZHcli2Kgtlo+0IIRzM3D5zyz2IYGW0i47h2D0LsFg0RXP7kX3ygN32VZvdsCyUUm5KqUbXeL6RUsrtkqdm2jyZEEKUU2xsR470/RJTSQHZs3vz+LcjjY7kcKwdWbwDdLnG852Bty9+obWeZ8NMQggBlA5vfnGIc2viOnbhQO/P8dAFPLlrE+fPyhSttmStLNprrRdf+aTWegnQ1T6RhBCicjre3oP/DehHfZ3DuVl9yc84aXQkh2GtLOrcxLpCCFHlMm4p4WW//ngXnSHtg74UZ58xOpJDsPYHP1UpFXflk0qpWEB+A0KIaik9qJh1se/hU3iC0+/1wXI+w+hINZ61snga+Fop9fKFm93uVUr9L/D1hdeEEKJauTgkum9cK34Of4tG+Ukcf/9edEHOVctW5JxIbXfDstBab6Z0aHIFjL3wUEAHrfUme4cTQoiKuHJI9MCuESwPfoWAnN0kfTAYiguMjlhjWT3voLVO1Vq/pLUecuHxN611alWEE0LUbhWdOOlaQ6IPHjWFJYHP0TxzE4c/HAmWEjsmdlxykloIUS1VZuKkaw2JrpRi0LhnWdJ4Gi3PrOLgJw+Blul4KkrKQghRLVVm4qTrDYnuZFLc8/A/WFZvNMHJSzj4xeNSGBXkbHQAIYS4lotHCRpdoYmTPM2eeJo9rxoS3cXJxJ1TZ7JiZjZ9Ds/jyJL6N745QFzG2nAfTkqph5VSryilbr/iten2jSaEqM1sPXESgJvZmdsfmcOv5u602Pk2nZJO3XzQWsLax1AfAt2ANOAdpdRbl7w22G6pboLMZyGE47D1xEkAXu6uhE/9gnVOccRmbEWdPFjuE+i1mbWyiNNaj9RazwA6AJ5KqcVKKVeq6UizMp+FEMKaRvU8yRn2P4xv0oStpnNMWDFeCsMKa2VhvvgPrXWx1noSkAD8CnjaMZcQQthVctEBikxgUYpiSyHr9y4xOlK1Zu0Ed7xSqo/WesXFJ7TWf1dKnQT+Y99oQghRceWdO+PSE+hOGqI2fkZB61G4+rayc8Kaydod3KMvLYpLnv9Ia+1iv1hCCGFfl55An3bLs7TOKyJrTn+KM+Wk97WU6z4LpZSTvYMIIURVu3gCfXyvB9jY4QM8itJJ+U9/dH6W0dGqHatloZTyApZWQRYhhLiMvadkvdQ9d/fnpzav45t3iOOzhkBxYZXst6awdp+FH7AKmF01cYQQwjiD7x/HIv9nCTy3mT/njgOLxehI1Ya1I4vfgde11t9XRRghhDCSUorB459hQd1xBJ5YzslFzxgdqdqwVhYZgH9VBBFCiOrA7Gyi7+Q3+M6lH00T53Bm1UyjI1UL1sriDqCvUmpaFWQRQohqwbuOmZjJH7JGxdJg3ctkJiwzOpLhlLYy8uKFK6E+1FpPqJpIthETE6Pj462PUimEENez6+hJ1Lx+tFQnUeNX4BYYbXQku1NKbdVaXzVqY3kmPyqpaUUhhBC2EN68KWn9PyNde5D36RBKziUbHckwMp+FEELcQLf24Wy5bRbOxedJ/XAgFGQbHckQ1i6dzVZKZV3jka2UkrtWhBC1wsA+vfk++P9odP4wyXNGQEmx0ZGqnLXhPry01nWv8fDSWtetqpBCCGG04SPH8WXDxwg4+zunFj5e62bak4+hhBCiHJxMisETp/O1eRB+Bz4n7dd3jI5UpaQshBCinLzcXOg46T1W0YH6v7/E+d3LjY5UZaQshBCiAgIbelJ35CfssQShvp1IScpeoyNVCSkLIYSooLhWARzuOZtsi5nMT4bA+XSjI9mdlIUQQlTCwG5xfN/6DTzyU0j5eDiUFBkdya4criyUUvcqpWZnZmYaHUUI4SDGrRjHuBXjrnp+7H3DmOvzBL5pm0j55gkDklUdhysLrfUyrfUkb29vo6MIIRycs5OJEROeZaHLQHz3fc65tR8aHcluHK4shBCiKnnXcSFmwjus1dF4/vo8BYfWGh3JLqQshBDiJrX09aZ40BySLL4Uzx+NTj9qdCSbk7IQQggb6BEVzNr271BUXMy5T4ZCQY7RkWxKykIIIWxk7L29mOP7P9TNPkz6/EkONSSIlIUQQtiIyaSYNG4Cc8wP0CDpB7JXzzA6ks1IWQghhA15u7vQffyr/Fd3oM7av1N0cI3RkWxCykIIIWwsxK8uuv/7HLb4UbhgDJw7bnSkmyZlIYQQdtCnfTCro96mpLiA9E/ug6I8oyPdFCkLIYSwIqcwh1O5p0hITajQehMG9uajRs/RICuR9K8frdEnvKUshBDiBhJSE9ifsZ8TOSeY+PPEChWGk0kxdtxUPna6jwYHvyFvw2z7BbUzKQshhLiB+JR4NKVHBEWWIuJT4iu0fgMPM9EPvs5qSxQuK19A//mHPWLanZSFEELcQIxvDAoFgIvJhRjfmApvo10zH072eIdkiw/nvxgF2adtHdPupCyEEOIGohpHEVI/BH9Pf+b0nkNU46hKbWdktwi+CHoNVZBN1uejoaTYtkHtTMpCCCGs8DR74ufhV+miAFBK8ZeR/XnLbRp1U7eQu+Il2wWsAlIWQghRRbzcXBg67gnmW3rhseU9ShJrzhzeUhZCCFGFWjepi9s9b7DT0pyibx+GGjJCrZSFEEJUsUFxLfmpzesUFFvI/nwkFOUbHckqKQshhDDAX4beydueT+GVkUjO0qeMjmOVlIUQQhjAzcWJceOn8BED8dz9BcXbvjI60g1JWQghhEGa+XjgP+gf/GFpg17+OKTsMTrSdUlZCCGEgfpG3sKa8NfJKHEvvWGvINvoSNckZSGEEAZ7fGAX3vR6BtfsY+R9O7VaDjgoZSGEEAZzc3Hi4QfHMMMyHPcD32PZ9KHRka4iZSGEENXArY09Cer/PCtL2qH/+yKc2Gp0pMsoXQ0Pd2whJiZGx8dXbHRIIYQw2vT5a5mybxw+Xu64TVsH7vWqdP9Kqa1a66tGS5QjCyGEqEaeH9yJ1zyexjnnJAWLp1Wb8xdSFkIIUY14uDoz7YGR/LtkOK4Hl2PZPMfoSICUhRBCVDtt/OrS9O6n+bUkCr3iBTiZYHQkxysLpdS9SqnZmZmZRkcRQohKG90xiJ+CX+KMxYv8+Q9CfpaheRyuLLTWy7TWk7y9vY2OIoQQlaaUYvqwLvzd9Smcs49T+N2jhp6/cLiyEEIIR+Ht7sKkB0Yzo3gY5n3foePnGpZFykIIIaqxqFvqUafnX/mtJALLT8/C6V2G5JCyEEKIam5yt2Dm+79IWokHhQseNGT8KCkLIYSo5kwmxf+OvIMXTH/B6VwSJcueqPLzF1IWQghRA/jWdWP4sJHMKBqM0+5vYPvnVbp/KQshhKgherX1JSv2MdaVhFLyw9OQklhl+5ayEEKIGuT5fmG8W+9ZMkrcKF44Bgpzq2S/UhZCCFGDuLk48croHjxVMg1T+kH0D1Uzf7eUhRBC1DCtfL3o1e9+3i0ehNoxHxLsP3+3lIUQQtRAozsEsid4Mn9Y2mJZ/iSk7rPr/qQshBCiBlJK8frQaF5xfYLMEjOWb8ZCUZ7d9idlIYQQNVQDDzMvDu/B4wWTMZ3ZCyuet9u+pCyEEKIG69SyIaFdBzOr+F7YOhf2LLHLfqQshBCihnvizlasbDKRnQRjWfooZCTZfB9SFkIIUcO5OJl4e0QsT+u/kFdkQX8zHkqKbLoPKQshhHAAgT51mDKwB38tmIA6uRV++btNty9lIYQQDmJgtD/ukYP5sqQnbHgHDq6y2balLIQQwoH8fWAYcz0ncUg1w7LkYcg6ZZPtSlkIIYQD8XR15p/DOzC14FGK8nJgySSwlNz0dqUshBDCwbRvVp+7e3RjeuGDcHQt/P7WTW9TykIIIRzQI91v5XDTAfxAZ/Sa/4NjG25qe1IWQgjhgJydTMwY3o6XLRM4bWqC/nYCnE+v9PakLIQQwkEF+tThmf4xTDw/FUt2CiydVunpWKUshBDCgQ1tH0BgWCdeKx4J+3+ETR9WajtSFkII4cCUUvzfoHCWuw1gg1MseuX/wMmECm9HykIIIRxcvTpm3ro/imnnJ5Bl8oZF46Agu0LbkLIQQohaoNOtDRnWJZIJOVPQGUmw/MkKnb+QshBCiFriqd6tyG0Sx38YBru+rtB0rFIWQghRS7g6OzFzeBTvFg0g0TUK/eNf4cyBcq0rZSGEELVIsK8Xz/cLZWzmRApwLT1/UY7pWKUshBCilnmgYzNCQ1rxaP4kSNkNP0+3uo6UhRBC1DJKKd4YGsk2cyyLXAfBlo8gcekN15GyEEKIWqiRlytvDI3g+cxBnPBoC0sfhYxj111eykIIIWqpnm18ub9jC4anT6TYUgLfPnTdZaUshBCiFnvx7raYG7bgb5ZJkLzlustJWQghRC3mbnZi5vBovs6LZX29e6+7nJSFEELUcmH+3jzWM5iHTg++7jJSFkIIIZh6R0tCbvG97utSFkIIIXB2MvHWfZHXfV3KQgghBAAtG3le9zUpCyGEEFZJWQghhLBKykIIIYRVUhZCCCGskrIQQghhlbPRAcpDKRUIvAOkAwe01q8bHEkIIWoVux9ZKKU+UUqlKqV2X/F8H6XUfqXUIaXUc1Y2Ew4s0lqPB6LtFlYIIcQ1VcWRxTzgPeCzi08opZyA94E7gWRgi1Lqe8AJeO2K9ccDfwCLlFLjgc+rILMQQohL2L0stNZrlVJBVzwdBxzSWh8BUEotAAZorV8D7rlyG0qpvwIvXdjWImDutfallJoETLrwZcGVRzMOwBvIdLB922K7ld1GRdaz9bLWlmkInC3n/moKef/adhv2ev8GX/NZrbXdH0AQsPuSr4cCH13y9QPAezdYPwxYBMwC/lXOfcZXxfdWlQ9gtqPt2xbbrew2KrKerZe1toy8f2vGvmvT+7dGnODWWu+mtGBqu2UOuG9bbLey26jIerZe1sjfpVHk/WvbbVTp+1ddaBK7uvAx1HKtddiFr28DXtZa33Xh6+cBdOnHULbaZ7zWOsZW2xOiKsn7V1Q3Rt1nsQUIVko1V0qZgeHA9zbex2wbb0+IqiTvX1Gt2P3IQik1H7iD0hN2KZSeqP5YKXU3MIPSK6A+0Vq/atcgQgghKq1KPoYSQghRs8lwH0IIIaySshBCCGFVrSgLpVQLpdTHF27oE6JGUUoNVErNUUotVEr1NjqPqJ1qbFlUZMwprfURrfVDxiQV4moVfP9+p7WeCEwG7jcirxA1tiwoHXOqz6VPXDLmVF+gLTBCKdW26qMJYdU8Kv7+nX7hdSGqXI0tC631WkqHLL9U2ZhTWutCYAEwoMrDCWFFRd6/qtQ/gZ+01tuqOqsQUIPL4jr8geOXfJ0M+CulfJRSs4Doi3eLC1ENXfP9CzwK9AKGKqUmGxFMiBoxNtTN0lqnUfp5rxA1jtb6HUon/xLCMI52ZHECuOWSrwMuPCdETSDvX1FtOVpZVMWYU0LYi7x/RbVVY8viwphTG4EQpVSyUuohrXUx8AjwX2Av8LXWeo+ROYW4Fnn/ippGxoYSQghhVY09shBCCFF1pCyEEEJYJWUhhBDCKikLIYQQVklZCCGEsErKQgghhFVSFkIYQCmVpJRqeOHfOUbnEcIaKQshhBBWSVkIYWdKqdFKqc1KqQSl1IcX5q0QokaRshDCjpRSbSid3e52rXUUUAKMMjSUEJVQK4YoF8JAPYH2wBalFIA7kGpoIiEqQcpCCPtSwKda68sm3VJKjTUmjhCVIx9DCWFfv1A6w11jAKVUA6VUM4MzCVFhcmQhhB1prROVUtOBn5VSJqAImGZwLCEqTIYoF0IIYZV8DCWEEMIqKQshhBBWSVkIIYSwSspCCCGEVVIWQgghrJKyEEIIYZWUhRBCCKukLIQQQlj1/wAIrWu62K5GPAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "loglog(ell, cl, label='CosmoSIS camb')\n",
    "loglog(ell_camb, cl_camb, label='Python camb')\n",
    "errorbar(planck_ell, planck_cl, planck_sigma, fmt='.', label='Planck')\n",
    "xlim(10, 800)\n",
    "ylim(1e-8, 2e-7)\n",
    "xlabel(\"ell\")\n",
    "ylabel(\"l^2 C_ell\")\n",
    "legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "textile-relief",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'cosmosis camb / python camb')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAu+0lEQVR4nO3dd3xW5f3/8dcnmwRCIIQZ9hCjICOAqyKOinuvKlXE0Tpbta38bLXVqvWr1VbrQgsCdbda0TrAgaggEmSDYGTIEsIIe2R8fn/cJxgR4Q7kzsl4Px+P87jvc51z7vtzxZgP13Wdc13m7oiIiEQrLuwARESkZlHiEBGRClHiEBGRClHiEBGRClHiEBGRCkkIO4Cq0KRJE2/Xrl3YYYiI1ChTp05d4+5Zu5fXicTRrl078vLywg5DRKRGMbMleypXV5WIiFSIEoeIiFSIEoeIiFSIEoeIiFSIEoeIiFSIEoeIiFSIEoeIiFSIEsd+yF+9mY8WFIQdhohIKGKaOMxsuJmtNrPZP3LczOwRM8s3s5lm1qvcscvM7Ktgu6xceW8zmxVc84iZWSzrsCcnPPQRlw3/nNJSrWUiInVPrFsczwID93L8ZKBzsF0NPAFgZo2BO4F+QF/gTjNrFFzzBHBVuev29vkxdd6TE9leVBLW14uIhCKmU464+wQza7eXU84ERnlkGcLPzCzDzFoAxwLj3H0dgJmNAwaa2Xgg3d0/C8pHAWcBb8ci/rkrNrJpexEJ8cb2olK27Sxh046iXce/+KaQZz5eSJvMNFplpJDToiH1kuJjEYqISLUR9lxVrYCl5faXBWV7K1+2h/IfMLOribRiaNOmzX4F93/vfsn4+Xsfy3hw7ILv7ScnxJGRmkjDeomkJMaTFB9HYnwciQlxJMXHkZRgeyiLK1dmPyhLSYwnLTmetOQEUpPiqZ+cQFpyAukpiSQlaJhKRKpW2IkjZtx9GDAMIDc3d78GI4aefDBX/aQDxaVOSkIc9ZLiqZcYT8N6iaTXS2TT9mKe/Ohrju/alA3bili0dgsbthZRuLWIwm072VFcys7iUopKStm6rWTX+7LXopJSdux675Tsx5hJWTxlW3q59xmpiTRtkEyLjHq0aJhC84YppKck7s+PQkRkl7ATx3Kgdbn97KBsOZHuqvLl44Py7D2cHxMHNW/AQTT40eMpifH84bScSvu+klLfLZlEksy2ohK27Chh685ituwoZsuOEjbvKGbjtiI27LYtW7+VuSsi77fs/OH4S/3kBJo3TIkkkvSUXUmldaNUOjZNo3l6CiHcbyAiNUjYiWMMcL2ZvUhkIHyDu680s3eBe8sNiP8UGOru68xso5kdDkwGfg48GkrkMRAfZ8THxZOSWDnjJDuLS1m9aTsrN0S2bzdsi7wv3M7KjdtZsKqA1Zt24OUaOmlJ8XTIqk/HrDQ6ZtWnc7P6dG2eTpvGqcTFKaGISIwTh5m9QKTl0MTMlhG5UyoRwN2fBN4CTgHyga3A4ODYOjO7G5gSfNRdZQPlwLVE7taqR2RQPCYD47VBUkIc2Y1SyW6U+qPnFJWUsnrTDpas3cLXBVv4evVmvi7YzJTF6/nv9BW7zktNiueg5g04uEV6ZGvegIOaN6CBur5E6hxzr/3PIuTm5roWcqq4rTuLyV+9mXkrNzJv5abgdSMbtxfvOqddZirdszM4rHUGh2U35JCWurNMpLYws6nunrt7edhdVVKNpSYl0D07g+7ZGbvK3J2VG7bvSiKzl29kyuJ1jJkRaZ3ExxldmjWgZ5sM+rVvzOEdMmmWnhJSDUQkFtTikEqxeuN2ZizbwIylhcxYVsi0bwrZvCPSMmmbmUq/9o3p1z6TIzpm0jKjXsjRikg0fqzFocQhMVFcUsq8lZuYvGgtny1cx5TF69iwLfLwZMesNI7pksUxXbI4vH2murZEqiklDiWOUJWWOvNXbeLT/DVM+GoNkxeuZUdxKUkJcfRt15hjujShf5emdGlWX7cDi1QTShxKHNXK9qISPl+0jgkLCpjwVQELVm0GILtRPU7MacaJOc3o264xCfF6Ml4kLEocShzV2soN2xg/v4D35q7i4/w17CwupWG9RI7r2pQTc5pxTJcs6ifrXg6RqqTEocRRY2zZUczHXxUwdu4qPvhyNYVbi0iKj+PITpkMPKQ5Aw9tTkZqUthhitR6ShxKHDVScUkpeUvWM27uKsbNXcU367aSEGcc0yWL0w9rwYk5zdUSEYkRJQ4ljhrP3ZmzYiNvzFjBGzNWsGLDdpIT4jj+4Kac3r0lA7o2rbTpWkREiUOJo5YpLXW++GY9b8xYwf9mrWTN5p2kpyRwVs9WnN+7NYe2StfdWSIHSIlDiaPWKi4pZdLCtfx76jLenv0tO4tL6dq8ARfktuasnq1onKbxEJH9ocShxFEnbNhaxJiZK3glbykzl20gMd44MacZl/ZryxEdM9UKEakAJQ4ljjpn3sqNvJK3jNemLWP91iI6Na3Pz49oy9k9W2lWX5EoKHEocdRZ24tKeHPmSkZNWszMZRtIS4rnnF7Z/PyItnRu9uMLdYnUdUocShwCTF9ayKhJi3lzxkp2lpRyRIdMBh/VjhMObqaFqkR2o8ShxCHlrN28g5fylvLcZ9+wvHAbHZqkMeQn7Tm3V7Zu6RUJKHEoccgeFJeU8vbsb3n644XMXLaBxmlJDDq8LYOOaEuT+slhhycSqv1OHGaWCfwROApw4BMiS7mujUGcMaHEIfvi7ny+aB1Pf7yQ9+atJjkhjnN6ZfOL/h1om5kWdngioTiQFQBfBCYA5wb7lwAvASdUXngi4TIz+nXIpF+HTPJXb+afnyzkP18s4+W8pZx5WEuuHdCJTk3rhx2mSLUQTYtjtrsfulvZLHfvFtPIKpFaHLI/Vm/czrAJC3lu8jdsLy7h1G4tuOG4zhzUXHdiSd3wYy2OaBY7GGtmF5lZXLBdALxb+SGKVC9N01P4/Wk5fPK7Afyif0c+/HI1J/1tAteMzmP28g1hhycSmh9tcZjZJiJjGgakAaXBoThgs7unV0mElUAtDqkMhVt3MvzTxYz4dBGbthdzXNem3HBcJ3q2aRR2aCIxobuqlDikkmzcXsSoiYt55pNFFG4tYsBBWdx60kEc0rJh2KGJVKoDShxm1h1oR7nBdHd/tTIDjCUlDomFLTuKGTlpMU99tJAN24o4rXsLbj6xCx2yNIgutcOB3I47HOgOzOG77ip39ysqPcoYUeKQWNqwrYinJyxk+KeL2FFcyvm9s7nx+M60zKgXdmgiB+RAEsdcd8+JWWRVQIlDqkLBph08Pj6f5z77BgwGHd6W6wd0opGmdZca6kDuqppkZjU6cYhUhawGydx5+iF8cGt/zurRkhGfLqL/Ax/y9ISF7CguCTs8kUoTTYujPzAG+BbYQeQuK3f37rEPr3KoxSFhmP/tJu59ax4fLSigTeNUfjewK6d0a641QaTGOJCuqnzgZmAW341x4O5LKjvIWFHikDBNWFDAPf+bx/xVm+jdthG3n3owvXQLr9QABzLlSIG7j4lBTCJ1wjFdsjiqUxNezlvKX8cu4JzHJ3LGYS0ZekpXWjTUALrUPNGMcUwzs+fN7GIzO6dsi+bDzWygmc03s3wzu20Px9ua2ftmNtPMxptZdrlj95vZ7GC7sFz58Wb2hZlNN7NPzKxTVDUVCVF8nHFx3zaM/82x3HBcJ96Z8y3H//Ujnhj/tcY/pMaJpqtqxB6K93k7rpnFAwuAE4FlwBTgYnefW+6cV4A33X2kmR0HDHb3QWZ2KvAr4GQgGRgPHO/uG81sAXCmu88zs2uBvu5++d5iUVeVVDffrN3K3f+by7i5q+jQJI07Ts/h2IOahh2WyPfsd1eVuw/ez+/sC+S7+8IggBeBM4G55c7JITJ+AvAh8N9y5RPcvRgoNrOZwEDgZSLToJRNd9IQWLGf8YmEpk1mKk//PJcP56/mrjfmcvmIKZyY04w7TsuhdePUsMMT2at9dlWZWYqZXWdmj5vZ8LItis9uBSwtt78sKCtvBlDW7XU20CBY/2MGMNDMUs2sCTAAaB2cdyXwlpktAwYBf/mRuK82szwzyysoKIgiXJGqN+Cgprzzq5/wu4Fd+TR/Dcc/9BF/f+8rdV9JtRbNGMdooDlwEvARkA1sqqTvvxXob2bTgP7AcqDE3ccCbwETgReASUDZ/0m/Bk5x92xgBPDQnj7Y3Ye5e66752ZlZVVSuCKVLzkhnl8e25H3b+nPiTnNePi9BZz6yCdMWbwu7NBE9iiaxNHJ3f8AbHH3kcCpQL8orlvOd60EiCSc5eVPcPcV7n6Ou/cEbg/KCoPXe9y9h7ufSOTZkQVmlgUc5u6Tg494CTgyilhEqr0WDevx2M96MeLyPmzbWcL5T05i6Kuz2LCtKOzQRL4nmsRR9ltbaGaHEhlXiGYUbwrQ2czam1kScBGRBwl3MbMmZlYWw1BgeFAeH3RZlU2w2B0YC6wHGppZl+CaE4F5UcQiUmMM6NqUcTcfw1U/ac9LU77hhIc+4s2ZK6gLM1lLzRBN4hhmZo2A3xP5wz8X+L99XRQMbF9PZNGnecDL7j7HzO4yszOC044F5gd3SjUD7gnKE4GPzWwuMAy41N2Lg8+8CviPmc0gMsbxm+iqKlJzpCYlcPupOYy5/miap6dw/fPTuOLZKSxbvzXs0ES0HodIdVdcUsqzExfz0LgFGDD0lIO5pF8bTV0iMbffkxya2b1mllFuv5GZ/bmS4xORH5EQH8eVP+nA2F8fQ882jfj9f2dzyTOTWbpOrQ8JRzRdVSeXDVgDuPt64JSYRSQie5TdKJXRQ/py79ndmLG0kIF/m8C/PluisQ+pctEkjngzSy7bMbN6RJ7mFpEqZmb8rF8b3lXrQ0IUTeJ4DnjfzIaY2RBgHDAytmGJyN6UtT7uOfvQXa2PV/KWqvUhVSLaNccHAicEu+Pc/d2YRlXJNDgutdmy9Vu55eUZTF60jpMPbc69Z3fTqoNSKfZ7PY7aQIlDaruSUufpjxfy17HzaZSaxIPnH8YxXTRjghyYA1k6VkSqufg44xf9O/Lf646iYb1Efj78c/44Zg7bizTnlVQ+JQ6RWuSQlg1544ajufzIdjw7cTGnP/oJc1ZsCDssqWWUOERqmZTEeP54xiGMuqIvG7YVcfbjExk9abEGzqXSRPMA4FFmNs7MFpjZQjNbZGYLqyI4Edl/x3TJ4u2bfsKRHTP5w+tzuO75LzRholSKaFoc/yQydfnRQB8gN3gVkWous34ywy/rw9CTuzJ2zipOe/RjZiwtDDssqeGiSRwb3P1td1/t7mvLtphHJiKVIi7OuKZ/R1665ghKS+G8JyfyzMcL1XUl+y2axPGhmT1gZkeYWa+yLeaRiUil6t22Ef+78WgGHNSUP/9vHleOzGP9lp1hhyU10D6f4zCzD/dQ7O5+XGxCqnx6jkPkO+7OsxMXc+9b88iqn8zjl/amR+uMsMOSakgPACpxiHzPzGWF/PJfX1CwaQd/OvMQLu7bJuyQpJo5kGnVG5rZQ2aWF2x/NbOGsQlTRKpK9+wM3rzhaPp1aMzQV2fxu3/P1AODEpVoxjiGA5uAC4JtIzAilkGJSNVolJbEs4P7ct2AjryUt5QLnprE8sJtYYcl1Vw0iaOju9/p7guD7U9Ah1gHJiJVIz7O+M1JXRk2qDeLCrZw+qOf8Gn+mrDDkmosmsSxzcyOLtsxs6MA/ZNEpJb56SHNef36o8hMS2LQPyfz1Edf65Zd2aNoEscvgMfMbLGZLQH+EZSJSC3TIas+/73uKE7u1oL73v6Sm1+eoXEP+YGEfZ3g7jOAw8wsPdjfGPOoRCQ0ackJ/OPinhzcvAEPjl3AwjVbGDaoN83SU8IOTaqJfSaOYNnYc4F2QIKZAeDud8U0MhEJjZlx/XGd6dysAb9+aTpn/OMThg3K5TA97yFE11X1OnAmUAxsKbeJSC130iHNefXaI0mMj+P8pybx+vTlYYck1cA+WxxAtrsPjHkkIlItdW2ezuvXHcUvn/uCm16czterN/PrE7tQ1vsgdU80LY6JZtYt5pGISLWVWT+Zfw3pxwW52TzyQT6/emk6O4o1aF5X/WiLw8xmAR6cMzhYg2MHYETmqupeNSGKSHWQlBDH/ed2p21mGg+8O58Vhdt4alAujdOSwg5NqtjeuqpOq7IoRKRGMDOuG9CJNo1TueWVGZzz+KeMGNyX9k3Swg5NqtCPdlW5+xJ3XwL8uex9+bKqC1FEqpvTD2vJC1f1Y+P2Ys5+/FM+X7Qu7JCkCkUzxnFI+R0ziwd6xyYcEakperdtzGvXHknjtCQufWYyb8xYEXZIUkV+NHGY2VAz2wR0N7ONwbYJWE3kFt19MrOBZjbfzPLN7LY9HG9rZu+b2UwzG29m2eWO3W9ms4PtwnLlZmb3BGugzzOzGytUYxGpNG0z03j1l0fSo3UGN744jRGfLgo7JKkCe+uqus/dGwAPuHt6sDVw90x3H7qvDw5aJo8BJwM5wMVmlrPbaQ8Co4KB9ruA+4JrTwV6AT2AfsCtZU+uA5cDrYGu7n4w8GLUtRWRSpeRmsSoIX35aU4z/vTGXO5/50vNcVXLRdNVdZCZnWJm0ZxbXl8gP5hRdyeRP/Bn7nZODvBB8P7DcsdzgAnuXuzuW4CZQNmzJL8E7nL3UgB3X13BuESkkqUkxvP4Jb35Wb82PDH+a259ZSZFJaVhhyUxEk0yeBy4BPjKzP5iZgdF+dmtgKXl9pcFZeXNAM4J3p8NNDCzzKB8oJmlmlkTYACRVgZAR+DCYFGpt82s856+3MyuLlt8qqCgIMqQRWR/xccZ95x1KL8+oQv/+WIZV43KY+vO4rDDkhjYZ+Jw9/fc/RIiXUeLgffMbKKZDTazxAP8/luB/mY2DegPLAdK3H0s8BYwEXgBmASUPW2UDGwPljN8mshCU3uKe5i757p7blZW1gGGKSLRMDNuOqEz957djQkLCrj46cms27Iz7LCkkkXV/RS0Ai4HrgSmAX8nkkjG7eWy5XzXSgDIDsp2cfcV7n6Ou/cEbg/KCoPXe9y9h7ufSOShwwXBZcuAV4P3rwF6EFGkmvlZvzY8cWlvvly5kfOemMgKrSpYq0Sz5vhrwMdAKnC6u5/h7i+5+w1A/b1cOgXobGbtzSwJuAgYs9tnNyk3djKUoPVgZvFBssLMuhNJDmOD8/5LpOsKIq2UBYhItXPSIc3515X9KNi0g/OfnMTiNZobtbaIpsXxiLvnBHdZrSx/IOgu2iN3LwauB94F5gEvu/scM7vLzM4ITjsWmG9mC4BmwD1BeSLwsZnNBYYBlwafB/AX4NxgSpT7iLSCRKQa6tOuMc9fdThbdxZzwVOT+GrVprBDkkpg+7ptzsxSgGuBo4nMXfUJ8IS7b499eJUjNzfX8/Lywg5DpM5asGoTlzwzmeKSUkYP6cehrRqGHZJEwcym7qmBEE2LYxSRp8cfJbJsbA4wunLDE5HarEuzBrxyzRGkJiVw8bDPmLpEU5TUZNEkjkPdfYi7fxhsV7HbNCQiIvvSrkkar/ziCJo0SObSZz7n0/w1YYck+ymaxPGFmR1etmNm/QD1+4hIhbXMqMdL1xxOm8apDH52Cu/NXRV2SLIfokkcvYks5rTYzBYTeaaij5nNMrOZMY1ORGqdpg1SeOmawzm4eQN+8a+pvDVr5b4vkmolmqVjtWysiFSqjNQk/nVlPwaPmMINL0zDHU7t3iLssCRK+0wcwfobIiKVqkFKIs9e0ZfBIz7nxhen4TindW8ZdlgShYpOXCgiUmnqJycwYnBferdpxE0vTteaHjWEEoeIhCqSPPrQu20jfvWSkkdNsLeFnN41s1+bWdeqDEhE6p605ARGXB5JHje9OI0xSh7V2t5aHJcB64E/mtkXZvaEmZ1pZlqVXkQqXVpyAs8O7kOfdo351YvT1PKoxva2AuC37v6su18E5BJ5grw3MNbM3jOz31ZVkCJSN6QmRbqtcts15lcvTWfsnG/DDkn2IKoxDncvdfdJ7n6Hux9FZKbb5fu6TkSkolKTEhh+eR+6tWrI9c9PY/x8LfJZ3ezX4Li7r3H35yo7GBERiAyYj7yiL52b1eea0VOZ9PXasEOScnRXlYhUSw3rJTJ6SD/aNE5lyMgpTF2yPuyQJKDEISLVVuO0JJ67sh/N0lO4fPjnzFq2IeyQhOhWAMw0s0eDO6ummtnfy1bnExGJtabpKTx3ZT8apiYyaPhk5n+rxaDCFk2L40VgNXAucB5QALwUy6BERMprmVGP5688nJSEeC555jO+Ltgcdkh1WjSJo4W73+3ui4Ltz0SWeRURqTJtMlN57qp+AFz6zGSWF24LOaK6K5rEMdbMLjKzuGC7gMg64iIiVapjVn1GD+nH5h3FDPrnZNZu3hF2SHXS3qYc2WRmG4GrgOeBncH2InB11YQnIvJ9B7dIZ/jlfVhRuI3LRnzOpu1FYYdU5+ztyfEG7p4evMa5e0Kwxbl7elUGKSJSXp92jXnikt58uXITV43KY3tRSdgh1SlR3Y5rZt3N7AwzO6dsi3VgIiJ7M6BrU/56wWFMXrSOG16YRnFJadgh1Rn7XMjJzIYD3YE5QNl/GQdejWFcIiL7dGaPVhRuLeLOMXO47dVZ/N+53YmLs7DDqvWiWTr2cHfPiXkkIiL74bIj27F+607+9t5XZNRL5PZTD8ZMySOWokkck8wsx93nxjwaEZH9cNPxnSncWsQznywis34yvzy2Y9gh1WrRJI5RRJLHt8AOwAB39+4xjUxEJEpmxh2n5bBuy07uf+dLmqUnc06v7LDDqrWiSRz/BAYBs/hujENEpFqJizMeOL87azbv4Lf/nklWg2R+0jkr7LBqpWjuqipw9zHBU+NLyraYRyYiUkHJCfE8Oag3nZrW5xejpzJ7uSZFjIVoEsc0M3vezC7W7bgiUt2lpyTy7OC+NKyXyOBnp7B03dawQ6p1okkc9YiMbfwUOD3YTotlUCIiB6J5wxRGXtGXHUUlXDbic9Zv2Rl2SLXKPhOHuw/ew3ZFNB9uZgPNbL6Z5ZvZbXs43tbM3jezmWY23syyyx2738xmB9uFe7j2ETPTFJkiskedmzXgmcv6sGz9NoaMnKKnyytRNOtxpJjZdWb2uJkNL9uiuC4eeAw4GcgBLjaz3Z8HeRAYFdyhdRdwX3DtqUAvoAfQD7jVzHZNc2JmuUCjaCooInVX3/aN+fuFPZi2tJAbXphGSamHHVKtEE1X1WigOXAS8BGQDUSzkkpfIN/dF7p72eSIZ+52Tg7wQfD+w3LHc4AJ7l7s7luAmcBA2JWQHgB+G0UMIlLHndytBXeelsO4uau4c8xs3JU8DlQ0iaOTu/8B2OLuI4FTibQC9qUVsLTc/rKgrLwZQNlA+9lAg2B1wRnAQDNLNbMmwACgdXDe9cAYd1+5ty83s6vNLM/M8goKCqIIV0Rqq8uPas81/Tvwr8++4bEP88MOp8aLJnGUzVlcaGaHAg2BppX0/bcC/c1sGtAfWA6UuPtY4C1gIvACMAkoMbOWwPnAo/v6YHcf5u657p6blaV7uUXqut+d1JWzerTkwbELeH368rDDqdGieQBwmJk1An4PjAHqA3dEcd1yvmslQKSL63v/tdx9BUGLw8zqA+e6e2Fw7B7gnuDY88ACoCfQCcgP5qJJNbN8d+8URTwiUofFxRn3n9edFRu285tXZtIyox592jUOO6wayWLV32dmCUT+2B9PJGFMAX7m7nPKndMEWOfupWZ2D5HWxh3BOEaGu681s+5EFpLq4e7Fu33HZnevv69YcnNzPS8vr/IqJyI11votOznniYkUbt3Ja9ceRbsmaWGHVG2Z2VR3z929PJq7qu41s4xy+43M7M/7ui74I389kWVm5wEvu/scM7vLzM4ITjsWmG9mC4isY35PUJ4IfGxmc4FhwKW7Jw0Rkf3RKC2JEZf3AeCKZ6dQuFXPeFTUPlscZjbN3XvuVvaFu/eKaWSVSC0OEdndlMXruOTpyfRok8HoIX1JTogPO6RqZ79bHEC8mSWX+6B6QPJezhcRqfb6tGvMA+d35/NF6xj6n1m6TbcCohkcfw5438xGBPuDgZGxC0lEpGqc2aMVi9ds5eH3FtCuSRo3Ht857JBqhH0mDne/38xmACcERXe7+7uxDUtEpGrceHwnlqzdwkPjFtA2M5Uze+z+uJnsLpoWB+7+DvBOjGMREalyZsZ953ZjWeE23aYbpWjGOEREarXkhHiGDepNdqN6XD0qj8VrtoQdUrWmxCEiAmSkJjE8uE13yMgpbNxetI8r6q4KJY7gGQ6tNS4itVK7Jmk8cWlvlqzdyo2aTfdHRfMA4HgzSzezxsAXwNNm9lDsQxMRqXqHd8jkT2cewvj5Bdz/zpdhh1MtRdPiaOjuG4nMKTXK3fvx3R1WIiK1ziX92vLzI9oybMJC/j11WdjhVDvRJI4EM2sBXAC8GeN4RESqhT+clsORHTP5f6/OYuqSdWGHU61EkzjuIjLfVL67TzGzDsBXsQ1LRCRcifFxPH5JL1pkpHDN6C9YUbgt7JCqjWjWHH/F3bu7+7XB/kJ3Pzf2oYmIhCsjNYlnfp7L9qISrhqVx9admmsV9pI4zOy3weujZvbI7lvVhSgiEp7OzRrw6MU9mbtyI795ZabmtGLvT47PC141rayI1GkDujbltoFdue/tL+nyfgNuOqFuz2n1o4nD3d8IXndNaGhmcUD94C4rEZE64+pjOjD/2008/N4CujSrz8ndWoQdUmiieY7j+eA5jjRgNjDXzH4T+9BERKoPM+Pec7rRs00GN788gzkrNoQdUmiiuasqJ2hhnAW8DbQHBsUyKBGR6iglMZ6nLu1Nw3qJXDN6Kuu31M3VA6NJHIlmlkgkcYxx9yJAo0MiUic1TU/hyUG9Wb1xBze+WDenJYkmcTwFLAbSgAlm1hbQGIeI1Fk9Wmdw91mH8PFXa3hw7Pyww6ly0TzH8Yi7t3L3UzxiCTCgCmITEam2LuzThov7tuGJ8V/z9qyVYYdTpaIZHG9oZg+ZWV6w/ZVI60NEpE774xk59Gidwa2vzOCrVZvCDqfKRNNVNRzYRGSuqguIdFON2OsVIiJ1QHJCPE9c2ot6SfFcM3pqnVnDI5rE0dHd7wymGlno7n8COsQ6MBGRmqBFw3o89rNeLFm3lVtenkFpHRgsjyZxbDOzo8t2zOwoQLN9iYgE+nXI5PZTDmbc3FU8Pj4/7HBibm9TjpT5JTDSzBoCBqwDLo9lUCIiNc3go9oxc1khfx23gENaNWTAQU3DDilmormrarq7HwZ0B7q5e093nxH70EREag4z475zutO1eTo3vTCNJWu3hB1SzERzV1WGmd0I/BH4s2bHFRHZs3pJkSfLzYxrRk+ttdOwRzPG8RbQDpgFTC23iYjIbtpkpvL3i3owf9Umhr46q1ZOwx7NGEeKu98c80hERGqJYw9qyi0nduHBsQvIbdeYQYe3DTukShVNi2O0mV1lZi3MrHHZFs2Hm9lAM5tvZvlmdtsejrc1s/fNbKaZjTez7HLH7jez2cF2Ybny54LPnG1mw4N5tEREqpVrj+3EsQdlcfcbc5m5rDDscCpVNIljJ/AAMInvuqn2ubiTmcUDjwEnAznAxWaWs9tpDwKj3L07kbXN7wuuPRXoBfQA+gG3mll6cM1zQFegG1APuDKKOoiIVKm4OOPhC3rQpH4S1z73BRu21p6HA6NJHLcAndy9nbu3D7ZoHgDsC+QHDw3uBF4EztztnBzgg+D9h+WO5wAT3L3Y3bcAM4GBAO7+VjBnlgOfA9mIiFRDjdKS+MclvVi1cTu3vDK91jwcGE3iyAe27sdntwKWlttfFpSVNwM4J3h/NtDAzDKD8oFmlmpmTYhMqti6/IVBF9Ug4J09fbmZXV02v1ZBQcF+hC8icuB6tWnE/zvlYN6bt5phHy8MO5xKEc3g+BZgupl9COwoK3T3Gyvh+28F/mFmlwMTgOVAibuPNbM+wESggEg3Wclu1z5OpFXy8Z4+2N2HAcMAcnNza0eaF5Ea6fIj25G3eD0PvDufnq0z6NchM+yQDkg0LY7/AvcQ+SNekdtxl/P9VkJ2ULaLu69w93PcvSdwe1BWGLze4+493P1EIk+sLyi7zszuBLIA3e0lItWemfGXc7vRpnEqN7wwjYJNO/Z9UTUWzZPjI8s2YAwwLXi/L1OAzmbW3sySgIuC63cxsyZmVhbDUCIz8WJm8UGXFWbWnchT62OD/SuBk4CL3b00mkqKiIStQUoij1/Siw3biriphq8cGM2T4+PNLD24BfcL4Gkze2hf17l7MXA98C4wD3jZ3eeY2V1mdkZw2rHAfDNbADQj0rIBSAQ+NrO5RLqbLg0+D+DJ4NxJZjbdzO6ItrIiImE6uEU6d595KBO/Xsvf3luw7wuqqWjGOBq6+8bgX/qj3P1OM5sZzYe7+1tEnjwvX3ZHuff/Bv69h+u2E7mzak+fGU3MIiLV0gV9WvP54nU8+kE+fdo15pguWWGHVGHRjHEkmFkLIos4vRnjeEREar27zzyUzk3rc/PL01m9aXvY4VRYNInjLiLdTV+7+xQz6wB8FduwRERqr3pJ8fzjZ73YtL2Ym1+qeYs/RTM4/oq7d3f3Xwb7C9393NiHJiJSex3UvAF/POMQPslfwxMffR12OBUSzeB4tpm9Zmarg+0/5eeUEhGR/XNRn9ac2r0FD41bQN7idWGHE7VouqpGELmNtmWwvRGUiYjIAYgs/tSNlhkp3PTidAq37gw7pKhEkziy3H1EMG9Usbs/S+ThOxEROUDpKYk8enFkPqvf/ntmjVi/I5rEsdbMLg0eyos3s0uBtbEOTESkrujROoPfDezK2LmrGP3ZkrDD2adoEscVRG7F/RZYCZwHDI5lUCIidc2Qo9sz4KAs/vzmPOas2BB2OHsVzV1VS9z9DHfPcvem7n6Wu39TFcGJiNQVcXHGg+cfRqO0RG54fhpbdlTf9cqjuatqpJlllNtvZGbDYxqViEgdlFk/mb9d2JNFa7fwh9dnhx3Oj4qmq6p72Yy1AO6+HugZs4hEROqwIzpmcsNxnXn1i+W8Pn35vi8IQTSJI87MGpXtBJMdar4oEZEYufG4TvRu24jfvzabpev2Zx292IomcfyVyEy0d5vZ3UTW5fi/2IYlIlJ3JcTH8bcLe+DAzS9Pr3ZTsEczOD6KyPKuq4LtHHcfHevARETqstaNU7n7rEOYsng9T4zPDzuc74mqy8nd5wJzYxyLiIiUc1aPVnz4ZQEPv/cVR3fOokfrjLBDAqLrqhIRkRCYGXefdSjN01O46cXqc4uuEoeISDXWsF4iD1/Yg6XrtvKnN+aEHQ6gxCEiUu31bd+Ya4/txMt5y3hr1sqww1HiEBGpCW46oTOHZTdk6KuzWLlhW6ixKHGIiNQAifFx/O2inhSVlIa+aqASh4hIDdG+SRp/PP0QJi1cyzOfLAwtDiUOEZEa5PzcbH6a04wH313A/G83hRKDEoeISA1iZtx7TjcapCRw88vT2VlcWuUxKHGIiNQwTeonc8/Z3ZizYiP/+OCrKv9+JQ4RkRpo4KHNObdXNo+N/5pp36yv0u9W4hARqaHuPCOHZg2SueXlGWzbWVJl36vEISJSQ6WnJPLA+YexcM0W7n/nyyr7XiUOEZEa7KhOTbj8yHY8O3Exn+avqZLvVOIQEanhfjewKx2y0rj1lRls2FYU8+9T4hARqeHqJcXz0AU9WL1pR5VMhBjTxGFmA81svpnlm9ltezje1szeN7OZZjbezLLLHbvfzGYH24Xlytub2eTgM18ys6RY1kFEpCbo0TqD647tyKtfLOed2d/G9LtiljjMLB54DDgZyAEuNrOc3U57EBjl7t2Bu4D7gmtPBXoBPYB+wK1mlh5ccz/wsLt3AtYDQ2JVBxGRmuT64zpzaKt0bn9tFms274jZ98SyxdEXyHf3he6+E3gROHO3c3KAD4L3H5Y7ngNMcPdid98CzAQGmpkBxwH/Ds4bCZwVuyqIiNQcSQlxPHRBDzbtKGboq7Nwj81EiLFMHK2ApeX2lwVl5c0gsp45wNlAAzPLDMoHmlmqmTUBBgCtgUyg0N2L9/KZAJjZ1WaWZ2Z5BQUFlVIhEZHqrkuzBtz60y6Mm7uK/05fHpPvCHtw/Fagv5lNA/oDy4ESdx8LvAVMBF4AJgEVerrF3Ye5e66752ZlZVVy2CIi1deQozvQu20j7nx9Dqs2bq/0z49l4lhOpJVQJjso28XdV7j7Oe7eE7g9KCsMXu9x9x7ufiJgwAJgLZBhZgk/9pkiInVdfJzx4PmH0bNNI4pjsG5HLBPHFKBzcBdUEnARMKb8CWbWxMzKYhgKDA/K44MuK8ysO9AdGOuRDrsPgfOCay4DXo9hHUREaqT2TdIYeUVfWmXUq/TPjlniCMYhrgfeBeYBL7v7HDO7y8zOCE47FphvZguAZsA9QXki8LGZzQWGAZeWG9f4HXCzmeUTGfP4Z6zqICIiP2SxGnWvTnJzcz0vLy/sMEREahQzm+ruubuXhz04LiIiNYwSh4iIVIgSh4iIVIgSh4iIVIgSh4iIVIgSh4iIVEiduB3XzAqAJftxaROgapbUqj5U57pBda4bDrTObd39B3M21YnEsb/MLG9P9zDXZqpz3aA61w2xqrO6qkREpEKUOEREpEKUOPZuWNgBhEB1rhtU57ohJnXWGIeIiFSIWhwiIlIhShwiIlIhShx7YGYDzWy+meWb2W1hx1OZzGy4ma02s9nlyhqb2Tgz+yp4bRSUm5k9EvwcZppZr/Ai3z9m1trMPjSzuWY2x8xuCsprc51TzOxzM5sR1PlPQXl7M5sc1O2lYIE1zCw52M8PjrcLtQIHIFgEbpqZvRns1+o6m9liM5tlZtPNLC8oi/nvthLHbswsHngMOBnIAS42s5xwo6pUzwIDdyu7DXjf3TsD7wf7EPkZdA62q4EnqijGylQM3OLuOcDhwHXBf8/aXOcdwHHufhjQAxhoZocD9wMPu3snYD0wJDh/CLA+KH84OK+muonIwnFl6kKdBwTLbJc9rxH7321311ZuA44A3i23PxQYGnZclVzHdsDscvvzgRbB+xbA/OD9U8DFezqvpm5Elho+sa7UGUgFvgD6EXmCOCEo3/V7TmSVziOC9wnBeRZ27PtR1+zgD+VxwJuA1YE6Lwaa7FYW899ttTh+qBWwtNz+sqCsNmvm7iuD998SWcYXatnPIuiO6AlMppbXOeiymQ6sBsYBXwOF/t0SzOXrtavOwfENRJZlrmn+BvwWKA32M6n9dXZgrJlNNbOrg7KY/24n7M9FUnu5u5tZrbtH28zqA/8BfuXuG81s17HaWGd3LwF6mFkG8BrQNdyIYsvMTgNWu/tUMzs25HCq0tHuvtzMmgLjzOzL8gdj9butFscPLQdal9vPDspqs1Vm1gIgeF0dlNeKn4WZJRJJGs+5+6tBca2ucxl3LwQ+JNJNk2FmZf9YLF+vXXUOjjcE1lZtpAfsKOAMM1sMvEiku+rv1O464+7Lg9fVRP6B0Jcq+N1W4vihKUDn4G6MJOAiYEzIMcXaGOCy4P1lRMYBysp/HtyNcTiwoVwTuEawSNPin8A8d3+o3KHaXOesoKWBmdUjMqYzj0gCOS84bfc6l/0szgM+8KATvKZw96Hunu3u7Yj8P/uBu19CLa6zmaWZWYOy98BPgdlUxe922IM71XEDTgEWEOkXvj3seCq5bi8AK4EiIn2cQ4j07b4PfAW8BzQOzjUid5h9DcwCcsOOfz/qezSRfuCZwPRgO6WW17k7MC2o82zgjqC8A/A5kA+8AiQH5SnBfn5wvEPYdTjA+h8LvFnb6xzUbUawzSn7W1UVv9uackRERCpEXVUiIlIhShwiIlIhShwiIlIhShwiIlIhShwiIlIhShwiIQtmOG0SvN8cdjwi+6LEISIiFaLEIVKFzOzSYK2M6Wb2VDCNv0iNosQhUkXM7GDgQuAod+8BlACXhBqUyH7Q7LgiVed4oDcwJZidtx7fTUAnUmMocYhUHQNGuvvQ7xWaXR5OOCL7R11VIlXnfeC8YO2EsrWh24Yck0iFqcUhUkXcfa6Z/Z7Iim1xRGYovi7ksEQqTLPjiohIhairSkREKkSJQ0REKkSJQ0REKkSJQ0REKkSJQ0REKkSJQ0REKkSJQ0REKuT/A0g5+k/tNtCuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(ell[:500], cl[:500]/cl_camb[:500])\n",
    "xlabel(\"ell\")\n",
    "ylabel(\"cosmosis camb / python camb\")"
   ]
  }
 ],
 "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}