scripts/ABC/WL_Cell_fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Toy model: quadratic function, ressembling weak-lensing power spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from os import environ\n",
    "import numpy as np\n",
    "import pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weak-lensing power spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "basedir = './dlogP_dlogp'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_Cell(path, ell_min, ell_max):\n",
    "    \"\"\"Read P_kappa file and return log_10(ell), log_10(C_ell) within ell-range\n",
    "    \"\"\"\n",
    "    \n",
    "    pk = np.loadtxt(path)\n",
    "    ell  = pk[:,0]\n",
    "    Cell = pk[:,1] / (ell * (ell+1) / (2 * np.pi))\n",
    "\n",
    "    ran     = (ell>=ell_min) & (ell<=ell_max)\n",
    "    \n",
    "    return np.log10(ell[ran]), np.log10(Cell[ran])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fitting range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ell_min = 10\n",
    "ell_max = 1e4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fiducial cosmological parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "Om_fid = 0.306\n",
    "s8_fid = 0.827"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot standard power spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "logell, logCell = read_Cell('{}/P_kappa'.format(basedir), ell_min, ell_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEICAYAAAB4YQKYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAm2klEQVR4nO3dd3hUZf738fc3jZBAQkuooQkiASRC6M2GIBasCCgWFOxYdt3Hturu/h5/u65YEJTqiiIouoodFlilKAECggRpkRpAOgGkJ/fzR8Z92DiBADM5M8nndV1z5cw595z53h7MJ6fdx5xziIiInK0IrwsQEZHSQYEiIiIBoUAREZGAUKCIiEhAKFBERCQgorwuwEvVqlVz9evX97oMEZGwsmjRop3OuaTC88t0oNSvX5/MzEyvyxARCStmtsHffB3yEhGRgFCgiIhIQChQREQkIBQoIiISEAoUEREJCAWKiIgEhAJFREQCokzfh3Km3p2/gYNH8kipUp46leNIqRJHYvlor8sSEfGUAuUMfLgoh+837v2veQmxUaRUiaNO5fKk+EImpUrBdJ3KcZSPifSmWBGREqJAOQMf39eJ3IPH2LTnIJt2H/T9PMSmPQfJ3n6Ab1bt4Mjx/P/6TFLFcqTWTCAtpRJpdSuRVqcSleNjPOqBiEjgKVDOUGJcNIlxiTSvnfibZc45duw/wqY9h8jxhc76XQdZlpPLsDVr+PUhmfWrxpGWUomWKZVIS6lEaq0EykVpT0ZEwpMCJQjMjOSEWJITYmldr/J/LTtw5DjLcnJZsmkvSzbtYd7aXUxZsgWAmMgImtZK4AJfwLSuV5mUKnFedEFE5LRZKD9T3szeB5r43lYC9jrn0vy0Ww/sB/KA48659OKsPz093YXC4JBbcw+xZONelmzay/eb9rIsJ5dDx/IAaJRcge6p1emeWp20OpWIiDCPqxWRss7MFvn7PRvSgXIiMxsK5Drn/uxn2Xog3Tm383TWGSqBUtjxvHxWbztAxtpdzFixjfnrdpOX76hWoRyXNk2me2p1OjWqRmy0Do+JSMkL60AxMwM2Ahc759b4Wb6eUhQoheUePMY3q7cz/cdtzFq1g/1HjlM+OpIujavRPbU6F5+XTNUK5bwuU0TKiHAPlK7AS0UdyjKzdcAewAGjnHOjT7KuwcBggLp167besMHvsP4h6+jxfOav28X0H7cx48dtbMk9TIRB63qVubRpdXq1qKnzLiISVCEbKGY2A6jhZ9FTzrlPfG3eALKdc0OLWEct59wWM0sGpgMPOudmn+q7w2UPpSjOOZZv2VcQLiu2sXzLPgC6NK5Gv7Z1ubRpdWKiNBiCiARWyAbKqZhZFLAZaO2cyylG++eAA865F0/VNtwDpbCcPQf5aPFm3l+4ic17D1GtQgw3tE6hb5sU6leL97o8ESklwjlQegJPOOe6FbE8Hohwzu33TU8H/uycm3qqdZe2QPlVXr5j9podTJq/kZkrt5OX7+jUqCr929aje6r2WkTk7BQVKOFwH0pfYNKJM8ysFjDWOdcLqA58XHDenihgYnHCpDSLjDAuapLMRU2S2bbvMJMXbuK9hZu4f+JiqsbHcEN6Hfq2qUsD7bWISACF/B5KMJXWPRR/8vIdc9bsYNKCjcxYUbDX0vGcqvRrW5cezWpor0VEii1sD3kFU1kKlBNt33eYDxblMGnBRnL2HKJWYiyDuzbkpjZ1NYiliJySAsWPshoov8rPd3yzejtvfPMTC9fvoWp8DAM7N2BAh3okxGo4fhHxT4HiR1kPlBMtWLebEV9nM2v1DiqWi+LWjvUY2KmBbpgUkd9QoPihQPmtrM25vP5NNl9l/Uy5qAj6ta3LoC4NqVWpvNeliUiIUKD4oUApWvb2A4yc9RNTvt+MGVx3QR3uufAcXRkmIgoUfxQop5az5yBjZq/lvYWbOJaXz+UtanLfhefQrNZvnwMjImWDAsUPBUrx7dh/hDe/Xcc78zZw4Mhxrji/Jr+/rIn2WETKIAWKHwqU05d76Bjj5qxl7Nx1HD2ez01tUnjoksYkJ8R6XZqIlBAFih8KlDO3Y/8RXvv3GibO30h0ZAQDO9fn7m7n6HJjkTJAgeKHAuXsrd/5C0Onr+azpVuoFBfN/Rc2YkCHenr4l0gppkDxQ4ESOFmbc/nb1JXMWbOTWomxPNz9XK5vVYdIPbJYpNQpKlA0gJMERPPaibxzZzsm3tWOpIrl+MOHP9Dzldn8a/nPlOU/WkTKEgWKBFTHRtWYcn8nXr+5FXn5jsHvLOKGkfPIXL/b69JEJMgUKBJwZkavFjWZ9khXnr+2BZt2H+SGkfMYMul7tuYe8ro8EQkSBYoETXRkBP3b1eWbxy5kyMWNmLr8Zy5+cRbD/72Gw8fyvC5PRAJMgSJBFxcTxaOXNWHmo93odm4SL/5rNd1fnsXULJ1fESlNFChSYlKqxDFyQGvevasd5aMjuWfCIgaMW8Cabfu9Lk1EAkCBIiWuU6NqfDmkC89dlcoPOXvp+eoc/vTZcnIPHfO6NBE5CyEdKGaWZmYZZrbEzDLNrG0R7Xqa2Sozyzazx0u6Tjl9UZER3N6pAd88dhE3tUnhre/Wc9GL3zBx/kby8nUYTCQchXSgAC8Af3LOpQHP+N7/FzOLBEYAlwOpQD8zSy3JIuXMVYmP4flrW/DZA505JymeJz9extXD57JQlxmLhJ1QDxQHJPimE4Etftq0BbKdc2udc0eB94DeJVSfBEjz2olMvrsDw/pdwK4DR7lx5Dx+N3kpew8e9bo0ESmmKK8LOIWHgWlm9iIF4dfRT5vawKYT3ucA7YpaoZkNBgYD1K1bN2CFytkzM65uWYtLmyYz/N/ZjJq9llmrt/Onq5vTq0UNzDSMi0go83wPxcxmmFmWn1dv4F7gEedcCvAIMM7fKvzMK/IgvHNutHMu3TmXnpSUFJhOSEDFxUTxh57n8ekDnaiRGMv9Excz+J1FbNt32OvSROQkQnpwSDPLBSo555wV/Hma65xLKNSmA/Ccc66H7/0TAM65/z3V+jU4ZOg7npfPuLnreGn6amKiIniyV1P6tknR3oqIh8J1cMgtQDff9MXAGj9tFgKNzayBmcUAfYFPS6g+CbKoyAju7nYOUx/uSrNaCTzx0TL6j5nP+p2/eF2aiBQS6oEyCBhqZkuB5/Gd+zCzWmb2JYBz7jjwADANWAFMds4t96heCZIG1eKZeFd7/ve6FmRtzqXHK7MZNesnjufle12aiPiE9CGvYNMhr/D0c+5hnp6SxYwV22hRO5G/XX8+qbUSTv1BEQmIcD3kJfIbNRJjGXNra4b3v4CtuYe4evhc/j5tpQacFPGYAkXCkplx5fm1mP5IN3qn1WbE1z/Ra9gclm7a63VpImWWAkXCWuX4GIb2acn4gW05dDSP69/4jhFfZ2v4FhEPKFCkVOh2bhJTH+pKj2Y1+Pu0VfQbk8HmvXqYl0hJUqBIqZEYF83w/hfw4o0tWb45l56vzOazpf5G6xGRYFCgSKliZtzQug5fPtSFRskVeHDS9zw6eQn7D2tofJFgU6BIqVSvajwf3N2Bhy5pzJTvN9Nr2BwWbdAIxiLBpECRUisqMoJHup/L5Ls74BzcOHIeL09frZshRYJEgSKlXnr9Knz1UBeuSavNqzPX0GfUPDbuOuh1WSKljgJFyoSKsdG8dFMar/ZNY832A/QaNod/LsqhLI8UIRJoChQpU3qn1earh7qQWiuB332wlAcnfa9n2YsEiAJFypw6leOYNKg9j/VowldZP9PrVZ2wFwkEBYqUSZERxv0XNeKDezoQEQF9RmXw2sw1usNe5CwoUKRMa1W3Ml8M6cIVLWoydPpq+o/JYGuu7rAXORMKFCnzEmKjebVvGi/e2JJlm3O5/NU5TFv+s9dliYQdBYoI//8O+88f7EydyuW5+51FPD1lmYbEFzkNChSREzRMqsBH93ZiUJcGTMjYyNXD57Lq5/1elyUSFhQoIoXEREXw1BWpjB/Ylt2/HOPq4XN5J2OD7lkROQUFikgRup2bxFcPdaF9w6r8cUoWd7+ziD2/HPW6LJGQFdKBYmZpZpZhZkvMLNPM2hbRbr2ZLfu1XUnXKaVXUsVy/OP2Njx9RVO+XrWdy1+dw7yfdnldlkhICulAAV4A/uScSwOe8b0vykXOuTTnXHqJVCZlRkSEcVeXhnx8XyfiYiLpPzaDv361kqPHNcikyIlCPVAckOCbTgT0tCTxTPPaiXw+pDN926QwctZPXPfGt/y044DXZYmEDAvlE41m1hSYBhgF4dfRObfBT7t1wB4KAmiUc270SdY5GBgMULdu3dYbNvxmdSKnNDXrZx7/6AeOHMvnmatS6dsmBTPzuiyREmFmi/wdDfI8UMxsBlDDz6KngEuAWc65f5pZH2Cwc+5SP+uo5ZzbYmbJwHTgQefc7FN9d3p6usvM1CkXOTPb9h3md5OXMjd7J5elVuev159PlfgYr8sSCbqQDZSTMbNcoJJzzlnBn3+5zrmEU3zmOeCAc+7FU61fgSJnKz/f8ea363hh6ioqxUUztE9LujRO8roskaAqKlBC/RzKFqCbb/piYE3hBmYWb2YVf50GLgOySqxCKdN+PWE/5f5OJJaPZsC4Bfzl8x91h72USaEeKIOAoWa2FHge37kPM6tlZl/62lQH5vraLAC+cM5N9aRaKbNSayXw2YOdua1DPcbNXcc1I75l9TbdYS9lS0gf8go2HfKSYPj3ym384cMf2H/4OE/2asqtHerphL2UKuF6yEsk7Fx8XnW+eqgrHc6pyrOfLmfgWwvZvv+w12WJBJ0CRSQIfr3D/k9XN+O7n3bR4+XZTM3a6nVZIkGlQBEJEjPjto71+WJIZ+pUjuOeCYt5dPIS9h3WM+yldFKgiARZo+SKfHRfR4Zc3IhPlmzh8lc0HpiUTgoUkRIQHRnBo5c14cN7OhATFUG/MRn8jy4vllJGgSJSgi6oW5kvhnRmQPt6jJ27jquHzyVrc67XZYkEhAJFpITFxUTxl2ua89Ydbdh78BjXvv4tI77O5nieRi+W8KZAEfHIhU2SmfZwVy5rVoO/T1tFn1Hz2LDrF6/LEjljChQRD1WOj2F4vwt4tW8a2dsPcPmrc5g4f6MeNyxhSYEi4jEzo3dabaY90pVWdSvz5MfLGPjWQrbmHvK6NJHTokARCRE1E8vz9sC2PHdVKvPW7uKyl2bz/kLtrUj4UKCIhJCICOP2Tg2Y9nBXUmsl8H/+uYxb31xAzp6DXpcmckoKFJEQVK9qPJMGtecv1zRn8YY99Hh5Nu9kbCA/X3srEroUKCIhKiLCGNC+XsG5lXqV+eOULPqPzdCVYBKyFCgiIa5O5TjeHtiWF64/n+Vb9tHjldmMm7uOPO2tSIhRoIiEATOjT5sUpj/SjY7nVOMvn/9In1HzyN5+wOvSRP5DgSISRmokxjLutnRe6tOS7O0H6DVsDm9885PuspeQoEARCTNmxnWt6jD90a5c1CSJv01dyXVvfMfKn/d5XZqUcSEdKGbW0szmmdkyM/vMzBKKaNfTzFaZWbaZPV7SdYp4IbliLCNvac3w/heQs+cQVw6by9+nrdQIxuKZkA4UYCzwuHOuBfAx8FjhBmYWCYwALgdSgX5mllqiVYp4xMy48vxazHi0G73TajPi65/o+cpsvs3e6XVpUgaFeqA0AWb7pqcD1/tp0xbIds6tdc4dBd4DepdQfSIhoUp8DEP7tOTdu9oBcPPY+Tw6eQm7fznqcWVSloR6oGQBV/umbwRS/LSpDWw64X2Ob55fZjbYzDLNLHPHjh0BK1QkFHRqVI2pD3flgYsa8emSLVwy9Bs+Wpyj4VukRHgeKGY2w8yy/Lx6AwOB+81sEVAR8PfnlvmZV+T/Pc650c65dOdcelJSUmA6IRJCYqMj+X2PJnwxpAsNqsXz6OSl3DJuPut36oZICa4orwtwzl16iiaXAZjZucAVfpbn8N97LnWALYGpTiR8NalRkQ/v6ci7Czbywlcr6fHKbIZc0pjBXRsSHen535JSCoX0vyozS/b9jACeBkb6abYQaGxmDcwsBugLfFpyVYqErl+Hb5nxu25cfF4yf5+2iqtem8vijXu8Lk1KoZAOFAqu2FoNrKRgr+MfAGZWy8y+BHDOHQceAKYBK4DJzrnlHtUrEpKqJ8Tyxi2tGXNrOrmHjnH9G9/xzCdZ7Dt8zOvSpBSxsnyyLj093WVmZnpdhkiJOnDkOEP/tYq3vltPUoVyPHtVM3q1qIGZv9ORIr9lZoucc+mF54f6HoqIBFiFclE8e1UzPrm/E8kJ5bh/4mLueGshG3fpmStydhQoImXU+XUqMeW+TjxzZSoL1+2m+8uzGPF1NkePa1wwOTOnFShm1tnM7jezhifMaxD4skSkJERFRjCwcwNm/u7C/5y0v2LYHBas2+11aRKGTncPJYmCO9OfNbOXzCwdeD7wZYlISaqRWHDS/s3b0zl4NI8+o+bx2AdLdae9nJZTnpQ3sz8CB51zQ33vI4HWQHMgEajknHs22IUGg07Ki/zWwaPHGTYzm7Fz1lIxNoonezXlhtZ1dNJe/uNsTsoPAN749Y1zLs85t8D32dhwDRMR8S8uJorHLz+Pz4d0pmFSBR778AduGp1B9vb9XpcmIa44gXLIOefv8o+3gVsCXI+IhIjzaiTwwd0d+N/rWrDq5/1c/uocXpy2SsPjS5GKFShmVrPwTN/IvscDX5KIhIqICKNf27rM/F03rjq/FsO/ztbw+FKk4gTKUOATM6t34kzfsCi6vlCkDKhWoRwv3ZTGhDvb4dDw+OLfKQeHdM59YGZxwCIzywCWUBBENwLPBbU6EQkpnRtXY9rDXXnt32sYNWstX6/czlNXpHJ9q9o6aS/Fu2zYOTceaABMBqKBw0A/59y7QaxNREJQbHQkj/U4jy+GdKFhUgV+/8FSbh47n3UaHr/M01heumxY5Izl5zsmLtjI36au5MjxfB68qBF3dzuHmCgNwlGaaSwvEQm4iAjjlvb1mPloN7o3rc7Q6au5YtgcMtfrTvuySIEiImctOSGWETe3+s+d9jeMnMcTHy0j95CGxy9LFCgiEjAXn1edfz3Slbs6N+D9hRu5ZOgsvly21euypIQoUEQkoOLLRfH0lal8+kBnaibGct+7i7l3wiJ27D/idWkSZAoUEQmK5rUT+fi+jvyhZxNmrtxO95dnMeX7zZTlC4FKOwWKiARNVGQE913YiC+HdKZBtXgefn8Jg97O5Ofcw16XJkGgQBGRoGuUXJEP7+nI01c0ZW72Trq/PIvJmZu0t1LKhHSgmFlLM5tnZsvM7DMzSyii3XpfmyVmphtLREJQZIRxV5eGTH2oK01rJvCHD3/g1jcXsHnvIa9LkwAJ6UABxgKPO+daAB8Dj52k7UXOuTR/N9uISOioXy2e9wa158+9m7Fowx4ue2kWEzI2kJ+vvZVwF+qB0gSY7ZueDlzvYS0iEiAREcatHeoz7eGuXFC3Mk9PyeLmsfPZuMvfkzIkXIR6oGQBV/umbwRSimjngH+Z2SIzG3yyFZrZYDPLNLPMHTt2BLBUETldKVXieOfOtvz1uhZkbc6lxyuzeXPuOu2thCnPx/IysxlADT+LngJWAcOAqsCnwBDnXFU/66jlnNviG1J/OvCgc2524XaFaSwvkdCxZe8hnvx4Gd+s2kGHhlUZ2qcltSqV97os8aOosbw8D5TiMrNzgQnOubanaPcccMA59+Kp1qlAEQktzjkmZ27iT5/9SFSE8X+vbcFVLWt5XZYUEpaDQ/r2ODCzCOBpYKSfNvFmVvHXaeAyCg6ViUiYMTNualOXL31D4z846XsefX8J+w5rTLBwENKBAvQzs9XASmAL8A8oOMRlZl/62lQH5prZUmAB8IVzbqon1YpIQNSvFs+H93TgoUsa88nSLVz+yhwWrNMIxqEubA55BYMOeYmEvkUb9vDI+0vI2XOQey88h4cuOVfPW/FYWB7yEhFpXa8yXz7UhRtbpzDi65+4/o3vyN5+wOuyxA8FioiEvArlovjbDecz8pZWbNpzkCtfm8OEjA0auiXEKFBEJGz0bF6TaQ93pU39Kjw9JYs7x2dqWPwQokARkbBSPSGW8Xe05dmrUpmbvZOer8xm5optXpclKFBEJAxFRBh3dGrA5w92JqliOe4cn8nzX67gWF6+16WVaQoUEQlb51avyCcPdGJA+3qMnr2W/mMy9KwVDylQRCSslYuK5C/XNOfVvmks37KPK4bNYe6anV6XVSYpUESkVOidVptPH+hElfgYBrw5n1dnrNEgkyVMgSIipUaj5IJDYNek1eblGau57R8L2HVAV4GVFAWKiJQqcTFRvNSnJc9f24L563ZzxbC5LNqgYVtKggJFREodM6N/u7p8dG9HYqIiuGlUBmPnrNWNkEGmQBGRUqt57UQ+e7AzF5+XzP98sYJ7JyzWyMVBpEARkVItsXw0owa05qleTZm+YhtXvTaX5VtyvS6rVFKgiEipZ2YM6tqQ9we358ixfK59/TsmZ27yuqxSR4EiImVGev0qfDGkM23qV+YPH/7A/3z+I3m6tDhgFCgiUqZUrVCO8Xe05faO9Rk7dx0D31qo8yoBokARkTInKjKC565uxvPXtuDb7J1cO+Jb1u/8xeuywp4CRUTKrP7t6vLOne3Y9ctReo/4lu+yNWTL2QiJQDGzG81suZnlm1l6oWVPmFm2ma0ysx5FfL6KmU03szW+n5VLpnIRCXcdzqnKp/d3JrliOQa8uYB3MjZ4XVLYColAAbKA64DZJ840s1SgL9AM6Am8bmaRfj7/ODDTOdcYmOl7LyJSLHWrxvHRfR3p2rgaf5ySxTOfZGko/DMQEoHinFvhnFvlZ1Fv4D3n3BHn3DogG2hbRLvxvunxwDVBKVRESq2KsdGMva0Ng7s25O15G7j9HwvYe/Co12WFlZAIlJOoDZx4sXiOb15h1Z1zWwF8P5OLWqGZDTazTDPL3LFjR0CLFZHwFhlhPNmrKX+/4XwWrtvDNSO+JXv7Aa/LChslFihmNsPMsvy8ep/sY37mndVF48650c65dOdcelJS0tmsSkRKqRvTU5g4qB0Hjhzn2te/5ZtV270uKSyUWKA45y51zjX38/rkJB/LAVJOeF8H2OKn3TYzqwng+6mtLyJnJb1+Fabc34nalcoz8K2FvDl3nQaXPIVQP+T1KdDXzMqZWQOgMbCgiHa3+aZvA04WUiIixVKnchz/vLcjlzatzp8//5FnP12uO+tPIiQCxcyuNbMcoAPwhZlNA3DOLQcmAz8CU4H7nXN5vs+MPeES478C3c1sDdDd915E5KzFl4ti5C2tGdSlAW/P28ADExdz+Fie12WFJCvLu3Dp6ekuMzPT6zJEJEyMnbOW//liBW0bVGHMrekklo/2uiRPmNki51x64fkhsYciIhIO7urSkFf7pvH9xj30GTmPrbmHvC4ppChQREROQ++02rx1R1s27z3E9a9/x5pt+70uKWQoUERETlOnRtV4b3B7juY5bhg5j8z1emY9KFBERM5I89qJfHxfR6rEx3Dz2PlMW/6z1yV5ToEiInKGUqrE8eE9HTivZgL3TljEu/PL9sCSChQRkbNQtUI5Jg1qR7dzk3jq4yxemr66zN4AqUARETlLcTFRjL41nRtb12HYzDU88dEyjpfB0YqjvC5ARKQ0iI6M4IUbzqd6QizDv85m54EjvNavFeVj/D1xo3TSHoqISICYGb/v0YS/9G7GzJXbuXlsBrmHys7z6hUoIiIBNqBDfV7v34plm3PpPyaD3b+UjeeqKFBERILg8hY1GX1rOtnbD9BvdAbb9x/2uqSgU6CIiATJRU2S+cftbdi4+yB9R2WU+qFaFCgiIkHUsVE13r6zLdv3H6HPqHls2n3Q65KCRoEiIhJkbepXYcJd7cg9eIw+o+axbucvXpcUFAoUEZESkJZSiUmD23PkeD59Rs0rlYNKKlBEREpIs1qJvD+4PQA3jc5g+ZZcjysKLAWKiEgJaly9IpPv7kBsVAT9RmewdNNer0sKGAWKiEgJa1Atnvfv7kBiXDQ3j53PwlIy/H1IBIqZ3Whmy80s/4TnxP+67AkzyzazVWbWo4jPP2dmm81sie/Vq2QqFxE5MylV4ph8dweSK5bj1nEL+C57p9clnbWQCBQgC7gOmH3iTDNLBfoCzYCewOtmVtTAOC8759J8ry+DWq2ISADUTCzPe3e3J6VKee54ayFfr9rudUlnJSQCxTm3wjm3ys+i3sB7zrkjzrl1QDbQtmSrExEJnuSKsbw3uAONkisw+O3MsH5QV0gEyknUBjad8D7HN8+fB8zsBzN708wqB780EZHAqBIfw8S72tOsViIPTFzMzBXbvC7pjJRYoJjZDDPL8vPqfbKP+Znn78k1bwDnAGnAVmDoSeoYbGaZZpa5Y8eO0+mCiEjQJMZFM35gW5rWTODeCYuZtTr8fj+VWKA45y51zjX38/rkJB/LAVJOeF8H2OJn3ducc3nOuXxgDCc5LOacG+2cS3fOpSclJZ1pd0REAi6xfDRvD2z7n8Nf34bZifpQP+T1KdDXzMqZWQOgMbCgcCMzq3nC22spOMkvIhJ2KsXFMOGudtSvGs+d4xeSsXaX1yUVW0gEiplda2Y5QAfgCzObBuCcWw5MBn4EpgL3O+fyfJ8Ze8Ilxi+Y2TIz+wG4CHikxDshIhIgVeILQqV2pfIMfGshmWFyn4o55++URNmQnp7uMjMzvS5DRMSv7fsOc9PoDHbsP8KEu9qRllLJ65IAMLNFzrn0wvNDYg9FRER+KzkhlomD2lElPoZbx80na3Noj/2lQBERCWE1E8szcVA7KsZGc8u4+azYus/rkoqkQBERCXF1KscxcVA7YqMiuXnsfFaH6ND3ChQRkTBQr2o8Ewe1IzLC6D9mPj/tOOB1Sb+hQBERCRMNkyowaVA7nHP0H5PB+hB78qMCRUQkjDRKrsi7g9px9Hg+/cdkhNQz6hUoIiJh5rwaCbxzZzsOHDlOvzEZbNl7yOuSAAWKiEhYal47kXfubEfuwWMMGDef3b8c9bokBYqISLhqmVKJsbelk7PnELf/YwH7Dx/ztB4FiohIGGvXsCqv39yK5Vv2MfjtRRw+ludZLQoUEZEwd0nT6gy9sSXz1u7iwUnfczwv35M6FCgiIqXANRfU5rmrUpn+4zYe/2gZ+fklP05jVIl/o4iIBMXtnRqw99AxXpmxhsTy0Tx9RVPM/D2nMDgUKCIipchDlzRm78FjjJu7jspx0TxwceMS+24FiohIKWJmPHNlKrmHjvHiv1aTGBfDgPb1SuS7FSgiIqVMRITxwg3ns+/QMZ75JIvE8tFc3bJW8L836N8gIiIlLjoyghE3t6JN/So8+v4Svlm1PejfqUARESmlYqMjGXtbOk1qVOSeCYuC/ihhBYqISCmWEBvN+IFtqZVY8Hz6YD6gKyQCxcxuNLPlZpZvZuknzK9qZl+b2QEzG36Sz1cxs+lmtsb3s3LJVC4iEvqqVSjH23e2Jb5cFAPGLWDDruAMex8SgQJkAdcBswvNPwz8Efj9KT7/ODDTOdcYmOl7LyIiPnUqx/HOnW3Jy8/nlnHz2bbvcMC/IyQCxTm3wjm3ys/8X5xzcykIlpPpDYz3TY8HrglshSIi4a9RckXGD2zLuckViYuJDPj6S8tlw9Wdc1sBnHNbzSy5qIZmNhgYDFC3bt0SKk9EJDScX6cS425vE5R1l1igmNkMoIafRU855z4pqTqcc6OB0QDp6eklP9iNiEgpVWKB4py7NIir32ZmNX17JzWB4F9wLSIi/yUkzqEEwKfAbb7p24AS2+MREZECIREoZnatmeUAHYAvzGzaCcvWAy8Bt5tZjpml+uaPPeES478C3c1sDdDd915EREpQSJyUd859DHxcxLL6Rcy/64TpXcAlQSlORESKJST2UEREJPwpUEREJCAUKCIiEhDmXNm9FcPMdgAbTtGsGrCzBMopCaWpL6D+hDr1J3SdbV/qOeeSCs8s04FSHGaW6ZxLP3XL0Fea+gLqT6hTf0JXsPqiQ14iIhIQChQREQkIBcqpjfa6gAAqTX0B9SfUqT+hKyh90TkUEREJCO2hiIhIQChQREQkIBQogJm9aWbbzSyriOVmZsPMLNvMfjCzViVd4+koRn8uNLNcM1viez1T0jUWl5mlmNnXZrbCzJab2UN+2oTN9ilmf8Ji+5hZrJktMLOlvr78yU+bcNo2xelPWGybE5lZpJl9b2af+1kW2O3jnCvzL6Ar0ArIKmJ5L+ArwID2wHyvaz7L/lwIfO51ncXsS02glW+6IrAaSA3X7VPM/oTF9vH9967gm44G5gPtw3jbFKc/YbFtCtX8KDDRX92B3j7aQwGcc7OB3Sdp0ht42xXIACr5HuQVkorRn7DhnNvqnFvsm94PrABqF2oWNtunmP0JC77/3gd8b6N9r8JX+YTTtilOf8KKmdUBrgDGFtEkoNtHgVI8tYFNJ7zPIUx/CZygg2/X/isza+Z1McVhZvWBCyj4y/FEYbl9TtIfCJPt4zucsoSCp6ROd86F9bYpRn8gTLaNzyvAH4D8IpYHdPsoUIrH/MwL579cFlMwFk9L4DVgirflnJqZVQD+CTzsnNtXeLGfj4T09jlFf8Jm+zjn8pxzaUAdoK2ZNS/UJKy2TTH6EzbbxsyuBLY75xadrJmfeWe8fRQoxZMDpJzwvg6wxaNazppzbt+vu/bOuS+BaDOr5nFZRTKzaAp++b7rnPvIT5Ow2j6n6k+4bR8A59xe4BugZ6FFYbVtflVUf8Js23QCrvY99fY94GIzm1CoTUC3jwKleD4FbvVdEdEeyHXObfW6qDNlZjXMzHzTbSn4d7DL26r889U5DljhnHupiGZhs32K059w2T5mlmRmlXzT5YFLgZWFmoXTtjllf8Jl2wA4555wztVxBU+97Qv82zl3S6FmAd0+IfEIYK+Z2SQKrt6oZgXPtn+WghNyOOdGAl9ScDVENnAQuMObSounGP25AbjXzI4Dh4C+znfJRwjqBAwAlvmObQM8CdSFsNw+xelPuGyfmsB4M4uk4BfrZOfc52Z2D4TltilOf8Jl2xQpmNtHQ6+IiEhA6JCXiIgEhAJFREQCQoEiIiIBoUAREZGAUKCIiEhAKFBERCQgFCgiIhIQChSREGNm8WY23HfnskjYUKCIhJ57gHJAZ68LETkdChSR0NOTggdvLfG4DpHTokARCSFmFgtEUvDEzVkelyNyWhQoIqGlMQWBstI5d8zrYkROh0YbFgktScC5FDyaVSSsaA9FJLTUouDhWxFmVtnrYkROhwJFJESYWRQF505qACOBPG8rEjk9eh6KiIgEhPZQREQkIBQoIiISEAoUEREJCAWKiIgEhAJFREQCQoEiIiIBoUAREZGA+H9I4CJxUofj8gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = plt.plot(logell, logCell)\n",
    "tx = plt.xlabel('$\\ell$')\n",
    "ty = plt.ylabel('$C_\\ell$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting functions\n",
    "We define a quadratic function $q(u) = c + a (u - u_0)^2$, which we will fit to the (ten-)logarithm of the weak-lensing power spectrum $C_\\ell$. We will then vary two parameters, that mimic cosmological parameter as follows:\n",
    "\n",
    "1. A tilt parameter `tilt`, corresponding to the matter density $\\Omega_{\\rm m}$. This is linearly related to the shift $u_0$.\n",
    "\n",
    "2. An amplitude parameter `ampl`, corresponding to the power-spectrum amplitude $\\sigma_8$. The logarithm of that parameter is linearly related to the offset $c$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add(ampl):\n",
    "    \"\"\"Return the additive constant of the quadratic function\n",
    "       from the amplitude fitting parameter\n",
    "       (mimics power-spectrum normalisation s8)\n",
    "    \"\"\"\n",
    "    \n",
    "    # This provides a best-fit amp=0.827, but the 10% increased\n",
    "    # spectrum (0.9097) gives a best-fit of 0.925\n",
    "    # Changing the prefactor of amp or lg(amp) does not help...\n",
    "    \n",
    "    c = np.log10(ampl)*2 - 6.11568527 + 0.1649\n",
    "    \n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def shift(tilt):\n",
    "    \"\"\"Return the shift parameter of the quadratic function\n",
    "       from the tilt parameter (mimics matter density)\n",
    "    \"\"\"\n",
    "    \n",
    "    u0 = tilt * 1.85132114 / 0.306\n",
    "    \n",
    "    return u0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadratic(u, *params):\n",
    "    \"\"\"Used to fit quadratic function varying all three parameters\n",
    "    \"\"\"\n",
    "    \n",
    "    (ampl, tilt, a) = np.array(params)\n",
    "    c  = add(ampl)\n",
    "    u0 = shift(tilt)\n",
    "    \n",
    "    return c + a * (u - u0)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadratic_c(u, *params):\n",
    "    \"\"\"Used to fit quadratic function varying constant coefficient only\n",
    "    \"\"\"\n",
    "    (ampl) = np.array(params)\n",
    "    c  = add(ampl)\n",
    "    u0 = shift(tilt_glob)\n",
    "    \n",
    "    return c + a_glob * (u - u0)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadratic_u0(u, *params):\n",
    "    \"\"\"Used to fit quadratic function varying shift parameter only\n",
    "    \"\"\"\n",
    "    (tilt) = np.array(params)\n",
    "    c = add(ampl_glob)\n",
    "    u0 = shift(tilt)\n",
    "    \n",
    "    return c + a_glob * (u - u0)**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit standard power spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "p0 = [1, 0, 0]\n",
    "popt, pcov = curve_fit(quadratic, logell, logCell + logell, p0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Best-fit parameter (ampl, tilt, $u_0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.82708472,  0.306     , -0.17586216])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAApf0lEQVR4nO3deXhU5dnH8e+dhBAQEBRQEGncKEUqW1RQRFxQi1RFsUAFpVh9kUKtu1hR3HCjYCsutWqhKIKKIm61WkFfCypBUVGEur1KcQGsiAiy3e8fTwIhJGEmmZkzk/l9rutckznnmZn7cMLcOc9q7o6IiEiscqIOQEREMosSh4iIxEWJQ0RE4qLEISIicVHiEBGRuORFHUAqNG3a1AsLC6MOQ0QkYyxYsGCluzer6FhWJI7CwkKKi4ujDkNEJGOY2f9VdkxVVSIiEhclDhERiYsSh4iIxCUr2jhEJH1t3LiRZcuWsX79+qhDyUoFBQW0atWKOnXqxPwaJQ4RidSyZcto2LAhhYWFmFnU4WQVd2fVqlUsW7aMffbZJ+bXqapKRCK1fv16dt99dyWNCJgZu+++e9x3e0ocUrF58+DGG8NjMl8jAkoaEarOv72qqmRH8+bBMcfAhg2Qnw///Cd065b414hIRtIdh+xozpyQADZvDo9z5iTnNSJpIjc3l44dO3LggQfSoUMHxo8fz5YtW6p8zSeffMLUqVNTFGF6UeLIFvFUI/XsGe4acnPDY8+eyXlNvHGJJEm9evVYuHAh7777Ls8//zzPPPMM11xzTZWvyebEgbvX+q1Lly6e1ebOda9Xzz03NzzOnRvba8aOja1sdV9Tnbik1nnvvffif1F1fj+rsMsuu2z3/MMPP/TddtvNt2zZ4h9//LF3797dO3Xq5J06dfJ//etf7u5+6KGHeqNGjbxDhw4+fvz4SstlgoquAVDslXynRv6lnoot6xPH2LHhyxnC49ixUUcUpGtcklJxJ44k/MFRPnG4uzdu3Ni/+OILX7t2ra9bt87d3ZcuXeql3yezZ8/2E088cWv5ysplgngThxrHs0FpNVJpw3Ws1UjJlq5xSXqrqD0tCR0xwndnGKA4YsQIFi5cSG5uLkuXLq2wfKzlagMljmzQrVvo5TRnTvhyTpfeTukal6S3FPzB8dFHH5Gbm0vz5s255ppr2GOPPXjrrbfYsmULBQUFFb5mwoQJMZWrDZQ4skW3bun5xZyucUn6SvIfHCtWrGDYsGGMGDECM2P16tW0atWKnJwcJk+ezObNmwFo2LAha9as2fq6ysrVRkocmWjevOz9Kz2bz122SfAfHOvWraNjx45s3LiRvLw8Bg8ezIUXXgjA8OHDOe2003jkkUc46qij2GWXXQA46KCDyMvLo0OHDgwZMqTScrWRldbj1WZFRUVeaxZyyuaBdtl87rXY4sWL+clPfhJ1GFmtomtgZgvcvaii8hrHkWmyeaBdNp+7SBpR4sg01R1oVxtk87mLpBG1cWSadOyJtGkTrFwJX30FP/wAZpVvOTnbP8/NhaZNoUmT8Lwq6XjuIllIiSMTpaIn0vr18MknIRmUbitWbP+8dN+qVTX/vPx82GMP2HPPqrcOHZQwRCKmxJHtNm+Gjz6Cd94J26JF4fHf/4aKJnnbfXdo3hyaNYP27cPPpVuzZlBQAGEseOXbli3bfi69W/nii23bp5/C66+HxFRR543ddoOf/jQkkdLtwAPDZ4tI0ilxZAv38KVcNjm88w689x6sWxfKmMF++4Uv5V/8An784/BXfrNmITE0bQp5KfyVqSipfPEFfPwxvPUW3HsvfP99KJubG+Lt0AEOOmhbQmnRYudVYCISFyWO2sodPvwQXnghbC+/HKqVSu2xR0gQw4aFx/btoV07SKe+53l526qoKrJlSzjHt97atv3rX/DQQ9vKNG26rXqrVy/o2jVUi4mkwJAhQ+jTpw/9+vWL+TUzZ86kTZs2tGvXDoCrrrqKHj16cOyxx8b0+ksuuYRnnnmG3r17s99++1G/fn3OPPNMJk2axHHHHUfLli2rdS5lKXHUJl9+CS++uC1ZfPpp2N+qFfTuDZ07b0sSzZpFG2si5OTAAQeErex/zG++gbff3j6ZjB0L118PDRrAUUfBcceF7YADdEciKbd582Zyc3MrPDZz5kz69OmzNXFce+21cb33n//8Z1asWEHdunW32z9p0iTat2+vxJH11qwJdxIvvBB6G73zTtjfuDEcfTRcdhkce2z2fTk2bgw9eoRt3jwoXTMhPz/cdbz9Njz5ZNjXuvW2JHLMMaH9RLLODTfcwN/+9jf23ntvmjVrRpcuXbj44ovp2bMn48aNo6ioiJUrV1JUVMQnn3zCJ598wuDBg1m7di0AEydO5LDDDsPdGTlyJC+++CL77LMPZQdYFxYWMnToUP7xj38wYsQI1qxZwz333MOGDRvYf//9mTJlCgsXLmTWrFm89NJLXH/99cyYMYPrrrtu613L/PnzOf/881m7di1169bln//8Jw0bNtz6GSeddBJr167l0EMPZdSoUSxevJgGDRpQWFhIcXExZ5xxBvXq1WPevHnUq1ev2v9eShyZ5sMPQ1XMc8/Bq6+GdoC6deGII+CXvwyJolOnUOcv2wYNbtkSOgIcfDA89ljoEPD88/CPf8Ajj4T2EjMoKtqWSA4/XP+Oqfa738HChYl9z44d4bbbKj28YMECpk2bxptvvsmmTZvo3LkzXbp0qfItmzdvzvPPP09BQQH//ve/GThwIMXFxTz++OMsWbKEd955hy+//JJ27doxdOjQra8rKCjglVdeAWDVqlWcc845AFx55ZXcd999jBw5kpNOOqnC6q0NGzbQv39/pk+fzsEHH8y33367w5f/rFmzaNCgAQtL/g3HjBkDQL9+/Zg4ceLWJFhTShxRi2XupS+/hOnTw1/Or70W9hUVwSWXhL+SDzsMavDXQ61W2Uyq++4L//M/Ydu0CebP35ZIbroJbrghNKwPHAiDB4d2kmy6a8si//u//0vfvn2pX78+EP5q35nKplB/+eWXGThwILm5ubRs2ZKjjz56u9f1799/68+LFi3iyiuv5JtvvuG7777j+OOPr/IzlyxZQosWLTj44IMBaNSoUVznmUhKHFGqau6lb7+FmTPhwQdDVdSWLeEvp1tugQEDYO+9o4w8c8QyaDAvb9vYmKuugtWrQwKZOhVuvx3Gjw/dfQcPhjPOCG1GkhxV3Bkkk1XyR0FeXt7WtcfXr1+/dX9VU6hX9l7AdhMfDhkyhJkzZ9KhQwcmTZrEnJ1MoePuVb53KmnKkSiVn3vphRfgiSegf//Q6+mss8J4ilGj4N134c03w12GkkZ8unUL/4axDhzcdVc4/XR4/HH4/HO4805o1Aguvzy0iRxzDPz1ryG5S8br0aMHjz/+OOvWrWPNmjU8Wdr+RWiXWLBgAQCPPvro1v2rV6+mRYsW5OTkMGXKlK1TqPfo0YNp06axefNmPv/8c2bPnl3p565Zs4YWLVqwceNGHnzwwa37y0/XXqpt27YsX76c+fPnb339pk2bYj7Pyt63OpQ4otSzJ9SpE3oHAdx6K5xySugZdfbZMHduaNO4/vrQVVZSb/fd4bzzwrX44AO4+urQW23o0JDcBw6Ep5+GjRujjlSqqXPnzvTv35+OHTty2mmnccQRR2w9dvHFF3PXXXdx2GGHsXLlyq37hw8fzuTJk+natStLly7deifRt29fDjjgAH76059y3nnnceSRR1b6uddddx2HHnoovXr1om3btlv3DxgwgFtvvZVOnTrx4Ycfbt2fn5/P9OnTGTlyJB06dKBXr17b3QXtzJAhQxg2bBgdO3ZkXenYrWrStOpRWbsWJk2Cm2+Gzz4LbRSnnbatgbtOnagjlMq4h7amKVNC29OqVaF78y9/CSNHhkGUErN0m1Z9zJgxNGjQgIsvvjjqUFJG06qnu2XLQpXH3nvDiBGhAfaBB8LgvClT4Gc/U9JId2ZhIOEdd8Dy5aF6sWfPUKXVpk0YU/Lqq1FHKZI0ahxPleJimDABHn44NHSfeipccEGod0+TBi8pI9aVBvPz4aSTwrZ8OUycCHfdBTNmhN5uF18cjqlbb8Yo7cIqldMdRzJt3hwaWHv0COMHnnwyVGV88EEYO3DYYUoa6ai0t9vo0eFx3rzYXteyZRih/tln8Kc/hYb1U0+Ftm3D3UjpvFqyg2yoMk9X1fm3T5vEYWYjzWyJmb1rZrdUcHxvM5ttZotLypwfRZwxWbMmfHG0aRO+OD77LHTpXLYsPO6zT9QRSlVqutJggwbhD4SlS8Md5m67wW9+E3pkjR4dxuXIVgUFBaxatUrJIwLuzqpVq7brThyLtGgcN7OjgN8DJ7r7D2bW3N2/KlemBdDC3d8ws4bAAuAUd39vZ++fssbxL7+EcePgL38JYwEOOyxUR51ySmpnlZWaSfTa5u5hvqxx42DWrPCegwbBhReqtxxhMN2yZcvi6iEkiVNQUECrVq2oU65ttarG8XRJHA8D97j7C3G85glgors/v7OySU8cq1eHL4UJE8ICSP36hYRx6KHJ+0xJrljbOOK1dGn4PZk0KfyunHJKqN5Ko15FIpAZiWMh8ARwArAeuNjd51dRvhB4GWjv7jsdhZW0xLF+fai7Hjs2dMns3x+uuy5MKihSlRUrQkP6hAmha/bQoTBmDOy1V9SRiQBp0h3XzF4ws0UVbCcTenc1AboClwAPWyVj682sATAD+F1VScPMzjWzYjMrXlF2HYpE2LQpjBxu0wYuugi6dAm9pqZNU9KQ2DRrBtdcEwZ4jhwJkyeH350rrgh3sCJpLF3uOP4O3OTuc0qefwh0dfcV5crVAZ4CnnP38bG+f8LuONzD/FG//z0sXgyHHAI33himMBepiY8+Cg3nU6eGxvQrr4Thw8PMxyIRSIs7jp2YCRwNYGZtgHxgZdkCJXcg9wGL40kaCTNnTqjrPvXUMA5jxowwyEtJQxJh333DhJYLFoQFty68MCyF+8ADFa/9LhKhdEkc9wP7mtkiYBpwlru7mbU0s2dKyhwODAaONrOFJVvvpEf25ptwwglh1bj//Ces27BoUUgg5WvT5s0LdyCx9vuX2iNR175z523Tu++2W5iRt3PnsP5KGtQOiAChH29t37p06eJxW73afcAAd3DfbTf3cePcv/++8vJz57rXq+eemxse586N/zMlMyXr2m/e7P7gg+6FheH38Jhj3IuLE/PeIjsBFHsl36npcseRfho0CAP3fv/7UP980UVVL5ZU00FjkrmSde1zcsLEie+/H9apWLgwLOB19tnw9deJ+QyRalDiqExOTljP+/rrw/oMO1O60lxu7vYrzUntl+xrX7cunH9+6IF1ySWhB9ZPfhIa0lV9JRFIi15VyZaykePJGjQm6S+V137hQjj33LDc7fHHh0kVNY2NJFjaDwBMtrRcj0OkJjZvDoNPr7gi/HzNNfC732lKfkmYTOiOKyLxyM0NAwffew+OOw4uvTTMwPz661FHJllAiUMkk+29dxiU+thjYRqTrl3ht78NMzSLJIkSh0hUEjnup2/fcPcxfHiYA6tdu7AyoUgSKHGIRKG6i0VVZdddQ9KYOxcaNw4z7556ahi4KpJAShwiUUjmuJ+uXeGNN8LdzLPPhruPqVMT9/6S9ZQ4RKKQ7LEfderA5ZfDO+9A+/Zwxhlw1llq+5CEUOIQiUK3bmFlweuuq/kKg1XZf3946SW4+uowYWKnTup5JTWmxCESlW7dYNSo5A8YzMsLi0S99BJs3AiHHx6qsTZvTu7nSq2lxCGSLbp3D6PO+/YNAwd79VLDuVSLEodIJqlpF94mTWD6dLj//lBlddBBYRyISByUOEQyRaK68JrBr34Vel7ts0+4AznvPPj++8TGK7WWEodIpkh0F942bcKYj0sugbvvDlO2v/VWIiKVWk6JQyRTJKMLb34+3HJLWHXwv/+FQw6BP/1J07VLlZQ4RDJFMrvwHnssvP12mDDx/PPh5JNh9erEvb/UKppWvSpaX0MyXby/w+5h2pILL4T99gvzXf34x8mOUtJQVdOq56U6mIxR2hC5YUO4nU/mIC2RZKjO77BZmK69Qwfo1y9UXT30EPTunZqYJSOoqqoyWkNcMl1Nfod79IDi4nDX0adP6AKcBbUTEhsljspoDXHJdDX9HW7dGl55BQYMCAMGBw6EtWuTEalkGFVVVaa0IVJtHJKpqvs7XL5d5MEHwxxXl10G778fBgwWFiYraskAahwXkW2qahf5+9/D3UdeHjz6qO7CazmtOS4isamqXeSEE2D+fGjePHTfvf12tXtkKSUOEdlmZ+0iBxwAr74KJ54Y1jb/9a/hhx+iiFQipMQhItvEMsiwUSN4/HG46qowWWLPnrB8ecpDleiojUNEqu+xx+DMM6FePfjFL2DQIHUkqSXUxiEiyXHqqWGCxFWr4M47w91HdWftlYyhxCEiNfPZZ2HEOYQG9QkToo1Hkk6JQ0RqpmdPqFsXcnLC9sgjYcbdLKgGz1YaACgiNVN2oOHhh4eqq8suCw3m48eHZCK1ihKHiNRct27bGsW7d4c99oDbboPPP4e//S3ckUitocQhIomVkxPuNPbaK6wuuGJF6L67667bymjJgoyWNonDzEYCI4BNwNPufmkl5XKBYuA/7t4nhSGKSKzM4OKLYc89w/rmRx4Jzz4LLVpoyYJaIC0qH83sKOBk4CB3PxAYV0Xx84HFKQlMRGpm0CB4+mn44AM47DBYskRLFtQCaZE4gPOAm9z9BwB3/6qiQmbWCjgRuDeFsYlITRx3XEgOa9eGxvNmzbRkQYZLl8TRBjjCzF4zs5fM7OBKyt0GXApsSVlkIlJzRUUwdy40bhzWNL/22uSsnS4pkbLEYWYvmNmiCraTCW0tTYCuwCXAw2alI4q2vr4P8JW7L4jx8841s2IzK16xYkWiT0dE4rX//vCvf0HbtnD55aH9Q0kjI6Wscdzdj63smJmdBzzmYeKs181sC9AUKPuNfzhwkpn1BgqARmb2gLsPquTz7gHugTBXVYJOQ0RqYo89QrXVaafB0KHw3XdhjXPJKOlSVTUTOBrAzNoA+cDKsgXcfZS7t3L3QmAA8GJlSUNE0ljDhvDUU3DKKWFq9ttuizoiiVO6JI77gX3NbBEwDTjL3d3MWprZMxHHJiKJlp8PDz8cJkm84AL4wx+ijkjikBbjONx9A7DD3YO7Lwd6V7B/DjAn6YGJSPLUqQPTpsEZZ4QxH5s3w6UVDt8KNGgwbaRF4hCRLFWnDkydGkabX3ZZSB6jRu1YToMG08pOE4eZ7RbD+2xx929qHo6IZJ28PHjggTCu44orQvK48srty1Q0aFCJIzKx3HEsL9msijK5QOuERCQi2ScvL0yGmJMDo0eHBHH11duOl66FXnrHoUGDkYolcSx2905VFTCzNxMUj4hkq9xcmDQpPI4ZA1u2hEez7aduVxtH5GJJHLFcIV1FEam53Fy4//7weO214c7juuu2JQ8ljLSw08Th7uvLPjezXYD17r65sjIiItWWkwN/+UtIHjfcEJLH2LHblqeVyMXSOJ5DGHB3BnAw8ANQ18xWAM8A97j7v5MapYhkl5ycsJJgbi7cdFNIHjffrOSRJmKpqpoNvACMAha5+xbY2tvqKOAmM3vc3R9IXpgiknVycuDOO8PjrbeG5DFunJJHGoglcRzr7hvL73T3r4EZwAwzq5PwyEREzGDixHDnMX48uIdR5koekYolcXQys+7AQuCl0rYNM2sMbHL37ypKLCIiCWEGf/xjeJwwIUzNftVVUUeV1WJJHIMJkw4acLaZvQ/8seT5vYT2DxGR5ClNGqtXh/EdTZpoVt0IxZI4LgOGAScQBvrtVfJ8BukzSaKI1HY5OXDvvfDNN2FW3d12C/NcVUXzWyVFLN1xvwfGA+PNrC5wANAMaAl0Tm54IiJl5OWFiRF794azzoJdd4U+fSouq/mtkmandwxlV+Jz9x/cfZG7z3b3BwnLuFJ+tT4RkaQpKIAnnoBOneD00+HllysuV9H8VpIQsVQ1zTazkWa23VxUZpYPHGlmk4GzkhKdiEhFGjaEZ5+FwkL4+c/hjTd2LFM6v1Vurua3SjALq7VWUcCsABhKGAC4D/ANUI+QdP4B3OHuC5MaZQ0VFRV5cXFx1GGISKItWwaHHw7r1sErr0CbNtsfVxtHtZnZAncvqvDYzhJHuTeqQ1gLfF0mTaOuxCFSiy1dCt27Q716IXnsvXfUEdUKVSWOWNo4RpvZRQDuvtHdP8+kpCEitVybNvDcc6G31XHHwcqVUUdU68XSxjEYuKv8TjP7tZlVsFSXiEiKdeoETz4Jn3wCP/sZfPtt1BHVarEkjnUlXXLLm0IF64SLiESiRw945BF48004+WRYr0m7kyWmxGFmLcrvdPcfgE2JD0lEpJr69IHJk0ODeP/+sElfUckQS+L4A/CEmf2o7E4zaw5sSUpUIiLVdcYZcPvtMGsW/PrXYSVBSahYRo4/Ymb1gQVm9iphssMc4HRgTFKjExGpjhEj4Ouvw7xWzZvDLbdEHVGtEstcVbj7ZDN7DOgLHAisBQa6u/q4ikh6Gj0avvgirOWx//5w7rlRR1RrxJQ4ANx9DfC3JMYiIpI4ZvCnP8HHH8Pw4WGU+XHHRR1VrRDLOI5DzOxCMzvazHLL7G9sZg2SG56ISA3k5cH06dCuXZjXatGiqsvPmwc33hgepVJaj0NEardGjeCpp+DQQ0Ovq9degz322LGcZtONWSy9qi4DlgDHA3sCxwCLgWtifL2ISLRatw7JY8UKOOkk+L6CoWmaTTdmO/3id/fv3X28ux8H9AZGEAb+vQZ8luT4REQSo0sXmDoV5s+HM8/csZuuZtONWVyTHO7wYrPc0jXI05kmORSRrcaPh4sugksvhZtv3v6YZtPdqqpJDnfaxmFmb7h7hSv9lSaNqsqIiKSVCy6ADz4IYzv23x/OOWfbsW7dsj5hxCKWxvGfmNnbVRw3YNcExSMiklxlu+med17opturV9RRZZRYEkfbGMqkfXWViMhWpd10u3eHfv1g7lw48MCoo8oYsUw58n+pCEREJKXKdtM98cTKu+nKDmIeOW5mF1awezWwIBFLx5rZSEKPrU3A0+5+aQVlGhPGjrQHHBjq7hqpIyLV07p1WMejR4/QTXf2bKhfP+qo0l484zCKgGHAXiXbuUBP4C9mtsOXfDzM7CjgZOAgdz8QGFdJ0T8Cf3f3tkAHwngSEZHqKyqqupuu7CCexLE70NndL3L3iwiJpBnQAxhSwzjOA24qWeMDd/+qfAEza1TyWfeVlNmgJWxFJCFOOQXGjYMZM+CKK6KOJu3FkzhaAxvKPN8I/Mjd1wE/1DCONsARZvaamb1kZgdXUGZfYAXwVzN708zuNbNdavi5IiLBBRfAsGFhbMeUKVFHk9ZibuMApgKvmtkThC64PwemlXx5v7ezF5vZC4QpS8r7fUkcTYCuwMHAw2a2r28/OjEP6AyMdPfXzOyPwOXA6Eo+71xCdRqtW7eO7QxFJHuZhQWg3n8/TMF+4IHQWcPTKhLXyHEz6wIcTkgcr7j7goQEYfZ3QlXVnJLnHwJd3X1FmTJ7Aq+6e2HJ8yOAy939xJ29v0aOi0jMvvoqtHuYwYIF0LRp1BFFoqYjx9cQejBt3VXmmLt7o5qHyEzgaGCOmbUhzMa7smwBd//CzD4zsx+7+xLCZIs7vdMREYlL8+bw2GNhjEf//vDcc2Hch2wVyySHDd29UZmtYZktEUkD4H5gXzNbBEwDznJ3N7OWZvZMmXIjgQdLRrJ3BMYm6PNFRLYpKoK774YXX4RRo6KOJu2kRRp19w2EGXfL719OmJG39PlCQm8uEZHkGjIEiotDb6suXWCAlh4qpfU0REQqM348HH44DB0Kb70VdTRpQ4lDRKQy+fnw6KPQpAn07Qtff71jmSxcbjYtqqpERNLWnnuGgYE9esAvfwlPPx0We4KsXW5WdxwiIjvTtSvccUfoYTW6zNCxLF1uVolDRCQW55wTthtvDHcgkLXLzaqqSkQkVrffDm+/DWedBW3bhmqpf/4z65abrdGa45lCI8dFJGH+85/QPbdRI3j9dWjcOOqIkqKqkeOqqhIRicdee4WeVh9/DIMGZeU07EocIiLx6t4dbrst9LC65pqoo0k5JQ4RkeoYPjyMLr/2Wpg1K+poUkqJQ0SkOszgrrvCvFaDB8NHH0UdUcoocYiIVFdBATzySEgiAwaEsRxZQIlDRKQmCgvh/vvDmuWXXx51NCmhxCEiUlOnngojRsCECVnR3qHEISKSCLfeCp06hQbzTz+NOpqkUuIQEUmEggKYPh02boSBA8NjLaXEISKSKAccAH/5C8ydC1dfHXU0SaPEISKSSAMGbJsM8bnnoo4mKZQ4REQS7bbboH37ML5j+fKoo0k4JQ4RkUSrXz+0d6xdC2ecEdbrqEWUOEREkqFdu7D405w5cP31UUeTUEocIiLJMmQInHlmmAhx9uyoo0kYJQ4RkWS64w5o0yZUWX31Vdg3b15oPJ83L9rYqkkrAIqIJFODBvDww3DIIaGx/KqroFevMK9Vfn5YQTDDVg5U4hARSbaDDoI//hGGDQsTIm7YEBrMN2wIbSBKHCIisoNzz4UXX4QZMyA3N+zLzw9rlWcYtXGIiKSCWRhVXlgITZrAqFEZWU0FShwiIqnTqFEY3/Hf/8KiRdC1a9QRVYsSh4hIKnXpAjfcADNnwuTJUUdTLUocIiKpdsEFcOSR8NvfwscfRx1N3JQ4RERSLTc33G2YhQGCGTYliRKHiEgUfvQjmDgRXnkFxo2LOpq4KHGIiERl0CDo1w9Gj4aFC6OOJmZKHCIiUTGDu++Gpk1DElm/PuqIYpI2icPMRprZEjN718xuqaTMBSXHF5nZQ2ZWkOo4RUQSavfd4a9/hXffhSuuiDqamKRF4jCzo4CTgYPc/UBghwo/M9sL+C1Q5O7tgVxgQEoDFRFJhuOPh9/8BiZMCIMC01xaJA7gPOAmd/8BwN2/qqRcHlDPzPKA+kDtW1pLRLLTLbfAj38cpmL/73+jjqZK6ZI42gBHmNlrZvaSmR1cvoC7/4dwJ/Ip8Dmw2t3/keI4RUSSo359eOAB+OILGDEi6miqlLLEYWYvlLRNlN9OJtxJNAG6ApcAD5uZlXt9E0J11j5AS2AXMxtUxeeda2bFZla8YsWKpJ2XiEjCFBWFadenToVp06KOplLm7lHHgJn9nVBVNafk+YdAV3dfUabM6cAJ7n52yfMzS8oM39n7FxUVeXFxcVJiFxFJqE2b4Igj4P334Z13oFWrSMIwswXuXlTRsXSpqpoJHA1gZm2AfGBluTKfAl3NrH7J3cgxwOJUBikiknR5eTBlSlirY8gQ2LIl6oh2kC6J435gXzNbBEwDznJ3N7OWZvYMgLu/BjwKvAG8Q4j9nqgCFhFJmv3339bD6vbbo45mB2lRVZVsqqoSkYzjDiedBM8/D2+8Ae3apfTjM6GqSkREyjKDe+8Na3gMGhSqrtKEEoeISLraY4+wauCbb8KYMVFHs5USh4hIOjv5ZDj7bLj5Zpg3L+poACUOEZH0N2FC6JY7dGhaTISoxCEiku4WLYJjjgljO669NupoyIs6ABERqcK8eSFpbNgQVg685ZawhkfnzpGFpDsOEZF0NmdOSBqbN4cuuvXrw69+FWkvKyUOEZF01rMn5OeHu426deHKK+Htt0NjeURUVSUiks66dQsjyOfMCUmkW7ewzOx110HfvtC+fcpD0shxEZFMs3JlGEleWAhz54b5rRJMI8dFRGqTpk1h4kSYPz901U0xJQ4RkUx0+umhqmr0aFiyJKUfrcQhIpKJzOCOO0Ivq1//OqXTrytxiIhkqhYtQlXVK6/AnXem7GOVOEREMtmZZ8IJJ8Dll8PHH6fkI5U4REQymRn8+c+QkwPnnBMGCSaZEoeISKZr3RpuvTWM97jvvqR/nBKHiEhtcM45cNRRcNFFsGxZUj9KiUNEpDbIyQmLPm3aBMOGJbXKSolDRKS22G8/GDsWnn4apk5N2scocYiI1CYjRoT5rH77W/jyy6R8hBKHiEhtkpsL998Pa9eGJJIEShwiIrVN27YwZgysWBESSIIpcYiI1EaXXAIvvgi77JLwt9Z6HCIitVFubtLeWnccIiISFyUOERGJixKHiIjERYlDRETiosQhIiJxUeIQEZG4KHGIiNRG8+bBjTeGxwTTOA4Rkdpm3jw45hjYsAHy88M6Hd26JeztdcchIlLbzJkTksbmzeFxzpyEvn1aJA4zm25mC0u2T8xsYSXlTjCzJWb2gZldnuIwRUQyQ8+e4U4jNzc89uyZ0LdPi6oqd+9f+rOZ/QFYXb6MmeUCdwC9gGXAfDOb5e7vpSxQEZFM0K1bqJ6aMyckjQRWU0GaJI5SZmbAL4CjKzh8CPCBu39UUnYacDKgxCEiUl63bglPGKXSoqqqjCOAL9393xUc2wv4rMzzZSX7REQkhVJ2x2FmLwB7VnDo9+7+RMnPA4GHKnuLCvZVuqiumZ0LnAvQunXrOCIVEZGqpCxxuPuxVR03szzgVKBLJUWWAXuXed4KWF7F590D3ANQVFSUvFXbRUSyTDpVVR0LvO/uyyo5Ph84wMz2MbN8YAAwK2XRiYgIkF6JYwDlqqnMrKWZPQPg7puAEcBzwGLgYXd/N+VRiohkubTpVeXuQyrYtxzoXeb5M8AzKQxLRETKMffaX/1vZiuA/6uiSFNgZYrCSQWdT3rT+aSv2nQuULPz+ZG7N6voQFYkjp0xs2J3L4o6jkTR+aQ3nU/6qk3nAsk7n3Rq4xARkQygxCEiInFR4gjuiTqABNP5pDedT/qqTecCSToftXGIiEhcdMchIiJxUeIQEZG4ZFXiMLP7zewrM1tUyXEzsz+VLBT1tpl1TnWMsYrhXHqa2eoyC2RdleoY42Fme5vZbDNbbGbvmtn5FZTJpOsTy/lkxDUyswIze93M3io5l2sqKJNJ1yaW88mIa1OWmeWa2Ztm9lQFxxJ7fdw9azagB9AZWFTJ8d7As4SZeLsCr0Udcw3OpSfwVNRxxnE+LYDOJT83BJYC7TL4+sRyPhlxjUr+vRuU/FwHeA3omsHXJpbzyYhrUy7mC4GpFcWd6OuTVXcc7v4y8HUVRU4G/ubBq0BjM2uRmujiE8O5ZBR3/9zd3yj5eQ1hPrLy661k0vWJ5XwyQsm/93clT+uUbOV71WTStYnlfDKKmbUCTgTuraRIQq9PViWOGNS2xaK6ldyOP2tmB0YdTKzMrBDoRPhLsKyMvD5VnA9kyDUqqQZZCHwFPO/uGX1tYjgfyJBrU+I24FJgSyXHE3p9lDi2F9diUWnuDcJcMx2A24GZ0YYTGzNrAMwAfufu35Y/XMFL0vr67OR8MuYauftmd+9IWAfnEDNrX65IRl2bGM4nY66NmfUBvnL3BVUVq2Bfta+PEsf24losKp25+7elt+MeZhWuY2ZNIw6rSmZWh/Al+6C7P1ZBkYy6Pjs7n0y8Ru7+DTAHOKHcoYy6NqUqO58MuzaHAyeZ2SfANOBoM3ugXJmEXh8lju3NAs4s6YHQFVjt7p9HHVR1mNmeZmYlPx9CuNaroo2qciWx3gcsdvfxlRTLmOsTy/lkyjUys2Zm1rjk53qULLpWrlgmXZudnk+mXBsAdx/l7q3cvZCwrtGL7j6oXLGEXp+0WY8jFczsIUJviaZmtgy4mtAwhrvfTVjrozfwAfA98KtoIt25GM6lH3CemW0C1gEDvKR7RZo6HBgMvFNS9wxwBdAaMu/6ENv5ZMo1agFMNrNcwhfow+7+lJkNg4y8NrGcT6Zcm0ol8/poyhEREYmLqqpERCQuShwiIhIXJQ4REYmLEoeIiMRFiUNEROKixCEiInFR4hARkbgocYikiJl9F0OZXcxsYsnoXpG0pMQhkl6GAXWB7lEHIlIZJQ6RFDOz0Wb2vpk9b2YPmdnFZQ6fQFj0aWE00YnsnBKHSAqZWRFwGmF9jlOBojLHCoBcwsqOL0USoEgMsmqSQ5E00B14wt3XAZjZk2WOHUBIHO+7+8YoghOJhRKHSGpVtKBOqWZAG8IynyJpS1VVIqn1CvBzMysoWR3wxDLHWhIWfsoxsyaRRCcSAyUOkRRy9/mERXXeAh4DioHVZpZHaNvYE7gb2BxZkCI7ofU4RFLMzBq4+3dmVh94GTjX3d+IOi6RWKmNQyT17jGzdkABMFlJQzKN7jhERCQuauMQEZG4KHGIiEhclDhERCQuShwiIhIXJQ4REYmLEoeIiMRFiUNEROLy/+0lU2405gh8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_1 = plt.plot(logell, logCell + logell, 'r.', label='Data')\n",
    "plt.plot(logell, quadratic(logell, *popt), 'r-', label='quadratic fit')\n",
    "plt.legend()\n",
    "xl_1   = plt.xlabel('$\\lg \\ell$')\n",
    "yl_1   = plt.ylabel('$\\lg [\\ell \\, C(\\ell)]$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit standard power spectrum with higher amplitude (small-scale normalization $\\sigma_8$)\n",
    "\n",
    "The normalisation is increased by 10%, thus $\\sigma_8^+ = 0.9097$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fit by varying all coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "logell_ps, logCell_ps = read_Cell('{}/normalization+dp/P_kappa'.format(basedir), ell_min, ell_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "p0 = [1, 0, 0]\n",
    "popt_ps, pcov_ps = curve_fit(quadratic, logell_ps, logCell_ps + logell_ps, p0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.9254462 ,  0.31677232, -0.18913766])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt_ps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fit by only varying amplitude"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "tilt_glob, a_glob  = popt[1], popt[2]\n",
    "p0  = [1]\n",
    "popt_ps_A, pcov_ps_A = curve_fit(quadratic_c, logell_ps, logCell_ps + logell_ps, p0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Best-fit ampl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.92487861])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt_ps_A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit standard power spectrum with higher matter density $\\Omega_{\\rm m}$ (tilt and amplitude parameter)\n",
    "\n",
    "Increase by 10% gives $\\Omega_{\\rm m}^+ = 0.3366$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fit by varying all parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "logell_pO, logCell_pO = read_Cell('{}/Omega_m+dp/P_kappa'.format(basedir), ell_min, ell_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "p0 = [1, 0, 0]\n",
    "popt_pO, pcov_pO = curve_fit(quadratic, logell_pO, logCell_pO + logell_pO, p0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.87065394,  0.32647783, -0.19338417])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt_pO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fit by varying only shift parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "ampl_glob, a_glob  = popt[0], popt[2]\n",
    "p0  = [0]\n",
    "popt_pO_x0, pcov_pO_x0 = curve_fit(quadratic_u0, logell_pO, logCell_pO + logell_pO, p0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.326074])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "popt_pO_x0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABkeElEQVR4nO2dd3xO1xvAvyebxN57lZoRW8yorUZbFK2iLUrRofSnVZWqokZ167D3Lmp0GLESNYOgdhBUIiSE7Dy/P04SEtkyOd/P537e99577r3PeW9yn3vOs5SIYDAYDAZDSrHIagEMBoPBkLMwisNgMBgMqcIoDoPBYDCkCqM4DAaDwZAqjOIwGAwGQ6qwymoBMoPChQtL+fLls1oMg8FgyDEcOnTopogUSWjfU6E4ypcvz8GDB7NaDIPBYMgxKKUuJbbPTFUZDAaDIVUYxWEwGAyGVGEUh8FgMBhSxVNh4zAYDGknPDwcHx8fQkJCsloUQwZgZ2dH6dKlsba2TvExRnEYDIYk8fHxIU+ePJQvXx6lVFaLY0hHRAR/f398fHyoUKFCio8zU1UGgyFJQkJCKFSokFEaTyBKKQoVKpTq0aRRHIYE8fCAyZP1Z0YeY8gZGKXx5JKWe2umqgyP4OEBrVtDWBjY2MC2beDsnP7HGAyGnIkZcRgewc1NK4DISP3p5pYxx4AZpRiSx8fHh27dulG5cmUqVarEu+++S1hYWIJtr127Ro8ePZI9Z6dOnQgICEiTPK6urkyfPj3Zdg4ODknuDwgI4Mcff0yTDFmNURxPCal5QLu46FGDpaX+dHHJmGNiRinjxulPozwM8RERXnrpJV544QXOnj3LmTNnCAoKYuzYsY+0jYiIoGTJkqxevTrZ827evJn8+fNngMQpxyiOdEApNUIpdVopdUIpNTWRNh2i25xTSo3JbBlzKql9QDs766mmzz9P+ZRTWo5J6yjFkANIp6Hk9u3bsbOz4/XXXwfA0tKSmTNnMnfuXO7fv8/8+fPp2bMnXbp0oV27dnh7e1OzZk0A7t+/z8svv4yjoyO9evWiUaNGsamHypcvz82bN/H29qZatWoMGjSIGjVq0K5dO4KDgwH49ddfadCgAbVr16Z79+7cv38/SVkvXryIs7MzDRo0YNy4cbHbg4KCaN26NXXr1qVWrVqsX78egDFjxnD+/HmcnJwYPXp0ou2yJSKS5QvQCtgK2EavF02gjSVwHqgI2ABHgeopOX+9evXkaWbSJBFLSxHQn5MmZbVEGnd3kVy5tEy5cul1Q/bj5MmTqTsgHW/sN998I++9994j252cnOTo0aMyb948KVWqlPj7+4uIyMWLF6VGjRoiIjJt2jQZPHiwiIgcP35cLC0t5cCBAyIiUq5cOfHz85OLFy+KpaWlHDlyREREevbsKYsWLRIRkZs3b8Zeb+zYsfLtt9+KiMj48eNl2rRpj8jUpUsXWbBggYiIfP/992Jvby8iIuHh4RIYGCgiIn5+flKpUiWJioqKI2tS7TKDhO4xcFASeaZmlxHHUGCKiIQCiIhvAm0aAudE5IKIhAHLgW6ZKGOOJS3TSJlBWkYphhxAOg4lRSRBr5+Ht7dt25aCBQs+0mbPnj307t0bgJo1a+Lo6JjgNSpUqICTkxMA9erVw9vbGwAvLy+aN29OrVq1WLJkCSdOnEhS1r1799KnTx8AXnvttTiyfvzxxzg6OtKmTRuuXr3KjRs3EuxTStplB7KLV1UVoLlS6gsgBBglIgfitSkFXHlo3QdolNgJlVKDgcEAZcuWTV9pcxgxD2g3N600stMD2tk5dfJ4eGTPfhgeIuZNJcbF7jHeVGrUqMGaNWvibLtz5w5XrlyhUqVKHDp0CHt7+wSP1S/NyWNraxv73dLSMnaqasCAAaxbt47atWszf/583FKgABNSckuWLMHPz49Dhw5hbW1N+fLlE4ybSGm77ECmjTiUUluVUl4JLN3QCqwA0BgYDaxUj96BhJyNE/3LEJFfRKS+iNQvUiTBlPJPFc7O8NFHOftha4zpOYR0HEq2bt2a+/fvs3DhQgAiIyP54IMPGDBgALlz507y2GbNmrFy5UoATp48yfHjx1N17bt371KiRAnCw8NZsmRJsu2bNm3K8uXLAeK0DwwMpGjRolhbW7Njxw4uXdLZyvPkycPdu3eTbZcdyTTFISJtRKRmAst69OhhbfTU2n4gCigc7xQ+QJmH1ksD1zJH+uxFdnVh9bjiweTdk/G4kjGCGWN6DiKd3lSUUvz222+sWrWKypUrU6VKFezs7Jg0aVKyx7799tv4+fnh6OjIl19+iaOjI/ny5UvxtT///HMaNWpE27ZtqVq1arLtv/nmG3744QcaNGhAYGBg7PZXX32VgwcPUr9+fZYsWRJ7rkKFCtG0aVNq1qzJ6NGjE22XHVEpHc5lqBBKDQFKisinSqkqwDagrDwknFLKCjgDtAauAgeAV0Qk6YlHoH79+vKkFHLKrEA7jyseuHm74VLeBecy+gIhIeDjA6VLg50dnD0L69dD//6wYd9x3pq+mchTnbGqvpF+dV6hYuEy7NsH330Ho0fDpUtw4YI+PjQUWraEc1dv4vTqSlaPG4CDbW5EoFw5yJsXeveGrl0hMBB++gl8feGHH7TisLaGzZvhuefSv++GuJw6dYpq1apltRipJjIykvDwcOzs7Dh//jytW7fmzJkz2NjYZLVo2Y6E7rFS6pCI1E+ofXaxccwF5iqlvIAwoL+IiFKqJDBbRDqJSIRSajjwJ9rDam5KlMaTRkJv3emlOO7cgZ07waLUQXpuaU3opdrI+q7MnnOMJs84snw5fPYZ9OkDwcFgbw9LlkDz5vBW92pEhtcCIGJnDebuhC5d4OpVfe4TJ7S8ISFQpAgUKAA3Q66zfd9dtleeioVdMUre68SVi7mwtYXz53Xbnj3h4EH4+OO4soaGwuHD0LAhXL+up68qVYKKFR8spUtrhwDD08n9+/dp1aoV4eHhiAizZs0ySiOdyBaKI9pLqm8C268BnR5a3wxszkTRsh2Pa3cMDtZKx8EBrl2D996DoUOhVSs9GujaFRq9ACFXZyDhNlDoNJ5+t5n4jiMXL+pzrFunH8xdu2pl4+AAM+ad4X/7+hNudw2bXOH8/cZ6mpV/oNFOnnxUlsm757Ou4jgiJRJ69uL1BlPpUXoktbT+wc0N6tQBR0etRCZPhrVrYf58uHgRVq2CChW0PAcPwurVum8xbN4MHTvC/v2wYsUDQ3ypUqn7zQw5kzx58piS0RlEtpiqymiepKkqSLlnUVSUbpsvH9SsCbdvQ6FCMH06jByp39SbN4dJk2DfPv1Avn07+mDbQFTV9dj1HMK2ftu4d9oZe3utMIoWhYTyoiU0vZVkP6540Hpha8Iiw7CxtGFbv20pOi6GXbu0AunfX6+3bQt+fvDVV1oJurlB4cLQqBG8/roepQCUKQONG+sptGLFUny5p5acOlVlSDmpnaoyiuMJY8cOENFz/5GR+sHZsyfMmqWni1xdISICzp2Dy5e1orC2hqlTtc3C2Vk/VG/be7DrcsqVQFpJi7JJrH1EBNy6pRUb6NFUrlx6pBIaClWrQsmSWnEcPgx79ui2Y8fC9u0PRiSNG+s2Bo1RHE8+RnEkwJOsOHbs0HaEvtETfQ0bQu7c+m3b318ri+rVYezUy/wyTcezFC784AH5zjt6qiknkNYRioeHVgxXr2rHgu7d9dRXwYLwySfa3rJkCRw6pLcD1KgBx4/rkZVIwiOspwWjOJ58cqpx3JBCDh/WnlSjR+v1uXP1lE2M4vjiCzhwQE9j7dkDGzfC2WAPFoS+jepeE5uyh1n/3myalM15AR1u3m6ERYYRKZGERYbh5u2WrOKI74UWE9AbHg7jx0OLFlqJtmkDVarAhAl6ai8iQttv8uWDBg2gbFl44QV4/nk93WcwPM1kl5QjhkSIiNBvyzFG302b9APv3j29Pm0a/PuvfpuuUQPatdNTL7duaTf6ypX1Azei0HGk1mIi8p9m5yW35C+cDSs5uZR3wcbSBktliY2lDS7lXZI9JrHYjzx54H//e2AjypMHpkzRiuWdd7QSKVpUd6VhQ21g799f20RcXLTh3ZA5JJee/ElgwIABsVl9Bw4cyMlob5KUxKvE5+FEjxmFGXFkQ0JD9fSInZ32IurVS48qmjeH4cNh8GA9RbVhA5QooV1kS5TQc/iDB2tvp4fLB7vY6AduzBRPsg/cbFrJybmMM9v6bUuVTSSlXmjFi2tFEkPlytpGUq+e7ka1atozq1kzrbxjnAi8vfWor1s3qFv36Z7Sym5ERERgZfV4j7jIyEgsM9mne/bs2bHfJ02axMfxfdGzAWbEkcXEf0m/dEm/6S5dqtc7dtTKo149OHpUjyaqVNExEsuWPXiAWVjAmjXw7rtxlQZEP3Drfs3nFq3ZVvfr5B+4mVnJKZU4l3Hmo+Yfpdhgn9bsF88+C19+qZUNgJWVnqKaMkXbPu7c0b//vn16erB+fShfXitxH5+09e1JIqOzG/z+++80atSIOnXq0KZNm9hkgK6urgwePJh27drRr18//Pz8aNu2LXXr1uWtt96iXLly3Lx5E4DFixfTsGFDnJyceOutt4iMHtY7ODjw6aef0qhRIzzidSCxVOsDBgxg6NChtGrViooVK7Jz507eeOMNqlWrxoABA2KPd3Bw4IMPPqBu3bq0bt0aPz+/R/rm4uLCwYMHGTNmDMHBwTg5OfHqq68+MpKYPn06rq6uABw6dIjatWvj7OzMDz/8ENsmMjKS0aNH06BBAxwdHfn5558f/8eH7JFWPaOX7JpW3d1dxM5ORCkRKyu9HhUl8sEHIv/8o9vcvfug/euv6/Z9+4ps2SISEpKKC6UmzXVa0mKn8JjwyHC5fve6RERGiIjIoT/myeLPX5Yla1xlybElsvTYUll2fJmERYSJiIjndU9Z6bVSVp1YJatPrJbVJ1bLxtMbJTwyXEQk09JOx6dRI5GePfV3Pz+RkSNFXFx06nobG5Hbt7NErAwhtWnV0ztdfkx68oe5detW7L3/9ddfZeTIkSKiU57XrVtX7t+/LyIiw4YNk0nRdQS2bNkigPj5+cnJkyelc+fOEham/86GDh0amxIdkBUrViQoS2Kp1vv37y+9evWSqKgoWbduneTJk0eOHTsmkZGRUrdu3di07YAsXrxYREQ+++wzGTZsWOzxq1atEhGRli1bxqZ/f7jv8dOwT5s2TcaPHy8iIrVq1RI3NzcRERk1alRsu59//lk+//xzEREJCQmRevXqyYULFx7pV2rTqpupqizg8mUdsHb6tDbSiugX9Zgo8GnT9Jz6wIGwfLmO5q5XT781z5ihvYBSRWrDzdOQTvd+/dqc+u0HLh/chm/lUvR0epaCwPp/1/P1P19zI+gGvvd88Q/2B8DnfR9KnbjMpmmD+LR5BBxHL9EEjgnE2tKaRccWMcNjRpxrKRRh43Tp0GGbh7HixAqKOxSPXUrnKc2Xbb8E4KTfSQ5fO8y5W+do/0z7dHMt9vDQow7Q01PffKNjY+bOhb17tVEdoF8/PQJ8801tYH8ayMjsBjH4+PjQq1cvrl+/TlhYGBUeGmZ37dqVXLlyATq1+m+//QZAhw4dKBD9z7Nt2zYOHTpEgwYNAAgODqZotB+3paUl3bt3T/C6Xl5efPLJJwQEBBAUFET79u1j93Xp0gWlFLVq1aJYsWLUio5krVGjBt7e3jg5OWFhYUGvXr0A6Nu3Ly+99NJj/xaBgYEEBATQsmVLQKd037JlCwB//fUXx44di7WfBAYGcvbs2Ti/V1owiiOTeNilc8oUWLAAfv897vx748bw/ffwyy96OiR3bp3eI29efVxMxHNqYx/SFG6eSL7zKIniwu0LFMldhHx2+fj7/N8M3zKcs/5nkZhkxafAqclLNCrdiCiJIjIqkupFqtOqfCuK2helqH1R7G3swc2N4fui6OUJWFgg772LDB4MgL21TpU9qskoBjgNiE2RLQjhkeFYWeg/3VblW6FQ/HfvP24E3eAfn384aHEwVnG89ftb7LmyB4AJuybQslxLujzbhZHOI5PsfnJBlko9UA4FC+opq5Il9XLrlo54nzdPu0QvXqyVfseOMGiQ9syytk7+FuRU0jGreqKMGDGCkSNH0rVrV9zc3GKnbIA4adZj/m7iIyL079+fyZMnP7LPzs4uUbtGUqnWY9KzW1hYxEnVbmFhQURERILnSygNe2JYWVkRFRUVux6Tcl0SqVkSs++7776Lo+DSA6M4MgF3d3j1VdiyRRuwP/oIPvxQz4lv3apdZrt00ak1XnhBG2Z/+imu0oghTbEMj1GQw/++P4uPLea473GO+x7Hy9eL++H3WdZ9Gb1r9qaIfRFqFKlBn5p9qFW0FhULVKSofVGKOeiQ7BervciL1V5M+OQuLhT43JYCAdFPmFY9oXDcjKAxo4jE6FmjJz1r9Ex0v1NxJ9yvuBNFFIJw3Pc4KGIVh8t8FyIlktrFauuleG3uXazJ8+1yp9jOr5S2ccQQFKSVfsWK2pC+ZYt2ZNiwAV58UUesDx+e+PlyOplR/yUwMJBS0W9SCxYsSLRdTGr1//3vf/z111/cjjYKtm7dmm7duvH+++9TtGhRbt26xd27dylXrlyS142far1UKvPXREVFsXr1anr37s3SpUtp1qxZku2tra0JDw/H2tqaYsWK4evri7+/Pw4ODmzcuJEOHTqQP39+8uXLx549e2jWrFmclO7t27dn1qxZPPfcc1hbW3PmzBlKlSqVaA2TlGIURwYgor2eChYEJyedfK9yZZ0nCnRUclSUzq30xRc66OyLL/QD6NQp/daaGGmJZQBSVDHpdvBtdnjvYOuFrbQs15JeNXtxL/we7/35HkVyF6FWsVoMqjuIWkVr0aRME0A/mNf2WpvCXyYBmTL4CfNKrVeYc2ROrKL9vc/vNCzVMHZ/neJ12H9tPwuOLiAoLAiAWha9CAtbTmQkhJbZwpYdzjg750/xNVu21C8LMXz/PRw7ptOg/PknNG2qty9aBCtX6iDERomWJMuZpLZAV1Lcv3+f0qVLx66PHDkSV1dXevbsSalSpWjcuDEXYxKpxWP8+PH06dOHFStW0LJlS0qUKEGePHkoXLgwEydOpF27dkRFRWFtbc0PP/yQrOKISbVerlw5atWqFaeeRkqwt7fnxIkT1KtXj3z58rFixYok2w8ePBhHR0fq1q3LkiVLYo32FSpUiJN2fd68ebzxxhvkzp07zuhi4MCBeHt7U7duXUSEIkWKsG7dulTJnCCJGT+epCWzjePBwSKFCon06/fovvBwkYULRapW1YbUKlVEZs8WiYhI2bndL7tLrgm2YumqJNcEW3G//HiWx6ioKPlk2yfS8NeGYvGZheCKOExykC92fRG7//rd6491jXTF3V0XTU+FxdX9srtM2jUpyd8qMipSzvmfk7Un18qPm3ZJrlwiFgUuCa6IhauFOM92lvE7xsvey3tjDfMpxddXZM8e/T0qSqRNG5E5c0TmzRMpWFD/HbRpI7JzZ6pOm2mkuuZ4NiIkJETCw/X9cnd3l9q1a2epPAkZ+rMDqTWOZ/lDPTOWzFAc338v0qmTfjCIiBw6pBVIfFau1L96rVoiy5enXGHE4u4u7pVsZFILJe6VbFL1AI2MihTP654yfe90+WTbJ7HbG/7aUJrOaSrjd4yX3Zd2x3o0ZTvS210nmUt9/kWY/Lhpl4zbPk4a/dooVrEu8NTeNzfv3ZTzt86n6ry3b4s8/7xItAOP3LqlveiKFdN/F2+8kc4dSQdysuI4c+aMODk5iaOjo9SvX1/279+fpfI8KYrDTFU9Bjdv6ukoCwtd98HGRkd0OzjoYDDQ01OzZ+tgvkGD9Bz3pk3QoYM+LtW4ueHsHYnzeQHLyBS5rGy/uJ05R+bw9/m/8buv/cbrlajHhFYTUErh/oY7lhY5oHBFZrjrRKOnWqyB5kBzJrSawK3gW2y/uJ2W5bT3yjKvZYzYMoJKBSrRrlI7OjzTgfaV2mNrZZvoefPn1zatGNav155yu3aBp+cDz6u7d3X3Onc2QYWPQ+XKlTly5EhWixFLUFBQVouQLpgAwDTi6akr1a2Nnt5/6y347bcHCQPv3tUZZ8uX1yks/vxTb7eygk6d0qg04IHLSoymSsBlJUqi2H1pNyER2uti7+W9/H3+b9o/05753eZz5f0rHBx8MNYTI0coDUhR3zOSgrkK0qN6D4rY6xr2Xap04buO31G9SHUWHVtEt+XdKD2zNHdDUz7v3bUrfPutjkgfMUJ3zdNTp7jv2lV7Z61erW1iBkN2wWTHTQVXrugYjKZN9UvvmDF6FFGlStx2S5boh8Dt29C+vY72bt78sS//gER8RU/4nmDxscUs9VrK5cDLrHl5DS9Ve4l7YfewtbKNdWF9LMLCtI/pzZuPLv7+Ol+KUnGWUBXFNct7lBAH7JQ15ywD2Wh9kX7h1Slo6cCB/PeYa3OC8ZXepHipZ9mvrrH4+p+Mb/05hXIXYp/PPladWMWnLT8l35GTHNm+lF2VLHnzpYk42DgQEBJAWGQYRXIXSZV7Y1Kk1uU5LDKMbRe2cfj6Yca2GAvA0I1DKe5QnL6OfalUsFKy5xDRmYxLloQ//tDZAyZNgjNndMqTjz/W3nmZPQIx2XGffFKbHTfL7Q+ZsaSXjaNZM5Fnn31gx3iY4GCRwED9fft2kW7dRDJrOvVG0A2pPau24IpYfmYpHRZ3kMVHF8vd0LvJH/ww9++LeHqKLFsm8umn2rrfqZNIw4YiFSuK5M2rJ+LjLTfskVn1kcsl7UWKFJH91fNJreGWcqBqXpECBWSjk73givxTxV4kTx5ZX9tOcEUOVrATsbGRNdWQIqORswX1+ZbWRPL/D/EumVukYkWZ272iOIyzlGuDXxH5+GOZPvVFwRUJuHxWRESm7pkquCJ3Qu6IiMicw3Ok1fxWsbaaPZf2yOKji2Oj1ZPD/bK75JqYSyw/s5RcE3OlyAEhvs0+IjJC2i1qJ8pVCa5IkzlNZNaBWeJ/3z/J89y+LXLunP4eEKCj02fMEKlZU+S551IkfrqTk20chpRhbBzpyIULuprc1KnaL/+777RN4+E3PhE9lfDhhzqw6/vvdRnWVq0yTq7AkEDWnFpDUFgQ7zR6hyK5i1C5UGXerPMmL9d4OTaGIlECArTfb/zl4kXdIdCdLF0aihYlqnAhLlUrgUOBohQpVJZL+eGV+4v5pNIAOlZ5nv+UL0PXtGF59zn0qtmLPDf/peLWMViP/QyK16bO3WvMO/8X5T7oCA7FaBcRwq3wYPJ9mg+UBS/du8dLvr7Q3Rd8fekTvfCyL/j58bqvL69vyAO+O8DXl/ejInndDvK6VobChWnXsBz21VrgMHsh1KiBsgjCQllgbamj7FadXMV8z/m8UusVAKbuncrRG0dZ8pL2d/e950s+23yxtonUujwnnN/Rkj/7/smVwCssPb6URccWMXTTUG4H3+aj5h8RFhmGiDxiD8mfXy+gg0C3b9dxP++9p6sbgs6F9eabejRSr14K/2gMhvQkMY3yJC1pHXHs2qWddxJzkzx4UKR5c/3i7egosm1bmi6TYs75n5Phm4ZL7i9yC65Io18bJZ+r6coVkSVLRIYNE2nVSqR48bgjBltb7eL18ssi48eLLF8uIYcPyPSdk+WPs3+IiPYewhWZvne6iIjcCbkjLea1kM1nNouISFhEmFwKuCSRUZEZ2X1NZKTI1asif/0lMnOmyMCBIs7Oj46GihTRyaOGDZOQH76RC26/aV9oEfl85+fSfUX32FN2WdpFav1YK3Z94s6JYjPBJsUjjkmTtKMX6M/o1EhxiIqKksPXDst/d/8TEZFlx5dJgSkF5J3N78iFW4/mDorhYc+8Dz4Q6dhRd71oUZ3j7PXXRa5dS/5nexyyesSRnTyR5s2bJ1evXo1df/PNN+XEiRMiIlKuXDnx8/OT27dvyw8//JBVIqYJ446bjopDRMQ/kZmFWbP0P26RIiK//JIGt9pU8u2+b0W5KrGeYC39f+sv+67se1RpREWJnD2rgwT69xepUOHBg9TBQWfmGzBA5MsvRTZs0G2jBf9yz5cy9/BcEdFuu4W+LCQfbf0o9tQLPBfI6ZunM7aTj0NUlFaSf/wh8tVXIm++KdK4sUiePA9+g9y59dTb1Kl6HjFakWw8vVGWHV8We6paP9aSZnObxcZ+TNk9JVaJJkRavIT3++yX3qt7i9UEK7H4zEJeXvWyHLh6IMljvvlGZNQo/T0gQOT990WsrfWtnTIl4SnU9CCnKo6Y+I3HISLeP/bDCQjjE6M44icjzAkYxZHOiuNh7t8XuR4dC3fhgsjo0fofOIY0xKYlGpwWHhkuy48vF68bXiIicvS/o/Lx1o/l2p2HXi8jI0WOHxf54QeRXr1ESpR48JAsVEjkxRf1W/mhQ49otq89vo4Ty9FkThN5be1rsesBwQHyRLB3rx5VKaWf7OXKPfiN8uRJUJHcD7svvkG+IqLvQ7FpxWT0X6NFRI8c3t3yruy9vDfOZdJy70VErgRekQ//+lDyTs4rVb+vmuJsv+fO6ZeWOXO0PS0j4z/SojhSEnSZUhJSHBs2bJCGDRuKk5OTtG7dWv77T4/kxo8fL4MGDZK2bdtKnz59xNfXV9q0aSN16tSRwYMHS9myZcXPz09ERBYtWiQNGjSQ2rVry+DBg2OVhL29vYwbN04aNmwou3fvjr3mqlWrxN7eXqpUqSK1a9eW+/fvx1EkMYqjV69eYmdnJ7Vr15ZRMZo+m2MURwYojqgokRUr9DOnY8eE26QpE3kCUeABwQEyw32GlJ1ZVnBF3v/j/bgHXbqkXz27dXsQdgwipUqJvPKKHgqdOPHI6+esA7Pk5VUvx64PXD9QOi5+0JlsG/T3uCQ0j3T9ur6hQ4eKVKsWV5F07KhHZP/8oxWzaEN3jKPBlcArkn9Kfpl9aLaIiPjf95dBGwbJ8RvHH0vMwJBAOfbfMRERuRt6V5rMaSK/HvpVgsMTiCIV/eLSvbuetRPRjhlRUSIHDugZyegs3ulCqtOqp8G5ICmyU1r1+COOhBTH0zDiMMbxZDh0SBsm9+zRSQhHjUq4XVpi09x2LSQsIpRICwgLD+WT30Zw4N4Z7obdpUW5Fnzb4Vu6PNsFzp7VASNr1uiC4qATYL3wgi6a3aKFDhh5yGp/7MYxfj30KzPaz8DG0oaIqAh87/nG5mz6pcsvcVxXYwzJTxwJpWotXhxeflkvAP/9pyPw3Nz0ElMKsEQJ6N4dy549cYhOMFU6b2lujr5JRJTOdnr65mlWnlhJ/9r9ATjjf4ZNZzbxep3XyW+XP0nR4npV56VWMZ2G++qdqwSHBzPo90GM3T6WEQ1HMLT+UArlflDsvEIF7ZQRwxtv6M/XXtN5serW1Wn5J07UhcEykzTnU0sFWZVW3aAxAYBJsHIlNGig62b88gscPgzPPZdw27TEprl4g1UkWESCTRSUDIikc5XOHBi4n531f6DbCk8sajvpQJExY/R78eTJ2rH/3DmYM0cXwq5QgZDIUDac3sCNIF0JzTvAm/lH53PK7xQAwxsOZ0f/HdhY6pJ26RXvkO1JSQnAGEXy449w8qRWJEuW6LazZ+ushaVLw7Bh4OaGpRDrDeVcxhn/D/0RESbvnsxPB39i1N+jCIvU9ULO+p/lSuCVRy4Z44k1bpz+fLjQ3LOFn+XQ4ENs77edeiXqMW7HOMp+XZbLgZcT7KKIDiBs1kyXsD13Tos6b55Orvndd4//M6aGtNSGTy0jRoxg+PDhHD9+nJ9//jk2xTikLq26p6cnnp6enD59OjY1e1Jp1Q3RJDYUeZKWtE5VBQSIjBsX146RFKmZ575w64K8+nM7UeOR17sh7hWtJWr2bJExY3TmQ9Dz8s2aaTvFpUuPnONOyB25EXRDRERO+p4UXJEf9mtvjrCIsESnOQyp4O5dnVSse3c9Bwk6sdTQoTpgJyLikamZdafWxR7ee3VvKTK1SGwMScw9SYknVgzHbxyPTTopoj2yLgdcTrT9zp0ihQuLrFql82JNnPh4P0F2tHE4OTnJwYMHRURkwIAB0rJlSxHRU1XTpk2Lbff222/LlClTRETkzz//jJ2qOnHihDzzzDNy44b+//H39xdvb+9ErxdD586dZfv27bHrCU1V3bx5U8qWLfsYPc58jI0jHRVHRuAb5CvvbH5HrCdYS66JuWTMV53kdstGD9xkLS11qtRZsx5Y4hMgNCI01p0zhp3eO59cW0V24O5dbRvp0eOBEilaVA51ayht+luI5TgdgDlp1wMtcNb/rGw6syl2vf4v9WXg+oFpztcYGBIoeSfnFdvPbeXDvz6U28G3H2nj5SXSp4/IvXt6/Y6Oi5QNG0QmTBAJS+WfSFZ7VSmlpFSpUrHLjBkzZN26dVKhQgVp1qyZjBo1KlHFcePGDXnuueekTp068t5770mJEiUkJLrm8vLly6V27dpSq1YtqVu3rnh4eIhI0opj9erVyRrHRUT69OkjNWrUMMbxnLxkF8URGRUpVb6rIhafWcigaS3Ep1ltiS1S/fzzInPnijxU0zg+X+z6Qvr/1j92/ddDv8q+K/syXnDDowQF6VTHL78sEbnsRECuOyBTWlrJQY+1CR4SGRUpU3ZPkcVHdc3pXXvCpOpnneX7jTtSdelLAZek32/9RLkqKTClgEzfOz3R0WVkpEj9+toD8J139J9bnTo6QUBKyWrF8Thkt7Tq2RWjOLKZ4giNCJU5h+fokcDx4/Ln+13lVHkH/dPXqCHy3Xc6z0Qix248vTHWe2T8jvHSe3XvFKfOMGQS9+7Jvz99IWeaVJUoCwsRCwuRrl11TElk4kGR3re9pdI3leT307+LiB5NXAp4dEryYR6eDvW87ikdFncQy88sE42xCQvTs5+rVun1tWv1TJuVlY73DA1Nvns5WXFkt7Tq2RWjOLKJ4oiMipRlx5dJpa8rCq7IqpeiKzfZ2Ii8+qrI7t3JRmz9euhXnePJ559MktoQS1oDM7y9RT7+WAdZgEilSiLTpyc6koyMioyNuJ/pMVMsPrNINJI8semth5XGh399KH+e+zNR8RYu1PmvevfW4q1Zk3yXcrLiMKQMoziygeL4+/zfUvebGoIrUmu4pWx+Bol6ppLItGki0XOgCeEb5CsdFneQVSf06+GdkDuy+cxmM8LIbNKjYFRIiMjSpdq5ISa1S//+Oj4kkReGSwGXYp0bRERmuM+Qb/Z9E7uenEH9dvBtqfB1BcEVabOwjRy+dviRa0ybpjOxRETobsWI4uGhRU6IkydPpjgw0ZDziIqKSrXiyDbuuEqpEUqp00qpE0qpqQnsL6OU2qGUOhXd5t2skDNJIiKIWr2K0bNexP/iCRaut+DIfy/Q8aetqNNndBBI4cJxDrl+9zruV3SB6oK5CnI39C7B4bo4eR7bPHSs3DHn1Mt4UkgoKCe12NpCnz6we7cOrHjjDR2H06iR9vGeO5d/zuxg8u7JeFzRvrhl85Xl7QZvx55i56Wd7L2yN3a9csMLSbp857fLz6lhp/i6/dccuX6Eur/U5dW1r/Jf0H+xbUaNgq1b9Tlq1tS1Yv79F9q00QkTE6o+YGdnh7+/djk2PFmICP7+/tjZ2aXquGxRj0Mp1QoYCzwvIqFKqaIi4huvTQmghIgcVkrlAQ4BL4jIyeTOn171OBIjJPguX/08gLd/OUL+Uxe5WL0kJfsMxvbNwTqILAlazGuB7z1fTg079fTEVmR3Ek53+/jnvXMHFi/W8SInTnDbDubUU3zX3Jblb29PMEguNCIUWytbrt+9TtmvyzLkmSmU9P4gfimWRwgMCeTLvV/y64EFDIk6TqdWBR9p/8cfuiLltm26dsxbb+kQltGjYfx4XbUSIDw8HB8fnzixEoYnBzs7O0qXLo21ddwg4GxfjwNYCbRJ5THrgbYpaZthU1UhIbL1m/ek8vtWgisyr2tZPWmcRMZD79ve8u6Wd+VemPaVPHztsJz1P5sx8hnSTlptHCkhKkoWfTtIVtRQEqGQe9bIvpebJel+fTf0rnz/z/dyzl8X6zh987QsP748yWlMd3cRO4cQsbQUscsdIZ1+6f9Ijq3oMAYR0faPPn30VFjNmtpxzPD0Qna3cQCewGfAP8BOoEEy7csDl4G8SbQZDBwEDqZ7MM79++I78wvp90ouwRV5ZrSt/L10YpLG7pg54l3eu8T2c1vZfmF7om0NTz4xQYPVRljIYidL7Y1lZyfy3nspypM+6s9RkvuL3OJ3L3Gb2cM2EYtCFySfaxlRrkqGbhz6SBLLgACdmf6NN0Q2b9Z1vAxPN9lCcQBbAa8Elm7Rn98CCmgIXCR6Gi2B8zigp6leSum1023Ecfeuti4WKyavvIRYf6pk7OxX5H7ovUQPCY8Ml16reskEtwmx25L6Zzc8PcSJrj5zRhvPLS21AnnnHREfn0SPjYyKlKP/HY1dH7F5hKw+sTru+ePZ+LfuuivvbXlPLD6zkJIzSsrak3HjTU6dEvH11ceNHSuyZYu25Q8a9CCY0PD0kC0UR1IL8Afg8tD6eaBIAu2sgT+Bkak5/2MrjoAAkS++kH8r5RfvfIi0aSOX/lwZm/I8IR5Of97vt34yeffkx5PB8HRw7px+7be01J5Yw4bpOiNJcDf0rjjOcoxNSxIVFRU7hZXQjNt+n/3iOMtRKn9bWUIj4gZyxCibmIw3w4frz5o1tWIxPD3kBMUxBJgQ/b0KcCX+iCN6NLIQ+Dq150+z4vD3F/n0UwkplE8+a4nYfGohPWclXvjZffPPMunzdvLO3J6Sa2Ku2GpvhqeA9LaJXLigqxtaWenYn6FDRS5dSjQHVGRUZKwS2HZhmzzz7TNyyi/xJ31YRFhsvEhQaJDMPjRbIqMi40xvKaW79McfOveVvb3IokXp0z1D9icnKA4bYHH0lNVh4Lno7SWBzdHfmwECHIu2iXgCnVJy/jQpjoAAkfz5xa0c8uwYB8EV6b26t1y/m7ABc+Xaz8VuLGL5KWI7Fhky+wW5E3In9dc15DzSI+4jMS5eFBk8WMTaWiKtreTXBpZS5gOLJOtc7PLeJZ2WdIpNQ3Ip4FKSRvQf9/8ouCLOs51l0Z/HH+nK33/r2A9nZ/3E+O239OueIfuSlOLIFnEcIhImIn1FpKaI1BWR7dHbr4lIp+jve0REiYijiDhFL5szTKh8+Zj5VjdcXoc7+Quz+ZXNLOu+jOIOxR9peif0Dq8dcyXUCiItIMICyl6/Tx7bPBkmniEbkR5xH4lRvjz8/DOcO4fn8/V47XAkp76NYtT2EPac/jvBQ5qXa86mVzZhZ2WHiNBteTe6r0y8vsSQ+kNY+MJCzvif4Y1/6tLr53F8OiEk1gv5xg24exe2bNHlBbp00ceFh6dfNw05jMQ0ypO0pGXE4e4uYpcvQFQrV7HLG/TIS2R4ZLhsPb81dn3S4rf0iGMckmss4r7551Rf05BDycgRx8OXuewu1UbayppqiIAEly6uk1Al4823wmuFbDy9UUR0NcPEUrL7BvnKa2tfE1yJU0ZY5IGHeViYyE8/abNLpUoPpq4y0nvZkDWQ3aeqMnpJi+JILr3DtL3TBFfiGMhjbBxGaTyFZNKTM8bG4bX8O5FatfQfaMuWidaKjW8T+fXQr2I30S7JUrd/nP1Dzt86LyIi98LuxUk3smqVvuTSpSLNm+vvXbpoR7AM1puGTMYojjSOOOK/RJ7yOyUnfE+IiM5kuvbkWpPDx5B1hIfrui2FCumMvIMHa3/aaBKq/X054LK47nCN/bu9cOtCon/DUVFR0m1ZN3l+yfPiG/TgvNFlKyQ8XIedxJRsT0lRKkPOISnFkS1sHNmR+BVH6zcMp/XC1oz+ezQAeW3z8mK1F02aEEPWYWUFQ4bomvQjRsDcubpW7MyZEBaWYO3vMvnKMN5lPEop7obexXmOM8M3D0/0Eq0rtObvC39T+6fabL2wFYDGjfW+q1dh/nyIKc+tVMrLJhtyNkZxJEH9huFU7LKCxo0Fa0trlnVfxrxu87JaLMOTgoeHriH/cMHxtFCgAHz9tU6m2LgxjBwJjo68cNE2ydrf9jb2TGo9idfrvA5AUFhQnISISilGNBrB/oH7yWeXj3aL2vG/v/8XW0+9WDGdu3HaNJ2/cdw4/ZJVu/bjdceQA0hsKPIkLWmN41jouVBwxaQHMaQ/GWVQj4oS2bhRpHJlEZBbrZzlp0Xvpaj298dbP5a8k/PKxtMbH4kVCQoNkkEbBkmxacVi69zHZ+xYka1bdcXBd9/VU1mGnAvGxpE2xREeGS5/nfvL2DEM6U9y3hePS2ioyIwZOgGVtbXIZ58lW2z8zM0z8u6Wd2PtInaf2z2icGJsHRGREXFqqQcEiFSpokvUxtg9OnTQ2xPCeGFlf5JSHGaqKgmsLKxoW6mtsWMY0h8XF5IsrvG42NjoKauzZ6FHD50nvUEDOHIk0UMqF6pMMftisXaRkMgQZrjPiNOmiH0RAOZ7zuf5pc/T77d+3A29S758sH+/nnmbORO+/FLX/WjSBC5ciHudmKz148bpz8edqTNkPkZxGAxZQXzvi/So95EQRYvC0qWwbp2O5GvQAD75BEJDE2zuUt4FG0sbLLDAUlnSr3Y/AKIkKk67/k79cW3pypLjS6jzcx0OXD1AvnxaD4aEaDt98+Zw/Tr07699rmLIyHhJQ+aQLQo5ZTQZXcjJYMgR3L4N778PCxZAjRowb55WJPHwuOKBm7cbLuVdYotLvbLmFUo4lGBG+7gjkD2X9/Dq2le5dvcaP3f+mTfqvAFoxfHMM7qOmVL6e1QUWFhkXJ0sQ/qSVCEnozgMhpyEh4d+RU+uBGBSbN6sy/1du6Zrybq6Qq5ciTaPjIpk1F+jKGJfhI+bf/zI/tvBt+m5qidVClXhNcfXHqlkOG+e9sBauhRKltTTWfv3P343DBmLURxGcRieBNLzVT0wUNeI/fVXePZZPURo0iRFh7pfcWfh0YVMazuNPLZ58LjiQeuFrQmLDMPG0oZhDYYx0nkkJfKUIDJSi1i0KJQpAz/9BF27wpIl4OCQNtENmUNSisPYOAyGnEJ6Ggfy5dMZC//+WxslmjXT01j37yd76P6r+/n7woMEiw8HGoZGhPLt/m+p/2t99l/dj6Ul7NihFcWPP8LUqbBxo77c5ctpF9+QtRjFYTDkFDLCE6tNGzh+HN5+WwcROjrCzp1JHvJe4/c4PvQ4eWzzECVRXA68jLWFNZbKElsrW2Z3mY2NpQ0t5rVg4dGF2NtrPQVa11WvDhcvQtu2EBHx+F0wZD5mqspgyEmkh40jMXbuhDffhPPn4Z13tE+tnV2Sh7hfcafFvBZ80uITbC1tYw3qN+/fpOeqnrh5u/FRs4+Y1HoSoCPMAwN1gLu/v/a8MmRPjI0jrYojI/9JDYbMILV/w/fuwccfw7ffQp06sGKFzn+VBKdvnqZKoSoopTjld4rKhSpjZWFFeGQ4I/8cSaPSjejr2DeuWFc8WPD7WSrZNcTmblWeeQaefz59u2J4PJJSHFke1Z0ZS5oix93dxb2SjUxqocS9ko0JcTXkPB4nrcnvv4sULCji4CCyZEmKDgkIDpDCUwvLwPUDE9y/5uQaOel7MjZrL5X+FIqckKqOd8XKSmTx4ozpiiFtkETkuFXm6rCcg8f2hbTuHUaYJdhEhrFt+0KczWuOISeRkDE9pX/DnTuDpye88gq8+ips365HIblzJ3pIPrt8fNfxO5yKOwH6pTQm60JoRCgj/xzJreBbvFD1BZ0o8eWXsAgtTK92I9j15Qf07atDTYYnkKz3cbpiSH+McTwR3MpDmKUuBRtmodcNhhzF4xrTy5TRLlFjx2p33YYN4cSJJA/pXbM3VQtXBWDUX6N4/4/3ERFsrWzZ/fpuqhSqwuJji1FKYWEbjG1BX9pXb4KjI1SooLPDT5iQ/l0xpC9GcSSCS4t+2FjZYonCxtoWlxb9slokgyF1pDWtycPp3q2sYOJE+PNP8PPTkeZz58bNIZIAIkJEVAQRURGxo44y+cqw+/Xd9KnVh4ioCGoUrcGfff+kcWlnLC2hY0cYMAAKF06/rhgyBmMcT4KEUi8YDE80SQUZ/vcf9O2rt73yio7my5MnydNFSRQWyoJzt85x5PoRetboiYgw3X06t0Nux3pbwQNddO2aPu3Fi9p119o6ozprSIqkjOPGxpEEzmWcjcIwPF0kZUwoXlyPPKZMgU8/1XlDVq7U3leJYKH0pMbUvVNZfXI1bSq2oUCuAoxuOjq2zdH/jlIodyFK5y1NZKT2rrK3h6NHoVUrfYkkMqIYsgAzVWUwGB6QnDHB0lLbPNzcIDhYB2R8/32yU1c/dPqBnQN2UiBXAQCu370OQERUBD1W9aDJnCac8juFpaWeJfviC5g+HTZtgg4ddOyHIftgpqoMBkNcUhowcfOmNkps2gQtW0KLFtpQkYwBYunxpQz+fTC7X99NnRJ1OHL9CB2XdCQ8KpyNfTbGGeW/957WS46O8PmcAxwL2mqmjjMJEwBoFIfBkDFERemn+3ff6XVbW+2JlYTyuHrnKjP3zWRKmylYWejZ8gu3L9BuUTuu3b3Gqp6reL7K84SGQs2a2lh+xDOSyEbTkefGYmNpw7Z+24zyyGBMkkODwZAxWFjoohsW0Y+S0FCdPDEJSuUtxfR207GysOJO6B3e3vQ2+e3ys/eNvVQrUo2fD/2sXXht9cDnjz/grZ/mEtniUyIlkrDIMNy83TK8a4bEMYrDYDA8Hi4ueqRhYaGrNi1cCF99lazdA+Afn39YcHQBx28cp5hDMdz6u7Gs+zKUUtwLu0fJkkK+fNCjRU3U6lWw8xNk4Z88a9UuRaI97FlsSD/MVJXBYHh8YuwiDRvq/Olr12r7x08/aaWSBH73/GJrmf8X9B/FHYpzL+weLea3oHnZ5nzV/ivCQi1o1+0WqvBpjmxsQOFCVri5QdmySYtkKg2mHTNVZTAYMhZnZ/joI/2kXrUKxo+H+fPhued0rfMkiFEaB64eoOI3FVl9cjW5rHPRomwLvvnnG/qu7YuFdRg7/yjIziXObN9mhb+/PrWPT+LnNbXNMw6jOAwGQ/piYaHL0a5cCUeOQP36+vNhEphDql6kOm/Ve4vnKjyHhbLgq/Zf8WWbL1nmtYznlz5PUNhdQBvLLSy00njuOR0wmBAmTUnGkW0Uh1JqhFLqtFLqhFJqahLtLJVSR5RSGzNTPoPBkEp69oS9e7Xdo2lTPRKBB3NI48bpz2jlYW9jz8wOMymYqyBREsW4HePoV7sf87vNZ8fFHQzeOBiA0qWhf39YsEDHJNrYJHx5k6Yk48gWkeNKqVZAN8BRREKVUkWTaP4ucArImynCGQyGtFOnDhw4AC+9BC+/rCPObWySTXV7wvcEM/fNpFy+cgyqN4hiDsVikydaWelihaBPuXevNq3cufNonitnZ6MwMoLsMuIYCkwRkVAAEfFNqJFSqjTwPDA7E2UzGAyPQ7FiOi3766/r1Ldbt+oEVEnMIdUqVotTw04xsO5AABqXbkz5/OWJjIrk420f43NHGzf++ENXEXzuOX0a3wSfHIb0JrsojipAc6XUP0qpnUqpBom0+xr4EIhK7oRKqcFKqYNKqYN+fn7pKKrBYEg1trYwZ4520921S883jRyZ5BxS2XxlUUpx9c5Vqv1Qje/++Y7zt8/z/f7vcZnvwpXAK7RvD7/+Cp99Bhcu6BLqN29mct+eQlKtOJRS9kopyzQct1Up5ZXA0g09ZVYAaAyMBlaqmFzMD47vDPiKyKGUXE9EfhGR+iJSv0iRIqkV12AwpDdKwfvvw+bNOkX7/Pl61JEMhXMXpke1HrSu2Joqharw12t/4XffD5cFLvjcvczAgdpUsnw5/PuvVh63buljPa54MHn3ZDyumECO9CRZxaGUslBKvaKU2qSU8gX+Ba5HG7GnKaWSLkgcjYi0EZGaCSzrAR9gbXTFwv3oEUX8rPxNga5KKW9gOfCcUmpxKvpqMBiyA+3bwz//gIODnmP6888km9ta2fJdp++oXqQ6AOdunWPNy2u4ef8mLvNduBx4GYCzZ3UGlFOntO3D44oHrRe2ZtyOcbRe2Nooj3QkJSOOHUAl4COguIiUEZGiQHNgHzBFKdU3qROkgHXAcwBKqSqADRBnwCkiH4lIaREpD/QGtovI417XYDBkBc8+C+7u8Mwzukzt4pS9A168fZE31r/B7ku7+fu1vwmJCOHi7YuAnvk6eBDWr4cvvwQ3bzfCIsNMmpIMICVeVW1EJDz+RhG5BawB1iilHrfUylxgrlLKCwgD+ouIKKVKArNFpNNjnt9gMGQ3iheHnTvhhRfgtdd0oOAHHyTe3sODCm5u7G0wC6cW/bC2tOb8O+fJZa2LdQRH3MfJSddE37YNvFb1R/lexKLqamzsQ3Ap75LxfXpaEJEUL0AzYBhQ8aFtFVJzjqxY6tWrJwaDIZsSEiLSo4cIiIwaJRIZ+Wgbd3eRXLlELC31p7u73Am5Ix0Xd5T9Pvtl7uG5UuHrCnLx9kUREfnoI5EqVUQsrSKlfO1LsuOMR+b26QkAOCiJPFNTaxwvAjQExiulvlJK1QcmJXOMwWAwJI6trbZsDxumqzcNGADh8SY5Esgf4h/sz7lb57hx7waOxRwJCAnAZb4LF29f5Isv9LTVooUWeB8ty7f/a0xkZFZ07skkVYpDRH4D3gB+ALzQdo4zGSCXwWB4mrC01DU9Pv8cFi2Cbt3g3r0H+xPIH1I+f3m83vaic5XO1CtZj2Xdl3En9A4uC1y4GHCBPHmgVy9o0gR++03rpacgp2umkBKvqnFKqdiJRxGJFJH9IjJXRGaKyPiMFdFgMDwVKAWffKIDM/78U3tcxQRlJJI/xMZS5xs5cv0IL618iVFNRhEUFoTLfBcCQwJRCqpV081/+UUHsSeFScOeMpJNq66UOgM4icj9eNsHAkVEZHIGypcumLTqBkMOY/166N0bypXTSqRcuSSbh0SEMHbbWMY0G8PVu1fZ6b2Tdxu/C2gXXaVg/36dASWx3FYmDXtcHjetenB8pRHNIsC4wxoMhvSnWzf4+2/tadWkCRw/nmRzOys7ZrSfQRH7IjgWcyS3dW7CI8M5cPUAFwLOoRSULKlL0Y4bp8uFxMekYU85KVIcSqkS8TeKzisVkf4iGQwGA9CsGezerb83b65TlaSA7Re3M3jjYFZ4raDPmj60WdiGq3euUrAgVK0KW7ZAnz6PKgaThj3lpERxzADWK6XijBWjM9gmmzPKYDAY0kzNmjpQsHhxaNcO1q1L9pA2Fduw5/U9vOr4Kst7LMc/2J/2i9sTanGLDRvgr7+gcmXo2hU8PR8cZ9Kwp5xkFYeIrEJ7UR1SSm1USk1USk0C9gLTM1pAg8HwlFOunM6d7uQEPXroAlHJ0LRsU5RSFMldhFpFa3H21lmeX/o898LuUbAgvP22dtpq00YnR4whppChURpJkyJ3XBFZAFQAVgLWQAjQR0SWZKBsBoPBoClUSNs8nJ31PNOyZSk67GLARbwDvPniuS/Yf3U/MzxmADrLyYsvavfcefMetDdJEVNGsilHlFIN0RHjnsASEYmM3p5fKeUgIkEZK6LBYDAAefJoA0XnztC3r7Zi903aP8flui3nQ4eQSzWl3mtbaViqIQBly8Lq1eDtrWfBrl8H7widFDEsMgwbSxu29duGcxkz9EiIlIw4XgOeBeoAC5VSnyil8gAKU1DJYDBkJg4OsGkTtGwJ/frp+rGJEe1fm+vTCdC6NVcO76Dh7IZ4+Xoxbe80RITy5XUZ2qZNYfhge0JDMUkRU0BKFMf/gNNAe6A40AadWv2zFB5vMBgM6Ye9PWzcqIMuXn8d5s5NuF08/9oK//5H1cJV2XB6Ax9u/RBXN1cA3nlHZ3o//Icj6reFWIgVNpY2JiliEiQ7VRUdw/EV8JVSyhaojM5ZVRKom7HiGQwGQwLkzg0bNujMum++qaP8Bg6M2ybGvzY6oq/5c6/T3NkZEeGM/xkm7JpAodyFeKfpOzRtChUrwocf9qBhpYrM/DbUTFMlQUpSjsRW4hORUBHxEpEd0YbxD+O3MRgMhkwhVy4dYd6xIwwaBD//HHd/Iv61Sinuh9/H3tqed/94lyXHtI/PoEFgZwf719flnzVGaSRFSlKOuKHrbqwXkcsPbbdBG837AztEZH7Gifl4mJQjBsMTTGgodO+ubR/ff6+zGSbDPz7/cOzGMZZ5LcP3ni9H3jqCtaU1y5bpHIuBgbpUiFVKKhY9oSSVciQlisMOnRH3VbRLbgCQCz1a+Qv4QUQ801HedMcoDoPhCSc0VNeL3bABvv0WRoxI0WF3Qu9w8OpBKhSoQIUCFQC4H51g6dgxaNw4owTO/iSlOFJi4wgBfgR+jK70VxidvyogXaU0GAyGtGJrC6tW6cSI77wDERHw/vvJHpbbOjdvbXqL4g7FqVa4GsMaDKN28dosXQqvvgodOmjHraJFM6EPOYhUDcREl5C9nkGyGAwGQ9qxsYEVK3SA4MiR2ptq1KgkD7GysGJVz1WEhIfQc3VPNpzewN439tKjRyVOnICZM7X9fft2bf8waJKdqoptqNTIBDYHAofMVJXBYMg2hIfrwMCVK2HKFPjf/1J02Cm/U9T/tT4FcxXk4KCDFHMoxurV0LMndOmi7fBPkxvQ46ZVj6E+MAQoFb0MBlyAX5VSHz6ukAaDwZAuWFvDkiV65DFmDEyblqLDCuUuhI2FDTeCbtBteTeCw4Pp0QMqVIDff4eJEzNY7hxEahRHIaCuiHwgIh+gFUkRoAUwIANkMxgMhrRhZQULF2qbx4cf6vJ/yVDUvihHhx5leY/lnPE/w783/wX0wKVdO+3te+fOo8c9jVUDU2PjKAuEPbQeDpQTkWClVGj6imUwGAyPSYzyuHMHhgyBvHm1IkmCsvnKUjZfWWoXq82nbp/ydfuvqV+/COvXw+3bcOaMriJoaanbP61VA1OjOJYC+5RS69F5qroAy5VS9sDJjBDOYDAYHgtra+1t1bEjvPaaVh6dOiXe3sMD3Ny44VSIzWc2k9cmL41LN6a/U38CA6FRI2jVCubPh9KlE64a+DQojhRPVYnI58Ag4Hb0MlhEPhOReyLyakYJaDAYDI9F7tzaSFG7tg4UTKySYMzwYdw4mnR/j/POyzl76yyDfh+Em7cbVavCpEnwzz+6CNS9e09v1cCUpFW/CzzseqUe2icikjcjBDMYDIZ0I29e+OMPXYK2c2fYsQPq1YvbJt7woaD7YVa9twrHnxxpt6gdR946wv/+V4OaNbXi6NIFtm7V01NublppPA2jDUhZBcA8IpL3oSXPQ4tRGgaDIWdQuLAuBlWwoI7sO3Uq7v4Ehg8FchXgrXpvIQhdl3XF/74/zz+vm+7YoWMNKe0BzSbrz6eEFMdx5GRMHIfBYIjl3Dlo1kwbz/fsgfLlH+yLtnHEHz7s8t5Fx6UdmddtHj2r98TPT9G7N+x1j0KNqEpEngtPXPGn9IrjMBgMhpzPM8/AX39pI0XbtvDffw/2JVJ0vEX5FpwZfobfz/zO2O1jKVoU/vwThvw0l/AwCyLv5Xuqij8ZxWEwGJ4+HB1h82a4dk0Hady+newhJfOUxN7anqt3rjLdfTrW1vBi01pEzd8Gi/7COqjCU1P8ySgOg8HwdOLsDOvWwenT8PzzEBSUZHOlFLOen4VSitF/j2bx0cW4VG7Ep5/dw/p2LSptO0y9Yk/GNFVyZBvFoZQaoZQ6rZQ6oZSamkib/Eqp1Uqpf5VSp5RST8ddMhgMGUPbtrBsmfaxfeklnZ49CZRS/Nz5Z+qXrM9r615j4q6JfPZeFZYstOHE4Ty88UYmyZ3FZAvFoZRqBXQDHEWkBjA9kabfAH+ISFWgNnAqkXYGg8GQMl56CebM0R5Xr7yiU7Inga2VLWteXoO9tT0z983EO8Cbnj31jNeSJbrg4JNOtlAcwFBgioiEAoiIb/wGSqm86LxYc6LbhJmaIAaDIV0YMAC+/hrWrtU1ZKOikmxeNl9ZDg4+SJREscBzATfv32TpUp0QcfJkuP6EF5/ILoqjCtBcKfWPUmqnUqpBAm0qAn7APKXUEaXU7Oh0JwmilBqslDqolDro5+eXUXIbDIYnhXffBVdXnU9k3Lhkm1ctXJUjbx2hfP7yVP6uMv6cwcNDF36yt39QSfBJJNMUh1Jqq1LKK4GlGzqCvQDQGBgNrFTqkcz3VkBdYJaI1AHuAWMSu56I/CIi9UWkfpEiRTKmUwaD4cni00/1iGPSJPj112Sbl89fnlYVWtGhUgc2ntlIsWI6qrxpU53VPZmBS44l0xSHiLQRkZoJLOsBH2CtaPYDUegStQ/jA/iIyD/R66vRisRgMBjSB6Xghx+gfXsYOlQHayRD2XxlKeZQjA/++oDlXsuJUPfo0gW2bIEvvsgEmbOA7DJVtQ54DkApVQWwAW4+3EBE/gOuKKWejd7UGpOV12AwpDcxGXVr1oQePeDo0WQP+bLNlzQs1ZC+a/vSZVkXJk4UevWC8eO1x++TRnZRHHOBikopL2A50F9ERClVUim1+aF2I4AlSqljgBMwKfNFNRgMTzx58sDGjZAvn47x8PFJsnmMp1Vu69yc9j9NUPhdfv4ZypbVTlvbtmWS3JlEtlAc0R5SfaOnruqKyPbo7ddEpNND7Tyj7RaOIvKCiCQf7mkwGAxpoXRpHV1+545WHgmV/3u4ed7S/N7nd24E3WDy7smEqFssWKBTYo0cCcHBmSR3JpAtFIfBYDBkSxwd9bTViRPw8ssQHv5om4dqx7Ys35I/+v7B85Wfp8I3FQgsvoG1a3XtKCurJ8dYnpoKgAaDwfD00b49/PST9rZ6+21dvzzG6TOB2rFtnNsQHB7MC8++gIWyoHNnHRzYvTtUrw5TpmRtd9IDM+IwGAyG5Bg4ED7+GGbPjvvkT6h2LJDLOhdB4UH0+60f52+d527ELezstIfv/v1xBik5ElOPw2AwGFJCVBT07atzWy1dqgM1EhhxxKRkv3D7AvV+qYe1hTX57PLxd49DtGqSl+BgCAjQmU3iHZKtMPU4DAaD4XGxsIB583T52QEDdO1yZ2f95P/880c0QMUCFVnWfRk3798kl1UuyhZzYM0auHlT51KMN0jJURjFYTAYDCnF1lYHZlSoAC+8oFOyJ1L8CaDDMx34vNXnHPc9zi+Hf6FCtdu89toDC3l0hdoch1EcBoPBkBoKFtRuulZW0LEj+D6SkzUOHzX/iHEtxtG8bHPq/VKPMv1cadLWlyLPnuWrVR7ZcpoqOYziMBgMhtRSsSL8/rsuO9u1a5IZDS2UBRNaTaB6ker0rN6TQrkKcbhpFW6+XJMR03bx9ZKclwDDKA6DwWBIC40a6QIc+/fDa68lG6ShlOL87fPM8JhBGHeRSAsivF5gxhf5clx8h1EcBoPBkFZefBGmT9d1PFJQwemVWq9w5c4VvWITgmXDX/A5VYqvv85YMdMbozgMBoPhcXj/fejfX9fy+O23JJu+VO0lxjQdQxRRlMtXjj++68QLL8CYMdq+nlMwisNgMBgeB6V0ZHmDBnrKyssryeYTn5tI6wqt8b3nS4k8xfnpl3Ds7ISpU+H8+UyS+TExisNgMBgeFzs7PdrIkwe6dYNbtxJtamlhyeKXFlOrWC1uh9ymz+b2dHOdC2j9kxMwisNgMBjSg1KltK3Dxwd69dKh4YlQ3KE4+97cR9MyTalepDrtWtuwdq0uQLhlC4SEZKLcacAoDoPBYEgvnJ1h1izYuhU+/DDJpkopwqPCsbG0wcbShm7dwNtb6NQpRXb2LMUoDoPBYEhP3ngDRoyAmTNh4cIkmyoU+3z2Mej3QSw6uoj+7vWpVjOMTZt0WpLsilEcBoPBkN7MmAGtWsHgwTrOIxGsLa1Z1n0ZlhaWTN4zGSsLS0aNDeToUe2sdflyJsqcCoziMBgMhvTG2hpWroQSJXSsx/XriTYtl78c87rN49TNUzQt05Q3ehfh3Xf1jFeLFjoZYnbDKA6DwWDICAoX1gkRAwJ0Fack5p5eqPoCwxsM5+dDP+N925vQVu9RtLw/t24lqXOyDKM4DAaDIaOoXRvmz9d1O4YNgyTqH01rN41Dgw9ROl9pzgZ68aLrAurW1anX797NPJFTglEcBoPBkJH07Aljx8KcOfDjj3pbAiUA7azsqFakGlYWVoxoOILvX3sHNzedub1cOfD0zBLpE8TUHDcYDIaMZsIEOHYM3n1Xr48enWDVQIC9l/fywooXGNdiHH0d+/LLtRnkyfsjly9b4uSUNeLHx4w4DAaDIaOxsIDFi6FKFR3fkUQJwKZlmzLAaQATd01k+fHl7L72F1hEMGwY+PtnjfjxMYrDYDAYMoO8eWH9erC01LYOC4tESwB+3/F7ni38LLMOzeKfYVtZu8qW//6DOnV0JvesxigOg8FgyCwqV9ZuukpBzZo6wjyBEoD2Nvas6LGC28G3GbZ5GHXrCh3e3M+VK/D331kgdzyM4jAYDIbMpEOHBzaPY8cSbeZYzJGfOv/EoLqD2Lk3FI+iA8hX/jzu7lkf26EkCfewJ4X69evLwYMH42wLDw/Hx8eHkOyeTczwCHZ2dpQuXRpra+usFsVgSBtRUdC5szaMu7tDvXqJNvXwgNatIcTOG+vgMqxeaUmhQnD4MAwfnnEiKqUOiUj9hPY9tV5VPj4+5MmTh/Lly6OUympxDClERPD398fHx4cKFSpktTgGQ9qwsIBFi7TRomdPOHQIChRIsKmbG4RUXYA0+4KIuf9w4Lhi55W/uLjpZd58E3LlylzR4SmeqgoJCaFQoUJGaeQwlFIUKlTIjBQNOZ9ChbS948oVeP31RIMDXVzA0iYcCpyHHr25WGoKe7yucNM/imvXMlfkGJ5axQEYpZFDMffN8MTQuDFMm6a9rb76KuE2pT2wfP4dsBCiKvxFxbK2rP+qDTbWFvTrp0uAZLbF4alWHAaDwZDlvPuuzmX1v//B3r2P7HbzdiNCwgCtHaa5T6Nu1SJ89502j3TvEcnOnZkrcrZRHEqpEUqp00qpE0qpqYm0eT96v5dSaplSyi6z5UwvQkJCaNiwIbVr16ZGjRqMHz8+q0VKNfPnz+faY46VQ0NDadOmDU5OTqxYsYKBAwdy8uRJACZNmpQeYhoM2RuldDqSChV05UA/vzi7Xcq7YGNpg6WyxNbSFoXi7/N/07jTOVT1tVhYSmLmkYxDRLJ8AVoBWwHb6PWiCbQpBVwEckWvrwQGpOT89erVk/icPHnykW2ZSVRUlNy9e1dERMLCwqRhw4bi4eGR7tcJDw9P93PG0LJlSzlw4ECqjokvj4eHh7Ro0SLBtvb29omeJ6vvn8GQ7hw5ImJrK9K2rUhERJxd7pfdZdKuSeJ+2V18g3xjt3+9bYk0bBQu//wj8t9/ImFh6ScOcFASeaZmF6+qocAUEQkFEBHfRNpZAbmUUuFAbiB9TEPvvZf+GcScnODrrxPdrZTCwcEB0K7B4eHhCc7de3t706FDBxo1asSRI0eoUqUKCxcuJHfu3EyYMIHff/+d4OBgmjRpws8//4xSChcXF5o0acLevXvp2rUrVapUYeLEiYSFhVGoUCGWLFlCsWLFcHV15eLFi1y/fp0zZ87w1VdfsW/fPrZs2UKpUqX4/fffsba25tChQ4wcOZKgoCAKFy7M/Pnz2bt3LwcPHuTVV18lV65ceHh4cPLkyUfalShR4hF5PvjgAwB8fX3p27cvfn5+ODk5sWbNGt58802mT5/O6tWrCQ4OxsnJiRo1arAkO4TLGgwZiZMTfPedLv70xRe6AHk0zmWccS4TN1DQzduNTrUb8I6HFdf/i6R6/VsMf6MIn32W8aJml6mqKkBzpdQ/SqmdSqkG8RuIyFVgOnAZuA4EishfiZ1QKTVYKXVQKXXQL97QL7sQGRmJk5MTRYsWpW3btjRq1CjBdqdPn2bw4MEcO3aMvHnz8mN0hs3hw4dz4MABvLy8CA4OZuPGjbHHBAQEsHPnTj744AOaNWvGvn37OHLkCL1792bq1AczgefPn2fTpk2sX7+evn370qpVK44fP06uXLnYtGkT4eHhjBgxgtWrV3Po0CHeeOMNxo4dS48ePahfvz5LlizB09MTKyurBNslJE8MRYsWZfbs2TRv3hxPT08qVaoUu2/KlCnkypULT09PozQMTw8DB8Jrr4Grq44qT4SgsCB6rOzBa7+9RqRE8O6OAYSW+YOKlcMzRcxMG3EopbYCxRPYNTZajgJAY6ABsFIpVTF6uBRzfAGgG1ABCABWKaX6isjihK4nIr8Av4AOAExSuCRGBhmJpaUlnp6eBAQE8OKLL+Ll5UXNmjUfaVemTBmaNm0KQN++ffn2228ZNWoUO3bsYOrUqdy/f59bt25Ro0YNunTpAkCvXr1ij/fx8aFXr15cv36dsLCwOPEPHTt2xNramlq1ahEZGUmHDh0AqFWrFt7e3pw+fRovLy/atm0LaGVXokSJR2RMrt3D8hgMhkRQSpf+O3wYXnkFjhyBUqUeaeZg48APnX6g95reTNkzhZ5VX2WHd1MmfmZFjxd1bIdFBg4LMk1xiEibxPYppYYCa6MVxX6lVBRQGHh4qNAGuCgiftHHrAWaAAkqjpxE/vz5cXFx4Y8//iBfvnyxD/8hQ4bQoUOHR6awlFKEhITw9ttvc/DgQcqUKYOrq2uc2AZ7e/vY7yNGjGDkyJF07doVNzc3XF1dY/fZ2toCYGFhgbW1dey1LCwsiIiIQESoUaMGHg/VDUiI5No9LI/BYEgCe3tYtQoaNIDevWH7dl2KNh69avZi/en1fLbzM/a9uY81y/LQqhU0bR5J1SqWLFum9VBGkF2mqtYBzwEopaoANsDNeG0uA42VUrmVfrq1Bk5lppDpiZ+fHwEBAQAEBwezdetWqlatSpkyZfD09MTT05MhQ4YAcPny5dgH8rJly2jWrFmskihcuDBBQUGsXr060WsFBgZSKvqtZcGCBamS89lnn8XPzy/2+uHh4Zw4cQKAPHnycDe6NFlS7dKKtbU14eGZM/Q2GLIV1arBL7/Anj3wySeJNvuh0w8Usy9G39/60rBJML0GXuPoEUtuhF7O0HxW2UVxzAUqKqW8gOVAfxERpVRJpdRmABH5B1gNHAaOo2X/JasEflyuX79Oq1atcHR0pEGDBrRt25bOnTsn2LZatWosWLAAR0dHbt26xdChQ8mfPz+DBg2iVq1avPDCCzRo8IhZKBZXV1d69uxJ8+bNKVy4cKrktLGxYfXq1fzvf/+jdu3aODk54e7uDsCAAQMYMmQITk5OREZGJtourQwePBhHR0deffXVxzqPwZAjeeUVGDIEpk6FDRsSbFIgVwHmvzCfofWHYmtlyxdfgH1Jby54Fs+w0QY8xUkOT506RbVq1bJIopTj7e1N586d8fLyympRshU55f4ZDI9FSAg0bQoXLmi7RzL52aIkCq/jFlhb67IfZ8/C88+n7dJJJTnMLiMOg8FgMMTHzk7bO0R0MsTQ0ESbbjqziTo/16FslQCeqRLOjOMfcMIiYzwSjeLI5pQvX96MNgyGp5mKFWH+fJ1B9+OPE21W1L4oJ3xPMHzzcCyUBV5/1cXqv8YZIlJ2CQA0GAwGQ2IUKwaNGulEiO3bQ7t2jzRpUKoB41qMw3WnK12f7crenzPONmhGHAaDwZCdiankdOCA9q995ZVH8lnF8HHzj2lYqiFDNg7h2t2My7luFIfBYDBkZ9zcICxMVw1UCgID4Y03Esylbm1pzaIXFxESEcJyr+UZJpJRHAaDwZCdcXEBGxvtJmVrC8OGwcaNOsI8AaoUqoLX216MdB6ZYSIZxZHF/Pbbbyil+Pfff9PtnG5ubrExIRs2bGDKlCkArFu3LjZleWpwcXEhvjtzUvz77784OTlRp04dzp8/T5MmTQDtWrx06dJUX99geKpxdta1yT//XH/OnAkdOsAHH0AiQbYVC1TMUJGM4shiYiLBly/PmGFl165dGTNmDJB2xZFa1q1bR7du3Thy5AiVKlWKDQQ0isNgSCPOzvDRR/pTKe1llTcv9OmjYz0yGaM4onFx0fcCIDxcry+OzoJ1/75eX7FCrwcG6vW1a/X6zZt6/fff9fp//6XsmkFBQezdu5c5c+bEURxubm60bNmSl19+mSpVqjBmzBiWLFlCw4YNqVWrFufPnwceRG43b96cKlWqxMmOG8P8+fMZPnw47u7ubNiwgdGjR+Pk5MT58+fjjCRu3rxJ+fLlAZ0CpXfv3jg6OtKrVy+Cg4Njz/fXX3/h7OxM3bp16dmzJ0FBQXGut3nzZr7++mtmz55Nq1atAGLTx48ZM4bdu3fj5OTEzJkzU/YjGQyGRylWDObNg+PHIfrFMDMxiiMLWbduHR06dKBKlSoULFiQw4cPx+47evQo33zzDcePH2fRokWcOXOG/fv3M3DgQL777rvYdt7e3uzcuZNNmzYxZMiQOIkOH6ZJkyZ07dqVadOmPZLCPD6zZs0id+7cHDt2jLFjx3Lo0CFAK5eJEyeydetWDh8+TP369fkqXp3kTp06MWTIEN5//3127NgRZ9+UKVNiU6i///77qf69DAbDQ3TqBCNGwDffwB9/ZOqlTRxHNG5uD75bW8ddz5077nq+fHHXCxeOu148oeTxCbBs2TLee+89AHr37s2yZcuoW7cuAA0aNIhNS16pUiXaRftt16pVK84D+eWXX8bCwoLKlStTsWLFdLGV7Nq1i3feeQcAR0dHHB0dAdi3bx8nT56MTfEeFhaGs7NzoucxGAwZzNSpsGMH9O8Px47pkUgmYBRHFuHv78/27dvx8vJCKUVkZCRKqdgiSzHpzkGnOH84/XlERETsvoRSrqcUKysroqKiAB4ZqSR0HhGhbdu2LFu2LMXXMBgMGYidHSxbBvXraxfdjRtjc6l7eOgXWhcXbRpJT8xUVRaxevVq+vXrx6VLl/D29ubKlStUqFCBPXv2pOo8q1atIioqivPnz3PhwgWeffbZRNs+nAYddDqTmGmoh9Oyt2jRIrbqnpeXF8eOHQOgcePG7N27l3PnzgFw//59zpw5k2JZ41/fYDCkAzVrwvTpsHkzfP898CBmcNw4/ZlMOZ1UYxRHFrFs2TJefPHFONu6d++eaq+jZ599lpYtW9KxY0d++ukn7OzsEm3bu3dvpk2bFusmO2rUKGbNmkWTJk24efNB+ZOhQ4cSFBSEo6MjU6dOpWHDhgAUKVKE+fPn06dPHxwdHWncuHGqpsYcHR2xsrKidu3axjhuMKQnw4bpNLijR8Px47i5QWgRDyKdJxNaxCPOVHp6YNKq52AGDBhA586d6dGjR1aLkuk8CffPYEhXfH3B0REKF2bOxO8ZeLATWIZBpA0/N93G4I6pm68yadUNBoPhSadoUViwAE6coPaij7CwDgOLSCysw/B3cEvXSxnjeA5mfkzgicFgMIDOnPv++9SfOZNu9jZseAZsLG1wKe+SrpcxIw6DwWB4kpg8GWrXZsXm3Hzl+CHb+m3DuUz6ulUZxWEwGAxPEra2sHQp1vdCeOenwziXapTulzCKw2AwGJ40qlfXyRDt7TMkl5VRHAaDwfAk8tZbsHq1Tn2RzhjFkQX4+/vj5OSEk5MTxYsXp1SpUjg5OeHg4MDbb78NPEhOCElntfXz86NRo0bUqVOH3bt306lTJwICAggICODHH3/MtD4ZDIZshlKxUeTpjfGqygIKFSqEp6cnAK6urjg4ODBq1KhE269bt47OnTtTvXr1R/Zt27aNqlWrsmDBAgCaN28O6OSHP/74Y6wiMhgMhvTCjDiicZnvwnzP+QCER4bjMt+Fxcd0XvX74fdxme/CCi+dVz0wJBCX+S6sPaXzqt+8fxOX+S78flrnVf8vKIV51ePxcAGmGBJKhx6Dp6cnH374IZs3b8bJyYng4GDKly/PzZs3GTNmDOfPn8fJyYnRo0enSR6DwWBICDPiyObEpENPKELcycmJCRMmcPDgQb6PzlETw5QpU/Dy8ood2RgMBkN6YRRHNG4D3GK/W1tax1nPbZ07zno+u3xx1gvnLhxnvbhDCvOqGwwGQw7ETFUZDAaDIVUYxZEDSEs6cpPC3GB4yvHw0FHk6Z1THaM4cgTx06GnhEKFCtG0aVNq1qxpjOMGw9NGBhfkMGnVDTkSc/8MhiSYPFkrjchIsLSEzz+Hjz5K1SmyfVp1pdQKpZRn9OKtlPJMpF0HpdRppdQ5pdSYTBbTYDAYcgYuLmBjo5WGjY1eT0eyhVeViPSK+a6UmgEExm+jlLIEfgDaAj7AAaXUBhFJOKTaYDAYnlacnWHbtgwrOp4tFEcMSikFvAw8l8DuhsA5EbkQ3XY50A1Is+IQEVQGheQbMo6nYXrVYHhsnJ3TXWHEkC2mqh6iOXBDRM4msK8UcOWhdZ/obQmilBqslDqolDro5+f3yH47Ozv8/f3NQyiHISL4+/snWVvdYDBkLJk24lBKbQUSiowbKyLro7/3AZYldooEtiX61BeRX4BfQBvH4+8vXbo0Pj4+JKRUDNkbOzs7SpcundViGAxPLZmmOESkTVL7lVJWwEtAvUSa+ABlHlovDVxLqzzW1tZUqFAhrYcbDAbDU0t2mqpqA/wrIj6J7D8AVFZKVVBK2QC9gQ2ZJp3BYDAYgOylOHoTb5pKKVVSKbUZQEQigOHAn8ApYKWInMh0KQ0Gg+EpJ9t4VYnIgAS2XQM6PbS+GdiciWIZDAaDIR5PReS4UsoPuJREk8LAzUwSJzMw/cnemP5kX56kvsDj9aeciBRJaMdToTiSQyl1MLHQ+pyI6U/2xvQn+/Ik9QUyrj/ZycZhMBgMhhyAURwGg8FgSBVGcWh+yWoB0hnTn+yN6U/25UnqC2RQf4yNw2AwGAypwow4DAaDwZAqjOIwGAwGQ6p4qhSHUmquUspXKeWVyH6llPo2ulDUMaVU3cyWMaWkoC8uSqnAhwpkfZrZMqYGpVQZpdQOpdQppdQJpdS7CbTJSfcnJf3JEfdIKWWnlNqvlDoa3ZfPEmiTk+5NSvqTI+7NwyilLJVSR5RSGxPYl773R0SemgVoAdQFvBLZ3wnYgs7E2xj4J6tlfoy+uAAbs1rOVPSnBFA3+nse4AxQPQffn5T0J0fco+jf2yH6uzXwD9A4B9+blPQnR9ybeDKPBJYmJHd635+nasQhIruAW0k06QYsFM0+IL9SqkTmSJc6UtCXHIWIXBeRw9Hf76LzkcWvt5KT7k9K+pMjiP69g6JXraOX+F41OenepKQ/OQqlVGngeWB2Ik3S9f48VYojBaSqWFQOwDl6OL5FKVUjq4VJKUqp8kAd9Jvgw+TI+5NEfyCH3KPoaRBPwBf4W0Ry9L1JQX8gh9ybaL4GPgSiEtmfrvfHKI64pKpYVDbnMDrXTG3gO2Bd1oqTMpRSDsAa4D0RuRN/dwKHZOv7k0x/csw9EpFIEXFC18FpqJSqGa9Jjro3KehPjrk3SqnOgK+IHEqqWQLb0nx/jOKIS7oWi8pKROROzHBcdFZha6VU4SwWK0mUUtboh+wSEVmbQJMcdX+S609OvEciEgC4AR3i7cpR9yaGxPqTw+5NU6CrUsobWA48p5RaHK9Nut4fozjisgHoF+2B0BgIFJHrWS1UWlBKFVdKqejvDdH32j9rpUqcaFnnAKdE5KtEmuWY+5OS/uSUe6SUKqKUyh/9PRfRRdfiNctJ9ybZ/uSUewMgIh+JSGkRKY+ua7RdRPrGa5au9yfb1OPIDJRSy9DeEoWVUj7AeLRhDBH5CV3roxNwDrgPvJ41kiZPCvrSAxiqlIoAgoHeEu1ekU1pCrwGHI+eewb4GCgLOe/+kLL+5JR7VAJYoJSyRD9AV4rIRqXUEMiR9yYl/ckp9yZRMvL+mJQjBoPBYEgVZqrKYDAYDKnCKA6DwWAwpAqjOAwGg8GQKoziMBgMBkOqMIrDYDAYDKnCKA6DwWAwpAqjOAwGg8GQKoziMBgyCaVUUAra2Culvo+O7jUYsiVGcRgM2YshgC3QLKsFMRgSwygOgyGTUUqNU0r9q5T6Wym1TCk16qHdHdBFnzyzRjqDIXmM4jAYMhGlVH2gO7o+x0tA/Yf22QGW6MqOO7NEQIMhBTxVSQ4NhmxAM2C9iAQDKKV+f2hfZbTi+FdEwrNCOIMhJRjFYTBkLgkV1ImhCFAFXebTYMi2mKkqgyFz2QN0UUrZRVcHfP6hfSXRhZ8slFIFskQ6gyEFGMVhMGQiInIAXVTnKLAWOAgEKqWs0LaN4sBPQGSWCWkwJIOpx2EwZDJKKQcRCVJK5QZ2AYNF5HBWy2UwpBRj4zAYMp9flFLVATtggVEahpyGGXEYDAaDIVUYG4fBYDAYUoVRHAaDwWBIFUZxGAwGgyFVGMVhMBgMhlRhFIfBYDAYUoVRHAaDwWBIFUZxGAwGgyFV/B9EXA8OUf8O9wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "d1, = plt.plot(logell, logCell + logell, 'r.')\n",
    "d2, = plt.plot(logell_ps, logCell_ps + logell, 'b.')\n",
    "d3, = plt.plot(logell_pO, logCell_pO + logell, 'g.')\n",
    "\n",
    "f1, = plt.plot(logell, quadratic(logell, *popt), 'r-')\n",
    "f2, = plt.plot(logell_ps, quadratic_c(logell_ps, *popt_ps_A), 'b:')\n",
    "f3, = plt.plot(logell_pO, quadratic_u0(logell_pO, *popt_pO_x0), 'g:')\n",
    "\n",
    "plt.plot(logell_ps, quadratic(logell_ps, *popt_ps), 'b--')\n",
    "plt.plot(logell_pO, quadratic(logell_pO, *popt_pO), 'g--')\n",
    "\n",
    "plot_data  = [d1, d2, d3]\n",
    "label_data = ['Original data', 'Larger amplitude', 'Larger tilt']\n",
    "leg_data   = plt.legend(plot_data, label_data, loc=1)\n",
    "\n",
    "plot_fit  = [f1, f2, f3]\n",
    "label_fit = ['3-parameter fit', 'Amplitude fit', 'Tilt fit']\n",
    "leg_fit   = plt.legend(plot_fit, label_fit, loc=3)\n",
    "\n",
    "plt.gca().add_artist(leg_data)\n",
    "\n",
    "text_x = plt.xlabel('$\\lg \\ell$')\n",
    "text_y = plt.ylabel('$\\lg(\\ell C_\\ell)$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABZeElEQVR4nO2dd3hURffHP5NCQg+E3mtoSQgdBKSDUgQbSFFQAUXwtSD+rID1VQFREQso7QWlKgKiYqFYIk2C1ITeISEQWnpyfn9MsinsJrvJJtmE+TzPfbbcmXvP3Zvs2Zkz53uUiGAwGAwGg724FbQBBoPBYChcGMdhMBgMBocwjsNgMBgMDmEch8FgMBgcwjgOg8FgMDiER0EbkB9UqFBB6tSpU9BmGAwGQ6Fh586dF0WkorV9t4TjqFOnDjt27ChoMwwGg6HQoJQ6YWufmaoyGAwGg0MYx2EwGAwGhzCOw2AwGAwOcUvEOAwGQ85JSEjg9OnTxMbGFrQphjzA29ubGjVq4OnpaXcf4zgMBkOWnD59mtKlS1OnTh2UUgVtjsGJiAiRkZGcPn2aunXr2t3PTFUZDIYsiY2NxdfX1ziNIohSCl9fX4dHk8ZxGAyGbDFOo+iSk3trHIfBKlWqgFJ6q1Il7/oYDIbCh3EcBqtcuGD9ubP7GAzZcfr0aQYOHEjDhg2pX78+Tz31FPHx8Vbbnj17lvvuuy/bY/bt25eoqKgc2TN16lSmT5+ebbtSpUpluT8qKopPPvkkRzYUNMZxGAoUM0oxZIWIcM899zBo0CAOHTpEWFgY169f5+WXX76pbWJiItWqVWPlypXZHnf9+vX4+PjkgcX2YxyHE1BKPamUClVK7VNKvWejzR0pbQ4rpV7IbxsLM45+QVeubP25s/uYUUoR5OhRaNYMPDz049GjOT7Ub7/9hre3Nw8//DAA7u7uzJw5k3nz5hEdHc2CBQu4//77GTBgAL179+b48eP4+/sDEB0dzeDBgwkMDGTIkCG0a9fOIj1Up04dLl68yPHjx2nSpAljxoyhWbNm9O7dm5iYGADmzp1LmzZtaN68Offeey/R0dFZ2nrs2DE6dOhAmzZtePXVVy3vX79+nR49etCyZUsCAgL47rvvAHjhhRc4cuQIQUFBTJo0yWY7l0RECnwDugG/AF4prytZaeMOHAHqAcWA3UBTe47fqlUrudWBjJur4Kp2GdLYv3+/Yx2aNhVxc9M31M1Nv84hH374oTz99NM3vR8UFCS7d++W+fPnS/Xq1SUyMlJERI4dOybNmjUTEZFp06bJ2LFjRURkz5494u7uLtu3bxcRkdq1a0tERIQcO3ZM3N3dZdeuXSIicv/998v//vc/ERG5ePGi5Xwvv/yyfPTRRyIiMmXKFJk2bdpNNg0YMEAWLlwoIiIff/yxlCxZUkREEhIS5MqVKyIiEhERIfXr15fk5OQMtmbVLj+wdo+BHWLjO9VVRhzjgHdEJA5ARMKttGkLHBaRoyISDywFBuajjYY8ICejFIOLExoKycn6eXKyfp1DRMTqqp/07/fq1Yvy5cvf1OaPP/7ggQceAMDf35/AwECr56hbty5BQUEAtGrViuPHjwOwd+9eOnfuTEBAAEuWLGHfvn1Z2vrnn38ydOhQAB588MEMtr700ksEBgbSs2dPzpw5wwUrw2t727kCruI4/IDOSqmtSqnNSqk2VtpUB06le3065T2rKKXGKqV2KKV2REREONncwoerfkGfP5823jh/vqCtMTiFRo3ALeWrxc1Nv84hzZo1u0nZ+urVq5w6dYr69esDULJkSat99Y/m7PHy8rI8d3d3JzExEYBRo0bx8ccfs2fPHqZMmWJXroM1J7dkyRIiIiLYuXMnISEhVK5c2eqx7G3nCuSb41BK/aKU2mtlG4jOYC8HtAcmAcvVzXfA2mJjm38ZIjJHRFqLSOuKFa1Kyt9S5McXdFhkGAO+GkBYZFjenAATTC8UrF0LjRuDu7t+XLs2x4fq0aMH0dHRLFq0CICkpCQmTpzIqFGjKFGiRJZ9O3XqxPLlywHYv38/e/bscejc165do2rVqiQkJLBkyZJs23fs2JGlS5cCZGh/5coVKlWqhKenJxs3buTECa1WXrp0aa5du5ZtO1ck3xyHiPQUEX8r23fo0cM3KVNr24BkoEKmQ5wGaqZ7XQM4mz/W33o46gTCIsPo8GUH1h9eT4cvOzjUz5HzmGB6IaBePdi3DxIT9WO9ejk+lFKKb7/9lhUrVtCwYUP8/Pzw9vbm7bffzrbvE088QUREBIGBgbz77rsEBgZStmxZu8/9xhtv0K5dO3r16kXjxo2zbf/hhx8ye/Zs2rRpw5UrVyzvDx8+nB07dtC6dWuWLFliOZavry8dO3bE39+fSZMm2Wzniih7h3N5aoRSjwPVRGSyUsoP+BWoJemMU0p5AGFAD+AMsB0YJiJZTzwCrVu3lqJUyKlKlbQvzcqVnT+CSHUCUbFR+Hj7EPxoMH6+fsTGwunTcPFixq1tn6Pc/UMbLu1pDX8/BUrwdPegY82OlPIqxcyZ0KAB/PQTfPZZ2ojhRsI1Np74hYQez1CuyjX+12knN87UoVIlLFu5cmmzHqD7pccF/nyLPAcOHKBJkyYFbYbDJCUlkZCQgLe3N0eOHKFHjx6EhYVRrFixgjbN5bB2j5VSO0WktbX2riJyOA+Yp5TaC8QDI0VElFLVgC9EpK+IJCqlJgA/oVdYzbPHaRRF8upX99Wr8NXaMzy75Gtizs6AqDpciq5A67AX2fHufzn0tx/9+9/cr+3FBUSViAK3RIjW04IJKHYdOk+9cg1ISEg7/pEj+ss+LjGOo5fOkJTUFFQCl2Muc9/rS4j5KeP6fA8P7awqV4b//Q+8vSF12rd8ebh+HbLJszLcokRHR9OtWzcSEhIQET799FPjNJyES4w48pqiNuLIza/umBgIC4MDB9K2YcNg0CDYvRtSFpdAqbNQ/giUDEfd9gH9uvswp9tafvkFKlQAX1/9WKECnIsP47Z5HbgccxlBUCjKFS9nGalYY8BXA1h/eD3Jkpx2XXE+dC4znMltPiY8HMv22mvagbz3Hnz+OZw5A3Fxuk+xYtp5eHrCd99BVBQ0bao3GzFTg4MU1hGHwX4K64jD4ACVK2ecqrJFXBzs2KG/XNu0gcuX9Rd+qqNxc4O6daFPH/26SRNY9uNJHvu7C1c4kcEJzOgTTFVfSLfK0EIZ/Ah+NNjq9JYtZvSZwV+n/8robHzcmPvof/Dztd7n+edh0KNhPLN+Ev/x+4CYc3U5d047DYBPPoENG9La16kD3bvDl1/q1zExULy47c/LYDDYhxlxFDE2boRNm2DLFvj7bz2tc/fd8M7cMCb+NJFGB+fRpllFmjQBPz899ZMZWzGO7AiL1OeY0WeG3e0dOU927RMT9VTYvn1pm4+PjquATmIG6NpVb1266DiKIWvMiKPo4+iIwziOQkxUFPz5p44BPPaYfq9tW9i5E1q0gNtv11uVpofpt7qdQ47AUSeQU+w9T6rTcGQ6LD3JyTBtmnaqv/8ON27o959/Ht59V4/CIiKMI7GGcRxFH+M4rFCUHMc//8CSJfDbbzomIQJly+rVTR4eOkm3alUoU0a3z+0XrqtgLSbiptzo26Ava4dlnSeQeRXaqVP6c9y0STvY3r3h0CE9AmvaNG1E0ru3/mxvdYzjKPo46jhcJXPcYIPERPj1V0jVV/v+e5g9W8cqpk7VU1PnzmmnATpJN9VpAEz8aSJRsVFISq6kIETFRjHxp4m2T5oTkTonCttZY0afGfh4+6BS8kAVCh9vH2b0mZFt38yr0Dw9oV07+L//084BoHRpeOcdqFULFi2CwYP16OPXX516GYYckJ08eVFg1KhRFlXf0aNHs3//fgC78lUyk17oMa8wjsMFiYvTDuKRR/Qv5J49dQ4EwIQJeqXRL7/A5Mn6l3FWAd8cfeEOGAAHD0JSkn4cMCB7o3PSxwH8fHUAvlzxcrgpN6ePmqpU0Y7khx/g0iX44w/9WbdqpffPng39+8PChXqK0OD6pEqH5IakpCQnWOIYX3zxBU2bNgVy5jjyA+M4XIwTJ6BiRf0ltWoV3HknfPNN2sqncuUyjiiyw8/Xj+DeyykX54ZbMpSLcyO49/Ksv3BzIlLnRGE7W6Q6j74N+jrkNBzV6fL0hI4dYcYMHVwHvQR6zx4YNUqPRPr10yMTQ8Gwdu1a2rVrR4sWLejZs6dFDHDq1KmMHTuW3r1789BDDxEREUGvXr1o2bIljz32GLVr1+bixYsALF68mLZt2xIUFMRjjz1mcRKlSpVi8uTJtGvXjuDg4AzntSW1PmrUKMaNG0e3bt2oV68emzdv5pFHHqFJkyaMGjXK0r9UqVJMnDiRli1b0qNHD6zp6HXt2pUdO3bwwgsvEBMTQ1BQEMOHD79pJDF9+nSmTp0KwM6dO2nevDkdOnRg9uzZljZJSUlMmjSJNm3aEBgYyOeff577Dx9cQ1Y9rzdXlVWPjxcpWzZNUrxyZZHkZJGJE0W+/14kNtZJJ2raVEIrKOk/FAmtoLKXuc6JLHZO+hw5otu5u+vHI0fsux4HCb0YKv2X9JfQi6G5Ok5yssjff4s895xI7doi3bql7Vu/XiRF2bvI4bCsuui/5fR/17khVZ48PZcuXbJIjs+dO1eeffZZEdGS5y1btpTo6GgRERk/fry8/fbbIiLyww8/CCARERGyf/9+6d+/v8THx4uIyLhx4yyS6IAsW7bMqi22pNZHjhwpQ4YMkeTkZFm9erWULl1a/v33X0lKSpKWLVtaZNsBWbx4sYiIvPbaazJ+/HhL/xUrVoiISJcuXSzy7+mvPbMM+7Rp02TKlCkiIhIQECCbNm0SEZHnnnvO0u7zzz+XN954Q0REYmNjpVWrVnL06NGbrquwyqrfUpw8qaeZateGdJI2XLigf9lOnw59+0I60c4MOCwmGBqK30Vh7dfgd1GyHw3kRKQuJ33yeHoLcq6hZQ2ldGxk2jQ4dgxWrNDvR0Zq06tWhREj4K+/jBRKXmuKnT59mj59+hAQEMC0adMySJ7fddddFE+Zv00vrX7HHXdQrlw5AH799Vd27txJmzZtCAoK4tdff+VoSlzO3d2de++91+p5s5JaHzBgAEopAgICqFy5MgEBAbi5udGsWTOLVLubmxtDhgwBYMSIEfzxxx+5/iyuXLlCVFQUXbp0ATJKum/YsIFFixYRFBREu3btiIyM5NChQ7k+p0kAzCeSknTCnVI6CPvZZ3oa6tw5x46TPpfhr9N/2Tdl06iR/mJOTrZP5jpVpM4RctInj6e3Mq8ouxxzmQ5fdnBKbEQpvUABtPTJ33/r+MeiRXrVW2AgfPGFTrw0OJ8nn3ySZ599lrvuuotNmzZZpmwgo8y62PDgIsLIkSP573//e9M+b29v3N3drfYbNWoUq1evpnnz5ixYsIBNmzZZ9qXKs7u5uWWQandzc7MZb7Emw24LDw8PkpPTVhWmSq6LjZolqftmzZpFn9S5bidhRhx5zPnz8PbbUL++/iUK8OKLetHR9987Nv+e/oswWZItX4TZ/op2osy1U3Fi3QZr5GhFWQqOyLcrBa1bw6xZWg7l88/14rKqVfX+3bshZZHMLUNe13+5cuUK1avrcjwLFy602S69tPqGDRu4fPkyoOXaV65cSXi4rhl36dIlu2TMHZVaz0xycrJl9dRXX31Fp06dsmzv6elJQorYW+XKlQkPDycyMpK4uDjWrVsHgI+PD2XLlrWMXtLb1adPHz799FPLMcLCwriRmsSUC4zjyANEdJ7F4MFQsya8/LJWh039EVOzppbDAMfqZOT4i9CJMtdOJY8dmjOX8NpLqVIwdqxOwqxRQ7/30kt6hXK3bnp6K1X0sSjjzPov0dHR1KhRw7K9//77TJ06lfvvv5/OnTtToULmCgxpTJkyhQ0bNtCyZUt++OEHqlatSunSpWnatClvvvkmvXv3JjAwkF69enHOjuG/o1LrmSlZsiT79u2jVatW/Pbbb0yePDnL9mPHjiUwMJDhw4fj6elpCdr3798/w/nnz5/P+PHj6dChg2WaDvTS3qZNm9KyZUv8/f157LHHnLLarMAD1/mx5VdwPLU8cEyMiK+vSPnyOtAdmruYrIXQi6FS/m0fUZMRpiJqMlL+bZ9cB30LDTkIpodeDJXy75YXt9fcpPy75e3+rJxZCz08XOSdd0Tq1NHHqlpV5OOPc3fM/CQnwXFXITY2VhISEkRE5K+//pLmzZsXqD3WAv2ugAmOFwCXLulkvNtu07EMb28ttnfmjA50+zkpQdvP14/g1b6Ui0UvrY2F4NW+hSoDPFfkIJieX0t4s6JiRZ0jcvgwrFuns9VTpeETEhwPDRns5+TJk5bls//5z3+YO3duQZtUNLDlUYrSllcjjvBwkRdeECldWv+SHDRIJN1qvbzB3V1CfdFLa33Rv75vFdzdMw4DCvG1p45O58/Xl3LvvSK7dxeoSTYpzCMOg32YEUc+ERKi4xTvvqtXR+3eDd9+m7bSJs9o1Ai/y256ae1l5weUXZo8DqanJ6/rp6cugrnrLnj1Vfj5Z2jeXCsZ79qVJ6c0GJyGWY7rAKdO6cJHvXtDQAA88QQ8+qiO6+Yba9fqKZrQUP3Fmd8rpOLjdeJC5vqxFy/q9+Pi0pYjpWxxKpmz7jeoKqXwVp4cdr/COs9jPJTQlPLupdjuc4N5xfYxpf6jVKneiG3qLIvP/cSUHm/gW8KXv0//zYp9K5i86ivK3juMXVcOsqV1RR59dymlgKjYKOKT4qlYoqJDyxttkaMlzzmkfHl4/XV45hn48EP44AMtDb97980FuwwGl8HWUKQobbmdqjpyRGTMGBFPTx3YTIm1FT2io0VCQkS+/lpk8mSRhx4S6dtXpG1bkXr1RMqUkZuixiAXSiKftkZOVispUrGibGtaVgImuMv2xmVEypWTdUElhanIVr+SIqVLy3fNvYWpyI663iLFismqJkjFScih8vp4X/kjPv+HHK9WQqRePZl3bz0p9aq7nB07TOSll2T6e3cLU5Gok4dEROS9P94TpiJXY6+KiMiX/3wp3RZ0k/hEnRX8x4k/ZPHuxZKYlJjtR5AaTFdTlV6AMFU5FFTPLZcvi+zdq59HRYkMGSKSkkRcYJipqqKPo1NVZsSRBUeP6l+DixfrFaNjxuj6DR6F/VOLispYOzZ1O3YsLeVZKb2etFIlkiv4cqJJVUqVq0RF31qc8IFh0Yt5pf4o7vTrx3kVzrhVPVl675cM8R9C6YsHqffLC3i+/BpUaU6La2eZf2QDtSfeCaUq0zsxlksJMZSdXBaUG/fcuME94eFwr64VOzRlY3A4RETwcHg4D68pDeEbITycZ5KTeNgbykxtCBUq0LttbUo2uZ1SXyyCZs1QbtdxU254uuvSgCv2r2BByAKGBQwD4L0/32P3hd0suUevdw+/EU5Zr7J4eXhlueQ5K/n2zNLtOV2C6uOTpo+1Z49eZLFsmVYSmDJF11sxGAocWx6lKG05HXFs2SJSvLjI00+LnDmTo0MUPKdOiSxZIjJ+vBZXqlIl44jBy0skIEBk8GCRKVNEli6V2H+2y/TN/5UfD/0oIiIXb1wUpiLT/5wuIiJXY6/K7fNvl/Vh60VEJD4xXk5EnZCk5KS8v56kJH0zNmwQmTlTZPRokQ4dbh4NVawo0rWryPjxEjv7Qzm66VvLUPGNzW/IvcvutRxywFcDJOCTABHRI45Sb5cSpuLQiMOZy3fTc+WKyFtv6aXdINKnj8i1a847vj0U9IjDlZawzp8/X86k+zJ49NFHZd++fSIiUrt2bYmIiJDLly/L7NmzC8rEHOHoiKPAv9TzY8vNVFVeC9c5S4BPRPRSnUOHRL78UmTkSJG6ddO+yUqVEmnXTmTUKJF33xVZs0a3TdTTN+/+8a7M+2eeiIgkJSeJ77u+8uIvL1oOvTBkoWvniyQnayf5448i778v8uijIu3b6+tO/Qzc3LQzee89kW3bLI5kXeg6+XrP15ZD+c3yE8/XPS25HxN/mmhxorbIK8eRytWrOhdk5Mi092JinH8ea+TEcTjz7zqnjiPBCXPKiYkZpzfTCxBmJtVxZBYjLAwYx+Fkx5GX5DQ5zUJSksiePSKzZ+vJ8KpV0765fH1F7r5b/yrfudPiIFL5IPgDeeXXVyyvb/vyNnnwmwctr6NionJzaa5DkyYiSqV9LsWKpT0vXVrHcDI5kuj4aAk+FSz9l/SXfeH7pPK0yjJpwyQREUlOTpanfnhK/jz5Z4bTOFMN1h4OH9aDqunTReLi8vZcjjqOXP9dZ8Ka41izZo20bdtWgoKCpEePHnL+/HkR0eq4Y8aMkV69esnQoUMlPDxcevbsKS1atJCxY8dKrVq1JCIiQkRE/ve//0mbNm2kefPmMnbsWIuTKFmypLz66qvStm1b+f333y3nXLFihZQsWVL8/PykefPmEh0dncGRpDqOIUOGiLe3tzRv3lyee+65XF17fmEcRwE6Dke+PHKcBX7ihMiHH4oMHJg2fwEi1auLDBsm8umnIvv2pSUKpPDp9k9l8IrBltejvxstdy6+0/I6NZBc5LCW+3HunMiyZSLjxmnHkt6R3HmnHpFt3aods4gkJiXKtTg9P3TqyinxecdHvtj5hYiIREZHypg1Y2TPhT35ellHj2qfByING4qsW5d353LEceTF4gJXklXPPOKw5jhuhRFHYQ/zuhSO6BtZgrApaQnihu0g7KFDuprTqlWwfbt+r359GDQIbr9db3XqZFi/+e+Ff5m7cy4z+sygmHsxEpMTCb8RTnxSPMXcizFnwJwMS1dTA8lFDmvKwFWqaCGxwYN1m/PnYcsWXYR80yad5g1apfDee3G//35KdewIQI0yNbg46SKJyVrvJ/RiKMv3Ladb3W68+MuLTGg7gf0R+3m4xcP4ePvk2WXVratFMtevh2ef1YW/BgyA1avTUl0KgpwuLnCU06dPM2TIEM6dO0d8fDx169a17Mssq/7tt98CtmXVAWJiYqhUqRKQtay6QWMSAAuIGX1m4BMLKkUlWSWDT6x+HxHYuxdee03rc/v5wQsv6Pf/+18IC9P6FV9+CSNHQt26xCbFsSZ0DReua491POo4C3Yv4EDEAQAmtJ3AxpEbKeZeTJ/vVkkSsEdIMdWRfPKJlrE9f15ro3fooLXRu3TRK8zGj4dNm3AX8PLQstkdamqZ9vHfj2f94fXcs/wenvv5OeKT4gE4FHmIU1dOWTXNEQVeW/TtC//+C++/r/9UUp1Gao36/CY3wpKO8OSTTzJhwgT27NnD559/bpEYB8dk1UNCQggJCSE0NNQizZ6VrLpBYxyHE3FE38jP14/gTfXTdKdi4J911fGbPl9/wQUEaMdRtizMnKlrym7frh1Iw4YAXIu7RvgNLQt97PIxBi4dyKoDqwC4s8GdREyKoHmV5nlyrYWGnCgDV64Mw4bpEV5EBCxdqmvJzp+vJW6rV9fZnxs3EhZ+gE7zOxEVG0WyJBOTEEMZrzJExUYBMHnTZFrNaUVSsi5LGpuY9gXnrGJHxYrpBMI339Svt2zRRcI++0zLeuUneV0bPpWCklUvXbo0165dy3WbQo+tOayitLlqcFyOHJHjzWrIt42Q+GIeaXPwPXvqWMW5cza7xiXGSbl3ysl/1v/H8t7m45uLbqzCFbh2TcdG7rtPr9MGuVSmmHzaGun2EOL+qo5Xub3mJv2X9BcRkUORh+T7sO8th2g9p7WM/m60iGQMvThzJdbevSJduuhjBgSI/PZb7o5X0KuqlFJSvXp1yzZjxgxZvXq11K1bVzp16iTPPfecdOnSRUR0jGPatGmWvhcuXJDu3btLixYt5Omnn5aqVatKbEpN5qVLl0rz5s0lICBAWrZsKcHBwSKS9SqulStXZhscFxEZOnSoNGvWzATHC/Pmco7j+nW9ZLZNG7Gs9OnXT2TevCxVEt/a8paM/Hak5fXcnXPl71N/54PBhpu4fl1k+XK5OvBOue6pv/nPlUTe7IwEvFDW6hdmUnKSvPP7O7J4t645XalqvDC0v1Bno9NXYiUni6xcqWujg8h//pNtF5sUdB5HbnA1WXVXxTgOV3Yce/aITJggUras/uibNROZNUvrTFghLjFO1oWus6wembJxijyw8gG7pDMM+UfYyRAZOaKUrGmEJCok2c1N5K67dE5Jku2kyOOXj0v9D+vL2tC1IiJyJfaKnIg64VTboqNFXnxR5JNPcn6Mwuw4wsLCJCgoSAIDA6V169aybdu2gjbJJTGOw9UcR0yMyOLFIp06iWV0MXy4yO+/37RkNjNzd87VGk+nt+aTsQYRyXHBqP5L+suRXb+JvPSSTrIAkfr1dbKFjZFkUnKSJeN+ZvBMcXvNTY5eOmrzPLnNF1m0SGTECMcSWwuz4zDYh3EcruI4wsJEnntOJ+KBSIMGItOmiaTMgVoj/Hq43LH4Dlmxb4WIaGmP9WHrzQgjv2naVGeZp2abN23q+DFiY0W++irtB4OXl0773rrV8oMhcxzgRNQJmb0tTapixl8z5MO/P8xw2NzGRaZNE/Hw0E7nm2/s67N//37LqNdQ9EhOTi68jgN4EggF9gHvWdlfE9gIHEhp85S9x843x5GQoCeWe/bUH62Hhw6k/vKLzSmLs1fPWrKQE5MSpeOXHWVRyKL8sddgHWcXjPr3X51smCp/0qqVnPvwLan+Rrkss6vv+vquDEmbRy4dcUpAfdcukRYtdP8hQ3RBsqw4evSoREREGOdRBElOTpaIiAg5evTmUW5WjkPp/QWLUqob8DLQT0TilFKVRCQ8U5uqQFUR+UcpVRrYCQwSkf3ZHb9169ayY8eOPLEd0PU/Fy+Gt9/W+RW1asHYsfDIIzqJLAtun3874TfCOTD+wK2TW+HqNGuWMWmwcWPn1He9ehUWLyZu1gd4HTzEZW+Y2xJm3AaJFctbXbYalxiHl4cX566do9YHtfD+4x2ub5gI5E6FNyEB3ntPqz+vWKELStlum8Dp06cz5EoYig7e3t7UqFEDT8+MScBKqZ0i0tpqJ1seJT83YDnQ08E+3wG97GmbZyOO2FiRzz4TqVNH/3xr0UJk1aqbdKHSc/zycXnqh6fkRvwNERH55+w/cijyUN7YZ8gZOYhxOEL/xf2ky8NKljbTwfQbHsiMDsiI2T1t9rkWd00+3vqxHI48LCJ6mmvpnqXZTmNmFxM5dSrt+apVWa4AN9xi4OpTVUAI8BqwFdgMtMmmfR3gJFAmizZjgR3Ajlq1ajn3E42O1npR1avrj7BdOy0WlMVQPnWYv+X4FvF6w0t+O5rLxfWGQkt6PSe/CcjCQO1Akry9tIb/2bPZHuO5n56TEm+VkIgbtmNmImL31FZUlFamL1dOB9DNrJTBJRwH8Auw18o2MOXxI0ABbYFjoKfRrBynFHqa6h57z+20Ece1azq6mPoz7vbbRX7+Ocv/soSkBBmyYoi8vul1y3vZ/bMbij6ZFWSPbv1JB8/d3UW8vXXixenTNvsnJSfJ7vO7La+fXP+krNy38qZ2jsREDhzQpU1ArybOIqXIcAuQlePIN8kREekpIv5Wtu+A00DqGo9tQDJQIfMxlFKewCpgiYh8k1+2c+WKjl/UqQOTJmk5kM2b9dazp0VcMCwyjAFfDSAsMoxz184B4OHmgZeHVwYRwQolbro0wy1GqjRH3wZ9CX40mLpte8OCBbqW/LBhMHu2FrKcMAFOn7b0S/0bO3zpMIGVAwG4Hn+dzSc2ExoZCugfg6kSJ47I4DRuDL//DqVKwZo1UKFC9n0Mtyi2PEp+bsDjwOspz/2AU2QacaBHI4uADxw9fo5HHJGRuva2j4/+Gda3r8hff1ltGrrrFyn/oru4TUZKvKLE+w0vOX/tfM7Oayhc5EVM5OhRXd3Qw0Pn/owbJ0dCNtqsc5GUnCRxibowx69Hf5UGHzWQAxEHcnRqa6MUM3V164ErjDiyYR5QTym1F1gKjBQRUUpVU0qtT2nTEXgQ6K6UCknZ+uaZRVeu6F98r7+uhe127NA61h063NQ0LDKMtst7c8kziWQ3iHYTVFw8Z66dyTPzDC7EgAF6FVZSkn4cMCD3x6xbF+bO1ZL6o0YhX3xBjVbdeHv5JapcSeZyzGU6fNmBsMgwANyUm0X52NPNEz9fP+r41AHg5JWTlhFITvjlF63CmxshRkPRwiWW4+Y1OV6OO3u2rnUREJBlszv+dwc/HfkJ0q2mdUuGvo36O7UGgcFF8fDIKEPr7q7VeJ3IIx/1pM3iX3nkH0h0g/92gpkdFd2b9Mvyb0xEaDmnJbXL1mb1A6vtOleVKmlOonJlmDEDRo8GHx+96rxHj9xfj8H1yWo5rquMOFySKm+MRwUGWK2XkJicyK9HfwXgo74fUSrRTUdmSKmtkeDu9BoEBhelUaO0QhipxaKczAvDP+GV+8rTdDz82ADe3AgHZis+j+ulZ5RsEBYZhhtu9G2oB+dJyUk264Okcv582kTV+fMwfDhs2wblykGvXjB5stP9oqGQYRxHFmRVL+GDvz+g5/96si98H36+fuy8bwPlE9x1bY0Ed4IH/+T0GgQGF8WeYlG5JDWYHlW9PIMfcGPg2NJUrtqAao88padSQ0Ju6hMWGcZt824j5EIIL/76ImGRYcwPmY/fx37sDd/r0PkDAnQ5mFGj4I03YOXKjPudUZTKUIiwFfwoSltOg+OZg4QHIg7IvvB9IqKVTL/Z/00GGQZn1iAwGKyR4W8sIUHXbfH11ZpaY8da9ENs1f7edGyTTN041fJ3e/TSUYelRH75JS1Yniq9lle1RQwFB64uOZLX5DTGkX6ut1LVBDyeq0NQlSC+H/a9ky00GHLB5cu6WuTs2VCyJEyZwt2+v7Dm+I8kS7KlmZtyo2+DvpaYyLW4azSc1ZB7m9zL7H6zHT7tiRMQFKTjH9OnZ9x3C3ytFHlMjCOHnDqTwNI9y0hOFi6c9eTre79m/sD5BW2WoShw9KjWxPLw0I9Hj+b8WOXKwQcf6OLj7dvDs8+y7I0D3Hu8RJa1v0sWK8nbPd7m4RYPAzof5Pz18xnykbIitcLu9On6MtK/byjamBFHFvxv9/94aPVD/PbQb3Sr2y0PLDPcsuSVkKIIrF+vi5AfOsTPjTx5qlciF2pmX/v75V9f5qNtH+Hh5sHVuKv4ePvYVS981Sp49FFd+/ybb6BTp9xfhqHgMSOOHDI0YCgbRmyga52uBW2KoagRGqqdBujH0FDnHFcp6NcP9u6FGTPoftaL3Z9C6KXh+JWpm2XXLnW6kCzJXIm9QrIkcynmUoZcEVvcey9s3aoHPvPmOecyDK6NcRxZ4OHmQa/6vYzcucH55PUS3mLF4NlncT98BM/BD1Dh3VnQpg3s2mWzy6yts4hNjEVIm4W4HHOZiT9NzPZ0jRpp5/HJJ/r12bNZL9k1q7AKN8ZxGAwFQT4s4QWgUiX46itYvVqv9GjTBl55BeLibmo6o88MfLx9LHERgLLeZZnRZ0aGILstfHzA2xtiY6F7d7jzTrh0yXrbrJa6G1wf4zgMhoKgXj0d00hM1I/16uXt+QYOhP37YcQIeOstaNVKJ2akIzVXpFzxcrgpN8oXL8/W0Vvx8/VjxDcjshx5pA+oe3vD889rDdB27eDAgby9NEP+YxyHwVBYyO1KrHLltALv999rLbb27eH//g9iYixNMqv2+vn6kZScROWSlfEt4Wv1sGGRYXT4sgPrD6+3xEQeeQQ2btRFD9u31/H69Dii2mtwPcyqKoOhsODMlVhXrugSAXPn6gDFvHlw2212df3r1F8s2r2Iab2mce76OTp82YHLMZcRBIWiXPG0FVwnT+rBjrc3/PlnWljH4PqYVVUGQ1HAmSuxypaFOXPg5591UKJTJ72ENzo6267bzmzj56M/AzDxp4lExUZZAuqCEBUbZZnWqlUL/vhDh1jc3PQIxJQuL/wYx2EwFBbyYiVWz56wZw888YROIgwM1MGJLHi6/dPsGbeH0l6lmdZ7Gl7uXpZ9VhMNS+rpKBEYOlRLa50/n3vTDQWHcRwGQ2Ehr1ZilS4NH38Mmzbp1127wlNPZTk0KOFZAoBLMZeIT4qnpGdJ3JRbhmmqzCgFjzyiE9xbtErg9reezjZHxOCi2BKxKkpbjkQOjxyR0Lb1pf8wJLRtfedUdTMY8pOcVCa8fl3XOweRFi1EwsKy7XIw4qAcjDgo/Zf0l/Vh6yUhKSHL9qs3Hhc3n5OC5w0pNXKoEQV1UTAih44Hx8PaNaBD1yNEeYNPLARvqo/f1sN5ZKHBkAfkJpi+bh2MHAnx8fD551qUKhuuxF6hwawGDGo0iLl3zbXaJnUF1qUIT/hqLVyvQrnn2/P34xuzlTbJXGDKTHflLSY47iBhkWF06HqEy96Q7AaXvaFD1yNmWG0oXOQmmN6/v67xERSkKzmNHp1t4Lysd1lm3TmLibfpwLi1H6WpwXRKXYCR3eDBXlxJOs+zP07MVlHXJA26DsZxWGHiTxOJ8gZJ+XTEDaK8sUt6wWBwGXIbTK9ZUydjvPyyXq7btm22I5YH/B+gcYXGADy34Tme+fGZDA4kQ3a61w1UxTB8vH2o/NdCnn46zc8ZXBvjOKxg+eNOXwo200oRg8HlyUkwPXOS4cmT8Oab8NNPEBGhJUvmzcu24IaIkJicSGJyYgatt8zZ6eWKl+OvR4Lx8S7PRx/pGTEraiiASRp0KWwFP4rSlpPgeGoFNbfX3KT8u+VNAM9wa9C0qa4mCPqxadO0fefOifToofcNGyZy9Wq2h0tKThIRkUORh2T53uWW961Vy3zvPX3oHj1Erlxx3iUZcgZZBMcL/Es9P7aclo41pWANtxzu7pKhBqy7e8b9iYkib76pnUqDBiL//GPXYcesGSPl3iknl6IvZdlu4UJ9yg4dRJKScnoRBmeQleMwU1VZ4Ofrx9pha7Nd7WEwFBmyi4u4u+uYx6ZNWuOqfXudAyJZT13N7jubzaM2U654OQDOXTtntd1DD8GaNfDsszfLk9hbmdCQ9xjHYTAY0rA3LtK5s1511asXPPmkljBxd7cpvujp7klA5QAAvtrzFQ1nNWTXOeu1Qfr2hfvu08+XL4d//rEupGgoOEweh8FgyDnJyVCtWtr6WKWgSZMsV1+duXqGmX/P5J2e7+Dh5mGzXVwc+PvD2XPJuA+7l+vVv7MqpGjIG0weh8FgyBvc3ODixbTXIjrpMAuql6nO9N7TLbXNn/j+CS7F3FzxyctLz4i5lzvJtXnLkL16GCKZhBQN+Y9xHAaDIXekj4uAHoW8/362cQ+Arae3snD3QvZc2GN1f/Xq8OumRDxq7YCVS2HbE1aFFA35i3EcBoMhd2SOi/TpAxMnakVDW0kZKfSq34vjTx2nS50uAJy/frOOSJv6Ddi5pSKeTX+EazUcmqYytc3zBuM4DAZD7khfBvfAAV3ub8oUXW2we/ds9UEqlqwIwPYz26n3YT1W7l95U5vAGg3ZvbEB/R77i+BHgykv9sU2jExJ3mAch8FgcC5ubjB1ql4StWsXtG6tH1OxUQK3acWmPNbqMbrX7W71sE0q+bFu+FqKXfOjaVN45518uBaDVVzGcSilnlRKhSql9iml3suinbtSapdSal1+2mcwGBzk/vt1vViloGNHWLFCvz9ggA6gJyXpxwEDAChZrCQz75hJ+eLlSZZkXvntFatTVzVq6FXAL76o1VCywsiU5A2218LlI0qpbsBAIFBE4pRSlbJo/hRwACiTL8YZDIac06IFbN8O99wDgwfD5MlpUu9gU7V3X/g+Zv49k9plazOm1ZgM+zw8YNEiHVJ59VXtf6ZMsX56I72eN7iE4wDGAe+ISByAiIRba6SUqgH0A94Cns0/8wwGQ46pXBl++w3GjYPXX9cVB69f16uubKj2BlQO4MD4A9QsUxOAqNgofLx9LPvd3WH+fP04dSrUqaPLhxjyB1eZqvIDOiultiqlNiul2tho9wHwPJCt+LJSaqxSaodSakdERIQTTTUYDA7j5QVffqmX6d64AcWKpRWXspGdXqtsLZRSnLl6hiazmzBr66wM+93d0w55//3WT2tkSvIGhx2HUqqkUso9B/1+UUrttbINRI98ygHtgUnAcpVei1n37w+Ei8hOe84nInNEpLWItK5YsaKj5hoMBmejFDzzjF515e0Nvr562FCvXpbdKpSowH1N7qNHvR437XNz04csUQKuXIHZs9PSR4xMSd6RreNQSrkppYYppb5XSoUDB4FzKUHsaUqphvacSER6ioi/le074DTwTYoo4zb0iKJCpkN0BO5SSh0HlgLdlVKLHbhWg8HgCvTpA1u3QqlSernuTz9l2dzLw4tZfWfRtGJTABb/u1hXEczEF1/AhAkwaRKEXtRO43LMZZIlmcsxl43zcCL2jDg2AvWBF4EqIlJTRCoBnYG/gXeUUiNyacdqoDuAUsoPKAZcTN9ARF4UkRoiUgd4APhNRHJ7XoPBUBA0agR//QUNGugytYvt+w147PIxHvnuEWYGz7xp37PPascxYwbc+eABLsdEIejhh5EpcS72BMd7ikhC5jdF5BKwClillPLMpR3zgHlKqb1APDBSREQpVQ34QkT65vL4BoPB1ahSBTZvhkGD4MEHdYbeRBtf7EePwoAB1A0N5c/2tQgaOuymJkrBRx/p2MeHHw7E6+qnxPV6HJQYmRInk+2II73TUEp1UkqNV0rVS/deXWuOxRFEJF5ERqRMXbUUkd9S3j9rzWmIyCYR6Z+bcxoMBhegbFn48Ueto/7cc3qeyVrh8XS5H22CT+A58B6uxV2j75K+bD+z3dJMKZg5U48+fE8+jE9yQ0uJWqOm6zwcDY5XBNoCU5RS7yulWgNvO98sg8Fwy+DlBUuXwvjxMH06jBoFCZl+i4aG3pT7ERkTyeFLh7lwI6OWiFL6MLt3ebL1qbXcWb8vfz1inIYzcSiPQ0S+VUqtAVoB/ug4h4k2GQyG3OHuDrNm6emrV1/VUu0rVkDJknp/o0ZpiYMpuR91fOqw94m9FHMvBuiqglVLVwW086hQAcon+1Hj97UsPghvvFFQF1f0sGdV1atKKcvEo4gkicg2EZknIjNFxEbOpsFgMDiAUvDKKzB3rl5p1b17Wq0PG5UJU53GrnO7aDCrAV/t+eqmQyYmammSGSa84TTsGXE8CARlflMpNRqoKCL/dbZRBoPhFmb0aKhYER54ADp10k4kVYHXBk0qNuHxVo/Tq16vDO8rBZ9/Dlev6hCKjw88+qjtU1epkqaiW7mykSyxhT0xjhgRibby/v8AsxzWYDA4n4ED4eef9bf4bbfBHuuFnlLx9vBmRp8ZVCxZkWRJZu7OuSQk6TiJu7te7XvHHTB2bJrWojWMDLt92OU4lFJVM7+ZoiuV6HyTDAaDAT3a+P13/bxzZ9iyxa5uvx37jbHrxmao61GsGKxaBT16pIVNDDnHHscxA/hOKVU7/ZspCrbZakYZDAZDjvH314mCVapA796wenW2XXrW68kfD//BA/4PZHi/RAk969U3ZYF/VNTNfY0Mu33Yk8exApgN7FRKrVNKvamUehv4E5ie1wYaDIZbnNq1dV2PoCCd77F8ebZdOtbqiFKKE1EnGLh0IBE3tNBpqgLesmU6bJK+vhTAln1h9F8ygNCLYSa+kQV25XGIyEKgLrAc8ARigaEisiQPbTMYDAaNr6+OeXToAEOHwtdf29XtWNQxtp/ZzumrpzO836GDlsrq0yetHIgRRbQfJalSkrYaKNUW6ASEAJtFJCnlfR8gUUSu57GNuaZ169ayY8eOgjbDYDDkluvXtbbV77/DwoUwIov1OSkyJTFHDlK8vl7CG1+7hmUJb1iYDqN4e8OitUe596c2XI65jKAlSm71bHOl1E4RaW1tnz0jjgeBRkALYJFS6hWlVGlAAV84z0yDwWDIhlKl4PvvoUsXeOgh7TxskSJTUjwuGQ4eZNGEzrT4vAUXruvlUn5+OuZx9Sr0v7MYl68kGFFEO7HHcfwfEAr0AaoAPdHS6q/Z2d9gMBicR8mSsG6dXiL18MMwb571dplkSuoeOE/jCo0p45VWdbpFC+2HnhjrRbmynih0EMSIImZNtlNVGRor5QU0RGtWVQNaiojLu2QzVWUwFEFiYrSy7oYNOtt89OiM+5s1yyhT0rixJYkwNjGWyzGXLRIloGMcbd8eydXieyhXxuuWnqaCrKeqss0cV0opSfEuKbkbe9PtW5q5jcFgMOQLxYvDd9/BPffAmDGQlASPPZa2f+1aPV0VGqq1rtKVqB21ehQ7z+1kz7g9eHt4A1DZww+14A8q1Ajhtx9L39JOIzvsKuSklHpSKVUr/ZtKqWJAF6XUQsCUiTcYDPmPtzd8+y306wePP65rx6aSKlOSmKgf05Wofab9M7zc+WWL0wCt8D7tXXci/m3F9Bf9MD+FbWOPVtUdwCPA10qpukAUUBztdDYAM0UkJK8MNBgMhizx8tJp4YMH6xKAycnw5JNZdmlXox3tarQDtECij7cPdcvVZfRoOHMGpk7Vfmby5HywvxDiaIzDE10LPEZEovLKKGdjYhwGwy1AfLwWRvz2W3j/fXjmmWy7JCYn0mR2E6qVrsamkZtQSiGiY+4LF+oYfL9++WC7C5KrGEd6Uir9nXOKVQaDweBMihXTKeFDh+oSgElJWhI3CzzcPFhx/wrKepVFpaSVKwVz5kCTJtCzZ34YXviw23EopZ618vYVYKeZqjIYDC6Bp6fOKh8xQpehTUqC//u/LLsEVQmyPJ/x1wwGNBqAn6+fpdulS7osiJ+JlVtwJA+jNfA4UD1lGwt0BeYqpZ53vmkGg8GQAzw9YckSPfJ44QWYNs2ubuE3wnn3z3f54p+Mec333af1FY12VRqOTFX5ovM2rgMopaYAK4HbgZ3Ae843z2AwGHKAhwcsWgQi8PzzesnU2LFZdqlUshI7xu6gRpkaGd6fNg1uv10rnWzebGTZwbERRy0gPt3rBKC2iMQAcU61ymAwGHJLqvPo21cv1V26NNsutcrWwk25cTH6IsO/GU7EjQhatdKhk1279CAmKSljnypVdFxEKf38VsARx/EV8LdSaopSairwF7BUKVUS2J8XxhkMBkOu8PTUJf86d4YHH4T16223PXpUZ5t7eBB2Rxt+CP2eveE637l/f5g1S+cQvv56xm63YtVAR5fjtgI6ogUO/xCRnXllmDMxy3ENhlucq1ehe3edCPjTT3ruKTOZJEqu+jekzO6DGZp8/LGOeaQfWaTW+EilqCQO5kodVyl1TSl1VSl1FdgIvAW8CWxKec9gMBhcmzJl4McfoU4dPXzYaeU3byZRxDL7DgOwNnQt96+4n4SkBCZM0E4jMRFSf4v6+oXB0AHgG3bLVA20pwJgaREpk24rnW4rk11/g8FgcAkqVNDFoMqXhzvugAMHMu5v1EiLIYJ+bNQIgLPXznIi6gTX49NKD73yip79WvbjSeTRDrg1Xk/55zuwZd+tUfzJoamqwoqZqjIYDBYOH9YVnDw84I8/9CgELIWfMogipuhbJSQl4OnuSbIko1BcvKho3S6eU+FXkUfbQ/kjRa74U24LORkMBkPRoUEDLcV+4wb06pWWoJGFKKKnuyeJyYmMXD2Sl397mYoVod74J5FkBYvXQ3T5W6r4k3EcBoPh1iMwUK+wOntWZ/ddvpxtF3flTknPkpT0LImI8PmoiZQeNQKu1IJVX91SxZ/MVJXBYLh1+flnHSxv1UqPQkqVyrK5iFg0rWISYjh19RStnnmTGyX2U67esSIzTQWFZKoqpeZHqFJqn1LKaha6UspHKbVSKXVQKXVAKdUhv+00GAxFiF69tLbV1q26IFRc1rnMqU7j2OVjNJndhJDzIeyc+Qr9bq9K8KPBEFk0nEZ2uITjUEp1AwYCgSLSDJhuo+mHwI8i0hhoDhyw0c5gMBjs45574Msv9ehj2DAd48iGyqUq07Z6Wxr5NsLP14+1w9ay5Vs/mjWD33/PB5sLGJeYqlJKLQfmiMgvWbQpA+wG6jlaptZMVRkMhmz58EN4+mkYNUo7Ejf7f1dfjL6IZ0IF2rbV4ZIdO6BWrez7uTKFYarKD+islNqqlNqslGpjpU09IAKYr5TapZT6IkXuxCpKqbFKqR1KqR0RERF5ZbfBYCgqPPWULv23YAG8+qrd3RbtXkTDWQ25kBjGd9/p2a5BgyA6Oq8MLXjyzXEopX5RSu21sg1Eq/SWA9oDk4DlSmVO5McDaAl8KiItgBvAC7bOJyJzRKS1iLSuWLFi3lyUwWAoWkyeDGPGwNtvw9y5dnXpWqcrQ/2HUqtsLRo3hq++gpAQGD06b00tSByqAJgbRMRmLS2l1Djgm5QpqG1KqWR0idr0Q4XTwGkR2ZryeiVZOA6DwWBwGKVg9mw4eRLGjdPzTX36ZNmlVtlafNLvEwBiE2Pp2iuJ6dNLUr16fhhcMOSb48iG1UB3tP6VH1AMuJi+gYicV0qdUko1EpFQoAdGlddgMDib9Iq6992ns8ubN8+2W1JyEn0W98HH24fVz6y2rMC6caPo1fBwlRjHPKCeUmovsBQYKSKilKqmlEqvg/wksEQp9S8QBLyd/6YaDIYiT+nSsG6dLgDVrx+cPp1tF3c3d4b5D2NEwAiL0/j2W52AHlbEJKxcYlVVXmNWVRkMhhzx779a16puXb3Otoz9uq6XYi5x9UJ52rTR+op//639UGGhMKyqMhgMBtcjMFBPW+3bB4MHQ0JCxv3pij/RrJl+DWw/s526H9bl37g1rFypdRWHD7+5emBhxTgOg8FgyIo+feCzz3QBqCeeyFipacAAXfwpKUk/DhgAgH8lf4Y0G0Lraq3p0gU++AC+/14v2ioKuEpw3GAwGFyX0aPh2DG9TLdePXjxRf1+puJPhIYCUNyzOHMGzAG0vtUDD19m9+7yKKX9zk3JBoUM4zgMBoPBHt54QzuPl17SNTyGDtV1O9KVm00t/pSeST9PYm3YWrZ9tJ2y3jpGUrkyhIdjeZ6q7F5YMI7DYDAY7MHNDebP1yusRo2C6tV1safMxZ8yMbDRQEoVK0UZr9IABAenOQ2ACxfyyX4nYhyHwWAw2IuXF6xeDbfdpnVFgoN14DwLOtfuTOfanQG4HHMZd4+UpVWVd0NCSbhU+BR1TXDcYDAYHKF8eV0EysMD7rwz4/AhCyKjI2k1pxULI56EO5+GCy2g8bf4+hW+JA/jOAwGg8FR6tXT01Lnz8Ndd9mlaFi+eHl61evF4n8XQ9uPofkC+GsS8Z1fJCyycDkP4zgMBoMhJ7RrB0uWwLZt8OCDaaurbKCU4uy1s1yPvw5KoN84qLSPazsHFLo65cZxGAwGQ065+26YPh2++UavusqGGX1m4OPto18Ui4VRXSk35LlCV6fcOA6DwWDIDc88AyNH6loe336bZVM/Xz+CHw2mrJcOkJcrL/w95i+8b/jx6af5YKuTMI7DYDAYcoNSOrO8TRs9ZbV3b5bN/Xz92DZmG/0b9ufv0X9T16cusz5O5okn4Lvv8snmXGIch8FgMOQWb2892ihdGgYOhEuXsmyeWqe8Zpma9FncB89eU2ndWg9cUuSuXBrjOAwGg8EZVK+uYx2nT8OQIZCYmG0Xbw9vmlZsSpNKDVm+XA9eBg+G2Nh8sDcXGMdhMBgMzqJDB/j0U/jlF3j++WybK6X4uO/HPNj8QerWhQULhJ074d1388HWXGAch8FgMDiTRx6BJ5+EmTNh0SK7u/169Fdev9Caz/93kWefzUP7nIBxHAaDweBsZsyAbt1g7Fid52EHpb1K4+nmyR0DoildWucUnjyZx3bmEFMB0GAwGPKCixf1Sqv4eNixA6pWzbaLiFjKzvbsqRPTt24tmJrlpgKgwWAw5DcVKmhBxKgouPdeiIvLtotSioSkBCasn0DQw/PZv//m2lGugHEcBoPBkFc0bw4LFmgV3fHj7fIASikOXjyIV5XDTJmiwyRffpn3pjqCkVU3GAyGvOT+++Hll+Gtt6BFC62om7mGR716luYebh78MPwHPN09SeoKf/4JEyZA69YQFFRgV5EBM+IwGAyGvOb117WzeOop6N7dap3y9Hi6ewJwJCoMn4cfoluvWEqVym+jbWMch8FgMOQ1bm6weDH4+cGJE1brlFtjX/g+Np3+kelfHqFBA9eJdRjHYTAYDPlBmTJajMot3deujTrlqdzd5G4O/+cwzSo14/p1nVW+ZEk+2JoNxnEYDAZDftGwIXzxhX6ulM065ekp41UGgEX7P+dg7GbGjYNjx/La0KwxjsNgMBjyk4cfhjff1PNO//lPhsC4LWITY/lw2/v4DZmHUjB8uF1SWHnGLZsAmJCQwOnTp4l1dTUxQ5HC29ubGjVq4OnpWdCmGAqS5GTo3x9+/RX++gtatcqyeZUqcCH6LFyvTJnS7ly9CpMnw2uv5Z2JWSUA3rKO49ixY5QuXRpfX19LpqbBkJeICJGRkVy7do26desWtDmGgiYyUi/P9fCAnTuhXDmbTVWFMOg9ETbMgBuV6PzwT4R9M4RDh7SSe15gMsetEBsba5yGIV9RSuHr62tGuQaNry8sXw6nTunpKxs/4sMiw+DRDtBwvX7s9Txbqz3Eut9P5ZnTyI5b1nEAxmkY8h3zN2fIQPv2MG2aXm31/vs37Q6LDKPDlx3A+zK4JevHZitZcs8SWjesSVKS7prfE0e3tOMwGAyGAuepp7SW1f/9n04TT8fEnyYSFRsFbimewU1wK36FhSELAXjzy38YdE8Cixfnr8ku4ziUUk8qpUKVUvuUUu/ZaPNMyv69SqmvlVLe+W2ns4iNjaVt27Y0b96cZs2aMWXKlII2yWHi4uLo2bMnQUFBLFu2jNGjR7N//34A3n77bZv9stqXW9LbkFtGjRrFypUrs2yzYMECzp4965TzGW5RlNJiVHXr6sqBERGWXTP6zMDH2weFHqkqFD7ePszoM4PDlw7z5oV21B7+HuPHw5Ej+WiziBT4BnQDfgG8Ul5XstKmOnAMKJ7yejkwyp7jt2rVSjKzf//+m97LT5KTk+XatWsiIhIfHy9t27aV4OBgp58nISHB6cdMJTg4WG6//Xar+0qWLGmzX1b7ckNiYqJTjzdy5EhZsWJFlm26dOki27dvd+i4Bf23Z3BRdu0S8fIS6dVLJN3fcujFUCn/bnlxe81Nyr9bXkIvhlr2zftnnuw5FCU+PiLt2onExzvPHGCH2PhOdZURxzjgHRGJAxCRcBvtPIDiSikPoATgnJ96Tz8NXbs6d3v66SxPqZSiVIr4TEJCAgkJCVbnv48fP07jxo0ZOXIkgYGB3HfffURHRwPw+uuv06ZNG/z9/Rk7dmyqg6Vr16689NJLdOnShQ8//JC1a9fSrl07WrRoQc+ePblw4QIAU6dOZeTIkfTu3Zs6derwzTff8PzzzxMQEMAdd9xBQkKCTfvDw8MZMWIEISEhBAUFceTIEbp27cqOHTt44YUXiImJISgoiOHDh2foZ23f+++/j7+/P/7+/nzwwQdWz/f1118TEBCAv78///d//2d5v1SpUkyePJl27doRHBxssQHgyy+/xM/Pj65duzJmzBgmTJiQ5T0RESZMmEDTpk3p168f4eFpf4bWPuuVK1eyY8cOhg8fTlBQEDExMTbvicGQLUFBMGsW/PyzFkRMwc/Xj+BHg+nboC/Bjwbj5+tn2fdwi4fxb1CWTz9LYuueCF5/PZ9steVR8nMDQoDXgK3AZqCNjXZPAdeBCGBJNsccC+wAdtSqVesmb5rhV99TT4l06eLc7amnsnLmIqJ/ITdv3lxKliwpzz//vNU2x44dE0D++OMPERF5+OGHZdq0aSIiEhkZaWk3YsQIWbNmjYjoX8Hjxo2z7Lt06ZIkJyeLiMjcuXPl2WefFRGRKVOmSMeOHSU+Pl5CQkKkePHisn79ehERGTRokHz77bdZ2r9x40bp16+f5XX6X9/2jjh27Ngh/v7+cv36dbl27Zo0bdpU/vnnnwztz5w5IzVr1pTw8HBJSEiQbt26WWwDZNmyZTfZcObMGaldu7ZERkZKfHy8dOrUScaPH5/l9axatUp69uwpiYmJcubMGSlbtqxlxJHVZ51+xGGrXXrMiMNgk+RkkQcfFFFK5Oef7e42ZMUQqTylhWz6Pc5pppDFiCPfZNWVUr8AVazsehk9kigHtAfaAMuVUvVSjE/tXw4YCNQFooAVSqkRImI1LCQic4A5oPM4sjTOxq/cvMbd3Z2QkBCioqK4++672bt3L/7+/je1q1mzJh07dgRgxIgRfPTRRzz33HNs3LiR9957j+joaC5dukSzZs0YkKK0OWTIEEv/06dPM2TIEM6dO0d8fHyGHII777wTT09PAgICSEpK4o477gAgICCA48eP5+HVa/744w/uvvtuSqaUOLvnnnv4/fffadGihaXN9u3b6dq1KxUrVgRg+PDhbNmyhUGDBuHu7s69995703G3bdtGly5dKF++PAD3338/YWFhWdqyZcsWhg4diru7O9WqVaN79+6WfVl91umxt53BYBWl4NNP4Z9/YNgw2LULqlfPttuIwBH0rHee21voxNLk5IySWM4m36aqRKSniPhb2b4DTgPfpDi6bUAyUCHTIXoCx0QkQkQSgG+A2/LL/rzEx8eHrl278uOPP3Lq1CmCgoIICgris88+A25ewqmUIjY2lieeeIKVK1eyZ88exowZkyE/oGS6WpNPPvkkEyZMYM+ePXz++ecZ2nl5eQHg5uaGp6en5Vxubm4k5oOmQbrfBjlq4+3tjbu7e46Oaw1r04XZfdaOtjMYsqRkSVixQhcdf+AByGLKOJX+fv0Z3XI0Sime+78Ehg3L2yW6rhLjWA10B1BK+QHFgIuZ2pwE2iulSij9390DOJCfRjqTiIgIoqKiAIiJieGXX36hcePG1KxZk5CQEEJCQnj88ccBOHnyJMHBwYCe6+/UqZPlC6lChQpcv349y9U/V65coXrKr5aFCxfm4VWl4enpaTNGkn7f7bffzurVq4mOjubGjRt8++23dO7cOUP7du3asXnzZi5evEhSUhJff/01Xbp0yfL8bdu2ZfPmzVy+fJnExERWrVqVrc233347S5cuJSkpiXPnzrFx40aALD/r0qVLc+3atWzbGQwO0aQJzJkDf/wBr7xid7fgU8HMK+HHsl8PsGBB3pnnKhUA5wHzlFJ7gXhgpIiIUqoa8IWI9BWRrUqplcA/QCKwi5SpqMLIuXPnGDlyJElJSSQnJzN48GD69+9vtW2TJk1YuHAhjz32GA0bNmTcuHGUKFGCMWPGEBAQQJ06dWjTpo3Nc02dOpX777+f6tWr0759e445KK05efJkWrduzV133WV3n7FjxxIYGEjLli1ZkkkHOvO+UaNG0bZtW0Avp00/TQVQtWpV/vvf/9KtWzdEhL59+zJw4MAsz1+9enVeeukl2rVrR7Vq1WjatClly5YFYM2aNezYsYPXM0US7777bn777TcCAgLw8/OzOCcfHx+bn/WoUaN4/PHHKV68OMHBwXbfE4MhW4YNg99/h/feg44dwY7/v9o+tQmqU48bbdx48kno1EkL8jqbW1ar6sCBAzRp0qSALLKf48eP079/f/bu3VvQphQ6rl+/TqlSpUhMTOTuu+/mkUce4e677y5oswrN357BBYiN1U7j6FEd97BT4+z0aQgMhOeeg5deytmpjVaV4ZZk6tSpBAUF4e/vT926dRk0aFBBm2QwOIa3t453iOja5XFxdnWrXDWB++dOpGKfuXlilqtMVRlsUKdOHTPayCHTp08vaBMMhtxTrx4sWAB3362HDzNmZNvFTbmR7HmVVtWylmvPKcZxGAwGgytz9Ci8/LJeqvv++3oOauTILLu4u7kz9668GW2AmaoyGAwG12bAADh4MG197ejRGfSsCgLjOAwGg8GVCQ3VGX2pJCbCI4/kv5Z6OozjMBgMBlemUaO0NHA3N6hcGdat0xnmBYRxHAWIu7u7JUs8KCgoXyQ+nEF6IcG+ffsSFRVFVFQUn3zyicPH2rRpk838FVtMmjSJZs2aMWnSJD777DMWLVoEGIlzQxFl7Vpo3Bjc3fXjn3/CHXfAxImwb1+BmGSC4wVI8eLFCQkJsbovVUzMLS8FZ5zA+vXrAZ1v8sknn/DEE0/k+Tk///xzIiIiLHIpqSxYsAB/f3+qVauW5zYYDPlGvXo3O4gFC3SQfOhQ2LZNL9vNR1z7WykfsaaMnvoDOjra+v7UlP6LF2/elxOOHz9OkyZNeOKJJ2jZsiWnTp1i2rRptGnThsDAwAzFnhYvXkzbtm0JCgriscceIykp6abjZSW7/swzz3D77bfTpEkTtm/fzj333EPDhg15JUXeICs59/TUqVOHixcv8sILL3DkyBGCgoKYNGnSTSOJCRMmsCDlA/vxxx9p3LgxnTp14ptvvrG0uXHjBo888ght2rShRYsWfPfddzed76677uLGjRu0a9eOZcuWMXXqVKZPn25V4txgKLJUrgzz58OePfDCC/l+euM4CpDUuhRBQUGWjObQ0FAeeughdu3aRWhoKIcOHWLbtm2EhISwc+dOtmzZwoEDB1i2bBl//vknISEhuLu73yTrAfrLevv27ezdu5eYmBjWrVtn2VesWDG2bNnC448/zsCBA5k9ezZ79+5lwYIFREZGWmwZO3Ys//77L2XKlMlyKuqdd96hfv36hISEMG3aNJvtYmNjGTNmDGvXruX333/n/Pnzln1vvfUW3bt3Z/v27WzcuJFJkyZx48aNDP3XrFljGamlVwC+7777aN26NUuWLCEkJITixYtn8+kbDIWcvn3hySfhww/hxx/z9dRmqiqFTZts7ytRIuv9FSpkvd8Wmaeqjh8/Tu3atWnfvj0AGzZsYMOGDRbtpuvXr3Po0CH+/fdfdu7cadFCiomJoVKlSjcdPyuJ71TdqYCAAJo1a0bVqlUBqFevHqdOncLHx8emnHtuOHjwIHXr1qVhioDOiBEjmDNnjuV616xZY0nci42N5eTJk0aew2CwxXvvwcaNOq/j33/1SCQfMI7DxUgvhy4ivPjiizz22GMZ2syaNYuRI0fy3//+1+ZxUiW+d+zYQc2aNZk6dapNOfX0sYL0curW5NztxcPDg+R0SwjTn9vWcUSEVatW0ahRI7vPYzDc0nh7w9dfQ+vWeonuunU6UTCPMVNVLkyfPn2YN28e169fB+DMmTOEh4fTo0cPVq5caSlteunSJU6cOJGhrzMkvq3Judsivbw4QO3atdm/fz9xcXFcuXKFX3/9FYDGjRtz7Ngxjhw5Yjlu+uudNWuWJRaza9cuh+zNbIPBcEvg7w/Tp8P69fDxx5a3KzQKQw0bgKoQRhVrJfRygXEcLkzv3r0ZNmwYHTp0ICAggPvuu49r167RtGlT3nzzTXr37k1gYCC9evXi3LlzGfqmlwIfNGhQjiS+U+XcAwMDuXTpEuPGjbPZ1tfXl44dO+Lv78+kSZOoWbMmgwcPJjAwkOHDh1um27y9vZkzZw79+vWjU6dO1K5d23KMV199lYSEBAIDA/H39+fVV191yN5UiXMTHDfccowfD/36waRJsGcPYZFhRA7qAA3Xw6MduJCYdfVLRzGy6garGDn3vMP87RnyhPBwCAwkrlwZ6o24yNn4KHATSFYQW47Q54Px8/Wz+3BGVt1gMBiKOpUqwcKFeB08xItrLmunAfrRO4qJP0102qmM4zBYxci5GwyFkD59uPz4KCZsg34ps1MKRfmSPszok70cu70Yx2EwGAxFiHIffEasf2Pmf6eoel1Rrng5gh91bJoqO4zjMBgMhqKElxfey1bhm1SMtRsqEPzwn051GmAch8FgMBQ9mjbFbeYHtGrQGb+StZx+eOM4DAaDoSjy2GOwcqWWvnAyxnEUIHktq/7BBx9kECZMlUAHKFWqFKCX3X711Vc2j2EkzA2GQopSeZZFbiRHCpC8llX/4IMPGDFiBCVSfnGkSqCnJ9VxDBs2zOoxjIS5wWDIjBlxpNB1QVcWhCwAICEpga4LurL438UARCdE03VBV5btXQbAldgrdF3QlW8OaEnwi9EX6bqgK2tD1wJw/vr5m09gB86UVf/oo484e/Ys3bp1o1u3bkCaBHp6XnjhBX7//XeCgoKYOXNmhn1GwtxgMFjDOI4CJC9l1f/zn/9QrVo1Nm7cyMaNG23a8M4779C5c2dCQkJ45plnMuwzEuYGg8EaZqoqhU2jNlmee7p7ZnhdwrNEhtdlvctmeF2hRIUMr6uUsk9RLK9l1Q0GgyEvMI7DxXCWrLrBYDDkFWaqyoXJjaw62CcznlMpciNhbjC4MEePQrNm4OGhH48ederhjeNwYXIjqw4wduxY7rzzTktw3BqBgYF4eHjQvHnzm4LjWWEkzA0GF2bAADh4EJKS9GNK5U9nYWTVDYZ8xvztGfIcDw/tNFJxd4eUyp724vKy6kqpZUqpkJTtuFIqxEa7O5RSoUqpw0qpF/LZTIPBYCgcNGoEqTlgbm76tRNxCcchIkNEJEhEgoBVwDeZ2yil3IHZwJ1AU2CoUqppvhpqMBgMhYG1a6FxYz3SaNxYv3YiLrWqSimlgMFAdyu72wKHReRoStulwEBgf07PJyKofCjsbjCkcitMDRtcgHr1YN++PDu8S4w40tEZuCAih6zsqw6cSvf6dMp7VlFKjVVK7VBK7YiIiLhpv7e3N5GRkeYf2ZBviAiRkZF4e3sXtCkGQ67ItxGHUuoXwFpm3Msi8l3K86HA17YOYeU9m9/6IjIHmAM6OJ55f40aNTh9+jTWnIrBkFd4e3tTo0aNgjbDYMgV+eY4RKRnVvuVUh7APUArG01OAzXTva4B5Fie1dPTk7p16+a0u8FgMNyyuNJUVU/goIictrF/O9BQKVVXKVUMeABYk2/WGQwGgwFwLcfxAJmmqZRS1ZRS6wFEJBGYAPwEHACWi0jeRX8MBoPBYBWXWVUlIqOsvHcW6Jvu9Xrg5qISBoPBYMg3bonMcaVUBHCzmFMaFYCLWewvbJjrcW3M9bguRelaIHfXU1tEKlrbcUs4juxQSu2wlVpfGDHX49qY63FditK1QN5djyvFOAwGg8FQCDCOw2AwGAwOYRyHZk5BG+BkzPW4NuZ6XJeidC2QR9djYhwGg8FgcAgz4jAYDAaDQxjHYTAYDAaHuKUch1JqnlIqXCm118Z+pZT6KKVQ1L9KqZb5baO92HEtXZVSV9IVyJqc3zY6glKqplJqo1LqgFJqn1LqKSttCtP9sed6CsU9Ukp5K6W2KaV2p1zLa1baFKZ7Y8/1FIp7kx6llLtSapdSap2Vfc69PyJyy2zA7UBLYK+N/X2BH9BKvO2BrQVtcy6upSuwrqDtdOB6qgItU56XBsKApoX4/thzPYXiHqV83qVSnnsCW4H2hfje2HM9heLeZLL5WeAra3Y7+/7cUiMOEdkCXMqiyUBgkWj+BnyUUlXzxzrHsONaChUick5E/kl5fg2tR5a53kphuj/2XE+hIOXzvp7y0jNly7yqpjDdG3uup1ChlKoB9AO+sNHEqffnlnIcduBQsahCQIeU4fgPSqlmBW2MvSil6gAt0L8E01Mo708W1wOF5B6lTIOEAOHAzyJSqO+NHdcDheTepPAB8DyQbGO/U++PcRwZcahYlIvzD1prpjkwC1hdsObYh1KqFLru/NMicjXzbitdXPr+ZHM9heYeiUiSiASh6+C0VUr5Z2pSqO6NHddTaO6NUqo/EC4iO7NqZuW9HN8f4zgy4tRiUQWJiFxNHY6LVhX2VEpVKGCzskQp5Yn+kl0iIt9YaVKo7k9211MY75GIRAGbgDsy7SpU9yYVW9dTyO5NR+AupdRxYCnQXSm1OFMbp94f4zgysgZ4KGUFQnvgioicK2ijcoJSqopSSqU8b4u+15EFa5VtUmz9EjggIu/baFZo7o8911NY7pFSqqJSyifleXFSiq5lalaY7k2211NY7g2AiLwoIjVEpA66rtFvIjIiUzOn3h+XqceRHyilvkavlqiglDoNTEEHxhCRz9C1PvoCh4Fo4OGCsTR77LiW+4BxSqlEIAZ4QFKWV7goHYEHgT0pc88ALwG1oPDdH+y7nsJyj6oCC5VS7ugv0OUisk4p9TgUyntjz/UUlntjk7y8P0ZyxGAwGAwOYaaqDAaDweAQxnEYDAaDwSGM4zAYDAaDQxjHYTAYDAaHMI7DYDAYDA5hHIfBYDAYHMI4DoPBYDA4hHEcBkM+oZS6bkebkkqpj1Oyew0Gl8Q4DoPBtXgc8AI6FbQhBoMtjOMwGPIZpdSrSqmDSqmflVJfK6WeS7f7DnTRp5CCsc5gyB7jOAyGfEQp1Rq4F12f4x6gdbp93oA7urLj5gIx0GCwg1tK5NBgcAE6Ad+JSAyAUmptun0N0Y7joIgkFIRxBoM9GMdhMOQv1grqpFIR8EOX+TQYXBYzVWUw5C9/AAOUUt4p1QH7pdtXDV34yU0pVa5ArDMY7MA4DoMhHxGR7eiiOruBb4AdwBWllAc6tlEF+AxIKjAjDYZsMPU4DIZ8RilVSkSuK6VKAFuAsSLyT0HbZTDYi4lxGAz5zxylVFPAG1honIahsGFGHAaDwWBwCBPjMBgMBoNDGMdhMBgMBocwjsNgMBgMDmEch8FgMBgcwjgOg8FgMDiEcRwGg8FgcAjjOAwGg8HgEP8PHue+u5P687UAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "msize = 8\n",
    "d1, = plt.plot(logell, logCell + logell, 'r.', markersize=msize)\n",
    "d2, = plt.plot(logell_ps, logCell_ps + logell, 'bs', markersize=msize-5)\n",
    "d3, = plt.plot(logell_pO, logCell_pO + logell, 'gD', markersize=msize-4)\n",
    "\n",
    "f1, = plt.plot(logell, quadratic(logell, *popt), 'r-')\n",
    "f2, = plt.plot(logell_ps, quadratic_c(logell_ps, *popt_ps_A), 'b--')\n",
    "f3, = plt.plot(logell_pO, quadratic_u0(logell_pO, *popt_pO_x0), 'g:')\n",
    "\n",
    "#plt.plot(logell_ps, quadratic(logell_ps, *popt_ps), 'b--')\n",
    "#plt.plot(logell_pO, quadratic(logell_pO, *popt_pO), 'g--')\n",
    "\n",
    "plot_data  = [d1, d2, d3]\n",
    "label_data = ['Original data', 'Larger amplitude', 'Larger tilt']\n",
    "leg_data   = plt.legend(plot_data, label_data, loc=1)\n",
    "\n",
    "plot_fit  = [f1, f2, f3]\n",
    "label_fit = ['3-param. fit to orig. data', 'Free amplitude fit', 'Free tilt fit']\n",
    "leg_fit   = plt.legend(plot_fit, label_fit, loc=3)\n",
    "\n",
    "plt.gca().add_artist(leg_data)\n",
    "\n",
    "text_x = plt.xlabel('$\\lg \\ell$')\n",
    "text_y = plt.ylabel('$\\lg(\\ell C_\\ell)$')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}