examples/3x2demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cosmology with LSST\n",
    "\n",
    "Our goal is to constrain cosmological parameters with LSST using galaxy shapes and galaxy positions. We know these observables are sensitive to the growth and expansion history of the Universe. While we cannot predict where a given galaxy will form in the Universe, we can predict, given a cosmological model, the statistics of how galaxy pairs are distributed, and of how galaxy shapes relate to galaxy positions and to one another. In this example, we are going to choose a cosmological model and we are going to plot the angular correlations of those observables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using the Core Cosmology Library (CCL)\n",
    "\n",
    "The first thing to do will be to be able to call the library from this notebook. For that, we are going to import it as we import any other module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyccl as ccl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also import some additional modules we will need. These are standard ones that come with python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now need to decide which cosmology we are going to be working with. Let's start with a flat LCDM cosmology:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cosmo_fid = ccl.Cosmology(Omega_c=0.222, Omega_b=0.045, h=0.71, sigma8=0.806, n_s=0.963)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what's inside this object. It contains all the parameters that CCL needs to perform its computations. Because we haven't provided all of them, CCL is making some assumptions about what values to adopt. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pyccl.Cosmology(Omega_c=0.222, Omega_b=0.045, h=0.71, n_s=0.963, sigma8=0.806, A_s=None, Omega_k=0.0, Omega_g=None, Neff=3.046, w0=-1.0, wa=0.0, bcm_log10Mc=14.079181246047625, bcm_etab=0.5, bcm_ks=55.0, m_nu=0.0, mnu_type=None, z_mg=None, df_mg=None, transfer_function='boltzmann_class', matter_power_spectrum='halofit', baryons_power_spectrum='nobaryons', mass_function='tinker10', halo_concentration='duffy2008', emulator_neutrinos='strict')\n"
     ]
    }
   ],
   "source": [
    "print(cosmo_fid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining our sample of galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this cosmology in place, we can make CCL give us some basic computations. To do this, we are going to have to decide at what redshift we need to evaluate the functions. For the observables we are interested in, we expect that LSST will roughly cover the range between 0<z<3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "nz = 25\n",
    "z = np.linspace(0.1,3.,nz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The galaxies that we are observing have some redshift distribution, and it is through these galaxies that we plan to constrain cosmology. To illustrate our example, we are going to assume we have the same sample to probe galaxy shapes and galaxy positions, and we'll use pre-defined redshift distributions for this. (We can't measure the shapes for all galaxies for which we have positions, so usually the first will be a fraction of the second.) You could, if you wanted, pass your own redshift distribution to CCL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate\n",
    "\n",
    "def pz_tot(z, alpha=1.24, beta=1.01, z0=0.51):\n",
    "    \"\"\"Smail et al. 1994 redshift dist\"\"\"\n",
    "    return np.power(z, alpha) * np.exp(-np.power(z/z0, beta))\n",
    "\n",
    "def photoz_dist(z_photo, z_true, sigma_z0=0.05):\n",
    "    return (\n",
    "        np.exp(-0.5 * np.power((z_photo - z_true) / sigma_z0 / (1 + z_true), 2)) / \n",
    "        np.sqrt(2.0 * np.pi) / \n",
    "        sigma_z0 / (1 + z_true))\n",
    "\n",
    "\n",
    "def pz_joint(z_photo, z_true):\n",
    "    return pz_tot(z_true) * photoz_dist(z_photo, z_true)\n",
    "\n",
    "\n",
    "denom, _ = scipy.integrate.dblquad(\n",
    "    pz_joint,\n",
    "    # range in z_true (second arg to pz_joint)\n",
    "    0, 5,\n",
    "    # range in z_photo (first argument to pz_joint) as a function of z_true\n",
    "    lambda x: 0.1, lambda x: 3.0)\n",
    "\n",
    "dndz = []\n",
    "for _z in z:\n",
    "    val, _ = scipy.integrate.quad(pz_joint, 0.1, 3.0, args=(_z,))\n",
    "    dndz.append(val / denom)\n",
    "dndz = np.array(dndz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the redshift distribution we have defined:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEYCAYAAAB/QtA+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4FFXWx/HvCQgICAiELQIBx1dxGTEEkX3fZAREWQQUVNxwkM1BEUYUeIkog+DIMiDKrqAs4oxoElRAAUciqAguhHUgQlhkX4Sc94/q+GZCQrZOV3XnfJ6nH0n3re5fUZLTVffWvaKqGGOMMf4S5nYAY4wxocUKizHGGL+ywmKMMcavrLAYY4zxKyssxhhj/MoKizHGGL+ywmKMMcavrLAYY4zxKyssxhhj/Kqw2wHcUL58eY2MjHQ7hjHGBJWEhIRDqhqeVbsCWVgiIyPZuHGj2zGMMSaoiMju7LSzS2HGGGP8ygqLMcYYv7LCYowxxq88VVhEpL+I7BSRsyKSICKNL9N2tohoBo9TgcxsjDHmv3mmsIhId2AyMA64DVgHrBSRaplsMhConO6xA1ic/2mNMcZkxjOFBRgCzFbVmaq6TVUHAEnAExk1VtVjqvpL6gO4FqgJzAxcZGOMMel5orCISBGgDhCb7qVYoEE23+YR4HtVXefPbMYYY3LGK/exlAcKAQfSPX8AaJXVxiJSGugKPHeZNo8CjwJUq5bZ1bWCY9++fSxYsIDSpUsTERFBlSpVqFKlChUqVCAszBPfN4wxQcorhSWVpvtZMnguI71xCtO8TN9YdQYwAyA6Ojo77xmy4uLi6NmzJ4cOHbrktcKFC1OpUiWqVKnyXwUn9c+33norFSpUcCG1MSZYeKWwHAIuApXSPV+BS89iMvIIsERVj/g7WCi5ePEiY8aMYfTo0dx44418+umnlClThv3797Nv3z7279//X3/+6aef+PTTT/n1119/f4+SJUsyadIkHnroIUTExb0xxniVJwqLqp4XkQSgNfBumpdaA0sut62I3A7cCgzKv4TB7+DBg/Tq1Yv4+Hj69OnD1KlTKV68OADXXHPNZbc9ffo0SUlJ7Nmzh7Fjx9KvXz+WL1/OzJkzqVQp/XcBY0xB56WL6ROBviLST0RqichkoAowHUBE5orI3Ay2exT4GVgduKjBZe3atdx22218/vnnzJo1i7feeuv3opIdxYsX59prr6V58+bExcUxefJk4uPjufnmm3nvvffyMbkxJhh5prCo6iKcs46RwGagEXCnqqZOelbN9/idiFwF9ADeUNUC3W+SkZSUFF5++WWaN29OiRIl2LBhQ54vYYWFhfHUU0+xadMmatasSdeuXenVqxdHjx71Y3JjTDCTgvj7ODo6WkN9duMjR47Qt29fPvjgA+69915mzZpFqVKl/PoZFy5cICYmhtGjR1OxYkVmzZpF27Zt/foZxhjvEJEEVY3Oqp1nzliM/3z11VdERUXx0Ucf8dprr7F48WK/FxVwRpD99a9/5csvv6R06dK0a9eO/v37c+qUzapjTEFmhSWEqCqvv/46DRs2RFVZu3YtAwYMyPfRW1FRUSQkJDB06FCmT5/OrbfeyhdffJGvn2mM8S4rLCHixIkT3HfffQwYMIA2bdqwadMm6tWrF7DPL1asGBMmTOCzzz7j4sWLNGnShGeffZZz584FLIMxxhussISAM2fOcMcdd/Duu+8SExPDihUrKFu2rCtZmjRpwrfffstDDz3E+PHjqVu3Lj/88IMrWYwx7rDCEgKmTZvG1q1bWb58Oc8++6zrU7JcddVVzJw5k3/+858cOHCA5s2bs337dlczGWMCx0aFBbmTJ09Ss2ZNbr31VuLi4tyOc4mtW7fSpEkTSpQowdq1a22eNmOCmI0KKyCmTJlCcnIyo0ePdjtKhm688UZiY2M5duwYLVu2JCkpye1Ixph8ZoUliB0/fpyXX36Z9u3bU79+fbfjZCoqKooPP/yQpKQkWrduzeHDh92OZIzJR1ZYgtjkyZM5cuSIZ89W0mrQoAErVqxg+/bttG3blmPHjrkdyRiTT6ywBKlff/2Vv/3tb3Ts2JHo6CwveXpCixYtWLJkCd988w0dOnSwGymNCVFWWILUxIkTOXbsWFCcraTVoUMHFi5cyPr16+ncuTNnz551O5Ixxs+ssAShw4cPM2nSJO69915uvfVWt+PkWNeuXXnzzTeJj4+nW7du/Pbbb25HMsb4kRWWIDRhwgROnjzJCy+84HaUXOvTpw9Tpkzhgw8+4P777+fixYtuRzLG+IknFvoy2Xfw4EFee+01evTowU033eR2nDxJnbBy2LBhFC9enDfeeMP1mzuNMXlnhSXIjB8/nrNnzzJq1Ci3o/jFX/7yF06ePMno0aMpWbIkkydPtiWPjQlyVliCSFJSElOnTuX+++/n+uuvdzuO37zwwgucPHmSiRMnUrJkScaNG+d2JGNMHlhhCSIxMTH89ttv/PWvf3U7il+JCBMmTODUqVPExMRQsmRJnnvuObdjGWNyyQpLkNi7dy//+Mc/ePDBB7n22mvdjuN3IsLUqVM5deoUI0aM4JprruGBBx5wO5YxJhc81VMqIv1FZKeInBWRBBFpnEX7IiIy2rfNORHZIyJPBSpvII0bNw5VZeTIkW5HyTdhYWG8+eabNGvWjMcee4yEhAS3IxljcsEzhUVEugOTgXHAbcA6YKWIXG463LeBdsCjwPVAV+DbfI4acLt27WLWrFn069eP6tWrux0nX11xxRUsXryYChUqcPfdd3Pw4EG3IxljcsgzhQUYAsxW1Zmquk1VBwBJwBMZNRaRNkAr4E5VjVPVXar6pap+FrjIgTFmzBjCwsIYMWKE21ECIjw8nGXLlpGcnGw3UBoThDxRWESkCFAHiE33UizQIJPNOgNfAUNE5D8i8rOIvCYiJfMxasBt376dOXPm8PjjjxMREeF2nICJiopi5syZrF69mqefftrtOMaYHPBK5315oBBwIN3zB3DOSjJSE2gEnAPuAcoAfweqAPembywij+JcMguqxaZGjx5NkSJFePbZZ92OEnC9e/cmISGBSZMmUadOHevMNyZIeOKMJY30y1lKBs+lCvO91tN3Cexj4M/APSJS8ZI3Vp2hqtGqGh0eHu7X0Pnlhx9+YMGCBTz55JNUqlTJ7TiueOWVV2jevDmPPvooobLqpzGhziuF5RBwEUj/27MCl57FpEoC9qlq2oU9tvn+GzynJJfx4osvcuWVVzJs2DC3o7imcOHCLFq0iIoVK1pnvjFBwhOFRVXPAwlA63QvtcYZHZaRL4Aq6fpU/sf3393+TRh4W7ZsYdGiRQwcOJBgOcPKL6md+YcOHaJr167WmW+Mx3misPhMBPqKSD8RqSUik3H6S6YDiMhcEZmbpv1C4DDwlojcJCINcYYrv6eqQf+1dtSoUVx11VUMHTrU7SieEBUVxRtvvMGaNWvs78QYj/NK5z2qukhEygEjgcrAFpyhxKlnH9XStT8pIq1wOuy/Ao4Cy4Gg7+XetGkTS5cuZdSoUZQtW9btOJ7Rq1cvEhISePXVV4mKiqJv375uRzLGZEBUM+sbD13R0dHq5Y7gHj168PHHH7Nr1y5Kly7tdhxPuXDhAm3btuWLL75g7dq11K1b1+1IxhQYIpKgqlmuhe6lS2EG5xfnypUr6dq1qxWVDKR25leqVIkuXbpw4EBmYzuMMW6xwuIxX331FcePH6d16/TjGEyq8uXLs2zZMg4fPmx35hvjQVZYPCYuLg4RoUWLFm5H8bTbbrvt9878IUOGuB3HGJOGZzrvjSM+Pp46depQrlw5t6N4Xs+ePUlISGDixInUqVPHOvON8Qg7Y/GQEydOsH79elq1ymwWG5Pe+PHjadmyJY8//jjffPON23GMMVhh8ZTVq1dz4cIF61/JgcKFC7Nw4ULKli1Lt27dOHHihNuRjCnwrLB4SHx8PFdeeSUNGmQ2obPJSIUKFXj77bfZvn07jz32GAVxCL0xXmKFxUPi4uJo3LgxxYoVcztK0GnatCmjR4/m7bffZubMmW7HMaZAs8LiEfv27WPr1q12GSwPhg8fTps2bXjqqaesv8UYF1lh8YhVq1YBWGHJg7CwMObNm0e5cuXo2rUrx48fdzuSMQWSFRaPiIuLIzw8nFtuucXtKEEttb8lMTGRRx991PpbjHGBFRYPUFXi4+Np1aoVYWF2SPKqSZMmjB07lkWLFjFjxgy34xhT4NhvMQ/4/vvv+eWXX+wymB8988wztG3bloEDB7Jp0ya34xhToFhh8YC4uDgAuzHSj1L7W8qXL0+3bt2sv8WYALLC4gFxcXFcf/31VK1a1e0oISU8PJx33nmHnTt38sgjj1h/izEBYoXFZefOnWP16tV2GSyfNGrUiLFjx7J48WKmT5/udhxjCgQrLC7bsGEDp0+ftsKSj4YNG0b79u0ZNGiQ9bcYEwBWWFwWFxdHoUKFaNq0qdtRQlZYWBhz584lPDycrl27cuzYMbcjGRPSPFVYRKS/iOwUkbMikiAijS/TtpmIaAaPGwKZOa/i4uKoV6+erRaZz8qXL88777zDrl27rL/FmHzmmcIiIt2BycA44DZgHbBSRKplselNQOU0j5/zM6c/HT16lI0bN9plsABp1KgR//u//8u7777LtGnT3I5jTMjyTGEBhgCzVXWmqm5T1QFAEvBEFtsdVNVf0jwu5n9U//j0009JSUmxYcYB9Je//IU777yTwYMHk5CQ4HYcY0KSJwqLiBQB6gCx6V6KBbKaQ36jiCSJyCoRaZ4vAfNJXFwcV111FfXq1XM7SoERFhbGnDlzqFChAt26dbP+FmPygScKC1AeKAQcSPf8AaBSJtukns3cA3QBfgRWiUiTjBqLyKMislFENiYnJ/sndR7Fx8fTrFkzrrjiCrejFCip/S27d++mX79+1t9ijJ95pbCkSv8vXDJ4zmmo+qOqTlfVBFVdr6r9gY+ApzNpP0NVo1U1Ojw83L+pc2HXrl1s377dLoO5pGHDhsTExPDee+8xZcoUt+MYE1K8UlgOARe59OykApeexVzOl8B1/gqVn1KncbGOe/cMHTqUDh06MHToUDZu3Oh2HGNChicKi6qeBxKA9L9lW+OMDsuu2jiXyDwvLi6OiIgIbrghqEZHh5TU/paKFSvSrVs3fv31V7cjGRMSPFFYfCYCfUWkn4jUEpHJQBVgOoCIzBWRuamNRWSQiHQWketE5CYRiQE6A6+7kj4HUlJSWLVqFa1atUJE3I5ToJUrV45Fixaxd+9eHn74YetvMcYPPFNYVHURMAgYCWwGGgF3qupuX5NqvkeqIsAE4Ftgra99B1VdGrDQubRp0yaOHDlil8E8on79+sTExLB06VJef93z30uM8TwpiN/QoqOj1c1r6i+99BLDhw/nl19+oWLFiq7lMP9PVenUqRMfffQR69atIzo62u1IxniOiCSoapb/ODxzxlKQxMfHc8stt1hR8RARYfbs2VSuXNn6W4zJIyssAXbmzBk+//xzuwzmQWXLlv29v+Whhx6y/hZjcskKS4CtXbuWc+fOWWHxqDvuuIPx48ezbNky/v73v7sdx5igZIUlwOLj4ylSpAiNG2c6cbNx2eDBg+nYsSNPP/00//73v92OY0zQscISYHFxcTRo0IASJUq4HcVkQkR46623qFKlCt27d+fo0aNuRzImqOS4sIjINhF5S0RsAZEcOnjwIJs3b7bLYEEgtb/lP//5Dw8++KD1txiTA7k5Y7keeAD4UESK+zlPSFu1ahVg07gEi3r16vHyyy/z/vvvM3nyZLfjGBM0cnsp7ARQH1jhm/L+EiIyTUS+zHWyEBQfH8/VV19NVFSU21FMNg0aNIhOnToxbNgw628xJptyW1iWAn8HWgDvikihDNoUA+wuMx9VJS4ujhYtWlCoUEZ/XcaLUvtbIiIi6NatG0eOHHE7kjGel+vOe1UdCMwG7gLm+StQqPrpp5/Yu3evXQYLQldffTWLFi1i//799O7dm5SUFLcjGeNpeR0V1g94F+guIjP8kCdkxcfHA9j6K0Hq9ttv57XXXmPlypWMHj3a7TjGeFrhvGysqiki0gsoDjwsIidUdah/ooWWuLg4atSowbXXXut2FJNLjz32GF9++SUvvvgidevWpUOHDm5HMsaT8nwfi6peAO4FPgMGicgLeX3PUHPhwgU+/fRTuwwW5ESEqVOnUrt2bXr37k1iYqLbkYzxJL/cIKmq54COOCs4/lVE7Kwlja+++orjx4/bZbAQcOWVV7J06VJEhHvuuYfTp0+7HckYz8lNYXkZ+Dr9k6p6CmiPsz7Ky4D9FvWJi4tDRGjRooXbUYwf1KhRgwULFvDtt9/y+OOP282TxqST48Kiqs+qaoarIanqMZzlhH8EIvKYLWTExcVRp04dypUr53YU4yft27fnhRdeYN68eUybNs3tOMZ4it/nClPVQzhnK58Cp/z9/sHmxIkTbNiwwS6DhaCRI0fSoUMHBg0axPr1692OY4xnZFlYROSYiHwiIq+ISHcRyXJYk6ruV9WWQIGfT2zNmjVcuHDBOu5DUFhYGPPmzaNq1arce++9HDhwwO1IxnhCds5YSgDNgKHAQuAnETkiInEiEiMi94hI9Yw2VLv4zObNmwFn3ikTeq6++mqWLl3K0aNH6d69OxcuXHA7kjGuy05hKQU0Bf4CvAfsAcrgTOcyDFgM7BCRZBFZKSJjRKSTiOS4j0VE+ovIThE5KyIJIpKtRUtEpJGIXBCRLTn9zPyWmJhI5cqVbZr8EHbrrbcyY8YMVq9ezfDhw92OY4zrsrxBUlVPA2t9DwBEZAjwCnAWp9CUAioDbYE2adodUNUq2QkiIt2ByUB/4HPff1eKyI2quucy210NzAVW4cEBA9u3b7ebIguA3r17s2HDBiZMmMDtt99O165d3Y5kjGtysx5LF5yiMgEoq6q1VDUCqA48DxwHBDgKVMrBWw8BZqvqTFXdpqoDgCTgiSy2mwXMATzZe5qYmGiFpYCYOHEi9evX58EHH2Tbtm1uxzHGNbkZFTYC2KKqz/hujARAVfeq6ljgJuA74ABQMztv6Jt6vw4Qm+6lWKDBZbbrj1O8xuZoDwLkzJkz7N+/3wpLAVGkSBHeffddSpQowd13383x48fdjmSMK3JTWGoBWzN7UVX3A3cDfwD6ZvM9ywOFcIpRWgfI5KxHRG4BRgG9VPViVh8gIo+KyEYR2ZicnJzNWHmzY8cOACssBUhERASLFy9m+/bttvKkKbByU1iO4hSNTKnqDmA10DmH753+X6Fk8BwiUhR4B3haVXdm641VZ6hqtKpGh4eH5zBW7qTOJfWHP1z2r8uEmKZNm/Lyyy+zdOlSXnnlFbfjGBNwuSkssUCUiDTJot0xILtf1Q8BF7n07KQCl57FgDNQ4EbgLd9osAs4/Ts3+X5uk8E2AZdaWOyMpeAZPHgw3bp1Y/jw4Xz44YduxzEmoHJTWGKAc8B7ItI8owYiUgJn6eJsvb+qngcScKaDSas1sC6DTfYBtwC10zymA9t9f85om4BLTEykdOnSlC1b1u0oJsBEhDfffJPatWvTo0cPtmzx3Eh4Y/JNbuYK+wmn76QUEC8i74tIVxGJFJGKItIM+BfOWcVXOXjriUBfEeknIrVEZDJQBadgICJzRWSuL8Nvqrol7QM4CJzz/Xwyp/uVH1KHGouI21GMC0qUKMGKFSsoWbIkf/rTnzh48KDbkYwJiFzNFaaqi3HmA9uBszTxO0AisB/nfpImwG/Aczl4z0XAIGAksBloBNypqrt9Tar5HkHDhhqbiIgIVqxYwcGDB7n77rs5e/as25GMyXd5WfP+c+AGnLOX94G9ODdMHgKWAw1UNUeXpFR1qqpGqmpRVa2jqmvSvNZMVZtdZtsXVPXmnO9J/rhw4QK7du2ywmKIjo5mzpw5rFu3jkceecRGipmQl51JKBtJJtdyVPWiqs5V1S6+glBCVSv6fr5kzZaCZO/evVy4cMEKiwGga9eujB49mvnz5xMTE+N2HGPyVXbOWNYASSIyXUTaiUiW08AYGxFmLjVy5Eh69uzJiBEjWLJkidtxjMk32Sksf8e5xPUoTqd8sojME5EuIlI8X9MFMbuHxaQnIsyaNYs77riD+++/n6+/LtAn9SaEZVlYVHWgqkYCt+MsOfwL0At4F6fILBOR+0WkTL4mDTKJiYkULVqUiAjPzYtpXFSsWDGWL19OeHg4d911F/v373c7kjF+l+3Oe1XdqKrDVbUWznxgo3CWIO4EzAYOiEisiDwuIpXzJW0QSUxMpEaNGoSF+X2RThPkKlasyAcffMDx48fp2LEjp0+fdjuSMX6V2+HG21R1rKpGATVw1mr5CmeNlqnAXhH5QkSGZmfFyVBk0+Wby/njH//IwoUL+frrr+nTpw8pKSluRzLGb/L8dVpVd6vqRFVthHND4xM497JE40yv/5Nv/ZYCQ1XtHhaTpbvuuotXXnmF9957j1GjRrkdxxi/8esIL1U9CPwD+IeIlAY64sx0XKCuBx08eJBTp05ZYTFZGjJkCNu2bWPs2LHccMMN9OrVy+1IxuRZvg0dVtVjwDzfo0CxocYmu0SEqVOnsn37dh5++GFq1qxJ/fr13Y5lTJ5k5wbJz0XkNRHpKyK3iEihQAQLZjbU2OREkSJFWLJkCddccw2dO3dm165dbkcyJk+yc4mqAfAkzhLAm4ETIvKliEzzTRgZJSJX5GvKIJOYmIiIEBkZ6XYUEyTKlSvHP//5T86dO0f79u05dOiQ25GMybXsXArrBUT5HrWBq4G6vkfqpEe/ichWnKnvv/Y9vlHVAjnjXmJiIlWrVqVo0aJuRzFB5IYbbmDFihW0bduWO++8k1WrVnHVVVe5HcuYHMvODZJvq+pfVLWlqpbDWbyrK866LHE4k04WwSk6DwOv46yHclxEvs235B5mI8JMbjVp0oRFixbx9ddf06VLF86dO+d2JGNyLDfrsexU1SWqOkJV26lqRZzp7DsDo/FN+4JzNnSTX9MGCbuHxeRFx44dmTVrFvHx8dx///1cvHjR7UjG5IhfRoWp6n9E5ChQDGfJ4AI7h9iJEydITk62wmLypE+fPhw6dIinn36asmXLMm3aNFswzgSNPBUWEbka516VLjjLCBcFBDgFLPE9ChQbamz8ZejQoRw6dIiXXnqJ8PBwxowZ43YkY7Ilx4VFRCriXPa6B2jqew8BjuJMTLkU+Lggd9yDDTU2/jFu3DgOHTrE2LFjKV++PAMHDnQ7kjFZylZhEZFqOGclXXCGH4vv8QvO6pFLgU9V9UI+5QwadsZi/ElEmDZtGkeOHGHQoEGUK1eO3r17ux3LmMvKzg2S/wZ2An/DWYd+DzAJaAxEqOoTqhrnj6IiIv1FZKeInBWRBBFpfJm2TUVknYgcFpEzIvKDiDyd1wx5lZiYSPny5SlVqpTbUUyIKFy4MAsWLKB58+b07duXf/3rX25HMuaysjMqLNr3321AV1WtqapDVfUL9ePi3SLSHZgMjANuwxmyvNJ3tpSRk8BrQBOcAQNjgRdFpL+/MuWGDTU2+SF1HZfatWtz77338vnnn7sdyZhMZXe4sQC1gMW+M4RVIvKKiNwnItf7KcsQYLaqzvRNyz8ASMKZLfkSqpqgqu+o6ve+IdDzgY9xzqRcY4XF5JdSpUqxcuVKqlWrxp/+9Ce+/bZA3iZmgkB2CsvDwBRgA3Aa58775sBQYD6wVUROpJtT7I85mVNMRIoAdYDYdC/F4vTpZOc9bvO1XZ3dz/W38+fPs2fPHissJt+Eh4cTGxtLyZIladu2LTt27HA7kjGXyM6d92+p6gBVbQiUwrnp8X7gVWANcBwogfNL/c84c4ptwjenWDZzlAcKAQfSPX8AqHS5DUXkPyJyDtgITFXV6Zm0e1RENorIxuTk5GzGypldu3aRkpJihcXkq+rVqxMbG8v58+dp3bo1v/zyi9uRjPkvObrzXh3bVHWBr5+luapeDVwHdAfGA/HAYZybJaMv83YZfkS6nyWD59Jr7Pucx4FBInJ/JtlnqGq0qkaHh4fnMFb22IgwEyg33ngjH374IQcOHKBdu3YcPXrU7UjG/C7L4cYicgxncskEnLOCjaqamLaN7+dEnPtYUreritMJnx2HgItcenZSgUvPYv6Lqu70/fE73z02L+DSGjB2D4sJpHr16rF06VLuuusuWrZsSVxcHOXKlXM7ljHZOmMpATTD6VNZiLPU8BERiRORGBG5R0Sqp99IVfeq6orshFDV8ziFq3W6l1rjjA7LrjCcu/9dkZiYSIkSJahYsaJbEUwB06ZNG95//322bt1K8+bNOXjwoNuRjMnWDZKlcDrWb0/zqA608D0AEJEj+M5o+P8zm305yDIRmOe7b+YLnEtbVYDpvvefC6CqD/h+HoBzf82Pvu2bAE8DU3PwmX6VmJhIzZo1bU4nE1Dt2rXjn//8Jx07dqRZs2asWrWKypUrux3LFGBZFhZVPQ2s9T0AEJEhwCvAWZwbJksBlYG2QJs07Q6oapXsBFHVRSJSDhjpe68twJ2qutvXJP39LIVw+nQigQs4l+KexVeI3JCYmMj//M//uPXxpgBr1aoVK1eupEOHDjRr1oxPPvmEiIgIt2OZAirH0+aLSBecojIBKKuqtVQ1Aucs5nmcUWKpc4dddkRXeqo6VVUjVbWoqtZR1TVpXmumqs3S/DxJVW9S1RKqWlpVo3zbp+R0n/whJSXF7mExrmratCkff/wxSUlJNGnShN27d2e9kTH5IMeFBRgBbFHVZ1T191WIfH0qY3GGI3+H0+le0z8xvW///v2cO3fOCotxVcOGDYmLi+Pw4cM0bdrU7nMxrshNYakFbM3sRVXdD9wN/AHom7tYwceGGhuvqFevHp988gknTpygadOm/Pzzz25HMgVMbgrLUZyikSlV3YFzB3zn3IQKRjbU2HhJVFQUn3zyCWfPnqVp06Zs27bN7UimAMlNYYkFokSkSRbtjgEF5ut7YmIihQsXplq1zObMNCawbr31Vj777DNSUlJo1qwZW7ZscTuSKSByU1higHPAeyLSPKPD+7T+AAAWQklEQVQGIlICqJ/L9w9KiYmJVK9encKF/bLaszF+cdNNN7F69WoKFy5Ms2bN2Lx5s9uRTAGQ41/8qvoTTt9JKSBeRN4Xka4iEikiFUWkGfAvnCHDX/kzrJfZiDDjVddffz2rV6+mePHitGjRgo0bN7odyYS4XJ1RqOpioBWwA7gLeAfnPpL9wCqcmxV/A57zT0zvs8JivOwPf/gDa9asoXTp0rRs2ZL169e7HcmEsFxfqlLVz4EbcM5e3gf24twweQhYDjRQ1ZxMxxK0jhw5wtGjR62wGE+LjIxkzZo1hIeH07p1a1uJ0uSbPPWBqOpFVZ2rql18NzaWUNWKvp+/9ldIr7OhxiZYVK1albVr13L99dfTsWNHpk51bQYkE8IKTOd6frKhxiaYVK5cmdWrV9OhQweefPJJhg4dysWLF92OZUKIFRY/SC0sNWsWmIkGTJArWbIky5YtY8CAAUycOJGuXbty+vRpt2OZEGGFxQ8SExOpXLkyxYsXdzuKMdlWqFAhXnvtNSZNmsTy5ctp3rw5Bw5cdvkjY7LFCosf2IgwE8wGDhzIsmXL+O6777jjjjvYujXTGZuMyRYrLH5ghcUEu06dOrF69WrOnDlDgwYN+OSTT9yOZIKYFZY8OnPmDPv27bPCYoJe3bp12bBhAxEREbRt25Y5c+a4HckEKSsseZQ6LbkVFhMKIiMj+eKLL2jatCl9+/Zl1KhRqKrbsUyQscKSR3YPiwk1ZcqU4cMPP+TBBx9k9OjRPPDAA5w7dy7rDY3xsRkT88juYTGhqEiRIsyaNYtrr72WkSNHsmfPHpYtW0bZsmXdjmaCgKfOWESkv4jsFJGzIpIgIo0v07aLiMSKSLKInBCRL0WkYyDzglNYSpcubf/gTMgREUaMGMGCBQvYsGED0dHRJCQkuB3LBAHPFBYR6Q5MBsYBtwHrgJUiktkCJ02BT4AOvvYfAssuV4zyQ+qIMBEJ5McaEzA9e/Zk9erVXLhwgQYNGjBlyhTrdzGX5ZnCAgwBZqvqTFXdpqoDgCTgiYwaq+pAVX1JVf+tqttV9UUggQCvWmlDjU1BcMcdd7Bp0yZatWrFn//8Z3r06MHx48fdjmU8yhOFRUSKAHVwVqdMKxZokIO3ugpn6eSAuHDhArt27bLCYgqEcuXK8cEHHxATE8OSJUuIjo7mm2++cTuW8SBPFBagPFAISD+fxAGgUnbeQESeBK4B5mXy+qMislFENiYnJ+cl6+/27t3Lb7/9ZoXFFBhhYWE8++yzfPLJJ5w8eZI77riDN954wy6Nmf/ilcKSKv3/nZLBc5cQkXuAV4Beqro7wzdWnaGq0aoaHR4envek2FBjU3A1adKEzZs306hRIx555BH69OnDqVOn3I5lPMIrheUQcJFLz04qcOlZzH/xFZV5wAOquiJ/4mXMhhqbgqxChQp89NFHvPjii8yfP5+6devy/fffux3LeIAnCouqnsfpeG+d7qXWOKPDMiQi3YD5QF9VfS//EmYsMTGRokWLEhEREeiPNsYTChUqxPPPP09cXByHDx/m9ttvZ+7cuW7HMi7zRGHxmQj0FZF+IlJLRCYDVYDpACIyV0R+/z9WRHoAC4BngTUiUsn3CNgNJYmJidSoUYOwMC/9NRoTeC1btmTz5s3UrVuXPn360K9fP86cOeN2LOMSz/xGVNVFwCBgJLAZaATcmabPpJrvkepxnJkDJuEMS059LA1UZhtqbMz/q1y5MvHx8YwYMYJZs2ZRr149tmzZ4nYs4wLPFBYAVZ2qqpGqWlRV66jqmjSvNVPVZul+lgwezTJ673zIaoXFmHQKFy7M2LFjWblyJb/88gtRUVGMHj2a8+fPux3NBJCnCkswSU5O5uTJk1ZYjMlAu3bt+P7777n33nsZNWoUdevWtelgChArLLm0fft2wIYaG5OZ8PBwFi5cyPLly0lOTqZevXo899xznD171u1oJp9ZYcklG2psTPZ06tSJ77//nj59+hATE8Ntt93G+vXr3Y5l8pEVllxKTExERIiMjHQ7ijGed/XVVzNr1iw+/vhjTp8+TcOGDRkyZAinT592O5rJB1ZYcikxMZGqVatStGhRt6MYEzTatGnDli1beOKJJ3j11Vf54x//yGeffeZ2LONnVlhyyUaEGZM7V111FVOmTPm9oDRv3pz+/ftz4sQJd4MZv7HCkktWWIzJm6ZNm/Ltt98yZMgQpk+fzs0338xHH33kdizjB1ZYcuHEiRMcPHjQCosxeVS8eHH+9re/8cUXX1C8eHHat29Px44d+fHHH92OZvLACksu7NixA7Chxsb4S/369dm0aRPjx4/ns88+4+abb+app57i8OHDbkczuWCFJRfsHhZj/K9YsWIMGzaM7du38/DDDzNlyhSuu+46Jk2aZHfuBxkrLLlg67AYk38qVKjA9OnT+eabb6hbty6DBw/m5ptv5v3337cFxYKEFZZcSExMpHz58pQuXdrtKMaErNTO/A8//JDChQvTuXPn32dRNt5mhSUXbESYMYEhIrRv355vvvmG119/nW+//ZaoqCgefvhhkpKS3I5nMmGFJRessBgTWFdccQVPPvkk27dvZ8iQIcybN4/rrruOMWPG2N37HmSFJYfOnz/Pnj17rLAY44IyZcowYcIEtm3bRtu2bXn++eepWbMmr7zyit1g6SFWWHJo9+7dpKSkWGExxkXXXnstS5YsYc2aNdxyyy0MGzaMyMhIxowZw6+//up2vALPCksO2YgwY7yjcePGxMXFsWHDBho2bMjzzz9P9erVGTFiBMnJyW7HK7CssORQ6j0sNl2+Md5Rr149VqxYwaZNm2jbti0xMTFERkYyZMgQ9u/f73a8AscKSw4lJiZSokQJKlas6HYUY0w6tWvXZvHixb+vXvnaa69Ro0YN+vfvz65du9yOV2B4qrCISH8R2SkiZ0UkQUQaX6ZtZRFZKCI/iMhFEZkdiIyJiYnUrFkTEQnExxljcqFWrVrMmTOHn376ib59+/LGG29w3XXX8dBDD/HTTz+5HS/keaawiEh3YDIwDrgNWAesFJFqmWxSFDgEvAR8GZCQ2FBjY4JJzZo1+cc//sGOHTt48skneeedd6hVqxadOnUiNjaWlJQUtyOGJM8UFmAIMFtVZ6rqNlUdACQBT2TUWFV3qepTqjobOBKIgCkpKezYscMKizFB5pprrmHSpEns2rWL4cOHs2HDBtq2bcsNN9zAq6++ytGjR92OGFI8UVhEpAhQB4hN91Is0CDwiTKWlJTE2bNnrbAYE6QqVKjA2LFj2bNnDwsWLKBChQoMGTKEiIgI+vXrx6ZNm9yOGBI8UViA8kAh4EC65w8AlfzxASLyqIhsFJGNuR2GWLlyZXbu3En37t39EckY45KiRYvSs2dPPv/8czZt2kTv3r15++23iYqKon79+syfP59z5865HTNoeaWwpEo/dalk8Fzu3lh1hqpGq2p0eHh4rt4jLCyMyMhIypYt649IxhgPqF27NjNmzGDfvn1MmjSJI0eOcP/993PNNdcwfPhwdu/e7XbEoOOVwnIIuMilZycVuPQsxhhj/K5MmTIMHDiQH374gbi4OBo1asTLL79MzZo16dixI++//z6//fab2zGDgicKi6qeBxKA1uleao0zOswYYwJCRGjVqhXLli37vbP/q6++onPnzlxzzTU8/fTTfP/9927H9DRPFBafiUBfEeknIrVEZDJQBZgOICJzRWRu2g1EpLaI1AZKAWV9P98Y8OTGmJBUtWpVxo4dy969e/nggw9o1KgRkydP5uabb6ZevXpMnz7d5ibLgHhpRTYR6Q8MAyoDW4DBqrrG99pnAKraLE37jMLvVtXIy31OdHS0bty40T+hjTEFSnJyMgsWLODNN9/ku+++o1ixYnTp0oUHH3yQFi1aEBbmpe/r/iUiCaoanWU7LxWWQLHCYozJK1Xl66+/5s0332ThwoX8+uuvVK9enT59+tC3b19q1KjhdkS/y25hCd3Saowx+UhEqFOnDlOmTCEpKYl33nmH66+/njFjxlCzZk1atmzJ22+/zdmzZ92OGnBWWIwxJo+KFStG9+7d+fjjj9m1axdjxoxhx44d9OzZk4iICAYPHszWrVvdjhkwVliMMcaPqlWrxsiRI0lMTCQ2NpaWLVsyZcoUbrrpJho2bMhbb73FqVOn3I6Zr6ywGGNMPggLC6N169YsXryYffv2MWHCBA4fPsxDDz1ElSpVeOKJJ/j666/djpkvrLAYY0w+Cw8PZ+jQoWzbto01a9bQqVMnZs+eTZ06dYiKimLatGkcO3bM7Zh+Y4XFGGMCRERo3Lgxc+fOZf/+/bz++uukpKTQv39/qlSpEjITYVphMcYYF1x99dU8+eSTbNq0iX//+9/07Nnz94kwGzRowIIFC4J2IkwrLMYY4yIRoW7dusycOZN9+/bx6quvcujQIXr37k3VqlV57rnn2LNnj9sxc8QKizHGeESZMmUYNGgQP/zwA7GxsTRo0IDx48dTo0YNOnfuTFxcXFCsemmFxRhjPCZ1RNny5cvZsWMHzzzzDOvWraNNmzbUqlWLSZMmeXqOMissxhjjYdWrV2fcuHHs3buX+fPnU65cOQYPHkxERASPPPIIGzZswGtTc1lhMcaYIFC0aFF69erFunXrSEhI4L777mPhwoXUr1+fWrVq8dJLL7Fv3z63YwJWWIwxJuhERUXxxhtvkJSUxKxZs6hQoQLDhw+nWrVqtGvXjkWLFrk6R5kVFmOMCVKlSpXioYceYs2aNfz8888899xzbN26lR49elC5cmWeeOIJvvzyy4BfKrNp840xJoSkpKTw6aef8tZbb7F06VLOnDlDrVq16Nu3L71796ZKlSq5fm+bNt8YYwqgsLAwWrZsyfz580lKSmLmzJmULVuWZ555hqpVq9KjR498z1A43z/BGGOMK0qXLk2/fv3o168fP/30E3Pnzg3IZTG7FGaMMSZbgvJSmIj0F5GdInJWRBJEpHEW7Zv62p0VkR0i8nigshpjjMmYZwqLiHQHJgPjgNuAdcBKEamWSfsawIe+drcBMcDfReSewCQ2xhiTEc8UFmAIMFtVZ6rqNlUdACQBT2TS/nFgv6oO8LWfCcwBng5QXmOMMRnwRGERkSJAHSA23UuxQINMNqufQfuPgWgRucK/CY0xxmSXJwoLUB4oBBxI9/wBoFIm21TKpH1h3/v9FxF5VEQ2isjG5OTkPMY1xhiTGa8UllTph6hJBs9l1T6j51HVGaoararR4eHheYhojDHmcrxSWA4BF7n07KQCl56VpPolk/YXgMN+TWeMMSbbPFFYVPU8kAC0TvdSa5xRXxlZD7TKoP1GVf3NvwmNMcZkl2dukPQNN54H9Ae+wBn19TBwk6ruFpG5AKr6gK99DWALMBP4B9AQmArcp6pLsvisZGB3uqfL45w5hZJQ26dQ2x8IvX0Ktf2B0NunvOxPdVXNsi/BM1O6qOoiESkHjAQq4xSNO1U1tQBUS9d+p4jcCbyKMyR5P/BUVkXFt+0lfzEisjE7d5QGk1Dbp1DbHwi9fQq1/YHQ26dA7I9nCguAqk7FOevI6LVmGTy3GojK51jGGGNywBN9LMYYY0KHFZb/N8PtAPkg1PYp1PYHQm+fQm1/IPT2Kd/3xzOd98YYY0KDnbEYY4zxKyssxhhj/KrAFJZQW+slJ/sjIs1ERDN43BDIzJcjIk1EZIWI7PNl65uNbW4RkdUicsa33fMiIlltFwg53R8RiczkGLULUOTLEpHhIvKViBwXkWQR+UBEbs7Gdp48RrnZnyA4Rk+KyLe+fTouIutFpEMW2+TL8SkQhSXU1nrJ6f6kcRPOPUKpj5/zM2cOlcS5d2kgcCarxiJSCojDmfKnLvAU8Bec5Re8IEf7k0Y7/vsYfeL/aLnSDOdWgAZAC5ypk+JFpGxmG3j8GDUjh/uThleP0X+AZ3BuwYjGybVcRP6YUeN8PT6qGvIP4EtgZrrnfgZiMmk/Hvg53XNvAOvd3pdc7k8znIk5y7udPZv7dxLom0WbJ4DjwJVpnhsJ7MM3KMUrj2zuT6TvGEW7nTeb+1QSZ36/u0LkGGVnf4LqGPkyHwEeC/TxCfkzllBb6yWX+5Nqo4gkicgqEWmeLwEDpz6wVlXTng18DFTB+QUQrJaKyEER+UJE7nU7zGVchXPF4+hl2gTTMcrO/qTy/DESkUIi0gOnYGY232K+HZ+QLywEYK2XAMvN/qSuxHkP0AX4EVglIk3yK2QAZHaMUl8LNidxVj/tBtwJrAIWiUhvV1NlbjKwGWcy2MwE0zHKzv54/hj5+kxOAueA6cDdqvpdJs3z7fh4akqXfJZva724JNv7o6o/4hSTVOtFJBLnH8ma/AgXIF4/RtmmqoeAv6V5aqOIlAeGAfPdSZUxEZkINAIaqerFLJp7/hhld3+C5Bj9CNQGyuB8kZwjIs1UdUsm7fPl+BSEM5ZQW+slN/uTkS+B6/wVygWZHSPI2d+Dl3nuGInIq8B9QAtV3ZFFc88foxzuT0Y8dYxU9byqblfVjao6HOcsbHAmzfPt+IR8YdEQW+sll/uTkdo4l8iC1XqgsYgUS/Nca5xZrne5ksj/PHWMRGQy0BPnl/AP2djE08coF/uTEU8dowyEAUUzeS3/jo/boxYCNDKiO3Ae6AfUwrmeehJnbQGAucDcNO1rAKeASb72/Xzb3+P2vuRyfwYBnXG+Wd2EM3xagS5u70uajCVx/pHWBk4Dz/v+XM33egywKk370jjfuN4BbsbpOzoODHV7X3K5P31wfsnVAq7HuUx5Hhjs9r748k3x/f22wPmWm/oomaZN0ByjXO6P14/RS0BjnI73W3z5U4D2gT4+rv9lBPAvvT9OFT6H842/SZrXPgM+S9e+KfC1r/1O4HG39yG3+4NzDXg7zv0UR4C1OGvduL4faTI2wyl26R+zfa/PBnal2+YWnD6iszjfGkfhkWGsOd0f3y+trThfaI4DG4Hebu9HmnwZ7YsCL6RpEzTHKDf7EwTHaDbOAobngINAPNDWjeNjk1AaY4zxq5DvYzHGGBNYVliMMcb4lRUWY4wxfmWFxRhjjF9ZYTHGGONXVliMMcb4lRUWY4wxfmWFxRhjjF9ZYTHGGONXVliMMcb4lRUWY1wmIuNFRLN4/Jj1OxnjDQVpoS9jvGofMCeT1+4GSuGs+2FMULBJKI3xIBERnOUQBuAUlfaqmp312I1xnZ2xGOMxIlIImAk8iLMEwl2qetLVUMbkgBUWYzxERK7AWT+9G/Av4F5VPetuKmNyxjrvjfEI3xKxS3GKymLgbisqJhhZYTHGA0SkJM4Zyp+AN4H7VPU3d1MZkztWWIxxmYiUAWJx1l+fDPRT1RR3UxmTe1ZYjHGRiIQDnwL1gbGqOkhtqKYJctZ5b4xLRCQCiANqAcNU9RWXIxnjF3YfizEuEZFVOJe/9uMUmIycVdXHA5fKmLyzwmKMC0QkDDgGlMyiaYKqRgcgkjF+Y4XFGGOMX1nnvTHGGL+ywmKMMcavrLAYY4zxKyssxhhj/MoKizHGGL+ywmKMMcavrLAYY4zxKyssxhhj/MoKizHGGL+ywmKMMcav/g9HiFTc6JottgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(z, dndz, 'k-')\n",
    "plt.xlabel(r'$z$', fontsize=22)\n",
    "plt.ylabel(r'$dN/dz$', fontsize=22)\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3x2 correlations with CCL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With two samples of galaxies, there are three types of correlations we can perform. We could measure and model the auto-correlations of galaxy positions (\"clustering\"), the auto-correlations of galaxy shapes (\"cosmic shear\") and the cross-correlation between positions and shapes. Normally, to gain information on the expansion history of the Universe, we would split the sample into different redshift bins with a sufficient number of galaxies to have a significant measurement in each one. But for now, let's just take the full redshift distribution.\n",
    "\n",
    "Correlations are expressed in different forms. Here, we are going to express them in terms of \"angular power spectra\", or C_ell. Imagine that we took the sphere of the sky and expanded any function of the coordinates of the sphere into a basis of spherical harmonics. Each harmonic would contribute to the expansion with a given amplitude. What we are going to plot is the square of that amplitude as a function of multipole index."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Galaxy positions\n",
    "\n",
    "To model galaxy positions we need to define a \"bias\" parameter. This parameter tells us how the galaxies are connected to the density field. To make it simple, we'll take a one-to-one relation. Galaxies are simply tracing the density field in this model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bias_gal = np.ones(z.size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know need to define a convenience function, called a \"tracer\" which will store this information. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "gal_pos = ccl.NumberCountsTracer(cosmo_fid, has_rsd=False, dndz=(z, dndz), bias=(z, bias_gal))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"False\" statements above control the modeling of potential contributions to the signal, like redshift-space distortions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "gal_shapes = ccl.WeakLensingTracer(cosmo_fid, dndz=(z, dndz))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Angular power spectra\n",
    "\n",
    "We are now ready to compute angular power spectra, C_ell. These are a function of multipole number, with high ell correponding to small scales on the sky and low ell, to large separations on the sky."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ell = np.arange(100, 5000)\n",
    "cls_auto_pos = ccl.angular_cl(cosmo_fid, gal_pos, gal_pos, ell)\n",
    "cls_auto_shape = ccl.angular_cl(cosmo_fid, gal_shapes, gal_shapes, ell)\n",
    "cls_pos_shape = ccl.angular_cl(cosmo_fid, gal_pos, gal_shapes, ell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know have three correlations we can plot, and which we can fit to our observables. Varying the cosmology will result in these curves changing shape and moving up and down. See our likelihood notebook for an example of how to make fits to your data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEcCAYAAAB53pugAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd41FX2+PH3JYRQEhIgdAhVaRakSC/SWyghEJIFRBRQispvWVfsqy7qWla/awMUUCEB0pQiRSnSqyIl9JIQAqFJC5B6f3/cEEJPSDKfKef1PHkgM5+ZucExZ+69556jtNYIIYQQ9qyQ1QMQQggh7kWClRBCCLsnwUoIIYTdk2AlhBDC7kmwEkIIYfckWAkhhLB7EqyEEELYPQlWQggh7J4EKyGEEHavsNUDcBa+vr66evXqVg9DCCEcytatW09rrcve6zoJVnmklPIH/GvXrs2WLVusHo4QQjgUpVRsTq6TZcA80lrP11qP9Pb2tnooQgjhtCRYCSGEsHsSrPJIKeWvlJpy/vx5q4cihBBOS4JVHskyoBBCFDxJsBBC2KXU1FTi4+O5evWq1UMReeDm5oaPjw++vr4UKnT/8yMJVkIIuxQfH4+XlxfVq1dHKWX1cMR90FqTmppKYmIi8fHx+Pn53fdzyTJgHuV5z+rwYZBPjkLc4urVq5QpU0YClQNTSlGkSBEqV65MUlJSnp5LglUe5XnPKiQEfH2hf3/47js4dSp/ByiEA5NA5RzysvyX9Rz5MA6RF2++CUOGwIYNMGwYVKgAbdrAf/4De/aA1laPUAghLCfBymrdusFXX0F8PGzZAq+/DklJ8M9/Qr16UKcOTJgAq1ZBWprVoxVCFIBnn32Wd9555473T5o0iWeeecaGI7I/Sssn9zzJVm5pxP79+/PviePiYMECmDcPli+H1FQoXRp69gR/f+jaFUqWzL/XE8LO7N69m3r16lk9DJtbuXIlgwcPJj4+3uqh5Ks7/fdUSm3VWje51+NlZpVHBXbOys8PRo+GxYvh9GkIDzeBauFCGDjQ7HP16AHffiv7XEIIpyfByhGULAmBgfD995CYaJYEx42D3bvhmWfMPleHDvD553DsmNWjFcLpVa9enffee4/69etTqlQpnnrqqazzYFOnTqV27dqULl2a3r17k5CQAJg07vHjx1OuXDm8vb155JFH2LlzJwDDhg3jtddeIykpie7du5OQkICnpyeenp4kJCTw1ltvMXjw4KzXnzdvHg0aNMDHx4f27duze/fuG8b20Ucf8cgjj+Dt7U1QUFDW2E6fPk2vXr3w8fGhdOnStGnThoyMDFv9s+WJnLNyNIULmwSMNm3go49g2zaIioLISBPAxo2D5s1NdmFAANSsafWIhcgfL75o3u8FqWFD+PTTHF06a9YslixZQokSJfD39+fdd9+lQ4cOTJw4kaVLl9KgQQMmTJjAoEGDWLVqFUuXLmXVqlXs27cPb29v9uzZg4+Pzw3PWaJECRYtWnTXZcB9+/YRHBzMjz/+SPv27fnvf/+Lv78/MTExFClSBIC5c+eyePFiihYtSqtWrZgxYwbPPvssH3/8MVWqVOFU5mrMhg0bHCbjUmZWjkwpeOwxeOcdiIkxX//+N6SkwD/+AbVqXb9/1y7JLBQiH40dO5aqVatSunRpXn31VcLCwpg1axbDhw+nUaNGeHh48N5777F+/XqOHDmCu7s7Fy9eZM+ePWitqVevHhUrVsz1686ZM4eePXvSuXNn3N3dmTBhAleuXGHdunVZ1zz//PNUqlSJ0qVL4+/vz7bMIO/u7s7x48eJjY3F3d2dNm3aOEywkpmVM6lXz3y98gocOQLR0WbG9eab8MYbJrMwMBAGDYKHHrJ6tELkTg5nPLZStWrVrL9Xq1aNhIQEEhISaNSoUdbtnp6elClThmPHjtGhQwfGjh3LmDFjiIuLo1+/fnz00UeUzGWiVEJCAtWqVcv6vlChQlStWpVj2bYAKlSokPX34sWLZy1F/uMf/+Ctt96iS5cuAIwcOZKXX345dz+4RWRmlUd2W3W9enUYPx7WrDH7WF99BVWqwHvvwcMPQ4MG8PbbsHev1SMVwiEdPXo06+9xcXFUqlSJSpUqERt7vZdgUlISZ86coXLlyoCZ8WzdupVdu3axb98+Pvzww1ue914znZtfQ2vN0aNHs17jbry8vPj44485dOgQ8+fP55NPPmHZsmX3fJw9kGCVRw5Rdb1iRXj2Wfj1Vzh+HL780mQTvvUW1K1r1unfew8OHbJ6pEI4jC+++IL4+HjOnj3LpEmTCAoKIiQkhOnTp7Nt2zaSk5N55ZVXaNasGdWrV2fz5s1s3LiR1NRUSpQoQdGiRXFzc7vlecuXL8+ZM2e40wfggQMHsnDhQpYtW0Zqaioff/wxHh4etGzZ8p5jXrBgAQcOHEBrTcmSJXFzc7vtGOyRBCtXU64cPPcc/PabOYj86adQvLhZOqxVC5o2NYkbcXFWj1QIuxYSEkKXLl2oWbMmNWvW5LXXXqNjx46888479O/fn4oVK3Lw4EFmz54NwIULFxgxYgSlSpWiWrVqlClThgkTJtzyvHXr1iU4OJiaNWvi4+OTtYR3TZ06dZg5cybjxo3D19eX+fPnM3/+/KzkirvZv38/nTp1wtPTkxYtWjB69Gjat2+fL/8eBU0OBeeTJk2a6C1btlg9jPsXFwdz58KcOaaSBkCLFhAUZM513cdGsBB5Yc+HgqtXr84333xDp06drB6Kw5BDwSJ/+PmZsk6bN8OBAzBpEly+bNKFq1SBzp1Nod0LF6weqRDCBUmwEreqVQsmTjRnWnbvhldfhYMHTaHd8uVNNuH8+SZFXgghbECClbi7unVN1uDBg7BuHQwfbhI1eveGSpVMSah16+QMl3ApR44ckSVAG5NgJXJGKbOH9cUXJqNw/nyzNDhjBrRqZWZjr79u2poIIUQ+k2CVR3Z7zqogubtDr14QFmZqFX73HTzwgNnnqlcPGjeG//7XBDUhhMgHEqzyyCHOWRUkLy8YOhSWLDGp8P/9r5mF/b//ZxIzunWD0FCTrCGEEPdJgpXIPxUrmuzBLVtMYsbEiWZZ8G9/M5Xhhw8357scpMqzEMJ+SLASBaNuXXj3XVMVY8UKU5MwPBzatzeV4F9/HfKzWaUQwqlJsBIFq1AhE6CmTTP7WzNnmoK6kybBgw9Cy5bw9dfw119Wj1SIPLu575TIPxKshO0UL26WBJcsMRUzPvgAzp835Z8qVIABA0yWYWqq1SMVQtgZCVbCGpUrw0svwc6dsHXr9XqFvXub+158EX7/Xc5vCSEACVbCakpBo0amoO6xYzBvHrRrZ1qaNG5s2pl8+CHcVMxTCKt98MEHVK5cGS8vL+rUqZPVaiMlJYWhQ4fi5eVFgwYNyF4z9P3336dWrVp4eXlRv359oqOjs+6bMWMGrVq1Yty4cXh7e1O3bt0b2necP3+ep59+mooVK1K5cmVee+010tPTAThw4ADt2rXD29sbX19fgoKCbPSvYDvSfFHYD3d38Pc3X2fPmsK6339vZmAvvwydOpk0+X79zJKicCkvvvhiVsfbgtKwYUM+zUGTx7179/L555+zefNmKlWqxJEjR0hPT2f16tXMmzePqKgopk+fzmuvvcbYsWPZsGEDALVq1WL16tVUqFCB8PBwBg8ezIEDB7I6Bm/cuJHAwEBOnz5NVFQUAQEBHD58mNKlS/Pkk09Svnx5Dhw4QFJSEr169aJq1aqMGjWK119/nS5durBixQpSUlJw6KLadyAzq5sopR5WSm3L9nVFKdXX6nG5nNKlTQ+udetg3z7TwmTvXhg82NQnHD4cVq6UNHhhCTc3N5KTk4mJiSE1NZXq1atTq1YtAFq3bk2PHj1wc3NjyJAh/Pnnn1mPGzBgAJUqVaJQoUIEBQXxwAMPsGnTpqz7y5Urx4svvoi7uztBQUHUqVOHhQsXkpiYyKJFi/j0008pUaIE5cqVY/z48VntR9zd3YmNjSUhIYGiRYvSunVr2/6D2IDMrG6itd4BNARQSnkBh4FfLB2Uq3vgAXjnHfjXv2D1alMxIzwcpk+HatVgyBDz9eCDVo9UFKCczHhspXbt2nz66ae89dZb7Nq1i65du/LJJ58At7aUv3r1KmlpaRQuXJjvv/+eTz75hCNHjgBw6dIlTp8+nXV95cqVb+gUXK1aNRISEoiNjSU1NTVrBgaQkZFB1apVAfjPf/7D66+/zuOPP06pUqX4+9//zvDhwwvyn8DmZGZ1d32AX7XWSVYPRGDS4Nu1u54GP2vW9TT4OnVM7cIvv4QzZ6weqXABISEhrFmzhtjYWJRS/POf/7zr9bGxsYwYMYLPP/+cM2fOcO7cOR566CGy9xQ8duzYDd/HxcVRqVIlqlatioeHB6dPn+bcuXOcO3eOCxcusGvXLsAEyKlTp5KQkMDkyZMZPXo0Bw4cKJgf3CIOFayUUm2VUvOUUseUUlopNew214xWSh1WSl1VSm1VSrXJw0sGAXPy8HhRUIoXh5CQG9PgL16EMWNMJY1+/SA6GpKTrR6pcEJ79+5l+fLlJCcnU7RoUYoVK3bP9vBJSUkopShbtiwA06dPZ+fOnTdcc/LkSf7v//6P1NRUwsPD2b17Nz169KBixYp06dKFv//971y4cIGMjAwOHjzIb7/9BkB4eDjx8fEAlCpVCqWUw7SrzymHClaAJ7ATeAG4cvOdSqkg4DNgEvAYsA5YpJTyy3bNzjt8Vb3puXyA5sCigvtxRL64lga/Y4dJdx87Ftavh4AAE7hGjzbfSxq8yCfJycm8/PLL+Pr6UqFCBU6ePMmkSZPu+pj69evz97//nRYtWlC+fHl27NhBq1atbrimWbNm7N+/H19fX1599VUiIiIoU6YMAN9//z0pKSnUr1+fUqVKERgYyPHMYtGbN2+mWbNmeHp60rt3bz777DNq1KhRMD+8RRy2rb1S6hIwVms9I9ttG4HtWusR2W7bD0RorSfm8vmHA+211kPvcs1IYCSAn59f49jY2Nz9EKLgpKXBL7/ADz+YGdbVq1C7ttnbGjzYlHwSds2e29oXhBkzZvDNN9+wZs0aq4dSIKStfSalVBGgMbD0pruWAi3v4ymDgNl3u0BrPUVr3URr3eTa1F7YicKFoXt3U/E9MdHsc1WpAm++aXpvtWkDU6ZImSchHITTBCvAF3ADEm+6PRGocOvld6aU8sVkBEoWoDMoWRKeesoU1I2NNQkZp0/DqFFmmXDAAHMYOSXF6pEKIe7AmYLVNTeva6rb3Hb3J9D6tNa6vNb6nkXqXLL5oiPz8zOtS2JiYPNmGDnSnNfq08fsfY0bZ2530OVx4biGDRvmtEuA+cGZgtVpIJ1bZ1HluHW2lW9cvvmio1IKmjSB//s/U8pp/nx44gmYOhUef9x0PP73v81MTAhhOacJVlrrFGAr0PmmuzpjsgILhMysnIC7O/TqZco7nThh9rLKlYPXXoPq1a+3OLlwweqRCuGyHCpYKaU8lVINlVINMWP3y/z+Wmr6J8AwpdQzSql6SqnPgErA1wU1JplZORkfHxgxAlatMo0j337bzLyeftqUeQoOhp9/NtmGQgibcahgBTQB/sj8Kgb8K/PvbwNorecALwKvAduA1kAPrbWs5Yjcq1HDdDTeu9ec0xo+HJYuhZ49zf7W+PHSxkQIG3GoYKW1Xqm1Vrf5Gpbtmi+11tW11h5a68Za61UFOSZZBnQBSkHz5vDFF3D8uDm31bq1+f5aG5MPPoCjR60eqRBOy6GClT2SZUAXU6QI9O0LkZFmf+urr0xq/Msvm6K67dubPa+zZ60eqRDMmDHDaSqwS7AS4n5lb2Oyfz+89ZaZeY0aBRUqmKA2dy5cuaUymBAilyRY5ZEsAwrAlHJ64w3Yswe2bDHntTZtgqAgk5jx5JNmv0sSM4S4LxKs8kiWAcUNlDL7WB9/bPawli0zFTJ++gm6djUln154ATZulMQMB1a9enXee++9rKKyTz31FFevXs26f+rUqdSuXZvSpUvTu3dvEhISANBaM378eMqVK4e3tzePPPLILZXXAc6ePUuVKlWYP38+YPpe1a5dm++///6245kxYwY1a9bEy8uLGjVqMGvWrBvunzBhAqVKlaJGjRosWnS9Nvf06dOpV68eXl5e1KxZk8mTJ2fdt3LlSqpUqcKkSZPw9fWlevXqNzxvcnIyEyZMwM/Pj/Lly/Pss89ypQBXEaT5ohAFxc0NOnQwX198YVLeQ0Nh8mRzGLl2bdPmJCTE9OMSd/Xii1DAXe1p2BBy2uNx1qxZLFmyhBIlSuDv78+7777Lu+++y/Lly5k4cSJLly6lQYMGTJgwgUGDBrFq1SqWLl3KqlWr2LdvH97e3uzZswcfH59bnrt06dJMmzaNoUOHsn37dl599VUaNmzI0KG31tVOSkri+eefZ/PmzdSpU4fjx49zNtue6caNG3nyySc5ffo0U6ZM4emnn+bYsWMopShXrhwLFiygZs2arFq1iu7du9O0aVMaNWoEwIkTJzh9+jTHjh1jw4YN9OjRgyZNmlCnTh3++c9/cujQIbZt24a7uzshISG8/fbbvPfee/f3j38PMrPKI1kGFDlStKhpWRIRYRIzpk0zCRnvvAN165pqGp98Ys50CYcwduxYqlatSunSpXn11VcJCwsDTBAbPnw4jRo1wsPDg/fee4/169dz5MgR3N3duXjxInv27EFrTb169W7o/ptdly5dGDBgAB07dmThwoU3zHpuVqhQIXbu3MmVK1eoWLEiDRo0yLqvWrVqjBgxAjc3N5588kmOHz9OYqIp6tOzZ09q1aqFUop27drRpUsXVq9efcNzv/POO3h4eNCuXTt69uzJ3Llz0VozdepU/vvf/1K6dGm8vLx45ZVXmD37rrW/80RmVnmktZ4PzG/SpMmIe158G2+/bf4MCIAGDcwqknByPj6msO5TT5ngNGeOmXH9/e8wYYIp+/S3v5kGkqVKWT1au2FHXe0BslrKw/X28wAJCQlZMxMAT09PypQpw7Fjx+jQoQNjx45lzJgxxMXF0a9fPz766CNKlix529cYOXIkn3/+Oa+88kpWX6ublShRgjlz5vDRRx/x9NNP06pVKz7++GPq1q0LmC7C1xQvXhwwy4oAixYt4l//+hf79u0jIyODy5cv8/DDD2ddX6pUKUqUKHHLz3nq1CkuX75M48aNs+7TWpOenp6zf7z7IDMri23dapLIHn7YrAS9/LLZl5ftDBdRqZI5XLx5s0nOeOMN0/n4WsWMXr1MTy4p9WR3jmY7V3et/TxApUqVyN7bLikpiTNnzlC5cmUAnn/+ebZu3cquXbvYt28fH3744W2fPz09nVGjRjF06FC++uqru7ap79q1K7/88gvHjx+nbt26jBhx78/OycnJ9O/fnwkTJpCYmMi5c+fo0aMH2Xsc/vXXXyQlJd3yc/r6+lKsWDF27drFuXPnOHfuHOfPn88KggVBgpXFfvrJfLj++mtTMOHjj6FZM1Mc/PnnTUFwSSBzEXXqmE8u+/aZjMIXXzTdj4cONbUK+/WD2bMh2y8PYZ0vvviC+Ph4zp49y6RJkwgKCgIgJCSE6dOns23bNpKTk3nllVdo1qwZ1atXZ/PmzWzcuJHU1FRKlChB0aJF79h+/lrn4WnTpjFhwgSGDh1625lLYmIi8+bNIykpCQ8PDzw9PXPU0j4lJYXk5GTKli1L4cKFWbRoEUuX3twOEN58801SUlJYvXo1CxYsYMCAARQqVIgRI0Ywfvx4Tp48CcCxY8dYsmRJjv/9ck1rLV/58NW4cWOdH86e1fr777Xu21frokW1Bq19fbV++mmtFy7U+urVfHkZ4SgyMrRet07rF17QulIl84YoVkzrAQO0jojQ+vJlq0dYYGJiYqwewh1Vq1ZNT5o0SderV097e3vroUOH6qSkpKz7v/rqK12zZk1dqlQp3bNnT3306FGttda//vqrfvjhh3WJEiV0mTJldEhIiL548eItz79lyxbt4+Oj9+/fr7XWOi0tTbds2VK/++67t1ybkJCg27Ztq0uWLKm9vb11u3bt9K5du7TWWk+fPl23atXqhuuBrOf9/PPPdbly5bS3t7cePHiwDgoK0q+++qrWWusVK1boypUr63fffVeXKVNGV61aVX///fdZz3PlyhU9ceJEXaNGDe3l5aXr1q2rP/vsszv+m93pvyewRefgd6zDtrW3F0opf8C/du3aI/bv35+vz52UBIsXQ1QULFhgVoK8vExpuoAA6NEDsi0nC2eXkQFr1pg9rvBwOHUKPD2hd28YNAi6dAEPD6tHmW/sua199erV+eabb+jUqZPVQykwK1euZPDgwcTHx+fL80lbe4vpAjxnVaIE9O8Ps2bByZMm8zkoCH79FQYONCtDQUGm8o8USXABhQpB27YmDT4hAX75xQSpxYtNwCpf3iRtLF4MqffsGyqEQ5Fg5SA8PKB7d9Mb8PhxWL7cFEVYsQICA6FsWXNc56efINvZROGsCheGTp3MG+LECfNJpm9fMw3v3h0qVjRdkJctk01P4RRkGTCfNGnSRG/ZssXmr5uWBr/9ZlaGoqLgzBmzVNinj5l1de7sVCtD4l6Sk2HJEvOG+Okns5bs62uSMwYMMIV23d2tHmWO2PMyoMg9WQZ0cYULQ8eOptD38eNmBWjAAFi4EPz9r68MLVoEKSlWj1YUOA8PsyQ4a5bZ04qIMJ9YwsLMnlaFCiYtXt4QwsFIsHIi7u6m/Ny335qVoYULzcpQdLRJxqhQAZ55xiwhFuDZPWEvihUzm56hoWbT88cfzRshIsL8Wa6cWUtesMDMyOyQrPw4h4yMjDw/hywD5lFBZgPml+RkU/B77lzz++rSJXMWNTjY7HM99phUznApyckmSyc83CwVnjtn1o579zYboF27mkBnscOHD+Pl5UWZMmVQ8gZ1SFprUlNTSUxMRGuNn5/fLdfkdBlQglU+sWrPKrcuXzYfpGfNMitBqammNN21eqq1alk9QmFTKSlmqh0RYabgZ8+aNNRevUzg6t7dsvMRqampxMfH31DNXDiewoUL4+3tja+vL4UK3bqYJ8HKxhwlWGV39qz5HTVrFqxaZW5r3twEraAgs0okXEhqqsnWiYgw2TqnTkHx4mbJMDDQHPDz9LR6lMLJSLCyMUcMVtkdPWr24END4c8/TXeLTp1MPdW+fc0qkXAhaWmwevX1wHXihKkc362bCVy9eoH0cBP5QIKVjTl6sMpu1y4TtEJD4cgRs33Ru7eZcXXrBkWKWD1CYVPp6bBunQlcERHmQHKRIibLsH9/c06idGmrRykclAQrG3OmYHWN1uZ3VGioObZz5oz5nRQUBIMHQ4sWkpjhcjIyYMMGUzYlMhJiY683mezf30zDy5e3epTCgUiwsjFnDFbZpaaajMJZs0xG4ZUrpkr83/5mvjJb5whXorXpcRMZaWZcBw6YklBt2pilwn79ILMthhB3IsHKxpw9WGV38aJJHJs501TzycgwjW4HDzal6uSDtQvS2rQzuRa4YmLM7S1amMDVv7/pjCzETSRY2YgjnLMqSMePmxZLM2fC77+bFaHOna8nZkjymIvavfv6UuG2bea2Jk1M0OrfHx54wNrxCbshwcrGXGlmdScxMWaZcNYss5VRvLhZCRo82GQWFi5s9QiFJQ4evB64Nm0ytz3yiAlagYFQv7614xOWkmBlYxKsrsvIgLVrTdCaOxf++suc2Ro0yASuJk0kMcNlxcWZVPiICJO9o7XZ8Ly2VPjoo/LmcDESrGxMgtXtJSebShkzZ14vQffggyZo/e1vULOm1SMUljl+3Gx+RkSYw8gZGeYNcS1wNW0qgcsFSLCyMQlW93bunFkJmjkTVq40t7VsaYLWwIGmk4VwUadOmTTTyMjrPbj8/ExL7MBAk6hxm1I9wvFJsLIxCVa5ExdnKmb88IM5hFy4sClDN3iwaW1iB3VUhVX++gvmzTOBa8kSU7+wYkWzARoYaFLjZQPUaUiwsjEJVvdHa9i+/XpiRkKCKe0UGGhmXO3bmwxD4aIuXDC9biIjTTfkK1fMFLxvX/Mm6dDBYZpJituTYGVjEqzyLj3dbF3MnGm2MS5eNK1MQkLMjOuRR2QLw6UlJZnuohERZgP00iXw8THlnvr3N2cmiha1epQilyRY2ZgEq/x15QrMn29mWz//bLYwHnrIBK2QEKha1eoRCktdvWpKqkRGmiXDaz25evUygat7d3N2Qtg9CVZ5oJQaAzyX+e0mYITW+q69dSVYFZzTp02fwJkzTbYzQLt2JnAFBpoP18KFXevJFRlpsgvPnDGBqnv3661NpG2A3ZJgdZ+UUr7ARqABkAzMAyZrrRfc7XESrGzj4EFTWHfmTNi3zxT/9vc3+1s9eoCHh9UjFJZKSzPN2SIjr7c28fCALl1M4PL3h1KlrB6lyEaC1X1SSpUDNgOPAknAImCi1nrz3R4nwcq2rtVQnTnTZBWePGlmWAMHmhlXq1aS6ezy0tNh/XqzxxUZCfHxJouwU6frFeLlvITlchqsHOp/Z6VUW6XUPKXUMaWUVkoNu801o5VSh5VSV5VSW5VSbXLzGlrrk8BHQBxwAthzr0AlbE8pUwnj00/h2DFz8LhnTxO82rY1Z0tfffV6PVXhgtzcoHVr8yaJjTWtTcaPN1PyESOgQgXo2BG+/NIcUBZ2zaGCFeAJ7AReAK7cfKdSKgj4DJgEPAasAxYppfyyXbPzDl9VM+8vBfQEqgNVgIeUUu0L9scSeVG4sGkKOXMmJCaaP+vVg/ffhwYNoFEj+OQT+X3k0goVgmbN4D//Ma1Mfv8dXn7ZfNIZM8a0MmnTxgS2o0etHq24DYddBlRKXQLGaq1nZLttI7Bdaz0i2237gQit9cQcPu8AoL3Wekzm9//A/Dv95zbXjgRGAvj5+TWOjY3Nw08k8tuJE6Zp5MyZsGWL+X3VsaNZJuzXT/bcBWY9OSbmemuTHTvM7Y8/fr3sk9QEK1BOuQx4N0qpIkBjYOlNdy0FWubiqY4CLZRSRZVSbkB7YO/tLtRaT9FaN9FaNylbtux9jFoUpAoV4IUXYPNm2LOn/FKjAAAgAElEQVTHLAseOABPPml6bgUHm/OmqalWj1RYRikz/X7jDXM6fe9emDTJ7He99BLUqmWm5v/+t7lPWMZpghXgC7gBiTfdnghUyOmTaK03YJIq/gC2AwcxGYG3pZTyV0pNOX/+fK4HLGynTh14+22TTbh2LQwbZo7p9OplDh6PGwcbN5oP2sKFPfggTJxopuKHDsFHH5mDxq+9ZqrDN2wIH3xg9sCETTnNMqBSqhJwDGirtV6d7bo3gWCtdYE2XpdsQMeTkmJKz82cac6VXr1qPkhfqwgv/QFFlvh4s1QYFmY+1YCpwhwcDAMGSHvsPHC5ZUDgNJDOrbOoctw62xIi64zWnDlmf2vaNNN5/e23zQfs5s3h889NQXDh4qpUMWvKGzaY6fm//23qgY0bZ6bmnTubN9C5c1aP1Gk5TbDSWqcAW4HON93VGZMVWCBkGdA5eHvDU0+Z7hRxcfDhh2amNW6cKfjds6f5UH35stUjFZarWRNeecXsce3caZYNDx2Cp582M6y+fc0nIHmz5CuHWgZUSnkCtTO/XQe8j9lPOqu1jstMXf8BGA2sBZ4FngYaaK0LdJFZlgGd044d1yvCx8eDp6dpsTR4sCn4LRXhBWA2OzdvNp9o5swx5yRKlIDevc1SYdeuZiovbuGUFSwyzzutuM1d32mth2VeMxp4CaiIOZM1Xmu9qgDH5A/4165de8T+/fsL6mWExTIyYPVqs78VHg7nz5tsw+BgE7gee0wqwotM6enmzRIWZtLhz541JZ769zeboW3bSnmVbJwyWNkzmVm5jqtXTcr7zJnXU98ffNCUegoKMpnQErgEYLJ4fvnFBK4ffzRtTvz8TNAaPBjq17d6hJaTYGVjEqxc09mz5sPznDmwcqWZgdWrZ4LWwIHm70IAZg/rp59Me+ylS80MrFEjGDLETNFdNKNQgpWNyDKguCYx0WQ3z51rCn9rbXpwXQtcDz5o9QiF3UhMNLOtH34wpZ/c3Exl+CFDTDNJF+rFJcHKxu53ZrVz506qVKmCjzRlcirHj5sZ19y5sGaNue3RR03gGjAAate+++OFC4mJMWvKM2eauoSenmZ/a8gQaN/e6bN4JFjZ2P0Gq8aNG7Njxw46d+5MYGAgffr0oXTp0gUwQmGV+PjrgWv9enNbo0bXA1eNGtaOT9iJjAwzJf/hB/OGuXDBFNi9tr/18MNWj7BA2CRYKaUqaq2lljX3H6w2b95MeHg44eHhHDlyhMKFC9OxY0cGDBhA3759KVOmTAGMVlglNvZ64Nq0ydzWtOn1wOXnd/fHCxdx5YopqzJzJixebJpKNmpkDgOGhIATfaC1VbDKAE5i6uhty/zzD621y2ze5Neeldaa33//PStwHTp0CDc3Nzp06EBgYCD9+vVDiuU6l8OHTRr83LmmkSSYYt/9+5uvWrWsHZ+wE6dOmf2tGTPgjz/Mea0+fUzg6tLF4ZcJbRWsumM66jbE9I+qBShMh93t3BjEdmZWmXBK+ZkNqLVm27ZtWYHrwIEDuLm50b59ewIDAwkICKBcuXL58lrCPhw4YAJXZOT1wPXoo9cDl2Q4CwC2bYPp080p9TNnTKmnJ580lZkdNIPHkj0rpVRxTPB6DBPAGgIPAUWBVK21R769mJ0pqNR1rTXbt2/PClz79u2jUKFCtG3blgEDBtC/f3/Ku2jKq7M6cgSiokzgWpdZKKxeveuB69FH5RyXy0tOhgULTOBatMjsd7VqZWZbAwc6VLO2fA1WSqlawBmtda6rNCqlCgH1gUe11rNy+3hHYYtzVlprdu7cmRW49uzZQ6FChejQoQODBg0iICCAUqVKFegYhG0lJEB0tNnnWrXK/E6qWfN64Hr8cQlcLu/4cZOUMX26adxWvLhpHPnUUw5RLSPfgpVSygc4APgAXbXWy/JniM7FikPBO3fuZM6cOYSFhXHw4EHc3d3p1q0bgwYNonfv3nh6etp0PKJgnTxpzpRGRpqCu2lpULWqqVXYv7/pWOHg2xciL7Q27UumTYPZs01V+Bo1YPhw81WpktUjvK38DFYvAp8AX2qtx+bghUsAV7TWGTkdrCOzh0PBWmu2bt3K7NmzmTNnDvHx8RQrVgx/f38GDRpE9+7dKVq0qCVjEwXjr79g/nwTuJYsMatC5ctDv34mcLVvD4ULWz1KYZnLl81a8rRpsGKF+RTTsyeMHAndutnVp5r8DFbLgDZAZa31PTv7KKUqA78Du4GeWuuknA3ZsdlLuaWMjAzWrl3L7NmzCQ8P59SpU5QsWZJ+/foRHBxMhw4dcHd3t3qYIh9dvGhqFEZGws8/m99TpUubhLF+/aBTJyhWzOpRCsscOADffGOWCU+eNNPxa7MtOzgrkZ/B6iRwQGvdMhcv/m/gZeBJrfXMnD7OkdlLsMouLS2N5cuXM3v2bKKiojh//jy+vr4EBgYSHBxM69atKWTn69kidy5fNjOtyEgz87pwwXSq6NbNBK6ePUGKpbiolBTzppg61dQmVMq8MUaONG8Mi6bi+RmskoEorXVwLl68KhALzNJaD8np4xyZPQar7JKTk1m8eDGzZ89m3rx5XL58mcqVKxMUFERwcDCNGzdGyU69U0lJMcV1o6NNwe8TJ8zvoyeeMIGrTx+73cYQBe3IEfj2W/N1/LjpMDp8uGkgaeOSKvk9s9qkte6VywHEYPauGufmcY7K3oNVdklJScyfP5/Zs2ezaNEiUlJSqF27NiEhIQQHB1O3bl2rhyjyWUaG2Xv/8UcTvK5trzZvbgJX374Oe0xH5EVamlk7njLFpMBrDZ07w4gR5tOMDbYM8jNY/QkU01rn6q2slPoRaKq1rpybxzkae0iwyItz584RFRVFaGgoK1asICMjg8cee4zg4GAGDRpE1apVrR6iyGdam9qp1wLXtUPI9etfD1yNG0tKvMs5etQkZHz7rfl7hQomaI0cCVWqFNjL5mew+hQYB7TTWq/JxQB+AAZorV0iDc2RZlZ3cvz4cebOnUtYWBgbN24EoE2bNoSEhBAYGIivr6/FIxQFIS7OBK4ffzRnudLTzR58374meLVpI5mFLiU93dQj/OorM+tSCnr3htGjoWPHfD+3lZ/BqgWwFtgKtNFaX83hAH4BGmmtXaISqzMEq+wOHjzI7NmzCQ0NJSYmhsKFC9OlSxeCg4Pp06cPXg50Ql7k3JkzZg8+OtrswV+9ajIL/f1N4OrSRTILXcrhw2aJ8Jtv4PRp09vmuedMead8Kqab3xUsIoAA4Gdg8L0qWSilvIFjwPbcZBE6MmcLVtdordmxYwehoaGEhYURFxeXdYYrJCSEbt264eHhtFW0XFpSksksjI42lX3OnTPFEbp2NYGrVy+QgikuIjnZpJh++SWsXQtFi8KgQWa21bRpnp46p8EKrfU9v4CywB4gHTgCBNzlWgX8kHntKzl5fmf4aty4sXZ26enpes2aNXrMmDG6bNmyGtA+Pj766aef1r/++qtOS0uzeoiigKSkaP3LL1qPHq11pUpag9aFC2vdqZPWn3+udXy81SMUNvPnn1o/+6zWJUqYN0LjxlofOXLfTwds0Tn4HZvjQrZKqYqYmdWjgMakpkcBm4FrPa3qAs8AjYFTQH2t9dkcvYCDc9aZ1Z2kpaWxbNkyQkNDiY6O5uLFi1SoUCErFf7xxx+XVHgnlZEBW7aYGVd0NOzda25//HEz4+rXD+rUsXaMwgYuXDD9tqKjzd7WfWYOFkjVdaWUO/Av4HmgOCZo3XIZcAHopXORkOHoXC1YZXflyhUWLlxIWFgYCxcuJDk5mZo1axIcHExISAj1pb+FU9uz53rg2rzZ3HYtszAgAB57TDILxZ0VaIsQpZQv8DTQHWgEXKuYegKYD0zSWsfm+okdmCsHq+zOnz9PdHQ0oaGhLFu2jIyMDB555BFCQkIYNGgQ1apVs3qIogAdPXo9Jf6338wsrFo1k1kYEGC6WNhRWTphB2zaz0opVQzI0Fon5/nJHJQEq1slJiZmpcKvX78egFatWhEcHMyAAQOkgaSTO33aZBZGRcEvv5g9+rJlr9cs7NgRJDdHWNJ80RU5+qFgWzl8+HBWKvzOnTtxc3OjU6dOhISE0LdvX0qWLGn1EEUBunjRHN2JijJFdy9eNP0Be/Y0gat7d4fqFyjykQQrG5OZVc7t2LGDsLAwwsLCOHLkCEWLFqVXr14EBwfTo0cPaWfi5JKTTT+u6GjTn+vUKTPD6tzZLBX6+4OcP3cdEqxsTIJV7mmt2bBhA2FhYcyZM4eTJ09SsmRJAgICCAkJ4YknnqCwlE5waunp5thOdLSZdcXFmQIJbduawNW3r6mmIZyXBCsbk2CVN2lpaaxYsYLQ0FCioqK4cOEC5cqVy0qFb968uaTCOzmt4Y8/rgeumBhze9Om11Pipcay85FgZWMSrPLP1atX+fnnnwkLC2P+/PkkJydTvXr1rFT4hx56yOohChvYu/d6SvymTea2evWup8Q3aiQp8c5AgpWNSbAqGBcuXCA6OpqwsDB+/fVX0tPTeeihh7JS4WvYuPeOsEZ8/I0p8enppsltQAAMGGBanUgfUcckwcrGJFgVvJMnTxIeHk5YWBhr164FoEWLFgQHBzNw4EDKly9v8QiFLVwrthsVZWoXpqRA5crQv78JXC1bSuByJBKsbEyClW3FxsZmpcJv376dQoUK0bFjR4KDgwkICMDb29vqIQobuHDBBK7wcJMan5xsuh9fC1ytWkngsncSrPJAKTUBeApTTup9rfXMez1GgpV1du3alZUKf+jQITw8POjRowchISH07NmTYtLTwiVcuGCqw4eHm6a3ycmmW/u1pcLWraV6hj2SYHWflFIPA98BLTF1DlcCXfU92qJIsLKe1ppNmzYRGhrKnDlzSExMxMvLi379+hEcHEynTp0kFd5FXLxoDh+Hh5saq1evmsa31wJXmzYSuOyFBKv7pJQaiOmKPCbz+8nACq317Ls9ToKVfUlPT2fFihWEhYURGRnJ+fPnKVu2LAMGDCAkJIQWLVpQSNaHXMKlSzcGritXoHz5GwOXfIaxTk6DlUP936qUaquUmqeUOqaU0kqpYbe5ZrRS6rBS6qpSaqtSqk0uX2YH0F4p5aOUKgU8AVTOh+ELG7pWzunbb78lMTGR6OhonnjiCaZNm0br1q2pUaMGL7/8Mn/++Sfygc25eXpCUBBERMDJkzBnjglQM2ZAhw4mOWPMGFizxhTeFfbJoWZWSqkeQGvgd+B7YLTWeka2+4OAmcBoYE3mn09h+mrFZV6z8w5P311rfTTzmhHAc8BfQDzwu9b6s7uNTWZWjuHixYv89NNPhIaGsnTpUtLT06lfvz7BwcEEBwdTq1Ytq4cobCQpyextzZ1rkjSuXjXVMoKCIDhYWpvYitMvAyqlLgFjbwpWG4HtWusR2W7bD0RorSfe5+t8A0RrrRfe7ToJVo7n1KlTREREEBoaypo1pvXa448/TnBwMEFBQVSsWNHiEQpbuXgR5s2DsDCTDp+WBg8+aDq3BwdL5YyC5HLBSilVBLgMBGutw7Nd9wXwkNa6XS6eu5zW+qRSqg4QDjTSWqfd5rqRwEgAPz+/xrGxLtXCy6nExcUxZ84cQkND2bZtG0op2rRpw4ABA+jfv78ELhdy5ow5wxUWBitXmjJQjz5qgtagQaY/l8g/rhisKgHHMMkRq7Jd9wbwN611jhttK6XWAT5AEvCs1nrrvR4jMyvnsXv3bubMmUN4eDgxMTEopWjdunVW4KpUqZLVQxQ2kpBgEjPCwmDjRnNbixYmcA0caBI1RN64crBqq7Vene26NzGzrQKZyEs/K+cWExNDeHg44eHh7Nq1C6UULVu2zApcVapUsXqIwkYOH4bZs03g2rHDHDZ+4gkICYHAQJCWbPfHFYNVvi0D3g+ZWTm/3bt3ExERQXh4ODt27AC4IXBVlV4WLmPXruuB6+BBKFrUdEAeMgS6dAF3d6tH6DhcLlhl3rYR+FNrPTLbbfuAyPtNsMjBOGRm5YL27t2bNePavn07AM2bN2fAgAEEBgbi5+dn8QiFLWhtlgd/+MEEr7NnoWxZs0w4ZAg0biwZhffilMFKKeUJ1M78dh3wPjAPOKu1jstMXf8Bk7K+FngWeBpooLUu0OwHmVm5rn379mXNuLZt2wZAs2bNsgJXNdmRdwkpKSYVfuZMk1mYkmKyCIcMgb/9TRIz7sRZg1V7YMVt7vpOaz0s85rRwEtARWAnMD57wkVBkWAlAPbv358VuP744w/ApMMHBgYSGBgoLU1cxLlzJjHjhx9gdeYOert2MHiwqZohdZavc8pgZY9kGVDcyYEDB4iMjCQ8PJytW01C6WOPPUZAQAABAQHUr1/f4hEKWzh8GGbNMoFr3z7w8IDevWHoUOjWTUo9SbCyMZlZibs5dOgQ0dHRREVFsW7dOgDq1q2bFbgaNWqEks0Np6Y1bN58fX/r9GlTXHfIEHjqKdMF2RVJsLIxCVYipxISEvjxxx+Jiopi5cqVpKenU61atazA1aJFC9ykJLhTS001+1vTppkiu2lp0KwZDB9uyj250jKhBCsbkWVAkRdnzpxh3rx5REVFsXTpUlJSUihfvjz9+vUjICCA9u3b4y550E4tMdEkZUyfblLiixUzFeGHD4f27Z2/eaQEKxuTmZXIqwsXLrBo0SIiIyP5+eefSUpKolSpUvj7+9O/f386d+4sjSSdmNawZYsJWqGhcP48VK8OTz4Jw4aZvzsjCVY2JsFK5KcrV67wyy+/EBUVxbx58/jrr78oUaIEPXr0ICAggJ49e+Ll5WX1MEUBuXIFfvzRBK5ffzWBrEMHs7cVEADFi1s9wvwjwcrGJFiJgpKamsrKlSuJiooiOjqaxMREihQpQpcuXQgICKB3796UKVPG6mGKAhIXB999Z/pvHTpk9rOGDIERI+CRR6weXd5JsLIR2bMStpSens6GDRuIjIwkKiqK2NhY3NzcaN++PQEBAfTt21cK7TqpjAxYtQq++cY0kkxONkkZI0eapIwSJawe4f2RYGVjMrMStqa15o8//iAqKorIyEj27NkDQIsWLQgICKBfv37STNJJnTljUuCnTIHdu8HLyxw4HjHCNI10JBKsbEyClbBaTExM1lmu33//HYBHH300KyW+QYMGcpbLyWgNa9fC1Kmm4/HVq9CkiZltDRpkgpi9k2BlYxKshD05cuRIVuBau3YtWmseeOCBrMDVpEkTCjl7TrSL+esvkwI/ZQrs3AmenqZ9yciRpqCuvZJgZSOyZyXs3YkTJ/jpp5+Iiopi+fLlpKWlUaVKFfr160e/fv1o06YNhV295o8T0Ro2bDBBa84ck1nYqBE895ypBm9ve1sSrGxMZlbCEfz1118sXLiQqKgoFi9ezJUrVyhTpgx9+vQhICCATp064eHhYfUwRT45d86c2frqKzPb8vY26e/PPQcPPmj16AwJVjYmwUo4mqSkJJYsWUJUVBTz58/nwoULeHl50bNnTwICAujevTuenp5WD1PkA61hzRr48kuTSZiWBp07w+jR0KuXtcV0JVjZmAQr4chSUlJYvnw5UVFR/Pjjj5w6dQoPDw+6du1KQEAA/v7+lC5d2uphinxw4oRJf588GeLjoWpVGDUKnnkGype3/XgkWNmYBCvhLNLT01m7di1RUVFERUVx9OhR3NzceOKJJ7LOclWsWNHqYYo8SkuD+fPhiy9g2TJwd4fAQBgzBlq2tF2HYwlWNibBSjgjrTVbt27NOsu1b98+lFI3nOWqWbOm1cMUebRnj9nXmjEDLlwwlTFGjzZntwo6IUOClY1INqBwFVprdu/enTXjutYJuWHDhjc0lJSzXI4rKck0ivziC9i+HXx8zPLgmDEFV0hXgpWNycxKuJrDhw/f0FBSa82DDz6Y1d6kadOmErgc1LXDxv/7H0RGmu/79IHnn4d27fJ3iVCClY1JsBKu7Pjx41lnuVasWHHDWa6AgABat24tZ7kc1NGjZolwyhRT5umRR0zQCgkxvbfySoKVjUmwEsI4e/YsCxYsICoqiiVLlnD16lV8fX2zznJ17NhRznI5oCtXzJmtzz6DHTugTBlTHWP0aKhS5f6fV4KVjUmwEuJWly5dYvHixURFRbFgwQIuXryIl5cXvXr1IiAggG7duslZLgejNfz2mwla8+aZJcG9e+F+ayZLsLIxCVZC3F1ycjLLli0jKiqKn376idOnT1O0aNEbznKVKlXK6mGKXDh82DSJHD/+/p9DgpWNSbASIufS0tJYs2ZNVoJGfHw87u7udOnShYEDB9K7d298fHysHqawAQlWNibBSoj7o7Vm8+bNhIeHM3fuXOLi4nB3d6dr165Zgcvb29vqYYoCIsHKRuSclRD5R2vNpk2bmDt3LnPnziU+Pp4iRYrQrVs3Bg4ciL+/PyVLlrR6mCIfSbCyMZlZCZG/MjIy2LhxI3PnziU8PJxjx47h4eFB9+7dGThwIL169cLLEboLiruSYGVjEqyEKDgZGRmsX78+K3AdP36cokWL0qNHj6zAVcLeGjWJHJFgZWMSrISwjYyMDNauXcvcuXOJiIjgxIkTFCtWjJ49ezJw4EB69uxJ8eLFrR6myCEJVjYmwUoI20tPT2fNmjXMnTuXyMhIEhMT8fT0pF+/fgQHB9OpUyfc3d2tHqa4CwlWNibBSghrpaen89tvvxEWFkZERATnzp3D19eXgQMHEhwcTMuWLSlUqJDVwxQ3kWBlYxKshLAfycnJLFmyhNDQUObNm8eVK1fw8/MjODiYkJAQHn74YSmyayckWNmYBCsh7NPFixeZN28eoaGhLFmyhPT0dBo0aEBwcDDBwcHSj8tiEqxsTIKVEPbv1KlTREREEBoaypo1awBo3rw5ISEhDBw4kPJW9HV3cTkNVi69gKuUmqeU+kspFXHT7T2UUnuVUvuVUqOtGp8QIn+VLVuW5557jtWrVxMbG8sHH3zAlStXeP7556lUqRJdunRh5syZJCUlWT1UcROXnlkppZ4APIEntdaBmbcVBnYDHYAzwBago9b6+N2eS2ZWQjiumJgYQkNDmTVrFkeOHKFEiRIEBgYydOhQ2rdvL4kZBUhmVjmgtV4BXLzp5seBGK31Ua31ZSAa6GXzwQkhbKZ+/fq8++67HDx4kFWrVhEcHEx0dDQdO3akevXqvPLKK+zevdvqYbo0uwxWSqm2mUt0x5RSWik17DbXjFZKHVZKXVVKbVVKtcmnl68EHM32fTxQOZ+eWwhhxwoVKkSbNm2YOnUqJ06cYPbs2Tz88MP85z//oX79+jRt2pT//e9/nDp1yuqhuhy7DFaYpbmdwAvAlZvvVEoFAZ8Bk4DHgHXAIqWUX7Zrdt7hq+o9Xvt2+ayuu1YqhIsqVqwYQUFBLFy4kPj4eD755BPS09Oz9rf69OlDZGQkycnJVg/VJdj9npVS6hIwVms9I9ttG4HtWusR2W7bD0RorSfm8vnbZz7/tT2rlsA/tdZ9Mr//N3BEaz31No8dCYwE8PPzaxwbG5vLn04I4Wh27NjBDz/8wMyZMzl+/Dg+Pj4EBQUxbNgwmjVrJue3cslp96yUUkWAxsDSm+5aCrTMh5fYBDRQSlVVShUD+gELb3eh1nqK1rqJ1rpJ2bJl8+GlhRD27tqy4NGjR1myZAk9e/bk+++/p0WLFjz00EN8/PHHnDx50uphOh2HC1aAL+AGJN50eyJQITdPpJT6FQgHeiil4pVSLbTWacB4YBmwA/hKa52Q92ELIZyJm5tbVqr7iRMnmDJlCiVLlmTChAlUrlyZgIAAFixYQFpamtVDdQqOGKyuuXn9Ut3mtrs/gdadtNZltdbFtdZVtNbrM2+fr7V+UGtdW2v9v7s9h1LKXyk15fz587kbvRDCaZQsWZIRI0awfv16du3axYsvvsjatWvx9/enatWqvPzyy+zbt8/qYTo0RwxWp4F0bp1FlePW2VaBywxsI6XtthACTBr8hx9+SHx8PNHR0TRt2pSPPvqIOnXq0KZNG6ZPn86lS5esHqbDcbhgpbVOAbYCnW+6qzMmK9CmZGYlhLgdd3d3+vbty7x58zh69Cjvv/8+J0+eZPjw4VSsWJFnnnmGdevWYe9JbvbCLrMBlVKeQO3Mb9cB7wPzgLNa67jM1PUfgNHAWuBZ4GmggdbakpQ8qWAhhLgXrTXr1q1j2rRpzJkzh6SkJOrXr8+oUaMYMmQIpUqVsnqINufQhWwz08lX3Oau77TWwzKvGQ28BFTEnMkar7VeZasx3kyClRAiNy5dusScOXOYPHkymzdvpmjRogwcOJBRo0bRokULl0mBd+hg5UiUUv6Af+3atUfs37/f6uEIIRzQH3/8wZQpU5g1axYXL17koYceYuTIkQwZMgQfHx+rh1egJFjZmMyshBB5denSJcLCwpgyZQpbtmzJqqIxatQopz1w7LSHgoUQwll5enoyYsQINm/ezNatWxkyZAgRERG0aNGChg0b8sUXX+CqyVwSrPJIsgGFEAWhUaNGTJ48mYSEBCZPnkzhwoUZO3YslSpV4plnnuH333+3eog2JcuA+USWAYUQBW3Lli1MnjyZ0NBQLl++TIsWLRgzZgyBgYF4eHhYPbz7IsuAQgjhZJo0acLUqVM5duwYn376KadPn2bw4MH4+fnx2muvcfTo0Xs/iYOSYJVHsgwohLA1Hx8fXnjhBfbs2cOSJUto3rw5kyZNokaNGvTv35/ly5c73WFjWQbMJ7IMKISw0pEjR/j666/55ptvOHPmDPXq1WP06NEMHTqUkiVLWj28O5JlQCGEcCHVq1fn/fffJz4+nhkzZuDp6cm4ceOoXLkyY8aMISYmxuoh5okEKyGEcCJFixblySefZNOmTWzcuJGAgAC+/fZbGjRoQOfOnVm4cCEZGRlWDzPXJFgJIYSTevzxx/nuu++Ij49n0qRJ7N69m4Me4WIAAAY8SURBVF69elG3bl0+//xzh6r+LsEqjyTBQghh73x9fZk4cSKHDx8mLCyM0qVLM27cOKpUqcKECRM4cuSI1UO8JwlWeST9rIQQjsLd3Z1BgwaxYcMG1q9fT7du3fj000+pVasWgYGBrFmzxm6zCCVYCSGEC2revDmzZ8/m8OHDvPTSSyxfvpw2bdrQtGlTZs6cSUpKitVDvIEEKyGEcGFVq1blvffe4+jRo3z99dckJSUxZMgQqlWrxjvvvMOpU6esHiIgwUoIIQRQokQJRo0axa5du1i8eDENGzbkjTfewM/Pj2effZa9e/daOj4JVkIIIbIUKlSIrl27smjRImJiYhgyZAgzZsygbt269OnTh9WrV1uyryXBKo8kG1AI4azq1avHlClTiI2N5Y033mDt2rW0bduW5s2bM3fuXNLS0mw2FglWeSTZgEIIZ1e+fHn+9a9/ERcXx5dffsnZs2cJCgrigQce4LPPPiMpKanAxyDBSgghRI4UL16c5557jj179hAdHU3lypV56aWXuHDhQoG/thSyzSdSyFYI4YoOHz5MjRo17vvxUshWCCFEgctLoMoNCVZCCCHsngQrIYQQdk+ClRBCCLsnwSqP5JyVEEIUPAlWeSTnrIQQouBJsBJCCGH3JFgJIYSwe3IoOJ8opU4BsXe42xu426aWL3A63wdlvXv93I762vnx3Pf7HLl9XG6uz8m18l52rte38r187bE+Wuuy97xSay1fBfwFTLnH/VusHqMVP7ejvnZ+PPf9PkduH5eb63NyrbyXnev1rXwv5/axsgxoG/OtHoBFrPy5C/K18+O57/c5cvu43Fyfk2vlvexcr2/lezlXj5VlQDuglNqic1AbSwh7J+9lUVBkZmUfplg9ACHyibyXRYGQmZUQQgi7JzMrIYQQdk+ClRBCCLsnwUoIIYTdk2AlhBDC7kmwslNKqapKqZVKqRil1J9KqQCrxyTE/VJKrct8H+9USr1h9XiE45FsQDullKoIlNdab1NKlQO2AnW01pctHpoQuaaUKqm1vqCUcgPWAM9prbdZPS7hOGRmZae01sev/c+stT4J/IWpuyaEw9FaX8j8a5HMLyFyRYJVAVFKtVVKzVNKHVNKaaXUsNtcM1opdVgpdVUptVUp1eYOz9UEcAeOFvCwhbhFfr2XlVIbgZPArzKrErklwargeAI7gReAKzffqZQKAj4DJgGPAeuARUopv5uuKwN8DzytZc1WWCNf3sta62ZAZaChUuqhgh60cC6yZ2UDSqlLwFit9Yxst20EtmutR2S7bT8QobWemPm9B/ALMFVr/YNtRy3Ere73vXzTc/wTSNdaf2SDIQsnITMrCyiligCNgaU33bUUaJl5jQJmAMslUAl7lcP3so9Syjfz70WBLsAeW45TOD4JVtbwBdyAxJtuTwQqZP69FRAE9FVKbcv8etiGYxQiJ3LyXi4NLFFKbQe2AL9prRfYbojCGRS2egAu7uY1WHXtNq31GuTDhHAcd3svH8LMvoS4b/LL0BqngXSuf/K8phy3fkIVwp7Je1nYhAQrC2itUzCHfDvfdFdnTCaVEA5B3svCVmQZsIAopTyB2pnfFgL8lFINgbNa6zjgE+AHpdQmYC3wLFAJ+NqK8QpxJ/JeFvZAUtcLiFKqPbDiNnd9p7UelnnNaOAloCLmHMt4rfUqW41RiJyQ97KwBxKshBBC2D3ZsxJCCGH3JFgJIYSwexKshBBC2D0JVkIIIeyeBCshhBB2T4KVEEIIuyfBSgghhN2TYCWEk1NKlc7s8HvM6rEIcb8kWAnh/Jpm/rnF0lEIkQcSrIRwfk0y/5RgJRyWBCshnJ8EK+HwJFgJ4fwkWAmHJ8FKCCemlCoHVAHitNanrB6PEPdLgpUQzu3/t3fHNgjFMBRF7REQ0FCw/05IlBQUSEwQir8BFN88nTNBuisnlmK5gghiBdlcARJBrCCbWBHB54sQrLsftf3ee1xrvfY+D3zLZAWhuvtSW6juQsW/EyvIZbmCGGIFubxXEUOsIJdYEcOCBYTq7mdVnarqsNZ6730e+IXJCgJ197WqzlV1EyoSmKwAGM9kBcB4YgXAeGIFwHhiBcB4YgXAeGIFwHhiBcB4YgXAeGIFwHgfqvbj+gxSTTEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ell, cls_auto_pos, 'r-', label='positions')\n",
    "plt.plot(ell, cls_auto_shape, 'k-', label='shapes')\n",
    "plt.plot(ell, cls_pos_shape, 'b-', label='pos x shape')\n",
    "plt.legend(loc='upper right', prop={'size':12})\n",
    "plt.xlabel(r'$l$', fontsize=22)\n",
    "plt.ylabel(r'$C_l$', fontsize=22)\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.gca().set_yscale('log')\n",
    "plt.gca().set_xscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}